diff --git a/diversity.cabal b/diversity.cabal --- a/diversity.cabal +++ b/diversity.cabal @@ -2,8 +2,8 @@ -- documentation, see http://haskell.org/cabal/users-guide/ name: diversity -version: 0.7.1.1 -synopsis: Return the diversity at each position by default for all sequences in a fasta file +version: 0.8.0.0 +synopsis: Quantify the diversity of a population description: Find the diversity of a collection of entities, mainly for use with fasta sequences. homepage: https://github.com/GregorySchwartz/diversity license: GPL-2 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 @@ -265,5 +265,4 @@ where opts = info (helper <*> options) ( fullDesc - <> progDesc "Return the diversity at each position for all sequences in a\ - \ fasta file" ) + <> progDesc "Quantify the diversity of a population" ) 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 @@ -12,6 +12,8 @@ , diversityOfMap , chao1 , chao2 + , chao1Var + , chao2Var , rarefactionCurve , rarefactionSampleCurve , rarefactionViable @@ -94,7 +96,9 @@ 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 +-- | Returns the number of b that appear in x number of a. Notice that this +-- function takes in a normal frequency map, as it converts it with +-- overlapSampleMap. overlapFreq :: (Ord a, Ord b) => Int -> Map.Map (a, b) Int -> Int overlapFreq x = Map.size . Map.filter (== x) . overlapSampleMap @@ -123,6 +127,44 @@ -- Saves time so don't have to recalculate t = fromIntegral . Map.size . Map.mapKeys fst $ samples +-- | Returns the chao1 estimator variance of a map of the species +-- and how many times each one appears. +chao1Var :: (Ord a) => Map.Map a Int -> Double +chao1Var sample + | f2 > 0 = f2 + * ( ((1 / 2) * ((n - 1) / n) * ((f1 / f2) ** 2)) + + ((((n - 1) / n) ** 2) * ((f1 / f2) ** 3)) + + ((1 / 4) * (((n - 1) / n) ** 2) * ((f1 / f2) ** 4)) + ) + | otherwise = (((n - 1) / n) * ((f1 * (f1 - 1)) / 2)) + + ((((n - 1) / n) ** 2) * ((f1 * (((2 * f1) - 1) ** 2)) / 4)) + + ((((n - 1) / n) ** 2) * ((f1 ** 4) / (4 * sest))) + where + sest = fromIntegral (Map.size sample) + chao1 sample + f1 = fromIntegral . abundanceFreq 1 $ sample + f2 = fromIntegral . abundanceFreq 2 $ sample + -- Saves time so don't have to recalculate + n = fromIntegral . Map.foldl' (+) 0 $ sample + +-- | Returns the chao2 estimator variance of a map of the sample labeled species +-- (sample, species) and how many times it appears. +chao2Var :: (Ord a, Ord b) => Map.Map (a, b) Int -> Double +chao2Var samples + | q2 > 0 = q2 + * ( ((1 / 2) * ((t - 1) / t) * ((q1 / q2) ** 2)) + + ((((t - 1) / t) ** 2) * ((q1 / q2) ** 3)) + + ((1 / 4) * ((t - 1 / t) ** 2) * ((q1 / q2) ** 4)) + ) + | otherwise = (((t - 1) / t) * ((q1 * (q1 - 1)) / 2)) + + ((((t - 1) / t) ** 2) * ((q1 * (((2 * q1) - 1) ** 2)) / 4)) + + ((((t - 1) / t) ** 2) * ((q1 ** 4) / (4 * sest))) + where + sest = fromIntegral (richness samples) + chao2 samples + 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 @@ -141,18 +183,18 @@ -> Integer -> Integer -> Map.Map (Sample, Fragment) Int - -> IO [(Int, (Double, Double))] + -> IO [(Int, Maybe (Double, Double))] 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)) + | n == 0 = return (fromIntegral n, Just (0, 0)) + | n == 1 = return (fromIntegral n, Just (1, 0)) + | n == n_total = return (fromIntegral n, Just (k, 0)) + | n > n_total = return (fromIntegral n, Just (estimation n, 0)) + | runs == 0 = return (fromIntegral n, Just (k - inner n, 0)) | otherwise = do -- Empirical version statTuple <- subsampleES runs @@ -189,17 +231,17 @@ -> Int -> Int -> Map.Map (Sample, Fragment) Int - -> IO [(Int, (Double, Double))] + -> IO [(Int, Maybe (Double, Double))] 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, (sobs, 0)) - | t > t_total = return (t, (estimation t, 0)) - | otherwise = return (t, (sobs - inner t, 0)) + | t == 0 = return (t, Just (0, 0)) + | t == t_total = return (t, Just (sobs, 0)) + | t > t_total = return (t, Just (estimation t, 0)) + | otherwise = return (t, Just (sobs - inner t, 0)) inner t = ( \x -> if fastBin then x / choose t_total t else x ) @@ -265,12 +307,14 @@ -> IO Int minRarefaction _ _ (-1) _ _ _ = return (-1) minRarefaction bySample fastBin threshold sample !oldRare !count = do - newRare <- fmap (fst . snd . head) + newRare <- fmap (maybe (-1) fst . snd . head) . rarefaction bySample count $ sample - if newRare - oldRare < fromIntegral threshold - then return count - else minRarefaction True fastBin threshold sample newRare (count + 1) + if newRare == (-1) + then return (-1) + else 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 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 @@ -13,6 +13,7 @@ -- Built in import Data.List +import Data.Maybe import qualified Data.Map.Strict as Map -- Local @@ -52,7 +53,8 @@ return (header ++ body) where header = "label,window,position,weight,\ - \additional_sampling,g_proportion,richness,S_est\n" + \additional_sampling,g_proportion,richness,\ + \S_est,S_est_var\n" mapLine (p, xs) = fmap (intercalate ",") . line p $ xs line p xs = do -- The minimum number of samples needed before any additional @@ -65,9 +67,12 @@ , show g , show . sobs $ xs , show . sest bySample $ xs + , show . var bySample $ xs ] sest True xs = sobs xs + chao2 xs sest False xs = sobs xs + chao1 xs + var True = chao2Var + var False = chao1Var sobs = fromIntegral . richness additionalSampling True = sampleG g additionalSampling False = individualG g @@ -104,30 +109,30 @@ , show . Map.foldl' (+) 0 $ xs , show . rarefactionViable - . map (snd . snd) + . map (maybe (-1) snd . snd) $ curve , intercalate "/" - . map (show . fst . snd) + . map (maybe "NA" (show . fst) . snd) $ curve , intercalate "/" - . map (show . snd . snd) + . map (maybe "NA" (show . snd) . snd) $ curve ] line True p xs = do curve <- getRarefactionCurve bySample xs return . intercalate "\n" - . map ( \(!x, (!y, !z)) + . map ( \(!x, !y) -> intercalate "," [ label , show window , show p , show . Map.foldl' (+) 0 $ xs , show . rarefactionViable - . map (snd . snd) + . map (maybe (-1) snd . snd) $ curve , show x - , show y - , show z + , maybe "NA" (show . fst) y + , maybe "NA" (show . snd) y ] ) $ curve getRarefactionCurve True = rarefactionSampleCurve fastBin start interval end diff --git a/src/src-lib/Math/Diversity/RandomSampling.hs b/src/src-lib/Math/Diversity/RandomSampling.hs --- a/src/src-lib/Math/Diversity/RandomSampling.hs +++ b/src/src-lib/Math/Diversity/RandomSampling.hs @@ -34,7 +34,12 @@ -- | Repeat the sampling to get a median and MAD value for the runs for the -- expected species counts -subsampleES :: (Eq a, Ord a) => Int -> Int -> Int -> [a] -> IO (Double, Double) +subsampleES :: (Eq a, Ord a) + => Int + -> Int + -> Int + -> [a] + -> IO (Maybe (Double, Double)) subsampleES runs size n xs = do allRuns <- mapM (const run) [1..runs] diff --git a/src/src-lib/Math/Diversity/Statistics.hs b/src/src-lib/Math/Diversity/Statistics.hs --- a/src/src-lib/Math/Diversity/Statistics.hs +++ b/src/src-lib/Math/Diversity/Statistics.hs @@ -12,23 +12,24 @@ -- | Get the median, adapted from -- http://rosettacode.org/wiki/Averages/Median#Haskell, but FIXED due to --- zero indexing +-- zero indexing. Also made sure that it is sorted median :: [Double] -> Maybe Double median xs | null xs = Nothing - | odd len = Just $ xs !! mid + | odd len = Just $ sorted !! mid | even len = Just meanMedian where - len = length xs + len = length sorted mid = len `div` 2 - meanMedian = (xs !! (mid - 1) + xs !! mid) / 2 + meanMedian = (sorted !! (mid - 1) + sorted !! mid) / 2 + sorted = sort xs -- | Get the median and median absolute deviation (MAD) of a list of numbers -medmad :: [Double] -> (Double, Double) -medmad xs = (sampleMedian, mad) +medmad :: [Double] -> Maybe (Double, Double) +medmad [] = Nothing +medmad xs = Just (sampleMedian, mad) where mad = fromJust . median - . sort . map (\x -> abs (x - sampleMedian)) $ xs - sampleMedian = fromJust . median . sort $ xs + sampleMedian = fromJust . median $ xs