diversity 0.7.1.1 → 0.8.0.0
raw patch · 6 files changed
+91/−37 lines, 6 files
Files
- diversity.cabal +2/−2
- src/src-exec/Main.hs +1/−2
- src/src-lib/Math/Diversity/Diversity.hs +60/−16
- src/src-lib/Math/Diversity/Print.hs +13/−8
- src/src-lib/Math/Diversity/RandomSampling.hs +6/−1
- src/src-lib/Math/Diversity/Statistics.hs +9/−8
diversity.cabal view
@@ -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
src/src-exec/Main.hs view
@@ -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" )
src/src-lib/Math/Diversity/Diversity.hs view
@@ -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
src/src-lib/Math/Diversity/Print.hs view
@@ -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
src/src-lib/Math/Diversity/RandomSampling.hs view
@@ -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]
src/src-lib/Math/Diversity/Statistics.hs view
@@ -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