diff --git a/diversity.cabal b/diversity.cabal --- a/diversity.cabal +++ b/diversity.cabal @@ -2,7 +2,7 @@ -- documentation, see http://haskell.org/cabal/users-guide/ name: diversity -version: 0.6.0.2 +version: 0.6.2.0 synopsis: Return the diversity at each position for all sequences in a fasta file description: Find the diversity of a collection of entities, mainly for use with fasta sequences. Produces a binary which works on fasta files to find the diversity of any order and rarefaction curves for a sliding window across all positions in the sequences. To analyze just a collection of entities, just use the whole sequences and list flag. homepage: https://github.com/GregorySchwartz/diversity @@ -20,7 +20,7 @@ hs-source-dirs: src/src-lib exposed-modules: Math.Diversity.Types, Math.Diversity.Diversity, Math.Diversity.GenerateDiversity, Math.Diversity.Print, Math.Diversity.Statistics, Math.Diversity.RandomSampling -- other-modules: - build-depends: base >=4.6 && <4.9, containers >=0.5 && <0.6, split >=0.2 && <0.3, parsec >=3.1 && <4.0, fasta >=0.6.1 && <0.7, math-functions >=0.1 && <0.2, scientific >= 0.3 && < 0.4, random-shuffle >=0.0 && < 0.1, MonadRandom >=0.3 && <0.4 + build-depends: base >=4.6 && <4.9, containers >=0.5 && <0.6, split >=0.2 && <0.3, parsec >=3.1 && <4.0, fasta >=0.6.1 && <0.7, math-functions >=0.1 && <0.2, scientific >= 0.3 && < 0.4, random-shuffle >=0.0 && < 0.1, MonadRandom >=0.3 && <0.4, data-ordlist >= 0.4 && < 0.5 executable diversity ghc-options: -O2 diff --git a/src/src-exec/Main.hs b/src/src-exec/Main.hs --- a/src/src-exec/Main.hs +++ b/src/src-exec/Main.hs @@ -5,11 +5,9 @@ window length (to split into fragments) by position. -} -{-# LANGUAGE BangPatterns #-} - -- Built-in import Data.List -import Control.Monad (forever) +import Control.Monad (forever, unless) import qualified Data.Map.Strict as Map import qualified System.IO as IO @@ -31,6 +29,7 @@ , inputFasta :: String , inputSampleField :: Int , inputSubsampling :: String + , inputG :: Double , fastBin :: Bool , runs :: Int , removeNBool :: Bool @@ -82,13 +81,25 @@ <*> strOption ( long "input-subsampling" <> short 'I' - <> metavar "INT INT" + <> metavar "INT INT (INT)" <> value "1 1" - <> help "The start point and interval of subsamples in the\ - \ rarefaction curve. For instance, '1 1' would be 1, 2, 3, ...\ - \ '2 6' would be 2, 8, 14, ... Note: input is a string so\ - \ use quotations around the entry and it always has the\ - \ number of subsamples overall as the last point" ) + <> help "The start point, interval, and optional endpoint of\ + \ subsamples in the rarefaction curve.\ + \ For instance, '1 1 4' would be 1, 2, 3, 4\ + \ '2 6 14' would be 2, 8, 14, ... Note: input is a string so\ + \ use quotations around the entry and it always includes the\ + \ number of subsamples overall in the result. Excluding the\ + \ endpoint results in the number of samples the endpoint, so\ + \ '1 1' would be 1, 2, 3, ..., N " ) + <*> option auto + ( long "input-g" + <> short 'g' + <> metavar "Double" + <> value 0.95 + <> help "Used for calculating the number of individuals\ + \ (or samples) needed before the proportion g of the total\ + \ number of estimated species is reached.\ + \ Sobs / Sest < g < 1" ) <*> switch ( long "fast-bin" <> short 'f' @@ -159,6 +170,12 @@ <> help "The csv file containing the diversities at each position.\ \ expects a string, so you need a string wven with std" ) +parseSampling :: (Num a, Read a) => String -> [a] +parseSampling = map read . parsing . words + where + parsing [x, y] = [x, y, "0"] + parsing x = x + pipesPositionMap :: Options -> IO PositionMap pipesPositionMap opts = do h <- if null . inputFasta $ opts @@ -188,33 +205,31 @@ window = inputWindow opts nFlag = removeNBool opts whole = wholeSeq opts - start = read . head . words . inputSubsampling $ opts - interval = read . last . words . inputSubsampling $ opts + [start, interval, end] = parseSampling . inputSubsampling $ opts howToOutput x = if std opts then putStrLn else writeFile x positionMap <- pipesPositionMap opts - if (null . output $ opts) - then return () - else howToOutput (output opts) - . printDiversity label order window - $ positionMap - if (null . outputRarefaction $ opts) - then return () - else do + unless (null . output $ opts) + $ howToOutput (output opts) + . printDiversity label order window + $ positionMap + unless (null . outputRarefaction $ opts) + $ do s <- printRarefaction (sample opts) (fastBin opts) (runs opts) start interval + end + (inputG opts) label window positionMap howToOutput (outputRarefaction opts) s - if (null . outputRarefactionCurve $ opts) - then return () - else do + unless (null . outputRarefactionCurve $ opts) + $ do s <- printRarefactionCurve (rarefactionDF opts) (sample opts) @@ -222,6 +237,7 @@ (runs opts) start interval + end label window positionMap diff --git a/src/src-lib/Math/Diversity/Diversity.hs b/src/src-lib/Math/Diversity/Diversity.hs --- a/src/src-lib/Math/Diversity/Diversity.hs +++ b/src/src-lib/Math/Diversity/Diversity.hs @@ -7,19 +7,27 @@ {-# LANGUAGE BangPatterns #-} module Math.Diversity.Diversity ( hamming + , richness , diversity , diversityOfMap + , chao1 + , chao2 , rarefactionCurve , rarefactionSampleCurve - , rarefactionViable ) where + , rarefactionViable + , individualG + , sampleG + , minRarefaction ) where -- Built-in import Data.List import qualified Data.Set as Set import qualified Data.Map.Strict as Map import Numeric.SpecFunctions (choose) -import Data.Function (on) +-- Cabal +import qualified Data.List.Ordered as LO + -- Local import Math.Diversity.RandomSampling import Math.Diversity.Types @@ -39,7 +47,7 @@ | otherwise = (fromInteger x / fromInteger y) * productDivision acc xs ys -- | Returns the diversity of a list of things -diversity :: (Ord b) => Double -> [b] -> Double +diversity :: (Ord a) => Double -> [a] -> Double diversity order sample | null sample = 0 | order == 1 = exp . h . speciesList $ sample @@ -57,8 +65,9 @@ ((Map.foldl' (+) 0 . speciesList $ sample) :: Double) speciesList = Map.fromListWith (+) . map (\x -> (x, 1)) --- | Returns the diversity of a map of things -diversityOfMap :: Double -> Map.Map (Sample, Fragment) Int -> Double +-- | Returns the diversity of a map of the species and how many times it +-- appears +diversityOfMap :: (Ord a) => Double -> Map.Map a Int -> Double diversityOfMap order sample | Map.null sample = 0 | order == 1 = exp . h $ sample @@ -70,11 +79,50 @@ pow = 1 / (1 - order) h = negate . Map.foldl' (+) 0 - . Map.map ( \x -> p_i (fromIntegral x) - * log (p_i (fromIntegral x))) - p_i x = (fromIntegral x :: Double) / - (fromIntegral (Map.foldl' (+) 0 sample) :: Double) + . Map.map ( \x -> p_i x * log (p_i x)) + p_i x = fromIntegral x / fromIntegral (Map.foldl' (+) 0 sample) +-- | Returns the richness of the observed data +richness :: (Ord a, Ord b) => Map.Map (a, b) c -> Int +richness = Map.size . Map.mapKeys snd + +-- | Returns the map of species with how many samples they appear in +overlapSampleMap :: (Ord a, Ord b) => Map.Map (a, b) Int -> Map.Map b Int +overlapSampleMap = Map.mapKeysWith (+) snd . Map.map (const 1) + +-- | Returns the number of a that appear x times +abundanceFreq :: (Ord a) => Int -> Map.Map a Int -> Int +abundanceFreq x = Map.size . Map.filter (== x) + +-- | Returns the number of b that appear in x number of a +overlapFreq :: (Ord a, Ord b) => Int -> Map.Map (a, b) Int -> Int +overlapFreq x = Map.size . Map.filter (== x) . overlapSampleMap + +-- | Returns the chao1 estimator of a map of the species and how many times +-- it appears +chao1 :: (Ord a) => Map.Map a Int -> Double +chao1 sample + | f2 > 0 = (f1 ** 2) / (2 * f2) + | otherwise = (f1 * (f1 - 1)) / (2 * (f2 + 1)) + where + f1 = fromIntegral . abundanceFreq 1 $ sample + f2 = fromIntegral . abundanceFreq 2 $ sample + +-- | Returns the chao2 estimator of a map of the sample labeled species +-- (sample, species) and how many times it appears. This will calculate the +-- overlap for you, so if you don't have the number of times it appears it +-- does not matter, you can set it to 1 and get the same result as it's all +-- about overlapping samples. +chao2 :: (Ord a, Ord b) => Map.Map (a, b) Int -> Double +chao2 samples + | q2 > 0 = ((t - 1) / t) * ((q1 ** 2) / (2 * q2)) + | otherwise = ((t - 1) / t) * ((q1 * (q1 - 1)) / (2 * (q2 + 1))) + where + q1 = fromIntegral . overlapFreq 1 $ samples + q2 = fromIntegral . overlapFreq 2 $ samples + -- Saves time so don't have to recalculate + t = fromIntegral . Map.size . Map.mapKeys fst $ samples + -- | Binomial for small or large numbers (slow but works for big numbers, -- fast but works for small numbers) specialBinomial :: Bool -> Integer -> Integer -> Integer -> Double @@ -91,25 +139,29 @@ -> Int -> Integer -> Integer + -> Integer -> Map.Map (Sample, Fragment) Int -> IO [(Int, (Double, Double))] -rarefactionCurve !fastBin !runs !start !interval !xs = - mapM rarefact $ [start,(start + interval)..(n_total - 1)] ++ [n_total] +rarefactionCurve !fastBin !runs !start !interval !end !sample = + mapM rarefact + . LO.nubSort + $ n_total : [start,(start + interval)..finish] where rarefact !n | n == 0 = return (fromIntegral n, (0, 0)) | n == 1 = return (fromIntegral n, (1, 0)) | n == n_total = return (fromIntegral n, (k, 0)) + | n > n_total = return (fromIntegral n, (estimation n, 0)) | runs == 0 = return (fromIntegral n, (k - inner n, 0)) - | otherwise = do + | otherwise = do -- Empirical version statTuple <- subsampleES runs (fromIntegral n_total) (fromIntegral n) . concatMap snd . Map.toAscList - . Map.mapWithKey (\(s, f) x -> replicate x f) - $ xs + . Map.mapWithKey (\(_, f) x -> replicate x f) + $ sample return (fromIntegral n, statTuple) inner n = ( \x -> if fastBin then x / choose (fromIntegral n_total) (fromIntegral n) @@ -117,9 +169,12 @@ . sum . map (\g -> specialBinomial fastBin n_total (fromIntegral g) n) $ grouped - n_total = fromIntegral . Map.foldl' (+) 0 $ xs - k = fromIntegral . Map.size $ xs - grouped = Map.elems xs + -- Unreadable unless I break the 80 column rule + estimation n = fromIntegral n_total + (chao1 sample * (1 - exp (((fromIntegral n - fromIntegral n_total) / (- fromIntegral n_total)) * (fromIntegral (abundanceFreq 1 sample) / chao1 sample)))) + finish = if end == 0 then n_total else end + n_total = fromIntegral . Map.foldl' (+) 0 $ sample + k = fromIntegral . Map.size $ sample + grouped = Map.elems sample -- | Each sample has a collection of species, return a list of these maps getSampleContents :: Map.Map (Sample, Fragment) Int -> [Set.Set Fragment] @@ -132,30 +187,37 @@ rarefactionSampleCurve :: Bool -> Int -> Int + -> Int -> Map.Map (Sample, Fragment) Int -> IO [(Int, (Double, Double))] -rarefactionSampleCurve !fastBin !start !interval !ls = - mapM rarefact $ [start,(start + interval)..(t_total - 1)] ++ [t_total] +rarefactionSampleCurve !fastBin !start !interval !end !ls = + mapM rarefact + . LO.nubSort + $ t_total : [start,(start + interval)..finish] where rarefact !t | t == 0 = return (t, (0, 0)) - | t == t_total = return (t, (richness, 0)) - | otherwise = return (t, (richness - inner t, 0)) - inner t = ( \x -> if fastBin - then x / choose t_total t - else x ) - . sum - . map ( \s -> specialBinomial - fastBin - (fromIntegral t_total) - (numHave s samples) - (fromIntegral t) ) - $ speciesList - numHave s = genericLength . filter (Set.member s) - richness = genericLength speciesList - speciesList = Map.keys . Map.mapKeys snd $ ls - t_total = genericLength samples - samples = getSampleContents ls + | t == t_total = return (t, (sobs, 0)) + | t > t_total = return (t, (estimation t, 0)) + | otherwise = return (t, (sobs - inner t, 0)) + inner t = ( \x -> if fastBin + then x / choose t_total t + else x ) + . sum + . map ( \s -> specialBinomial + fastBin + (fromIntegral t_total) + (numHave s samples) + (fromIntegral t) ) + $ speciesList + -- Unreadable unless I break the 80 column rule + estimation t = sobs + (chao2 ls * (1 - exp ((- (fromIntegral t - fromIntegral t_total) * fromIntegral (overlapFreq 1 ls)) / (fromIntegral (overlapFreq 1 ls) + (fromIntegral t_total * chao2 ls))))) + finish = if end == 0 then t_total else end + numHave s = genericLength . filter (Set.member s) + sobs = fromIntegral $ richness ls + speciesList = Map.keys . Map.mapKeys snd $ ls + t_total = genericLength samples + samples = getSampleContents ls -- | Calculates the percent of the curve that is above 95% of height of the -- curve @@ -163,3 +225,54 @@ rarefactionViable !xs = (genericLength valid / genericLength xs) * 100 where valid = dropWhile (< (0.95 * last xs)) xs + +-- | Returns the number of individuals needed to get the proportion g of +-- the estimated total richness of the assemblage. Sobs / Sest < g < 1 +individualG :: Double -> Map.Map (Sample, Fragment) Int -> Double +individualG g sample = ((sobs * f1) / (2 * f2)) * log (f0 / ((1 - g) * sest)) + where + sest = sobs + f0 + sobs = fromIntegral $ richness sample + f2 = fromIntegral . abundanceFreq 2 $ sample + f1 = fromIntegral . abundanceFreq 1 $ sample + f0 = chao1 sample + +-- | Returns the number of samples needed to get the proportion g of +-- the estimated total richness of the assemblage. Sobs / Sest < g < 1 +sampleG :: Double -> Map.Map (Sample, Fragment) Int -> Double +sampleG g sample = log (1 - ((t / (t - 1)) * ((2 * q2) / (q1 ** 2)) * ((g * sest) - sobs))) + / log (1 - ((2 * q2) / (((t - 1) * q1) + (2 * q2)))) + where + t = fromIntegral . Map.size . Map.mapKeys fst $ sample + sest = sobs + q0 + sobs = fromIntegral $ richness sample + q2 = fromIntegral . overlapFreq 2 $ sample + q1 = fromIntegral . overlapFreq 1 $ sample + q0 = chao2 sample + +-- | Returns the number of samples needed before a new sample returns less +-- than x new species. Warning, goes forever until threshold is met! +minRarefaction :: Bool + -> Bool + -> Int + -> Map.Map (Sample, Fragment) Int + -> Double + -> Int + -> IO Int +minRarefaction _ _ (-1) _ _ _ = return (-1) +minRarefaction bySample fastBin threshold sample !oldRare !count = do + newRare <- fmap (fst . snd . head) + . rarefaction bySample count + $ sample + if newRare - oldRare < fromIntegral threshold + then return count + else minRarefaction True fastBin threshold sample newRare (count + 1) + where + rarefaction True x = rarefactionSampleCurve fastBin x 1 x + rarefaction False x = rarefactionCurve + fastBin + 0 + (fromIntegral x) + 1 + (fromIntegral x) + diff --git a/src/src-lib/Math/Diversity/Print.hs b/src/src-lib/Math/Diversity/Print.hs --- a/src/src-lib/Math/Diversity/Print.hs +++ b/src/src-lib/Math/Diversity/Print.hs @@ -45,28 +45,46 @@ -> Int -> Int -> Int + -> Int + -> Double -> Label -> Window -> PositionMap -> IO String printRarefaction - sample fastBin runs start interval label window positionMap = do + bySample fastBin runs start interval end g label window positionMap = do body <- fmap unlines . mapM mapLine . Map.toAscList $ positionMap return (header ++ body) where - header = "label,window,position,weight,percent_above\n" + header = "label,window,position,weight,percent_above,\ + \additional_sampling,g_proportion,richness,S_est\n" mapLine (p, xs) = fmap (intercalate ",") . line p $ xs line p xs = do - curve <- getRarefactionCurve sample xs + curve <- getRarefactionCurve bySample xs + -- The minimum number of samples needed before any additional + -- sampling returns less than the threshold (min) number of species return [ label , show window , show p , show . Map.foldl' (+) 0 $ xs , show . rarefactionViable . map (snd . snd) $ curve + , show . additionalSampling bySample $ xs + , show g + , show . sobs $ xs + , show . sest bySample $ xs ] - getRarefactionCurve True = rarefactionSampleCurve fastBin start interval + sest True xs = sobs xs + chao2 xs + sest False xs = sobs xs + chao1 xs + sobs = fromIntegral . richness + additionalSampling True = sampleG g + additionalSampling False = individualG g + getRarefactionCurve True = rarefactionSampleCurve fastBin start interval end getRarefactionCurve False = rarefactionCurve - fastBin runs (fromIntegral start) (fromIntegral interval) + fastBin + runs + (fromIntegral start) + (fromIntegral interval) + (fromIntegral end) -- Return the results of the rarefaction analysis of the entire curve in -- string form for saving to a file @@ -76,12 +94,13 @@ -> Int -> Int -> Int + -> Int -> Label -> Window -> PositionMap -> IO String printRarefactionCurve - asDF sample fastBin runs start interval label window positionMap = do + asDF bySample fastBin runs start interval end label window positionMap = do body <- fmap unlines . mapM mapLine . Map.toAscList $ positionMap return (header asDF ++ body) @@ -92,7 +111,7 @@ \expected_richness,mad\n" mapLine (!p, !xs) = line asDF p xs line False p xs = do - curve <- getRarefactionCurve sample xs + curve <- getRarefactionCurve bySample xs return . intercalate "," $ [ label , show window , show p @@ -105,7 +124,7 @@ $ curve ] line True p xs = do - curve <- getRarefactionCurve sample xs + curve <- getRarefactionCurve bySample xs return . intercalate "\n" . map ( \(!x, (!y, !z)) -> intercalate "," [ label @@ -117,6 +136,10 @@ , show z ] ) $ curve - getRarefactionCurve True = rarefactionSampleCurve fastBin start interval + getRarefactionCurve True = rarefactionSampleCurve fastBin start interval end getRarefactionCurve False = rarefactionCurve - fastBin runs (fromIntegral start) (fromIntegral interval) + fastBin + runs + (fromIntegral start) + (fromIntegral interval) + (fromIntegral end)