diff --git a/diversity.cabal b/diversity.cabal --- a/diversity.cabal +++ b/diversity.cabal @@ -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 diff --git a/src/src-exec/Main.hs b/src/src-exec/Main.hs --- a/src/src-exec/Main.hs +++ b/src/src-exec/Main.hs @@ -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" ) diff --git a/src/src-lib/Math/Diversity/Diversity.hs b/src/src-lib/Math/Diversity/Diversity.hs --- a/src/src-lib/Math/Diversity/Diversity.hs +++ b/src/src-lib/Math/Diversity/Diversity.hs @@ -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 diff --git a/src/src-lib/Math/Diversity/GenerateDiversity.hs b/src/src-lib/Math/Diversity/GenerateDiversity.hs --- a/src/src-lib/Math/Diversity/GenerateDiversity.hs +++ b/src/src-lib/Math/Diversity/GenerateDiversity.hs @@ -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) diff --git a/src/src-lib/Math/Diversity/Print.hs b/src/src-lib/Math/Diversity/Print.hs --- a/src/src-lib/Math/Diversity/Print.hs +++ b/src/src-lib/Math/Diversity/Print.hs @@ -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 diff --git a/src/src-lib/Math/Diversity/Types.hs b/src/src-lib/Math/Diversity/Types.hs --- a/src/src-lib/Math/Diversity/Types.hs +++ b/src/src-lib/Math/Diversity/Types.hs @@ -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