diversity 0.5.0.2 → 0.6.0.0
raw patch · 6 files changed
+162/−89 lines, 6 filesdep +pipesdep ~fastaPVP ok
version bump matches the API change (PVP)
Dependencies added: pipes
Dependency ranges changed: fasta
API changes (from Hackage documentation)
+ Math.Diversity.Diversity: diversityOfMap :: Double -> Map (Sample, Fragment) Int -> Double
+ Math.Diversity.Types: type FrequencyMap = Map (Sample, Fragment) Int
- Math.Diversity.Diversity: rarefactionCurve :: (Eq a, Ord a) => Bool -> Int -> Integer -> Integer -> [a] -> IO [(Int, (Double, Double))]
+ Math.Diversity.Diversity: rarefactionCurve :: Bool -> Int -> Integer -> Integer -> Map (Sample, Fragment) Int -> IO [(Int, (Double, Double))]
- Math.Diversity.Diversity: rarefactionSampleCurve :: (Ord a, Ord b) => Bool -> Int -> Int -> [(a, b)] -> IO [(Int, (Double, Double))]
+ Math.Diversity.Diversity: rarefactionSampleCurve :: Bool -> Int -> Int -> Map (Sample, Fragment) Int -> IO [(Int, (Double, Double))]
- Math.Diversity.GenerateDiversity: fragmentPos :: Bool -> Int -> [(Position, Fragment)] -> [(Position, Fragment)]
+ Math.Diversity.GenerateDiversity: fragmentPos :: Bool -> Int -> [(Position, Char)] -> [(Position, Fragment)]
- Math.Diversity.GenerateDiversity: generatePositionMap :: Bool -> Int -> Bool -> Window -> [FastaSequence] -> PositionMap
+ Math.Diversity.GenerateDiversity: generatePositionMap :: Bool -> Int -> Bool -> Window -> FastaSequence -> PositionMap
- Math.Diversity.Types: type Fragment = Seq Char
+ Math.Diversity.Types: type Fragment = String
- Math.Diversity.Types: type PositionMap = Map Position [(Sample, Fragment)]
+ Math.Diversity.Types: type PositionMap = Map Position FrequencyMap
Files
- diversity.cabal +4/−4
- src/src-exec/Main.hs +59/−20
- src/src-lib/Math/Diversity/Diversity.hs +58/−25
- src/src-lib/Math/Diversity/GenerateDiversity.hs +27/−24
- src/src-lib/Math/Diversity/Print.hs +9/−11
- src/src-lib/Math/Diversity/Types.hs +5/−5
diversity.cabal view
@@ -2,7 +2,7 @@ -- documentation, see http://haskell.org/cabal/users-guide/ name: diversity-version: 0.5.0.2+version: 0.6.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@@ -16,11 +16,11 @@ cabal-version: >=1.8 library- ghc-options: -O2+ -- ghc-options: -O2 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.5.1.5 && <0.6, 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.5.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@@ -28,4 +28,4 @@ hs-source-dirs: src/src-exec main-is: Main.hs -- other-modules: - build-depends: diversity, base >=4.6 && <4.9, optparse-applicative >=0.10 && <0.12, fasta >=0.5.1.5 && <0.6+ build-depends: diversity, base >=4.6 && <4.9, optparse-applicative >=0.10 && <0.12, fasta >=0.5.3 && <0.6, pipes >= 4.1 && < 4.2, containers >=0.5 && <0.6
src/src-exec/Main.hs view
@@ -5,14 +5,22 @@ window length (to split into fragments) by position. -} +{-# LANGUAGE BangPatterns #-}+ -- Built-in import Data.List+import Control.Monad (forever)+import qualified Data.Map.Strict as Map+import qualified System.IO as IO -- Cabal import Options.Applicative-import Data.Fasta.String.Parse+import Data.Fasta.String+import Pipes+import qualified Pipes.Prelude as P -- Local+import Math.Diversity.Types import Math.Diversity.GenerateDiversity import Math.Diversity.Print @@ -25,7 +33,7 @@ , inputSubsampling :: String , fastBin :: Bool , runs :: Int- , removeN :: Bool+ , removeNBool :: Bool , wholeSeq :: Bool , list :: Bool , sample :: Bool@@ -151,32 +159,64 @@ <> help "The csv file containing the diversities at each position.\ \ expects a string, so you need a string wven with std" ) +pipesFasta :: (MonadIO m) => IO.Handle -> Pipe String FastaSequence m ()+pipesFasta h = do+ first <- await+ getRest first ""+ where+ getRest x !acc = do+ eof <- liftIO $ IO.hIsEOF h+ if eof+ then yield FastaSequence { fastaHeader = tail x+ , fastaSeq = filter+ (`notElem` "\n\r ")+ acc }+ else do+ y <- await+ if take 1 y == ">"+ then do+ yield FastaSequence { fastaHeader = tail x+ , fastaSeq = filter+ (`notElem` "\n\r ")+ acc }+ getRest y ""+ else getRest x (acc ++ y)++pipesPositionMap :: Options -> IO PositionMap+pipesPositionMap opts = do+ h <- if null . inputFasta $ opts+ then return IO.stdin+ else IO.openFile (inputFasta opts) IO.ReadMode+ runEffect+ $ P.fold (Map.unionWith (Map.unionWith (+))) Map.empty id+ $ P.fromHandle h+ >-> toFastaSequence (list opts) h+ >-> P.map ( generatePositionMap+ (sample opts)+ (inputSampleField opts)+ (wholeSeq opts)+ (inputWindow opts)+ . removals )+ where+ toFastaSequence True _ = P.map ( \x -> FastaSequence { fastaHeader = ""+ , fastaSeq = x } )+ toFastaSequence False h = pipesFasta h+ removals = if removeNBool opts then removeN else id+ generateDiversity :: Options -> IO () generateDiversity opts = do- 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+ let contents = parseFasta label = inputLabel opts order = inputOrder opts window = inputWindow opts- nFlag = removeN opts+ nFlag = removeNBool 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- (sample opts)- (inputSampleField opts)- whole- window- fastaList- howToOutput x = if std opts then putStrLn else writeFile x + positionMap <- pipesPositionMap opts+ if (null . output $ opts) then return () else howToOutput (output opts)@@ -216,5 +256,4 @@ opts = info (helper <*> options) ( fullDesc <> progDesc "Return the diversity at each position for all sequences in a\- \ fasta file"- <> header "Diversity, Gregory W. Schwartz" )+ \ fasta file" )
src/src-lib/Math/Diversity/Diversity.hs view
@@ -8,6 +8,7 @@ module Math.Diversity.Diversity ( hamming , diversity+ , diversityOfMap , rarefactionCurve , rarefactionSampleCurve , rarefactionViable ) where@@ -15,11 +16,13 @@ -- 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) -- Local import Math.Diversity.RandomSampling+import Math.Diversity.Types -- | Takes two strings, returns Hamming distance hamming :: String -> String -> Int@@ -39,15 +42,39 @@ diversity :: (Ord b) => Double -> [b] -> Double diversity order sample | null sample = 0- | order == 1 = exp . h $ speciesList- | otherwise = (sum . map ((** order) . p_i) $ speciesList) ** pow+ | order == 1 = exp . h . speciesList $ sample+ | otherwise = ( Map.foldl' (+) 0+ . Map.map ((** order) . p_i)+ . speciesList+ $ sample )+ ** pow where pow = 1 / (1 - order)- 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+ h = negate+ . Map.foldl' (+) 0+ . Map.map (\x -> p_i x * log (p_i x))+ p_i x = (x :: Double) /+ ((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+diversityOfMap order sample+ | Map.null sample = 0+ | order == 1 = exp . h $ sample+ | otherwise = ( Map.foldl' (+) 0+ . Map.map ((** order) . p_i)+ $ sample )+ ** pow+ where+ 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)+ -- | Binomial for small or large numbers (slow but works for big numbers, -- fast but works for small numbers) specialBinomial :: Bool -> Integer -> Integer -> Integer -> Double@@ -60,12 +87,11 @@ (fromIntegral n) -- | Returns the rarefaction curve for each position in a list-rarefactionCurve :: (Eq a, Ord a)- => Bool+rarefactionCurve :: Bool -> Int -> Integer -> Integer- -> [a]+ -> 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]@@ -76,31 +102,37 @@ | 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+ statTuple <- subsampleES+ runs+ (fromIntegral n_total)+ (fromIntegral n)+ . concatMap snd+ . Map.toAscList+ . Map.mapWithKey (\(s, f) x -> replicate x f)+ $ xs return (fromIntegral n, statTuple) inner n = ( \x -> if fastBin then x / choose (fromIntegral n_total) (fromIntegral n) else x ) . sum- . map (\g -> specialBinomial fastBin n_total g n)+ . map (\g -> specialBinomial fastBin n_total (fromIntegral g) n) $ grouped- n_total = genericLength xs- k = genericLength grouped- grouped = map genericLength . group . sort $ xs+ n_total = fromIntegral . Map.foldl' (+) 0 $ xs+ k = fromIntegral . Map.size $ xs+ grouped = Map.elems 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)+getSampleContents :: Map.Map (Sample, Fragment) Int -> [Set.Set Fragment]+getSampleContents = Map.elems+ . Map.fromListWith Set.union+ . map (\(x, y) -> (x, Set.singleton y))+ . Map.keys -- | Returns the rarefaction curve for each position in a list-rarefactionSampleCurve :: (Ord a, Ord b)- => Bool+rarefactionSampleCurve :: Bool -> Int -> Int- -> [(a, b)]+ -> Map.Map (Sample, Fragment) Int -> IO [(Int, (Double, Double))] rarefactionSampleCurve !fastBin !start !interval !ls = mapM rarefact $ [start,(start + interval)..(t_total - 1)] ++ [t_total]@@ -119,13 +151,14 @@ (numHave s samples) (fromIntegral t) ) $ speciesList- numHave s = sum . map (\x -> if Set.member s x then 1 else 0)+ numHave s = genericLength . filter (Set.member s) richness = genericLength speciesList- speciesList = nub . map snd $ 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+-- | 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
src/src-lib/Math/Diversity/GenerateDiversity.hs view
@@ -11,10 +11,11 @@ , generatePositionMap ) where -- Built in-import qualified Data.Map as M+import qualified Data.Map.Strict as Map import Data.List import Data.Fasta.String import qualified Data.Sequence as Seq+import Debug.Trace -- Cabal import qualified Data.List.Split as Split@@ -28,34 +29,36 @@ -- | 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, Fragment)] -> [(Position, Fragment)]-fragmentPos whole win xs | whole && null xs = error "Empty line in file!!"- | whole = [combine xs]- | length xs < win = []- | otherwise = combine (take win xs)- : fragmentPos whole win (tail xs)+-- entry. Is in tail recursive form+fragmentPos :: Bool -> Int -> [(Position, Char)] -> [(Position, Fragment)]+fragmentPos whole win ls = fragmentPosLoop ls [] where- combine = foldl1' (\(!x, !ys) (_, y) -> (x, ys Seq.>< y))+ fragmentPosLoop [] !acc = acc+ fragmentPosLoop !xs !acc+ | whole && null xs = error "Empty line in file!!"+ | whole = [combine xs]+ | length xs < win = acc+ | otherwise = fragmentPosLoop+ (tail xs)+ (combine (take win xs) : acc)+ combine ((!x, !y):ys) = (x, y:(map snd ys)) --- | Generate the PositionMap from a list of FastaSequences+-- | Generate the frequency from a FastaSequence generatePositionMap :: Bool -> Int -> Bool -> Window- -> [FastaSequence]+ -> FastaSequence -> PositionMap-generatePositionMap sample sampleField whole win =- M.fromListWith (++) . posSeqList+generatePositionMap !sample !sampleField !whole !win = posSeqList where- 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)+ posSeqList !x = Map.fromList+ . map (\(!p, !f) -> (p, Map.singleton (sampleIt sample x f) 1))+ . fragmentPos whole win+ . filter (noGaps . snd)+ . zip [1..]+ . fastaSeq+ $ x+ noGaps y = y /= '-' && y /= '.'+ sampleIt True !s !f = (getSample sampleField s, f)+ sampleIt False _ !f = ("Sample", f)
src/src-lib/Math/Diversity/Print.hs view
@@ -13,7 +13,7 @@ -- Built in import Data.List-import qualified Data.Map as M+import qualified Data.Map.Strict as Map -- Local import Math.Diversity.Types@@ -27,15 +27,15 @@ header = "label,order,window,position,weight,diversity\n" body = unlines . map mapLine- . M.toAscList+ . Map.toAscList $ positionMap mapLine (p, xs) = intercalate "," . line p $ xs line p xs = [ label , show order , show window , show p- , show . length $ xs- , show . diversity order $ xs+ , show . Map.foldl' (+) 0 $ xs+ , show . diversityOfMap order $ xs ] -- Return the results of the rarefaction analysis in string form for saving@@ -51,7 +51,7 @@ -> IO String printRarefaction sample fastBin runs start interval label window positionMap = do- body <- fmap unlines . mapM mapLine . M.toAscList $ positionMap+ body <- fmap unlines . mapM mapLine . Map.toAscList $ positionMap return (header ++ body) where header = "label,window,position,weight,percent_above\n"@@ -61,7 +61,7 @@ return [ label , show window , show p- , show . length $ xs+ , show . Map.foldl' (+) 0 $ xs , show . rarefactionViable . map (snd . snd) $ curve ] getRarefactionCurve True = rarefactionSampleCurve fastBin start interval@@ -82,7 +82,7 @@ -> IO String printRarefactionCurve asDF sample fastBin runs start interval label window positionMap = do- body <- fmap unlines . mapM mapLine . M.toAscList $ positionMap+ body <- fmap unlines . mapM mapLine . Map.toAscList $ positionMap return (header asDF ++ body) where@@ -96,7 +96,7 @@ return . intercalate "," $ [ label , show window , show p- , show . length $ xs+ , show . Map.foldl' (+) 0 $ xs , intercalate "/" . map (show . fst . snd) $ curve@@ -111,9 +111,7 @@ -> intercalate "," [ label , show window , show p- , show- . length- $ xs+ , show . Map.foldl' (+) 0 $ xs , show x , show y , show z
src/src-lib/Math/Diversity/Types.hs view
@@ -6,11 +6,10 @@ module Math.Diversity.Types where -import qualified Data.Map as M-import qualified Data.Sequence as Seq+import qualified Data.Map.Strict as Map -- Basic-type Fragment = Seq.Seq Char+type Fragment = String type Sample = String type Position = Int type Diversity = Double@@ -19,8 +18,9 @@ type Window = Int -- Advanced+type FrequencyMap = Map.Map (Sample, Fragment) Int -- | At each position we have a collection of fragments to find the -- diversity of-type PositionMap = M.Map Position [(Sample, Fragment)]+type PositionMap = Map.Map Position FrequencyMap -- | At each position we have a diversity-type DiversityMap = M.Map Position Diversity+type DiversityMap = Map.Map Position Diversity