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