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.2.0.5 +version: 0.3.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,12 +16,14 @@ cabal-version: >=1.8 library + ghc-options: -O2 hs-source-dirs: src/src-lib - exposed-modules: Diversity.Types, Diversity.Diversity, Diversity.GenerateDiversity, Diversity.Print + 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 executable diversity + ghc-options: -O2 -- Directories containing source files. hs-source-dirs: src/src-exec main-is: Main.hs 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 @@ -13,8 +13,8 @@ import Data.Fasta.String.Parse -- Local -import Diversity.GenerateDiversity -import Diversity.Print +import Math.Diversity.GenerateDiversity +import Math.Diversity.Print -- Command line arguments data Options = Options { inputLabel :: String diff --git a/src/src-lib/Diversity/Diversity.hs b/src/src-lib/Diversity/Diversity.hs deleted file mode 100644 --- a/src/src-lib/Diversity/Diversity.hs +++ /dev/null @@ -1,67 +0,0 @@ --- Diversity module. --- By G.W. Schwartz --- -{- | Collection of functions pertaining to finding the diversity of samples. --} - -module Diversity.Diversity ( hamming - , diversity - , rarefactionCurve - , rarefactionViable ) where - --- Built-in -import Data.List -import Data.Ratio -import Numeric.SpecFunctions (choose) - --- | Takes two strings, returns Hamming distance -hamming :: String -> String -> Int -hamming xs ys = length $ filter not $ zipWith (==) xs ys - --- | Returns the diversity of a list of things -diversity :: (Ord b) => Double -> [b] -> Double -diversity order sample - | length sample == 0 = 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))) - 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) -specialBinomial :: Bool -> Integer -> Integer -> Integer -> Double -specialBinomial False n_total g n = fromRational - $ product [(n_total - g - n + 1)..(n_total - g)] - % product [(n_total - n + 1)..n_total] -specialBinomial True n_total g n = choose - (fromIntegral n_total - fromIntegral g) - (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] - where - rarefact n - | n == 0 = 0 - | n == 1 = 1 - | n == n_total = k - | otherwise = k - inner n - inner n = ( \x -> if fastBin - then x / choose (fromIntegral n_total) (fromIntegral n) - else x ) - . sum - . map (\g -> specialBinomial fastBin n_total g n) - $ grouped - n_total = genericLength xs - k = genericLength grouped - grouped = map genericLength . group . sort $ xs - --- | 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/Diversity/GenerateDiversity.hs b/src/src-lib/Diversity/GenerateDiversity.hs deleted file mode 100644 --- a/src/src-lib/Diversity/GenerateDiversity.hs +++ /dev/null @@ -1,42 +0,0 @@ --- GenerateDiversity module. --- By G.W. Schwartz --- -{- | Collection of functions for the collection of fragments for the -diversity calculations. --} - -module Diversity.GenerateDiversity ( fragmentPos - , generatePositionMap ) where - --- Built in -import qualified Data.Map as M -import Data.List -import Data.Fasta.String - --- Local -import Diversity.Types - --- | 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 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) - where - combine = foldl1' (\(x, ys) (_, y) -> (x, ys ++ y)) - --- | Generate the PositionMap from a list of FastaSequences -generatePositionMap :: Bool -> Window -> [FastaSequence] -> PositionMap -generatePositionMap whole win = M.fromListWith (++) . posSeqList - where - posSeqList = map toList . concatMap (\x -> fragmentPos whole win - . map (\(p, f) -> (p, [f])) - . filter (\(_, f) -> noGaps f) - . zip [1..] - . fastaSeq - $ x) - toList (x, y) = (x, [y]) - noGaps y = y /= '-' && y /= '.' diff --git a/src/src-lib/Diversity/Print.hs b/src/src-lib/Diversity/Print.hs deleted file mode 100644 --- a/src/src-lib/Diversity/Print.hs +++ /dev/null @@ -1,73 +0,0 @@ --- Print module --- By G.W. Schwartz --- -{- | Collection of functions for the printing of data (converting data -structures into strings for use with writing to output files). --} - -module Diversity.Print ( printDiversity - , printRarefaction - , printRarefactionCurve ) where - --- Built in -import Data.List -import qualified Data.Map as M - --- Local -import Diversity.Types -import Diversity.Diversity - --- Return the results of the diversity analysis in string form for saving --- to a file -printDiversity :: Label -> Order -> Window -> PositionMap -> String -printDiversity label order window positionMap = header ++ body - where - header = "label,order,window,position,weight,diversity\n" - body = unlines - . map mapLine - . M.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 - ] - --- 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 - 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 . rarefactionCurve fastBin $ xs - ] - --- 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 - where - header = "label,window,position,weight,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 - ] diff --git a/src/src-lib/Diversity/Types.hs b/src/src-lib/Diversity/Types.hs deleted file mode 100644 --- a/src/src-lib/Diversity/Types.hs +++ /dev/null @@ -1,24 +0,0 @@ --- Types module. --- By G.W. Schwartz --- -{- | Collects all application specific types. --} - -module Diversity.Types where - -import qualified Data.Map as M - --- Basic -type Fragment = String -type Position = Int -type Diversity = Double -type Order = Double -type Label = String -type Window = Int - --- Advanced --- | At each position we have a collection of fragments to find the --- diversity of -type PositionMap = M.Map Position [Fragment] --- | At each position we have a diversity -type DiversityMap = M.Map Position Diversity diff --git a/src/src-lib/Math/Diversity/Diversity.hs b/src/src-lib/Math/Diversity/Diversity.hs new file mode 100644 --- /dev/null +++ b/src/src-lib/Math/Diversity/Diversity.hs @@ -0,0 +1,67 @@ +-- Diversity module. +-- By G.W. Schwartz +-- +{- | Collection of functions pertaining to finding the diversity of samples. +-} + +module Math.Diversity.Diversity ( hamming + , diversity + , rarefactionCurve + , rarefactionViable ) where + +-- Built-in +import Data.List +import Data.Ratio +import Numeric.SpecFunctions (choose) + +-- | Takes two strings, returns Hamming distance +hamming :: String -> String -> Int +hamming xs ys = length $ filter not $ zipWith (==) xs ys + +-- | Returns the diversity of a list of things +diversity :: (Ord b) => Double -> [b] -> Double +diversity order sample + | length sample == 0 = 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))) + 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) +specialBinomial :: Bool -> Integer -> Integer -> Integer -> Double +specialBinomial False n_total g n = fromRational + $ product [(n_total - g - n + 1)..(n_total - g)] + % product [(n_total - n + 1)..n_total] +specialBinomial True n_total g n = choose + (fromIntegral n_total - fromIntegral g) + (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] + where + rarefact n + | n == 0 = 0 + | n == 1 = 1 + | n == n_total = k + | otherwise = k - inner n + inner n = ( \x -> if fastBin + then x / choose (fromIntegral n_total) (fromIntegral n) + else x ) + . sum + . map (\g -> specialBinomial fastBin n_total g n) + $ grouped + n_total = genericLength xs + k = genericLength grouped + grouped = map genericLength . group . sort $ xs + +-- | 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 new file mode 100644 --- /dev/null +++ b/src/src-lib/Math/Diversity/GenerateDiversity.hs @@ -0,0 +1,42 @@ +-- GenerateDiversity module. +-- By G.W. Schwartz +-- +{- | Collection of functions for the collection of fragments for the +diversity calculations. +-} + +module Math.Diversity.GenerateDiversity ( fragmentPos + , generatePositionMap ) where + +-- Built in +import qualified Data.Map as M +import Data.List +import Data.Fasta.String + +-- Local +import Math.Diversity.Types + +-- | 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 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) + where + combine = foldl1' (\(x, ys) (_, y) -> (x, ys ++ y)) + +-- | Generate the PositionMap from a list of FastaSequences +generatePositionMap :: Bool -> Window -> [FastaSequence] -> PositionMap +generatePositionMap whole win = M.fromListWith (++) . posSeqList + where + posSeqList = map toList . concatMap (\x -> fragmentPos whole win + . map (\(p, f) -> (p, [f])) + . filter (\(_, f) -> noGaps f) + . zip [1..] + . fastaSeq + $ x) + toList (x, y) = (x, [y]) + noGaps y = y /= '-' && y /= '.' diff --git a/src/src-lib/Math/Diversity/Print.hs b/src/src-lib/Math/Diversity/Print.hs new file mode 100644 --- /dev/null +++ b/src/src-lib/Math/Diversity/Print.hs @@ -0,0 +1,73 @@ +-- Print module +-- By G.W. Schwartz +-- +{- | Collection of functions for the printing of data (converting data +structures into strings for use with writing to output files). +-} + +module Math.Diversity.Print ( printDiversity + , printRarefaction + , printRarefactionCurve ) where + +-- Built in +import Data.List +import qualified Data.Map as M + +-- Local +import Math.Diversity.Types +import Math.Diversity.Diversity + +-- Return the results of the diversity analysis in string form for saving +-- to a file +printDiversity :: Label -> Order -> Window -> PositionMap -> String +printDiversity label order window positionMap = header ++ body + where + header = "label,order,window,position,weight,diversity\n" + body = unlines + . map mapLine + . M.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 + ] + +-- 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 + 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 . rarefactionCurve fastBin $ xs + ] + +-- 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 + where + header = "label,window,position,weight,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 + ] diff --git a/src/src-lib/Math/Diversity/Types.hs b/src/src-lib/Math/Diversity/Types.hs new file mode 100644 --- /dev/null +++ b/src/src-lib/Math/Diversity/Types.hs @@ -0,0 +1,24 @@ +-- Types module. +-- By G.W. Schwartz +-- +{- | Collects all application specific types. +-} + +module Math.Diversity.Types where + +import qualified Data.Map as M + +-- Basic +type Fragment = String +type Position = Int +type Diversity = Double +type Order = Double +type Label = String +type Window = Int + +-- Advanced +-- | At each position we have a collection of fragments to find the +-- diversity of +type PositionMap = M.Map Position [Fragment] +-- | At each position we have a diversity +type DiversityMap = M.Map Position Diversity