packages feed

tlynx 0.1.0 → 0.2.1

raw patch · 34 files changed

+2495/−1632 lines, 34 filesdep +aesondep +comonaddep +gnuplotdep −MonadRandomdep −QuickCheckdep −elynx-seqdep ~arraydep ~basedep ~bytestring

Dependencies added: aeson, comonad, gnuplot, scientific

Dependencies removed: MonadRandom, QuickCheck, elynx-seq, math-functions, quickcheck-instances

Dependency ranges changed: array, base, bytestring, containers, elynx-tools, elynx-tree, lifted-async, megaparsec, monad-logger, optparse-applicative, primitive, statistics, text, vector

Files

ChangeLog.md view
@@ -1,10 +1,13 @@ -# Table of Contents--1.  [Changelog for ELynx](#org263dcac)+# Changelog for ELynx  -<a id="org263dcac"></a>+## 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. 
README.md view
@@ -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) 
app/Main.hs view
@@ -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
− app/Options.hs
@@ -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"]
− src/Compare/Compare.hs
@@ -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
− src/Compare/Options.hs
@@ -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)"
− src/Connect/Connect.hs
@@ -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'
− src/Connect/Options.hs
@@ -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"
− src/Distance/Distance.hs
@@ -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
− src/Distance/Options.hs
@@ -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."-  ]
− src/Examine/Examine.hs
@@ -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
− src/Examine/Options.hs
@@ -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)"
− src/Shuffle/Options.hs
@@ -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"
− src/Shuffle/Shuffle.hs
@@ -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
− src/Simulate/Options.hs
@@ -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 ]
− src/Simulate/Simulate.hs
@@ -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
+ src/TLynx/Coalesce/Coalesce.hs view
@@ -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
+ src/TLynx/Coalesce/Options.hs view
@@ -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."+  ]
+ src/TLynx/Compare/Compare.hs view
@@ -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)
+ src/TLynx/Compare/Options.hs view
@@ -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"
+ src/TLynx/Connect/Connect.hs view
@@ -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'
+ src/TLynx/Connect/Options.hs view
@@ -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"
+ src/TLynx/Distance/Distance.hs view
@@ -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
+ src/TLynx/Distance/Options.hs view
@@ -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."+  ]
+ src/TLynx/Examine/Examine.hs view
@@ -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
+ src/TLynx/Examine/Options.hs view
@@ -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"
+ src/TLynx/Options.hs view
@@ -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
+ src/TLynx/Parsers.hs view
@@ -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 " ++)
+ src/TLynx/Shuffle/Options.hs view
@@ -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"
+ src/TLynx/Shuffle/Shuffle.hs view
@@ -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
+ src/TLynx/Simulate/Options.hs view
@@ -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+  ]
+ src/TLynx/Simulate/Simulate.hs view
@@ -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
+ src/TLynx/TLynx.hs view
@@ -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
tlynx.cabal view
@@ -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