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 phylogenetic trees -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,30 +14,7 @@ module Main where -import Control.Monad.Trans.Reader - -import Options - -import Compare.Compare -import Connect.Connect -import Distance.Distance -import Examine.Examine -import Shuffle.Shuffle -import Simulate.Simulate - -import ELynx.Tools.Logger -import ELynx.Tools.Options +import TLynx.TLynx main :: IO () -main = do - Arguments g c <- parseArguments - let fn = outFileBaseName g - lvl = verbosity g - lf = (++ ".log") <$> fn - case c of - Distance a -> runReaderT (eLynxWrapper lvl lf distanceDescription $ distance fn) a - Examine a -> runReaderT (eLynxWrapper lvl lf examineDescription $ examine fn) a - Simulate a -> runReaderT (eLynxWrapper lvl lf simulateDescription $ simulate fn) a - Compare a -> runReaderT (eLynxWrapper lvl lf compareDescription $ compareCmd fn) a - Connect a -> runReaderT (eLynxWrapper lvl lf connectDescription $ connectCmd fn) a - Shuffle a -> runReaderT (eLynxWrapper lvl lf shuffleDescription $ shuffleCmd fn) a +main = rTLynx diff --git a/app/Options.hs b/app/Options.hs deleted file mode 100644 --- a/app/Options.hs +++ /dev/null @@ -1,118 +0,0 @@ -{- | -Module : Options -Description : TLynx 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 - , distanceDescription - , examineDescription - , simulateDescription - , compareDescription - , connectDescription - , shuffleDescription - ) where - -import Options.Applicative -import Options.Applicative.Help.Pretty - -import Compare.Options -import Connect.Options -import Distance.Options -import Examine.Options -import Shuffle.Options -import Simulate.Options - -import ELynx.Tools.Options - -data CommandArguments = - Distance DistanceArguments - | Examine ExamineArguments - | Simulate SimulateArguments - | Compare CompareArguments - | Connect ConnectArguments - | Shuffle ShuffleArguments - -distanceDescription :: String -distanceDescription = "Compute distances between many phylogenetic trees." - -examineDescription :: String -examineDescription = "Compute summary statistics of phylogenetic trees." - -simulateDescription :: String -simulateDescription = "Simulate phylogenetic trees using birth and death processes." - -compareDescription :: String -compareDescription = "Compare two phylogenetic trees (compute distances and branch-wise differences)." - -connectDescription :: String -connectDescription = "Connect two phylogenetic trees in all ways (possibly honoring constraints)." - -shuffleDescription :: String -shuffleDescription = "Shuffle a phylogenetic tree (keep coalescent times, but shuffle topology and leaves)." - -distanceCommand :: Mod CommandFields CommandArguments -distanceCommand = command "distance" $ - info (Distance <$> distanceArguments) - ( fullDesc - <> progDesc distanceDescription - <> footerDoc (Just $ pretty distanceFooter) ) - -examineCommand :: Mod CommandFields CommandArguments -examineCommand = command "examine" $ - info (Examine <$> examineArguments) - $ progDesc examineDescription - -simulateCommand :: Mod CommandFields CommandArguments -simulateCommand = command "simulate" $ - info (Simulate <$> simulateArguments) - ( fullDesc - <> progDesc simulateDescription - <> footerDoc (Just $ pretty simulateFooter) ) - -compareCommand :: Mod CommandFields CommandArguments -compareCommand = command "compare" $ - info (Compare <$> compareArguments) - $ progDesc compareDescription - -connectCommand :: Mod CommandFields CommandArguments -connectCommand = command "connect" $ - info (Connect <$> connectArguments) - $ progDesc connectDescription - -shuffleCommand :: Mod CommandFields CommandArguments -shuffleCommand = command "shuffle" $ - info (Shuffle <$> shuffleArguments) - $ progDesc shuffleDescription - -commandArguments :: Parser CommandArguments -commandArguments = hsubparser $ - distanceCommand - <> examineCommand - <> simulateCommand - <> compareCommand - <> connectCommand - <> shuffleCommand - -parseArguments :: IO (Arguments CommandArguments) -parseArguments = parseArgumentsWith desc ftr commandArguments - -desc :: [String] -desc = [ "Compare, examine, and simulate phylogenetic trees." ] - -ftr :: [String] -ftr = "File formats:" : fs - where - toListItem = (" - " ++) - fs = map toListItem ["Newick"] diff --git a/src/Compare/Compare.hs b/src/Compare/Compare.hs deleted file mode 100644 --- a/src/Compare/Compare.hs +++ /dev/null @@ -1,199 +0,0 @@ -{-# LANGUAGE OverloadedStrings #-} -{-# LANGUAGE TemplateHaskell #-} - -{- | -Module : Compare.Compare -Description : Compare two phylogenies -Copyright : (c) Dominik Schrempf 2019 -License : GPL-3 - -Maintainer : dominik.schrempf@gmail.com -Stability : unstable -Portability : portable - -Creation date: Thu Sep 19 15:01:52 2019. - --} - -module Compare.Compare - ( compareCmd - ) 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.ByteString.Lazy.Char8 as L -import Data.List -import qualified Data.Map as M -import Data.Monoid -import qualified Data.Set as S -import qualified Data.Text as T -import qualified Data.Text.Encoding as E -import qualified Data.Text.IO as T - -import Data.Tree -import System.IO -import Text.Printf - -import Compare.Options - -import ELynx.Data.Tree.Bipartition -import ELynx.Data.Tree.BranchSupportTree as BS -import ELynx.Data.Tree.Distance -import ELynx.Data.Tree.MeasurableTree -import ELynx.Data.Tree.NamedTree -import ELynx.Data.Tree.PhyloTree -import ELynx.Data.Tree.Tree -import ELynx.Export.Tree.Newick (toNewick) -import ELynx.Import.Tree.Newick -import ELynx.Tools.InputOutput - -treesOneFile :: FilePath - -> Compare (Tree (PhyloLabel L.ByteString), Tree (PhyloLabel L.ByteString)) -treesOneFile tf = do - a <- lift ask - let nw = if argsNewickIqTree a then manyNewickIqTree else manyNewick - $(logInfo) $ T.pack $ "Parse file '" ++ tf ++ "'." - ts <- liftIO $ parseFileWith nw tf - let n = length ts - case compare n 2 of - LT -> error "Not enough trees in file." - GT -> error "Too many trees in file." - EQ -> return (head ts, head . tail $ ts) - -treesTwoFiles :: FilePath -> FilePath - -> Compare (Tree (PhyloLabel L.ByteString), Tree (PhyloLabel L.ByteString)) -treesTwoFiles tf1 tf2 = do - a <- lift ask - let nw = if argsNewickIqTree a then oneNewickIqTree else oneNewick - $(logInfo) $ T.pack $ "Parse first tree file '" ++ tf1 ++ "'." - t1 <- liftIO $ parseFileWith nw tf1 - $(logInfo) $ T.pack $ "Parse second tree file '" ++ tf2 ++ "'." - t2 <- liftIO $ parseFileWith nw tf2 - return (t1, t2) - --- | More detailed comparison of two trees. -compareCmd :: Maybe FilePath -> Compare () -compareCmd outFile = do - a <- lift ask - -- Determine output handle (stdout or file). - let outFn = (++ ".out") <$> outFile - outH <- outHandle "results" outFn - -- liftIO $ hPutStrLn outH "" - - -- Read input. - let inFiles = argsInFiles a - nFiles = length inFiles - (t1, t2) <- case nFiles of - 1 -> treesOneFile (head inFiles) - 2 -> treesTwoFiles (head inFiles) (head . tail $ inFiles) - _ -> error "Need two input files with one tree each or one input file with two trees." - -- liftIO $ hPutStrLn outH "" - - liftIO $ hPutStrLn outH "Tree 1:" - liftIO $ L.hPutStrLn outH $ toNewick t1 - -- liftIO $ hPrint outH t1 - liftIO $ hPutStrLn outH "Tree 2:" - liftIO $ L.hPutStrLn outH $ toNewick t2 - -- liftIO $ hPrint outH t2 - liftIO $ hPutStrLn outH "" - - -- Check input. - let lvs1 = leaves t1 - lvs2 = leaves t2 - lfns1 = map getName lvs1 - lfns2 = map getName lvs2 - s1 = S.fromList lfns1 - s2 = S.fromList lfns2 - if s1 == s2 - then liftIO $ hPutStrLn outH "Trees have the same set of leaf names." - else liftIO $ hPutStrLn outH "Trees do not have the same set of leaf names." - liftIO $ hPutStrLn outH "" - - -- Distances. - let formatD str val = T.justifyLeft 25 ' ' str <> " " <> val - liftIO $ hPutStrLn outH "Distances." - liftIO $ T.hPutStrLn outH $ formatD "Symmetric" - (T.pack $ show $ symmetric t1 t2) - liftIO $ T.hPutStrLn outH $ formatD "Branch score" - (T.pack $ show $ branchScoreWith getName getLen t1 t2) - let t1' = BS.normalize t1 - t2' = BS.normalize t2 - $(logDebug) "Trees with normalized branch support values:" - $(logDebug) $ E.decodeUtf8 $ L.toStrict $ toNewick t1' - $(logDebug) $ E.decodeUtf8 $ L.toStrict $ toNewick t2' - liftIO $ T.hPutStrLn outH $ formatD "Incompatible split (0.10)" - (T.pack $ show $ incompatibleSplits (collapse 0.1 t1') (collapse 0.1 t2')) - liftIO $ T.hPutStrLn outH $ formatD "Incompatible split (0.50)" - (T.pack $ show $ incompatibleSplits (collapse 0.5 t1') (collapse 0.5 t2')) - -- liftIO $ T.hPutStrLn outH $ formatD "Incompatible split (0.60)" - -- (T.pack $ show $ incompatibleSplits (collapse 0.6 t1') (collapse 0.6 t2')) - -- liftIO $ T.hPutStrLn outH $ formatD "Incompatible split (0.70)" - -- (T.pack $ show $ incompatibleSplits (collapse 0.7 t1') (collapse 0.7 t2')) - liftIO $ T.hPutStrLn outH $ formatD "Incompatible split (0.80)" - (T.pack $ show $ incompatibleSplits (collapse 0.8 t1') (collapse 0.8 t2')) - liftIO $ T.hPutStrLn outH $ formatD "Incompatible split (0.90)" - (T.pack $ show $ incompatibleSplits (collapse 0.9 t1') (collapse 0.9 t2')) - -- liftIO $ T.hPutStrLn outH $ formatD "Incompatible split (1.01)" - -- (T.pack $ show $ incompatibleSplits (collapse 1.01 t1') (collapse 1.01 t2')) - -- liftIO $ L.hPutStrLn outH $ toNewick (collapse 1.01 t1') - - -- Bipartitions. - when (argsBipartitions a) (do - let bp1 = bipartitions (fmap getName t1) - bp2 = bipartitions (fmap getName t2) - bp1Only = bp1 S.\\ bp2 - bp2Only = bp2 S.\\ bp1 - unless (S.null bp1Only) - (do - liftIO $ hPutStrLn outH "" - liftIO $ hPutStrLn outH "Bipartitions in Tree 1 that are not in Tree 2." - -- let bp1Strs = map (bphuman L.unpack . bpmap getName) (S.toList bp1Only) - forM_ bp1Only (liftIO . hPutStrLn outH . bphuman L.unpack)) - -- let bp1Strs = map (bphuman L.unpack) (S.toList bp1Only) - -- liftIO $ hPutStrLn outH $ intercalate "\n" bp1Strs) - unless (S.null bp2Only) - (do - liftIO $ hPutStrLn outH "" - liftIO $ hPutStrLn outH "Bipartitions in Tree 2 that are not in Tree 1." - forM_ bp2Only (liftIO . hPutStrLn outH . bphuman L.unpack)) - - -- Common bipartitions and their respective differences in branch lengths. - liftIO $ hPutStrLn outH "" - let bpCommon = bp1 `S.intersection` bp2 - if S.null bpCommon - then liftIO $ hPutStrLn outH "There are no common bipartitions." - else do - let bpToBrLen1 = M.map getSum $ bipartitionToBranchLength getName (Sum . getLen) t1 - bpToBrLen2 = M.map getSum $ bipartitionToBranchLength getName (Sum . getLen) t2 - liftIO $ hPutStrLn outH "Common bipartitions and their respective differences in branch lengths." - -- Header. - liftIO $ hPutStrLn outH header - forM_ bpCommon (liftIO . hPutStrLn outH . getCommonBpStr L.unpack bpToBrLen1 bpToBrLen2)) - - liftIO $ hClose outH - -- TODO. Plot using haskell-chart. - -- See https://github.com/timbod7/haskell-chart/wiki. - -header :: String -header = intercalate " " $ cols ++ ["Bipartition"] - where cols = map (T.unpack . T.justifyRight 12 ' ') - [ "Length 1", "Length 2", "Delta", "Relative [%]"] - -getCommonBpStr :: (Ord a, Fractional b, PrintfArg b) - => (a -> String) - -> M.Map (Bipartition a) b -> M.Map (Bipartition a) b - -> Bipartition a -> String -getCommonBpStr f m1 m2 p = intercalate " " - [ printf "% 12.7f" l1 - , printf "% 12.7f" l2 - , printf "% 12.7f" d - , printf "% 12.7f" rd - , s ] - where l1 = m1 M.! p - l2 = m2 M.! p - d = l1 - l2 - rd = 2 * d / (l1 + l2) - s = bphuman f p diff --git a/src/Compare/Options.hs b/src/Compare/Options.hs deleted file mode 100644 --- a/src/Compare/Options.hs +++ /dev/null @@ -1,65 +0,0 @@ -{- | -Module : Compare.Options -Description : Options for the compare subcommand -Copyright : (c) Dominik Schrempf 2019 -License : GPL-3 - -Maintainer : dominik.schrempf@gmail.com -Stability : unstable -Portability : portable - -Creation date: Thu Sep 19 15:02:21 2019. - --} - -module Compare.Options - ( CompareArguments (..) - , Compare - , compareArguments - ) where - -import Control.Monad.Logger -import Control.Monad.Trans.Reader -import Options.Applicative - --- | Arguments of compare command. -data CompareArguments = CompareArguments - { argsNormalize :: Bool - , argsBipartitions :: Bool - , argsNewickIqTree :: Bool - , argsInFiles :: [FilePath] } - - --- | Logger and reader data type. -type Compare = LoggingT (ReaderT CompareArguments IO) - --- | Parse arguments of compare command. -compareArguments :: Parser CompareArguments -compareArguments = CompareArguments <$> - normalize - <*> bipartitions - <*> newickIqTree - <*> some file - -normalize :: Parser Bool -normalize = switch $ - long "normalize" - <> short 'n' - <> help "Normalize trees before comparison" - -bipartitions :: Parser Bool -bipartitions = switch $ - long "bipartitions" - <> short 'b' - <> help "Print common and missing bipartitions" - -file :: Parser FilePath -file = strArgument $ - metavar "NAME" - <> help "Tree file" - -newickIqTree :: Parser Bool -newickIqTree = switch $ - long "newick-iqtree" - <> short 'i' - <> help "Use IQ-TREE Newick format (internal node labels are branch support values)" diff --git a/src/Connect/Connect.hs b/src/Connect/Connect.hs deleted file mode 100644 --- a/src/Connect/Connect.hs +++ /dev/null @@ -1,107 +0,0 @@ -{-# LANGUAGE OverloadedStrings #-} -{-# LANGUAGE TemplateHaskell #-} - -{- | -Module : Connect.Connect -Description : Connect two phylogenies -Copyright : (c) Dominik Schrempf 2019 -License : GPL-3 - -Maintainer : dominik.schrempf@gmail.com -Stability : unstable -Portability : portable - -Creation date: Thu Sep 19 15:01:52 2019. - --} - -module Connect.Connect - ( connectCmd - ) where - -import Control.Monad.IO.Class -import Control.Monad.Logger -import Control.Monad.Trans.Class -import Control.Monad.Trans.Reader -import qualified Data.ByteString.Lazy.Char8 as L -import Data.Tree -import System.IO - -import Connect.Options - -import ELynx.Data.Tree.Bipartition (bipartition, compatible) -import ELynx.Data.Tree.NamedTree (getName) -import ELynx.Data.Tree.PhyloTree (PhyloLabel (PhyloLabel)) -import ELynx.Data.Tree.Subset (Subset, smap) -import ELynx.Data.Tree.Tree (clades, connect) -import ELynx.Export.Tree.Newick (toNewick) -import ELynx.Import.Tree.Newick (manyNewick, oneNewick) -import ELynx.Tools.InputOutput (outHandle, parseFileWith) -import ELynx.Tools.Text (fromBs, tShow) - --- TODO: Write a proper documentation. - --- | Connect two trees honoring possible constraints. -connectCmd :: Maybe FilePath -> Connect () -connectCmd outFile = do - -- Determine output handle (stdout or file). - a <- lift ask - let outFn = (++ ".out") <$> outFile - outH <- outHandle "results" outFn - - -- Do we have constraints or not? - let cs = constraints a - l = inFileA a - r = inFileB a - case cs of - Nothing -> connectOnly outH l r - Just c -> connectAndFilter outH c l r - - liftIO $ hClose outH - -connectTrees :: Tree (PhyloLabel L.ByteString) - -> Tree (PhyloLabel L.ByteString) - -> [Tree (PhyloLabel L.ByteString)] -connectTrees = connect (PhyloLabel "" Nothing Nothing) - -type Constraint a = Subset a - -compatibleAll :: (Show a, Ord a) => Tree a -> [Constraint a] -> Bool -compatibleAll (Node _ [l, r]) cs = all (compatible (bipartition l)) cs && - all (compatible (bipartition r)) cs -compatibleAll _ _ = error "Tree is not bifurcating." - -compatibleWith :: (Show b, Ord b) => (a -> b) -> [Constraint a] -> Tree a -> Bool -compatibleWith f cs t = compatibleAll (fmap f t) (map (smap f) cs) - -parseTrees :: FilePath -> FilePath - -> Connect (Tree (PhyloLabel L.ByteString), Tree (PhyloLabel L.ByteString)) -parseTrees l r = do - tl <- liftIO $ parseFileWith oneNewick l - tr <- liftIO $ parseFileWith oneNewick r - $(logInfo) "Tree 1:" - $(logInfo) $ fromBs $ toNewick tl - $(logInfo) "Tree 2:" - $(logInfo) $ fromBs $ toNewick tr - return (tl, tr) - -connectOnly :: Handle -> FilePath -> FilePath -> Connect () -connectOnly h l r = do - (tl, tr) <- parseTrees l r - let ts = connectTrees tl tr - $(logInfo) $ "Connected trees: " <> tShow (length ts) - liftIO $ L.hPutStr h $ L.unlines $ map toNewick ts - -connectAndFilter :: Handle -> FilePath -> FilePath -> FilePath -> Connect () -connectAndFilter h c l r = do - cts <- liftIO $ parseFileWith manyNewick c - $(logInfo) "Constraints:" - $(logInfo) $ fromBs $ L.intercalate "\n" $ map toNewick cts - (tl, tr) <- parseTrees l r - let ts = connectTrees tl tr - cs = concatMap clades cts :: [Constraint (PhyloLabel L.ByteString)] - -- Only collect trees that are compatible with the constraints. - ts' = filter (compatibleWith getName cs) ts - $(logInfo) $ "Connected trees: " <> tShow (length ts) - $(logInfo) $ "Compatible trees: " <> tShow (length ts') - liftIO $ L.hPutStr h $ L.unlines $ map toNewick ts' diff --git a/src/Connect/Options.hs b/src/Connect/Options.hs deleted file mode 100644 --- a/src/Connect/Options.hs +++ /dev/null @@ -1,64 +0,0 @@ -{- | -Module : Connect.Options -Description : Options for the connect subcommand -Copyright : (c) Dominik Schrempf 2019 -License : GPL-3 - -Maintainer : dominik.schrempf@gmail.com -Stability : unstable -Portability : portable - -Creation date: Thu Sep 19 15:02:21 2019. - --} - -module Connect.Options - ( ConnectArguments (..) - , Connect - , connectArguments - ) where - -import Control.Monad.Logger -import Control.Monad.Trans.Reader -import Options.Applicative - --- | Arguments of connect command. -data ConnectArguments = ConnectArguments - { newickIqTreeFlag :: Bool - , constraints :: Maybe FilePath - , inFileA :: FilePath - , inFileB :: FilePath } - --- | Logger and reader data type. -type Connect = LoggingT (ReaderT ConnectArguments IO) - --- | Parse arguments of connect command. -connectArguments :: Parser ConnectArguments -connectArguments = ConnectArguments - <$> newickIqTree - <*> constraintsFile - <*> fileA - <*> fileB - -newickIqTree :: Parser Bool -newickIqTree = switch $ - long "newick-iqtree" - <> short 'i' - <> help "Use IQ-TREE Newick format (internal node labels are branch support values)" - -constraintsFile :: Parser (Maybe FilePath) -constraintsFile = optional $ strOption $ - metavar "CONSTRAINTS" - <> short 'c' - <> long "contraints" - <> help "File containing one or more Newick trees to be used as constraints" - -fileA :: Parser FilePath -fileA = strArgument $ - metavar "TREE-FILE-A" - <> help "File containing the first Newick tree" - -fileB :: Parser FilePath -fileB = strArgument $ - metavar "TREE-FILE-B" - <> help "File containing the second Newick tree" diff --git a/src/Distance/Distance.hs b/src/Distance/Distance.hs deleted file mode 100644 --- a/src/Distance/Distance.hs +++ /dev/null @@ -1,157 +0,0 @@ -{-# LANGUAGE OverloadedStrings #-} -{-# LANGUAGE TemplateHaskell #-} - -{- | -Description : Compute distances between trees -Copyright : (c) Dominik Schrempf 2019 -License : GPL-3 - -Maintainer : dominik.schrempf@gmail.com -Stability : unstable -Portability : portable - -Creation date: Wed May 29 18:09:39 2019. - -- Symmetric (Robinson-Foulds) distance. -- Incompatible splits distance. - --} - -module Distance.Distance - ( distance - ) -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.ByteString.Lazy.Char8 as L -import Data.Maybe -import qualified Data.Text as T -import qualified Data.Text.IO 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 Statistics.Sample -import System.IO -import Text.Printf - -import Distance.Options - -import ELynx.Data.Tree.BranchSupportTree as B -import ELynx.Data.Tree.Distance -import ELynx.Data.Tree.MeasurableTree as M -import ELynx.Data.Tree.PhyloTree -import ELynx.Export.Tree.Newick -import ELynx.Import.Tree.Newick -import ELynx.Tools.ByteString (alignLeft, alignRight) -import ELynx.Tools.InputOutput -import ELynx.Tools.Logger - -pf :: String -pf = "%.3f" - -header :: Int -> DistanceMeasure -> L.ByteString -header n d = alignLeft (n+2) "Tree 1" - <> alignLeft (n+2) "Tree 2" - <> alignRight 20 (L.pack $ show d) - -showTriplet :: (PrintfArg a) => Int -> [String] -> (Int, Int, a) -> L.ByteString -showTriplet n args (i, j, d) = i' <> j' <> d' - where i' = alignLeft (n+2) $ L.pack (args !! i) - j' = alignLeft (n+2) $ L.pack (args !! j) - d' = alignRight 20 $ L.pack (printf pf d) - --- | Compute distance functions between phylogenetic trees. -distance :: Maybe FilePath -> Distance () -distance outFileBN = do - a <- lift ask - let nw = if argsNewickIqTree a then manyNewickIqTree else manyNewick - -- Determine output handle (stdout or file). - let outFile = (++ ".out") <$> outFileBN - outH <- outHandle "results" outFile - -- Master tree (in case it is given). - let mname = argsMasterTreeFile a - mtree <- case mname of - Nothing -> return Nothing - Just f -> do $(logInfo) $ T.pack $ "Read master tree from file: " <> f <> "." - ts <- liftIO $ parseFileWith nw f - let n = length ts - when (n > 1) (error "More than one tree found in master file.") - $(logInfo) "Compute distances between all trees and master tree." - $(logInfo) $ T.pack $ "Trees are numbered from 0 to " ++ show (n-1) ++ "." - return $ Just (head ts) - let tfps = argsInFiles a - (trees, names) <- - if length tfps <= 1 - then - do ts <- if null tfps - then do $(logInfo) "Read trees from standard input." - liftIO $ parseIOWith nw - else do let f = head tfps - $(logInfo) $ T.pack $ "Read trees from file: " <> f <> "." - liftIO $ parseFileWith nw f - let n = length ts - when (n < 1) (error "Not enough trees found in file.") - when (isNothing mtree) $ $(logInfo) "Compute pairwise distances between trees in the same file." - $(logInfo) $ T.pack $ "Trees are numbered from 0 to " ++ show (n-1) ++ "." - return (ts, take n (map show [0 :: Int ..])) - else - do $(logInfo) "Read trees from files." - ts <- liftIO $ mapM (parseFileWith newick) tfps - when (length ts <= 1) (error "Not enough trees found in files.") - when (isNothing mtree) $ $(logInfo) "Compute pairwise distances between trees from different files." - $(logInfo) "Trees are named according to their file names." - return (ts, tfps) - $(logDebug) "The trees are:" - $(logDebug) $ LT.toStrict $ LT.decodeUtf8 $ L.unlines $ map toNewick trees - case outFile of - Nothing -> logNewSection "Write results to standard output." - Just f -> logNewSection $ T.pack $ "Write results to file " <> f <> "." - let n = maximum $ 6 : map length names - dist = argsDistance a - case dist of - Symmetric -> $(logInfo) "Use symmetric (Robinson-Foulds) distance." - IncompatibleSplit val -> do - $(logInfo) "Use incompatible split distance." - $(logInfo) $ T.pack $ "Collapse nodes with support less than " ++ show val ++ "." - BranchScore -> $(logInfo) "Use branch score distance." - when (argsNormalize a) $ $(logInfo) "Normalize trees before calculation of distances." - let distanceMeasure :: Tree (PhyloLabel L.ByteString) -> Tree (PhyloLabel L.ByteString) -> Double - distanceMeasure = case dist of - Symmetric -> \t1 t2 -> fromIntegral $ symmetric t1 t2 - IncompatibleSplit _ -> \t1 t2 -> fromIntegral $ incompatibleSplits t1 t2 - BranchScore -> branchScore - normalizeF = if argsNormalize a then M.normalize else id - collapseF = case dist of - -- For the incompatible split distance we have to collapse branches with - -- support lower than the given value. Before doing so, we normalize the - -- branch support values. - IncompatibleSplit val -> collapse val . B.normalize - _ -> id - trees' = map (collapseF . normalizeF) trees - $(logDebug) "The prepared trees are:" - $(logDebug) $ LT.toStrict $ LT.decodeUtf8 $ L.unlines $ map toNewick trees' - let dsTriplets = case mtree of - Nothing -> pairwise distanceMeasure trees' - Just t -> [ (0, i, distanceMeasure t t') | ( i, t' ) <- zip [1..] trees' ] - ds = map (\(_, _, x) -> x) dsTriplets - dsVec = V.fromList ds - liftIO $ hPutStrLn outH $ "Summary statistics of " ++ show dist ++ " Distance:" - liftIO $ T.hPutStrLn outH $ T.justifyLeft 10 ' ' "Mean: " <> T.pack (printf pf (mean dsVec)) - liftIO $ T.hPutStrLn outH $ T.justifyLeft 10 ' ' "Variance: " <> T.pack (printf pf (variance dsVec)) - -- L.putStrLn $ L.unlines $ map toNewick ts - -- L.putStrLn $ L.unlines $ map toNewick tsN - -- L.putStrLn $ L.unlines $ map toNewick tsC - lift $ unless (argsSummaryStatistics a) ( - do - lift $ hPutStrLn outH "" - lift $ L.hPutStrLn outH $ header n dist - case mname of - Nothing -> lift $ L.hPutStr outH $ L.unlines (map (showTriplet n names) dsTriplets) - Just mn -> lift $ L.hPutStr outH $ L.unlines (map (showTriplet n (mn : names)) dsTriplets) - ) - liftIO $ hClose outH diff --git a/src/Distance/Options.hs b/src/Distance/Options.hs deleted file mode 100644 --- a/src/Distance/Options.hs +++ /dev/null @@ -1,147 +0,0 @@ -{- | -Module : Distance.Options -Description : Options of tree-dist -Copyright : (c) Dominik Schrempf 2019 -License : GPL-3 - -Maintainer : dominik.schrempf@gmail.com -Stability : unstable -Portability : portable - -Creation date: Thu Aug 29 13:02:22 2019. - --} - -module Distance.Options - ( DistanceArguments (..) - , DistanceMeasure (..) - , Distance - , distanceArguments - , distanceFooter - ) where - -import Control.Monad.Logger -import Control.Monad.Trans.Reader -import Data.List -import Data.Void -import Options.Applicative -import Text.Megaparsec (Parsec, eof, try) -import Text.Megaparsec.Char (char, string) -import Text.Megaparsec.Char.Lexer (float) -import Text.Printf - -import ELynx.Tools.Options - --- | Supported distance measures. -data DistanceMeasure = - Symmetric -- ^ Symmetric distance. - | IncompatibleSplit Double -- ^ Incompatible split distance; collapse nodes - -- with branch support below given value. - | BranchScore -- ^ Branch score distance. - -instance Show DistanceMeasure where - show Symmetric = "Symmetric" - show (IncompatibleSplit c) = "Incompatible Split (" ++ printf "%.1f" c ++ ")" - show BranchScore = "Branch Score" - --- | Arguments needed to compute distance measures. -data DistanceArguments = DistanceArguments - { argsDistance :: DistanceMeasure - , argsNormalize :: Bool - , argsSummaryStatistics :: Bool - , argsMasterTreeFile :: Maybe FilePath - , argsNewickIqTree :: Bool - , argsInFiles :: [FilePath] - } - --- | Logger and reader data type. -type Distance = LoggingT (ReaderT DistanceArguments IO) - --- | COmmand line parser. -distanceArguments :: Parser DistanceArguments -distanceArguments = DistanceArguments <$> - distanceOpt - <*> normalizeSwitch - <*> summaryStatisticsSwitch - <*> masterTreeFile - <*> newickIqTree - <*> many inFilesArg - -masterTreeFile :: Parser (Maybe FilePath) -masterTreeFile = optional $ strOption $ - long "master-tree-file" <> - short 'm' <> - metavar "MASTER-TREE-File" <> - help "Compare all trees to the tree in the master tree file." - -inFilesArg :: Parser FilePath -inFilesArg = strArgument $ - metavar "INPUT-FILES" <> - help "Read tree(s) from INPUT-FILES; if more files are given, one tree is expected per file" - -symmetric :: Parsec Void String DistanceMeasure -symmetric = do - _ <- string "symmetric" - _ <- eof - pure Symmetric - -incompatibleSplit :: Parsec Void String DistanceMeasure -incompatibleSplit = do - _ <- string "incompatible-split" - _ <- char '[' - f <- float - _ <- char ']' - _ <- eof - if (0 < f) && (f < 1) - then pure $ IncompatibleSplit f - else error "Branch support has to be between 0 and 1." - -branchScore :: Parsec Void String DistanceMeasure -branchScore = do - _ <- string "branch-score" - _ <- eof - pure BranchScore - -distanceParser :: Parsec Void String DistanceMeasure -distanceParser = try symmetric - <|> try incompatibleSplit - -- Try first the normalized one, since the normal branch score - -- parser also succeeds in this case. - <|> branchScore - -distanceOpt :: Parser DistanceMeasure -distanceOpt = option (megaReadM distanceParser) $ - long "distance" <> - short 'd' <> - metavar "MEASURE" <> - help "Type of distance to calculate (available distance measures are listed below)" - -summaryStatisticsSwitch :: Parser Bool -summaryStatisticsSwitch = switch $ - long "summary-statistics" <> - short 's' <> - help "Report summary statistics only" - -normalizeSwitch :: Parser Bool -normalizeSwitch = switch $ - long "normalize" <> - short 'n' <> - help "Normalize trees before distance calculation; only affect distances depending on branch lengths" - -newickIqTree :: Parser Bool -newickIqTree = switch $ - long "newick-iqtree" - <> short 'i' - <> help "Use IQ-TREE Newick format (internal node labels are branch support values)" - --- | Information about provided distance types. -distanceFooter :: String -distanceFooter = intercalate "\n" - [ "Available distance measures:" - , " symmetric Symmetric distance (Robinson-Foulds distance)." - , " incompatible-split[VAL] Incompatible split distance. Collapse branches" - , " with support less than VAL before distance calculation;" - , " in this way, only well supported difference contribute" - , " to the distance measure." - , " branch-score Branch score distance." - ] diff --git a/src/Examine/Examine.hs b/src/Examine/Examine.hs deleted file mode 100644 --- a/src/Examine/Examine.hs +++ /dev/null @@ -1,67 +0,0 @@ -{-# LANGUAGE OverloadedStrings #-} -{-# LANGUAGE TemplateHaskell #-} - -{- | -Description : Analyze trees -Copyright : (c) Dominik Schrempf 2019 -License : GPL-3 - -Maintainer : dominik.schrempf@gmail.com -Stability : unstable -Portability : portable - -Creation date: Fri May 24 13:47:56 2019. - --} - -module Examine.Examine - ( examine - ) -where - -import Control.Monad (unless) -import Control.Monad.IO.Class -import Control.Monad.Logger -import Control.Monad.Trans.Class -import Control.Monad.Trans.Reader -import qualified Data.ByteString.Lazy.Char8 as L -import Data.List (nub, (\\)) -import qualified Data.Text as T -import Data.Tree -import System.IO (Handle, hPutStrLn) - -import Examine.Options - -import ELynx.Data.Tree.MeasurableTree -import ELynx.Data.Tree.NamedTree -import ELynx.Data.Tree.PhyloTree -import ELynx.Data.Tree.Tree -import ELynx.Import.Tree.Newick -import ELynx.Tools.InputOutput - -readTrees :: Maybe FilePath -> Examine [Tree (PhyloLabel L.ByteString)] -readTrees mfp = do - case mfp of - Nothing -> $(logInfo) "Read tree(s) from standard input." - Just fp -> $(logInfo) $ T.pack $ "Read tree(s) from file " <> fp <> "." - a <- lift ask - let nw = if argsNewickIqTree a then manyNewickIqTree else manyNewick - liftIO $ parseFileOrIOWith nw mfp - -examineTree :: (Measurable a, Named a) - => Handle -> Tree a -> IO () -examineTree h t = do - L.hPutStrLn h $ summarize t - unless (null dups) (hPutStrLn h $ "Duplicate leaves: " ++ show dups) - where lvs = map getName $ leaves t - dups = lvs \\ nub lvs - --- | Examine phylogenetic trees. -examine :: Maybe FilePath -> Examine () -examine outFn = do - a <- lift ask - let inFn = argsInFile a - trs <- readTrees inFn - let outFilePath = (++ ".out") <$> outFn - outH <- outHandle "results" outFilePath - liftIO $ mapM_ (examineTree outH) trs diff --git a/src/Examine/Options.hs b/src/Examine/Options.hs deleted file mode 100644 --- a/src/Examine/Options.hs +++ /dev/null @@ -1,48 +0,0 @@ -{- | -Module : Examine.Options -Description : Tree analysis options -Copyright : (c) Dominik Schrempf 2019 -License : GPL-3 - -Maintainer : dominik.schrempf@gmail.com -Stability : unstable -Portability : portable - -Creation date: Thu Aug 29 08:16:45 2019. - --} - -module Examine.Options - ( ExamineArguments (..) - , Examine - , examineArguments - ) where - -import Control.Monad.Logger -import Control.Monad.Trans.Reader -import Options.Applicative - --- | Arguments needed to examine phylogenetic trees. -data ExamineArguments = ExamineArguments - { argsInFile :: Maybe FilePath - , argsNewickIqTree :: Bool } - --- | Logger and reader data type. -type Examine = LoggingT (ReaderT ExamineArguments IO) - --- | Command line parser. -examineArguments :: Parser ExamineArguments -examineArguments = ExamineArguments <$> - optional inFile - <*> newickIqTree - -inFile :: Parser FilePath -inFile = strArgument $ - metavar "INPUT-FILE" <> - help "Read trees from INPUT-FILE" - -newickIqTree :: Parser Bool -newickIqTree = switch $ - long "newick-iqtree" - <> short 'i' - <> help "Use IQ-TREE Newick format (internal node labels are branch support values)" diff --git a/src/Shuffle/Options.hs b/src/Shuffle/Options.hs deleted file mode 100644 --- a/src/Shuffle/Options.hs +++ /dev/null @@ -1,68 +0,0 @@ -{- | -Module : Shuffle.Options -Description : Options for the connect subcommand -Copyright : (c) Dominik Schrempf 2019 -License : GPL-3 - -Maintainer : dominik.schrempf@gmail.com -Stability : unstable -Portability : portable - -Creation date: Thu Sep 19 15:02:21 2019. - --} - -module Shuffle.Options - ( ShuffleArguments (..) - , Shuffle - , shuffleArguments - ) where - -import Control.Monad.Logger -import Control.Monad.Trans.Reader --- import Data.Word -import Options.Applicative - --- import ELynx.Tools.Options - --- | Arguments of shuffle command. -data ShuffleArguments = ShuffleArguments - { newickIqTreeFlag :: Bool - , nReplicates :: Int - , inFile :: FilePath } - -- -- TODO. Too complicated at the moment with MWC, use MonadRandom, but then - -- -- seed option would not match other commands, so leave the option out for - -- -- now. Maybe, when working on this, try using splitmix consistently. - -- - -- -- | Seed of NRG, random if 'Nothing'. - -- , argsSeed :: Maybe [Word32] } - --- | Logger and reader data type. -type Shuffle = LoggingT (ReaderT ShuffleArguments IO) - --- | Parse arguments of shuffle command. -shuffleArguments :: Parser ShuffleArguments -shuffleArguments = ShuffleArguments - <$> newickIqTree - <*> n - <*> file - -- <*> seedOpt - -newickIqTree :: Parser Bool -newickIqTree = switch $ - long "newick-iqtree" - <> short 'i' - <> help "Use IQ-TREE Newick format (internal node labels are branch support values)" - -n :: Parser Int -n = option auto $ - long "replicates" - <> short 'n' - <> metavar "N" - <> value 1 - <> help "Number of trees to generate" - -file :: Parser FilePath -file = strArgument $ - metavar "TREE-FILE" - <> help "File containing a Newick tree" diff --git a/src/Shuffle/Shuffle.hs b/src/Shuffle/Shuffle.hs deleted file mode 100644 --- a/src/Shuffle/Shuffle.hs +++ /dev/null @@ -1,142 +0,0 @@ -{-# LANGUAGE OverloadedStrings #-} -{-# LANGUAGE TemplateHaskell #-} - -{- | -Module : Shuffle.Shuffle -Description : Shuffle a phylogeny. -Copyright : (c) Dominik Schrempf 2019 -License : GPL-3 - -Maintainer : dominik.schrempf@gmail.com -Stability : unstable -Portability : portable - -Creation date: Thu Sep 19 15:01:52 2019. - -The coalescent times are unaffected. The topology and the leaf order is -shuffled. Branch support values are ignored and lost. - --} - -module Shuffle.Shuffle - ( shuffleCmd - ) where - -import Control.Monad.IO.Class -import Control.Monad.Logger -import Control.Monad.Random -import Control.Monad.Trans.Reader -import Data.Array (elems) -import Data.Array.ST (newListArray, readArray, - runSTArray, writeArray) -import qualified Data.ByteString.Lazy.Char8 as L -import Data.List (filter) -import Data.Maybe (isNothing) -import Data.Tree (Tree (Node), flatten, - rootLabel) -import System.IO - -import Shuffle.Options - -import ELynx.Data.Tree.MeasurableTree (distancesOriginLeaves, height, - rootHeight) -import ELynx.Data.Tree.NamedTree (getName) -import ELynx.Data.Tree.PhyloTree (PhyloLabel, brLen) -import ELynx.Data.Tree.Tree (leaves) -import ELynx.Export.Tree.Newick (toNewick) -import ELynx.Import.Tree.Newick (oneNewick) -import ELynx.Simulate.PointProcess (PointProcess (PointProcess), - toReconstructedTree) -import ELynx.Tools.Definitions (eps) -import ELynx.Tools.InputOutput (outHandle, parseFileWith) -import ELynx.Tools.Text (fromBs, tShow) - --- | Shuffle a tree. Get all coalescent times, shuffle them. Get all leaves, --- shuffle them. Connect the shuffled leaves with the shuffled coalescent times. --- The shuffled tree has a new topology while keeping the same set of coalescent --- times and leaves. -shuffleCmd :: Maybe FilePath -> Shuffle () -shuffleCmd outFile = do - a <- lift ask - - let outFn = (++ ".tree") <$> outFile - h <- outHandle "results" outFn - - t <- liftIO $ parseFileWith oneNewick (inFile a) - $(logInfo) "Input tree:" - $(logInfo) $ fromBs $ toNewick t - - -- Check if all branches have a given length. However, the length of the stem is not important. - let r = rootLabel t - r' = r {brLen = Just 0} - t' = t {rootLabel = r'} - when (isNothing $ traverse brLen t') - (do - $(logDebug) $ tShow t' - error "Not all branches have a given length.") - - -- Check if tree is ultrametric enough. - let dh = sum $ map (height t -) (distancesOriginLeaves t) - $(logDebug) $ "Distance in branch length to being ultrametric: " <> tShow dh - when (dh > 2e-4) - (error "Tree is not ultrametric.") - when (dh > eps && dh < 2e-4) $ - $(logInfo) "Tree is nearly ultrametric, ignore branch length differences smaller than 2e-4." - when (dh < eps) $ - $(logInfo) "Tree is ultrametric." - - let cs = filter (>0) $ flatten $ mapTree rootHeight t - ls = map getName $ leaves t - $(logDebug) $ "Number of coalescent times: " <> tShow (length cs) - $(logDebug) $ "Number of leaves: " <> tShow (length ls) - $(logDebug) "The coalescent times are: " - $(logDebug) $ tShow cs - - ts <- liftIO $ shuffle (nReplicates a) (height t) cs ls - liftIO $ L.hPutStr h $ L.unlines $ map toNewick ts - - liftIO $ hClose h - -shuffle :: MonadRandom m - => Int -- How many? - -> Double -- Stem length. - -> [Double] -- Coalescent times. - -> [L.ByteString] -- Leave names. - -> m [Tree (PhyloLabel L.ByteString)] -shuffle n o cs ls = do - css <- grabble cs n (length cs) - lss <- grabble ls n (length ls) - return [toReconstructedTree "" (PointProcess names times o) | (times, names) <- zip css lss] - --- TODO: This seems to be 'extend' from comonad! -mapTree :: (Tree a -> b) -> Tree a -> Tree b -mapTree f t@(Node _ cs) = Node (f t) (map (mapTree f) cs) - --- | From https://wiki.haskell.org/Random_shuffle. --- --- @grabble xs m n@ is /O(m*n')/, where @n' = min n (length xs)@. Choose @n@ --- elements from @xs@, without replacement, and that @m@ times. -grabble :: MonadRandom m => [a] -> Int -> Int -> m [[a]] -grabble xs m n = do - swapss <- replicateM m $ forM [0 .. min (maxIx - 1) n] $ \i -> do - j <- getRandomR (i, maxIx) - return (i, j) - return $ map (take n . swapElems xs) swapss - where - maxIx = length xs - 1 - --- grabbleOnce :: MonadRandom m => [a] -> Int -> m [a] --- grabbleOnce xs n = head `liftM` grabble xs 1 n - -swapElems :: [a] -> [(Int, Int)] -> [a] -swapElems xs swaps = elems $ runSTArray (do - arr <- newListArray (0, maxIx) xs - mapM_ (swap arr) swaps - return arr) - where - maxIx = length xs - 1 - swap arr (i,j) = do - vi <- readArray arr i - vj <- readArray arr j - writeArray arr i vj - writeArray arr j vi diff --git a/src/Simulate/Options.hs b/src/Simulate/Options.hs deleted file mode 100644 --- a/src/Simulate/Options.hs +++ /dev/null @@ -1,162 +0,0 @@ -{- | -Module : Simulate.Options -Description : Argument parser for seq-ana -Copyright : (c) Dominik Schrempf 2019 -License : GPL-3 - -Maintainer : dominik.schrempf@gmail.com -Stability : unstable -Portability : portable - -Creation date: Fri May 3 11:51:07 2019. - --} - -module Simulate.Options - ( SimulateArguments (..) - , Simulate - , simulateArguments - , reportSimulateArguments - , simulateFooter - ) where - -import Control.Monad.Logger -import Control.Monad.Trans.Reader -import Data.List -import Data.Word -import Options.Applicative - -import ELynx.Tools.Options - --- | Arguments need to simulate phylogenetic trees. -data SimulateArguments = SimulateArguments - { argsNTrees :: Int -- ^ Simulated trees. - , argsNLeaves :: Int -- ^ Number of leaves. - , argsHeight :: Maybe Double -- ^ Tree height (time to origin or MRCA). - , argsConditionMRCA :: Bool -- ^ False: condition on origin; True: - -- condition on MRCA. - , argsLambda :: Double -- ^ Birth rate. - , argsMu :: Double -- ^ Death rate. - , argsRho :: Double -- ^ Smapling rate. - , argsSubSample :: Bool -- ^ Perform actual sub-sampling. - , argsSumStat :: Bool -- ^ Only print summary statistics? - , argsSeed :: Maybe [Word32] -- ^ Seed of NRG, random if 'Nothing'. - } - --- | Logger and reader data type. -type Simulate = LoggingT (ReaderT SimulateArguments IO) - --- | Print useful information about the provided arguments. -reportSimulateArguments :: SimulateArguments -> String -reportSimulateArguments a = - intercalate "\n" [ "Number of simulated trees: " ++ show (argsNTrees a) - , "Number of leaves per tree: " ++ show (argsNLeaves a) - , "Height of trees: " ++ hStr - , "Birth rate: " ++ show (argsLambda a) - , "Death rate: " ++ show (argsMu a) - , "Sampling probability: " ++ show (argsRho a) - , "Perform sub-sampling: " ++ show (argsSubSample a) - , "Summary statistics only: " ++ show (argsSumStat a) - , "Seed: " ++ sStr ] - where hStr = case argsHeight a of Nothing -> "Random height of origin" - Just h -> show h ++ ", conditioned on " - ++ if argsConditionMRCA a then "MRCA" else "origin" - sStr = case argsSeed a of Nothing -> "Random" - Just i -> show i - --- | Command line parser. -simulateArguments :: Parser SimulateArguments -simulateArguments = SimulateArguments - <$> nTreeOpt - <*> nLeavesOpt - <*> treeHeightOpt - <*> conditionMRCAOpt - <*> lambdaOpt - <*> muOpt - <*> rhoOpt - <*> subSampleOpt - <*> sumStatOpt - <*> seedOpt - -nTreeOpt :: Parser Int -nTreeOpt = option auto $ - long "nTrees" - <> short 't' - <> metavar "INT" - <> value 10 - <> showDefault - <> help "Number of trees" - -nLeavesOpt :: Parser Int -nLeavesOpt = option auto $ - long "nLeaves" - <> short 'n' - <> metavar "INT" - <> value 5 - <> showDefault - <> help "Number of leaves per tree" - -treeHeightOpt :: Parser (Maybe Double) -treeHeightOpt = optional $ option auto $ - long "height" - <> short 'H' - <> metavar "DOUBLE" - <> help "Fix tree height (no default)" - -conditionMRCAOpt :: Parser Bool -conditionMRCAOpt = switch $ - long "condition-on-mrca" - <> short 'M' - <> showDefault - <> help "Do not condition on height of origin but on height of MRCA" - -lambdaOpt :: Parser Double -lambdaOpt = option auto $ - long "lambda" - <> short 'l' - <> metavar "DOUBLE" - <> value 1.0 - <> showDefault - <> help "Birth rate lambda" - -muOpt :: Parser Double -muOpt = option auto $ - long "mu" - <> short 'm' - <> metavar "DOUBLE" - <> value 0.9 - <> showDefault - <> help "Death rate mu" - -rhoOpt :: Parser Double -rhoOpt = option auto $ - long "rho" - <> short 'r' - <> metavar "DOUBLE" - <> value 1.0 - <> help "Sampling probability rho (default: 1.0)" - -subSampleOpt :: Parser Bool -subSampleOpt = switch $ - long "sub-sample" - <> short 'u' - <> showDefault - <> help "Perform sub-sampling; see below." - -sumStatOpt :: Parser Bool -sumStatOpt = switch $ - long "summary-statistics" - <> short 's' - <> showDefault - <> help "Only output number of children for each branch" - -citation :: String -citation = "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" - --- | And a footer. -simulateFooter :: String -simulateFooter = intercalate "\n" - [ "Height of Trees: if no tree height is given, the heights will be randomly drawn from the expected distribution given the number of leaves, the birth and the death rate." - , "Summary statistics only: only print (NumberOfExtantChildren BranchLength) pairs for each branch of each tree. The trees are separated by a newline character." - , "Sub-sampling: simulate one big tree with n'=round(n/rho), n'>=n, leaves, and randomly sample sub-trees with n leaves. Hence, with rho=1.0, the same tree is reported over and over again." - , citation ] diff --git a/src/Simulate/Simulate.hs b/src/Simulate/Simulate.hs deleted file mode 100644 --- a/src/Simulate/Simulate.hs +++ /dev/null @@ -1,113 +0,0 @@ -{-# LANGUAGE OverloadedStrings #-} -{-# LANGUAGE TemplateHaskell #-} -{-# LANGUAGE TupleSections #-} - -{- | - Description : Simulate reconstructed trees - Copyright : (c) Dominik Schrempf 2018 - License : GPL-3 - - Maintainer : dominik.schrempf@gmail.com - Stability : unstable - Portability : portable - -Creation date: Tue Feb 27 17:27:16 2018. - -See Gernhard, T. (2008). The conditioned reconstructed process. Journal of -Theoretical Biology, 253(4), 769–778. http://doi.org/10.1016/j.jtbi.2008.04.005. - -Mon Feb 4 14:26:11 CET 2019: Adding sampling probability rho. See Article -(Stadler2009) Stadler, T. On incomplete sampling under birth–death models and -connections to the sampling-based coalescent Journal of Theoretical Biology, -Elsevier BV, 2009, 261, 58-66 - --} - -module Simulate.Simulate - ( simulate - ) -where - -import Control.Concurrent (getNumCapabilities) -import Control.Concurrent.Async.Lifted.Safe (mapConcurrently) -import Control.Monad -import Control.Monad.IO.Class -import Control.Monad.Logger -import Control.Monad.Trans.Class -import Control.Monad.Trans.Reader -import Control.Parallel.Strategies -import qualified Data.ByteString.Lazy.Char8 as L -import Data.Maybe -import qualified Data.Sequence as Seq -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 Simulate.Options - -import ELynx.Data.Tree.MeasurableTree -import ELynx.Data.Tree.PhyloTree (PhyloLabel) -import ELynx.Data.Tree.SumStat (formatNChildSumStat, - toNChildSumStat) -import ELynx.Data.Tree.Tree -import ELynx.Export.Tree.Newick (toNewick) -import ELynx.Simulate.PointProcess (simulateNReconstructedTrees, - simulateReconstructedTree) -import ELynx.Tools.Concurrent -import ELynx.Tools.InputOutput -import ELynx.Tools.Logger - --- | Simulate phylogenetic trees. -simulate :: Maybe FilePath -> Simulate () -simulate outFile = do - a <- lift ask - when (isNothing (argsHeight a) && argsConditionMRCA a) $ - error "Cannot condition on MRCA (-M) when height is not given (-H)." - let s = argsSumStat a - nCap <- liftIO getNumCapabilities - logNewSection "Arguments" - $(logInfo) $ T.pack $ reportSimulateArguments a - logNewSection "Simulation" - $(logInfo) $ T.pack $ "Number of used cores: " <> show nCap - trs <- if argsSubSample a - then simulateAndSubSampleNTreesConcurrently nCap - else simulateNTreesConcurrently nCap - let ls = if s - then parMap rpar (formatNChildSumStat . toNChildSumStat) trs - else parMap rpar toNewick trs - let outFile' = (++ ".tree") <$> outFile - let res = L.unlines ls - out "simulated trees" res outFile' - -simulateNTreesConcurrently :: Int -> Simulate [Tree (PhyloLabel Int)] -simulateNTreesConcurrently c = do - (SimulateArguments nT nL h cM l m r _ _ s) <- lift ask - let l' = l * r - m' = m - l * (1.0 - r) - gs <- liftIO $ getNGen c s - let chunks = getChunks c nT - timeSpec = fmap (, cM) h - trss <- liftIO $ mapConcurrently - (\(n, g) -> simulateNReconstructedTrees n nL timeSpec l' m' g) - (zip chunks gs) - return $ concat trss - -simulateAndSubSampleNTreesConcurrently :: Int -> Simulate [Tree (PhyloLabel Int)] -simulateAndSubSampleNTreesConcurrently c = do - (SimulateArguments nT nL h cM l m r _ _ s) <- lift ask - let nLeavesBigTree = (round $ fromIntegral nL / r) :: Int - gs <- liftIO $ getNGen c s - let chunks = getChunks c nT - timeSpec = fmap (, cM) h - tr <- liftIO $ simulateReconstructedTree nLeavesBigTree timeSpec l m (head gs) - logNewSection $ T.pack $ "Simulate one big tree with " <> show nLeavesBigTree <> " leaves." - -- Log the base tree. - $(logInfo) $ LT.toStrict $ LT.decodeUtf8 $ toNewick tr - logNewSection $ T.pack $ "Sub sample " <> show nT <> " trees with " <> show nL <> " leaves." - let lvs = Seq.fromList $ leaves tr - trss <- liftIO $ mapConcurrently - (\(nSamples, g) -> nSubSamples nSamples lvs nL tr g) - (zip chunks gs) - let trs = catMaybes $ concat trss - return $ map prune trs diff --git a/src/TLynx/Coalesce/Coalesce.hs b/src/TLynx/Coalesce/Coalesce.hs new file mode 100644 --- /dev/null +++ b/src/TLynx/Coalesce/Coalesce.hs @@ -0,0 +1,111 @@ +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE TupleSections #-} + +{- | + Description : Simulate reconstructed trees using the coalescent process + Copyright : (c) Dominik Schrempf 2020 + License : GPL-3.0-or-later + + Maintainer : dominik.schrempf@gmail.com + Stability : unstable + Portability : portable + +Creation date: Tue Feb 3 17:00:00 2020. + +-} + +module TLynx.Coalesce.Coalesce + ( coalesce + ) +where + +import Control.Concurrent ( getNumCapabilities ) +import Control.Concurrent.Async.Lifted.Safe + ( mapConcurrently ) +import Control.Monad +import Control.Monad.IO.Class +import Control.Monad.Trans.Reader ( ask ) +import Control.Monad.Logger +import Control.Parallel.Strategies +import qualified Data.ByteString.Lazy.Char8 as L +import Data.Maybe +import qualified Data.Sequence as Seq +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 System.Random.MWC ( initialize ) + +import TLynx.Coalesce.Options + +import ELynx.Data.Tree +import ELynx.Export.Tree.Newick ( toNewick ) +import ELynx.Simulate.Coalescent ( simulate ) + +import ELynx.Tools + +-- | Simulate phylogenetic trees. +coalesce :: ELynx CoalesceArguments () +coalesce = do + l <- local <$> ask + let s = argsSumStat l + logNewSection "Arguments" + $(logInfo) $ T.pack $ reportCoalesceArguments l + logNewSection "Simulation" + c <- liftIO getNumCapabilities + $(logInfo) $ T.pack $ "Number of used cores: " <> show c + trs <- case argsRho l of + Nothing -> simulateNTreesConcurrently + Just _ -> simulateAndSubSampleNTreesConcurrently + let ls = if s + then parMap rpar (formatNChildSumStat . toNChildSumStat) trs + else parMap rpar toNewick trs + let res = L.unlines ls + out "simulated trees" res ".tree" + +simulateNTreesConcurrently :: ELynx CoalesceArguments [Tree (PhyloLabel Int)] +simulateNTreesConcurrently = do + (CoalesceArguments nT nL _ _ (Fixed s)) <- local <$> ask + c <- liftIO getNumCapabilities + gs <- liftIO $ initialize s >>= \g -> splitGen c g + let chunks = getChunks c nT + trss <- liftIO $ mapConcurrently (\(n, g) -> replicateM n $ simulate nL g) + (zip chunks gs) + return $ concat trss + +simulateAndSubSampleNTreesConcurrently + :: ELynx CoalesceArguments [Tree (PhyloLabel Int)] +simulateAndSubSampleNTreesConcurrently = do + (CoalesceArguments nT nL mR _ (Fixed s)) <- local <$> ask + c <- liftIO getNumCapabilities + let + r = fromMaybe + (error + "cimulateAndSubSampleNTreesConcurrently: no sampling probability given." + ) + mR + let nLeavesBigTree = (round $ fromIntegral nL / r) :: Int + gs <- liftIO $ initialize s >>= \g -> splitGen c g + let chunks = getChunks c nT + tr <- liftIO $ simulate nLeavesBigTree (head gs) + logNewSection + $ T.pack + $ "Simulate one big tree with " + <> show nLeavesBigTree + <> " leaves." + -- Log the base tree. + $(logInfo) $ LT.toStrict $ LT.decodeUtf8 $ toNewick tr + logNewSection + $ T.pack + $ "Sub sample " + <> show nT + <> " trees with " + <> show nL + <> " leaves." + let lvs = Seq.fromList $ leaves tr + trss <- liftIO $ mapConcurrently + (\(nSamples, g) -> nSubSamples nSamples lvs nL tr g) + (zip chunks gs) + let trs = catMaybes $ concat trss + return $ map prune trs diff --git a/src/TLynx/Coalesce/Options.hs b/src/TLynx/Coalesce/Options.hs new file mode 100644 --- /dev/null +++ b/src/TLynx/Coalesce/Options.hs @@ -0,0 +1,121 @@ +{-# LANGUAGE DeriveGeneric #-} + +{- | +Module : TLynx.Coalesce.Options +Description : Argument parser for @slynx coalesce@ +Copyright : (c) Dominik Schrempf 2020 +License : GPL-3.0-or-later + +Maintainer : dominik.schrempf@gmail.com +Stability : unstable +Portability : portable + +Creation date: Fri May 3 11:51:07 2019. + +-} + +module TLynx.Coalesce.Options + ( CoalesceArguments(..) + , coalesceArguments + , reportCoalesceArguments + , coalesceFooter + ) +where + +import Data.List +import Options.Applicative + +import ELynx.Tools + +-- | Arguments need to simulate phylogenetic trees using the coalescent process. +data CoalesceArguments = CoalesceArguments + { argsNTrees :: Int -- ^ Number of simulated trees. + , argsNLeaves :: Int -- ^ Number of leaves per tree. + , argsRho :: Maybe Double -- ^ Perform random sub-sampling with given sampling rate. + , argsSumStat :: Bool -- ^ Only print summary statistics? + , argsSeed :: Seed -- ^ Seed of NRG, random if 'Nothing'. + } + deriving (Eq, Show, Generic) + +instance Reproducible CoalesceArguments where + inFiles _ = [] + outSuffixes _ = [".tree"] + getSeed = Just . argsSeed + setSeed a s = a { argsSeed = Fixed s } + parser = coalesceArguments + cmdName = "coalesce" + cmdDsc = + [ "Simulate phylogenetic trees using the coalescent processes (see also the 'simulate' command for simulations using the birth and death process)." + ] + cmdFtr = coalesceFooter + +instance FromJSON CoalesceArguments + +instance ToJSON CoalesceArguments + +-- | Print useful information about the provided arguments. +reportCoalesceArguments :: CoalesceArguments -> String +reportCoalesceArguments a = intercalate + "\n" + [ "Number of simulated trees: " ++ show (argsNTrees a) + , "Number of leaves per tree: " ++ show (argsNLeaves a) + , "Sub-sampling: " ++ ssStr + , "Summary statistics only: " ++ show (argsSumStat a) + ] + where + ssStr = case argsRho a of + Nothing -> "No" + Just r -> "Yes, with probability " <> show r + +-- | Command line parser. +coalesceArguments :: Parser CoalesceArguments +coalesceArguments = + CoalesceArguments + <$> nTreeOpt + <*> nLeavesOpt + <*> rhoOpt + <*> sumStatOpt + <*> seedOpt + +nTreeOpt :: Parser Int +nTreeOpt = + option auto + $ long "nTrees" + <> short 't' + <> metavar "INT" + <> value 10 + <> showDefault + <> help "Number of trees" + +nLeavesOpt :: Parser Int +nLeavesOpt = + option auto + $ long "nLeaves" + <> short 'n' + <> metavar "INT" + <> value 5 + <> showDefault + <> help "Number of leaves per tree" + +rhoOpt :: Parser (Maybe Double) +rhoOpt = + optional + $ option auto + $ long "sub-sample" + <> short 'r' + <> metavar "DOUBLE" + <> showDefault + <> help + "Perform random sub-sampling with sampling probability rho (see below)" + +sumStatOpt :: Parser Bool +sumStatOpt = + switch $ long "summary-statistics" <> short 's' <> showDefault <> help + "Only output number of children for each branch" + +-- | And a footer. +coalesceFooter :: [String] +coalesceFooter = + [ "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." + , "Summary statistics: only print (NumberOfExtantChildren BranchLength) pairs for each branch of each tree. The trees are separated by a newline character." + ] diff --git a/src/TLynx/Compare/Compare.hs b/src/TLynx/Compare/Compare.hs new file mode 100644 --- /dev/null +++ b/src/TLynx/Compare/Compare.hs @@ -0,0 +1,274 @@ +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE TemplateHaskell #-} + +{- | +Module : TLynx.Compare.Compare +Description : Compare two phylogenies +Copyright : (c) Dominik Schrempf 2020 +License : GPL-3.0-or-later + +Maintainer : dominik.schrempf@gmail.com +Stability : unstable +Portability : portable + +Creation date: Thu Sep 19 15:01:52 2019. + +-} + +module TLynx.Compare.Compare + ( compareCmd + ) +where + +import Control.Monad +import Control.Monad.IO.Class +import Control.Monad.Logger +import Control.Monad.Trans.Reader ( ask ) +import qualified Data.ByteString.Lazy.Char8 as L +import Data.List ( intercalate ) +import qualified Data.Map as M +import Data.Monoid +import qualified Data.Set as S +import qualified Data.Text as T +import qualified Data.Text.Encoding as E +import qualified Data.Text.IO as T +import Graphics.Gnuplot.Simple + +import Data.Tree +import System.IO +import Text.Printf + +import TLynx.Compare.Options + +import ELynx.Data.Tree +import ELynx.Export.Tree.Newick ( toNewick ) +import ELynx.Import.Tree.Newick +import ELynx.Tools ( parseFileWith + , outHandle + , ELynx + , Arguments(..) + , GlobalArguments(..) + ) + +treesOneFile + :: FilePath + -> ELynx + CompareArguments + (Tree (PhyloLabel L.ByteString), Tree (PhyloLabel L.ByteString)) +treesOneFile tf = do + nwF <- argsNewickFormat . local <$> ask + $(logInfo) $ T.pack $ "Parse file '" ++ tf ++ "'." + ts <- liftIO $ parseFileWith (manyNewick nwF) tf + let n = length ts + case compare n 2 of + LT -> error "Not enough trees in file." + GT -> error "Too many trees in file." + EQ -> return (head ts, head . tail $ ts) + +treesTwoFiles + :: FilePath + -> FilePath + -> ELynx + CompareArguments + (Tree (PhyloLabel L.ByteString), Tree (PhyloLabel L.ByteString)) +treesTwoFiles tf1 tf2 = do + nwF <- argsNewickFormat . local <$> ask + $(logInfo) $ T.pack $ "Parse first tree file '" ++ tf1 ++ "'." + t1 <- liftIO $ parseFileWith (oneNewick nwF) tf1 + $(logInfo) $ T.pack $ "Parse second tree file '" ++ tf2 ++ "'." + t2 <- liftIO $ parseFileWith (oneNewick nwF) tf2 + return (t1, t2) + +-- | More detailed comparison of two trees. +compareCmd :: ELynx CompareArguments () +compareCmd = do + l <- local <$> ask + -- Determine output handle (stdout or file). + outH <- outHandle "results" ".out" + + -- Read input. + let inFiles = argsInFiles l + nFiles = length inFiles + (tr1, tr2) <- case nFiles of + 1 -> treesOneFile (head inFiles) + 2 -> treesTwoFiles (head inFiles) (head . tail $ inFiles) + _ -> + error + "Need two input files with one tree each or one input file with two trees." + + liftIO $ hPutStrLn outH "Tree 1:" + liftIO $ L.hPutStrLn outH $ toNewick tr1 + liftIO $ hPutStrLn outH "Tree 2:" + liftIO $ L.hPutStrLn outH $ toNewick tr2 + liftIO $ hPutStrLn outH "" + + -- Intersect trees. + (t1, t2) <- if argsIntersect l + then do + let [x, y] = intersectWith getName extendBranchLength [tr1, tr2] + liftIO $ hPutStrLn outH "Intersected trees are:" + liftIO $ L.hPutStrLn outH $ toNewick x + liftIO $ L.hPutStrLn outH $ toNewick y + return (x, y) + else return (tr1, tr2) + + -- Check input (moved to library functions). + -- let lvs1 = leaves t1 + -- lvs2 = leaves t2 + -- lfns1 = map getName lvs1 + -- lfns2 = map getName lvs2 + -- s1 = S.fromList lfns1 + -- s2 = S.fromList lfns2 + -- if s1 == s2 + -- then liftIO $ hPutStrLn outH "Trees have the same set of leaf names." + -- else error "Trees do not have the same set of leaf names." + -- liftIO $ hPutStrLn outH "" + + -- Distances. + let formatD str val = T.justifyLeft 25 ' ' str <> " " <> val + liftIO $ hPutStrLn outH "Distances." + liftIO $ T.hPutStrLn outH $ formatD "Symmetric" + (T.pack $ show $ symmetric t1 t2) + liftIO $ T.hPutStrLn outH $ formatD + "Branch score" + (T.pack $ show $ branchScoreWith getName getLen t1 t2) + let t1' = normalizeBranchSupport t1 + t2' = normalizeBranchSupport t2 + $(logDebug) "Trees with normalized branch support values:" + $(logDebug) $ E.decodeUtf8 $ L.toStrict $ toNewick t1' + $(logDebug) $ E.decodeUtf8 $ L.toStrict $ toNewick t2' + liftIO $ T.hPutStrLn outH $ formatD + "Incompatible split" + (T.pack $ show $ incompatibleSplits t1' t2') + liftIO $ T.hPutStrLn outH $ formatD + "Incompatible split (0.10)" + (T.pack $ show $ incompatibleSplits (collapse 0.1 t1') (collapse 0.1 t2')) + liftIO $ T.hPutStrLn outH $ formatD + "Incompatible split (0.50)" + (T.pack $ show $ incompatibleSplits (collapse 0.5 t1') (collapse 0.5 t2')) + -- liftIO $ T.hPutStrLn outH $ formatD "Incompatible split (0.60)" + -- (T.pack $ show $ incompatibleSplits (collapse 0.6 t1') (collapse 0.6 t2')) + -- liftIO $ T.hPutStrLn outH $ formatD "Incompatible split (0.70)" + -- (T.pack $ show $ incompatibleSplits (collapse 0.7 t1') (collapse 0.7 t2')) + liftIO $ T.hPutStrLn outH $ formatD + "Incompatible split (0.80)" + (T.pack $ show $ incompatibleSplits (collapse 0.8 t1') (collapse 0.8 t2')) + liftIO $ T.hPutStrLn outH $ formatD + "Incompatible split (0.90)" + (T.pack $ show $ incompatibleSplits (collapse 0.9 t1') (collapse 0.9 t2')) + -- liftIO $ T.hPutStrLn outH $ formatD "Incompatible split (1.01)" + -- (T.pack $ show $ incompatibleSplits (collapse 1.01 t1') (collapse 1.01 t2')) + -- liftIO $ L.hPutStrLn outH $ toNewick (collapse 1.01 t1') + + -- Bipartitions. + when + (argsBipartitions l) + (do + let bp1 = bipartitions (fmap getName t1) + bp2 = bipartitions (fmap getName t2) + bp1Only = bp1 S.\\ bp2 + bp2Only = bp2 S.\\ bp1 + unless + (S.null bp1Only) + (do + liftIO $ hPutStrLn outH "" + liftIO + $ hPutStrLn outH "Bipartitions in Tree 1 that are not in Tree 2." + -- let bp1Strs = map (bphuman L.unpack . bpmap getName) (S.toList bp1Only) + forM_ bp1Only (liftIO . hPutStrLn outH . bphuman L.unpack) + ) + -- let bp1Strs = map (bphuman L.unpack) (S.toList bp1Only) + -- liftIO $ hPutStrLn outH $ intercalate "\n" bp1Strs) + unless + (S.null bp2Only) + (do + liftIO $ hPutStrLn outH "" + liftIO + $ hPutStrLn outH "Bipartitions in Tree 2 that are not in Tree 1." + forM_ bp2Only (liftIO . hPutStrLn outH . bphuman L.unpack) + ) + + -- Common bipartitions and their respective differences in branch lengths. + liftIO $ hPutStrLn outH "" + let bpCommon = bp1 `S.intersection` bp2 + if S.null bpCommon + then do + liftIO $ hPutStrLn outH "There are no common bipartitions." + liftIO $ hPutStrLn outH "No plots have been generated." + else do + let + bpToBrLen1 = + M.map getSum $ bipartitionToBranchLength getName (Sum . getLen) t1 + bpToBrLen2 = + M.map getSum $ bipartitionToBranchLength getName (Sum . getLen) t2 + liftIO $ hPutStrLn + outH + "Common bipartitions and their respective differences in branch lengths." + -- Header. + liftIO $ hPutStrLn outH header + forM_ + bpCommon + ( liftIO + . hPutStrLn outH + . getCommonBpStr L.unpack bpToBrLen1 bpToBrLen2 + ) + + -- XXX: This circumvents the extension checking, and hash creation for + -- elynx files. + bn <- outFileBaseName . global <$> ask + case bn of + Nothing -> + $(logInfo) "No output file name provided. Do not generate plots." + Just fn -> do + let compareCommonBps = + [ (bpToBrLen1 M.! b, bpToBrLen2 M.! b) + | b <- S.toList bpCommon + ] + liftIO $ epspdfPlot fn (plotBps compareCommonBps) + $(logInfo) + "Comparison of branch lengths plot generated (EPS and PDF)" + ) + + liftIO $ hClose outH + +header :: String +header = intercalate " " $ cols ++ ["Bipartition"] + where + cols = map (T.unpack . T.justifyRight 12 ' ') + ["Length 1", "Length 2", "Delta", "Relative [%]"] + +getCommonBpStr + :: (Ord a, Fractional b, PrintfArg b) + => (a -> String) + -> M.Map (Bipartition a) b + -> M.Map (Bipartition a) b + -> Bipartition a + -> String +getCommonBpStr f m1 m2 p = intercalate + " " + [ printf "% 12.7f" l1 + , printf "% 12.7f" l2 + , printf "% 12.7f" d + , printf "% 12.7f" rd + , s + ] + where + l1 = m1 M.! p + l2 = m2 M.! p + d = l1 - l2 + rd = 2 * d / (l1 + l2) + s = bphuman f p + +plotBps :: [(Double, Double)] -> [Attribute] -> IO () +plotBps xs as = plotPathsStyle as' [(ps1, xs), (ps2, line)] + where + as' = + as + ++ [ Title "Comparison of branch lengths of common branches" + , XLabel "Branch lengths, tree 1" + , YLabel "Branch lengths, tree 2" + ] + ps1 = PlotStyle Points (DefaultStyle 1) + m = maximum $ map fst xs ++ map snd xs + line = [(0, 0), (m, m)] + ps2 = PlotStyle Lines (DefaultStyle 1) diff --git a/src/TLynx/Compare/Options.hs b/src/TLynx/Compare/Options.hs new file mode 100644 --- /dev/null +++ b/src/TLynx/Compare/Options.hs @@ -0,0 +1,83 @@ +{-# LANGUAGE DeriveGeneric #-} + +{- | +Module : TLynx.Compare.Options +Description : Options for the compare subcommand +Copyright : (c) Dominik Schrempf 2020 +License : GPL-3.0-or-later + +Maintainer : dominik.schrempf@gmail.com +Stability : unstable +Portability : portable + +Creation date: Thu Sep 19 15:02:21 2019. + +-} + +module TLynx.Compare.Options + ( CompareArguments(..) + , compareArguments + ) +where + +import Options.Applicative + +import ELynx.Tools + +import TLynx.Parsers + +-- | Arguments of compare command. +data CompareArguments = CompareArguments + { argsNormalize :: Bool + , argsBipartitions :: Bool + , argsIntersect :: Bool + , argsNewickFormat :: NewickFormat + , argsInFiles :: [FilePath] } + deriving (Eq, Show, Generic) + +instance Reproducible CompareArguments where + inFiles = argsInFiles +-- XXX: The plots are not checked with the ELynx file, because they are not +-- always created and this is hard with the actual setup. + outSuffixes _ = [".out"] +-- outSuffixes a = ".out" : if argsBipartitions a then [".eps", ".pdf"] else [] + getSeed _ = Nothing + setSeed a _ = a + parser = compareArguments + cmdName = "compare" + cmdDsc = + [ "Compare two phylogenetic trees (compute distances and branch-wise differences)." + ] + +instance FromJSON CompareArguments + +instance ToJSON CompareArguments + +-- | Parse arguments of compare command. +compareArguments :: Parser CompareArguments +compareArguments = + CompareArguments + <$> normalize + <*> bipartitions + <*> intersect + <*> newickFormat + <*> some file + +normalize :: Parser Bool +normalize = switch $ long "normalize" <> short 'n' <> help + "Normalize trees before comparison" + +bipartitions :: Parser Bool +bipartitions = switch $ long "bipartitions" <> short 'b' <> help + "Print and plot common and missing bipartitions" + +intersect :: Parser Bool +intersect = + switch + $ long "intersect" + <> short 't' + <> help + "Compare intersections; i.e., before comparison, drop leaves that are not present in the other tree" + +file :: Parser FilePath +file = strArgument $ metavar "NAME" <> help "Tree file" diff --git a/src/TLynx/Connect/Connect.hs b/src/TLynx/Connect/Connect.hs new file mode 100644 --- /dev/null +++ b/src/TLynx/Connect/Connect.hs @@ -0,0 +1,116 @@ +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE TemplateHaskell #-} + +{- | +Module : TLynx.Connect.Connect +Description : Connect two phylogenies +Copyright : (c) Dominik Schrempf 2020 +License : GPL-3.0-or-later + +Maintainer : dominik.schrempf@gmail.com +Stability : unstable +Portability : portable + +Creation date: Thu Sep 19 15:01:52 2019. + +-} + +module TLynx.Connect.Connect + ( connectCmd + ) +where + +import Control.Monad.IO.Class +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 Data.Tree +import System.IO + +import TLynx.Connect.Options + +import ELynx.Data.Tree +import ELynx.Export.Tree.Newick ( toNewick ) +import ELynx.Import.Tree.Newick ( manyNewick + , oneNewick + ) +import ELynx.Tools ( outHandle + , parseFileWith + , ELynx + , Arguments(..) + , fromBs + , tShow + ) + +-- | Connect two trees honoring possible constraints. See 'connect'. +connectCmd :: ELynx ConnectArguments () +connectCmd = do + lArgs <- local <$> ask + outH <- outHandle "results" ".out" + + -- Do we have constraints or not? + let cs = constraints lArgs + l = inFileA lArgs + r = inFileB lArgs + case cs of + Nothing -> connectOnly outH l r + Just c -> connectAndFilter outH c l r + + liftIO $ hClose outH + +connectTrees + :: Tree (PhyloLabel L.ByteString) + -> Tree (PhyloLabel L.ByteString) + -> [Tree (PhyloLabel L.ByteString)] +connectTrees = connect (PhyloLabel "" Nothing Nothing) + +type Constraint a = S.Set a + +compatibleAll :: (Show a, Ord a) => Tree a -> [Constraint a] -> Bool +compatibleAll (Node _ [l, r]) cs = + all (bpcompatible (bipartition l)) cs && all (bpcompatible (bipartition r)) cs +compatibleAll _ _ = error "Tree is not bifurcating." + +compatibleWith + :: (Show b, Ord b) => (a -> b) -> [Constraint a] -> Tree a -> Bool +compatibleWith f cs t = compatibleAll (fmap f t) (map (S.map f) cs) + +parseTrees + :: FilePath + -> FilePath + -> ELynx + ConnectArguments + (Tree (PhyloLabel L.ByteString), Tree (PhyloLabel L.ByteString)) +parseTrees l r = do + nwF <- nwFormat . local <$> ask + tl <- liftIO $ parseFileWith (oneNewick nwF) l + tr <- liftIO $ parseFileWith (oneNewick nwF) r + $(logInfo) "Tree 1:" + $(logInfo) $ fromBs $ toNewick tl + $(logInfo) "Tree 2:" + $(logInfo) $ fromBs $ toNewick tr + return (tl, tr) + +connectOnly :: Handle -> FilePath -> FilePath -> ELynx ConnectArguments () +connectOnly h l r = do + (tl, tr) <- parseTrees l r + let ts = connectTrees tl tr + $(logInfo) $ "Connected trees: " <> tShow (length ts) + liftIO $ L.hPutStr h $ L.unlines $ map toNewick ts + +connectAndFilter + :: Handle -> FilePath -> FilePath -> FilePath -> ELynx ConnectArguments () +connectAndFilter h c l r = do + nwF <- nwFormat . local <$> ask + cts <- liftIO $ parseFileWith (manyNewick nwF) c + $(logInfo) "Constraints:" + $(logInfo) $ fromBs $ L.intercalate "\n" $ map toNewick cts + (tl, tr) <- parseTrees l r + let ts = connectTrees tl tr + cs = concatMap clades cts :: [Constraint (PhyloLabel L.ByteString)] + -- Only collect trees that are compatible with the constraints. + ts' = filter (compatibleWith getName cs) ts + $(logInfo) $ "Connected trees: " <> tShow (length ts) + $(logInfo) $ "Compatible trees: " <> tShow (length ts') + liftIO $ L.hPutStr h $ L.unlines $ map toNewick ts' diff --git a/src/TLynx/Connect/Options.hs b/src/TLynx/Connect/Options.hs new file mode 100644 --- /dev/null +++ b/src/TLynx/Connect/Options.hs @@ -0,0 +1,72 @@ +{-# LANGUAGE DeriveGeneric #-} + +{- | +Module : TLynx.Connect.Options +Description : Options for the connect subcommand +Copyright : (c) Dominik Schrempf 2020 +License : GPL-3.0-or-later + +Maintainer : dominik.schrempf@gmail.com +Stability : unstable +Portability : portable + +Creation date: Thu Sep 19 15:02:21 2019. + +-} + +module TLynx.Connect.Options + ( ConnectArguments(..) + , connectArguments + ) +where + +import Options.Applicative + +import ELynx.Tools + +import TLynx.Parsers + +-- | Arguments of connect command. +data ConnectArguments = ConnectArguments + { nwFormat :: NewickFormat + , constraints :: Maybe FilePath + , inFileA :: FilePath + , inFileB :: FilePath } + deriving (Eq, Show, Generic) + +instance Reproducible ConnectArguments where + inFiles a = [inFileA a, inFileB a] + outSuffixes _ = [".out"] + getSeed _ = Nothing + setSeed a _ = a + parser = connectArguments + cmdName = "connect" + cmdDsc = + [ "Connect two phylogenetic trees in all ways (possibly honoring constraints)." + ] + +instance FromJSON ConnectArguments + +instance ToJSON ConnectArguments + +-- | Parse arguments of connect command. +connectArguments :: Parser ConnectArguments +connectArguments = + ConnectArguments <$> newickFormat <*> constraintsFile <*> fileA <*> fileB + +constraintsFile :: Parser (Maybe FilePath) +constraintsFile = + optional + $ strOption + $ metavar "CONSTRAINTS" + <> short 'c' + <> long "contraints" + <> help "File containing one or more Newick trees to be used as constraints" + +fileA :: Parser FilePath +fileA = strArgument $ metavar "TREE-FILE-A" <> help + "File containing the first Newick tree" + +fileB :: Parser FilePath +fileB = strArgument $ metavar "TREE-FILE-B" <> help + "File containing the second Newick tree" diff --git a/src/TLynx/Distance/Distance.hs b/src/TLynx/Distance/Distance.hs new file mode 100644 --- /dev/null +++ b/src/TLynx/Distance/Distance.hs @@ -0,0 +1,212 @@ +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE TemplateHaskell #-} + +{- | +Description : Compute distances between trees +Copyright : (c) Dominik Schrempf 2020 +License : GPL-3.0-or-later + +Maintainer : dominik.schrempf@gmail.com +Stability : unstable +Portability : portable + +Creation date: Wed May 29 18:09:39 2019. + +- Symmetric (Robinson-Foulds) distance. +- Incompatible splits distance. + +-} + +module TLynx.Distance.Distance + ( distance + ) +where + +import Control.Monad ( unless + , when + ) +import Control.Monad.IO.Class ( liftIO ) +import Control.Monad.Logger ( logDebug + , logInfo + ) +import Control.Monad.Trans.Class ( lift ) +import Control.Monad.Trans.Reader ( ask ) +import qualified Data.ByteString.Lazy.Char8 as L +import Data.List ( sort ) +import Data.Maybe ( isNothing ) +import qualified Data.Text as T +import qualified Data.Text.IO as T +import qualified Data.Text.Lazy as LT +import qualified Data.Text.Lazy.Encoding as LT +import Data.Tree ( Tree ) +import qualified Data.Vector.Unboxed as V +import Statistics.Sample ( mean + , variance + ) +import System.IO ( hClose + , hPutStrLn + ) +import Text.Printf ( PrintfArg + , printf + ) + +import TLynx.Distance.Options + +import ELynx.Data.Tree +import ELynx.Export.Tree.Newick +import ELynx.Import.Tree.Newick +import ELynx.Tools ( alignLeft + , alignRight + , tShow + , outHandle + , parseFileWith + , ELynx + , Arguments(..) + ) + +median :: Ord a => [a] -> a +median xs = sort xs !! l2 where l2 = length xs `div` 2 + +pf :: String +pf = "%.3f" + +header :: Int -> Int -> DistanceMeasure -> L.ByteString +header n m d = + alignLeft (n + 2) "Tree 1" <> alignLeft (n + 2) "Tree 2" <> alignRight + (m + 2) + (L.pack $ show d) + +showTriplet + :: (PrintfArg a) => Int -> Int -> [String] -> (Int, Int, a) -> L.ByteString +showTriplet n m args (i, j, d) = i' <> j' <> d' + where + i' = alignLeft (n + 2) $ L.pack (args !! i) + j' = alignLeft (n + 2) $ L.pack (args !! j) + d' = alignRight (m + 2) $ L.pack (printf pf d) + +-- | Compute distance functions between phylogenetic trees. +distance :: ELynx DistanceArguments () +distance = do + l <- local <$> ask + let nwFormat = argsNewickFormat l + -- Determine output handle (stdout or file). + outH <- outHandle "results" ".out" + -- Master tree (in case it is given). + let mname = argsMasterTreeFile l + mtree <- case mname of + Nothing -> return Nothing + Just f -> do + $(logInfo) $ T.pack $ "Read master tree from file: " <> f <> "." + t <- liftIO $ parseFileWith (oneNewick nwFormat) f + $(logInfo) "Compute distances between all trees and master tree." + return $ Just t + let tfps = argsInFiles l + (trees, names) <- case tfps of + [] -> error "No tree input files given." + [tf] -> do + $(logInfo) "Read trees from single file." + ts <- liftIO $ parseFileWith (manyNewick nwFormat) tf + $(logInfo) $ tShow (length ts) <> " trees found in file." + $(logInfo) "Trees are indexed with integers." + return (ts, map show [0 .. length ts - 1]) + _ -> do + $(logInfo) "Read trees from files." + ts <- liftIO $ mapM (parseFileWith (oneNewick nwFormat)) tfps + $(logInfo) "Trees are named according to their file names." + return (ts, tfps) + + when (null trees) (error "Not enough trees found in files.") + when (isNothing mtree && length trees == 1) + (error "Not enough trees found in files.") + -- when (isNothing mtree) $ $(logInfo) + -- "Compute pairwise distances between trees from different files." + $(logDebug) "The trees are:" + $(logDebug) $ LT.toStrict $ LT.decodeUtf8 $ L.unlines $ map toNewick trees + + -- Set the distance measure. + let dist = argsDistance l + case argsDistance l of + Symmetric -> $(logInfo) "Use symmetric (Robinson-Foulds) distance." + IncompatibleSplit val -> do + $(logInfo) "Use incompatible split distance." + $(logInfo) + $ T.pack + $ "Collapse nodes with support less than " + ++ show val + ++ "." + BranchScore -> $(logInfo) "Use branch score distance." + let distanceMeasure' + :: Tree (PhyloLabel L.ByteString) + -> Tree (PhyloLabel L.ByteString) + -> Double + distanceMeasure' = case dist of + Symmetric -> \t1 t2 -> fromIntegral $ symmetric t1 t2 + IncompatibleSplit _ -> + \t1 t2 -> fromIntegral $ incompatibleSplits t1 t2 + BranchScore -> branchScore + -- Possibly intersect trees before distance calculation. + when (argsIntersect l) + $ $(logInfo) "Intersect trees before calculation of distances." + let distanceMeasure = if argsIntersect l + then + (\t1 t2 -> + let [t1', t2'] = intersectWith getName extendBranchLength [t1, t2] + in distanceMeasure' t1' t2' + ) + else distanceMeasure' + + -- Possibly normalize trees. + when (argsNormalize l) + $ $(logInfo) "Normalize trees before calculation of distances." + let normalizeF = if argsNormalize l then normalizeBranchLength else id + + -- Possibly collapse unsupported nodes. + let collapseF = case dist of + -- For the incompatible split distance we have to collapse branches with + -- support lower than the given value. Before doing so, we normalize the + -- branch support values. + IncompatibleSplit val -> collapse val . normalizeBranchSupport + _ -> id + + -- The trees can be prepared now. + let trees' = map (collapseF . normalizeF) trees + + $(logDebug) "The prepared trees are:" + $(logDebug) $ LT.toStrict $ LT.decodeUtf8 $ L.unlines $ map toNewick trees' + let dsTriplets = case mtree of + Nothing -> pairwise distanceMeasure trees' + Just t -> + [ (0, i, distanceMeasure t t') | (i, t') <- zip [1 ..] trees' ] + ds = map (\(_, _, x) -> x) dsTriplets + dsVec = V.fromList ds + + liftIO + $ hPutStrLn outH + $ "Summary statistics of " + ++ show dist + ++ " Distance:" + liftIO $ T.hPutStrLn outH $ T.justifyLeft 10 ' ' "Mean: " <> T.pack + (printf pf (mean dsVec)) + liftIO $ T.hPutStrLn outH $ T.justifyLeft 10 ' ' "Median: " <> T.pack + (printf pf (median ds)) + liftIO $ T.hPutStrLn outH $ T.justifyLeft 10 ' ' "Variance: " <> T.pack + (printf pf (variance dsVec)) + -- L.putStrLn $ L.unlines $ map toNewick ts + -- L.putStrLn $ L.unlines $ map toNewick tsN + -- L.putStrLn $ L.unlines $ map toNewick tsC + + lift $ unless + (argsSummaryStatistics l) + (do + let n = maximum $ 6 : map length names + m = length $ show dist + lift $ hPutStrLn outH "" + lift $ L.hPutStrLn outH $ header n m dist + case mname of + Nothing -> lift $ L.hPutStr outH $ L.unlines + (map (showTriplet n m names) dsTriplets) + Just mn -> lift $ L.hPutStr outH $ L.unlines + (map (showTriplet n m (mn : names)) dsTriplets) + ) + + liftIO $ hClose outH diff --git a/src/TLynx/Distance/Options.hs b/src/TLynx/Distance/Options.hs new file mode 100644 --- /dev/null +++ b/src/TLynx/Distance/Options.hs @@ -0,0 +1,184 @@ +{-# LANGUAGE DeriveGeneric #-} + +{- | +Module : TLynx.Distance.Options +Description : Options of tree-dist +Copyright : (c) Dominik Schrempf 2020 +License : GPL-3.0-or-later + +Maintainer : dominik.schrempf@gmail.com +Stability : unstable +Portability : portable + +Creation date: Thu Aug 29 13:02:22 2019. + +-} + +module TLynx.Distance.Options + ( DistanceArguments(..) + , DistanceMeasure(..) + , distanceArguments + , distanceFooter + ) +where + +import Data.Void +import Data.Scientific ( toRealFloat ) +import Options.Applicative +import Text.Megaparsec ( Parsec + , eof + , try + ) +import Text.Megaparsec.Char ( char + , string + ) +import Text.Megaparsec.Char.Lexer ( scientific ) +import Text.Printf + +import ELynx.Tools + +import TLynx.Parsers + +-- | Supported distance measures. +data DistanceMeasure = + Symmetric -- ^ Symmetric distance. + | IncompatibleSplit Double -- ^ Incompatible split distance; collapse nodes + -- with branch support below given value. + | BranchScore -- ^ Branch score distance. + deriving (Eq, Generic) + +instance FromJSON DistanceMeasure + +instance ToJSON DistanceMeasure + +instance Show DistanceMeasure where + show Symmetric = "Symmetric" + show (IncompatibleSplit c) = "Incompatible Split (" ++ printf "%.2f" c ++ ")" + show BranchScore = "Branch Score" + +-- | Arguments needed to compute distance measures. +data DistanceArguments = DistanceArguments + { argsDistance :: DistanceMeasure + , argsNormalize :: Bool + , argsIntersect :: Bool + , argsSummaryStatistics :: Bool + , argsMasterTreeFile :: Maybe FilePath + , argsNewickFormat :: NewickFormat + , argsInFiles :: [FilePath] + } + deriving (Eq, Show, Generic) + +instance Reproducible DistanceArguments where + inFiles a = case argsMasterTreeFile a of + Nothing -> argsInFiles a + Just f -> f : argsInFiles a + outSuffixes _ = [".out"] + getSeed _ = Nothing + setSeed = const + parser = distanceArguments + cmdName = "distance" + cmdDsc = ["Compute distances between many phylogenetic trees."] + cmdFtr = distanceFooter + +instance FromJSON DistanceArguments + +instance ToJSON DistanceArguments + +-- | COmmand line parser. +distanceArguments :: Parser DistanceArguments +distanceArguments = + DistanceArguments + <$> distanceOpt + <*> normalizeSwitch + <*> intersectSwitch + <*> summaryStatisticsSwitch + <*> masterTreeFile + <*> newickFormat + <*> many inFilesArg + +masterTreeFile :: Parser (Maybe FilePath) +masterTreeFile = + optional + $ strOption + $ long "master-tree-file" + <> short 'm' + <> metavar "MASTER-TREE-File" + <> help "Compare all trees to the tree in the master tree file." + +inFilesArg :: Parser FilePath +inFilesArg = + strArgument + $ metavar "INPUT-FILES" + <> help + "Read tree(s) from INPUT-FILES; if more files are given, one tree is expected per file" + +symmetric :: Parsec Void String DistanceMeasure +symmetric = do + _ <- string "symmetric" + _ <- eof + pure Symmetric + +incompatibleSplit :: Parsec Void String DistanceMeasure +incompatibleSplit = do + _ <- string "incompatible-split" + _ <- char '[' + f <- toRealFloat <$> scientific + _ <- char ']' + _ <- eof + if (0 <= f) && (f <= 1) + then pure $ IncompatibleSplit f + else error "Branch support has to be in [0, 1]." + +branchScore :: Parsec Void String DistanceMeasure +branchScore = do + _ <- string "branch-score" + _ <- eof + pure BranchScore + +distanceParser :: Parsec Void String DistanceMeasure +distanceParser = + -- Try first the normalized one, since the normal branch score + -- parser also succeeds in this case. + try symmetric <|> try incompatibleSplit <|> branchScore + +distanceOpt :: Parser DistanceMeasure +distanceOpt = + option (megaReadM distanceParser) + $ long "distance" + <> short 'd' + <> metavar "MEASURE" + <> help + "Type of distance to calculate (available distance measures are listed below)" + +summaryStatisticsSwitch :: Parser Bool +summaryStatisticsSwitch = + switch $ long "summary-statistics" <> short 's' <> help + "Report summary statistics only" + +normalizeSwitch :: Parser Bool +normalizeSwitch = + switch + $ long "normalize" + <> short 'n' + <> help + "Normalize trees before distance calculation; only affect distances depending on branch lengths" + +intersectSwitch :: Parser Bool +intersectSwitch = + switch + $ long "intersect" + <> short 't' + <> help + "Compare intersections; i.e., before comparison, drop leaves that are not present in the other tree" + +-- | Information about provided distance types. +distanceFooter :: [String] +distanceFooter = + [ "Distance measures:" + , " symmetric Symmetric distance (Robinson-Foulds distance)." + , " incompatible-split[VAL] Incompatible split distance. Collapse branches" + , " with support less than VAL before distance calculation;" + , " in this way, only well supported difference contribute" + , " to the distance measure." + , " branch-score Branch score distance." + ] diff --git a/src/TLynx/Examine/Examine.hs b/src/TLynx/Examine/Examine.hs new file mode 100644 --- /dev/null +++ b/src/TLynx/Examine/Examine.hs @@ -0,0 +1,65 @@ +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE TemplateHaskell #-} + +{- | +Description : Analyze trees +Copyright : (c) Dominik Schrempf 2020 +License : GPL-3.0-or-later + +Maintainer : dominik.schrempf@gmail.com +Stability : unstable +Portability : portable + +Creation date: Fri May 24 13:47:56 2019. + +-} + +module TLynx.Examine.Examine + ( examine + ) +where + +import Control.Monad ( unless ) +import Control.Monad.IO.Class +import Control.Monad.Logger +import Control.Monad.Trans.Reader ( ask ) +import qualified Data.ByteString.Lazy.Char8 as L +import Data.List ( nub + , (\\) + ) +import qualified Data.Text as T +import Data.Tree +import System.IO ( Handle + , hPutStrLn + ) + +import TLynx.Examine.Options + +import ELynx.Data.Tree +import ELynx.Import.Tree.Newick +import ELynx.Tools + +readTrees :: FilePath -> ELynx ExamineArguments [Tree (PhyloLabel L.ByteString)] +readTrees fp = do + $(logInfo) $ T.pack $ "Read tree(s) from file " <> fp <> "." + nf <- argsNewickFormat . local <$> ask + liftIO $ parseFileWith (manyNewick nf) fp + +examineTree :: (Measurable a, Named a) => Handle -> Tree a -> IO () +examineTree h t = do + L.hPutStrLn h $ summarize t + unless (null dups) $ hPutStrLn h "" >> hPutStrLn + h + ("Duplicate leaves: " ++ show dups) + where + lvs = map getName $ leaves t + dups = lvs \\ nub lvs + +-- | Examine phylogenetic trees. +examine :: ELynx ExamineArguments () +examine = do + l <- local <$> ask + let inFn = argsInFile l + trs <- readTrees inFn + outH <- outHandle "results" ".out" + liftIO $ mapM_ (examineTree outH) trs diff --git a/src/TLynx/Examine/Options.hs b/src/TLynx/Examine/Options.hs new file mode 100644 --- /dev/null +++ b/src/TLynx/Examine/Options.hs @@ -0,0 +1,54 @@ +{-# LANGUAGE DeriveGeneric #-} + +{- | +Module : TLynx.Examine.Options +Description : Tree analysis options +Copyright : (c) Dominik Schrempf 2020 +License : GPL-3.0-or-later + +Maintainer : dominik.schrempf@gmail.com +Stability : unstable +Portability : portable + +Creation date: Thu Aug 29 08:16:45 2019. + +-} + +module TLynx.Examine.Options + ( ExamineArguments(..) + , examineArguments + ) +where + +import Options.Applicative + +import ELynx.Tools + +import TLynx.Parsers + +-- | Arguments needed to examine phylogenetic trees. +data ExamineArguments = ExamineArguments + { argsInFile :: FilePath + , argsNewickFormat :: NewickFormat } + deriving (Eq, Show, Generic) + +instance Reproducible ExamineArguments where + inFiles = pure . argsInFile + outSuffixes _ = [".out"] + getSeed _ = Nothing + setSeed = const + parser = examineArguments + cmdName = "examine" + cmdDsc = ["Compute summary statistics of phylogenetic trees."] + +instance FromJSON ExamineArguments + +instance ToJSON ExamineArguments + +-- | Command line parser. +examineArguments :: Parser ExamineArguments +examineArguments = ExamineArguments <$> inFile <*> newickFormat + +inFile :: Parser FilePath +inFile = + strArgument $ metavar "INPUT-FILE" <> help "Read trees from INPUT-FILE" diff --git a/src/TLynx/Options.hs b/src/TLynx/Options.hs new file mode 100644 --- /dev/null +++ b/src/TLynx/Options.hs @@ -0,0 +1,124 @@ +{-# LANGUAGE DeriveGeneric #-} + +{- | +Module : TLynx.Options +Description : TLynx 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 TLynx.Options + ( Arguments(..) + , CommandArguments(..) + , parseArguments + ) +where + +import Options.Applicative + +import ELynx.Tools + +import TLynx.Coalesce.Options +import TLynx.Compare.Options +import TLynx.Connect.Options +import TLynx.Distance.Options +import TLynx.Examine.Options +import TLynx.Shuffle.Options +import TLynx.Simulate.Options +import TLynx.Parsers ( newickHelp ) + +-- | The different TLynx commands and their arguments. +data CommandArguments = + Coalesce CoalesceArguments + | Compare CompareArguments + | Connect ConnectArguments + | Distance DistanceArguments + | Examine ExamineArguments + | Shuffle ShuffleArguments + | Simulate SimulateArguments + deriving (Eq, Show, Generic) + +instance Reproducible CommandArguments where + inFiles (Coalesce a) = inFiles a + inFiles (Compare a) = inFiles a + inFiles (Connect a) = inFiles a + inFiles (Distance a) = inFiles a + inFiles (Examine a) = inFiles a + inFiles (Shuffle a) = inFiles a + inFiles (Simulate a) = inFiles a + + outSuffixes (Coalesce a) = outSuffixes a + outSuffixes (Compare a) = outSuffixes a + outSuffixes (Connect a) = outSuffixes a + outSuffixes (Distance a) = outSuffixes a + outSuffixes (Examine a) = outSuffixes a + outSuffixes (Shuffle a) = outSuffixes a + outSuffixes (Simulate a) = outSuffixes a + + getSeed (Coalesce a) = getSeed a + getSeed (Compare a) = getSeed a + getSeed (Connect a) = getSeed a + getSeed (Distance a) = getSeed a + getSeed (Examine a) = getSeed a + getSeed (Shuffle a) = getSeed a + getSeed (Simulate a) = getSeed a + + setSeed (Coalesce a) = Coalesce . setSeed a + setSeed (Compare a) = Compare . setSeed a + setSeed (Connect a) = Connect . setSeed a + setSeed (Distance a) = Distance . setSeed a + setSeed (Examine a) = Examine . setSeed a + setSeed (Shuffle a) = Shuffle . setSeed a + setSeed (Simulate a) = Simulate . setSeed a + + parser = commandArguments + + cmdName = "tlynx" + + cmdDsc = ["Compare, examine, and simulate phylogenetic trees."] + + cmdFtr = "Available tree file formats:" : indent newickHelp + where indent = map (" " ++) + +instance FromJSON CommandArguments + +instance ToJSON CommandArguments + +coalesceCommand :: Mod CommandFields CommandArguments +coalesceCommand = createCommandReproducible Coalesce + +compareCommand :: Mod CommandFields CommandArguments +compareCommand = createCommandReproducible Compare + +connectCommand :: Mod CommandFields CommandArguments +connectCommand = createCommandReproducible Connect + +distanceCommand :: Mod CommandFields CommandArguments +distanceCommand = createCommandReproducible Distance + +examineCommand :: Mod CommandFields CommandArguments +examineCommand = createCommandReproducible Examine + +shuffleCommand :: Mod CommandFields CommandArguments +shuffleCommand = createCommandReproducible Shuffle + +simulateCommand :: Mod CommandFields CommandArguments +simulateCommand = createCommandReproducible Simulate + +commandArguments :: Parser CommandArguments +commandArguments = + hsubparser + $ coalesceCommand + <> compareCommand + <> connectCommand + <> distanceCommand + <> examineCommand + <> shuffleCommand + <> simulateCommand diff --git a/src/TLynx/Parsers.hs b/src/TLynx/Parsers.hs new file mode 100644 --- /dev/null +++ b/src/TLynx/Parsers.hs @@ -0,0 +1,43 @@ +{- | +Module : TLynx.Parsers +Description : Common parsers +Copyright : (c) Dominik Schrempf 2020 +License : GPL-3.0-or-later + +Maintainer : dominik.schrempf@gmail.com +Stability : unstable +Portability : portable + +Creation date: Wed Apr 22 13:34:35 2020. + +-} + +module TLynx.Parsers + ( NewickFormat + , newickFormat + , newickHelp + ) +where + +import Options.Applicative + +import ELynx.Import.Tree.Newick ( NewickFormat(..) + , description + ) +import ELynx.Tools ( allValues ) + +-- | Parse 'NewickFormat'. +newickFormat :: Parser NewickFormat +newickFormat = + option auto + $ long "newick-format" + <> short 'f' + <> metavar "FORMAT" + <> value Standard + <> showDefault + <> help "Newick tree format; see 'tlynx --help'" + +-- | Help for different 'NewickFormat's. +newickHelp :: [String] +newickHelp = map (toListItem . description) (allValues :: [NewickFormat]) + where toListItem = ("- Newick " ++) diff --git a/src/TLynx/Shuffle/Options.hs b/src/TLynx/Shuffle/Options.hs new file mode 100644 --- /dev/null +++ b/src/TLynx/Shuffle/Options.hs @@ -0,0 +1,63 @@ +{-# LANGUAGE DeriveGeneric #-} + +{- | +Module : TLynx.Shuffle.Options +Description : Options for the connect subcommand +Copyright : (c) Dominik Schrempf 2020 +License : GPL-3.0-or-later + +Maintainer : dominik.schrempf@gmail.com +Stability : unstable +Portability : portable + +Creation date: Thu Sep 19 15:02:21 2019. + +-} + +module TLynx.Shuffle.Options + ( ShuffleArguments(..) + , shuffleArguments + ) +where + +import Options.Applicative + +import ELynx.Tools + +import TLynx.Parsers + +-- | Arguments of shuffle command. +data ShuffleArguments = ShuffleArguments + { nwFormat :: NewickFormat + , nReplicates :: Int + , inFile :: FilePath + , argsSeed :: Seed } + deriving (Eq, Show, Generic) + +instance Reproducible ShuffleArguments where + inFiles = pure . inFile + outSuffixes _ = [".tree"] + getSeed = Just . argsSeed + setSeed a s = a { argsSeed = Fixed s } + parser = shuffleArguments + cmdName = "shuffle" + cmdDsc = + [ "Shuffle a phylogenetic tree (keep coalescent times, but shuffle topology and leaves)." + ] + +instance FromJSON ShuffleArguments + +instance ToJSON ShuffleArguments + +-- | Parse arguments of shuffle command. +shuffleArguments :: Parser ShuffleArguments +shuffleArguments = ShuffleArguments <$> newickFormat <*> n <*> file <*> seedOpt + +n :: Parser Int +n = + option auto $ long "replicates" <> short 'n' <> metavar "N" <> value 1 <> help + "Number of trees to generate" + +file :: Parser FilePath +file = + strArgument $ metavar "TREE-FILE" <> help "File containing a Newick tree" diff --git a/src/TLynx/Shuffle/Shuffle.hs b/src/TLynx/Shuffle/Shuffle.hs new file mode 100644 --- /dev/null +++ b/src/TLynx/Shuffle/Shuffle.hs @@ -0,0 +1,160 @@ +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE TemplateHaskell #-} + +{- | +Module : TLynx.Shuffle.Shuffle +Description : Shuffle a phylogeny +Copyright : (c) Dominik Schrempf 2020 +License : GPL-3.0-or-later + +Maintainer : dominik.schrempf@gmail.com +Stability : unstable +Portability : portable + +Creation date: Thu Sep 19 15:01:52 2019. + +The coalescent times are unaffected. The topology and the leaf order is +shuffled. Branch support values are ignored and lost. + +-} + +module TLynx.Shuffle.Shuffle + ( shuffleCmd + ) +where + +import Control.Comonad ( extend ) +import Control.Monad ( forM + , replicateM + , when + ) +import Control.Monad.IO.Class ( liftIO ) +import Control.Monad.Logger ( logDebug + , logInfo + ) +import Control.Monad.Trans.Reader ( ask ) +import Control.Monad.Primitive ( PrimMonad + , PrimState + ) +import Data.Array ( elems ) +import Data.Array.ST ( newListArray + , readArray + , runSTArray + , writeArray + ) +import qualified Data.ByteString.Lazy.Char8 as L +import Data.Maybe ( isNothing ) +import Data.Tree ( Tree + , flatten + , rootLabel + ) +import System.IO ( hClose ) +import System.Random.MWC ( Gen + , initialize + , uniformR + ) + +import TLynx.Shuffle.Options + +import ELynx.Data.Tree +import ELynx.Export.Tree.Newick ( toNewick ) +import ELynx.Import.Tree.Newick ( oneNewick ) +import ELynx.Simulate.PointProcess ( PointProcess(PointProcess) + , toReconstructedTree + ) +import ELynx.Tools + +-- | Shuffle a tree. Get all coalescent times, shuffle them. Get all leaves, +-- shuffle them. Connect the shuffled leaves with the shuffled coalescent times. +-- The shuffled tree has a new topology while keeping the same set of coalescent +-- times and leaves. +shuffleCmd :: ELynx ShuffleArguments () +shuffleCmd = do + l <- local <$> ask + h <- outHandle "results" ".tree" + let nwF = nwFormat l + t <- liftIO $ parseFileWith (oneNewick nwF) (inFile l) + $(logInfo) "Input tree:" + $(logInfo) $ fromBs $ toNewick t + + -- Check if all branches have a given length. However, the length of the stem is not important. + let r = rootLabel t + r' = r { brLen = Just 0 } + t' = t { rootLabel = r' } + when + (isNothing $ traverse brLen t') + (do + $(logDebug) $ tShow t' + error "Not all branches have a given length." + ) + + -- Check if tree is ultrametric enough. + let dh = sum $ map (height t -) (distancesOriginLeaves t) + $(logDebug) $ "Distance in branch length to being ultrametric: " <> tShow dh + when (dh > 2e-4) (error "Tree is not ultrametric.") + when (dh > eps && dh < 2e-4) + $ $(logInfo) + "Tree is nearly ultrametric, ignore branch length differences smaller than 2e-4." + when (dh < eps) $ $(logInfo) "Tree is ultrametric." + + let cs = filter (> 0) $ flatten $ extend rootHeight t + ls = map getName $ leaves t + $(logDebug) $ "Number of coalescent times: " <> tShow (length cs) + $(logDebug) $ "Number of leaves: " <> tShow (length ls) + $(logDebug) "The coalescent times are: " + $(logDebug) $ tShow cs + + gen <- case argsSeed l of + Random -> error "Seed not available; please contact maintainer." + Fixed s -> liftIO $ initialize s + + ts <- liftIO $ shuffle (nReplicates l) (height t) cs ls gen + liftIO $ L.hPutStr h $ L.unlines $ map toNewick ts + + liftIO $ hClose h + +shuffle + :: PrimMonad m + => Int -- How many? + -> Double -- Stem length. + -> [Double] -- Coalescent times. + -> [L.ByteString] -- Leave names. + -> Gen (PrimState m) + -> m [Tree (PhyloLabel L.ByteString)] +shuffle n o cs ls gen = do + css <- grabble cs n (length cs) gen + lss <- grabble ls n (length ls) gen + return + [ toReconstructedTree "" (PointProcess names times o) + | (times, names) <- zip css lss + ] + +-- | From https://wiki.haskell.org/Random_shuffle. +-- +-- @grabble xs m n@ is /O(m*n')/, where @n' = min n (length xs)@. Choose @n@ +-- elements from @xs@, without replacement, and that @m@ times. +grabble :: PrimMonad m => [a] -> Int -> Int -> Gen (PrimState m) -> m [[a]] +grabble xs m n gen = do + swapss <- replicateM m $ forM [0 .. min (maxIx - 1) n] $ \i -> do + j <- uniformR (i, maxIx) gen + return (i, j) + return $ map (take n . swapElems xs) swapss + where maxIx = length xs - 1 + +-- grabbleOnce :: MonadRandom m => [a] -> Int -> m [a] +-- grabbleOnce xs n = head `liftM` grabble xs 1 n + +swapElems :: [a] -> [(Int, Int)] -> [a] +swapElems xs swaps = elems $ runSTArray + (do + arr <- newListArray (0, maxIx) xs + mapM_ (swap arr) swaps + return arr + ) + where + maxIx = length xs - 1 + swap arr (i, j) = do + vi <- readArray arr i + vj <- readArray arr j + writeArray arr i vj + writeArray arr j vi diff --git a/src/TLynx/Simulate/Options.hs b/src/TLynx/Simulate/Options.hs new file mode 100644 --- /dev/null +++ b/src/TLynx/Simulate/Options.hs @@ -0,0 +1,176 @@ +{-# LANGUAGE DeriveGeneric #-} + +{- | +Module : TLynx.Simulate.Options +Description : Argument parser for seq-ana +Copyright : (c) Dominik Schrempf 2020 +License : GPL-3.0-or-later + +Maintainer : dominik.schrempf@gmail.com +Stability : unstable +Portability : portable + +Creation date: Fri May 3 11:51:07 2019. + +-} + +module TLynx.Simulate.Options + ( SimulateArguments(..) + , simulateArguments + , reportSimulateArguments + , simulateFooter + ) +where + +import Data.List +import Options.Applicative + +import ELynx.Tools + +-- | Arguments need to simulate phylogenetic trees using birth and death processes. +data SimulateArguments = SimulateArguments + { argsNTrees :: Int -- ^ Simulated trees. + , argsNLeaves :: Int -- ^ Number of leaves. + , argsHeight :: Maybe Double -- ^ Tree height (time to origin or MRCA). + , argsConditionMRCA :: Bool -- ^ False: condition on origin; True: + -- condition on MRCA. + , argsLambda :: Double -- ^ Birth rate. + , argsMu :: Double -- ^ Death rate. + , argsRho :: Double -- ^ Smapling rate. + , argsSubSample :: Bool -- ^ Perform actual sub-sampling. + , argsSumStat :: Bool -- ^ Only print summary statistics? + , argsSeed :: Seed -- ^ Seed of NRG, random if 'Nothing'. + } + deriving (Eq, Show, Generic) + +instance Reproducible SimulateArguments where + inFiles _ = [] + outSuffixes _ = [".tree"] + getSeed = Just . argsSeed + setSeed a s = a { argsSeed = Fixed s } + parser = simulateArguments + cmdName = "simulate" + cmdDsc = + [ "Simulate phylogenetic trees using birth and death processes (see also the 'coalesce' command for simulations using the coalescent process)." + ] + cmdFtr = simulateFooter + +instance FromJSON SimulateArguments + +instance ToJSON SimulateArguments + +-- | Print useful information about the provided arguments. +reportSimulateArguments :: SimulateArguments -> String +reportSimulateArguments a = intercalate + "\n" + [ "Number of simulated trees: " ++ show (argsNTrees a) + , "Number of leaves per tree: " ++ show (argsNLeaves a) + , "Height of trees: " ++ hStr + , "Birth rate: " ++ show (argsLambda a) + , "Death rate: " ++ show (argsMu a) + , "Sampling probability: " ++ show (argsRho a) + , "Perform sub-sampling: " ++ show (argsSubSample a) + , "Summary statistics only: " ++ show (argsSumStat a) + ] + where + hStr = case argsHeight a of + Nothing -> "Random height of origin" + Just h -> show h ++ ", conditioned on " ++ if argsConditionMRCA a + then "MRCA" + else "origin" + +-- | Command line parser. +simulateArguments :: Parser SimulateArguments +simulateArguments = + SimulateArguments + <$> nTreeOpt + <*> nLeavesOpt + <*> treeHeightOpt + <*> conditionMRCAOpt + <*> lambdaOpt + <*> muOpt + <*> rhoOpt + <*> subSampleOpt + <*> sumStatOpt + <*> seedOpt + +nTreeOpt :: Parser Int +nTreeOpt = + option auto + $ long "nTrees" + <> short 't' + <> metavar "INT" + <> value 10 + <> showDefault + <> help "Number of trees" + +nLeavesOpt :: Parser Int +nLeavesOpt = + option auto + $ long "nLeaves" + <> short 'n' + <> metavar "INT" + <> value 5 + <> showDefault + <> help "Number of leaves per tree" + +treeHeightOpt :: Parser (Maybe Double) +treeHeightOpt = + optional + $ option auto + $ long "height" + <> short 'H' + <> metavar "DOUBLE" + <> help "Fix tree height (no default)" + +conditionMRCAOpt :: Parser Bool +conditionMRCAOpt = + switch $ long "condition-on-mrca" <> short 'M' <> showDefault <> help + "Do not condition on height of origin but on height of MRCA" + +lambdaOpt :: Parser Double +lambdaOpt = + option auto + $ long "lambda" + <> short 'l' + <> metavar "DOUBLE" + <> value 1.0 + <> showDefault + <> help "Birth rate lambda" + +muOpt :: Parser Double +muOpt = + option auto + $ long "mu" + <> short 'm' + <> metavar "DOUBLE" + <> value 0.9 + <> showDefault + <> help "Death rate mu" + +rhoOpt :: Parser Double +rhoOpt = + option auto $ long "rho" <> short 'r' <> metavar "DOUBLE" <> value 1.0 <> help + "Sampling probability rho (default: 1.0)" + +subSampleOpt :: Parser Bool +subSampleOpt = switch $ long "sub-sample" <> short 'u' <> showDefault <> help + "Perform sub-sampling; see below." + +sumStatOpt :: Parser Bool +sumStatOpt = + switch $ long "summary-statistics" <> short 's' <> showDefault <> help + "Only output number of children for each branch" + +citation :: String +citation = + "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" + +-- | And a footer. +simulateFooter :: [String] +simulateFooter = + [ "Height of Trees: if no tree height is given, the heights will be randomly drawn from the expected distribution given the number of leaves, the birth and the death rate." + , "Summary statistics only: only print (NumberOfExtantChildren BranchLength) pairs for each branch of each tree. The trees are separated by a newline character." + , "Sub-sampling: simulate one big tree with n'=round(n/rho), n'>=n, leaves, and randomly sample sub-trees with n leaves. Hence, with rho=1.0, the same tree is reported over and over again." + , citation + ] diff --git a/src/TLynx/Simulate/Simulate.hs b/src/TLynx/Simulate/Simulate.hs new file mode 100644 --- /dev/null +++ b/src/TLynx/Simulate/Simulate.hs @@ -0,0 +1,140 @@ +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE TupleSections #-} + +{- | + Description : Simulate reconstructed trees + Copyright : (c) Dominik Schrempf 2018 + License : GPL-3.0-or-later + + Maintainer : dominik.schrempf@gmail.com + Stability : unstable + Portability : portable + +Creation date: Tue Feb 27 17:27:16 2018. + +See Gernhard, T. (2008). The conditioned reconstructed process. Journal of +Theoretical Biology, 253(4), 769–778. http://doi.org/10.1016/j.jtbi.2008.04.005. + +Mon Feb 4 14:26:11 CET 2019: Adding sampling probability rho. See Article +(Stadler2009) Stadler, T. On incomplete sampling under birth–death models and +connections to the sampling-based coalescent Journal of Theoretical Biology, +Elsevier BV, 2009, 261, 58-66 + +-} + +module TLynx.Simulate.Simulate + ( simulate + ) +where +import Control.Monad.Trans.Reader ( ask ) +import Control.Concurrent ( getNumCapabilities ) +import Control.Concurrent.Async.Lifted.Safe + ( mapConcurrently ) +import Control.Monad +import Control.Monad.IO.Class +import Control.Monad.Logger +import Control.Parallel.Strategies +import qualified Data.ByteString.Lazy.Char8 as L +import Data.Maybe +import qualified Data.Sequence as Seq +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 System.Random.MWC ( GenIO + , initialize + ) + +import TLynx.Simulate.Options + +import ELynx.Data.Tree +import ELynx.Export.Tree.Newick ( toNewick ) +import ELynx.Simulate.PointProcess ( TimeSpec + , simulateNReconstructedTrees + , simulateReconstructedTree + ) +import ELynx.Tools + +-- | Simulate phylogenetic trees. +simulate :: ELynx SimulateArguments () +simulate = do + l <- local <$> ask + let + SimulateArguments nTrees nLeaves tHeight mrca lambda mu rho subS sumS (Fixed s) + = l + -- error "simulate: seed not available; please contact maintainer." + when (isNothing tHeight && mrca) + $ error "Cannot condition on MRCA (-M) when height is not given (-H)." + c <- liftIO getNumCapabilities + logNewSection "Arguments" + $(logInfo) $ T.pack $ reportSimulateArguments l + logNewSection "Simulation" + $(logInfo) $ T.pack $ "Number of used cores: " <> show c + gs <- liftIO $ initialize s >>= \gen -> splitGen c gen + let chunks = getChunks c nTrees + timeSpec = fmap (, mrca) tHeight + trs <- if subS + then simulateAndSubSampleNTreesConcurrently nLeaves + lambda + mu + rho + timeSpec + chunks + gs + else simulateNTreesConcurrently nLeaves lambda mu rho timeSpec chunks gs + let ls = if sumS + then parMap rpar (formatNChildSumStat . toNChildSumStat) trs + else parMap rpar toNewick trs + let res = L.unlines ls + out "simulated trees" res ".tree" + +simulateNTreesConcurrently + :: Int + -> Double + -> Double + -> Double + -> TimeSpec + -> [Int] + -> [GenIO] + -> ELynx SimulateArguments [Tree (PhyloLabel Int)] +simulateNTreesConcurrently nLeaves l m r timeSpec chunks gs = do + let l' = l * r + m' = m - l * (1.0 - r) + trss <- liftIO $ mapConcurrently + (\(n, g) -> simulateNReconstructedTrees n nLeaves timeSpec l' m' g) + (zip chunks gs) + return $ concat trss + +simulateAndSubSampleNTreesConcurrently + :: Int + -> Double + -> Double + -> Double + -> TimeSpec + -> [Int] + -> [GenIO] + -> ELynx SimulateArguments [Tree (PhyloLabel Int)] +simulateAndSubSampleNTreesConcurrently nLeaves l m r timeSpec chunks gs = do + let nLeavesBigTree = (round $ fromIntegral nLeaves / r) :: Int + logNewSection + $ T.pack + $ "Simulate one big tree with " + <> show nLeavesBigTree + <> " leaves." + tr <- liftIO $ simulateReconstructedTree nLeavesBigTree timeSpec l m (head gs) + -- Log the base tree. + $(logInfo) $ LT.toStrict $ LT.decodeUtf8 $ toNewick tr + logNewSection + $ T.pack + $ "Sub sample " + <> show (sum chunks) + <> " trees with " + <> show nLeaves + <> " leaves." + let lvs = Seq.fromList $ leaves tr + trss <- liftIO $ mapConcurrently + (\(nSamples, g) -> nSubSamples nSamples lvs nLeaves tr g) + (zip chunks gs) + let trs = catMaybes $ concat trss + return $ map prune trs diff --git a/src/TLynx/TLynx.hs b/src/TLynx/TLynx.hs new file mode 100644 --- /dev/null +++ b/src/TLynx/TLynx.hs @@ -0,0 +1,54 @@ +{- | +Module : TLynx.TLynx +Description : TLynx 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:36:43 2020. + +-} + +module TLynx.TLynx + ( tlynx + , rTLynx + ) +where + +import TLynx.Options + +import TLynx.Coalesce.Coalesce +import TLynx.Compare.Compare +import TLynx.Connect.Connect +import TLynx.Distance.Distance +import TLynx.Examine.Examine +import TLynx.Shuffle.Shuffle +import TLynx.Simulate.Simulate + +import ELynx.Tools + +-- TODO: Use a class here (e.g., elynx-wrappable) which defines the extractor function. +-- | Run TLynx with given arguments. +tlynx :: Arguments CommandArguments -> IO () +tlynx c = case local c of + Coalesce _ -> + eLynxWrapper c (\(Arguments g (Coalesce l)) -> Arguments g l) coalesce + Compare _ -> + eLynxWrapper c (\(Arguments g (Compare l)) -> Arguments g l) compareCmd + Connect _ -> + eLynxWrapper c (\(Arguments g (Connect l)) -> Arguments g l) connectCmd + Distance _ -> + eLynxWrapper c (\(Arguments g (Distance l)) -> Arguments g l) distance + Examine _ -> + eLynxWrapper c (\(Arguments g (Examine l)) -> Arguments g l) examine + Shuffle _ -> + eLynxWrapper c (\(Arguments g (Shuffle l)) -> Arguments g l) shuffleCmd + Simulate _ -> + eLynxWrapper c (\(Arguments g (Simulate l)) -> Arguments g l) simulate + +-- | Run TLynx, parse arguments from command line. +rTLynx :: IO () +rTLynx = parseArguments >>= tlynx diff --git a/tlynx.cabal b/tlynx.cabal --- a/tlynx.cabal +++ b/tlynx.cabal @@ -1,9 +1,9 @@ cabal-version: 1.12 name: tlynx -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,76 +23,77 @@ library exposed-modules: - Compare.Compare - Compare.Options - Connect.Connect - Connect.Options - Distance.Distance - Distance.Options - Examine.Examine - Examine.Options - Shuffle.Options - Shuffle.Shuffle - Simulate.Options - Simulate.Simulate + TLynx.Coalesce.Coalesce + TLynx.Coalesce.Options + TLynx.Compare.Compare + TLynx.Compare.Options + TLynx.Connect.Connect + TLynx.Connect.Options + TLynx.Distance.Distance + TLynx.Distance.Options + TLynx.Examine.Examine + TLynx.Examine.Options + TLynx.Options + TLynx.Parsers + TLynx.Shuffle.Options + TLynx.Shuffle.Shuffle + TLynx.Simulate.Options + TLynx.Simulate.Simulate + TLynx.TLynx hs-source-dirs: src other-modules: Paths_tlynx default-language: Haskell2010 - ghc-options: -Wall -O2 + ghc-options: -Wall -fllvm build-depends: - MonadRandom >=0.5.1.1 && <0.6, - QuickCheck >=2.13.2 && <2.14, - array >=0.5.3.0 && <0.6, - base >=4.12.0.0 && <4.13, - bytestring >=0.10.8.2 && <0.11, - containers >=0.6.0.1 && <0.7, - elynx-seq >=0.1.0 && <0.2, - elynx-tools >=0.1.0 && <0.2, - elynx-tree >=0.1.0 && <0.2, - lifted-async >=0.10.0.4 && <0.11, - math-functions >=0.3.2.1 && <0.4, - megaparsec >=7.0.5 && <7.1, - monad-logger >=0.3.30 && <0.4, + aeson >=1.4.7.1 && <1.5, + array >=0.5.4.0 && <0.6, + base >=4.13.0.0 && <4.14, + bytestring >=0.10.10.0 && <0.11, + comonad >=5.0.6 && <5.1, + containers >=0.6.2.1 && <0.7, + elynx-tools >=0.2.1 && <0.3, + elynx-tree >=0.2.1 && <0.3, + gnuplot >=0.5.6.1 && <0.6, + lifted-async >=0.10.0.6 && <0.11, + 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, + optparse-applicative >=0.15.1.0 && <0.16, parallel >=3.2.2.0 && <3.3, - primitive >=0.6.4.0 && <0.7, - quickcheck-instances >=0.3.22 && <0.4, - statistics >=0.15.0.0 && <0.16, - text >=1.2.3.1 && <1.3, + primitive >=0.7.0.1 && <0.8, + scientific >=0.3.6.2 && <0.4, + statistics >=0.15.2.0 && <0.16, + text >=1.2.4.0 && <1.3, transformers >=0.5.6.2 && <0.6, - vector >=0.12.0.3 && <0.13 + vector >=0.12.1.2 && <0.13 executable tlynx main-is: Main.hs hs-source-dirs: app other-modules: - Options Paths_tlynx default-language: Haskell2010 - ghc-options: -Wall -O2 -eventlog -threaded -rtsopts - -with-rtsopts=-N + ghc-options: -Wall -fllvm -threaded -rtsopts -with-rtsopts=-N build-depends: - MonadRandom >=0.5.1.1 && <0.6, - QuickCheck >=2.13.2 && <2.14, - array >=0.5.3.0 && <0.6, - base >=4.12.0.0 && <4.13, - bytestring >=0.10.8.2 && <0.11, - containers >=0.6.0.1 && <0.7, - elynx-tools >=0.1.0 && <0.2, - elynx-tree >=0.1.0 && <0.2, - lifted-async >=0.10.0.4 && <0.11, - math-functions >=0.3.2.1 && <0.4, - megaparsec >=7.0.5 && <7.1, - monad-logger >=0.3.30 && <0.4, + aeson >=1.4.7.1 && <1.5, + array >=0.5.4.0 && <0.6, + base >=4.13.0.0 && <4.14, + bytestring >=0.10.10.0 && <0.11, + comonad >=5.0.6 && <5.1, + containers >=0.6.2.1 && <0.7, + elynx-tools >=0.2.1 && <0.3, + gnuplot >=0.5.6.1 && <0.6, + lifted-async >=0.10.0.6 && <0.11, + 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, + optparse-applicative >=0.15.1.0 && <0.16, parallel >=3.2.2.0 && <3.3, - primitive >=0.6.4.0 && <0.7, - quickcheck-instances >=0.3.22 && <0.4, - statistics >=0.15.0.0 && <0.16, - text >=1.2.3.1 && <1.3, + primitive >=0.7.0.1 && <0.8, + scientific >=0.3.6.2 && <0.4, + statistics >=0.15.2.0 && <0.16, + text >=1.2.4.0 && <1.3, tlynx -any, transformers >=0.5.6.2 && <0.6, - vector >=0.12.0.3 && <0.13 + vector >=0.12.1.2 && <0.13