diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -1,10 +1,13 @@ -# Table of Contents - -1. [Changelog for ELynx](#org263dcac) +# Changelog for ELynx - +## Version 0.2.1 -# Changelog for ELynx +- Validation and repetition of previous analyses is finally possible with the + new `elynx` binary. +- A library `elynx-markov` for running Markov processes along phylogenetic trees + has been split off `elynx-seq`. This library performs the computations when + executing `slynx simulate ...`. +- Many other small improvements. diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -2,30 +2,31 @@ # The ELynx Suite -Version: 0.1.0. +Version: 0.2.1. 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 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. +A Haskell library and tool set for computational biology. The goal of ELynx is +reproducible research. Evolutionary sequences and phylogenetic trees can be +read, viewed, modified and simulated. The command line with all arguments is +logged consistently, and automatically. Data integrity is verified using SHA256 +sums so that validation of past analyses is possible without the need to +recompute the result. -The Elynx Suite consists of three library packages and two executables providing +The Elynx Suite consists of four library packages and three executables providing a range of sub commands. The library packages are: +- **elynx-markov:** Simulate multi sequence alignments along phylogenetic trees. - **elynx-seq:** Handle evolutionary sequences and multi sequence alignments. -- **elynx-tree:** Handle phylogenetic trees. - **elynx-tools:** Tools for the provided executables. +- **elynx-tree:** Handle phylogenetic trees. 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. +- **tlynx:** Analyze, modify, and simulate phylogenetic trees. +- **elynx:** Validate and redo past analyses. **ELynx is actively developed. We happily receive comments, ideas, feature requests, and pull requests!** @@ -67,10 +68,12 @@ slynx --help - ELynx Suite version 0.1.0. Developed by Dominik Schrempf. Compiled on January - 30, 2020, at 09:18 am, UTC. + ELynx Suite version 0.2.1. + Developed by Dominik Schrempf. + Compiled on April 27, 2020, at 06:31 am, UTC. - Usage: slynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME] COMMAND + Usage: slynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME] + [-f|--force] COMMAND Analyze, and simulate multi sequence alignments. Available options: @@ -80,21 +83,22 @@ Debug (default: Info) -o,--output-file-basename NAME Specify base name of output file + -f,--force Ignore previous analysis and overwrite existing + output files. Available commands: concatenate Concatenate sequences found in input files. - examine Examine sequences. If data is a multi sequence - alignment, additionally analyze columns. + examine Examine sequences. If data is a multi sequence alignment, additionally analyze columns. + filter-columns Filter columns of multi sequence alignments. 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: + Available sequence file formats: - FASTA - Alphabet types: + Available alphabets: - DNA (nucleotides) - DNAX (nucleotides; including gaps) - DNAI (nucleotides; including gaps, and IUPAC codes) @@ -103,19 +107,20 @@ - ProteinS (amino acids; including gaps, and translation stops) - ProteinI (amino acids; including gaps, translation stops, and IUPAC codes) - 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. + ELynx + ----- + A Haskell library and tool set for computational biology. The goal of ELynx is + reproducible research. Evolutionary sequences and phylogenetic trees can be + read, viewed, modified and simulated. The command line with all arguments is + logged consistently, and automatically. Data integrity is verified using SHA256 + sums so that validation of past analyses is possible without the need to + recompute the result. slynx Analyze, modify, and simulate evolutionary sequences. tlynx Analyze, modify, and simulate phylogenetic trees. + elynx Validate and redo past analyses. - Get help for specific commands: + Get help for sub commands: slynx examine --help @@ -125,13 +130,36 @@ slynx concatenate --help + ELynx Suite version 0.2.1. + Developed by Dominik Schrempf. + Compiled on April 27, 2020, at 06:31 am, UTC. + Usage: slynx concatenate (-a|--alphabet NAME) INPUT-FILE Concatenate sequences found in input files. Available options: + -h,--help Show this help text + -V,--version Show version -a,--alphabet NAME Specify alphabet type NAME INPUT-FILE Read sequences from INPUT-FILE -h,--help Show this help text + + + ELynx + ----- + A Haskell library and tool set for computational biology. The goal of ELynx is + reproducible research. Evolutionary sequences and phylogenetic trees can be + read, viewed, modified and simulated. The command line with all arguments is + logged consistently, and automatically. Data integrity is verified using SHA256 + sums so that validation of past analyses is possible without the need to + recompute the result. + + slynx Analyze, modify, and simulate evolutionary sequences. + tlynx Analyze, modify, and simulate phylogenetic trees. + elynx Validate and redo past analyses. + + Get help for sub commands: + slynx examine --help ## Examine @@ -140,15 +168,37 @@ slynx examine --help - Usage: slynx examine (-a|--alphabet NAME) [INPUT-FILE] [--per-site] - Examine sequences. If data is a multi sequence alignment, additionally analyze - columns. + ELynx Suite version 0.2.1. + Developed by Dominik Schrempf. + Compiled on April 27, 2020, at 06:31 am, UTC. + Usage: slynx examine (-a|--alphabet NAME) INPUT-FILE [--per-site] + Examine sequences. If data is a multi sequence alignment, additionally analyze columns. + Available options: + -h,--help Show this help text + -V,--version Show version -a,--alphabet NAME Specify alphabet type NAME INPUT-FILE Read sequences from INPUT-FILE --per-site Report per site summary statistics -h,--help Show this help text + + + ELynx + ----- + A Haskell library and tool set for computational biology. The goal of ELynx is + reproducible research. Evolutionary sequences and phylogenetic trees can be + read, viewed, modified and simulated. The command line with all arguments is + logged consistently, and automatically. Data integrity is verified using SHA256 + sums so that validation of past analyses is possible without the need to + recompute the result. + + slynx Analyze, modify, and simulate evolutionary sequences. + tlynx Analyze, modify, and simulate phylogenetic trees. + elynx Validate and redo past analyses. + + Get help for sub commands: + slynx examine --help ## Filter @@ -157,12 +207,17 @@ slynx filter-rows --help - Usage: slynx filter-rows (-a|--alphabet NAME) [INPUT-FILE] - [--longer-than LENGTH] [--shorter-than LENGTH] - [--standard-characters] + ELynx Suite version 0.2.1. + Developed by Dominik Schrempf. + Compiled on April 27, 2020, at 06:31 am, UTC. + + 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: + -h,--help Show this help text + -V,--version Show version -a,--alphabet NAME Specify alphabet type NAME INPUT-FILE Read sequences from INPUT-FILE --longer-than LENGTH Only keep sequences longer than LENGTH @@ -170,21 +225,61 @@ --standard-characters Only keep sequences containing at least one standard (i.e., non-IUPAC) character -h,--help Show this help text + + + ELynx + ----- + A Haskell library and tool set for computational biology. The goal of ELynx is + reproducible research. Evolutionary sequences and phylogenetic trees can be + read, viewed, modified and simulated. The command line with all arguments is + logged consistently, and automatically. Data integrity is verified using SHA256 + sums so that validation of past analyses is possible without the need to + recompute the result. + + slynx Analyze, modify, and simulate evolutionary sequences. + tlynx Analyze, modify, and simulate phylogenetic trees. + elynx Validate and redo past analyses. + + Get help for sub commands: + slynx examine --help Filter columns of multi sequence alignments with `filter-columns`. slynx filter-columns --help - Usage: slynx filter-columns (-a|--alphabet NAME) [INPUT-FILE] + ELynx Suite version 0.2.1. + Developed by Dominik Schrempf. + Compiled on April 27, 2020, at 06:31 am, UTC. + + Usage: slynx filter-columns (-a|--alphabet NAME) INPUT-FILE [--standard-chars DOUBLE] - Filter columns of multi-sequence alignments. + Filter columns of multi sequence alignments. Available options: + -h,--help Show this help text + -V,--version Show version -a,--alphabet NAME Specify alphabet type NAME INPUT-FILE Read sequences from INPUT-FILE --standard-chars DOUBLE Keep columns with a proportion standard (non-IUPAC) characters larger than DOUBLE in [0,1] -h,--help Show this help text + + + ELynx + ----- + A Haskell library and tool set for computational biology. The goal of ELynx is + reproducible research. Evolutionary sequences and phylogenetic trees can be + read, viewed, modified and simulated. The command line with all arguments is + logged consistently, and automatically. Data integrity is verified using SHA256 + sums so that validation of past analyses is possible without the need to + recompute the result. + + slynx Analyze, modify, and simulate evolutionary sequences. + tlynx Analyze, modify, and simulate phylogenetic trees. + elynx Validate and redo past analyses. + + Get help for sub commands: + slynx examine --help ## Simulate @@ -193,16 +288,22 @@ slynx simulate --help - 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,...]"] + ELynx Suite version 0.2.1. + Developed by Dominik Schrempf. + Compiled on April 27, 2020, at 06:31 am, UTC. + + 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 + -h,--help Show this help text + -V,--version Show version + -t,--tree-file Name Read tree from Newick file NAME -s,--substitution-model MODEL Set the phylogenetic substitution model; available models are shown below (mutually exclusive with -m @@ -224,17 +325,24 @@ Substitution models: -s "MODEL[PARAMETER,PARAMETER,...]{STATIONARY_DISTRIBUTION}" - Supported DNA models: JC, HKY. + Supported DNA models: JC, F81, HKY, GTR4. For example, -s HKY[KAPPA]{DOUBLE,DOUBLE,DOUBLE,DOUBLE} - Supported Protein models: Poisson, Poisson-Custom, LG, LG-Custom, WAG, WAG-Custom. + -s GTR4[e_AC,e_AG,e_AT,e_CG,e_CT,e_GT]{DOUBLE,DOUBLE,DOUBLE,DOUBLE} + where the 'e_XY' are the exchangeabilities from nucleotide X to Y. + Supported Protein models: Poisson, Poisson-Custom, LG, LG-Custom, WAG, WAG-Custom, GTR20. MODEL-Custom means that only the exchangeabilities of MODEL are used, and a custom stationary distribution is provided. For example, + -s LG -s LG-Custom{...} + -s GTR20[e_AR,e_AN,...]{...} + the 'e_XY' are the exchangeabilities from amino acid X to Y (alphabetical order). + Notes: The F81 model for DNA is equivalent to the Poisson-Custom for proteins. + The GTR4 model for DNA is equivalent to the GTR20 for proteins. Mixture models: - -m "MIXTURE(SUBSTITUTION_MODEL_1,SUBSTITUTION_MODEL_2)" + -m "MIXTURE(SUBSTITUTION_MODEL_1,SUBSTITUTION_MODEL_2[PARAMETERS]{STATIONARY_DISTRIBUTION},...)" For example, -m "MIXTURE(JC,HKY[6.0]{0.3,0.2,0.2,0.3})" Mixture weights have to be provided with the -w option. @@ -249,6 +357,22 @@ LG exchangeabilities with stationary distributions given in FILE. -m "EDM(LG-Custom)" -e FILE For special mixture models, mixture weights are optional. + + ELynx + ----- + A Haskell library and tool set for computational biology. The goal of ELynx is + reproducible research. Evolutionary sequences and phylogenetic trees can be + read, viewed, modified and simulated. The command line with all arguments is + logged consistently, and automatically. Data integrity is verified using SHA256 + sums so that validation of past analyses is possible without the need to + recompute the result. + + slynx Analyze, modify, and simulate evolutionary sequences. + tlynx Analyze, modify, and simulate phylogenetic trees. + elynx Validate and redo past analyses. + + Get help for sub commands: + slynx examine --help ## Sub-sample @@ -257,12 +381,18 @@ slynx sub-sample --help - Usage: slynx sub-sample (-a|--alphabet NAME) [INPUT-FILE] + ELynx Suite version 0.2.1. + Developed by Dominik Schrempf. + Compiled on April 27, 2020, at 06:31 am, UTC. + + Usage: slynx sub-sample (-a|--alphabet NAME) INPUT-FILE (-n|--number-of-sites INT) (-m|--number-of-alignments INT) [-S|--seed [INT]] Sub-sample columns from multi sequence alignments. Available options: + -h,--help Show this help text + -V,--version Show version -a,--alphabet NAME Specify alphabet type NAME INPUT-FILE Read sequences from INPUT-FILE -n,--number-of-sites INT Number of sites randomly drawn with replacement @@ -272,8 +402,23 @@ 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. + 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. + + ELynx + ----- + A Haskell library and tool set for computational biology. The goal of ELynx is + reproducible research. Evolutionary sequences and phylogenetic trees can be + read, viewed, modified and simulated. The command line with all arguments is + logged consistently, and automatically. Data integrity is verified using SHA256 + sums so that validation of past analyses is possible without the need to + recompute the result. + + slynx Analyze, modify, and simulate evolutionary sequences. + tlynx Analyze, modify, and simulate phylogenetic trees. + elynx Validate and redo past analyses. + + Get help for sub commands: + slynx examine --help ## Translate @@ -282,17 +427,40 @@ slynx translate --help - Usage: slynx translate (-a|--alphabet NAME) [INPUT-FILE] - (-r|--reading-frame INT) (-u|--universal-code CODE) + ELynx Suite version 0.2.1. + Developed by Dominik Schrempf. + Compiled on April 27, 2020, at 06:31 am, UTC. + + 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: + -h,--help Show this help text + -V,--version Show version -a,--alphabet NAME Specify alphabet type NAME INPUT-FILE Read sequences from INPUT-FILE -r,--reading-frame INT Reading frame [0|1|2]. -u,--universal-code CODE universal code; one of: Standard, VertebrateMitochondrial. -h,--help Show this help text + + + ELynx + ----- + A Haskell library and tool set for computational biology. The goal of ELynx is + reproducible research. Evolutionary sequences and phylogenetic trees can be + read, viewed, modified and simulated. The command line with all arguments is + logged consistently, and automatically. Data integrity is verified using SHA256 + sums so that validation of past analyses is possible without the need to + recompute the result. + + slynx Analyze, modify, and simulate evolutionary sequences. + tlynx Analyze, modify, and simulate phylogenetic trees. + elynx Validate and redo past analyses. + + Get help for sub commands: + slynx examine --help # TLynx @@ -301,10 +469,12 @@ tlynx --help - ELynx Suite version 0.1.0. Developed by Dominik Schrempf. Compiled on January - 30, 2020, at 09:18 am, UTC. + ELynx Suite version 0.2.1. + Developed by Dominik Schrempf. + Compiled on April 27, 2020, at 06:31 am, UTC. - Usage: tlynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME] COMMAND + Usage: tlynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME] + [-f|--force] COMMAND Compare, examine, and simulate phylogenetic trees. Available options: @@ -314,35 +484,37 @@ Debug (default: Info) -o,--output-file-basename NAME Specify base name of output file + -f,--force Ignore previous analysis and overwrite existing + output files. Available commands: + coalesce Simulate phylogenetic trees using the coalescent processes (see also the 'simulate' command for simulations using the birth and death process). + compare Compare two phylogenetic trees (compute distances and branch-wise differences). + connect Connect two phylogenetic trees in all ways (possibly honoring constraints). 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). + shuffle Shuffle a phylogenetic tree (keep coalescent times, but shuffle topology and leaves). + simulate Simulate phylogenetic trees using birth and death processes (see also the 'coalesce' command for simulations using the coalescent process). - File formats: - - Newick + Available tree file formats: + - Newick Standard: Branch support values are stored in square brackets after branch lengths. + - Newick IqTree: Branch support values are stored as node names after the closing bracket of forests. + - Newick RevBayes Key-value pairs is provided in square brackets after node names as well as branch lengths. XXX: Key value pairs are IGNORED at the moment. - 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. + ELynx + ----- + A Haskell library and tool set for computational biology. The goal of ELynx is + reproducible research. Evolutionary sequences and phylogenetic trees can be + read, viewed, modified and simulated. The command line with all arguments is + logged consistently, and automatically. Data integrity is verified using SHA256 + sums so that validation of past analyses is possible without the need to + recompute the result. slynx Analyze, modify, and simulate evolutionary sequences. tlynx Analyze, modify, and simulate phylogenetic trees. + elynx Validate and redo past analyses. - Get help for specific commands: + Get help for sub commands: slynx examine --help @@ -352,18 +524,43 @@ tlynx compare --help - Usage: tlynx compare [-n|--normalize] [-b|--bipartitions] [-i|--newick-iqtree] - NAME - Compare two phylogenetic trees (compute distances and branch-wise - differences). + ELynx Suite version 0.2.1. + Developed by Dominik Schrempf. + Compiled on April 27, 2020, at 06:31 am, UTC. + Usage: tlynx compare [-n|--normalize] [-b|--bipartitions] [-t|--intersect] + [-f|--newick-format FORMAT] NAME + Compare two phylogenetic trees (compute distances and branch-wise differences). + Available options: + -h,--help Show this help text + -V,--version Show version -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) + -b,--bipartitions Print and plot common and missing bipartitions + -t,--intersect Compare intersections; i.e., before comparison, drop + leaves that are not present in the other tree + -f,--newick-format FORMAT + Newick tree format; see 'tlynx + --help' (default: Standard) NAME Tree file -h,--help Show this help text + + + ELynx + ----- + A Haskell library and tool set for computational biology. The goal of ELynx is + reproducible research. Evolutionary sequences and phylogenetic trees can be + read, viewed, modified and simulated. The command line with all arguments is + logged consistently, and automatically. Data integrity is verified using SHA256 + sums so that validation of past analyses is possible without the need to + recompute the result. + + slynx Analyze, modify, and simulate evolutionary sequences. + tlynx Analyze, modify, and simulate phylogenetic trees. + elynx Validate and redo past analyses. + + Get help for sub commands: + slynx examine --help ## Examine @@ -372,14 +569,38 @@ tlynx examine --help - Usage: tlynx examine [INPUT-FILE] [-i|--newick-iqtree] + ELynx Suite version 0.2.1. + Developed by Dominik Schrempf. + Compiled on April 27, 2020, at 06:31 am, UTC. + + Usage: tlynx examine INPUT-FILE [-f|--newick-format FORMAT] Compute summary statistics of phylogenetic trees. Available options: + -h,--help Show this help text + -V,--version Show version INPUT-FILE Read trees from INPUT-FILE - -i,--newick-iqtree Use IQ-TREE Newick format (internal node labels are - branch support values) + -f,--newick-format FORMAT + Newick tree format; see 'tlynx + --help' (default: Standard) -h,--help Show this help text + + + ELynx + ----- + A Haskell library and tool set for computational biology. The goal of ELynx is + reproducible research. Evolutionary sequences and phylogenetic trees can be + read, viewed, modified and simulated. The command line with all arguments is + logged consistently, and automatically. Data integrity is verified using SHA256 + sums so that validation of past analyses is possible without the need to + recompute the result. + + slynx Analyze, modify, and simulate evolutionary sequences. + tlynx Analyze, modify, and simulate phylogenetic trees. + elynx Validate and redo past analyses. + + Get help for sub commands: + slynx examine --help ## Simulate @@ -388,13 +609,19 @@ tlynx simulate --help - 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] + ELynx Suite version 0.2.1. + Developed by Dominik Schrempf. + Compiled on April 27, 2020, at 06:31 am, UTC. + + 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 phylogenetic trees using birth and death processes. + Simulate phylogenetic trees using birth and death processes (see also the 'coalesce' command for simulations using the coalescent process). Available options: + -h,--help Show this help text + -V,--version Show version -t,--nTrees INT Number of trees (default: 10) -n,--nLeaves INT Number of leaves per tree (default: 5) -H,--height DOUBLE Fix tree height (no default) @@ -413,15 +640,75 @@ 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 + ----- + A Haskell library and tool set for computational biology. The goal of ELynx is + reproducible research. Evolutionary sequences and phylogenetic trees can be + read, viewed, modified and simulated. The command line with all arguments is + logged consistently, and automatically. Data integrity is verified using SHA256 + sums so that validation of past analyses is possible without the need to + recompute the result. + + slynx Analyze, modify, and simulate evolutionary sequences. + tlynx Analyze, modify, and simulate phylogenetic trees. + elynx Validate and redo past analyses. + + Get help for sub commands: + slynx examine --help # ELynx -Documentation of the library can be found on [Hackage](https://hackage.haskell.org/): +Validate and (optionally) redo past ELynx analyses. + elynx --help + + ELynx Suite version 0.2.1. + Developed by Dominik Schrempf. + Compiled on April 27, 2020, at 06:31 am, UTC. + + Usage: elynx COMMAND + Validate and redo past ELynx analyses + + Available options: + -h,--help Show this help text + -V,--version Show version + + Available commands: + validate Validate an ELynx analysis + redo Redo an ELynx analysis + + + ELynx + ----- + A Haskell library and tool set for computational biology. The goal of ELynx is + reproducible research. Evolutionary sequences and phylogenetic trees can be + read, viewed, modified and simulated. The command line with all arguments is + logged consistently, and automatically. Data integrity is verified using SHA256 + sums so that validation of past analyses is possible without the need to + recompute the result. + + slynx Analyze, modify, and simulate evolutionary sequences. + tlynx Analyze, modify, and simulate phylogenetic trees. + elynx Validate and redo past analyses. + + Get help for sub commands: + slynx examine --help + + +# Library documentation + +Documentation of the libraries can be found on [Hackage](https://hackage.haskell.org/): + +- [elynx-markov](https://hackage.haskell.org/package/elynx-markov) - [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) +- [elynx-tree](https://hackage.haskell.org/package/elynx-tree) + +Documentation of the executables is also available: + +- [elynx](https://hackage.haskell.org/package/elynx) - [slynx](https://hackage.haskell.org/package/slynx) - [tlynx](https://hackage.haskell.org/package/tlynx) diff --git a/app/Main.hs b/app/Main.hs --- a/app/Main.hs +++ b/app/Main.hs @@ -1,8 +1,8 @@ {- | Module : Main Description : Work with molecular sequence data -Copyright : (c) Dominik Schrempf 2019 -License : GPL-3 +Copyright : (c) Dominik Schrempf 2020 +License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable @@ -14,31 +14,7 @@ module Main where -import Control.Monad.Trans.Reader - -import Options - -import Concatenate.Concatenate -import Examine.Examine -import Filter.Filter -import Simulate.Simulate -import SubSample.SubSample -import Translate.Translate - -import ELynx.Tools.Logger -import ELynx.Tools.Options +import SLynx.SLynx main :: IO () -main = do - Arguments g c <- parseArguments - let fn = outFileBaseName g - lvl = verbosity g - lf = (++ ".log") <$> fn - case c of - Concatenate a -> runReaderT (eLynxWrapper lvl lf concatenateDescription $ concatenateCmd fn) a - Examine a -> runReaderT (eLynxWrapper lvl lf examineDescription $ examineCmd fn) a - FilterRows a -> runReaderT (eLynxWrapper lvl lf filterRowsDescription $ filterRowsCmd fn) a - FilterCols a -> runReaderT (eLynxWrapper lvl lf filterColumnsDescription $ filterColsCmd fn) a - Simulate a -> runReaderT (eLynxWrapper lvl lf simulateDescription $ simulateCmd fn) a - SubSample a -> runReaderT (eLynxWrapper lvl lf subSampleDescription $ subSampleCmd fn) a - Translate a -> runReaderT (eLynxWrapper lvl lf translateDescription $ translateCmd fn) a +main = rSLynx diff --git a/app/Options.hs b/app/Options.hs deleted file mode 100644 --- a/app/Options.hs +++ /dev/null @@ -1,126 +0,0 @@ -{- | -Module : Options -Description : SLynx general options -Copyright : (c) Dominik Schrempf 2019 -License : GPL-3 - -Maintainer : dominik.schrempf@gmail.com -Stability : unstable -Portability : portable - -Creation date: Sat Sep 7 18:55:03 2019. - --} - -module Options - ( Arguments (..) - , CommandArguments (..) - , parseArguments - , concatenateDescription - , examineDescription - , filterRowsDescription - , filterColumnsDescription - , simulateDescription - , subSampleDescription - , translateDescription - ) where - -import Options.Applicative -import Options.Applicative.Help.Pretty - -import Concatenate.Options -import Examine.Options -import Filter.Options -import Simulate.Options -import SubSample.Options -import Translate.Options - -import ELynx.Data.Alphabet.Alphabet -import ELynx.Tools.Options --- TODO: Removed reproducibility stuff for now, have to work on that. --- import ELynx.Tools.Reproduction - -data CommandArguments = - Concatenate ConcatenateArguments - | Examine ExamineArguments - | FilterRows FilterRowsArguments - | FilterCols FilterColsArguments - | Simulate SimulateArguments - | SubSample SubSampleArguments - | Translate TranslateArguments - --- instance Reproducible CommandArguments where --- inFiles (Concatenate a) = inFiles a - -concatenateDescription, examineDescription, filterRowsDescription, filterColumnsDescription, simulateDescription, subSampleDescription, translateDescription :: String -concatenateDescription = "Concatenate sequences found in input files." -examineDescription = "Examine sequences. If data is a multi sequence alignment, additionally analyze columns." -filterRowsDescription = "Filter rows (or sequences) found in input files." -filterColumnsDescription = "Filter columns of multi-sequence alignments." -simulateDescription = "Simulate multi sequence alignments." -subSampleDescription = "Sub-sample columns from multi sequence alignments." -translateDescription = "Translate from DNA to Protein or DNAX to ProteinX." - -concatenateCommand :: Mod CommandFields CommandArguments -concatenateCommand = command "concatenate" $ - info (Concatenate <$> concatenateArguments) - $ progDesc concatenateDescription - -examineCommand :: Mod CommandFields CommandArguments -examineCommand = command "examine" $ - info (Examine <$> examineArguments) - ( fullDesc - <> progDesc examineDescription ) - -filterRowsCommand :: Mod CommandFields CommandArguments -filterRowsCommand = command "filter-rows" $ - info (FilterRows <$> filterRowsArguments) - $ progDesc filterRowsDescription - -filterColumnsCommand :: Mod CommandFields CommandArguments -filterColumnsCommand = command "filter-columns" $ - info (FilterCols <$> filterColsArguments) - $ progDesc filterColumnsDescription - -simulateCommand :: Mod CommandFields CommandArguments -simulateCommand = command "simulate" $ - info (Simulate <$> simulateArguments) - (fullDesc - <> progDesc simulateDescription - <> footerDoc (Just $ pretty simulateFooter) ) - -subSampleCommand :: Mod CommandFields CommandArguments -subSampleCommand = command "sub-sample" $ - info (SubSample <$> subSampleArguments) - ( fullDesc - <> progDesc subSampleDescription - <> footer "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." ) - -translateCommand :: Mod CommandFields CommandArguments -translateCommand = command "translate" $ - info (Translate <$> translateArguments) - $ progDesc translateDescription - -commandArguments :: Parser CommandArguments -commandArguments = hsubparser $ - concatenateCommand - <> examineCommand - <> filterRowsCommand - <> filterColumnsCommand - <> simulateCommand - <> subSampleCommand - <> translateCommand - -parseArguments :: IO (Arguments CommandArguments) -parseArguments = parseArgumentsWith desc ftr commandArguments - -desc :: [String] -desc = [ "Analyze, and simulate multi sequence alignments." ] - -ftr :: [String] -ftr = [ "File formats:" ] ++ fs ++ - [ "", "Alphabet types:" ] ++ as - where - toListItem = (" - " ++) - fs = map toListItem ["FASTA"] - as = map (toListItem . description) [(minBound :: Alphabet) ..] diff --git a/slynx.cabal b/slynx.cabal --- a/slynx.cabal +++ b/slynx.cabal @@ -1,9 +1,9 @@ cabal-version: 1.12 name: slynx -version: 0.1.0 +version: 0.2.1 license: GPL-3 license-file: LICENSE -copyright: Dominik Schrempf (2019) +copyright: Dominik Schrempf (2020) maintainer: dominik.schrempf@gmail.com author: Dominik Schrempf homepage: https://github.com/dschrempf/elynx#readme @@ -23,82 +23,67 @@ library exposed-modules: - Concatenate.Concatenate - Concatenate.Options - Examine.Examine - Examine.Options - Filter.Filter - Filter.Options - Simulate.Options - Simulate.PhyloModel - Simulate.Simulate - SubSample.Options - SubSample.SubSample - Tools - Translate.Options - Translate.Translate + SLynx.Concatenate.Concatenate + SLynx.Concatenate.Options + SLynx.Examine.Examine + SLynx.Examine.Options + SLynx.Filter.Filter + SLynx.Filter.Options + SLynx.Options + SLynx.Simulate.Options + SLynx.Simulate.PhyloModel + SLynx.Simulate.Simulate + SLynx.SLynx + SLynx.SubSample.Options + SLynx.SubSample.SubSample + SLynx.Tools + SLynx.Translate.Options + SLynx.Translate.Translate hs-source-dirs: src other-modules: Paths_slynx default-language: Haskell2010 - ghc-options: -Wall -O2 + ghc-options: -Wall -fllvm 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-seq >=0.1.0 && <0.2, - elynx-tools >=0.1.0 && <0.2, - elynx-tree >=0.1.0 && <0.2, + base >=4.13.0.0 && <4.14, + bytestring >=0.10.10.0 && <0.11, + containers >=0.6.2.1 && <0.7, + elynx-markov >=0.2.1 && <0.3, + elynx-seq >=0.2.1 && <0.3, + elynx-tools >=0.2.1 && <0.3, + elynx-tree >=0.2.1 && <0.3, hmatrix >=0.20.0.0 && <0.21, - integration >=0.2.1 && <0.3, - math-functions >=0.3.2.1 && <0.4, - matrices >=0.5.0 && <0.6, - megaparsec >=7.0.5 && <7.1, - monad-logger >=0.3.30 && <0.4, + megaparsec >=8.0.0 && <8.1, + monad-logger >=0.3.32 && <0.4, mwc-random >=0.14.0.0 && <0.15, - optparse-applicative >=0.14.3.0 && <0.15, - parallel >=3.2.2.0 && <3.3, - primitive >=0.6.4.0 && <0.7, - statistics >=0.15.0.0 && <0.16, - text >=1.2.3.1 && <1.3, + optparse-applicative >=0.15.1.0 && <0.16, + scientific >=0.3.6.2 && <0.4, + text >=1.2.4.0 && <1.3, transformers >=0.5.6.2 && <0.6, - vector >=0.12.0.3 && <0.13, - vector-th-unbox >=0.2.1.6 && <0.3, - word8 >=0.1.3 && <0.2 + vector >=0.12.1.2 && <0.13 executable slynx main-is: Main.hs hs-source-dirs: app other-modules: - Options Paths_slynx default-language: Haskell2010 - ghc-options: -Wall -O2 -eventlog -threaded -rtsopts - -with-rtsopts=-N + ghc-options: -Wall -fllvm -threaded -rtsopts -with-rtsopts=-N 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-seq >=0.1.0 && <0.2, - elynx-tools >=0.1.0 && <0.2, + base >=4.13.0.0 && <4.14, + bytestring >=0.10.10.0 && <0.11, + containers >=0.6.2.1 && <0.7, + elynx-seq >=0.2.1 && <0.3, + elynx-tools >=0.2.1 && <0.3, hmatrix >=0.20.0.0 && <0.21, - integration >=0.2.1 && <0.3, - math-functions >=0.3.2.1 && <0.4, - matrices >=0.5.0 && <0.6, - megaparsec >=7.0.5 && <7.1, - monad-logger >=0.3.30 && <0.4, + megaparsec >=8.0.0 && <8.1, + monad-logger >=0.3.32 && <0.4, mwc-random >=0.14.0.0 && <0.15, - optparse-applicative >=0.14.3.0 && <0.15, - parallel >=3.2.2.0 && <3.3, - primitive >=0.6.4.0 && <0.7, + optparse-applicative >=0.15.1.0 && <0.16, + scientific >=0.3.6.2 && <0.4, slynx -any, - statistics >=0.15.0.0 && <0.16, - text >=1.2.3.1 && <1.3, + text >=1.2.4.0 && <1.3, transformers >=0.5.6.2 && <0.6, - vector >=0.12.0.3 && <0.13, - vector-th-unbox >=0.2.1.6 && <0.3, - word8 >=0.1.3 && <0.2 + vector >=0.12.1.2 && <0.13 diff --git a/src/Concatenate/Concatenate.hs b/src/Concatenate/Concatenate.hs deleted file mode 100644 --- a/src/Concatenate/Concatenate.hs +++ /dev/null @@ -1,43 +0,0 @@ -{-# LANGUAGE OverloadedStrings #-} -{-# LANGUAGE TemplateHaskell #-} - -{- | -Module : Analyze.Analyze -Description : Parse sequence file formats and analyze them -Copyright : (c) Dominik Schrempf 2018 -License : GPL-3 - -Maintainer : dominik.schrempf@gmail.com -Stability : unstable -Portability : portable - -Creation date: Fri Oct 5 08:41:05 2018. - --} - -module Concatenate.Concatenate - ( concatenateCmd - ) -where - -import Control.Monad.Logger -import Control.Monad.Trans.Class -import Control.Monad.Trans.Reader - -import Concatenate.Options -import Tools - -import qualified ELynx.Data.Sequence.Sequence as S -import ELynx.Export.Sequence.Fasta -import ELynx.Tools.InputOutput - --- | Concatenate sequences. -concatenateCmd :: Maybe FilePath -> Concatenate () -concatenateCmd outFileBaseName = do - $(logInfo) "Command: Concatenate sequences." - ConcatenateArguments al fps <- lift ask - sss <- mapM (readSeqs al . Just) fps - let result = sequencesToFasta $ S.concatSequences sss - let outFilePath = (++ ".fasta") <$> outFileBaseName - out "concatenated multi sequence alignment " result outFilePath - diff --git a/src/Concatenate/Options.hs b/src/Concatenate/Options.hs deleted file mode 100644 --- a/src/Concatenate/Options.hs +++ /dev/null @@ -1,52 +0,0 @@ -{- | -Module : Concatenate.Options -Description : ELynxSeq argument parsing -Copyright : (c) Dominik Schrempf 2018 -License : GPL-3 - -Maintainer : dominik.schrempf@gmail.com -Stability : unstable -Portability : portable - -Creation date: Sun Oct 7 17:29:45 2018. - --} - -module Concatenate.Options - ( ConcatenateArguments (..) - , Concatenate - , concatenateArguments - ) where - -import Control.Applicative -import Control.Monad.Logger -import Control.Monad.Trans.Reader -import Options.Applicative - -import Tools - -import ELynx.Data.Alphabet.Alphabet -import ELynx.Tools.Reproduction - --- | Arguments needed to concatenate multi sequence alignments. -data ConcatenateArguments = ConcatenateArguments - { ccAlphabet :: Alphabet - , ccInFiles :: [FilePath] } - -instance Reproducible ConcatenateArguments where - inFiles = ccInFiles - parser _ = concatenateArguments - --- | Logger and Reader type. -type Concatenate = LoggingT (ReaderT ConcatenateArguments IO) - --- | Command line parser. -concatenateArguments :: Parser ConcatenateArguments -concatenateArguments = ConcatenateArguments - <$> alphabetOpt - <*> some inFileArg - -inFileArg :: Parser FilePath -inFileArg = strArgument $ - metavar "INPUT-FILE" <> - help "Read sequences from INPUT-FILE" diff --git a/src/Examine/Examine.hs b/src/Examine/Examine.hs deleted file mode 100644 --- a/src/Examine/Examine.hs +++ /dev/null @@ -1,118 +0,0 @@ -{-# LANGUAGE OverloadedStrings #-} -{-# LANGUAGE TemplateHaskell #-} - -{- | - -Module : Analyze.Analyze -Description : Parse sequence file formats and analyze them -Copyright : (c) Dominik Schrempf 2018 -License : GPL-3 - -Maintainer : dominik.schrempf@gmail.com -Stability : unstable -Portability : portable - -Creation date: Fri Oct 5 08:41:05 2018. - --} - -module Examine.Examine - ( examineCmd - ) - where - -import Control.Monad.Logger -import Control.Monad.Trans.Class -import Control.Monad.Trans.Reader -import qualified Data.ByteString.Lazy.Char8 as L -import qualified Data.Set as S -import Text.Printf - -import Examine.Options -import Tools - -import qualified ELynx.Data.Alphabet.Alphabet as A -import qualified ELynx.Data.Alphabet.Character as C -import qualified ELynx.Data.Sequence.Alignment as M -import qualified ELynx.Data.Sequence.Sequence as Seq -import ELynx.Tools.InputOutput - -examineAlignment :: Bool -> M.Alignment -> L.ByteString -examineAlignment perSiteFlag a = - L.unlines [ L.pack "Sequences have equal length (multi sequence alignment, or single sequence)." - , L.pack $ "Total number of columns in alignment: " - ++ show (M.length a) - , L.pack $ "Number of columns without gaps: " - ++ show (M.length aNoGaps) - , L.pack $ "Number of columns with standard characters only: " - ++ show (M.length aOnlyStd) - , L.empty - , L.pack $ "Total number of characters: " ++ show nTot - , L.pack $ "Standard (i.e., not extended IUPAC) characters: " - ++ show (nTot - nIUPAC - nGaps - nUnknowns) - , L.pack $ "Extended IUPAC characters: " ++ show nIUPAC - , L.pack $ "Gaps: " ++ show nGaps - , L.pack $ "Unknowns: " ++ show nUnknowns - , L.pack $ "Percentage of standard characters: " - ++ printf "%.3f" (1.0 - percentageIUPAC - percentageGaps - percentageUnknowns) - , L.pack $ "Percentage of extended IUPAC characters: " - ++ printf "%.3f" percentageIUPAC - , L.pack $ "Percentage of gaps: " - ++ printf "%.3f" percentageGaps - , L.pack $ "Percentage of unknowns: " - ++ printf "%.3f" percentageUnknowns - , L.empty - , L.pack "Distribution of characters:" - -- Holy crap. - , L.pack $ concatMap ((: " ") . C.toChar) $ S.toList $ - A.std $ A.alphabetSpec $ M.alphabet a - , L.pack $ unwords $ map (printf "%.3f") charFreqs - , L.empty - , L.pack "Mean effective number of states (measured using entropy):" - , L.pack "Across whole alignment: " - <> L.pack (printf "%.3f" kEffMean) - , L.pack "Across columns without gaps: " - <> L.pack (printf "%.3f" kEffMeanNoGaps) - , L.pack "Across columns without extended IUPAC characters: " - <> L.pack (printf "%.3f" kEffMeanOnlyStd) - ] - <> perSiteBS - where - nTot = M.length a * M.nSequences a - nIUPAC = M.countIUPACChars a - nGaps = M.countGaps a - nUnknowns = M.countUnknowns a - percentageIUPAC = fromIntegral nIUPAC / fromIntegral nTot :: Double - percentageGaps = fromIntegral nGaps / fromIntegral nTot :: Double - percentageUnknowns = fromIntegral nUnknowns / fromIntegral nTot :: Double - aNoGaps = M.filterColsNoGaps a - aOnlyStd = M.filterColsOnlyStd aNoGaps - charFreqsPerSite = M.toFrequencyData a - charFreqs = M.distribution charFreqsPerSite - kEffs = M.kEffEntropy charFreqsPerSite - kEffsNoGaps = M.kEffEntropy . M.toFrequencyData $ aNoGaps - kEffsOnlyStd = M.kEffEntropy . M.toFrequencyData $ aOnlyStd - kEffMean = sum kEffs / fromIntegral (length kEffs) - kEffMeanNoGaps = sum kEffsNoGaps / fromIntegral (length kEffsNoGaps) - kEffMeanOnlyStd = sum kEffsOnlyStd / fromIntegral (length kEffsOnlyStd) - perSiteBS = if perSiteFlag - then L.unlines [ L.pack "Effective number of used states per site:" - , L.pack . show $ kEffs - ] - else L.empty - -examine :: Bool -> [Seq.Sequence] -> L.ByteString -examine perSiteFlag ss = Seq.summarizeSequences ss <> - case M.fromSequences ss of - Left _ -> L.empty - Right a -> L.pack "\n" <> examineAlignment perSiteFlag a - --- | Examine sequences. -examineCmd :: Maybe FilePath -> Examine () -examineCmd outFileBaseName = do - $(logInfo) "Command: Examine sequences." - ExamineArguments al inFile perSiteFlag <- lift ask - ss <- readSeqs al inFile - let result = examine perSiteFlag ss - let outFilePath = (++ ".out") <$> outFileBaseName - out "result of examination" result outFilePath diff --git a/src/Examine/Options.hs b/src/Examine/Options.hs deleted file mode 100644 --- a/src/Examine/Options.hs +++ /dev/null @@ -1,63 +0,0 @@ -{- | -Module : Examine.Options -Description : ELynxSeq argument parsing -Copyright : (c) Dominik Schrempf 2018 -License : GPL-3 - -Maintainer : dominik.schrempf@gmail.com -Stability : unstable -Portability : portable - -Creation date: Sun Oct 7 17:29:45 2018. - --} - -module Examine.Options - ( ExamineArguments (..) - , Examine - , examineArguments - ) where - -import Control.Applicative (optional) -import Control.Monad.Logger -import Control.Monad.Trans.Reader -import Data.Maybe (fromJust) -import Options.Applicative - -import Tools - -import ELynx.Data.Alphabet.Alphabet -import ELynx.Tools.Reproduction - --- | Arguments needed to examine sequences. -data ExamineArguments = ExamineArguments - { exAlphabet :: Alphabet - , exInFile :: Maybe FilePath - , exPerSite :: Bool } - -instance Reproducible ExamineArguments where - -- TODO: How do I combine optional input files with reproducibility? The - -- answer is: use checksums! Also other input files have to be checked with - -- check sums. StdIO can also be checked in the same way. - inFiles = pure . fromJust . exInFile - parser _ = examineArguments - --- | Logger and Reader type. -type Examine = LoggingT (ReaderT ExamineArguments IO) - --- | Command line parser. -examineArguments :: Parser ExamineArguments -examineArguments = ExamineArguments - <$> alphabetOpt - <*> optional filePathArg - <*> examinePerSiteOpt - -examinePerSiteOpt :: Parser Bool -examinePerSiteOpt = switch $ - long "per-site" <> - help "Report per site summary statistics" - -filePathArg :: Parser FilePath -filePathArg = strArgument $ - metavar "INPUT-FILE" <> - help "Read sequences from INPUT-FILE" diff --git a/src/Filter/Filter.hs b/src/Filter/Filter.hs deleted file mode 100644 --- a/src/Filter/Filter.hs +++ /dev/null @@ -1,80 +0,0 @@ -{-# LANGUAGE OverloadedStrings #-} -{-# LANGUAGE TemplateHaskell #-} - -{- | -Module : Analyze.Analyze -Description : Parse sequence file formats and analyze them -Copyright : (c) Dominik Schrempf 2018 -License : GPL-3 - -Maintainer : dominik.schrempf@gmail.com -Stability : unstable -Portability : portable - -Creation date: Fri Oct 5 08:41:05 2018. - --} - -module Filter.Filter - ( filterRowsCmd - , filterColsCmd - ) - where - -import Control.Monad (when) -import Control.Monad.Logger -import Control.Monad.Trans.Class -import Control.Monad.Trans.Reader -import qualified Data.ByteString.Lazy.Char8 as L -import Data.Maybe (fromMaybe) -import qualified Data.Text as T - -import Filter.Options -import Tools - -import qualified ELynx.Data.Sequence.Alignment as M -import qualified ELynx.Data.Sequence.Sequence as S -import ELynx.Export.Sequence.Fasta -import ELynx.Tools.InputOutput -import ELynx.Tools.Misc - -filterRows :: Maybe Int -> Maybe Int -> Bool -> [S.Sequence] -> L.ByteString -filterRows ml ms std ss = sequencesToFasta $ compose filters ss - where filters' = map (fromMaybe id) [S.filterLongerThan <$> ml, S.filterShorterThan <$> ms] - filters = if std then S.filterStandard : filters' else filters' - --- | Filter sequences. -filterRowsCmd :: Maybe FilePath -> FilterRows () -filterRowsCmd outFileBaseName = do - $(logInfo) "Command: Filter sequences of a list of sequences." - FilterRowsArguments al inFile long short std <- lift ask - maybe (return ()) - (\val -> $(logInfo) $ T.pack $ " Keep sequences longer than " <> show val <> ".") long - maybe (return ()) - (\val -> $(logInfo) $ T.pack $ " Keep sequences shorter than " <> show val <> ".") short - when std $ - $(logInfo) " Keep sequences containing at least one standard (i.e., non-IUPAC) character." - ss <- readSeqs al inFile - let result = filterRows long short std ss - let outFilePath = (++ ".fasta") <$> outFileBaseName - out "filtered sequences" result outFilePath - -filterCols :: Maybe Double -> [S.Sequence] -> L.ByteString -filterCols ms ss = sequencesToFasta . M.toSequences $ compose filters a - where a = either error id (M.fromSequences ss) - filters = map (fromMaybe id) [ M.filterColsStd <$> ms ] - --- | Filter columns. -filterColsCmd :: Maybe FilePath -> FilterCols () -filterColsCmd outFileBaseName = do - $(logInfo) "Command: Filter columns of a multi sequence alignment." - FilterColsArguments al inFile standard <- lift ask - case standard of - Nothing -> return () - Just p -> $(logInfo) $ T.pack $ - " Keep columns with a proportion of standard (non-IUPAC) characters larger than " - ++ show p ++ "." - ss <- readSeqs al inFile - let result = filterCols standard ss - let outFilePath = (++ ".fasta") <$> outFileBaseName - out "filtered sequences" result outFilePath diff --git a/src/Filter/Options.hs b/src/Filter/Options.hs deleted file mode 100644 --- a/src/Filter/Options.hs +++ /dev/null @@ -1,95 +0,0 @@ -{- | -Module : Filter.Options -Description : ELynxSeq argument parsing -Copyright : (c) Dominik Schrempf 2018 -License : GPL-3 - -Maintainer : dominik.schrempf@gmail.com -Stability : unstable -Portability : portable - -Creation date: Sun Oct 7 17:29:45 2018. - --} - -module Filter.Options - ( FilterRowsArguments (..) - , FilterColsArguments (..) - , FilterRows - , FilterCols - , filterRowsArguments - , filterColsArguments - ) where - -import Control.Applicative -import Control.Monad.Logger -import Control.Monad.Trans.Reader -import Options.Applicative - -import Tools - -import ELynx.Data.Alphabet.Alphabet - --- | Arguments needed for filtering sequences. -data FilterRowsArguments = FilterRowsArguments - { frAlphabet :: Alphabet - , frInFile :: Maybe FilePath - , frLonger :: Maybe Int - , frShorter :: Maybe Int - , frStandard :: Bool } - --- | Arguments needed for filtering columns of a multi sequence alignment. -data FilterColsArguments = FilterColsArguments - { fcAlphabet :: Alphabet - , fcInFile :: Maybe FilePath - , fcStandard :: Maybe Double } - --- | Logger and Reader for sequence filtering. -type FilterRows = LoggingT (ReaderT FilterRowsArguments IO) - --- | Logger and Reader for column filtering. -type FilterCols = LoggingT (ReaderT FilterColsArguments IO) - --- | Command line parser. -filterRowsArguments :: Parser FilterRowsArguments -filterRowsArguments = FilterRowsArguments - <$> alphabetOpt - <*> optional inFileArg - <*> filterLongerThanOpt - <*> filterShorterThanOpt - <*> filterStandardChars - -filterLongerThanOpt :: Parser (Maybe Int) -filterLongerThanOpt = optional $ option auto $ - long "longer-than" <> - metavar "LENGTH" <> - help "Only keep sequences longer than LENGTH" - -filterShorterThanOpt :: Parser (Maybe Int) -filterShorterThanOpt = optional $ option auto $ - long "shorter-than" <> - metavar "LENGTH" <> - help "Only keep sequences shorter than LENGTH" - -filterStandardChars :: Parser Bool -filterStandardChars = switch $ - long "standard-characters" <> - help "Only keep sequences containing at least one standard (i.e., non-IUPAC) character" - --- | Command line parser. -filterColsArguments :: Parser FilterColsArguments -filterColsArguments = FilterColsArguments - <$> alphabetOpt - <*> optional inFileArg - <*> filterStandardOpt - -filterStandardOpt :: Parser (Maybe Double) -filterStandardOpt = optional $ option auto $ - long "standard-chars" <> - metavar "DOUBLE" <> - help "Keep columns with a proportion standard (non-IUPAC) characters larger than DOUBLE in [0,1]" - -inFileArg :: Parser FilePath -inFileArg = strArgument $ - metavar "INPUT-FILE" <> - help "Read sequences from INPUT-FILE" diff --git a/src/SLynx/Concatenate/Concatenate.hs b/src/SLynx/Concatenate/Concatenate.hs new file mode 100644 --- /dev/null +++ b/src/SLynx/Concatenate/Concatenate.hs @@ -0,0 +1,40 @@ +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE TemplateHaskell #-} + +{- | +Module : Analyze.Analyze +Description : Parse sequence file formats and analyze them +Copyright : (c) Dominik Schrempf 2018 +License : GPL-3.0-or-later + +Maintainer : dominik.schrempf@gmail.com +Stability : unstable +Portability : portable + +Creation date: Fri Oct 5 08:41:05 2018. + +-} + +module SLynx.Concatenate.Concatenate + ( concatenateCmd + ) +where + +import Control.Monad.Logger +import Control.Monad.Trans.Reader ( ask ) + +import SLynx.Concatenate.Options +import SLynx.Tools + +import qualified ELynx.Data.Sequence.Sequence as S +import ELynx.Export.Sequence.Fasta +import ELynx.Tools + +-- | Concatenate sequences. +concatenateCmd :: ELynx ConcatenateArguments () +concatenateCmd = do + (ConcatenateArguments al fps) <- local <$> ask + $(logInfo) "Command: Concatenate sequences." + sss <- mapM (readSeqs al) fps + let result = sequencesToFasta $ S.concatSequences sss + out "concatenated multi sequence alignment " result ".fasta" diff --git a/src/SLynx/Concatenate/Options.hs b/src/SLynx/Concatenate/Options.hs new file mode 100644 --- /dev/null +++ b/src/SLynx/Concatenate/Options.hs @@ -0,0 +1,56 @@ +{-# LANGUAGE DeriveGeneric #-} + +{- | +Module : SLynx.Concatenate.Options +Description : ELynxSeq argument parsing +Copyright : (c) Dominik Schrempf 2018 +License : GPL-3.0-or-later + +Maintainer : dominik.schrempf@gmail.com +Stability : unstable +Portability : portable + +Creation date: Sun Oct 7 17:29:45 2018. + +-} + +module SLynx.Concatenate.Options + ( ConcatenateArguments(..) + , concatenateArguments + ) +where + +import Control.Applicative +import Options.Applicative + +import SLynx.Tools + +import ELynx.Data.Alphabet.Alphabet +import ELynx.Tools + +-- | Arguments needed to concatenate multi sequence alignments. +data ConcatenateArguments = ConcatenateArguments + { ccAlphabet :: Alphabet + , ccInFiles :: [FilePath] } + deriving (Eq, Show, Generic) + +instance Reproducible ConcatenateArguments where + inFiles = ccInFiles + outSuffixes _ = [".fasta"] + getSeed _ = Nothing + setSeed = const + parser = concatenateArguments + cmdName = "concatenate" + cmdDsc = ["Concatenate sequences found in input files."] + +instance FromJSON ConcatenateArguments + +instance ToJSON ConcatenateArguments + +-- | Command line parser. +concatenateArguments :: Parser ConcatenateArguments +concatenateArguments = ConcatenateArguments <$> alphabetOpt <*> some inFileArg + +inFileArg :: Parser FilePath +inFileArg = + strArgument $ metavar "INPUT-FILE" <> help "Read sequences from INPUT-FILE" diff --git a/src/SLynx/Examine/Examine.hs b/src/SLynx/Examine/Examine.hs new file mode 100644 --- /dev/null +++ b/src/SLynx/Examine/Examine.hs @@ -0,0 +1,138 @@ +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE TemplateHaskell #-} + +{- | + +Module : Analyze.Analyze +Description : Parse sequence file formats and analyze them +Copyright : (c) Dominik Schrempf 2018 +License : GPL-3.0-or-later + +Maintainer : dominik.schrempf@gmail.com +Stability : unstable +Portability : portable + +Creation date: Fri Oct 5 08:41:05 2018. + +-} + +module SLynx.Examine.Examine + ( examineCmd + ) +where + +import Control.Monad.Logger +import Control.Monad.Trans.Reader ( ask ) +import qualified Data.ByteString.Lazy.Char8 as L +import qualified Data.Set as S +import Text.Printf + +import SLynx.Examine.Options +import SLynx.Tools + +import qualified ELynx.Data.Alphabet.Alphabet as A +import qualified ELynx.Data.Alphabet.Character as C +import qualified ELynx.Data.Sequence.Alignment as M +import qualified ELynx.Data.Sequence.Sequence as Seq +import ELynx.Tools + +pRow :: String -> String -> L.ByteString +pRow name val = alignLeft 50 n <> alignRight 10 v + where + n = L.pack name + v = L.pack val + +examineAlignment :: Bool -> M.Alignment -> L.ByteString +examineAlignment perSiteFlag a = + L.unlines + [ L.pack + "Sequences have equal length (multi sequence alignment, or single sequence)." + , pRow "Total number of columns in alignment:" $ show (M.length a) + , pRow "Number of columns without gaps:" $ show (M.length aNoGaps) + , pRow "Number of columns with standard characters only:" + $ show (M.length aOnlyStd) + , L.empty + , pRow "Total number of characters:" $ show nTot + , pRow "Standard (i.e., not extended IUPAC) characters:" + $ show (nTot - nIUPAC - nGaps - nUnknowns) + , pRow "Extended IUPAC characters:" $ show nIUPAC + , pRow "Gaps:" $ show nGaps + , pRow "Unknowns:" $ show nUnknowns + , L.empty + , pRow "Percentage of standard characters:" + $ printf "%2.2f" (100.0 - percentIUPAC - percentGaps - percentUnknowns) + , pRow "Percentage of extended IUPAC characters:" + $ printf "%2.2f" percentIUPAC + , pRow "Percentage of gaps:" $ printf "%2.2f" percentGaps + , pRow "Percentage of unknowns:" $ printf "%2.2f" percentUnknowns + , L.empty + , L.pack "Distribution of characters:" + -- Holy crap. + , L.pack + $ concatMap ((: " ") . C.toChar) + $ S.toList + $ A.std + $ A.alphabetSpec + $ M.alphabet a + , L.pack $ unwords $ map (printf "%.3f") charFreqs + , L.empty + , L.pack "Mean effective number of states (measured using entropy):" + , pRow "Across whole alignment:" $ printf "%.3f" kEffMean + , pRow "Across columns without gaps:" $ printf "%.3f" kEffMeanNoGaps + , pRow "Across columns without extended IUPAC characters:" + $ printf "%.3f" kEffMeanOnlyStd + , L.empty + , L.pack "Mean effective number of states (measured using homoplasy):" + , pRow "Across whole alignment:" $ printf "%.3f" kEffMeanHomo + , pRow "Across columns without gaps:" $ printf "%.3f" kEffMeanNoGapsHomo + , pRow "Across columns without extended IUPAC characters:" + $ printf "%.3f" kEffMeanOnlyStdHomo + ] + <> perSiteBS + where + nTot = M.length a * M.nSequences a + nIUPAC = M.countIUPACChars a + nGaps = M.countGaps a + nUnknowns = M.countUnknowns a + percentIUPAC = 100 * fromIntegral nIUPAC / fromIntegral nTot :: Double + percentGaps = 100 * fromIntegral nGaps / fromIntegral nTot :: Double + percentUnknowns = 100 * fromIntegral nUnknowns / fromIntegral nTot :: Double + aNoGaps = M.filterColsNoGaps a + aOnlyStd = M.filterColsOnlyStd aNoGaps + charFreqsPerSite = M.toFrequencyData a + charFreqs = M.distribution charFreqsPerSite + kEffs = M.kEffEntropy charFreqsPerSite + kEffsNoGaps = M.kEffEntropy . M.toFrequencyData $ aNoGaps + kEffsOnlyStd = M.kEffEntropy . M.toFrequencyData $ aOnlyStd + kEffMean = sum kEffs / fromIntegral (length kEffs) + kEffMeanNoGaps = sum kEffsNoGaps / fromIntegral (length kEffsNoGaps) + kEffMeanOnlyStd = sum kEffsOnlyStd / fromIntegral (length kEffsOnlyStd) + kEffsHomo = M.kEffHomoplasy charFreqsPerSite + kEffsNoGapsHomo = M.kEffHomoplasy . M.toFrequencyData $ aNoGaps + kEffsOnlyStdHomo = M.kEffHomoplasy . M.toFrequencyData $ aOnlyStd + kEffMeanHomo = sum kEffsHomo / fromIntegral (length kEffsHomo) + kEffMeanNoGapsHomo = + sum kEffsNoGapsHomo / fromIntegral (length kEffsNoGapsHomo) + kEffMeanOnlyStdHomo = + sum kEffsOnlyStdHomo / fromIntegral (length kEffsOnlyStdHomo) + perSiteBS = if perSiteFlag + then L.unlines + [ L.pack "Effective number of used states per site:" + , L.pack . show $ kEffs + ] + else L.empty + +examine :: Bool -> [Seq.Sequence] -> L.ByteString +examine perSiteFlag ss = + Seq.summarizeSequences ss <> case M.fromSequences ss of + Left _ -> L.empty + Right a -> L.pack "\n" <> examineAlignment perSiteFlag a + +-- | Examine sequences. +examineCmd :: ELynx ExamineArguments () +examineCmd = do + (ExamineArguments al inFile perSiteFlag) <- local <$> ask + $(logInfo) "Command: Examine sequences." + ss <- readSeqs al inFile + let result = examine perSiteFlag ss + out "result of examination" result ".out" diff --git a/src/SLynx/Examine/Options.hs b/src/SLynx/Examine/Options.hs new file mode 100644 --- /dev/null +++ b/src/SLynx/Examine/Options.hs @@ -0,0 +1,64 @@ +{-# LANGUAGE DeriveGeneric #-} + +{- | +Module : SLynx.Examine.Options +Description : ELynxSeq argument parsing +Copyright : (c) Dominik Schrempf 2018 +License : GPL-3.0-or-later + +Maintainer : dominik.schrempf@gmail.com +Stability : unstable +Portability : portable + +Creation date: Sun Oct 7 17:29:45 2018. + +-} + +module SLynx.Examine.Options + ( ExamineArguments(..) + , examineArguments + ) +where + +import Options.Applicative + +import SLynx.Tools + +import ELynx.Data.Alphabet.Alphabet +import ELynx.Tools + +-- | Arguments needed to examine sequences. +data ExamineArguments = ExamineArguments + { exAlphabet :: Alphabet + , exInFile :: FilePath + , exPerSite :: Bool } + deriving (Eq, Show, Generic) + +instance Reproducible ExamineArguments where + inFiles = pure . exInFile + outSuffixes _ = [".out"] + getSeed _ = Nothing +-- XXX: Probably throw error when seed is set. + setSeed = const + parser = examineArguments + cmdName = "examine" + cmdDsc = + [ "Examine sequences. If data is a multi sequence alignment, additionally analyze columns." + ] + +instance FromJSON ExamineArguments + +instance ToJSON ExamineArguments + +-- | Command line parser. +examineArguments :: Parser ExamineArguments +examineArguments = + ExamineArguments <$> alphabetOpt <*> filePathArg <*> examinePerSiteOpt + +examinePerSiteOpt :: Parser Bool +examinePerSiteOpt = + switch $ long "per-site" <> help "Report per site summary statistics" + +filePathArg :: Parser FilePath +filePathArg = + strArgument $ metavar "INPUT-FILE" <> help "Read sequences from INPUT-FILE" diff --git a/src/SLynx/Filter/Filter.hs b/src/SLynx/Filter/Filter.hs new file mode 100644 --- /dev/null +++ b/src/SLynx/Filter/Filter.hs @@ -0,0 +1,91 @@ +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE TemplateHaskell #-} + +{- | +Module : Analyze.Analyze +Description : Parse sequence file formats and analyze them +Copyright : (c) Dominik Schrempf 2018 +License : GPL-3.0-or-later + +Maintainer : dominik.schrempf@gmail.com +Stability : unstable +Portability : portable + +Creation date: Fri Oct 5 08:41:05 2018. + +-} + +module SLynx.Filter.Filter + ( filterRowsCmd + , filterColsCmd + ) +where + +import Control.Monad ( when ) +import Control.Monad.Logger +import Control.Monad.Trans.Reader ( ask ) +import qualified Data.ByteString.Lazy.Char8 as L +import Data.Maybe ( fromMaybe ) +import qualified Data.Text as T + +import SLynx.Filter.Options +import SLynx.Tools + +import qualified ELynx.Data.Sequence.Alignment as M +import qualified ELynx.Data.Sequence.Sequence as S +import ELynx.Export.Sequence.Fasta +import ELynx.Tools + +filterRows :: Maybe Int -> Maybe Int -> Bool -> [S.Sequence] -> L.ByteString +filterRows ml ms std ss = sequencesToFasta $ compose filters ss + where + filters' = + map (fromMaybe id) [S.filterLongerThan <$> ml, S.filterShorterThan <$> ms] + filters = if std then S.filterStandard : filters' else filters' + +-- | Filter sequences. +filterRowsCmd :: ELynx FilterRowsArguments () +filterRowsCmd = do + (FilterRowsArguments al inFile long short std) <- local <$> ask + $(logInfo) "Command: Filter sequences of a list of sequences." + maybe + (return ()) + (\val -> + $(logInfo) $ T.pack $ " Keep sequences longer than " <> show val <> "." + ) + long + maybe + (return ()) + (\val -> + $(logInfo) $ T.pack $ " Keep sequences shorter than " <> show val <> "." + ) + short + when std + $ $(logInfo) + " Keep sequences containing at least one standard (i.e., non-IUPAC) character." + ss <- readSeqs al inFile + let result = filterRows long short std ss + out "filtered sequences" result ".fasta" + +filterCols :: Maybe Double -> [S.Sequence] -> L.ByteString +filterCols ms ss = sequencesToFasta . M.toSequences $ compose filters a + where + a = either error id (M.fromSequences ss) + filters = map (fromMaybe id) [M.filterColsStd <$> ms] + +-- | Filter columns. +filterColsCmd :: ELynx FilterColsArguments () +filterColsCmd = do + (FilterColsArguments al inFile standard) <- local <$> ask + $(logInfo) "Command: Filter columns of a multi sequence alignment." + case standard of + Nothing -> return () + Just p -> + $(logInfo) + $ T.pack + $ " Keep columns with a proportion of standard (non-IUPAC) characters larger than " + ++ show p + ++ "." + ss <- readSeqs al inFile + let result = filterCols standard ss + out "filtered sequences" result ".fasta" diff --git a/src/SLynx/Filter/Options.hs b/src/SLynx/Filter/Options.hs new file mode 100644 --- /dev/null +++ b/src/SLynx/Filter/Options.hs @@ -0,0 +1,118 @@ +{-# LANGUAGE DeriveGeneric #-} + +{- | +Module : SLynx.Filter.Options +Description : ELynxSeq argument parsing +Copyright : (c) Dominik Schrempf 2018 +License : GPL-3.0-or-later + +Maintainer : dominik.schrempf@gmail.com +Stability : unstable +Portability : portable + +Creation date: Sun Oct 7 17:29:45 2018. + +-} + +module SLynx.Filter.Options + ( FilterRowsArguments(..) + , FilterColsArguments(..) + , filterRowsArguments + , filterColsArguments + ) +where + +import Control.Applicative +import Options.Applicative + +import SLynx.Tools + +import ELynx.Data.Alphabet.Alphabet +import ELynx.Tools + +-- | Arguments needed for filtering sequences. +data FilterRowsArguments = FilterRowsArguments + { frAlphabet :: Alphabet + , frInFile :: FilePath + , frLonger :: Maybe Int + , frShorter :: Maybe Int + , frStandard :: Bool } + deriving (Eq, Show, Generic) + +instance Reproducible FilterRowsArguments where + inFiles = pure . frInFile + outSuffixes _ = [".fasta"] + getSeed _ = Nothing + setSeed = const + parser = filterRowsArguments + cmdName = "filter-rows" + cmdDsc = ["Filter rows (or sequences) found in input files."] + +instance FromJSON FilterRowsArguments + +instance ToJSON FilterRowsArguments + +-- | Arguments needed for filtering columns of a multi sequence alignment. +data FilterColsArguments = FilterColsArguments + { fcAlphabet :: Alphabet + , fcInFile :: FilePath + , fcStandard :: Maybe Double } + deriving (Eq, Show, Generic) + +instance Reproducible FilterColsArguments where + inFiles = pure . fcInFile + outSuffixes _ = [".fasta"] + getSeed _ = Nothing + setSeed = const + parser = filterColsArguments + cmdName = "filter-columns" + cmdDsc = ["Filter columns of multi sequence alignments."] + +instance FromJSON FilterColsArguments + +instance ToJSON FilterColsArguments + +-- | Command line parser. +filterRowsArguments :: Parser FilterRowsArguments +filterRowsArguments = + FilterRowsArguments + <$> alphabetOpt + <*> inFileArg + <*> filterLongerThanOpt + <*> filterShorterThanOpt + <*> filterStandardChars + +filterLongerThanOpt :: Parser (Maybe Int) +filterLongerThanOpt = + optional $ option auto $ long "longer-than" <> metavar "LENGTH" <> help + "Only keep sequences longer than LENGTH" + +filterShorterThanOpt :: Parser (Maybe Int) +filterShorterThanOpt = + optional $ option auto $ long "shorter-than" <> metavar "LENGTH" <> help + "Only keep sequences shorter than LENGTH" + +filterStandardChars :: Parser Bool +filterStandardChars = + switch + $ long "standard-characters" + <> help + "Only keep sequences containing at least one standard (i.e., non-IUPAC) character" + +-- | Command line parser. +filterColsArguments :: Parser FilterColsArguments +filterColsArguments = + FilterColsArguments <$> alphabetOpt <*> inFileArg <*> filterStandardOpt + +filterStandardOpt :: Parser (Maybe Double) +filterStandardOpt = + optional + $ option auto + $ long "standard-chars" + <> metavar "DOUBLE" + <> help + "Keep columns with a proportion standard (non-IUPAC) characters larger than DOUBLE in [0,1]" + +inFileArg :: Parser FilePath +inFileArg = + strArgument $ metavar "INPUT-FILE" <> help "Read sequences from INPUT-FILE" diff --git a/src/SLynx/Options.hs b/src/SLynx/Options.hs new file mode 100644 --- /dev/null +++ b/src/SLynx/Options.hs @@ -0,0 +1,130 @@ +{-# LANGUAGE DeriveGeneric #-} + +{- | +Module : SLynx.Options +Description : SLynx general options +Copyright : (c) Dominik Schrempf 2020 +License : GPL-3.0-or-later + +Maintainer : dominik.schrempf@gmail.com +Stability : unstable +Portability : portable + +Creation date: Sat Sep 7 18:55:03 2019. + +-} + +module SLynx.Options + ( Arguments(..) + , CommandArguments(..) + , parseArguments + ) +where + +import Options.Applicative + +import ELynx.Data.Alphabet.Alphabet +import ELynx.Tools + +import SLynx.Concatenate.Options +import SLynx.Examine.Options +import SLynx.Filter.Options +import SLynx.Simulate.Options +import SLynx.SubSample.Options +import SLynx.Translate.Options + +-- | The different SLynx commands and their arguments. +data CommandArguments = + Concatenate ConcatenateArguments + | Examine ExamineArguments + | FilterCols FilterColsArguments + | FilterRows FilterRowsArguments + | Simulate SimulateArguments + | SubSample SubSampleArguments + | Translate TranslateArguments + deriving (Eq, Show, Generic) + +instance Reproducible CommandArguments where + inFiles (Concatenate a) = inFiles a + inFiles (Examine a) = inFiles a + inFiles (FilterCols a) = inFiles a + inFiles (FilterRows a) = inFiles a + inFiles (Simulate a) = inFiles a + inFiles (SubSample a) = inFiles a + inFiles (Translate a) = inFiles a + + outSuffixes (Concatenate a) = outSuffixes a + outSuffixes (Examine a) = outSuffixes a + outSuffixes (FilterCols a) = outSuffixes a + outSuffixes (FilterRows a) = outSuffixes a + outSuffixes (Simulate a) = outSuffixes a + outSuffixes (SubSample a) = outSuffixes a + outSuffixes (Translate a) = outSuffixes a + + getSeed (Concatenate a) = getSeed a + getSeed (Examine a) = getSeed a + getSeed (FilterCols a) = getSeed a + getSeed (FilterRows a) = getSeed a + getSeed (Simulate a) = getSeed a + getSeed (SubSample a) = getSeed a + getSeed (Translate a) = getSeed a + + setSeed (Concatenate a) = Concatenate . setSeed a + setSeed (Examine a) = Examine . setSeed a + setSeed (FilterCols a) = FilterCols . setSeed a + setSeed (FilterRows a) = FilterRows . setSeed a + setSeed (Simulate a) = Simulate . setSeed a + setSeed (SubSample a) = SubSample . setSeed a + setSeed (Translate a) = Translate . setSeed a + + parser = commandArguments + + cmdName = "slynx" + + cmdDsc = ["Analyze, and simulate multi sequence alignments."] + + cmdFtr = + ["Available sequence file formats:"] + ++ fs + ++ ["", "Available alphabets:"] + ++ as + where + toListItem = (" - " ++) + fs = map toListItem ["FASTA"] + as = map (toListItem . alphabetDescription) (allValues :: [Alphabet]) + +instance FromJSON CommandArguments + +instance ToJSON CommandArguments + +concatenateCommand :: Mod CommandFields CommandArguments +concatenateCommand = createCommandReproducible Concatenate + +examineCommand :: Mod CommandFields CommandArguments +examineCommand = createCommandReproducible Examine + +filterColumnsCommand :: Mod CommandFields CommandArguments +filterColumnsCommand = createCommandReproducible FilterCols + +filterRowsCommand :: Mod CommandFields CommandArguments +filterRowsCommand = createCommandReproducible FilterRows + +simulateCommand :: Mod CommandFields CommandArguments +simulateCommand = createCommandReproducible Simulate + +subSampleCommand :: Mod CommandFields CommandArguments +subSampleCommand = createCommandReproducible SubSample + +translateCommand :: Mod CommandFields CommandArguments +translateCommand = createCommandReproducible Translate + +commandArguments :: Parser CommandArguments +commandArguments = + hsubparser + $ concatenateCommand + <> examineCommand + <> filterColumnsCommand + <> filterRowsCommand + <> simulateCommand + <> subSampleCommand + <> translateCommand diff --git a/src/SLynx/SLynx.hs b/src/SLynx/SLynx.hs new file mode 100644 --- /dev/null +++ b/src/SLynx/SLynx.hs @@ -0,0 +1,61 @@ +{- | +Module : SLynx.SLynx +Description : SLynx module +Copyright : (c) Dominik Schrempf 2020 +License : GPL-3.0-or-later + +Maintainer : dominik.schrempf@gmail.com +Stability : unstable +Portability : portable + +Creation date: Thu Apr 23 16:38:55 2020. + +-} + +module SLynx.SLynx + ( slynx + , rSLynx + ) +where + +import SLynx.Options + +import SLynx.Concatenate.Concatenate +import SLynx.Examine.Examine +import SLynx.Filter.Filter +import SLynx.Simulate.Simulate +import SLynx.SubSample.SubSample +import SLynx.Translate.Translate + +import ELynx.Tools + +-- TODO: Use a class here (e.g., elynx-wrappable) which defines the extractor function. +-- | Run SLynx with given arguments. +slynx :: Arguments CommandArguments -> IO () +slynx c = case local c of + Concatenate _ -> eLynxWrapper + c + (\(Arguments g (Concatenate l)) -> Arguments g l) + concatenateCmd + Examine _ -> + eLynxWrapper c (\(Arguments g (Examine l)) -> Arguments g l) examineCmd + FilterCols _ -> eLynxWrapper + c + (\(Arguments g (FilterCols l)) -> Arguments g l) + filterColsCmd + FilterRows _ -> eLynxWrapper + c + (\(Arguments g (FilterRows l)) -> Arguments g l) + filterRowsCmd + Simulate _ -> + eLynxWrapper c (\(Arguments g (Simulate l)) -> Arguments g l) simulateCmd + SubSample _ -> eLynxWrapper c + (\(Arguments g (SubSample l)) -> Arguments g l) + subSampleCmd + Translate _ -> eLynxWrapper c + (\(Arguments g (Translate l)) -> Arguments g l) + translateCmd + +-- | Run SLynx, parse arguments from command line. +rSLynx :: IO () +rSLynx = parseArguments >>= slynx diff --git a/src/SLynx/Simulate/Options.hs b/src/SLynx/Simulate/Options.hs new file mode 100644 --- /dev/null +++ b/src/SLynx/Simulate/Options.hs @@ -0,0 +1,212 @@ +{-# LANGUAGE DeriveGeneric #-} + +{- | +Module : SLynx.Simulate.Options +Description : ELynxSim argument parsing +Copyright : (c) Dominik Schrempf 2018 +License : GPL-3.0-or-later + +Maintainer : dominik.schrempf@gmail.com +Stability : unstable +Portability : portable + +Creation date: Sun Oct 7 17:29:45 2018. + +Available options: + -h,--help Show this help text + -v,--version Show version + -t,--tree-file NAME Specify tree file NAME + -s,--substitution-model MODEL + Set the phylogenetic substitution model; available + models are shown below + -m,--mixture-model MODEL Set the phylogenetic mixture model; available models + are shown below + -l,--length NUMBER Set alignment length to NUMBER + -e,--edm-file NAME empirical distribution model file NAME in Phylobayes + format + -w,--mixture-model-weights [DOUBLE,DOUBLE,...] + weights of mixture model components + -g,--gamma-rate-heterogeneity (NCAT, SHAPE) + number of gamma rate categories and shape parameter + -e,--seed [INT] Set seed for the random number generator; list of 32 + bit integers with up to 256 elements (default: [0]) + -q,--quiet Be quiet + -o,--output-file NAME Specify output file NAME + + +-} + + +module SLynx.Simulate.Options + ( GammaRateHeterogeneityParams + , SimulateArguments(..) + , simulateArguments + , simulateFooter + ) +where + +import Data.Maybe ( maybeToList + , fromMaybe + ) +import Options.Applicative + +import ELynx.Tools + +-- | Number of gamma rate categories and alpha parameter. +type GammaRateHeterogeneityParams = (Int, Double) + +-- | Arguments needed to simulate sequences. +data SimulateArguments = SimulateArguments + { argsTreeFile :: FilePath + , argsSubstitutionModelString :: Maybe String + , argsMixtureModelString :: Maybe String + , argsEDMFile :: Maybe FilePath + , argsSiteprofilesFiles :: Maybe [FilePath] + , argsMixtureWeights :: Maybe [Double] + , argsGammaParams :: Maybe GammaRateHeterogeneityParams + , argsLength :: Int + , argsSeed :: Seed + } + deriving (Eq, Show, Generic) + +instance Reproducible SimulateArguments where + inFiles a = + argsTreeFile a + : (maybeToList (argsEDMFile a) ++ fromMaybe [] (argsSiteprofilesFiles a)) + outSuffixes _ = [".model.gz", ".fasta"] + getSeed = Just . argsSeed + setSeed a s = a { argsSeed = Fixed s } + parser = simulateArguments + cmdName = "simulate" + cmdDsc = ["Simulate multi sequence alignments."] + cmdFtr = simulateFooter + +instance FromJSON SimulateArguments + +instance ToJSON SimulateArguments + +-- | Sub command parser. +simulateArguments :: Parser SimulateArguments +simulateArguments = + SimulateArguments + <$> treeFileOpt + <*> phyloSubstitutionModelOpt + <*> phyloMixtureModelOpt + <*> maybeEDMFileOpt + <*> maybeSiteprofilesFilesOpt + <*> maybeMixtureWeights + <*> maybeGammaParams + <*> lengthOpt + <*> seedOpt + +treeFileOpt :: Parser FilePath +treeFileOpt = + strOption $ long "tree-file" <> short 't' <> metavar "Name" <> help + "Read tree from Newick file NAME" + +phyloSubstitutionModelOpt :: Parser (Maybe String) +phyloSubstitutionModelOpt = + optional + $ strOption + $ long "substitution-model" + <> short 's' + <> metavar "MODEL" + <> help + "Set the phylogenetic substitution model; available models are shown below (mutually exclusive with -m option)" + +phyloMixtureModelOpt :: Parser (Maybe String) +phyloMixtureModelOpt = optional $ strOption + ( long "mixture-model" + <> short 'm' + <> metavar "MODEL" + <> help + "Set the phylogenetic mixture model; available models are shown below (mutually exclusive with -s option)" + ) + +maybeEDMFileOpt :: Parser (Maybe FilePath) +maybeEDMFileOpt = optional $ strOption + (long "edm-file" <> short 'e' <> metavar "NAME" <> help + "Empirical distribution model file NAME in Phylobayes format" + ) + +-- fn :: Parsec Void String FilePath +-- fn = space *> takeWhile1P () <* space + +-- fns :: Parsec Void String [FilePath] +-- fns = do +-- many string + +maybeSiteprofilesFilesOpt :: Parser (Maybe [FilePath]) +maybeSiteprofilesFilesOpt = optional $ words <$> strOption + (long "siteprofile-files" <> short 'p' <> metavar "NAMES" <> help + "File names of site profiles in Phylobayes format" + ) + +maybeMixtureWeights :: Parser (Maybe [Double]) +maybeMixtureWeights = optional $ option + auto + ( long "mixture-model-weights" + <> short 'w' + <> metavar "\"[DOUBLE,DOUBLE,...]\"" + <> help "Weights of mixture model components" + ) + +maybeGammaParams :: Parser (Maybe GammaRateHeterogeneityParams) +maybeGammaParams = optional $ option + auto + ( long "gamma-rate-heterogeneity" + <> short 'g' + <> metavar "\"(NCAT,SHAPE)\"" + <> help "Number of gamma rate categories and shape parameter" + ) + +lengthOpt :: Parser Int +lengthOpt = option + auto + (long "length" <> short 'l' <> metavar "NUMBER" <> help + "Set alignment length to NUMBER" + ) + +-- | The model specification is somewhat complicated, so we need to provide +-- additional help. +simulateFooter :: [String] +simulateFooter = sms ++ mms + where + sms = + [ "Substitution models:" + , "-s \"MODEL[PARAMETER,PARAMETER,...]{STATIONARY_DISTRIBUTION}\"" + , " Supported DNA models: JC, F81, HKY, GTR4." + , " For example," + , " -s HKY[KAPPA]{DOUBLE,DOUBLE,DOUBLE,DOUBLE}" + , " -s GTR4[e_AC,e_AG,e_AT,e_CG,e_CT,e_GT]{DOUBLE,DOUBLE,DOUBLE,DOUBLE}" + , " where the 'e_XY' are the exchangeabilities from nucleotide X to Y." + , " Supported Protein models: Poisson, Poisson-Custom, LG, LG-Custom, WAG, WAG-Custom, GTR20." + , " MODEL-Custom means that only the exchangeabilities of MODEL are used," + , " and a custom stationary distribution is provided." + , " For example," + , " -s LG" + , " -s LG-Custom{...}" + , " -s GTR20[e_AR,e_AN,...]{...}" + , " the 'e_XY' are the exchangeabilities from amino acid X to Y (alphabetical order)." + , " Notes: The F81 model for DNA is equivalent to the Poisson-Custom for proteins." + , " The GTR4 model for DNA is equivalent to the GTR20 for proteins." + ] + mms = + [ "" + , "Mixture models:" + , "-m \"MIXTURE(SUBSTITUTION_MODEL_1,SUBSTITUTION_MODEL_2[PARAMETERS]{STATIONARY_DISTRIBUTION},...)\"" + , " For example," + , " -m \"MIXTURE(JC,HKY[6.0]{0.3,0.2,0.2,0.3})\"" + , "Mixture weights have to be provided with the -w option." + , "" + , "Special mixture models:" + , "-m CXX" + , " where XX is 10, 20, 30, 40, 50, or 60; CXX models, Quang et al., 2008." + , "-m \"EDM(EXCHANGEABILITIES)\"" + , " Arbitrary empirical distribution mixture (EDM) models." + , " Stationary distributions have to be provided with the -e option." + , " For example," + , " LG exchangeabilities with stationary distributions given in FILE." + , " -m \"EDM(LG-Custom)\" -e FILE" + , "For special mixture models, mixture weights are optional." + ] diff --git a/src/SLynx/Simulate/PhyloModel.hs b/src/SLynx/Simulate/PhyloModel.hs new file mode 100644 --- /dev/null +++ b/src/SLynx/Simulate/PhyloModel.hs @@ -0,0 +1,238 @@ +{- | +Module : SLynx.Simulate.PhyloModel +Description : Parse and interpret the model string +Copyright : (c) Dominik Schrempf 2020 +License : GPL-3.0-or-later + +Maintainer : dominik.schrempf@gmail.com +Stability : unstable +Portability : portable + +Creation date: Fri Feb 1 13:32:16 2019. + +-} + +module SLynx.Simulate.PhyloModel + ( getPhyloModel + ) +where + +import Control.Monad ( when ) +import qualified Data.ByteString.Lazy.Char8 as L +import Data.Either ( rights ) +import Data.List.NonEmpty ( fromList ) +import Data.Maybe +import Data.Scientific hiding ( scientific ) +import Data.Void +import Data.Word ( Word8 ) +import Numeric.LinearAlgebra ( norm_1 + , size + , vector + ) +import Text.Megaparsec +import Text.Megaparsec.Byte +import Text.Megaparsec.Byte.Lexer + +import ELynx.Data.MarkovProcess.AminoAcid +import ELynx.Data.MarkovProcess.CXXModels +import qualified ELynx.Data.MarkovProcess.MixtureModel + as M +import ELynx.Data.MarkovProcess.Nucleotide +import qualified ELynx.Data.MarkovProcess.PhyloModel + as P +import ELynx.Data.MarkovProcess.RateMatrix +import qualified ELynx.Data.MarkovProcess.SubstitutionModel + as S +import ELynx.Import.MarkovProcess.EDMModelPhylobayes + ( EDMComponent ) +import ELynx.Tools + +type Parser = Parsec Void L.ByteString + +bs :: String -> L.ByteString +bs = L.pack + +nNuc :: Int +-- nNuc = length (alphabet :: [Nucleotide]) +nNuc = 4 + +nAA :: Int +-- nAA = length (alphabet :: [AminoAcid]) +nAA = 20 + +-- Model parameters between square brackets. +paramsStart :: Word8 +paramsStart = c2w '[' + +paramsEnd :: Word8 +paramsEnd = c2w ']' + +-- Stationary distribution between curly brackets. +sdStart :: Word8 +sdStart = c2w '{' + +sdEnd :: Word8 +sdEnd = c2w '}' + +-- Mixture model components between round brackets. +mmStart :: Word8 +mmStart = c2w '(' + +mmEnd :: Word8 +mmEnd = c2w ')' + +separator :: Word8 +separator = c2w ',' + +name :: Parser String +name = L.unpack <$> takeWhile1P + (Just "Model name") + (`notElem` [paramsStart, paramsEnd, sdStart, sdEnd, mmStart, mmEnd, separator] + ) + +params :: Parser [Double] +params = map toRealFloat <$> between (char paramsStart) + (char paramsEnd) + (sepBy1 scientific (char separator)) + +stationaryDistribution :: Parser StationaryDistribution +stationaryDistribution = do + f <- vector . map toRealFloat <$> between + (char sdStart) + (char sdEnd) + (sepBy1 scientific (char separator)) + if nearlyEq (norm_1 f) 1.0 + then return f + else + error + $ "Sum of stationary distribution is " + ++ show (norm_1 f) + ++ " but should be 1.0." + +assertLength :: StationaryDistribution -> Int -> a -> a +assertLength d n r = if size d /= n + then + error + $ "Length of stationary distribution is " + ++ show (size d) + ++ " but should be " + ++ show n + ++ "." + else r + +-- This is the main function that connects the model string, the parameters and +-- the stationary distribution. It should check that the model is valid. +assembleSubstitutionModel + :: String + -> Maybe S.Params + -> Maybe StationaryDistribution + -> Either String S.SubstitutionModel +-- DNA models. +assembleSubstitutionModel "JC" Nothing Nothing = Right jc +assembleSubstitutionModel "F81" Nothing (Just d) = + Right $ assertLength d nNuc $ f81 d +assembleSubstitutionModel "HKY" (Just [k]) (Just d) = + Right $ assertLength d nNuc $ hky k d +assembleSubstitutionModel "GTR4" (Just es) (Just d) = + Right $ assertLength d nNuc $ gtr4 es d +-- Protein models. +assembleSubstitutionModel "Poisson" Nothing Nothing = Right poisson +assembleSubstitutionModel "Poisson-Custom" Nothing (Just d) = + Right $ assertLength d nAA $ poissonCustom Nothing d +assembleSubstitutionModel "LG" Nothing Nothing = Right lg +assembleSubstitutionModel "LG-Custom" Nothing (Just d) = + Right $ assertLength d nAA $ lgCustom Nothing d +assembleSubstitutionModel "WAG" Nothing Nothing = Right wag +assembleSubstitutionModel "WAG-Custom" Nothing (Just d) = + Right $ assertLength d nAA $ wagCustom Nothing d +assembleSubstitutionModel "GTR20" (Just es) (Just d) = + Right $ assertLength d nAA $ gtr20 es d +-- Ohterwisse, we cannot assemble the model. +assembleSubstitutionModel n mps mf = Left $ unlines + [ "Cannot assemble substitution model." + , "Name: " ++ show n + , "Parameters: " ++ show mps + , "Stationary distribution: " ++ show mf + ] + +parseSubstitutionModel :: Parser S.SubstitutionModel +parseSubstitutionModel = do + n <- name + mps <- optional params + mf <- optional stationaryDistribution + let esm = assembleSubstitutionModel n mps mf + case esm of + Left err -> fail err + Right sm -> return sm + +edmModel :: [EDMComponent] -> Maybe [M.Weight] -> Parser M.MixtureModel +edmModel cs mws = do + _ <- chunk (bs "EDM") + _ <- char mmStart + n <- name + mps <- optional params + _ <- char mmEnd + when (null cs) $ error "edmModel: no EDM components given." + let sms = map (\c -> assembleSubstitutionModel n mps (Just $ snd c)) cs + edmName = "EDM" ++ show (length cs) + ws = fromMaybe (map fst cs) mws + errs = [ e | (Left e) <- sms ] + when (length sms /= length ws) + $ error "edmModel: number of substitution models and weights differs." + if not $ null errs + then fail $ head errs + else return + $ M.fromSubstitutionModels edmName (fromList ws) (fromList $ rights sms) + +cxxModel :: Maybe [M.Weight] -> Parser M.MixtureModel +cxxModel mws = do + _ <- char (c2w 'C') + n <- decimal :: Parser Int + return $ cxx n mws + +standardMixtureModel :: [M.Weight] -> Parser M.MixtureModel +standardMixtureModel ws = do + _ <- chunk (bs "MIXTURE") + _ <- char mmStart + sms <- parseSubstitutionModel `sepBy1` char separator + _ <- char mmEnd + -- XXX: The use of `Data.List.NonEmpty.fromList` leads to uninformative error messages. + return $ M.fromSubstitutionModels "MIXTURE" (fromList ws) (fromList sms) + +mixtureModel + :: Maybe [EDMComponent] -> Maybe [M.Weight] -> Parser M.MixtureModel +mixtureModel Nothing Nothing = + try (cxxModel Nothing) <|> fail "No weights provided." +mixtureModel Nothing mws@(Just ws) = + try (cxxModel mws) <|> standardMixtureModel ws +mixtureModel (Just cs) mws = edmModel cs mws + +-- | Parse the phylogenetic model string. The argument list is somewhat long, +-- but models can have many parameters and we have to check for redundant +-- parameters. +-- +-- @ +-- getPhyloModel maybeSubstitutionModelString maybeMixtureModelString maybeEDMComponents +-- @ +getPhyloModel + :: Maybe String + -> Maybe String + -> Maybe [M.Weight] + -> Maybe [EDMComponent] + -> Either String P.PhyloModel +getPhyloModel Nothing Nothing _ _ = Left "No model was given. See help." +getPhyloModel (Just _) (Just _) _ _ = + Left "Both, substitution and mixture model string given; use only one." +getPhyloModel (Just s) Nothing Nothing Nothing = + Right $ P.SubstitutionModel $ parseStringWith "Substitution model string" + parseSubstitutionModel + s +getPhyloModel (Just _) Nothing (Just _) _ = + Left "Weights given; but cannot be used with substitution model." +getPhyloModel (Just _) Nothing _ (Just _) = + Left + "Empirical distribution mixture model components given; but cannot be used with substitution model." +getPhyloModel Nothing (Just m) mws mcs = + Right $ P.MixtureModel $ parseStringWith "Mixture model string" + (mixtureModel mcs mws) + m diff --git a/src/SLynx/Simulate/Simulate.hs b/src/SLynx/Simulate/Simulate.hs new file mode 100644 --- /dev/null +++ b/src/SLynx/Simulate/Simulate.hs @@ -0,0 +1,228 @@ +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE TemplateHaskell #-} + +{- | +Module : SLynx.Simulate.Simulate +Description : Simulate multiple sequence alignments +Copyright : (c) Dominik Schrempf 2020 +License : GPL-3.0-or-later + +Maintainer : dominik.schrempf@gmail.com +Stability : unstable +Portability : portable + +Creation date: Mon Jan 28 14:12:52 2019. + +-} + +module SLynx.Simulate.Simulate + ( simulateCmd + ) +where + +import Control.Applicative ( (<|>) ) +import Control.Concurrent +import Control.Concurrent.Async +import Control.Monad ( unless + , when + ) +import Control.Monad.IO.Class +import Control.Monad.Logger +import Control.Monad.Trans.Class +import Control.Monad.Trans.Reader ( ask ) +import qualified Data.ByteString.Lazy as L +import qualified Data.ByteString.Lazy.Char8 as LC +import Data.List.NonEmpty ( toList ) +import Data.Maybe +import qualified Data.Set as Set +import qualified Data.Text as T +import qualified Data.Text.Lazy as LT +import qualified Data.Text.Lazy.Encoding as LT +import Data.Tree +import qualified Data.Vector.Unboxed as V +import Numeric.LinearAlgebra hiding ( (<>) + , toList + ) +import System.Random.MWC + +import SLynx.Simulate.Options +import SLynx.Simulate.PhyloModel + +import ELynx.Data.Alphabet.Alphabet as A +import ELynx.Data.MarkovProcess.GammaRateHeterogeneity +import qualified ELynx.Data.MarkovProcess.MixtureModel + as M +import qualified ELynx.Data.MarkovProcess.PhyloModel + as P +import qualified ELynx.Data.MarkovProcess.SubstitutionModel + as SM +import qualified ELynx.Data.Sequence.Alignment as A +import qualified ELynx.Data.Sequence.Sequence as Seq + hiding ( name ) +import ELynx.Data.Tree +import ELynx.Export.Sequence.Fasta +import ELynx.Import.MarkovProcess.EDMModelPhylobayes + hiding ( Parser ) +import ELynx.Import.MarkovProcess.SiteprofilesPhylobayes +import ELynx.Import.Tree.Newick + hiding ( name ) +import ELynx.Simulate.MarkovProcessAlongTree + +import ELynx.Tools + +-- Simulate a 'Alignment' for a given phylogenetic model, +-- phylogenetic tree, and alignment length. +simulateAlignment + :: (Measurable a, Named a) + => P.PhyloModel + -> Tree a + -> Int + -> GenIO + -> IO A.Alignment +simulateAlignment pm t n g = do + c <- getNumCapabilities + gs <- splitGen c g + let chunks = getChunks c n + leafStatesS <- case pm of + -- XXX @performace: This parallelization is not very intelligent, because + -- the matrices exponentiation is done in all threads. So ten threads will + -- exponentiate the same matrix ten times. + P.SubstitutionModel sm -> mapConcurrently + (\(num, gen) -> simulateAndFlatten num d e t gen) + (zip chunks gs) + where + d = SM.stationaryDistribution sm + e = SM.exchangeabilityMatrix sm + -- P.MixtureModel mm -> mapConcurrently + -- (\(num, gen) -> simulateAndFlattenNSitesAlongTreeMixtureModel num ws ds es t gen) (zip chunks gs) + P.MixtureModel mm -> simulateAndFlattenMixtureModelPar n ws ds es t g + where + ws = vector . toList $ M.getWeights mm + ds = map SM.stationaryDistribution $ toList $ M.getSubstitutionModels mm + es = map SM.exchangeabilityMatrix $ toList $ M.getSubstitutionModels mm + -- XXX @performace. The horizontal concatenation might be slow. If so, + -- 'concatenateSeqs' or 'concatenateAlignments' can be used, which directly + -- appends vectors. + let leafStates = horizontalConcat leafStatesS + leafNames = map getName $ leaves t + code = P.getAlphabet pm + -- XXX: Probably use type safe stuff here? + alph = A.all $ alphabetSpec code + sequences = + [ Seq.Sequence sName "" code (V.fromList $ map (`Set.elemAt` alph) ss) + | (sName, ss) <- zip leafNames leafStates + ] + return $ either error id $ A.fromSequences sequences + +-- Summarize EDM components; line to be printed to screen or log. +summarizeEDMComponents :: [EDMComponent] -> L.ByteString +summarizeEDMComponents cs = + LC.pack + $ "Empiricial distribution mixture model with " + ++ show (length cs) + ++ " components." + +-- XXX. Maybe provide human readable model file. But then, why is this +-- necessary. A human readable summary is reported anyways, and for Protein +-- models the exchangeabilities are too many. +reportModel :: P.PhyloModel -> ELynx SimulateArguments () +reportModel m = do + bn <- outFileBaseName . global <$> ask + case bn of + Nothing -> + $(logInfo) + "No output file provided; omit writing machine-readable phylogenetic model." + Just _ -> + out "model definition (machine readable)" (bShow m <> "\n") ".model.gz" + +-- | Simulate sequences. +simulateCmd :: ELynx SimulateArguments () +simulateCmd = do + l <- local <$> ask + let treeFile = argsTreeFile l + + $(logInfo) "" + $(logInfo) $ T.pack $ "Read tree from file '" ++ treeFile ++ "'." + tree <- liftIO $ parseFileWith (newick Standard) treeFile + $(logInfo) $ LT.toStrict $ LT.decodeUtf8 $ summarize tree + + let edmFile = argsEDMFile l + let sProfileFiles = argsSiteprofilesFiles l + $(logInfo) "" + $(logDebug) "Read EDM file or siteprofile files." + when (isJust edmFile && isJust sProfileFiles) + $ error "Got both: --edm-file and --siteprofile-files." + edmCs <- case edmFile of + Nothing -> return Nothing + Just edmF -> do + $(logInfo) "Read EDM file." + liftIO $ Just <$> parseFileWith phylobayes edmF + maybe (return ()) + ($(logInfo) . LT.toStrict . LT.decodeUtf8 . summarizeEDMComponents) + edmCs + sProfiles <- case sProfileFiles of + Nothing -> return Nothing + Just fns -> do + $(logInfo) + $ T.pack + $ "Read siteprofiles from " + ++ show (length fns) + ++ " file(s)." + $(logDebug) $ T.pack $ "The file names are:" ++ show fns + xs <- liftIO $ mapM (parseFileWith siteprofiles) fns + return $ Just $ concat xs + maybe (return ()) + ($(logInfo) . LT.toStrict . LT.decodeUtf8 . summarizeEDMComponents) + sProfiles + let edmCsOrSiteprofiles = edmCs <|> sProfiles + + $(logInfo) "Read model string." + let ms = argsSubstitutionModelString l + mm = argsMixtureModelString l + mws = argsMixtureWeights l + eitherPhyloModel' = getPhyloModel ms mm mws edmCsOrSiteprofiles + phyloModel' <- case eitherPhyloModel' of + Left err -> lift $ error err + Right pm -> return pm + + let maybeGammaParams = argsGammaParams l + phyloModel <- case maybeGammaParams of + Nothing -> do + $(logInfo) $ LT.toStrict $ LT.decodeUtf8 $ LC.unlines $ P.summarize + phyloModel' + return phyloModel' + Just (n, alpha) -> do + $(logInfo) + $ LT.toStrict + $ LT.decodeUtf8 + $ LC.intercalate "\n" + $ P.summarize phyloModel' + $(logInfo) "" + $(logInfo) + $ LT.toStrict + $ LT.decodeUtf8 + $ LC.intercalate "\n" + $ summarizeGammaRateHeterogeneity n alpha + return $ expand n alpha phyloModel' + + -- -- XXX: Do not report possibly huge empirical distribution mixture models + -- -- for now, because it takes too long and uses too much disk space :). + unless + (isJust sProfiles) + (do + $(logInfo) "" + reportModel phyloModel + $(logInfo) "" + ) + + $(logInfo) "Simulate alignment." + let alignmentLength = argsLength l + $(logInfo) $ T.pack $ "Length: " <> show alignmentLength <> "." + gen <- case argsSeed l of + Random -> + error "simulateCmd: seed not available; please contact maintainer." + Fixed s -> liftIO $ initialize s + alignment <- liftIO $ simulateAlignment phyloModel tree alignmentLength gen + let output = (sequencesToFasta . A.toSequences) alignment + $(logInfo) "" + out "simulated multi sequence alignment" output ".fasta" diff --git a/src/SLynx/SubSample/Options.hs b/src/SLynx/SubSample/Options.hs new file mode 100644 --- /dev/null +++ b/src/SLynx/SubSample/Options.hs @@ -0,0 +1,99 @@ +{-# LANGUAGE DeriveGeneric #-} + +{- | +Module : SLynx.SubSample.Options +Description : ELynxSeq argument parsing +Copyright : (c) Dominik Schrempf 2018 +License : GPL-3.0-or-later + +Maintainer : dominik.schrempf@gmail.com +Stability : unstable +Portability : portable + +Creation date: Sun Oct 7 17:29:45 2018. + +-} + +module SLynx.SubSample.Options + ( SubSampleArguments(..) + , subSampleArguments + , getOutSuffixes + ) +where + +import Options.Applicative +import qualified Data.Text as T +import qualified Data.Text.Lazy as LT +import qualified Data.Text.Lazy.Builder as LT +import qualified Data.Text.Lazy.Builder.Int as LT + +import SLynx.Tools + +import ELynx.Data.Alphabet.Alphabet +import ELynx.Tools + +-- | Data structure holding the Command line arguments. +data SubSampleArguments = SubSampleArguments + { ssAlphabet :: Alphabet + , ssInFile :: FilePath + , ssNSites :: Int + , ssNAlignments :: Int + , ssMbSeed :: Seed } + deriving (Eq, Show, Generic) + +-- | Get a given number of output file suffixes. +-- +-- > getOutSuffixes 11 "fasta" +-- +-- Will result in @.00.fasta@ up to @.10.fasta@. +getOutSuffixes :: Int -> String -> [String] +getOutSuffixes n suffix = + [ "." ++ digitStr i ++ "." ++ suffix | i <- [0 .. n - 1] ] + where + nDigits = ceiling $ logBase (10 :: Double) (fromIntegral n) + digitStr i = T.unpack + $ T.justifyRight nDigits '0' (LT.toStrict $ LT.toLazyText $ LT.decimal i) + + +instance Reproducible SubSampleArguments where + inFiles = pure . ssInFile + outSuffixes a = getOutSuffixes (ssNAlignments a) "fasta" + getSeed = Just . ssMbSeed + setSeed a s = a { ssMbSeed = Fixed s } + parser = subSampleArguments + cmdName = "sub-sample" + cmdDsc = ["Sub-sample columns from multi sequence alignments."] + cmdFtr = + [ "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." + ] + +instance FromJSON SubSampleArguments + +instance ToJSON SubSampleArguments + +-- | Sub command parser. +subSampleArguments :: Parser SubSampleArguments +subSampleArguments = + SubSampleArguments + <$> alphabetOpt + <*> filePathArg + <*> subSampleNSitesOpt + <*> subSampleNAlignmentsOpt + <*> seedOpt + +subSampleNSitesOpt :: Parser Int +subSampleNSitesOpt = + option auto $ long "number-of-sites" <> short 'n' <> metavar "INT" <> help + "Number of sites randomly drawn with replacement" + +subSampleNAlignmentsOpt :: Parser Int +subSampleNAlignmentsOpt = + option auto + $ long "number-of-alignments" + <> short 'm' + <> metavar "INT" + <> help "Number of multi sequence alignments to be created" + +filePathArg :: Parser FilePath +filePathArg = + strArgument $ metavar "INPUT-FILE" <> help "Read sequences from INPUT-FILE" diff --git a/src/SLynx/SubSample/SubSample.hs b/src/SLynx/SubSample/SubSample.hs new file mode 100644 --- /dev/null +++ b/src/SLynx/SubSample/SubSample.hs @@ -0,0 +1,54 @@ +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE TemplateHaskell #-} + +{- | +Module : Analyze.Analyze +Description : Parse sequence file formats and analyze them +Copyright : (c) Dominik Schrempf 2018 +License : GPL-3.0-or-later + +Maintainer : dominik.schrempf@gmail.com +Stability : unstable +Portability : portable + +Creation date: Fri Oct 5 08:41:05 2018. + +-} + +module SLynx.SubSample.SubSample + ( subSampleCmd + ) +where + +import Control.Monad +import Control.Monad.IO.Class +import Control.Monad.Logger +import Control.Monad.Trans.Reader ( ask ) +import qualified Data.Text as T +import System.Random.MWC + +import SLynx.SubSample.Options +import SLynx.Tools + +import qualified ELynx.Data.Sequence.Alignment as M +import ELynx.Export.Sequence.Fasta +import ELynx.Tools + +-- | Sub sample sequences. +subSampleCmd :: ELynx SubSampleArguments () +subSampleCmd = do + (SubSampleArguments al inFile nSites nAlignments (Fixed s)) <- local <$> ask + $(logInfo) "Command: Sub sample from a multi sequence alignment." + $(logInfo) $ T.pack $ " Sample " <> show nSites <> " sites." + $(logInfo) + $ T.pack + $ " Sample " + <> show nAlignments + <> " multi sequence alignments." + ss <- readSeqs al inFile + gen <- liftIO $ initialize s + let a = either error id (M.fromSequences ss) + samples <- liftIO $ replicateM nAlignments $ M.randomSubSample nSites a gen + let results = map (sequencesToFasta . M.toSequences) samples + sfxs = getOutSuffixes nAlignments "fasta" + zipWithM_ (out "sub sampled multi sequence alignments") results sfxs diff --git a/src/SLynx/Tools.hs b/src/SLynx/Tools.hs new file mode 100644 --- /dev/null +++ b/src/SLynx/Tools.hs @@ -0,0 +1,74 @@ +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE TemplateHaskell #-} + +{- | +Module : SLynx.Tools +Description : Common tools for sequence lynx +Copyright : (c) Dominik Schrempf 2020 +License : GPL-3.0-or-later + +Maintainer : dominik.schrempf@gmail.com +Stability : unstable +Portability : portable + +Creation date: Sat Sep 7 06:24:22 2019. + +-} + +module SLynx.Tools + ( -- * SLynx.Tools + readSeqs + -- * Options + , alphabetOpt + ) +where + +import Control.Monad.IO.Class +import Control.Monad.Logger +import qualified Data.Text as T +import Options.Applicative + +import ELynx.Data.Alphabet.Alphabet +import ELynx.Data.Sequence.Sequence +import ELynx.Import.Sequence.Fasta + hiding ( Parser ) +import ELynx.Tools + +-- -- | Read sequences of given alphabet from file or standard input. +-- readSeqs +-- :: (MonadIO m, MonadLogger m) => Alphabet -> Maybe FilePath -> m [Sequence] +-- readSeqs a mfp = do +-- case mfp of +-- Nothing -> +-- $(logInfo) +-- $ T.pack +-- $ "Read sequences from standard input; alphabet " +-- <> show a +-- <> "." +-- Just fp -> +-- $(logInfo) +-- $ T.pack +-- $ "Read sequences from file " +-- <> fp +-- <> "; alphabet " +-- <> show a +-- <> "." +-- liftIO $ parseFileOrIOWith (fasta a) mfp + +-- | Read sequences of given alphabet from file or standard input. +readSeqs :: (MonadIO m, MonadLogger m) => Alphabet -> FilePath -> m [Sequence] +readSeqs a fp = do + $(logInfo) + $ T.pack + $ "Read sequences from file " + <> fp + <> "; alphabet " + <> show a + <> "." + liftIO $ parseFileWith (fasta a) fp + +-- | Command line option to specify the alphabet. Used by various commands. +alphabetOpt :: Parser Alphabet +alphabetOpt = + option auto $ long "alphabet" <> short 'a' <> metavar "NAME" <> help + "Specify alphabet type NAME" diff --git a/src/SLynx/Translate/Options.hs b/src/SLynx/Translate/Options.hs new file mode 100644 --- /dev/null +++ b/src/SLynx/Translate/Options.hs @@ -0,0 +1,78 @@ +{-# LANGUAGE DeriveGeneric #-} + +{- | +Module : SLynx.Translate.Options +Description : ELynxSeq argument parsing +Copyright : (c) Dominik Schrempf 2018 +License : GPL-3.0-or-later + +Maintainer : dominik.schrempf@gmail.com +Stability : unstable +Portability : portable + +Creation date: Sun Oct 7 17:29:45 2018. + +-} + +module SLynx.Translate.Options + ( TranslateArguments(..) + , translateArguments + ) +where + +import Data.List +import Options.Applicative + +import SLynx.Tools + +import ELynx.Data.Alphabet.Alphabet +import ELynx.Data.Character.Codon +import ELynx.Tools + +-- | Arguments needed to translate sequences. +data TranslateArguments = TranslateArguments + { trAlphabet :: Alphabet + , trInFile :: FilePath + , trReadingFrame :: Int + , trUniversalCode :: UniversalCode } + deriving (Eq, Show, Generic) + +instance Reproducible TranslateArguments where + inFiles = pure . trInFile + outSuffixes _ = [".fasta"] + getSeed _ = Nothing + setSeed = const + parser = translateArguments + cmdName = "translate" + cmdDsc = ["Translate from DNA to Protein or DNAX to ProteinX."] + +instance FromJSON TranslateArguments + +instance ToJSON TranslateArguments + +-- | Command line parser. +translateArguments :: Parser TranslateArguments +translateArguments = + TranslateArguments + <$> alphabetOpt + <*> inFileArg + <*> readingFrameOpt + <*> universalCodeOpt + +readingFrameOpt :: Parser Int +readingFrameOpt = + option auto $ long "reading-frame" <> short 'r' <> metavar "INT" <> help + "Reading frame [0|1|2]." + +universalCodeOpt :: Parser UniversalCode +universalCodeOpt = + option auto $ long "universal-code" <> short 'u' <> metavar "CODE" <> help + ("universal code; one of: " ++ codeStr ++ ".") + where + codes = allValues :: [UniversalCode] + codeWords = map show codes + codeStr = intercalate ", " codeWords + +inFileArg :: Parser FilePath +inFileArg = + strArgument $ metavar "INPUT-FILE" <> help "Read sequences from INPUT-FILE" diff --git a/src/SLynx/Translate/Translate.hs b/src/SLynx/Translate/Translate.hs new file mode 100644 --- /dev/null +++ b/src/SLynx/Translate/Translate.hs @@ -0,0 +1,49 @@ +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE TemplateHaskell #-} + +{- | +Module : Analyze.Analyze +Description : Parse sequence file formats and analyze them +Copyright : (c) Dominik Schrempf 2018 +License : GPL-3.0-or-later + +Maintainer : dominik.schrempf@gmail.com +Stability : unstable +Portability : portable + +Creation date: Fri Oct 5 08:41:05 2018. + +-} + +module SLynx.Translate.Translate + ( translateCmd + ) +where + +import Control.Monad.Logger +import Control.Monad.Trans.Reader ( ask ) +import qualified Data.Text as T + +import SLynx.Tools +import SLynx.Translate.Options + +import ELynx.Data.Character.Codon +import ELynx.Data.Sequence.Sequence +import ELynx.Data.Sequence.Translate +import ELynx.Export.Sequence.Fasta +import ELynx.Tools + +translateSeqs :: Int -> UniversalCode -> [Sequence] -> [Sequence] +translateSeqs rf uc = map (translateSeq uc rf) + +-- | Translate sequences. +translateCmd :: ELynx TranslateArguments () +translateCmd = do + (TranslateArguments al inFile rf uc) <- local <$> ask + $(logInfo) "Command: Translate sequences to amino acids." + $(logInfo) $ T.pack $ " Universal code: " <> show uc <> "." + $(logInfo) $ T.pack $ " Reading frame: " <> show rf <> "." + $(logInfo) "" + ss <- readSeqs al inFile + let result = sequencesToFasta $ translateSeqs rf uc ss + out "translated sequences" result "fasta" diff --git a/src/Simulate/Options.hs b/src/Simulate/Options.hs deleted file mode 100644 --- a/src/Simulate/Options.hs +++ /dev/null @@ -1,188 +0,0 @@ -{- | -Module : Simulate.Options -Description : ELynxSim argument parsing -Copyright : (c) Dominik Schrempf 2018 -License : GPL-3 - -Maintainer : dominik.schrempf@gmail.com -Stability : unstable -Portability : portable - -Creation date: Sun Oct 7 17:29:45 2018. - -Available options: - -h,--help Show this help text - -v,--version Show version - -t,--tree-file NAME Specify tree file NAME - -s,--substitution-model MODEL - Set the phylogenetic substitution model; available - models are shown below - -m,--mixture-model MODEL Set the phylogenetic mixture model; available models - are shown below - -l,--length NUMBER Set alignment length to NUMBER - -e,--edm-file NAME empirical distribution model file NAME in Phylobayes - format - -w,--mixture-model-weights [DOUBLE,DOUBLE,...] - weights of mixture model components - -g,--gamma-rate-heterogeneity (NCAT, SHAPE) - number of gamma rate categories and shape parameter - -e,--seed [INT] Set seed for the random number generator; list of 32 - bit integers with up to 256 elements (default: [0]) - -q,--quiet Be quiet - -o,--output-file NAME Specify output file NAME - - --} - - -module Simulate.Options - ( GammaRateHeterogeneityParams - , SimulateArguments (..) - , Simulate - , simulateArguments - , simulateFooter - ) where - -import Control.Monad.Logger -import Control.Monad.Trans.Reader -import Data.List -import Data.String (words) --- import Data.Void -import Data.Word -import Options.Applicative --- import Text.Megaparsec (Parsec, takeWhile1P) --- import Text.Megaparsec.Char (space) - -import ELynx.Tools.Options - --- | Number of gamma rate categories and alpha parameter. -type GammaRateHeterogeneityParams = (Int, Double) - --- | Arguments needed to simulate sequences. -data SimulateArguments = SimulateArguments - { argsTreeFile :: FilePath - , argsSubstitutionModelString :: Maybe String - , argsMixtureModelString :: Maybe String - , argsEDMFile :: Maybe FilePath - , argsSiteprofilesFiles :: Maybe [FilePath] - , argsMixtureWeights :: Maybe [Double] - , argsGammaParams :: Maybe GammaRateHeterogeneityParams - , argsLength :: Int - , argsMaybeSeed :: Maybe [Word32] - } - --- | A simulation can log stuff and read the necessary arguments. -type Simulate = LoggingT (ReaderT SimulateArguments IO) - --- | Sub command parser. -simulateArguments :: Parser SimulateArguments -simulateArguments = SimulateArguments - <$> treeFileOpt - <*> phyloSubstitutionModelOpt - <*> phyloMixtureModelOpt - <*> maybeEDMFileOpt - <*> maybeSiteprofilesFilesOpt - <*> maybeMixtureWeights - <*> maybeGammaParams - <*> lengthOpt - <*> seedOpt - -treeFileOpt :: Parser FilePath -treeFileOpt = strOption $ - long "tree-file" - <> short 't' - <> metavar "Name" - <> help "Read trees from file NAME" - -phyloSubstitutionModelOpt :: Parser (Maybe String) -phyloSubstitutionModelOpt = optional $ strOption $ - long "substitution-model" - <> short 's' - <> metavar "MODEL" - <> help "Set the phylogenetic substitution model; available models are shown below (mutually exclusive with -m option)" - -phyloMixtureModelOpt :: Parser (Maybe String) -phyloMixtureModelOpt = optional $ strOption - ( long "mixture-model" - <> short 'm' - <> metavar "MODEL" - <> help "Set the phylogenetic mixture model; available models are shown below (mutually exclusive with -s option)" ) - -maybeEDMFileOpt :: Parser (Maybe FilePath) -maybeEDMFileOpt = optional $ strOption - ( long "edm-file" - <> short 'e' - <> metavar "NAME" - <> help "Empirical distribution model file NAME in Phylobayes format" ) - --- fn :: Parsec Void String FilePath --- fn = space *> takeWhile1P () <* space - --- fns :: Parsec Void String [FilePath] --- fns = do --- many string - -maybeSiteprofilesFilesOpt :: Parser (Maybe [FilePath]) -maybeSiteprofilesFilesOpt = optional $ words <$> strOption - ( long "siteprofile-files" - <> short 'p' - <> metavar "NAMES" - <> help "File names of site profiles in Phylobayes format" ) - -maybeMixtureWeights :: Parser (Maybe [Double]) -maybeMixtureWeights = optional $ option auto - ( long "mixture-model-weights" - <> short 'w' - <> metavar "\"[DOUBLE,DOUBLE,...]\"" - <> help "Weights of mixture model components" ) - -maybeGammaParams :: Parser (Maybe GammaRateHeterogeneityParams) -maybeGammaParams = optional $ option auto - ( long "gamma-rate-heterogeneity" - <> short 'g' - <> metavar "\"(NCAT,SHAPE)\"" - <> help "Number of gamma rate categories and shape parameter" ) - -lengthOpt :: Parser Int -lengthOpt = option auto - ( long "length" - <> short 'l' - <> metavar "NUMBER" - <> help "Set alignment length to NUMBER" ) - --- | The model specification is somewhat complicated, so we need to provide --- additional help. -simulateFooter :: String -simulateFooter = intercalate "\n" $ sms ++ mms - where - sms = - [ "Substitution models:" - , "-s \"MODEL[PARAMETER,PARAMETER,...]{STATIONARY_DISTRIBUTION}\"" - , " Supported DNA models: JC, HKY." - , " For example," - , " -s HKY[KAPPA]{DOUBLE,DOUBLE,DOUBLE,DOUBLE}" - , " Supported Protein models: Poisson, Poisson-Custom, LG, LG-Custom, WAG, WAG-Custom." - , " MODEL-Custom means that only the exchangeabilities of MODEL are used," - , " and a custom stationary distribution is provided." - , " For example," - , " -s LG-Custom{...}" - ] - mms = - [ "" - , "Mixture models:" - , "-m \"MIXTURE(SUBSTITUTION_MODEL_1,SUBSTITUTION_MODEL_2)\"" - , " For example," - , " -m \"MIXTURE(JC,HKY[6.0]{0.3,0.2,0.2,0.3})\"" - , "Mixture weights have to be provided with the -w option." - , "" - , "Special mixture models:" - , "-m CXX" - , " where XX is 10, 20, 30, 40, 50, or 60; CXX models, Quang et al., 2008." - , "-m \"EDM(EXCHANGEABILITIES)\"" - , " Arbitrary empirical distribution mixture (EDM) models." - , " Stationary distributions have to be provided with the -e option." - , " For example," - , " LG exchangeabilities with stationary distributions given in FILE." - , " -m \"EDM(LG-Custom)\" -e FILE" - , "For special mixture models, mixture weights are optional." - ] diff --git a/src/Simulate/PhyloModel.hs b/src/Simulate/PhyloModel.hs deleted file mode 100644 --- a/src/Simulate/PhyloModel.hs +++ /dev/null @@ -1,188 +0,0 @@ -{- | -Module : Simulate.PhyloModel -Description : Parse and interpret the model string -Copyright : (c) Dominik Schrempf 2019 -License : GPL-3 - -Maintainer : dominik.schrempf@gmail.com -Stability : unstable -Portability : portable - -Creation date: Fri Feb 1 13:32:16 2019. - --} - -module Simulate.PhyloModel - ( getPhyloModel - ) where - -import qualified Data.ByteString.Lazy.Char8 as L -import Data.Maybe -import Data.Void -import Data.Word (Word8) -import Numeric.LinearAlgebra (norm_1, size, - vector) -import Text.Megaparsec -import Text.Megaparsec.Byte -import Text.Megaparsec.Byte.Lexer - -import ELynx.Data.MarkovProcess.AminoAcid -import ELynx.Data.MarkovProcess.CXXModels -import qualified ELynx.Data.MarkovProcess.MixtureModel as M -import ELynx.Data.MarkovProcess.Nucleotide -import qualified ELynx.Data.MarkovProcess.PhyloModel as P -import ELynx.Data.MarkovProcess.RateMatrix -import qualified ELynx.Data.MarkovProcess.SubstitutionModel as S -import ELynx.Import.MarkovProcess.EDMModelPhylobayes (EDMComponent) -import ELynx.Tools.ByteString -import ELynx.Tools.Equality -import ELynx.Tools.InputOutput - -type Parser = Parsec Void L.ByteString - -bs :: String -> L.ByteString -bs = L.pack - -nNuc :: Int --- nNuc = length (alphabet :: [Nucleotide]) -nNuc = 4 - -nAA :: Int --- nAA = length (alphabet :: [AminoAcid]) -nAA = 20 - --- Model parameters between square brackets. -paramsStart :: Word8 -paramsStart = c2w '[' - -paramsEnd :: Word8 -paramsEnd = c2w ']' - --- Stationary distribution between curly brackets. -sdStart :: Word8 -sdStart = c2w '{' - -sdEnd :: Word8 -sdEnd = c2w '}' - --- Mixture model components between round brackets. -mmStart :: Word8 -mmStart = c2w '(' - -mmEnd :: Word8 -mmEnd = c2w ')' - -separator :: Word8 -separator = c2w ',' - -name :: Parser String -name = L.unpack <$> - takeWhile1P (Just "Model name") (`notElem` [paramsStart, paramsEnd, sdStart, sdEnd, mmStart, mmEnd, separator]) - -params :: Parser [Double] -params = between (char paramsStart) (char paramsEnd) (sepBy1 float (char separator)) - -stationaryDistribution :: Parser StationaryDistribution -stationaryDistribution = do - f <- vector <$> between (char sdStart) (char sdEnd) (sepBy1 float (char separator)) - if nearlyEq (norm_1 f) 1.0 - then return f - else error $ "Sum of stationary distribution is " ++ show (norm_1 f) - ++ " but should be 1.0." - -assertLength :: StationaryDistribution -> Int -> a -> a -assertLength d n r = if size d /= n - then error $ "Length of stationary distribution is " ++ show (size d) - ++ " but should be " ++ show n ++ "." - else r - --- This is the main function that connects the model string, the parameters and --- the stationary distribution. It should check that the model is valid. -assembleSubstitutionModel :: String -> Maybe S.Params -> Maybe StationaryDistribution - -> Either String S.SubstitutionModel --- DNA models. -assembleSubstitutionModel "JC" Nothing Nothing = Right jc -assembleSubstitutionModel "HKY" (Just [k]) (Just d) = Right $ assertLength d nNuc $ hky k d --- Protein models. -assembleSubstitutionModel "LG" Nothing Nothing = Right lg -assembleSubstitutionModel "LG-Custom" Nothing (Just d) = Right $ assertLength d nAA $ lgCustom Nothing d -assembleSubstitutionModel "WAG" Nothing Nothing = Right wag -assembleSubstitutionModel "WAG-Custom" Nothing (Just d) = Right $ assertLength d nAA $ wagCustom Nothing d -assembleSubstitutionModel "Poisson" Nothing Nothing = Right poisson -assembleSubstitutionModel "Poisson-Custom" Nothing (Just d) = Right $ assertLength d nAA $ poissonCustom Nothing d --- Ohterwisse, we cannot assemble the model. -assembleSubstitutionModel n mps mf = Left $ unlines - [ "Cannot assemble substitution model. " - , "Name: " ++ show n - , "Parameters: " ++ show mps - , "Stationary distribution: " ++ show mf ] - -parseSubstitutionModel :: Parser S.SubstitutionModel -parseSubstitutionModel = do - n <- name - mps <- optional params - mf <- optional stationaryDistribution - let esm = assembleSubstitutionModel n mps mf - case esm of - Left err -> fail err - Right sm -> return sm - -edmModel :: [EDMComponent] -> Maybe [M.Weight] -> Parser M.MixtureModel -edmModel cs mws = do - _ <- chunk (bs "EDM") - _ <- char mmStart - n <- name - mps <- optional params - _ <- char mmEnd - let sms = map (\c -> assembleSubstitutionModel n mps (Just $ snd c)) cs - edmName = "EDM" ++ show (length cs) - ws = fromMaybe (map fst cs) mws - errs = [ e | (Left e) <- sms ] - if not $ null errs - then fail $ head errs - else return $ M.MixtureModel edmName - [ M.Component w sm | (w, Right sm) <- zip ws sms ] - -cxxModel :: Maybe [M.Weight] -> Parser M.MixtureModel -cxxModel mws = do - _ <- char (c2w 'C') - n <- decimal :: Parser Int - case cxx n mws of - Nothing -> fail "Only 10, 20, 30, 40, 50, and 60 components are supported." - Just m -> return m - -standardMixtureModel :: [M.Weight] -> Parser M.MixtureModel -standardMixtureModel ws = do - _ <- chunk (bs "MIXTURE") - _ <- char mmStart - sms <- parseSubstitutionModel `sepBy1` char separator - _ <- char mmEnd - return $ M.MixtureModel "MIXTURE" - [ M.Component w sm | (w, sm) <- zip ws sms] - -mixtureModel :: Maybe [EDMComponent] -> Maybe [M.Weight] -> Parser M.MixtureModel -mixtureModel Nothing Nothing = try (cxxModel Nothing) <|> fail "No weights provided." -mixtureModel Nothing mws@(Just ws) = try (cxxModel mws) <|> standardMixtureModel ws -mixtureModel (Just cs) mws = edmModel cs mws - --- | Parse the phylogenetic model string. The argument list is somewhat long, --- but models can have many parameters and we have to check for redundant --- parameters. --- --- @ --- getPhyloModel maybeSubstitutionModelString maybeMixtureModelString maybeEDMComponents --- @ -getPhyloModel :: Maybe String -> Maybe String -> Maybe [M.Weight] -> Maybe [EDMComponent] - -> Either String P.PhyloModel -getPhyloModel Nothing Nothing _ _ = - Left "No model was given. See help." -getPhyloModel (Just _) (Just _) _ _ = - Left "Both, substitution and mixture model string given; use only one." -getPhyloModel (Just s) Nothing Nothing Nothing = - Right $ P.SubstitutionModel $ parseStringWith "Substitution model string" parseSubstitutionModel s -getPhyloModel (Just _) Nothing (Just _) _ = - Left "Weights given; but cannot be used with substitution model." -getPhyloModel (Just _) Nothing _ (Just _) = - Left "Empirical distribution mixture model components given; but cannot be used with substitution model." -getPhyloModel Nothing (Just m) mws mcs = - Right $ P.MixtureModel $ parseStringWith "Mixture model string" (mixtureModel mcs mws) m diff --git a/src/Simulate/Simulate.hs b/src/Simulate/Simulate.hs deleted file mode 100644 --- a/src/Simulate/Simulate.hs +++ /dev/null @@ -1,192 +0,0 @@ -{-# LANGUAGE OverloadedStrings #-} -{-# LANGUAGE TemplateHaskell #-} - -{- | -Module : Simulate.Simulate -Description : Simulate multiple sequence alignments -Copyright : (c) Dominik Schrempf 2019 -License : GPL-3 - -Maintainer : dominik.schrempf@gmail.com -Stability : unstable -Portability : portable - -Creation date: Mon Jan 28 14:12:52 2019. - --} - -module Simulate.Simulate - ( simulateCmd ) -where - -import Control.Applicative ((<|>)) -import Control.Concurrent -import Control.Concurrent.Async -import Control.Monad (unless, - when) -import Control.Monad.IO.Class -import Control.Monad.Logger -import Control.Monad.Trans.Class -import Control.Monad.Trans.Reader -import qualified Data.ByteString.Lazy as L -import qualified Data.ByteString.Lazy.Char8 as LC -import Data.Maybe -import qualified Data.Set as Set -import qualified Data.Text as T -import qualified Data.Text.Lazy as LT -import qualified Data.Text.Lazy.Encoding as LT -import Data.Tree -import qualified Data.Vector.Unboxed as V -import Numeric.LinearAlgebra hiding ((<>)) -import System.Random.MWC - -import Simulate.Options -import Simulate.PhyloModel - -import ELynx.Data.Alphabet.Alphabet as A -import ELynx.Data.MarkovProcess.GammaRateHeterogeneity -import qualified ELynx.Data.MarkovProcess.MixtureModel as M -import qualified ELynx.Data.MarkovProcess.PhyloModel as P -import qualified ELynx.Data.MarkovProcess.SubstitutionModel as SM -import qualified ELynx.Data.Sequence.Alignment as A -import qualified ELynx.Data.Sequence.Sequence as Seq hiding - (name) -import ELynx.Data.Tree.MeasurableTree -import ELynx.Data.Tree.NamedTree -import ELynx.Data.Tree.Tree -import ELynx.Export.Sequence.Fasta -import ELynx.Import.MarkovProcess.EDMModelPhylobayes hiding - (Parser) -import ELynx.Import.MarkovProcess.SiteprofilesPhylobayes hiding - (Parser) -import ELynx.Import.Tree.Newick hiding (name) -import ELynx.Simulate.MarkovProcessAlongTree - -import ELynx.Tools.ByteString -import ELynx.Tools.Concurrent -import ELynx.Tools.InputOutput -import ELynx.Tools.Misc - --- Simulate a 'Alignment' for a given phylogenetic model, --- phylogenetic tree, and alignment length. -simulateAlignment :: (Measurable a, Named a) - => P.PhyloModel -> Tree a -> Int -> GenIO - -> IO A.Alignment -simulateAlignment pm t n g = do - c <- getNumCapabilities - gs <- splitGen c g - let chunks = getChunks c n - leafStatesS <- case pm of - -- TODO: This parallelization is not very intelligent, because the matrices - -- exponentiation is done in all threads. So ten threads will exponentiate - -- the same matrix ten times. - P.SubstitutionModel sm -> mapConcurrently - (\(num, gen) -> simulateAndFlatten num d e t gen) (zip chunks gs) - where d = SM.stationaryDistribution sm - e = SM.exchangeabilityMatrix sm - -- P.MixtureModel mm -> mapConcurrently - -- (\(num, gen) -> simulateAndFlattenNSitesAlongTreeMixtureModel num ws ds es t gen) (zip chunks gs) - P.MixtureModel mm -> simulateAndFlattenMixtureModelPar n ws ds es t g - where - ws = vector $ M.getWeights mm - ds = map SM.stationaryDistribution $ M.getSubstitutionModels mm - es = map SM.exchangeabilityMatrix $ M.getSubstitutionModels mm - -- XXX: The horizontal concatenation might be slow. If so, 'concatenateSeqs' - -- or 'concatenateAlignments' can be used, which directly appends vectors. - let leafStates = horizontalConcat leafStatesS - leafNames = map getName $ leaves t - code = P.getAlphabet pm - -- XXX: Probably use type safe stuff here? - alph = A.all $ alphabetSpec code - sequences = [ Seq.Sequence sName code (V.fromList $ map (`Set.elemAt` alph) ss) | - (sName, ss) <- zip leafNames leafStates ] - return $ either error id $ A.fromSequences sequences - --- Summarize EDM components; line to be printed to screen or log. -summarizeEDMComponents :: [EDMComponent] -> L.ByteString -summarizeEDMComponents cs = LC.pack - $ "Empiricial distribution mixture model with " - ++ show (length cs) ++ " components." - --- XXX. Maybe provide human readable model file. But then, why is this --- necessary. A human readable summary is reported anyways, and for Protein --- models the exchangeabilities are too many. -reportModel :: Maybe FilePath -> P.PhyloModel -> Simulate () -reportModel Nothing _ = $(logInfo) "No output file provided; omit detailed report of phylogenetic model." -reportModel (Just outFn) m = do - let modelFn = outFn <> ".model" - out "model definition (machine readable)" (bShow m <> "\n") (Just modelFn) - --- | Simulate sequences. -simulateCmd :: Maybe FilePath -> Simulate () -simulateCmd outFn = do - a <- lift ask - let treeFile = argsTreeFile a - - $(logInfo) "" - $(logInfo) $ T.pack $ "Read tree from file '" ++ treeFile ++ "'." - tree <- liftIO $ parseFileWith newick treeFile - $(logInfo) $ LT.toStrict $ LT.decodeUtf8 $ summarize tree - - let edmFile = argsEDMFile a - let sProfileFiles = argsSiteprofilesFiles a - $(logInfo) "" - $(logDebug) "Read EDM file or siteprofile files." - when (isJust edmFile && isJust sProfileFiles) $ error "Got both: --edm-file and --siteprofile-files." - edmCs <- case edmFile of - Nothing -> return Nothing - Just edmF -> do - $(logInfo) "Read EDM file." - liftIO $ Just <$> parseFileWith phylobayes edmF - maybe (return ()) ($(logInfo) . LT.toStrict . LT.decodeUtf8 . summarizeEDMComponents) edmCs - sProfiles <- case sProfileFiles of - Nothing -> return Nothing - Just fns -> do - $(logInfo) $ T.pack $ "Read siteprofiles from " ++ show (length fns) ++ " file(s)." - $(logDebug) $ T.pack $ "The file names are:" ++ show fns - xs <- liftIO $ mapM (parseFileWith siteprofiles) fns - return $ Just $ concat xs - maybe (return ()) ($(logInfo) . LT.toStrict . LT.decodeUtf8 . summarizeEDMComponents) sProfiles - let edmCsOrSiteprofiles = edmCs <|> sProfiles - - $(logInfo) "Read model string." - let ms = argsSubstitutionModelString a - mm = argsMixtureModelString a - mws = argsMixtureWeights a - eitherPhyloModel' = getPhyloModel ms mm mws edmCsOrSiteprofiles - phyloModel' <- case eitherPhyloModel' of - Left err -> lift $ error err - Right pm -> return pm - - let maybeGammaParams = argsGammaParams a - phyloModel <- case maybeGammaParams of - Nothing -> do - $(logInfo) $ LT.toStrict $ LT.decodeUtf8 - $ LC.unlines $ P.summarize phyloModel' - return phyloModel' - Just (n, alpha) -> do - $(logInfo) $ LT.toStrict $ LT.decodeUtf8 - $ LC.intercalate "\n" $ P.summarize phyloModel' ++ summarizeGammaRateHeterogeneity n alpha - return $ expand n alpha phyloModel' - - -- -- XXX: Do not report possibly huge empirical distribution mixture models - -- -- for now, because it takes too long and uses too much disk space :). - unless (isJust sProfiles) ( - do reportModel outFn phyloModel - $(logInfo) "" - ) - - $(logInfo) "Simulate alignment." - let alignmentLength = argsLength a - $(logInfo) $ T.pack $ "Length: " <> show alignmentLength <> "." - let maybeSeed = argsMaybeSeed a - gen <- case maybeSeed of - Nothing -> $(logInfo) "Seed: random" - >> liftIO createSystemRandom - Just s -> $(logInfo) (T.pack ("Seed: " <> show s <> ".")) - >> liftIO (initialize (V.fromList s)) - alignment <- liftIO $ simulateAlignment phyloModel tree alignmentLength gen - let output = (sequencesToFasta . A.toSequences) alignment - outFile = (<> ".fasta") <$> outFn - $(logInfo) "" - out "simulated multi sequence alignment" output outFile diff --git a/src/SubSample/Options.hs b/src/SubSample/Options.hs deleted file mode 100644 --- a/src/SubSample/Options.hs +++ /dev/null @@ -1,69 +0,0 @@ -{- | -Module : SubSample.Options -Description : ELynxSeq argument parsing -Copyright : (c) Dominik Schrempf 2018 -License : GPL-3 - -Maintainer : dominik.schrempf@gmail.com -Stability : unstable -Portability : portable - -Creation date: Sun Oct 7 17:29:45 2018. - --} - -module SubSample.Options - ( SubSampleArguments (..) - , SubSample - , subSampleArguments - ) where - -import Control.Applicative -import Control.Monad.Logger -import Control.Monad.Trans.Reader -import Data.Word -import Options.Applicative - -import Tools - -import ELynx.Data.Alphabet.Alphabet -import ELynx.Tools.Options - --- | Data structure holding the Command line arguments. -data SubSampleArguments = SubSampleArguments - { ssAlphabet :: Alphabet - , ssInFile :: Maybe FilePath - , ssNSites :: Int - , ssNAlignments :: Int - , ssMbSeed :: Maybe [Word32] } - --- | The sub sample command can log stuff and read necessary arguments. -type SubSample = LoggingT (ReaderT SubSampleArguments IO) - --- | Sub command parser. -subSampleArguments :: Parser SubSampleArguments -subSampleArguments = SubSampleArguments - <$> alphabetOpt - <*> optional filePathArg - <*> subSampleNSitesOpt - <*> subSampleNAlignmentsOpt - <*> seedOpt - -subSampleNSitesOpt :: Parser Int -subSampleNSitesOpt = option auto $ - long "number-of-sites" <> - short 'n' <> - metavar "INT" <> - help "Number of sites randomly drawn with replacement" - -subSampleNAlignmentsOpt :: Parser Int -subSampleNAlignmentsOpt = option auto $ - long "number-of-alignments" <> - short 'm' <> - metavar "INT" <> - help "Number of multi sequence alignments to be created" - -filePathArg :: Parser FilePath -filePathArg = strArgument $ - metavar "INPUT-FILE" <> - help "Read sequences from INPUT-FILE" diff --git a/src/SubSample/SubSample.hs b/src/SubSample/SubSample.hs deleted file mode 100644 --- a/src/SubSample/SubSample.hs +++ /dev/null @@ -1,70 +0,0 @@ -{-# LANGUAGE OverloadedStrings #-} -{-# LANGUAGE TemplateHaskell #-} - -{- | -Module : Analyze.Analyze -Description : Parse sequence file formats and analyze them -Copyright : (c) Dominik Schrempf 2018 -License : GPL-3 - -Maintainer : dominik.schrempf@gmail.com -Stability : unstable -Portability : portable - -Creation date: Fri Oct 5 08:41:05 2018. - --} - -module SubSample.SubSample - ( subSampleCmd - ) -where - -import Control.Monad -import Control.Monad.IO.Class -import Control.Monad.Logger -import Control.Monad.Trans.Class -import Control.Monad.Trans.Reader -import qualified Data.Text as T -import qualified Data.Text.Lazy as LT -import qualified Data.Text.Lazy.Builder as LT -import qualified Data.Text.Lazy.Builder.Int as LT -import qualified Data.Vector as V -import System.Random.MWC - -import SubSample.Options -import Tools - -import qualified ELynx.Data.Sequence.Alignment as M -import ELynx.Export.Sequence.Fasta -import ELynx.Tools.InputOutput - --- | Get a given number of output file names with provided suffix. --- --- > getOutFilePaths "BasePath" 11 "fasta" --- --- Will result in @BasePath.00.fasta@ up to @BasePath.10.fasta@. -getOutFilePaths :: String -> Int -> String -> [FilePath] -getOutFilePaths file n suffix = [ file ++ "." ++ digitStr i ++ "." ++ suffix - | i <- [0 .. n-1] ] - where nDigits = ceiling $ logBase (10 :: Double) (fromIntegral n) - digitStr i = T.unpack $ T.justifyRight nDigits '0' (LT.toStrict $ LT.toLazyText $ LT.decimal i) - --- | Sub sample sequences. -subSampleCmd :: Maybe FilePath -- ^ Output file base name - -> SubSample () -subSampleCmd outFileBaseName = do - SubSampleArguments al inFile nSites nAlignments seed <- lift ask - $(logInfo) "Command: Sub sample from a multi sequence alignment." - $(logInfo) $ T.pack $ " Sample " <> show nSites <> " sites." - $(logInfo) $ T.pack $ " Sample " <> show nAlignments <> " multi sequence alignments." - ss <- readSeqs al inFile - gen <- liftIO $ maybe createSystemRandom (initialize . V.fromList) seed - let a = either error id (M.fromSequences ss) - samples <- lift $ replicateM nAlignments $ M.randomSubSample nSites a gen - let results = map (sequencesToFasta . M.toSequences) samples - outFilePaths <- case outFileBaseName of - Nothing -> return $ repeat Nothing - Just fn -> return $ Just <$> getOutFilePaths fn nAlignments "fasta" - zipWithM_ (out "sub sampled multi sequence alignments") results outFilePaths - diff --git a/src/Tools.hs b/src/Tools.hs deleted file mode 100644 --- a/src/Tools.hs +++ /dev/null @@ -1,53 +0,0 @@ -{-# LANGUAGE OverloadedStrings #-} -{-# LANGUAGE TemplateHaskell #-} - -{- | -Module : Tools -Description : Common tools for sequence lynx -Copyright : (c) Dominik Schrempf 2019 -License : GPL-3 - -Maintainer : dominik.schrempf@gmail.com -Stability : unstable -Portability : portable - -Creation date: Sat Sep 7 06:24:22 2019. - --} - -module Tools - ( -- * Tools - readSeqs - -- * Options - , alphabetOpt - ) where - -import Control.Monad.IO.Class -import Control.Monad.Logger -import qualified Data.Text as T -import Options.Applicative - -import ELynx.Data.Alphabet.Alphabet -import ELynx.Data.Sequence.Sequence -import ELynx.Import.Sequence.Fasta hiding (Parser) -import ELynx.Tools.InputOutput - --- | Read sequences of given alphabet from file or standard input. -readSeqs :: (MonadIO m, MonadLogger m) => Alphabet -> Maybe FilePath -> m [Sequence] -readSeqs a mfp = do - case mfp of - Nothing -> $(logInfo) $ T.pack - $ "Read sequences from standard input; alphabet " - <> show a <> "." - Just fp -> $(logInfo) $ T.pack - $ "Read sequences from file " - <> fp <> "; alphabet " <> show a <> "." - liftIO $ parseFileOrIOWith (fasta a) mfp - --- | Command line option to specify the alphabet. Used by various commands. -alphabetOpt :: Parser Alphabet -alphabetOpt = option auto $ - long "alphabet" <> - short 'a' <> - metavar "NAME" <> - help "Specify alphabet type NAME" diff --git a/src/Translate/Options.hs b/src/Translate/Options.hs deleted file mode 100644 --- a/src/Translate/Options.hs +++ /dev/null @@ -1,71 +0,0 @@ -{- | -Module : Translate.Options -Description : ELynxSeq argument parsing -Copyright : (c) Dominik Schrempf 2018 -License : GPL-3 - -Maintainer : dominik.schrempf@gmail.com -Stability : unstable -Portability : portable - -Creation date: Sun Oct 7 17:29:45 2018. - --} - -module Translate.Options - ( TranslateArguments (..) - , Translate - , translateArguments - ) where - -import Control.Applicative -import Control.Monad.Logger -import Control.Monad.Trans.Reader -import Data.List -import Options.Applicative - -import Tools - -import ELynx.Data.Alphabet.Alphabet -import ELynx.Data.Character.Codon -import ELynx.Tools.Misc - --- | Arguments needed to translate sequences. -data TranslateArguments = TranslateArguments - { trAlphabet :: Alphabet - , trInFile :: Maybe FilePath - , trReadingFrame :: Int - , trUniversalCode :: UniversalCode } - --- | Logger and Reader type. -type Translate = LoggingT (ReaderT TranslateArguments IO) - --- | Command line parser. -translateArguments :: Parser TranslateArguments -translateArguments = TranslateArguments <$> - alphabetOpt <*> - optional inFileArg <*> - readingFrameOpt <*> - universalCodeOpt - -readingFrameOpt :: Parser Int -readingFrameOpt = option auto $ - long "reading-frame" <> - short 'r' <> - metavar "INT" <> - help "Reading frame [0|1|2]." - -universalCodeOpt :: Parser UniversalCode -universalCodeOpt = option auto $ - long "universal-code" <> - short 'u' <> - metavar "CODE" <> - help ("universal code; one of: " ++ codeStr ++ ".") - where codes = allValues :: [UniversalCode] - codeWords = map show codes - codeStr = intercalate ", " codeWords - -inFileArg :: Parser FilePath -inFileArg = strArgument $ - metavar "INPUT-FILE" <> - help "Read sequences from INPUT-FILE" diff --git a/src/Translate/Translate.hs b/src/Translate/Translate.hs deleted file mode 100644 --- a/src/Translate/Translate.hs +++ /dev/null @@ -1,51 +0,0 @@ -{-# LANGUAGE OverloadedStrings #-} -{-# LANGUAGE TemplateHaskell #-} - -{- | -Module : Analyze.Analyze -Description : Parse sequence file formats and analyze them -Copyright : (c) Dominik Schrempf 2018 -License : GPL-3 - -Maintainer : dominik.schrempf@gmail.com -Stability : unstable -Portability : portable - -Creation date: Fri Oct 5 08:41:05 2018. - --} - -module Translate.Translate - ( translateCmd - ) - where - -import Control.Monad.Logger -import Control.Monad.Trans.Class -import Control.Monad.Trans.Reader -import qualified Data.Text as T - -import Tools -import Translate.Options - -import ELynx.Data.Character.Codon -import ELynx.Data.Sequence.Sequence -import ELynx.Data.Sequence.Translate -import ELynx.Export.Sequence.Fasta -import ELynx.Tools.InputOutput - -translateSeqs :: Int -> UniversalCode -> [Sequence] -> [Sequence] -translateSeqs rf uc = map (translateSeq uc rf) - --- | Translate sequences. -translateCmd :: Maybe FilePath -> Translate () -translateCmd outFileBaseName = do - $(logInfo) "Command: Translate sequences to amino acids." - TranslateArguments al inFile rf uc <- lift ask - $(logInfo) $ T.pack $ " Universal code: " <> show uc <> "." - $(logInfo) $ T.pack $ " Reading frame: " <> show rf <> "." - $(logInfo) "" - ss <- readSeqs al inFile - let result = sequencesToFasta $ translateSeqs rf uc ss - let outFilePath = (++ ".fasta") <$> outFileBaseName - out "translated sequences" result outFilePath