diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -1,10 +1,10 @@ # Table of Contents -1. [Changelog for ELynx](#org9916021) +1. [Changelog for ELynx](#org263dcac) - + # Changelog for ELynx diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -2,33 +2,33 @@ # The ELynx Suite -Version: 0.0.1. +Version: 0.1.0. Reproducible evolution made easy. The ELynx Suite is a Haskell library and a tool set for computational biology. The goal of the ELynx Suite is reproducible research. Evolutionary sequences and phylogenetic trees can be read, viewed, modified and simulated. Exact specification of all options is necessary, and nothing is assumed about the data -(e.g., the type of code). The command line with all arguments is consistently, -and automatically logged. The work overhead in the beginning usually pays off in -the end. +(e.g., the type of the genetic code). The command line with all arguments is +consistently, and automatically logged. The work overhead in the beginning +usually pays off in the end. The Elynx Suite consists of three library packages and two executables providing a range of sub commands. The library packages are: -- **elynx-seq:** Handle evolutionary sequences and multi sequence alignments; -- **elynx-tree:** Handle phylogenetic trees; -- **elynx-tools:** Tools for the provided executables; +- **elynx-seq:** Handle evolutionary sequences and multi sequence alignments. +- **elynx-tree:** Handle phylogenetic trees. +- **elynx-tools:** Tools for the provided executables. The executables are: -- **SLynx:** Analyze, modify, and simulate evolutionary sequences (FASTA format); -- **TLynx:** Analyze, modify, and simulate phylogenetic trees (Newick format). +- **slynx:** Analyze, modify, and simulate evolutionary sequences (FASTA format). +- **tlynx:** Analyze, modify, and simulate phylogenetic trees (Newick format). -**ELynx is still under development. We happily receive comments, ideas, feature -requests, or pull requests!** +**ELynx is actively developed. We happily receive comments, ideas, feature +requests, and pull requests!** # Installation @@ -57,8 +57,8 @@ stack install - The binaries are installed into `~/.local/bin/` which has to be added [PATH](https://en.wikipedia.org/wiki/PATH_(variable)). - Then, they can be used directly. + The binaries are installed into `~/.local/bin/` which has to be added to the + [PATH](https://en.wikipedia.org/wiki/PATH_(variable)) environment variable. Then, they can be used directly. # SLynx @@ -67,8 +67,8 @@ slynx --help - ELynx Suite version 0.5.1. Developed by Dominik Schrempf. Compiled on September - 9, 2019, at 10:48 am, UTC. + ELynx Suite version 0.1.0. Developed by Dominik Schrempf. Compiled on January + 30, 2020, at 09:18 am, UTC. Usage: slynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME] COMMAND Analyze, and simulate multi sequence alignments. @@ -82,17 +82,14 @@ Specify base name of output file Available commands: - concatenate - examine If data is a multi sequence alignment, additionally - analyze columns. - filter-rows - filter-columns - simulate - sub-sample Create a given number of multi sequence alignments, - each of which containing a given number of random - sites drawn from the original multi sequence - alignment. - translate + concatenate Concatenate sequences found in input files. + examine Examine sequences. If data is a multi sequence + alignment, additionally analyze columns. + filter-rows Filter rows (or sequences) found in input files. + filter-columns Filter columns of multi-sequence alignments. + simulate Simulate multi sequence alignments. + sub-sample Sub-sample columns from multi sequence alignments. + translate Translate from DNA to Protein or DNAX to ProteinX. File formats: - FASTA @@ -106,15 +103,20 @@ - ProteinS (amino acids; including gaps, and translation stops) - ProteinI (amino acids; including gaps, translation stops, and IUPAC codes) - The ELynx Suite. + The ELynx Suite + --------------- A Haskell library and a tool set for computational biology. The goal of the ELynx Suite is reproducible research. Evolutionary sequences and phylogenetic trees can be read, viewed, modified and simulated. Exact specification of all options is necessary, and nothing is assumed about the data (e.g., the type of code). The command line with all arguments is consistently, and automatically - logged. The work overhead in the beginning usually pays off in the end. + logged. + slynx Analyze, modify, and simulate evolutionary sequences. tlynx Analyze, modify, and simulate phylogenetic trees. + + Get help for specific commands: + slynx examine --help ## Concatenate @@ -123,9 +125,8 @@ slynx concatenate --help - Concatenate sequences found in input files. - Usage: slynx concatenate (-a|--alphabet NAME) INPUT-FILE + Concatenate sequences found in input files. Available options: -a,--alphabet NAME Specify alphabet type NAME @@ -139,10 +140,9 @@ slynx examine --help - Examine sequences. - Usage: slynx examine (-a|--alphabet NAME) [INPUT-FILE] [--per-site] - If data is a multi sequence alignment, additionally analyze columns. + Examine sequences. If data is a multi sequence alignment, additionally analyze + columns. Available options: -a,--alphabet NAME Specify alphabet type NAME @@ -157,26 +157,27 @@ slynx filter-rows --help - Filter rows (or sequences) found in input files. - Usage: slynx filter-rows (-a|--alphabet NAME) [INPUT-FILE] [--longer-than LENGTH] [--shorter-than LENGTH] + [--standard-characters] + Filter rows (or sequences) found in input files. Available options: -a,--alphabet NAME Specify alphabet type NAME INPUT-FILE Read sequences from INPUT-FILE --longer-than LENGTH Only keep sequences longer than LENGTH --shorter-than LENGTH Only keep sequences shorter than LENGTH + --standard-characters Only keep sequences containing at least one standard + (i.e., non-IUPAC) character -h,--help Show this help text Filter columns of multi sequence alignments with `filter-columns`. slynx filter-columns --help - Filter columns of multi-sequence alignments. - Usage: slynx filter-columns (-a|--alphabet NAME) [INPUT-FILE] [--standard-chars DOUBLE] + Filter columns of multi-sequence alignments. Available options: -a,--alphabet NAME Specify alphabet type NAME @@ -192,13 +193,13 @@ slynx simulate --help - Simulate multi sequence alignments. - Usage: slynx simulate (-t|--tree-file Name) [-s|--substitution-model MODEL] [-m|--mixture-model MODEL] [-e|--edm-file NAME] + [-p|--siteprofile-files NAMES] [-w|--mixture-model-weights "[DOUBLE,DOUBLE,...]"] [-g|--gamma-rate-heterogeneity "(NCAT,SHAPE)"] (-l|--length NUMBER) [-S|--seed [INT]] + Simulate multi sequence alignments. Available options: -t,--tree-file Name Read trees from file NAME @@ -210,6 +211,8 @@ are shown below (mutually exclusive with -s option) -e,--edm-file NAME Empirical distribution model file NAME in Phylobayes format + -p,--siteprofile-files NAMES + File names of site profiles in Phylobayes format -w,--mixture-model-weights "[DOUBLE,DOUBLE,...]" Weights of mixture model components -g,--gamma-rate-heterogeneity "(NCAT,SHAPE)" @@ -257,8 +260,7 @@ Usage: slynx sub-sample (-a|--alphabet NAME) [INPUT-FILE] (-n|--number-of-sites INT) (-m|--number-of-alignments INT) [-S|--seed [INT]] - Create a given number of multi sequence alignments, each of which containing a - given number of random sites drawn from the original multi sequence alignment. + Sub-sample columns from multi sequence alignments. Available options: -a,--alphabet NAME Specify alphabet type NAME @@ -269,6 +271,9 @@ -S,--seed [INT] Seed for random number generator; list of 32 bit integers with up to 256 elements (default: random) -h,--help Show this help text + + Create a given number of multi sequence alignments, each of which contains a + given number of random sites drawn from the original multi sequence alignment. ## Translate @@ -277,10 +282,9 @@ slynx translate --help - Translate from DNA to Protein or DNAX to ProteinX. - Usage: slynx translate (-a|--alphabet NAME) [INPUT-FILE] (-r|--reading-frame INT) (-u|--universal-code CODE) + Translate from DNA to Protein or DNAX to ProteinX. Available options: -a,--alphabet NAME Specify alphabet type NAME @@ -297,8 +301,8 @@ tlynx --help - ELynx Suite version 0.5.1. Developed by Dominik Schrempf. Compiled on September - 9, 2019, at 10:48 am, UTC. + ELynx Suite version 0.1.0. Developed by Dominik Schrempf. Compiled on January + 30, 2020, at 09:18 am, UTC. Usage: tlynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME] COMMAND Compare, examine, and simulate phylogenetic trees. @@ -312,26 +316,34 @@ Specify base name of output file Available commands: - compare - examine - simulate Simulate reconstructed trees using the point process. - See Gernhard, T. (2008). The conditioned - reconstructed process. Journal of Theoretical - Biology, 253(4), 769–778. - http://doi.org/10.1016/j.jtbi.2008.04.005 + distance Compute distances between many phylogenetic trees. + examine Compute summary statistics of phylogenetic trees. + simulate Simulate phylogenetic trees using birth and death + processes. + compare Compare two phylogenetic trees (compute distances and + branch-wise differences). + connect Connect two phylogenetic trees in all ways (possibly + honoring constraints). + shuffle Shuffle a phylogenetic tree (keep coalescent times, + but shuffle topology and leaves). File formats: - Newick - The ELynx Suite. + The ELynx Suite + --------------- A Haskell library and a tool set for computational biology. The goal of the ELynx Suite is reproducible research. Evolutionary sequences and phylogenetic trees can be read, viewed, modified and simulated. Exact specification of all options is necessary, and nothing is assumed about the data (e.g., the type of code). The command line with all arguments is consistently, and automatically - logged. The work overhead in the beginning usually pays off in the end. + logged. + slynx Analyze, modify, and simulate evolutionary sequences. tlynx Analyze, modify, and simulate phylogenetic trees. + + Get help for specific commands: + slynx examine --help ## Compare @@ -340,24 +352,18 @@ tlynx compare --help - Compute distances between phylogenetic trees. - - Usage: tlynx compare (-d|--distance MEASURE) [-s|--summary-statistics] - [INPUT-FILES] + Usage: tlynx compare [-n|--normalize] [-b|--bipartitions] [-i|--newick-iqtree] + NAME + Compare two phylogenetic trees (compute distances and branch-wise + differences). Available options: - -d,--distance MEASURE Type of distance to calculate (available distance - measures are listed below) - -s,--summary-statistics Report summary statistics only - INPUT-FILES Read tree(s) from INPUT-FILES; if more files are - given, one tree is expected per file + -n,--normalize Normalize trees before comparison + -b,--bipartitions Print common and missing bipartitions + -i,--newick-iqtree Use IQ-TREE Newick format (internal node labels are + branch support values) + NAME Tree file -h,--help Show this help text - - Available distance measures: - Symmetric distance: -d symmetric - Incompatible split distance: -d incompatible-split[VAL] - Collapse branches with support less than VAL before distance calculation; - in this way, only well supported difference contribute to the distance measure. ## Examine @@ -366,12 +372,13 @@ tlynx examine --help - Compute summary statistics of phylogenetic trees. - - Usage: tlynx examine [INPUT-FILE] + Usage: tlynx examine [INPUT-FILE] [-i|--newick-iqtree] + Compute summary statistics of phylogenetic trees. Available options: INPUT-FILE Read trees from INPUT-FILE + -i,--newick-iqtree Use IQ-TREE Newick format (internal node labels are + branch support values) -h,--help Show this help text @@ -381,15 +388,11 @@ tlynx simulate --help - Simulate phylogenetic trees using birth and death processes. - Usage: tlynx simulate [-t|--nTrees INT] [-n|--nLeaves INT] [-H|--height DOUBLE] [-M|--condition-on-mrca] [-l|--lambda DOUBLE] [-m|--mu DOUBLE] [-r|--rho DOUBLE] [-u|--sub-sample] [-s|--summary-statistics] [-S|--seed [INT]] - Simulate reconstructed trees using the point process. See Gernhard, T. (2008). - The conditioned reconstructed process. Journal of Theoretical Biology, 253(4), - 769–778. http://doi.org/10.1016/j.jtbi.2008.04.005 + Simulate phylogenetic trees using birth and death processes. Available options: -t,--nTrees INT Number of trees (default: 10) @@ -409,9 +412,16 @@ Height of Trees: if no tree height is given, the heights will be randomly drawn from the expected distribution given the number of leaves, the birth and the death rate. Summary statistics only: only print (NumberOfExtantChildren BranchLength) pairs for each branch of each tree. The trees are separated by a newline character. Sub-sampling: simulate one big tree with n'=round(n/rho), n'>=n, leaves, and randomly sample sub-trees with n leaves. Hence, with rho=1.0, the same tree is reported over and over again. + Gernhard, T. (2008). The conditioned reconstructed process. Journal of Theoretical Biology, 253(4), 769–778. http://doi.org/10.1016/j.jtbi.2008.04.005 # ELynx -Documentation of the library can be found on Hackage. +Documentation of the library can be found on [Hackage](https://hackage.haskell.org/): + +- [elynx-seq](https://hackage.haskell.org/package/elynx-seq) +- [elynx-tree](https://hackage.haskell.org/package/elynx-tree) +- [elynx-tools](https://hackage.haskell.org/package/elynx-tools) +- [slynx](https://hackage.haskell.org/package/slynx) +- [tlynx](https://hackage.haskell.org/package/tlynx) diff --git a/elynx-seq.cabal b/elynx-seq.cabal --- a/elynx-seq.cabal +++ b/elynx-seq.cabal @@ -1,6 +1,6 @@ cabal-version: 1.12 name: elynx-seq -version: 0.0.1 +version: 0.1.0 license: GPL-3 license-file: LICENSE copyright: Dominik Schrempf (2019) @@ -46,13 +46,14 @@ ELynx.Data.MarkovProcess.PhyloModel ELynx.Data.MarkovProcess.RateMatrix ELynx.Data.MarkovProcess.SubstitutionModel + ELynx.Data.Sequence.Alignment ELynx.Data.Sequence.Defaults - ELynx.Data.Sequence.MultiSequenceAlignment ELynx.Data.Sequence.Sequence ELynx.Data.Sequence.Translate ELynx.Export.Sequence.CountsFile ELynx.Export.Sequence.Fasta ELynx.Import.MarkovProcess.EDMModelPhylobayes + ELynx.Import.MarkovProcess.SiteprofilesPhylobayes ELynx.Import.Sequence.Fasta ELynx.Simulate.MarkovProcess ELynx.Simulate.MarkovProcessAlongTree @@ -60,19 +61,18 @@ other-modules: Paths_elynx_seq default-language: Haskell2010 - ghc-options: -Wall + ghc-options: -Wall -O2 build-depends: async >=2.2.2 && <2.3, base >=4.12.0.0 && <4.13, bytestring >=0.10.8.2 && <0.11, containers >=0.6.0.1 && <0.7, data-memocombinators >=0.5.1 && <0.6, - elynx-tools >=0.0.1 && <0.1, - elynx-tree >=0.0.1 && <0.1, + elynx-tools >=0.1.0 && <0.2, + elynx-tree >=0.1.0 && <0.2, hmatrix >=0.20.0.0 && <0.21, integration >=0.2.1 && <0.3, - lens >=4.17.1 && <4.18, - math-functions >=0.3.1.0 && <0.4, + math-functions >=0.3.2.1 && <0.4, matrices >=0.5.0 && <0.6, megaparsec >=7.0.5 && <7.1, mwc-random >=0.14.0.0 && <0.15, @@ -93,17 +93,19 @@ ELynx.Data.Alphabet.DistributionDiversitySpec ELynx.Data.MarkovProcess.AminoAcidSpec ELynx.Data.MarkovProcess.NucleotideSpec - ELynx.Data.Sequence.MultiSequenceAlignmentSpec + ELynx.Data.Sequence.AlignmentSpec ELynx.Data.Sequence.SequenceSpec ELynx.Data.Sequence.TranslateSpec ELynx.Export.Sequence.FastaSpec ELynx.Import.MarkovProcess.EDMModelPhylobayesSpec + ELynx.Import.MarkovProcess.SiteprofilesPhylobayesSpec ELynx.Import.Sequence.FastaSpec ELynx.Simulate.MarkovProcessAlongTreeSpec Files Paths_elynx_seq default-language: Haskell2010 - ghc-options: -Wall -eventlog -threaded -rtsopts -with-rtsopts=-N + ghc-options: -Wall -O2 -eventlog -threaded -rtsopts + -with-rtsopts=-N build-depends: async >=2.2.2 && <2.3, base >=4.12.0.0 && <4.13, @@ -111,14 +113,13 @@ containers >=0.6.0.1 && <0.7, data-memocombinators >=0.5.1 && <0.6, elynx-seq -any, - elynx-tools >=0.0.1 && <0.1, - elynx-tree >=0.0.1 && <0.1, + elynx-tools >=0.1.0 && <0.2, + elynx-tree >=0.1.0 && <0.2, hmatrix >=0.20.0.0 && <0.21, hspec >=2.7.1 && <2.8, hspec-megaparsec >=2.0.1 && <2.1, integration >=0.2.1 && <0.3, - lens >=4.17.1 && <4.18, - math-functions >=0.3.1.0 && <0.4, + math-functions >=0.3.2.1 && <0.4, matrices >=0.5.0 && <0.6, megaparsec >=7.0.5 && <7.1, mwc-random >=0.14.0.0 && <0.15, diff --git a/src/ELynx/Data/Alphabet/Alphabet.hs b/src/ELynx/Data/Alphabet/Alphabet.hs --- a/src/ELynx/Data/Alphabet/Alphabet.hs +++ b/src/ELynx/Data/Alphabet/Alphabet.hs @@ -34,7 +34,7 @@ Alphabet (..) , AlphabetSpec (..) , alphabetSpec - , alphabetNameVerbose + , description , isStd , isGap , isUnknown @@ -53,14 +53,14 @@ deriving (Show, Read, Eq, Ord, Enum, Bounded) -- | Verbose alphabet name. -alphabetNameVerbose :: Alphabet -> String -alphabetNameVerbose DNA = "DNA (nucleotides)" -alphabetNameVerbose DNAX = "DNAX (nucleotides; including gaps)" -alphabetNameVerbose DNAI = "DNAI (nucleotides; including gaps, and IUPAC codes)" -alphabetNameVerbose Protein = "Protein (amino acids)" -alphabetNameVerbose ProteinX = "ProteinX (amino acids; including gaps)" -alphabetNameVerbose ProteinS = "ProteinS (amino acids; including gaps, and translation stops)" -alphabetNameVerbose ProteinI = "ProteinI (amino acids; including gaps, translation stops, and IUPAC codes)" +description :: Alphabet -> String +description DNA = "DNA (nucleotides)" +description DNAX = "DNAX (nucleotides; including gaps)" +description DNAI = "DNAI (nucleotides; including gaps, and IUPAC codes)" +description Protein = "Protein (amino acids)" +description ProteinX = "ProteinX (amino acids; including gaps)" +description ProteinS = "ProteinS (amino acids; including gaps, and translation stops)" +description ProteinI = "ProteinI (amino acids; including gaps, translation stops, and IUPAC codes)" -- | Alphabet specification. 'S.Set' is used because it provides fast lookups. data AlphabetSpec = AlphabetSpec { diff --git a/src/ELynx/Data/Alphabet/DistributionDiversity.hs b/src/ELynx/Data/Alphabet/DistributionDiversity.hs --- a/src/ELynx/Data/Alphabet/DistributionDiversity.hs +++ b/src/ELynx/Data/Alphabet/DistributionDiversity.hs @@ -84,6 +84,9 @@ else fromIntegral value / fromIntegral divisor -- | For a given code vector of characters, calculate frequency of characters. +-- The input vector has arbitrary length (most often the number of sequences in +-- an alignment), the length of the output vector is the number of characters in +-- the alphabet. frequencyCharacters :: AlphabetSpec -> V.Vector Character -> V.Vector Double frequencyCharacters alph d = V.map (`saveDivision` s) counts where diff --git a/src/ELynx/Data/Character/BoundaryMutation.hs b/src/ELynx/Data/Character/BoundaryMutation.hs --- a/src/ELynx/Data/Character/BoundaryMutation.hs +++ b/src/ELynx/Data/Character/BoundaryMutation.hs @@ -12,7 +12,8 @@ * Changelog -TODO: BM states can not be read and written with a single letter, like characters. +TODO: BM states can not be read and written with a single letter, like +characters. Check status of implementation. -} @@ -34,7 +35,6 @@ , neighbors ) where -import Control.Lens import qualified Data.ByteString.Lazy.Char8 as L import Numeric.SpecFunctions (choose) @@ -140,12 +140,12 @@ | i >= stateSpaceSize n = error $ "Index " ++ show i ++ "out of bounds when population size is " ++ show n ++ "." | i < nAlleles = Bnd n (toEnum i) - | otherwise = Ply n (i' - p^._1 + 1) (p^._2) (p^._3) + | otherwise = Ply n (i' - j + 1) x k where i' = i - nAlleles l = [ (enumCombination a b * (n-1), a, b) | a <- [minBound .. pred maxBound] , b <- [succ a ..]] - p = last $ takeWhile (\e -> e^._1 <= i') l + (j, x, k) = last $ takeWhile (\(e, _, _) -> e <= i') l -- | Convert 'State' to a number 'Int' for the given population size 'PopulationSize'. -- Back conversion can be done with 'fromIndexWith', with the same population size. diff --git a/src/ELynx/Data/MarkovProcess/GammaRateHeterogeneity.hs b/src/ELynx/Data/MarkovProcess/GammaRateHeterogeneity.hs --- a/src/ELynx/Data/MarkovProcess/GammaRateHeterogeneity.hs +++ b/src/ELynx/Data/MarkovProcess/GammaRateHeterogeneity.hs @@ -22,7 +22,6 @@ , expand ) where -import Control.Lens import qualified Data.ByteString.Lazy.Char8 as L import Numeric.Integration.TanhSinh import Statistics.Distribution @@ -65,13 +64,13 @@ expandSubstitutionModel :: Int -> Double -> S.SubstitutionModel -> M.MixtureModel expandSubstitutionModel n alpha sm = M.fromSubstitutionModels name (repeat 1.0) sms where - name = sm ^. S.name <> getName n alpha + name = S.name sm <> getName n alpha sms = splitSubstitutionModel n alpha sm expandMixtureModel :: Int -> Double -> M.MixtureModel -> M.MixtureModel expandMixtureModel n alpha mm = M.concatenate name renamedMMs where - name = mm ^. M.name <> getName n alpha + name = M.name mm <> getName n alpha means = getMeans n alpha scaledMMs = map (`M.scale` mm) means names = map (("; gamma rate category " ++) . show) [1 :: Int ..] diff --git a/src/ELynx/Data/MarkovProcess/MixtureModel.hs b/src/ELynx/Data/MarkovProcess/MixtureModel.hs --- a/src/ELynx/Data/MarkovProcess/MixtureModel.hs +++ b/src/ELynx/Data/MarkovProcess/MixtureModel.hs @@ -1,5 +1,3 @@ -{-# LANGUAGE TemplateHaskell #-} - {- | Module : ELynx.Data.MarkovProcess.MixtureModel Description : Mixture models are a set of substitution models with weights @@ -21,7 +19,7 @@ Weight , Component (Component) , MixtureModel (MixtureModel) - -- * Lenses and other accessors + -- * Getters , name , getAlphabet , getWeights @@ -40,7 +38,6 @@ , summarize ) where -import Control.Lens import qualified Data.ByteString.Builder as L import qualified Data.ByteString.Lazy.Char8 as L @@ -53,37 +50,32 @@ -- | A mixture model component has a weight and a substitution model. data Component = Component - { _weight :: Weight - , _substModel :: S.SubstitutionModel + { weight :: Weight + , substModel :: S.SubstitutionModel } deriving (Show, Read) -makeLenses ''Component - -- | A mixture model with its components. data MixtureModel = MixtureModel - { _name :: S.Name - , _components :: [Component] + { name :: S.Name -- ^ Name + , components :: [Component] } deriving (Show, Read) -makeLenses ''MixtureModel - --- | Get alphabet used with mixture model. Throws error if components use different --- 'Alphabet's. +-- | Get alphabet used with mixture model. Throws error if components use +-- different 'Alphabet's. getAlphabet :: MixtureModel -> Alphabet getAlphabet mm = if isValid mm - -- then S.alphabet . substModel $ head (components mm) - then head $ mm ^.. components . traverse . substModel . S.alphabet + then S.alphabet . substModel $ head (components mm) else error "Mixture model is invalid." -- | Get weights. getWeights :: MixtureModel -> [Weight] -getWeights m = m ^.. components . traverse . weight +getWeights = map weight . components -- | Get substitution models. getSubstitutionModels :: MixtureModel -> [S.SubstitutionModel] -getSubstitutionModels m = m ^.. components . traverse . substModel +getSubstitutionModels = map substModel . components -- | Create a mixture model from a list of substitution models. @@ -93,11 +85,17 @@ -- | Concatenate mixture models. concatenate :: S.Name -> [MixtureModel] -> MixtureModel -concatenate n mms = MixtureModel n $ concatMap (view components) mms +concatenate n mms = MixtureModel n $ concatMap components mms +scaleComponent :: Double -> Component -> Component +scaleComponent s c = c {substModel = s'} + where s' = S.scale s $ substModel c + -- | Scale all substitution models of the mixture model. scale :: Double -> MixtureModel -> MixtureModel -scale s = over (components . traverse . substModel) (S.scale s) +scale s m = m {components = cs'} + where cs = components m + cs' = map (scaleComponent s) cs -- | Globally normalize a mixture model so that on average one event happens per -- unit time. @@ -107,9 +105,15 @@ weights = getWeights mm scales = map S.totalRate $ getSubstitutionModels mm +appendNameComponent :: S.Name -> Component -> Component +appendNameComponent n c = c {substModel = s'} + where s' = S.appendName n $ substModel c + -- | Append byte string to all substitution models of mixture model. appendName :: S.Name -> MixtureModel -> MixtureModel -appendName n = over (components . traverse . substModel) (S.appendName n) +appendName n m = m {components = cs'} + where cs = components m + cs' = map (appendNameComponent n) cs -- | Checks if a mixture model is valid. -- @@ -118,19 +122,24 @@ -- and providing an algebraic way of creating mixture models (empty and -- addComponent which performs necessary checks). isValid :: MixtureModel -> Bool -isValid mm = not (null $ mm ^. components) - && allEqual alphabets - where alphabets = mm ^.. components . traverse . substModel . S.alphabet +isValid m = not (null $ components m) && allEqual alphabets + where alphabets = map (S.alphabet . substModel) $ components m -- | Summarize a mixture model component; lines to be printed to screen or log. summarizeComponent :: Component -> [L.ByteString] -summarizeComponent mmc = - L.pack "Weight: " <> (L.toLazyByteString . L.doubleDec $ mmc ^. weight) - : S.summarize (mmc ^. substModel) +summarizeComponent c = + L.pack "Weight: " <> (L.toLazyByteString . L.doubleDec $ weight c) + : S.summarize (substModel c) -- | Summarize a mixture model; lines to be printed to screen or log. summarize :: MixtureModel -> [L.ByteString] -summarize mm = - L.pack ("Mixture model: " ++ mm ^. name ++ ".") - : concat [ L.pack ("Component " ++ show i ++ ":") : summarizeComponent c - | (i, c) <- zip [1 :: Int ..] (mm ^. components) ] +summarize m = + [ L.pack $ "Mixture model: " ++ name m ++ "." + , L.pack $ "Number of components: " ++ show n ++ "." ] + ++ detail + where + n = length $ components m + detail = if n <= 100 + then concat [ L.pack ("Component " ++ show i ++ ":") : summarizeComponent c + | (i, c) <- zip [1 :: Int ..] (components m) ] + else [] diff --git a/src/ELynx/Data/MarkovProcess/Nucleotide.hs b/src/ELynx/Data/MarkovProcess/Nucleotide.hs --- a/src/ELynx/Data/MarkovProcess/Nucleotide.hs +++ b/src/ELynx/Data/MarkovProcess/Nucleotide.hs @@ -34,6 +34,11 @@ -- data PhyloModel = MixtureModel | SubstitutionModel +-- +-- I think it could simply be: +-- data PhyloModel = [(Weight, SubstitutionModel)] +-- + -- data MixtureModel = [(Weight, SubstitutionModel)] -- data SubstitutionModel = SMDNA DNASubstitutionModel | SMAA AASubstitutionModel diff --git a/src/ELynx/Data/MarkovProcess/PhyloModel.hs b/src/ELynx/Data/MarkovProcess/PhyloModel.hs --- a/src/ELynx/Data/MarkovProcess/PhyloModel.hs +++ b/src/ELynx/Data/MarkovProcess/PhyloModel.hs @@ -24,7 +24,6 @@ , summarize ) where -import Control.Lens import qualified Data.ByteString.Lazy.Char8 as L import ELynx.Data.Alphabet.Alphabet @@ -39,7 +38,7 @@ -- | Extract code from phylogenetic model. getAlphabet :: PhyloModel -> Alphabet getAlphabet (MixtureModel mm) = M.getAlphabet mm -getAlphabet (SubstitutionModel sm) = sm ^. S.alphabet +getAlphabet (SubstitutionModel sm) = S.alphabet sm -- | Summarize a phylogenetic model; lines to be printed to screen or log. summarize :: PhyloModel -> [L.ByteString] diff --git a/src/ELynx/Data/MarkovProcess/SubstitutionModel.hs b/src/ELynx/Data/MarkovProcess/SubstitutionModel.hs --- a/src/ELynx/Data/MarkovProcess/SubstitutionModel.hs +++ b/src/ELynx/Data/MarkovProcess/SubstitutionModel.hs @@ -1,5 +1,3 @@ -{-# LANGUAGE TemplateHaskell #-} - {- | Module : ELynx.Data.MarkovProcess.SubstitutionModel Description : Data type describing substitution model @@ -26,7 +24,7 @@ , name , stationaryDistribution , exchangeabilityMatrix - , getRateMatrix + , rateMatrix , totalRate -- * Building substitution models , substitutionModel @@ -39,7 +37,6 @@ , summarize ) where -import Control.Lens import qualified Data.ByteString.Lazy.Char8 as L import qualified Numeric.LinearAlgebra as LinAlg @@ -56,25 +53,26 @@ type Params = [Double] -- | Complete definition of a substitution model. Create instances with --- 'substitutionModel'. +-- 'substitutionModel'. A substitution model has an alphabet, a name, and a list +-- of parameters (e.g., the kappa value for the HKY model). Further, the +-- transition rate matrix is defined by a stationary distribution and a set of +-- exchangeabilities. data SubstitutionModel = SubstitutionModel - { _alphabet :: Alphabet - , _name :: Name - , _params :: Params - , _stationaryDistribution :: R.StationaryDistribution - , _exchangeabilityMatrix :: R.ExchangeabilityMatrix + { alphabet :: Alphabet -- ^ Alphabet + , name :: Name -- ^ Name + , params :: Params -- ^ List of parameters + , stationaryDistribution :: R.StationaryDistribution -- ^ Stationary distribution + , exchangeabilityMatrix :: R.ExchangeabilityMatrix -- ^ Exchangeability matrix } deriving (Show, Read) -makeLenses ''SubstitutionModel - -- | Calculate rate matrix from substitution model. -getRateMatrix :: SubstitutionModel -> R.RateMatrix -getRateMatrix sm = R.fromExchangeabilityMatrix (sm^.exchangeabilityMatrix) (sm^.stationaryDistribution) +rateMatrix :: SubstitutionModel -> R.RateMatrix +rateMatrix sm = R.fromExchangeabilityMatrix (exchangeabilityMatrix sm) (stationaryDistribution sm) -- | Get scale of substitution model. totalRate :: SubstitutionModel -> Double -totalRate sm = R.totalRate (sm^.stationaryDistribution) (getRateMatrix sm) +totalRate sm = R.totalRate (stationaryDistribution sm) (rateMatrix sm) -- | Create normalized 'SubstitutionModel'. See 'normalize'. substitutionModel :: Alphabet -> Name -> Params @@ -90,7 +88,8 @@ -- | Scale the rate of a substitution model by given factor. scale :: Double -> SubstitutionModel -> SubstitutionModel -scale r = over exchangeabilityMatrix (LinAlg.scale r) +scale r sm = sm {exchangeabilityMatrix = em'} + where em' = LinAlg.scale r $ exchangeabilityMatrix sm -- | Normalize a substitution model, so that, on average, one substitution -- happens per unit time. @@ -100,19 +99,20 @@ -- | Abbend to name. appendName :: Name -> SubstitutionModel -> SubstitutionModel -appendName n = over name (<> n) +appendName n sm = sm {name = n'} + where n' = name sm <> n -- | Summarize a substitution model; lines to be printed to screen or log. summarize :: SubstitutionModel -> [L.ByteString] summarize sm = map L.pack $ - (show (sm^.alphabet) ++ " substitution model: " ++ sm^.name ++ ".") : - [ "Parameters: " ++ show (sm^.params) ++ "." | not (null (sm^.params))] ++ - case sm^.alphabet of - DNA -> [ "Stationary distribution: " ++ dispv precision (sm^.stationaryDistribution) ++ "." - , "Exchangeability matrix:\n" ++ dispmi 2 precision (sm^.exchangeabilityMatrix) ++ "." + (show (alphabet sm) ++ " substitution model: " ++ name sm ++ ".") : + [ "Parameters: " ++ show (params sm) ++ "." | not (null (params sm))] ++ + case alphabet sm of + DNA -> [ "Stationary distribution: " ++ dispv precision (stationaryDistribution sm) ++ "." + , "Exchangeability matrix:\n" ++ dispmi 2 precision (exchangeabilityMatrix sm) ++ "." , "Scale: " ++ show (roundN precision $ totalRate sm) ++ "." ] - Protein -> [ "Stationary distribution: " ++ dispv precision (sm^.stationaryDistribution) ++ "." + Protein -> [ "Stationary distribution: " ++ dispv precision (stationaryDistribution sm) ++ "." , "Scale: " ++ show (roundN precision $ totalRate sm) ++ "." ] _ -> error "Extended character sets are not supported with substitution models." diff --git a/src/ELynx/Data/Sequence/Alignment.hs b/src/ELynx/Data/Sequence/Alignment.hs new file mode 100644 --- /dev/null +++ b/src/ELynx/Data/Sequence/Alignment.hs @@ -0,0 +1,244 @@ +{- | +Module : ELynx.Data.Sequence.Alignment +Description : Multi sequence alignment related types and functions +Copyright : (c) Dominik Schrempf 2018 +License : GPL-3 + +Maintainer : dominik.schrempf@gmail.com +Stability : unstable + +Portability : portable + +Creation date: Thu Oct 4 18:40:18 2018. + +This module is to be imported qualified. + +-} + + +module ELynx.Data.Sequence.Alignment + ( Alignment (..) + , length + , nSequences + -- | * Input, output + , fromSequences + , toSequences + , summarize + -- | * Manipulation + , join + , concat + , concatAlignments + , filterColsOnlyStd + , filterColsStd + , filterColsNoGaps + -- | * Analysis + , FrequencyData + , distribution + , toFrequencyData + , kEffEntropy + , kEffHomoplasy + , countIUPACChars + , countGaps + , countUnknowns + -- | * Sub sample + , subSample + , randomSubSample + ) where + +import Control.Monad hiding (join) +import Control.Monad.Primitive +import qualified Data.ByteString.Lazy.Char8 as L +import Data.List hiding (concat, + length) +import qualified Data.Matrix.Unboxed as M +import qualified Data.Vector.Unboxed as V +import Prelude hiding (concat, + length) +import System.Random.MWC + +import qualified ELynx.Data.Alphabet.Alphabet as A +import ELynx.Data.Alphabet.Character +import qualified ELynx.Data.Alphabet.DistributionDiversity as D +import ELynx.Data.Sequence.Defaults +import qualified ELynx.Data.Sequence.Sequence as S +import ELynx.Tools.Concurrent +import ELynx.Tools.Definitions +import ELynx.Tools.Equality +import ELynx.Tools.Matrix + +-- | A collection of sequences. +data Alignment = Alignment + { names :: [S.Name] + , alphabet :: A.Alphabet + , matrix :: M.Matrix Character + } + deriving (Show, Eq) + +-- | Number of sites. +length :: Alignment -> Int +length = M.cols . matrix + +-- | Number of sequences. +nSequences :: Alignment -> Int +nSequences = M.rows . matrix + +-- | Create 'Alignment' from a list of 'S.Sequence's. +fromSequences :: [S.Sequence] -> Either String Alignment +fromSequences ss + | S.equalLength ss && allEqual (map S.alphabet ss) = + Right $ Alignment ns a d + | S.equalLength ss = + Left "Sequences do not have equal codes." + | otherwise = + Left "Sequences do not have equal lengths." + where + ns = map S.name ss + a = S.alphabet $ head ss + bss = map S.characters ss + d = M.fromRows bss + +-- | Conversion to list of 'S.Sequence's. +toSequences :: Alignment -> [S.Sequence] +toSequences (Alignment ns a d) = zipWith (\n r -> S.Sequence n a r) ns rows + where + rows = M.toRows d + +-- -- | Show a 'Alignment', untrimmed. +-- summarizeSeq :: Alignment -> Int -> L.ByteString +-- summarizeSeq m i = +-- L.unwords [ alignLeft nameWidth $ names m !! i +-- , summarizeByteString summaryLength $ +-- L.pack $ V.toList $ V.map toChar $ M.takeRow (matrix m) i ] + +header :: Alignment -> L.ByteString +header a = L.unlines $ + [ L.pack "Multi sequence alignment." + , L.pack $ "Code: " ++ A.description (alphabet a) ++ "." + , L.pack $ "Length: " ++ show (length a) ++ "." ] + ++ reportLengthSummary ++ reportNumberSummary + where reportLengthSummary = + [ L.pack $ "For each sequence, the " + ++ show summaryLength ++ " first bases are shown." + | length a > summaryLength ] + reportNumberSummary = + [ L.pack $ show summaryNSequences ++ " out of " ++ + show (nSequences a) ++ " sequences are shown." + | nSequences a > summaryNSequences ] + +-- | Similar to 'S.summarizeSequenceList' but with different Header. +summarize :: Alignment -> L.ByteString +summarize a = header a <> S.body (toSequences a) + +-- | Join two 'Alignment's vertically. That is, add more sequences +-- to an alignment. See also 'concat'. +join :: Alignment + -> Alignment + -> Alignment +-- top bottom. +join t b + | length t == length b && + al == alphabet b + = Alignment ns al (tD === bD) + | otherwise = error "join: Multi sequence alignments do not have equal length." + where + ns = names t ++ names b + tD = matrix t + bD = matrix b + al = alphabet t + +-- | Concatenate two 'Alignment's horizontally. That is, add more +-- sites to an alignment. See also 'join'. +concat :: Alignment + -> Alignment + -> Alignment +-- left right. +concat l r + | nSequences l == nSequences r && + al == alphabet r + = Alignment (names l) al (lD ||| rD) + | otherwise = error "concat: Multi sequence alignments do not have equal length." + where + lD = matrix l + rD = matrix r + al = alphabet l + +-- | Concatenate a list of 'Alignment's horizontally. See +-- 'concat'. +concatAlignments :: [Alignment] -> Alignment +concatAlignments [] = error "concatAlignments: Nothing to concatenate." +concatAlignments [a] = a +concatAlignments as = foldl' concat (head as) (tail as) + +-- Only keep columns from alignment that satisfy given predicate. +filterColsWith :: (V.Vector Character -> Bool) -> Alignment -> Alignment +filterColsWith p a = a {matrix = m'} + where m' = M.fromColumns . filter p . M.toColumns $ matrix a + +-- | Only keep columns with standard characters. Alignment columns with IUPAC +-- characters are removed. +filterColsOnlyStd :: Alignment -> Alignment +filterColsOnlyStd a = filterColsWith (V.all $ A.isStd (alphabet a)) a + +-- | Filter columns with proportion of standard character larger than given number. +filterColsStd :: Double -> Alignment -> Alignment +filterColsStd prop a = filterColsWith + (\col -> prop * n <= fromIntegral (V.length (V.filter (A.isStd al) col))) a + where al = alphabet a + n = fromIntegral $ nSequences a + +-- | Only keep columns without gaps or unknown characters. +filterColsNoGaps :: Alignment -> Alignment +filterColsNoGaps a = filterColsWith (V.all $ not . A.isGap (alphabet a)) a + +-- | Frequency data; do not store the actual characters, but their frequencies. +-- The matrix is of size @N x K@, where @N@ is the number of sites, and @K@ is +-- the number of characters. +type FrequencyData = M.Matrix Double + +-- | Calculcate frequency of characters at each site of a multi sequence alignment. +toFrequencyData :: Alignment -> FrequencyData +toFrequencyData a = fMapColParChunk 100 (D.frequencyCharacters spec) (matrix a) + where spec = A.alphabetSpec (alphabet a) + +-- | Calculate the distribution of characters. +distribution :: FrequencyData -> [Double] +distribution fd = map (/ fromIntegral nSites) $ V.toList $ + foldl1 (V.zipWith (+)) (M.toColumns fd) + where nSites = M.cols fd + +-- | Diversity analysis. See 'kEffEntropy'. +kEffEntropy :: FrequencyData -> [Double] +kEffEntropy fd = parMapChunk chunksize D.kEffEntropy (M.toColumns fd) + +-- | Diversity analysis. See 'kEffEntropy'. +kEffHomoplasy :: FrequencyData -> [Double] +kEffHomoplasy fd = parMapChunk chunksize D.kEffHomoplasy (M.toColumns fd) + +-- | Count the number of standard (i.e., not extended IUPAC) characters in the +-- alignment. +countIUPACChars :: Alignment -> Int +countIUPACChars a = V.length . V.filter (A.isIUPAC (alphabet a)) $ allChars + where allChars = M.flatten $ matrix a + +-- | Count the number of gaps in the alignment. +countGaps :: Alignment -> Int +countGaps a = V.length . V.filter (A.isGap (alphabet a)) $ allChars + where allChars = M.flatten $ matrix a + +-- | Count the number of unknown characters in the alignment. +countUnknowns :: Alignment -> Int +countUnknowns a = V.length . V.filter (A.isUnknown (alphabet a)) $ allChars + where allChars = M.flatten $ matrix a + +-- | Sample the given sites from a multi sequence alignment. +subSample :: [Int] -> Alignment -> Alignment +subSample is a = a {matrix = m'} + where m' = subSampleMatrix is $ matrix a + +-- | Randomly sample a given number of sites of the multi sequence alignment. +randomSubSample :: PrimMonad m + => Int -> Alignment -> Gen (PrimState m) -> m Alignment +randomSubSample n a g = do + let l = length a + is <- replicateM n $ uniformR (0, l-1) g + return $ subSample is a diff --git a/src/ELynx/Data/Sequence/Defaults.hs b/src/ELynx/Data/Sequence/Defaults.hs --- a/src/ELynx/Data/Sequence/Defaults.hs +++ b/src/ELynx/Data/Sequence/Defaults.hs @@ -14,24 +14,24 @@ module ELynx.Data.Sequence.Defaults - ( defSequenceNameWidth - , defSequenceSummaryLength - , defSequenceListSummaryNumber - , defFieldWidth + ( nameWidth + , summaryLength + , summaryNSequences + , fieldWidth ) where -- | Space reserved for sequence names when printing them. -defSequenceNameWidth :: Int -defSequenceNameWidth = 23 +nameWidth :: Int +nameWidth = 23 -- | The length shown when summarizing sequences. -defSequenceSummaryLength :: Int -defSequenceSummaryLength = 60 +summaryLength :: Int +summaryLength = 60 -- | How many sequences are shown in summary. -defSequenceListSummaryNumber :: Int -defSequenceListSummaryNumber = 200 +summaryNSequences :: Int +summaryNSequences = 200 -- | Field width for tables. -defFieldWidth :: Int -defFieldWidth = 13 +fieldWidth :: Int +fieldWidth = 13 diff --git a/src/ELynx/Data/Sequence/MultiSequenceAlignment.hs b/src/ELynx/Data/Sequence/MultiSequenceAlignment.hs deleted file mode 100644 --- a/src/ELynx/Data/Sequence/MultiSequenceAlignment.hs +++ /dev/null @@ -1,253 +0,0 @@ -{-# LANGUAGE TemplateHaskell #-} - -{- | -Module : ELynx.Data.MultiSequenceAlignment -Description : Multi sequence alignment related types and functions -Copyright : (c) Dominik Schrempf 2018 -License : GPL-3 - -Maintainer : dominik.schrempf@gmail.com -Stability : unstable - -Portability : portable - -Creation date: Thu Oct 4 18:40:18 2018. - --} - - -module ELynx.Data.Sequence.MultiSequenceAlignment - ( MultiSequenceAlignment (MultiSequenceAlignment) - , msaLength - , msaNSequences - -- | * Input, output - , fromSequenceList - , toSequenceList - , showMSA - , summarizeMSA - -- | * Manipulation - , msaJoin - , msaConcatenate - , msasConcatenate - , filterColumnsOnlyStd - , filterColumnsStd - , filterColumnsNoGaps - -- | * Analysis - , FrequencyData - , toFrequencyData - , kEffEntropy - , kEffHomoplasy - , countIUPACChars - , countGaps - , countUnknowns - -- | * Sub sample - , subSample - , randomSubSample - ) where - -import Control.Lens -import Control.Monad -import Control.Monad.Primitive -import qualified Data.ByteString.Lazy.Char8 as L -import Data.List -import qualified Data.Matrix.Unboxed as M -import qualified Data.Vector.Unboxed as V -import System.Random.MWC - -import qualified ELynx.Data.Alphabet.Alphabet as A -import ELynx.Data.Alphabet.Character -import qualified ELynx.Data.Alphabet.DistributionDiversity as D -import ELynx.Data.Sequence.Defaults -import qualified ELynx.Data.Sequence.Sequence as S -import ELynx.Tools.ByteString -import ELynx.Tools.Concurrent -import ELynx.Tools.Equality -import ELynx.Tools.Matrix - --- | A collection of sequences. -data MultiSequenceAlignment = MultiSequenceAlignment - { _names :: [S.SequenceName] - , _alphabet :: A.Alphabet - , _matrix :: M.Matrix Character - } - deriving (Read, Show, Eq) - -makeLenses ''MultiSequenceAlignment - --- | Number of sites. -msaLength :: MultiSequenceAlignment -> Int -msaLength = M.cols . view matrix - --- | Number of sequences. -msaNSequences :: MultiSequenceAlignment -> Int -msaNSequences = M.rows . view matrix - --- | Create 'MultiSequenceAlignment' from a list of 'S.Sequence's. -fromSequenceList :: [S.Sequence] -> Either String MultiSequenceAlignment -fromSequenceList ss - | S.equalLength ss && allEqual (map (view S.alphabet) ss) = - Right $ MultiSequenceAlignment ns a d - | S.equalLength ss = - Left "Sequences do not have equal codes." - | otherwise = - Left "Sequences do not have equal lengths." - where - ns = map (view S.name) ss - a = head ss ^. S.alphabet - bss = map (view S.characters) ss - d = M.fromRows bss - --- | Conversion to list of 'S.Sequence's. -toSequenceList :: MultiSequenceAlignment -> [S.Sequence] -toSequenceList (MultiSequenceAlignment ns a d) = zipWith (\n r -> S.Sequence n a r) ns rows - where - rows = M.toRows d - -msaHeader :: L.ByteString -msaHeader = L.unwords [ alignLeft defSequenceNameWidth (L.pack "Name") - , L.pack "Sequence" ] - --- | Show a 'S.Sequence', untrimmed. -showSequence :: MultiSequenceAlignment -> Int -> L.ByteString -showSequence m i = - L.unwords [ alignLeft defSequenceNameWidth $ (m ^. names) !! i - , S.fromCharacters $ M.takeRow (m ^. matrix) i ] - --- | Show a 'S.Sequence', untrimmed. -summarizeSequence :: MultiSequenceAlignment -> Int -> L.ByteString -summarizeSequence m i = - L.unwords [ alignLeft defSequenceNameWidth $ (m ^. names) !! i - , summarizeByteString defSequenceSummaryLength $ - S.fromCharacters $ M.takeRow (m ^. matrix) i ] - --- | Show a 'MultiSequenceAlignment' in text form. -showMSA :: MultiSequenceAlignment -> L.ByteString -showMSA msa = L.unlines $ msaHeader : - map (showSequence msa) [0 .. (msaNSequences msa - 1)] - -summarizeMSAHeader :: MultiSequenceAlignment -> L.ByteString -summarizeMSAHeader msa = L.unlines $ - [ L.pack "Multi sequence alignment." - , L.pack $ "Code: " ++ A.alphabetNameVerbose (msa^.alphabet) ++ "." - , L.pack $ "Length: " ++ show (msaLength msa) ++ "." ] - ++ reportLengthSummary ++ reportNumberSummary - where reportLengthSummary = - [ L.pack $ "For each sequence, the " - ++ show defSequenceSummaryLength ++ " first bases are shown." - | msaLength msa > defSequenceSummaryLength ] - reportNumberSummary = - [ L.pack $ show defSequenceListSummaryNumber ++ " out of " ++ - show (msaNSequences msa) ++ " sequences are shown." - | msaNSequences msa > defSequenceListSummaryNumber ] - --- | Similar to 'S.summarizeSequenceList' but with different Header. -summarizeMSA :: MultiSequenceAlignment -> L.ByteString -summarizeMSA msa = L.unlines $ summarizeMSAHeader msa : - map (summarizeSequence msa) [0 .. n - 1] - where n = min (msaNSequences msa) defSequenceListSummaryNumber - --- | Join two 'MultiSequenceAlignment's vertically. That is, add more sequences --- to an alignment. See also 'msaConcatenate'. -msaJoin :: MultiSequenceAlignment - -> MultiSequenceAlignment - -> MultiSequenceAlignment --- top bottom. -msaJoin t b - | msaLength t == msaLength b - && t^.alphabet == b^.alphabet - = MultiSequenceAlignment ns a (tD === bD) - | otherwise = error "msaJoin: Multi sequence alignments do not have equal length." - where - ns = t^.names ++ b^.names - tD = t^.matrix - bD = t^.matrix - a = t^.alphabet - --- | Concatenate two 'MultiSequenceAlignment's horizontally. That is, add more --- sites to an alignment. See also 'msaJoin'. -msaConcatenate :: MultiSequenceAlignment - -> MultiSequenceAlignment - -> MultiSequenceAlignment --- left right. -msaConcatenate l r - | msaNSequences l == msaNSequences r - && l^.alphabet == r^.alphabet - = MultiSequenceAlignment (l ^. names) a (lD ||| rD) - | otherwise = error "msaConcatenate: Multi sequence alignments do not have equal length." - where - lD = l^.matrix - rD = r^.matrix - a = l^.alphabet - --- | Concatenate a list of 'MultiSequenceAlignment's horizontally. See --- 'msaConcatenate'. -msasConcatenate :: [MultiSequenceAlignment] -> MultiSequenceAlignment -msasConcatenate [] = error "msasConcatenate: Nothing to concatenate." -msasConcatenate [msa] = msa -msasConcatenate msas = foldl' msaConcatenate (head msas) (tail msas) - --- Only keep columns from alignment that satisfy given predicate. -filterColumnsWith :: (V.Vector Character -> Bool) -> MultiSequenceAlignment -> MultiSequenceAlignment -filterColumnsWith p = over matrix (M.fromColumns . filter p . M.toColumns) - --- | Only keep columns with standard characters. Alignment columns with IUPAC --- characters are removed. -filterColumnsOnlyStd :: MultiSequenceAlignment -> MultiSequenceAlignment -filterColumnsOnlyStd msa = filterColumnsWith (V.all $ A.isStd (msa^.alphabet)) msa - --- | Filter columns with proportion of standard character larger than given number. -filterColumnsStd :: Double -> MultiSequenceAlignment -> MultiSequenceAlignment -filterColumnsStd prop msa = filterColumnsWith - (\col -> prop * nSeqs <= fromIntegral (V.length (V.filter (A.isStd a) col))) - msa - where a = msa^.alphabet - nSeqs = fromIntegral $ msaNSequences msa - --- | Only keep columns without gaps or unknown characters. -filterColumnsNoGaps :: MultiSequenceAlignment -> MultiSequenceAlignment -filterColumnsNoGaps msa = filterColumnsWith (V.all $ not . A.isGap (msa^.alphabet)) msa - --- | Frequency data; do not store the actual characters, but only their --- frequencies. -type FrequencyData = M.Matrix Double - --- | Calculcate frequency of characters in multi sequence alignment. -toFrequencyData :: MultiSequenceAlignment -> FrequencyData -toFrequencyData msa = fMapColParChunk 100 (D.frequencyCharacters spec) (msa^.matrix) - where spec = A.alphabetSpec (msa^.alphabet) - --- | Diversity analysis. See 'kEffEntropy'. -kEffEntropy :: FrequencyData -> [Double] -kEffEntropy fd = parMapChunk 500 D.kEffEntropy (M.toColumns fd) - --- | Diversity analysis. See 'kEffEntropy'. -kEffHomoplasy :: FrequencyData -> [Double] -kEffHomoplasy fd = parMapChunk 500 D.kEffHomoplasy (M.toColumns fd) - --- | Count the number of standard (i.e., not extended IUPAC) characters in the --- alignment. -countIUPACChars :: MultiSequenceAlignment -> Int -countIUPACChars msa = V.length . V.filter (A.isIUPAC (msa^.alphabet)) $ allChars - where allChars = M.flatten $ msa^.matrix - --- | Count the number of gaps in the alignment. -countGaps :: MultiSequenceAlignment -> Int -countGaps msa = V.length . V.filter (A.isGap (msa^.alphabet)) $ allChars - where allChars = M.flatten $ msa^.matrix - --- | Count the number of unknown characters in the alignment. -countUnknowns :: MultiSequenceAlignment -> Int -countUnknowns msa = V.length . V.filter (A.isUnknown (msa^.alphabet)) $ allChars - where allChars = M.flatten $ msa^.matrix - --- | Sample the given sites from a multi sequence alignment. -subSample :: [Int] -> MultiSequenceAlignment -> MultiSequenceAlignment -subSample is = over matrix (subSampleMatrix is) - --- | Randomly sample a given number of sites of the multi sequence alignment. -randomSubSample :: PrimMonad m - => Int -> MultiSequenceAlignment -> Gen (PrimState m) -> m MultiSequenceAlignment -randomSubSample n msa g = do - let l = msaLength msa - is <- replicateM n $ uniformR (0, l-1) g - return $ subSample is msa diff --git a/src/ELynx/Data/Sequence/Sequence.hs b/src/ELynx/Data/Sequence/Sequence.hs --- a/src/ELynx/Data/Sequence/Sequence.hs +++ b/src/ELynx/Data/Sequence/Sequence.hs @@ -1,4 +1,4 @@ -{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE OverloadedStrings #-} {- | Module : ELynx.Data.Sequence @@ -12,46 +12,44 @@ Creation date: Thu Oct 4 18:54:51 2018. +This module is to be imported qualified. + -} module ELynx.Data.Sequence.Sequence ( -- * Types - SequenceName - , SequenceCharacters - , Sequence (Sequence) - -- * Lenses - , name - , alphabet - , characters + Name + , Characters + , Sequence (..) -- * Input - , toCharacters + , fromByteString -- * Output - , fromCharacters - , showSequence - , showSequenceList - , sequenceListHeader - , summarizeSequence - , summarizeSequenceList - , summarizeSequenceListBody + , toByteString + , header + , summarize + , summarizeSequences + , body -- * Analysis - , lengthSequence + , length , equalLength , longest -- * Manipulation - , trimSequence - , concatenate - , concatenateSeqs + , trim + , concat + , concatSequences -- * Filtering , filterShorterThan , filterLongerThan + , filterStandard ) where -import Control.Lens import Control.Parallel.Strategies import qualified Data.ByteString.Lazy.Char8 as L import Data.List (maximumBy) import Data.Ord (comparing) import qualified Data.Vector.Unboxed as V +import Prelude hiding (concat, length) +import qualified Prelude as Pr (length) import qualified Text.Printf as P import qualified ELynx.Data.Alphabet.Alphabet as A @@ -60,122 +58,116 @@ import ELynx.Tools.ByteString import ELynx.Tools.Equality --- | For now, 'SequenceName's are just 'L.ByteString's. -type SequenceName = L.ByteString +-- | For now, 'Name's are just 'L.ByteString's. +type Name = L.ByteString -- | The vector of characters of a sequence. -type SequenceCharacters = V.Vector Character - --- | Sequences have a name, a code and hopefully a lot of data. -data Sequence = Sequence { _name :: SequenceName - , _alphabet :: A.Alphabet - , _characters :: SequenceCharacters } - deriving (Eq) - -makeLenses ''Sequence +type Characters = V.Vector Character -- | Convert byte string to sequence characters. -toCharacters :: L.ByteString -> SequenceCharacters -toCharacters = V.fromList . map fromChar . L.unpack +fromByteString :: L.ByteString -> Characters +fromByteString = V.fromList . map fromChar . L.unpack -- | Convert sequence characters to byte string. -fromCharacters :: SequenceCharacters -> L.ByteString -fromCharacters = L.pack . map toChar . V.toList +toByteString :: Characters -> L.ByteString +toByteString = L.pack . map toChar . V.toList -showInfo :: Sequence -> L.ByteString -showInfo s = L.unwords [ alignLeft defSequenceNameWidth (s^.name) - , alignRight defFieldWidth (L.pack $ show $ s^.alphabet) - , alignRight defFieldWidth (L.pack . show $ len) - , alignRight defFieldWidth (L.pack $ P.printf "%.3f" pGaps) ] - where len = lengthSequence s +-- | Sequences have a name, a code and hopefully a lot of data. +data Sequence = Sequence { name :: Name + , alphabet :: A.Alphabet + , characters :: Characters } + deriving (Show, Eq) + +getInfo :: Sequence -> L.ByteString +getInfo s = L.unwords [ alignLeft nameWidth (name s) + , alignRight fieldWidth (L.pack $ show $ alphabet s) + , alignRight fieldWidth (L.pack . show $ len) + , alignRight fieldWidth (L.pack $ P.printf "%.3f" pGaps) ] + where len = length s nGaps = countGaps s pGaps = fromIntegral nGaps / fromIntegral len :: Double -instance Show Sequence where - show s = L.unpack $ showSequence s - --- | Show a 'Sequence', untrimmed. -showSequence :: Sequence -> L.ByteString -showSequence s = L.unwords [showInfo s, fromCharacters $ s^.characters] - --- | Show a list of 'Sequence's, untrimmed. -showSequenceList :: [Sequence] -> L.ByteString -showSequenceList = L.unlines . map showSequence - --- | Header printed before 'Sequence' list. -sequenceListHeader :: L.ByteString -sequenceListHeader = L.unwords [ alignLeft defSequenceNameWidth (L.pack "Name") - , alignRight defFieldWidth (L.pack "Code") - , alignRight defFieldWidth (L.pack "Length") - , alignRight defFieldWidth (L.pack "Gaps [%]") - , L.pack "Sequence" ] - -- | Trim and show a 'Sequence'. -summarizeSequence :: Sequence -> L.ByteString -summarizeSequence s = L.unwords [ showInfo s - , summarizeByteString defSequenceSummaryLength - (fromCharacters $ s^.characters) ] +summarize :: Sequence -> L.ByteString +summarize s = L.unwords [ getInfo s , summarizeByteString summaryLength $ toByteString (characters s) ] -- | Trim and show a list of 'Sequence's. -summarizeSequenceList :: [Sequence] -> L.ByteString -summarizeSequenceList ss = summarizeSequenceListHeader ss <> - summarizeSequenceListBody (take defSequenceListSummaryNumber ss) +summarizeSequences :: [Sequence] -> L.ByteString +summarizeSequences ss = header ss <> body (take summaryNSequences ss) -summarizeSequenceListHeader :: [Sequence] -> L.ByteString -summarizeSequenceListHeader ss = L.unlines $ +-- | Header printed before 'Sequence' list. +tableHeader :: L.ByteString +tableHeader = L.unwords [ alignLeft nameWidth "Name" + , alignRight fieldWidth "Code" + , alignRight fieldWidth "Length" + , alignRight fieldWidth "Gaps [%]" + , "Sequence" ] + +-- | A short description of the sequence. +header :: [Sequence] -> L.ByteString +header ss = L.unlines $ reportIfSubsetIsShown ++ - [ L.pack $ "For each sequence, the " ++ show defSequenceSummaryLength ++ " first bases are shown." - , L.pack $ "List contains " ++ show (length ss) ++ " sequences." - , L.pack "" - , sequenceListHeader ] - where l = length ss - s = show defSequenceListSummaryNumber ++ " out of " ++ - show (length ss) ++ " sequences are shown." + [ L.pack $ "For each sequence, the " ++ show summaryLength ++ " first bases are shown." + , L.pack $ "List contains " ++ show (Pr.length ss) ++ " sequences." + , "" + , tableHeader ] + where l = Pr.length ss + s = show summaryNSequences ++ " out of " ++ + show (Pr.length ss) ++ " sequences are shown." reportIfSubsetIsShown - | l > defSequenceListSummaryNumber = [L.pack s] + | l > summaryNSequences = [L.pack s] | otherwise = [] -- | Trim and show a list of 'Sequence's. -summarizeSequenceListBody :: [Sequence] -> L.ByteString -summarizeSequenceListBody ss = L.unlines (map summarizeSequence ss `using` parListChunk 5 rdeepseq) +body :: [Sequence] -> L.ByteString +body ss = L.unlines (map summarize ss `using` parListChunk 5 rdeepseq) -- | Calculate length of 'Sequence'. -lengthSequence :: Sequence -> Int -lengthSequence s = fromIntegral $ V.length $ s ^. characters +length :: Sequence -> Int +length = fromIntegral . V.length . characters -- | Check if all 'Sequence's have equal length. equalLength :: [Sequence] -> Bool -equalLength = allEqual . map lengthSequence +equalLength = allEqual . map length -- | Find the longest 'Sequence' in a list. longest :: [Sequence] -> Sequence -longest = maximumBy (comparing lengthSequence) +longest = maximumBy (comparing length) -- | Count number of gaps or unknown characters in sequence. countGaps :: Sequence -> Int -countGaps s = V.length . V.filter (A.isGap (s^.alphabet)) $ s^.characters +countGaps s = V.length . V.filter (A.isGap $ alphabet s) $ characters s -- | Trim to given length. -trimSequence :: Int -> Sequence -> Sequence -trimSequence n = over characters (V.take $ fromIntegral n) +trim :: Int -> Sequence -> Sequence +trim n (Sequence nm a cs) = Sequence nm a (V.take (fromIntegral n) cs) --- | Concatenate two sequences. 'SequenceName's have to match. -concatenate :: Sequence -> Sequence -> Sequence -concatenate (Sequence i c cs) (Sequence j k ks) +-- | Concatenate two sequences. 'Name's have to match. +concat :: Sequence -> Sequence -> Sequence +concat (Sequence i c cs) (Sequence j k ks) | i == j && c == k = Sequence i c (cs <> ks) | otherwise = error $ "concatenate: Sequences do not have equal names: " ++ L.unpack i ++ ", " ++ L.unpack j ++ "." --- | Concatenate a list of sequences, see 'concatenate'. -concatenateSeqs :: [[Sequence]] -> [Sequence] -concatenateSeqs [] = error "concatenateSeqs: Nothing to concatenate." -concatenateSeqs [ss] = ss -concatenateSeqs sss = foldl1 (zipWith concatenate) sss +-- | Concatenate a list of sequences, see 'concat'. +concatSequences :: [[Sequence]] -> [Sequence] +concatSequences [] = error "concatenateSequences: Nothing to concatenate." +concatSequences [ss] = ss +concatSequences sss = foldl1 (zipWith concat) sss -- | Only take 'Sequence's that are shorter than a given number. filterShorterThan :: Int -> [Sequence] -> [Sequence] -filterShorterThan n = filter (\x -> lengthSequence x < n) +filterShorterThan n = filter (\x -> length x < n) -- | Only take 'Sequence's that are longer than a given number. filterLongerThan :: Int -> [Sequence] -> [Sequence] -filterLongerThan n = filter (\x -> lengthSequence x > n) +filterLongerThan n = filter (\x -> length x > n) + +-- | Only take 'Sequence's that contain at least on non-IUPAC character. +filterStandard :: [Sequence] -> [Sequence] +filterStandard = filter (\s -> anyStandard (alphabet s) s) + +-- Are all characters IUPAC characters? +anyStandard :: A.Alphabet -> Sequence -> Bool +anyStandard a s = V.any (A.isStd a) cs + where cs = characters s diff --git a/src/ELynx/Data/Sequence/Translate.hs b/src/ELynx/Data/Sequence/Translate.hs --- a/src/ELynx/Data/Sequence/Translate.hs +++ b/src/ELynx/Data/Sequence/Translate.hs @@ -37,8 +37,8 @@ -- Translate from DNA to Protein with given reading frame (0, 1, 2). translateVecWith :: (V.Unbox a, Ord a, V.Unbox b) => (Codon a -> b) -> Int -> V.Vector a -> V.Vector b -translateVecWith f rf cs | rf > 2 = error "translate: reading frame is larger than 2." - | rf < 0 = error "translate: reading frame is negative." +translateVecWith f rf cs | rf > 2 = error "translateVecWith: reading frame is larger than 2." + | rf < 0 = error "translateVecWith: reading frame is negative." | otherwise = aas where codons = map unsafeFromVec $ chop 3 $ V.drop rf cs aas = V.fromList $ map f codons diff --git a/src/ELynx/Export/Sequence/Fasta.hs b/src/ELynx/Export/Sequence/Fasta.hs --- a/src/ELynx/Export/Sequence/Fasta.hs +++ b/src/ELynx/Export/Sequence/Fasta.hs @@ -21,7 +21,6 @@ , sequencesToFasta ) where -import Control.Lens import qualified Data.ByteString.Lazy.Char8 as L import ELynx.Data.Sequence.Sequence @@ -31,7 +30,7 @@ -- | Convert a 'Sequence' to Fasta format. sequenceToFasta :: Sequence -> L.ByteString -sequenceToFasta s = L.unlines [ fastaHeader $ s^.name , fromCharacters $ s^.characters ] +sequenceToFasta s = L.unlines [ fastaHeader $ name s, toByteString $ characters s ] -- | Convert a list 'Sequence's to Fasta format. A newline is added between any -- two 'Sequence's. diff --git a/src/ELynx/Import/MarkovProcess/SiteprofilesPhylobayes.hs b/src/ELynx/Import/MarkovProcess/SiteprofilesPhylobayes.hs new file mode 100644 --- /dev/null +++ b/src/ELynx/Import/MarkovProcess/SiteprofilesPhylobayes.hs @@ -0,0 +1,77 @@ +{- | +Module : ELynx.Import.MarkovProcess.SiteprofilesPhylobayes +Description : Import site profiles in Phylobayes format +Copyright : (c) Dominik Schrempf 2019 +License : GPL-3 + +Maintainer : dominik.schrempf@gmail.com +Stability : unstable +Portability : portable + +Creation date: Tue Jan 29 12:12:55 2019. + +For now I just try to go with a huge empirical distribution mixture model. Let's +see if performance is good enough. + +-} + +module ELynx.Import.MarkovProcess.SiteprofilesPhylobayes + ( Parser + , EDMComponent + , siteprofiles + ) where + +import Control.Monad +import qualified Data.ByteString.Lazy.Char8 as L +import Data.List (nub) +import qualified Data.Vector.Storable as V +import Data.Void +import Text.Megaparsec +import Text.Megaparsec.Byte +import Text.Megaparsec.Byte.Lexer + +import ELynx.Data.MarkovProcess.EDMModel +import ELynx.Tools.ByteString (c2w) + +-- | Shortcut. +type Parser = Parsec Void L.ByteString + +-- | Parse stationary distributions from Phylobayes format. +siteprofiles :: Parser [EDMComponent] +siteprofiles = do + _ <- headerLines + cs <- many dataLine + _ <- many newline *> eof + "phylobayes siteprofiles" + let ls = map length cs + nLs = length $ nub ls + when (nLs /= 1) (error "The site profiles have a different number of entries.") + return cs + +horizontalSpace :: Parser () +horizontalSpace = skipMany $ char (c2w ' ') <|> tab + +line :: Parser () +line = do + _ <- many $ noneOf [c2w '\n'] + pure () + +-- For now, just ignore the header. +headerLines :: Parser () +headerLines = do + _ <- line + _ <- many newline + "headerLine" + pure () + +dataLine :: Parser EDMComponent +dataLine = do + -- Ignore site number. + _ <- decimal :: Parser Integer + _ <- horizontalSpace + -- Also ignore additional white space on line. + vals <- float `sepEndBy1` horizontalSpace + _ <- many newline + "dataLine" + -- Set the weight to 1.0 for all sites. + return (1.0, V.fromList vals) diff --git a/src/ELynx/Import/Sequence/Fasta.hs b/src/ELynx/Import/Sequence/Fasta.hs --- a/src/ELynx/Import/Sequence/Fasta.hs +++ b/src/ELynx/Import/Sequence/Fasta.hs @@ -20,7 +20,8 @@ module ELynx.Import.Sequence.Fasta - ( fastaSequence + ( Parser + , fastaSequence , fasta ) where @@ -72,7 +73,7 @@ let !alph = S.map toWord (A.all . alphabetSpec $ a) lns <- some (sequenceLine alph) _ <- many eol - return $ Sequence hd a (toCharacters $ L.concat lns) + return $ Sequence hd a (fromByteString $ L.concat lns) -- | Parse a Fasta file with given 'Alphabet'. fasta :: Alphabet -> Parser [Sequence] diff --git a/src/ELynx/Simulate/MarkovProcess.hs b/src/ELynx/Simulate/MarkovProcess.hs --- a/src/ELynx/Simulate/MarkovProcess.hs +++ b/src/ELynx/Simulate/MarkovProcess.hs @@ -37,8 +37,8 @@ probMatrix :: RateMatrix -> Double -> ProbMatrix probMatrix q t | t == 0 = if rows q == cols q then ident (rows q) - else error "Matrix is not square." - | t < 0 = error "Time is negative." + else error "probMatrix: Matrix is not square." + | t < 0 = error "probMatrix: Time is negative." | otherwise = expm $ scale t q -- | Move from a given state to a new one according to a transition probability @@ -48,7 +48,6 @@ -- This function is the bottleneck of the simulator and takes up most of the -- computation time. However, I was not able to find a faster implementation -- than the one from Data.Distribution. --- jump :: (PrimMonad m) => State -> ProbMatrix -> Gen (PrimState m) -> m State jump i p = categorical (p ! i) diff --git a/src/ELynx/Simulate/MarkovProcessAlongTree.hs b/src/ELynx/Simulate/MarkovProcessAlongTree.hs --- a/src/ELynx/Simulate/MarkovProcessAlongTree.hs +++ b/src/ELynx/Simulate/MarkovProcessAlongTree.hs @@ -1,5 +1,5 @@ {- | - Description : Functions to work with transition probability matrices on rooted trees + Description : Work with transition probability matrices on rooted trees Copyright : (c) Dominik Schrempf 2017 License : GPLv3 @@ -17,27 +17,30 @@ module ELynx.Simulate.MarkovProcessAlongTree ( -- * Single rate matrix. - simulateNSitesAlongTree - , simulateAndFlattenNSitesAlongTree + simulate + , simulateAndFlatten -- * Mixture models. - , simulateNSitesAlongTreeMixtureModel - , simulateAndFlattenNSitesAlongTreeMixtureModel + , simulateMixtureModel + , simulateAndFlattenMixtureModel + , simulateAndFlattenMixtureModelPar ) where import Control.Monad import Control.Monad.Primitive +import Control.Parallel.Strategies import Data.Tree import Numeric.LinearAlgebra -import System.Random.MWC -import System.Random.MWC.Distributions +import System.Random.MWC (Gen, GenIO) +import System.Random.MWC.Distributions (categorical) import ELynx.Data.MarkovProcess.RateMatrix import ELynx.Data.Tree.MeasurableTree import ELynx.Simulate.MarkovProcess +import ELynx.Tools.Concurrent -measureableTreeToProbTree :: (Measurable a) => RateMatrix -> Tree a -> Tree ProbMatrix -measureableTreeToProbTree q = fmap (probMatrix q . getLen) +toProbTree :: (Measurable a) => RateMatrix -> Tree a -> Tree ProbMatrix +toProbTree q = fmap (probMatrix q . getLen) getRootStates :: PrimMonad m => Int -> StationaryDistribution -> Gen (PrimState m) -> m [State] @@ -46,48 +49,48 @@ -- | Simulate a number of sites for a given substitution model. Only the states -- at the leafs are retained. The states at internal nodes are removed. This has -- a lower memory footprint. -simulateAndFlattenNSitesAlongTree :: (PrimMonad m, Measurable a) +simulateAndFlatten :: (PrimMonad m, Measurable a) => Int -> StationaryDistribution -> ExchangeabilityMatrix -> Tree a -> Gen (PrimState m) -> m [[State]] -simulateAndFlattenNSitesAlongTree n d e t g = do +simulateAndFlatten n d e t g = do let q = fromExchangeabilityMatrix e d - pt = measureableTreeToProbTree q t + pt = toProbTree q t is <- getRootStates n d g - simulateAndFlattenAlongProbTree is pt g + simulateAndFlatten' is pt g -- This is the heart of the simulation. Take a tree and a list of root states. -- Recursively jump down the branches to the leafs. Forget states at internal. -simulateAndFlattenAlongProbTree :: (PrimMonad m) +simulateAndFlatten' :: (PrimMonad m) => [State] -> Tree ProbMatrix -> Gen (PrimState m) -> m [[State]] -simulateAndFlattenAlongProbTree is (Node p f) g = do +simulateAndFlatten' is (Node p f) g = do is' <- mapM (\i -> jump i p g) is if null f then return [is'] - else concat <$> sequence [simulateAndFlattenAlongProbTree is' t g | t <- f] + else concat <$> sequence [simulateAndFlatten' is' t g | t <- f] -- | Simulate a number of sites for a given substitution model. Keep states at -- internal nodes. The result is a tree with the list of simulated states as -- node labels. -simulateNSitesAlongTree :: (PrimMonad m, Measurable a) +simulate :: (PrimMonad m, Measurable a) => Int -> StationaryDistribution -> ExchangeabilityMatrix -> Tree a -> Gen (PrimState m) -> m (Tree [State]) -simulateNSitesAlongTree n d e t g = do +simulate n d e t g = do let q = fromExchangeabilityMatrix e d - pt = measureableTreeToProbTree q t + pt = toProbTree q t is <- getRootStates n d g - simulateAlongProbTree is pt g + simulate' is pt g -- This is the heart of the simulation. Take a tree and a list of root states. -- Recursively jump down the branches to the leafs. -simulateAlongProbTree :: (PrimMonad m) +simulate' :: (PrimMonad m) => [State] -> Tree ProbMatrix -> Gen (PrimState m) -> m (Tree [State]) -simulateAlongProbTree is (Node p f) g = do +simulate' is (Node p f) g = do is' <- mapM (\i -> jump i p g) is - f' <- sequence [simulateAlongProbTree is' t g | t <- f] + f' <- sequence [simulate' is' t g | t <- f] return $ Node is' f' -measureableTreeToProbTreeMixtureModel :: (Measurable a) +toProbTreeMixtureModel :: (Measurable a) => [RateMatrix] -> Tree a -> Tree [ProbMatrix] -measureableTreeToProbTreeMixtureModel qs = - fmap (\a -> [probMatrix q . getLen $ a | q <- qs]) +toProbTreeMixtureModel qs = + fmap (\a -> [probMatrix q . getLen $ a | q <- qs] `using` parList rpar) getComponentsAndRootStates :: PrimMonad m => Int -> Vector R -> [StationaryDistribution] -> Gen (PrimState m) -> m ([Int], [State]) @@ -98,41 +101,54 @@ -- | Simulate a number of sites for a given set of substitution models with -- corresponding weights. Forget states at internal nodes. See also --- 'simulateAndFlattenNSitesAlongTree'. -simulateAndFlattenNSitesAlongTreeMixtureModel :: (PrimMonad m, Measurable a) +-- 'simulateAndFlatten'. +simulateAndFlattenMixtureModel :: (PrimMonad m, Measurable a) => Int -> Vector R -> [StationaryDistribution] -> [ExchangeabilityMatrix] -> Tree a -> Gen (PrimState m) -> m [[State]] -simulateAndFlattenNSitesAlongTreeMixtureModel n ws ds es t g = do +simulateAndFlattenMixtureModel n ws ds es t g = do let qs = zipWith fromExchangeabilityMatrix es ds - pt = measureableTreeToProbTreeMixtureModel qs t + pt = toProbTreeMixtureModel qs t (cs, is) <- getComponentsAndRootStates n ws ds g - simulateAndFlattenAlongProbTreeMixtureModel is cs pt g + simulateAndFlattenMixtureModel' is cs pt g -simulateAndFlattenAlongProbTreeMixtureModel :: (PrimMonad m) +simulateAndFlattenMixtureModel' :: (PrimMonad m) => [State] -> [Int] -> Tree [ProbMatrix] -> Gen (PrimState m) -> m [[State]] -simulateAndFlattenAlongProbTreeMixtureModel is cs (Node ps f) g +simulateAndFlattenMixtureModel' is cs (Node ps f) g = do is' <- sequence [ jump i (ps !! c) g | (i, c) <- zip is cs ] if null f then return [is'] - else concat <$> sequence [ simulateAndFlattenAlongProbTreeMixtureModel is' cs t g | t <- f ] + else concat <$> sequence [ simulateAndFlattenMixtureModel' is' cs t g | t <- f ] +-- | See 'simulateAndFlattenMixtureModel', parallel version; +-- needs to be run in IO monad. +simulateAndFlattenMixtureModelPar + :: Measurable a + => Int -> Vector R -> [StationaryDistribution] -> [ExchangeabilityMatrix] -> Tree a + -> GenIO -> IO [[[State]]] +simulateAndFlattenMixtureModelPar n ws ds es t g = do + let qs = zipWith fromExchangeabilityMatrix es ds + pt = toProbTreeMixtureModel qs t + parComp n (\n' g' -> do + (cs, is) <- getComponentsAndRootStates n' ws ds g' + simulateAndFlattenMixtureModel' is cs pt g') g + -- | Simulate a number of sites for a given set of substitution models with -- corresponding weights. Keep states at internal nodes. See also --- 'simulateNSitesAlongTree'. -simulateNSitesAlongTreeMixtureModel :: (PrimMonad m, Measurable a) +-- 'simulate'. +simulateMixtureModel :: (PrimMonad m, Measurable a) => Int -> Vector R -> [StationaryDistribution] -> [ExchangeabilityMatrix] -> Tree a -> Gen (PrimState m) -> m (Tree [State]) -simulateNSitesAlongTreeMixtureModel n ws ds es t g = do +simulateMixtureModel n ws ds es t g = do let qs = zipWith fromExchangeabilityMatrix es ds - pt = measureableTreeToProbTreeMixtureModel qs t + pt = toProbTreeMixtureModel qs t (cs, is) <- getComponentsAndRootStates n ws ds g - simulateAlongProbTreeMixtureModel is cs pt g + simulateMixtureModel' is cs pt g -- See 'simulateAlongProbTree', only we have a number of mixture components. The -- starting states and the components for each site have to be provided. -simulateAlongProbTreeMixtureModel :: (PrimMonad m) +simulateMixtureModel' :: (PrimMonad m) => [State] -> [Int] -> Tree [ProbMatrix] -> Gen (PrimState m) -> m (Tree [State]) -simulateAlongProbTreeMixtureModel is cs (Node ps f) g = do +simulateMixtureModel' is cs (Node ps f) g = do is' <- sequence [ jump i (ps !! c) g | (i, c) <- zip is cs ] - f' <- sequence [ simulateAlongProbTreeMixtureModel is' cs t g | t <- f ] + f' <- sequence [ simulateMixtureModel' is' cs t g | t <- f ] return $ Node is' f' diff --git a/test/ELynx/Data/MarkovProcess/AminoAcidSpec.hs b/test/ELynx/Data/MarkovProcess/AminoAcidSpec.hs --- a/test/ELynx/Data/MarkovProcess/AminoAcidSpec.hs +++ b/test/ELynx/Data/MarkovProcess/AminoAcidSpec.hs @@ -14,7 +14,6 @@ module ELynx.Data.MarkovProcess.AminoAcidSpec (spec) where -import Control.Lens import Numeric.LinearAlgebra import Test.Hspec @@ -117,13 +116,13 @@ statDistUniform = vector $ replicate 20 0.05 statDistLG :: R.StationaryDistribution -statDistLG = lg ^. S.stationaryDistribution +statDistLG = S.stationaryDistribution lg exchLG :: R.ExchangeabilityMatrix -exchLG = lg ^. S.exchangeabilityMatrix +exchLG = S.exchangeabilityMatrix lg rmLG :: R.RateMatrix -rmLG = S.getRateMatrix lg +rmLG = S.rateMatrix lg spec :: Spec spec = do @@ -142,14 +141,14 @@ describe "lgCustom" $ it "stationary distribution can be recovered" $ do - let f = R.getStationaryDistribution $ S.getRateMatrix $ lgCustom Nothing statDistUniform + let f = R.getStationaryDistribution $ S.rateMatrix $ lgCustom Nothing statDistUniform f `nearlyEqVec` statDistUniform `shouldBe` True describe "poisson" $ it "stationary distribution is uniform 1/20" $ - R.getStationaryDistribution (S.getRateMatrix poisson) `nearlyEqVec` statDistUniform `shouldBe` True + R.getStationaryDistribution (S.rateMatrix poisson) `nearlyEqVec` statDistUniform `shouldBe` True describe "poissonCustom" $ it "stationary distribution can be recovered" $ do - let f = R.getStationaryDistribution $ S.getRateMatrix $ poissonCustom Nothing statDistLGPython + let f = R.getStationaryDistribution $ S.rateMatrix $ poissonCustom Nothing statDistLGPython f `nearlyEqVec` statDistLGPython `shouldBe` True diff --git a/test/ELynx/Data/MarkovProcess/NucleotideSpec.hs b/test/ELynx/Data/MarkovProcess/NucleotideSpec.hs --- a/test/ELynx/Data/MarkovProcess/NucleotideSpec.hs +++ b/test/ELynx/Data/MarkovProcess/NucleotideSpec.hs @@ -32,5 +32,5 @@ spec = describe "getStationaryDistribution" $ it "extracts the stationary distribution from a rate matrix" $ do - let sd = getStationaryDistribution (getRateMatrix hkyModel) + let sd = getStationaryDistribution (rateMatrix hkyModel) sd `nearlyEqVec` stationaryDist `shouldBe` True diff --git a/test/ELynx/Data/Sequence/AlignmentSpec.hs b/test/ELynx/Data/Sequence/AlignmentSpec.hs new file mode 100644 --- /dev/null +++ b/test/ELynx/Data/Sequence/AlignmentSpec.hs @@ -0,0 +1,40 @@ +{- | +Module : ELynx.Data.Sequence.AlignmentSpec +Copyright : (c) Dominik Schrempf 2018 +License : GPL-3 + +Maintainer : dominik.schrempf@gmail.com +Stability : unstable +Portability : portable + +Creation date: Fri Oct 5 14:25:42 2018. + +-} + +module ELynx.Data.Sequence.AlignmentSpec + (spec) where + +import qualified Data.ByteString.Lazy.Char8 as L +import qualified Data.Matrix.Unboxed as M +import Test.Hspec + +import ELynx.Data.Alphabet.Alphabet +import ELynx.Data.Alphabet.Character +import ELynx.Data.Sequence.Alignment +import ELynx.Import.Sequence.Fasta +import ELynx.Tools.InputOutput +import Files + +ssData :: M.Matrix Character +ssData = M.fromLists $ map (reverse . map fromChar) [ "AAA", "GAA", "TAA" ] + +ssA :: Alignment +ssA = Alignment (map L.pack ["SEQUENCE_1", "SEQUENCE_2", "SEQUENCE_3"]) DNAI ssData + +spec :: Spec +spec = + describe "subSample" $ + it "correctly sub sample an Alignment" $ do + a <- either error id . fromSequences <$> parseFileWith (fasta DNAI) fastaNucleotideIUPACFN + let ss = subSample [0,3,5] a + ss `shouldBe` ssA diff --git a/test/ELynx/Data/Sequence/MultiSequenceAlignmentSpec.hs b/test/ELynx/Data/Sequence/MultiSequenceAlignmentSpec.hs deleted file mode 100644 --- a/test/ELynx/Data/Sequence/MultiSequenceAlignmentSpec.hs +++ /dev/null @@ -1,40 +0,0 @@ -{- | -Module : ELynx.Data.Sequence.MultiSequenceAlignmentSpec -Copyright : (c) Dominik Schrempf 2018 -License : GPL-3 - -Maintainer : dominik.schrempf@gmail.com -Stability : unstable -Portability : portable - -Creation date: Fri Oct 5 14:25:42 2018. - --} - -module ELynx.Data.Sequence.MultiSequenceAlignmentSpec - (spec) where - -import qualified Data.ByteString.Lazy.Char8 as L -import qualified Data.Matrix.Unboxed as M -import Test.Hspec - -import ELynx.Data.Alphabet.Alphabet -import ELynx.Data.Alphabet.Character -import ELynx.Data.Sequence.MultiSequenceAlignment -import ELynx.Import.Sequence.Fasta -import ELynx.Tools.InputOutput -import Files - -ssData :: M.Matrix Character -ssData = M.fromLists $ map (reverse . map fromChar) [ "AAA", "GAA", "TAA" ] - -ssMSA :: MultiSequenceAlignment -ssMSA = MultiSequenceAlignment (map L.pack ["SEQUENCE_1", "SEQUENCE_2", "SEQUENCE_3"]) DNAI ssData - -spec :: Spec -spec = - describe "subSample" $ - it "correctly sub sample an MSA" $ do - msa <- either error id . fromSequenceList <$> parseFileWith (fasta DNAI) fastaNucleotideIUPACFN - let ss = subSample [0,3,5] msa - ss `shouldBe` ssMSA diff --git a/test/ELynx/Import/MarkovProcess/SiteprofilesPhylobayesSpec.hs b/test/ELynx/Import/MarkovProcess/SiteprofilesPhylobayesSpec.hs new file mode 100644 --- /dev/null +++ b/test/ELynx/Import/MarkovProcess/SiteprofilesPhylobayesSpec.hs @@ -0,0 +1,46 @@ +{- | +Module : ELynx.Import.MarkovProcess.SiteprofilesPhylobayesSpec +Copyright : (c) Dominik Schrempf 2019 +License : GPL-3 + +Maintainer : dominik.schrempf@gmail.com +Stability : unstable +Portability : portable + +Creation date: Wed Sep 18 09:51:35 2019. + +-} + +module ELynx.Import.MarkovProcess.SiteprofilesPhylobayesSpec + ( spec + ) where + +import qualified Data.Vector.Storable as V +import Test.Hspec + +import ELynx.Data.MarkovProcess.EDMModel +import ELynx.Import.MarkovProcess.SiteprofilesPhylobayes +import ELynx.Tools.Equality +import ELynx.Tools.InputOutput + +fn :: FilePath +fn = "data/HSSPMany.siteprofiles" + +getProfiles :: IO [EDMComponent] +getProfiles = parseFileWith siteprofiles fn + +firstProfile :: V.Vector Double +firstProfile = V.fromList + [ 0.0267009, 0.013874, 0.022432, 0.0440502, 0.0485174, 0.0407515 + , 0.0170806, 0.348043, 0.0371379, 0.0715536, 0.0454168, 0.0342213 + , 0.0146872, 0.0239681, 0.0768379, 0.0210387, 0.0123336 + , 0.0512678, 0.0149093, 0.0351776 ] + +spec :: Spec +spec = describe "import phylobayes siteprofiles" $ + it "parses a text file with siteprofiles in phylobayes format" $ do + profiles <- getProfiles + length profiles `shouldBe` 701 + map fst profiles `shouldBe` replicate 701 1.0 + map (V.sum . snd) profiles `shouldSatisfy` nearlyEqListWith 1e-5 (replicate 701 1.0) + snd (head profiles) `shouldBe` firstProfile diff --git a/test/ELynx/Import/Sequence/FastaSpec.hs b/test/ELynx/Import/Sequence/FastaSpec.hs --- a/test/ELynx/Import/Sequence/FastaSpec.hs +++ b/test/ELynx/Import/Sequence/FastaSpec.hs @@ -16,7 +16,7 @@ import Data.Either import ELynx.Data.Alphabet.Alphabet -import ELynx.Data.Sequence.MultiSequenceAlignment +import qualified ELynx.Data.Sequence.Alignment as M import ELynx.Import.Sequence.Fasta import ELynx.Tools.InputOutput import Files @@ -24,27 +24,27 @@ spec :: Spec spec = - describe "fastaFileMSA" $ do + describe "fastaFileAlignment" $ do it "parses a fasta file with nucleotide sequences with equal length" $ do - msa <- either error id . fromSequenceList <$> parseFileWith (fasta DNA) fastaNucleotideFN - msaNSequences msa `shouldBe` 3 - msaLength msa `shouldBe` 40 + a <- either error id . M.fromSequences <$> parseFileWith (fasta DNA) fastaNucleotideFN + M.nSequences a `shouldBe` 3 + M.length a `shouldBe` 40 it "parses a fasta file with nucleotide IUPAC sequences with equal length" $ do - msa <- either error id . fromSequenceList <$> parseFileWith (fasta DNAI) fastaNucleotideIUPACFN - msaNSequences msa `shouldBe` 3 - msaLength msa `shouldBe` 40 + a <- either error id . M.fromSequences <$> parseFileWith (fasta DNAI) fastaNucleotideIUPACFN + M.nSequences a `shouldBe` 3 + M.length a `shouldBe` 40 it "should not parse erroneous files" $ do - emsa <- runParserOnFile (fasta DNAI) fastaErroneousFN - emsa `shouldSatisfy` isLeft + ea <- runParserOnFile (fasta DNAI) fastaErroneousFN + ea `shouldSatisfy` isLeft it "parses a fasta file with amino acid sequences with equal length" $ do - msa <- either error id . fromSequenceList <$> parseFileWith (fasta Protein) fastaAminoAcidFN - msaNSequences msa `shouldBe` 2 - msaLength msa `shouldBe` 237 + a <- either error id . M.fromSequences <$> parseFileWith (fasta Protein) fastaAminoAcidFN + M.nSequences a `shouldBe` 2 + M.length a `shouldBe` 237 it "should not parse erroneous files" $ do - msa <- runParserOnFile (fasta ProteinI) fastaErroneousFN - msa `shouldSatisfy` isLeft + a <- runParserOnFile (fasta ProteinI) fastaErroneousFN + a `shouldSatisfy` isLeft diff --git a/test/ELynx/Simulate/MarkovProcessAlongTreeSpec.hs b/test/ELynx/Simulate/MarkovProcessAlongTreeSpec.hs --- a/test/ELynx/Simulate/MarkovProcessAlongTreeSpec.hs +++ b/test/ELynx/Simulate/MarkovProcessAlongTreeSpec.hs @@ -14,7 +14,6 @@ module ELynx.Simulate.MarkovProcessAlongTreeSpec (spec) where -import Control.Lens import Data.Tree import System.Random.MWC import Test.Hspec @@ -28,14 +27,14 @@ -- sampleNewickText :: T.Text -- sampleNewickText = T.pack "(Aeropyrum0:0.5478645225,(((((((((Arabidopsi:0.0701001024,Oryza_sati:0.0765988261):0.0309636193,Gymnosperm:0.0520325624):0.0338982245,Physcomitr:0.0768008916):0.0895714685,(Chlamydomo:0.1136227755,Dunaliella:0.1406347323):0.1117340620):0.0818876186,Rhodophyta:0.3405656487):0.0363527066,((((((Babesia_bo:0.1646969208,Theileria0:0.1519889486):0.1908081096,Plasmodium:0.3250696762):0.0637865908,(Toxoplasma:0.1153570425,Eimeria000:0.1671916078):0.0980136930):0.0518956330,Cryptospor:0.3175062809):0.1607708388,Ciliophora:0.5687502950):0.0624078848,(Phytophtho:0.2016424948,((Thalassios:0.1202730781,Phaeodacty:0.1290341329):0.1772775509,Phaeophyce:0.1989260715):0.0312359673):0.1154768302):0.0311952864):0.0149160316,(((((((((Candida_al:0.1027755272,Saccharomy:0.1190206560):0.1333487870,Neurospora:0.1977309079):0.0522926266,Schizosacc:0.2019603227):0.0567441011,(Cryptococc:0.1948614959,Ustilago_m:0.1564451295):0.0775729694):0.0323959951,Glomus_int:0.1573670796):0.0194701292,Chytridiom:0.2228415254):0.0384370601,Encephalit:1.4622174644):0.0416231688,(((Drosophila:0.2160627753,(Mammalians:0.1080484094,Tunicates0:0.1739253014):0.0289624371):0.0346633757,Hydrozoa00:0.2058137032):0.0480963050,Monosiga_b:0.3020637584):0.0654894239):0.0380915725,(Dictyostel:0.3453588998,Mastigamoe:0.3844779231):0.0478795653):0.0129578395):1.7592083381,((Archaeoglo:0.5402784445,Methanococ:0.4088567459):0.0993669265,Pyrococcus:0.4058713829):0.1734405968):0.2193511807,Pyrobaculu:0.7507718047):0.1646616482,Sulfolobus:0.5404967897);" -testTree :: Tree PhyloIntLabel -testTree = Node (PhyloLabel 0 Nothing 0.0) [ Node (PhyloLabel 1 Nothing 1.0) [], Node (PhyloLabel 2 Nothing 1.0) [] ] +testTree :: Tree (PhyloLabel Int) +testTree = Node (PhyloLabel 0 Nothing (Just 0.0)) [ Node (PhyloLabel 1 Nothing (Just 1.0)) [], Node (PhyloLabel 2 Nothing (Just 1.0)) [] ] -- testStateTree :: Tree State -- testStateTree = Node 0 [Node 3 [], Node 2 []] -nullTree :: Tree PhyloIntLabel -nullTree = Node (PhyloLabel 0 Nothing 0.0) [] +nullTree :: Tree (PhyloLabel Int) +nullTree = Node (PhyloLabel 0 Nothing (Just 0.0)) [] nullStateTree :: Tree State nullStateTree = fmap (const 0) nullTree @@ -48,13 +47,13 @@ describe "simulateNSitesAlongTree" $ do it "simulates one site along an easy tree" $ do gen <- create - tr <- simulateNSitesAlongTree 1 d e nullTree gen + tr <- simulate 1 d e nullTree gen fmap head tr `shouldBe` nullStateTree it "simulates some sites along a harder tree" $ do gen <- create - tr <- simulateNSitesAlongTree 10 d e testTree gen + tr <- simulate 10 d e testTree gen (length . rootLabel $ tr) `shouldBe` 10 - where d = jc ^. S.stationaryDistribution - e = jc ^. S.exchangeabilityMatrix + where d = S.stationaryDistribution jc + e = S.exchangeabilityMatrix jc