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.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 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 @@ -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 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 @@ -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 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 @@ -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) diff --git a/src/src-lib/Math/Diversity/RandomSampling.hs b/src/src-lib/Math/Diversity/RandomSampling.hs new file mode 100644 --- /dev/null +++ b/src/src-lib/Math/Diversity/RandomSampling.hs @@ -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 diff --git a/src/src-lib/Math/Diversity/Statistics.hs b/src/src-lib/Math/Diversity/Statistics.hs new file mode 100644 --- /dev/null +++ b/src/src-lib/Math/Diversity/Statistics.hs @@ -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