diversity 0.4.0.2 → 0.5.0.0

raw patch · 6 files changed

+201/−81 lines, 6 filesdep +MonadRandomdep +random-shufflePVP ok

version bump matches the API change (PVP)

Dependencies added: MonadRandom, random-shuffle

API changes (from Hackage documentation)

Files

• diversity.cabal +3/−3
• src/src-exec/Main.hs +37/−21
• src/src-lib/Math/Diversity/Diversity.hs +24/−17
• src/src-lib/Math/Diversity/Print.hs +56/−40
• src/src-lib/Math/Diversity/RandomSampling.hs +47/−0
• src/src-lib/Math/Diversity/Statistics.hs +34/−0

@@ -2,7 +2,7 @@ -- documentation, see http://haskell.org/cabal/users-guide/  name:                diversity-version:             0.4.0.2+version:             0.5.0.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@@ -18,9 +18,9 @@ library   ghc-options: -O2   hs-source-dirs:      src/src-lib-  exposed-modules:     Math.Diversity.Types, Math.Diversity.Diversity, Math.Diversity.GenerateDiversity, Math.Diversity.Print+  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.8, containers >=0.5 && <0.6, split >=0.2 && <0.3, parsec >=3.1 && <4.0, fasta >=0.5.1.3 && <0.6, math-functions >=0.1 && <0.2, scientific >= 0.3 && < 0.4+  build-depends:       base >=4.6 && <4.8, containers >=0.5 && <0.6, split >=0.2 && <0.3, parsec >=3.1 && <4.0, fasta >=0.5.1.3 && <0.6, math-functions >=0.1 && <0.2, scientific >= 0.3 && < 0.4, random-shuffle >=0.0 && < 0.1, MonadRandom >=0.3 && <0.4  executable diversity   ghc-options: -O2

@@ -24,6 +24,7 @@                        , inputSampleField       :: Int                        , inputSubsampling       :: String                        , fastBin                :: Bool+                       , runs                   :: Int                        , removeN                :: Bool                        , wholeSeq               :: Bool                        , list                   :: Bool@@ -87,6 +88,17 @@                  \ coefficient for the rarefaction analysis. This method\                  \ results in NaNs for larger numbers, so in that case you\                  \ you should use the slower, more accurate default method" )+      <*> option auto+          ( long "input-runs"+         <> short 'R'+         <> metavar "INT"+         <> value 0+         <> help "The number of runs for empirical resampling rarefaction.\+                 \ This method does not compute the theoretical, it reports the\+                 \ actual median and median absolute deviation (MAD) values of\+                 \ this many runs. If this value is not 0, empirical\+                 \ rarefaction is automatically enabled (individual based only,\+                 \ not for sample based)" )       <*> switch           ( long "remove-N"          <> short 'n'@@ -141,11 +153,11 @@  generateDiversity :: Options -> IO () generateDiversity opts = do-    contentsFile    <- if (null . inputFasta $ opts)+    contentsFile    <- if null . inputFasta $ opts                         then getContents                         else readFile . inputFasta $ opts     let lineIt       = unlines . (:) ">" . intersperse ">" . lines-        contents     = if (list opts) then lineIt contentsFile else contentsFile+        contents     = if list opts then lineIt contentsFile else contentsFile         label        = inputLabel opts         order        = inputOrder opts         window       = inputWindow opts@@ -172,27 +184,31 @@            $ positionMap     if (null . outputRarefaction $ opts)         then return ()-        else howToOutput (outputRarefaction opts)-           $ printRarefaction-             (sample opts)-             (fastBin opts)-             start-             interval-             label-             window-             positionMap+        else do+            s <- printRarefaction+                 (sample opts)+                 (fastBin opts)+                 (runs opts)+                 start+                 interval+                 label+                 window+                 positionMap+            howToOutput (outputRarefaction opts) s     if (null . outputRarefactionCurve $ opts)         then return ()-        else howToOutput (outputRarefactionCurve opts)-           $ printRarefactionCurve-             (rarefactionDF opts)-             (sample opts)-             (fastBin opts)-             start-             interval-             label-             window-             positionMap+        else do+            s <- printRarefactionCurve+                 (rarefactionDF opts)+                 (sample opts)+                 (fastBin opts)+                 (runs opts)+                 start+                 interval+                 label+                 window+                 positionMap+            howToOutput (outputRarefactionCurve opts) s  main :: IO () main = execParser opts >>= generateDiversity

@@ -14,11 +14,13 @@  -- Built-in import Data.List-import Data.Ratio import qualified Data.Set as Set import Numeric.SpecFunctions (choose) import Data.Function (on) +-- Local+import Math.Diversity.RandomSampling+ -- | Takes two strings, returns Hamming distance hamming :: String -> String -> Int hamming xs ys = length $ filter not $ zipWith (==) xs ys@@ -60,18 +62,23 @@ -- | Returns the rarefaction curve for each position in a list rarefactionCurve :: (Eq a, Ord a)                  => Bool+                 -> Int                  -> Integer                  -> Integer                  -> [a]-                 -> [(Int, Double)]-rarefactionCurve fastBin start interval xs =-    map rarefact $ [start,(start + interval)..(n_total - 1)] ++ [n_total]+                 -> IO [(Int, (Double, Double))]+rarefactionCurve !fastBin !runs !start !interval !xs =+        mapM rarefact $ [start,(start + interval)..(n_total - 1)] ++ [n_total]   where-    rarefact n-        | n == 0       = (fromIntegral n, 0)-        | n == 1       = (fromIntegral n, 1)-        | n == n_total = (fromIntegral n, k)-        | otherwise    = (fromIntegral n, k - inner n)+    rarefact !n+        | n == 0       = return (fromIntegral n, (0, 0))+        | n == 1       = return (fromIntegral n, (1, 0))+        | n == n_total = return (fromIntegral n, (k, 0))+        | runs == 0    = return (fromIntegral n, (k - inner n, 0))+        | otherwise    = do+            statTuple <-+                subsampleES runs (fromIntegral n_total) (fromIntegral n) xs+            return (fromIntegral n, statTuple)     inner n = ( \x -> if fastBin                         then x / choose (fromIntegral n_total) (fromIntegral n)                         else x )@@ -94,14 +101,14 @@                        -> Int                        -> Int                        -> [(a, b)]-                       -> [(Int, Double)]-rarefactionSampleCurve fastBin start interval ls =-    map rarefact $ [start,(start + interval)..(t_total - 1)] ++ [t_total]+                       -> IO [(Int, (Double, Double))]+rarefactionSampleCurve !fastBin !start !interval !ls =+    mapM rarefact $ [start,(start + interval)..(t_total - 1)] ++ [t_total]   where-    rarefact t-        | t == 0       = (t, 0)-        | t == t_total = (t, richness)-        | otherwise    = (t, richness - inner t)+    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 )@@ -120,6 +127,6 @@  -- | Calculates the percent of the curve that is above 95% of height of the curve rarefactionViable :: [Double] -> Double-rarefactionViable xs = (genericLength valid / genericLength xs) * 100+rarefactionViable !xs = (genericLength valid / genericLength xs) * 100   where     valid = dropWhile (< (0.95 * last xs)) xs

@@ -44,28 +44,29 @@                  -> Bool                  -> Int                  -> Int+                 -> Int                  -> Label                  -> Window                  -> PositionMap-                 -> String-printRarefaction sample fastBin start interval label window positionMap =-    header ++ body+                 -> IO String+printRarefaction+    sample fastBin runs start interval label window positionMap = do+    body <- fmap unlines . mapM mapLine . M.toAscList $ positionMap+    return (header ++ body)   where     header           = "label,window,position,weight,percent_above\n"-    body             = unlines-                     . map mapLine-                     . M.toAscList-                     $ positionMap-    mapLine (p, xs) = intercalate "," . line p $ xs-    line p xs  = [ label-                 , show window-                 , show p-                 , show . length $ xs-                 , show-                 . rarefactionViable . map snd . getRarefactionCurve sample $ xs-                 ]+    mapLine (p, xs)  = fmap (intercalate ",") . line p $ xs+    line p xs        = do+        curve <- getRarefactionCurve sample xs+        return [ label+               , show window+               , show p+               , show . length $ xs+               , show . rarefactionViable . map (snd . snd) $ curve+               ]     getRarefactionCurve True = rarefactionSampleCurve fastBin start interval-    getRarefactionCurve False = rarefactionCurve fastBin (fromIntegral start) (fromIntegral interval)+    getRarefactionCurve False = rarefactionCurve+        fastBin runs (fromIntegral start) (fromIntegral interval)  -- Return the results of the rarefaction analysis of the entire curve in -- string form for saving to a file@@ -74,35 +75,50 @@                       -> Bool                       -> Int                       -> Int+                      -> Int                       -> Label                       -> Window                       -> PositionMap-                      -> String-printRarefactionCurve asDF sample fastBin start interval label window positionMap =-    header asDF ++ body+                      -> IO String+printRarefactionCurve+    asDF sample fastBin runs start interval label window positionMap = do+    body <- fmap unlines . mapM mapLine . M.toAscList $ positionMap++    return (header asDF ++ body)   where-    header False     = "label,window,position,weight,curve\n"-    header True      = "label,window,position,weight,subsample,vertical_curve\n"-    body             = unlines-                     . map mapLine-                     . M.toAscList-                     $ positionMap+    header False      = "label,window,position,weight,expected_richness,\+                        \mad\n"+    header True       = "label,window,position,weight,subsample,\+                        \expected_richness,mad\n"     mapLine (!p, !xs) = line asDF p xs-    line False p xs   = intercalate "," [ label+    line False p xs   = do+        curve <- getRarefactionCurve sample xs+        return . intercalate "," $ [ label+                                   , show window+                                   , show p+                                   , show . length $ xs+                                   , intercalate "/"+                                   . map (show . fst . snd)+                                   $ curve+                                   , intercalate "/"+                                   . map (show . snd . snd)+                                   $ curve+                                   ]+    line True p xs    = do+        curve <- getRarefactionCurve sample xs+        return . intercalate "\n"+               . map ( \(!x, (!y, !z))+                     -> intercalate "," [ label                                         , show window                                         , show p-                                        , show . length $ xs-                                        , intercalate "/" . map (show . snd) . getRarefactionCurve sample $ xs-                                        ]-    line True p xs    = intercalate "\n"-                      . map ( \(!x, !y) -> intercalate "," [ label-                                                           , show window-                                                           , show p-                                                           , show . length $ xs-                                                           , show x-                                                           , show y-                                                           ] )-                      . getRarefactionCurve sample $ xs+                                        , show+                                        . length+                                        $ xs+                                        , show x+                                        , show y+                                        , show z+                                        ] )+               $ curve     getRarefactionCurve True = rarefactionSampleCurve fastBin start interval-    getRarefactionCurve False =-        rarefactionCurve fastBin (fromIntegral start) (fromIntegral interval)+    getRarefactionCurve False = rarefactionCurve+        fastBin runs (fromIntegral start) (fromIntegral interval)

@@ -0,0 +1,47 @@+-- RandomSampling module.+-- By Gregory W. Schwartz+--+{- | Collection of functions pertaining to getting random samples from+- a population+-}++module Math.Diversity.RandomSampling ( subsample+                                     , subsampleSpecies+                                     , subsampleES ) where++-- Built-in+import qualified Data.Set as Set++-- Cabal+import Control.Monad.Random+import System.Random.Shuffle++-- Local+import Math.Diversity.Statistics++-- | Gets a random sampling by getting either the requested amount either+-- directly or through the inverse (for speed)+subsample :: (Eq a, Ord a) => Int -> Int -> StdGen -> [a] -> [a]+subsample size n gen xs+    | n > size `div` 2 = drop (size - n) shuffled+    | otherwise = take n shuffled+  where+    shuffled  = shuffle' xs size gen++-- | Get the number of types of entities in a subsample+subsampleSpecies :: (Eq a, Ord a) => Int -> Int -> StdGen -> [a] -> Int+subsampleSpecies size n gen = Set.size . Set.fromList . subsample size n gen++-- | 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 runs size n xs = do+    allRuns <- mapM (const run) [1..runs]++    let statTuple = medmad . map fromIntegral $ allRuns++    return statTuple+  where+    run = do+        gen <- newStdGen+        return $ subsampleSpecies size n gen xs

@@ -0,0 +1,34 @@+-- Statistics module.+-- By Gregory W. Schwartz+--+{- | Collection of functions pertaining statistics needed in the library+-}++module Math.Diversity.Statistics ( median, medmad ) where++-- Built-in+import Data.List+import Data.Maybe++-- | Get the median, adapted from+-- http://rosettacode.org/wiki/Averages/Median#Haskell, but FIXED due to+-- zero indexing+median :: [Double] -> Maybe Double+median xs | null xs  = Nothing+          | odd  len = Just $ xs !! mid+          | even len = Just meanMedian+  where+    len        = length xs+    mid        = len `div` 2+    meanMedian = (xs !! (mid - 1) + xs !! mid) / 2++-- | Get the median and median absolute deviation (MAD) of a list of numbers+medmad :: [Double] -> (Double, Double)+medmad xs = (sampleMedian, mad)+  where+    mad          = fromJust+                 . median+                 . sort+                 . map (\x -> abs (x - sampleMedian))+                 $ xs+    sampleMedian = fromJust . median . sort $ xs