diversity 0.3.0.0 → 0.3.4.0

raw patch · 6 files changed

+195/−40 lines, 6 filesdep ~fastaPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependency ranges changed: fasta

API changes (from Hackage documentation)

Files

• diversity.cabal +3/−3
• src/src-exec/Main.hs +54/−4
• src/src-lib/Math/Diversity/Diversity.hs +59/−10
• src/src-lib/Math/Diversity/GenerateDiversity.hs +27/−8
• src/src-lib/Math/Diversity/Print.hs +48/−13
• src/src-lib/Math/Diversity/Types.hs +4/−2

@@ -2,7 +2,7 @@ -- documentation, see http://haskell.org/cabal/users-guide/  name:                diversity-version:             0.3.0.0+version:             0.3.4.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   -- 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.2 && <0.6, math-functions >=0.1 && <0.2+  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  executable diversity   ghc-options: -O2@@ -28,4 +28,4 @@   hs-source-dirs:      src/src-exec   main-is:             Main.hs   -- other-modules:       -  build-depends:       diversity, base >=4.6 && <4.8, optparse-applicative >=0.10 && <0.12, fasta >=0.5.1.2 && <0.6+  build-depends:       diversity, base >=4.6 && <4.8, optparse-applicative >=0.10 && <0.12, fasta >=0.5.1.3 && <0.6

@@ -21,10 +21,14 @@                        , inputOrder             :: Double                        , inputWindow            :: Int                        , inputFasta             :: String+                       , inputSampleField       :: Int+                       , inputSubsampling       :: String                        , fastBin                :: Bool                        , removeN                :: Bool                        , wholeSeq               :: Bool                        , list                   :: Bool+                       , sample                 :: Bool+                       , rarefactionDF          :: Bool                        , outputRarefaction      :: String                        , outputRarefactionCurve :: String                        , output                 :: String@@ -58,6 +62,21 @@          <> metavar "FILE"          <> value ""          <> help "The fasta file containing the germlines and clones" )+      <*> option auto+          ( long "input-sample-field"+         <> short 'S'+         <> metavar "INT"+         <> value 1+         <> help "The index for the sample ID in the header separated by '|'\+                 \ (1 indexed)" )+      <*> strOption+          ( long "input-subsampling"+         <> short 'I'+         <> metavar "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, ..." )       <*> switch           ( long "fast-bin"          <> short 'f'@@ -77,6 +96,15 @@          <> short 'L'          <> help "Analyze a diversity of species in a list separated by lines\                  \ instead of a fasta file" )+      <*> switch+          ( long "sample"+         <> short 's'+         <> help "Whether to use sample based rarefaction (requires sample ID\+                 \ field from input-sample-field)" )+      <*> switch+          ( long "rarefaction-df"+         <> short 'd'+         <> help "Whether to output the rarefaction curve as a data frame" )       <*> strOption           ( long "output-rarefaction"          <> short 'O'@@ -108,10 +136,17 @@         window       = inputWindow opts         nFlag        = removeN opts         whole        = wholeSeq opts+        start        = read . head . words . inputSubsampling $ opts+        interval     = read . last . words . inputSubsampling $ opts          fastaListN   = parseFasta contents         fastaList    = if nFlag then removeNs fastaListN else fastaListN-        positionMap  = generatePositionMap whole window fastaList+        positionMap  = generatePositionMap+                       (sample opts)+                       (inputSampleField opts)+                       whole+                       window+                       fastaList      if (null . output $ opts)         then return ()@@ -120,12 +155,27 @@            $ positionMap     if (null . outputRarefaction $ opts)         then return ()-        else writeFile (outputRarefaction opts) $-            printRarefaction (fastBin opts) label window positionMap+        else writeFile (outputRarefaction opts)+           $ printRarefaction+             (sample opts)+             (fastBin opts)+             start+             interval+             label+             window+             positionMap     if (null . outputRarefactionCurve $ opts)         then return ()         else writeFile (outputRarefactionCurve opts) $-            printRarefactionCurve (fastBin opts) label window positionMap+            printRarefactionCurve+            (rarefactionDF opts)+            (sample opts)+            (fastBin opts)+            start+            interval+            label+            window+            positionMap  main :: IO () main = execParser opts >>= generateDiversity

@@ -4,15 +4,21 @@ {- | Collection of functions pertaining to finding the diversity of samples. -} +{-# LANGUAGE BangPatterns #-}+ module Math.Diversity.Diversity ( hamming                                 , diversity                                 , rarefactionCurve+                                , rarefactionSampleCurve                                 , rarefactionViable ) where  -- Built-in import Data.List import Data.Ratio+import qualified Data.Set as Set+import qualified Data.Foldable as F import Numeric.SpecFunctions (choose)+import Data.Function (on)  -- | Takes two strings, returns Hamming distance hamming :: String -> String -> Int@@ -21,17 +27,18 @@ -- | Returns the diversity of a list of things diversity :: (Ord b) => Double -> [b] -> Double diversity order sample-    | length sample == 0 = 0+    | null sample        = 0     | order == 1         = exp . h $ speciesList     | otherwise          = (sum . map ((** order) . p_i) $ speciesList) ** pow   where     pow          = 1 / (1 - order)-    h            = negate . sum . map (\x -> (p_i x) * (log (p_i x)))+    h            = negate . sum . map (\x -> p_i x * log (p_i x))     p_i x        = ((fromIntegral . length $ x) :: Double) /                    ((fromIntegral . length $ sample) :: Double)     speciesList  = group . sort $ sample --- | Binomial for large numbers (slow but works for big numbers)+-- | Binomial for small or large numbers (slow but works for big numbers,+-- fast but works for small numbers) specialBinomial :: Bool -> Integer -> Integer -> Integer -> Double specialBinomial False n_total g n = fromRational     $ product [(n_total - g - n + 1)..(n_total - g)]@@ -41,14 +48,20 @@                                    (fromIntegral n)  -- | Returns the rarefaction curve for each position in a list-rarefactionCurve :: (Eq a, Ord a) => Bool -> [a] -> [Double]-rarefactionCurve fastBin xs = map rarefact [1..n_total]+rarefactionCurve :: (Eq a, Ord a)+                 => Bool+                 -> Integer+                 -> Integer+                 -> [a]+                 -> [(Int, Double)]+rarefactionCurve fastBin start interval xs =+    map rarefact [start,(start + interval)..n_total]   where     rarefact n-        | n == 0       = 0-        | n == 1       = 1-        | n == n_total = k-        | otherwise    = k - inner n+        | n == 0       = (fromIntegral n, 0)+        | n == 1       = (fromIntegral n, 1)+        | n == n_total = (fromIntegral n, k)+        | otherwise    = (fromIntegral n, k - inner n)     inner n = ( \x -> if fastBin                         then x / choose (fromIntegral n_total) (fromIntegral n)                         else x )@@ -59,9 +72,45 @@     k       = genericLength grouped     grouped = map genericLength . group . sort $ xs +-- | Each sample has a collection of species, return a list of these maps+getSampleContents :: (Ord a, Ord b) => [(a, b)] -> [Set.Set b]+getSampleContents = map (Set.fromList . map snd)+                  . groupBy ((==) `on` fst)+                  . sortBy (compare `on` fst)+                  . map (\(!x, !y) -> (x, y))++-- | Returns the rarefaction curve for each position in a list+rarefactionSampleCurve :: (Ord a, Ord b)+                       => Bool+                       -> Int+                       -> Int+                       -> [(a, b)]+                       -> [(Int, Double)]+rarefactionSampleCurve fastBin start interval ls =+    map rarefact [start,(start + interval)..t_total]+  where+    rarefact t+        | t == 0       = (t, 0)+        | t == t_total = (t, richness)+        | otherwise    = (t, richness - inner t)+    inner t = ( \x -> if fastBin+                        then x / choose t_total t+                        else x )+            . F.sum+            . Set.map ( \s -> specialBinomial+                              fastBin+                              (fromIntegral t_total)+                              (numHave s samples)+                              (fromIntegral t) )+            $ speciesList+    numHave s = sum . map (\x -> if Set.member s x then 1 else 0)+    richness = fromIntegral . Set.size $ speciesList+    speciesList = Set.fromList . map snd $ ls+    t_total = genericLength samples+    samples = getSampleContents ls+ -- | 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   where     valid = dropWhile (< (0.95 * last xs)) xs-

@@ -5,6 +5,8 @@ diversity calculations. -} +{-# LANGUAGE BangPatterns #-}+ module Math.Diversity.GenerateDiversity ( fragmentPos                                         , generatePositionMap ) where @@ -12,31 +14,48 @@ import qualified Data.Map as M import Data.List import Data.Fasta.String+import qualified Data.Sequence as Seq +-- Cabal+import qualified Data.List.Split as Split+ -- Local import Math.Diversity.Types +-- | Get the sample ID of a sequence+getSample :: Int -> FastaSequence -> Sample+getSample x = (!! (x - 1)) . Split.splitOn "|" . fastaHeader+ -- | Generates fragment list from string of "win" length. This version -- differs from normal as it takes a tuple with the position as the first -- entry-fragmentPos :: Bool -> Int -> [(Position, String)] -> [(Position, String)]+fragmentPos :: Bool -> Int -> [(Position, Fragment)] -> [(Position, Fragment)] fragmentPos whole win xs | whole && null xs = error "Empty line in file!!"-                         | whole            = combine xs : []+                         | whole            = [combine xs]                          | length xs < win  = []                          | otherwise        = combine (take win xs)                                             : fragmentPos whole win (tail xs)   where-    combine = foldl1' (\(x, ys) (_, y) -> (x, ys ++ y))+    combine = foldl1' (\(!x, !ys) (_, y) -> (x, ys Seq.>< y))  -- | Generate the PositionMap from a list of FastaSequences-generatePositionMap :: Bool -> Window -> [FastaSequence] -> PositionMap-generatePositionMap whole win = M.fromListWith (++) . posSeqList+generatePositionMap :: Bool+                    -> Int+                    -> Bool+                    -> Window+                    -> [FastaSequence]+                    -> PositionMap+generatePositionMap sample sampleField whole win =+    M.fromListWith (++) . posSeqList   where-    posSeqList    = map toList . concatMap (\x -> fragmentPos whole win-                                           . map (\(p, f) -> (p, [f]))-                                           . filter (\(_, f) -> noGaps f)+    posSeqList    = map toList . concatMap (\x -> map (\(!p, !f) -> (p, sampleIt sample x f))+                                           . fragmentPos whole win+                                           . map (\(!p, !f) -> (p, Seq.singleton f))+                                           . filter (\(_, !f) -> noGaps f)                                            . zip [1..]                                            . fastaSeq                                            $ x)     toList (x, y) = (x, [y])     noGaps y = y /= '-' && y /= '.'+    sampleIt True s f = (getSample sampleField s, f)+    sampleIt False _ f = ("Sample", f)

@@ -5,6 +5,8 @@ structures into strings for use with writing to output files). -} +{-# LANGUAGE BangPatterns #-}+ module Math.Diversity.Print ( printDiversity                             , printRarefaction                             , printRarefactionCurve ) where@@ -38,8 +40,16 @@  -- Return the results of the rarefaction analysis in string form for saving -- to a file-printRarefaction :: Bool -> Label -> Window -> PositionMap -> String-printRarefaction fastBin label window positionMap = header ++ body+printRarefaction :: Bool+                 -> Bool+                 -> Int+                 -> Int+                 -> Label+                 -> Window+                 -> PositionMap+                 -> String+printRarefaction sample fastBin start interval label window positionMap =+    header ++ body   where     header           = "label,window,position,weight,percent_above\n"     body             = unlines@@ -51,23 +61,48 @@                  , show window                  , show p                  , show . length $ xs-                 , show . rarefactionViable . rarefactionCurve fastBin $ xs+                 , show+                 . rarefactionViable . map snd . getRarefactionCurve sample $ xs                  ]+    getRarefactionCurve True = rarefactionSampleCurve fastBin start interval+    getRarefactionCurve False = rarefactionCurve fastBin (fromIntegral start) (fromIntegral interval)  -- Return the results of the rarefaction analysis of the entire curve in -- string form for saving to a file-printRarefactionCurve :: Bool -> Label -> Window -> PositionMap -> String-printRarefactionCurve fastBin label window positionMap = header ++ body+printRarefactionCurve :: Bool+                      -> Bool+                      -> Bool+                      -> Int+                      -> Int+                      -> Label+                      -> Window+                      -> PositionMap+                      -> String+printRarefactionCurve asDF sample fastBin start interval label window positionMap =+    header asDF ++ body   where-    header           = "label,window,position,weight,curve\n"+    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-    mapLine (p, xs) = intercalate "," . line p $ xs-    line p xs  = [ label-                 , show window-                 , show p-                 , show . length $ xs-                 , intercalate "/" . map show . rarefactionCurve fastBin $ xs-                 ]+    mapLine (!p, !xs) = line asDF p xs+    line False p xs   = 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+    getRarefactionCurve True = rarefactionSampleCurve fastBin start interval+    getRarefactionCurve False =+        rarefactionCurve fastBin (fromIntegral start) (fromIntegral interval)

@@ -7,9 +7,11 @@ module Math.Diversity.Types where  import qualified Data.Map as M+import qualified Data.Sequence as Seq  -- Basic-type Fragment  = String+type Fragment  = Seq.Seq Char+type Sample    = String type Position  = Int type Diversity = Double type Order     = Double@@ -19,6 +21,6 @@ -- Advanced -- | At each position we have a collection of fragments to find the -- diversity of-type PositionMap     = M.Map Position [Fragment]+type PositionMap     = M.Map Position [(Sample, Fragment)] -- | At each position we have a diversity type DiversityMap    = M.Map Position Diversity