diversity 0.6.0.2 → 0.6.2.0

raw patch · 4 files changed

+222/−70 lines, 4 filesdep +data-ordlistPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: data-ordlist

API changes (from Hackage documentation)

Files

• diversity.cabal +2/−2
• src/src-exec/Main.hs +38/−22
• src/src-lib/Math/Diversity/Diversity.hs +149/−36
• src/src-lib/Math/Diversity/Print.hs +33/−10

@@ -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

@@ -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

@@ -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)+

@@ -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)