tlynx 0.2.2 → 0.3.0
raw patch · 21 files changed
+1517/−1630 lines, 21 filesdep ~aesondep ~arraydep ~base
Dependency ranges changed: aeson, array, base, bytestring, comonad, containers, elynx-tools, elynx-tree, gnuplot, lifted-async, megaparsec, monad-logger, mwc-random, optparse-applicative, parallel, primitive, scientific, statistics, text, transformers, vector
Files
- ChangeLog.md +11/−0
- README.md +68/−274
- app/Main.hs +12/−15
- src/TLynx/Coalesce/Coalesce.hs +105/−71
- src/TLynx/Coalesce/Options.hs +70/−66
- src/TLynx/Compare/Compare.hs +204/−197
- src/TLynx/Compare/Options.hs +46/−42
- src/TLynx/Connect/Connect.hs +95/−70
- src/TLynx/Connect/Options.hs +38/−36
- src/TLynx/Distance/Distance.hs +160/−137
- src/TLynx/Distance/Options.hs +89/−87
- src/TLynx/Examine/Examine.hs +45/−45
- src/TLynx/Examine/Options.hs +21/−25
- src/TLynx/Options.hs +55/−59
- src/TLynx/Parsers.hs +32/−31
- src/TLynx/Shuffle/Options.hs +27/−29
- src/TLynx/Shuffle/Shuffle.hs +52/−124
- src/TLynx/Simulate/Options.hs +116/−98
- src/TLynx/Simulate/Simulate.hs +166/−107
- src/TLynx/TLynx.hs +23/−27
- tlynx.cabal +82/−90
ChangeLog.md view
@@ -2,6 +2,17 @@ # Changelog for ELynx +## Version 0.3.1; in development+++## Version 0.3.0++- **`elynx-nexus`:** library to import and export Nexus files.+- **`elynx-tree`:** major refactor and big cleanup; use rose trees with branch+ labels.+- **`elynx-tree`:** provide zippers.++ ## Version 0.2.2 - Validation and repetition of previous analyses is finally possible with the
README.md view
@@ -2,9 +2,11 @@ # The ELynx Suite -Version: 0.2.2.+Version: 0.3.0. Reproducible evolution made easy. +<p align="center"><img src="https://travis-ci.org/dschrempf/elynx.svg?branch=master"/></p>+ 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@@ -12,11 +14,12 @@ sums so that validation of past analyses is possible without the need to recompute the result. -The Elynx Suite consists of four library packages and three executables providing-a range of sub commands.+The Elynx Suite consists of library packages and executables providing a range+of sub commands. The library packages are: +- **elynx-nexus:** Nexus file support. - **elynx-markov:** Simulate multi sequence alignments along phylogenetic trees. - **elynx-seq:** Handle evolutionary sequences and multi sequence alignments. - **elynx-tools:** Tools for the provided executables.@@ -43,7 +46,7 @@ 2. Clone the ELynx repository. - git clone clone https://github.com/dschrempf/elynx+ git clone https://github.com/dschrempf/elynx 3. Navigate to the newly created `elynx` folder and build the binaries. This will take a while.@@ -62,15 +65,34 @@ [PATH](https://en.wikipedia.org/wiki/PATH_(variable)) environment variable. Then, they can be used directly. +# Documentation++Documentation is available on [Hackage](https://hackage.haskell.org/).++Libraries:++- [elynx-nexus](https://hackage.haskell.org/package/elynx-nexus)+- [elynx-markov](https://hackage.haskell.org/package/elynx-markov)+- [elynx-seq](https://hackage.haskell.org/package/elynx-seq)+- [elynx-tools](https://hackage.haskell.org/package/elynx-tools)+- [elynx-tree](https://hackage.haskell.org/package/elynx-tree)++Executables:++- [elynx](https://hackage.haskell.org/package/elynx)+- [slynx](https://hackage.haskell.org/package/slynx)+- [tlynx](https://hackage.haskell.org/package/tlynx)++ # SLynx Handle evolutionary sequences. - slynx --help+ slynx --help | head -n -16 - ELynx Suite version 0.2.1.+ ELynx Suite version 0.3.0. Developed by Dominik Schrempf.- Compiled on April 27, 2020, at 12:08 pm, UTC.+ Compiled on July 30, 2020, at 12:40 pm, UTC. Usage: slynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME] [-f|--force] COMMAND@@ -95,6 +117,7 @@ sub-sample Sub-sample columns from multi sequence alignments. translate Translate from DNA to Protein or DNAX to ProteinX. + Available sequence file formats: - FASTA @@ -105,23 +128,6 @@ - Protein (amino acids) - ProteinX (amino acids; including gaps) - ProteinS (amino acids; including gaps, and translation stops)- - ProteinI (amino acids; including gaps, translation stops, and IUPAC codes)- - 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 ## Concatenate@@ -130,9 +136,9 @@ slynx concatenate --help - ELynx Suite version 0.2.1.+ ELynx Suite version 0.3.0. Developed by Dominik Schrempf.- Compiled on April 27, 2020, at 12:08 pm, UTC.+ Compiled on July 30, 2020, at 12:40 pm, UTC. Usage: slynx concatenate (-a|--alphabet NAME) INPUT-FILE Concatenate sequences found in input files.@@ -143,23 +149,6 @@ -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@@ -168,9 +157,9 @@ slynx examine --help - ELynx Suite version 0.2.1.+ ELynx Suite version 0.3.0. Developed by Dominik Schrempf.- Compiled on April 27, 2020, at 12:08 pm, UTC.+ Compiled on July 30, 2020, at 12:40 pm, UTC. Usage: slynx examine (-a|--alphabet NAME) INPUT-FILE [--per-site] Examine sequences. If data is a multi sequence alignment, additionally analyze columns.@@ -182,23 +171,6 @@ 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@@ -207,9 +179,9 @@ slynx filter-rows --help - ELynx Suite version 0.2.1.+ ELynx Suite version 0.3.0. Developed by Dominik Schrempf.- Compiled on April 27, 2020, at 12:08 pm, UTC.+ Compiled on July 30, 2020, at 12:40 pm, UTC. Usage: slynx filter-rows (-a|--alphabet NAME) INPUT-FILE [--longer-than LENGTH] [--shorter-than LENGTH] [--standard-characters]@@ -225,31 +197,14 @@ --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 - ELynx Suite version 0.2.1.+ ELynx Suite version 0.3.0. Developed by Dominik Schrempf.- Compiled on April 27, 2020, at 12:08 pm, UTC.+ Compiled on July 30, 2020, at 12:40 pm, UTC. Usage: slynx filter-columns (-a|--alphabet NAME) INPUT-FILE [--standard-chars DOUBLE]@@ -263,23 +218,6 @@ --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@@ -288,9 +226,9 @@ slynx simulate --help - ELynx Suite version 0.2.1.+ ELynx Suite version 0.3.0. Developed by Dominik Schrempf.- Compiled on April 27, 2020, at 12:08 pm, UTC.+ Compiled on July 30, 2020, at 12:40 pm, UTC. Usage: slynx simulate (-t|--tree-file Name) [-s|--substitution-model MODEL] [-m|--mixture-model MODEL] [-e|--edm-file NAME] @@ -352,27 +290,13 @@ where XX is 10, 20, 30, 40, 50, or 60; CXX models, Quang et al., 2008. -m "EDM(EXCHANGEABILITIES)" Arbitrary empirical distribution mixture (EDM) models.- Stationary distributions have to be provided with the -e option.+ Stationary distributions have to be provided with the -e or -p option. For example, LG exchangeabilities with stationary distributions given in FILE. -m "EDM(LG-Custom)" -e FILE+ LG exchangeabilities with site profiles (Phylobayes) given in FILES.+ -m "EDM(LG-Custom)" -p FILES 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@@ -381,9 +305,9 @@ slynx sub-sample --help - ELynx Suite version 0.2.1.+ ELynx Suite version 0.3.0. Developed by Dominik Schrempf.- Compiled on April 27, 2020, at 12:08 pm, UTC.+ Compiled on July 30, 2020, at 12:40 pm, UTC. Usage: slynx sub-sample (-a|--alphabet NAME) INPUT-FILE (-n|--number-of-sites INT)@@ -403,22 +327,6 @@ -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.- - 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@@ -427,9 +335,9 @@ slynx translate --help - ELynx Suite version 0.2.1.+ ELynx Suite version 0.3.0. Developed by Dominik Schrempf.- Compiled on April 27, 2020, at 12:08 pm, UTC.+ Compiled on July 30, 2020, at 12:40 pm, UTC. Usage: slynx translate (-a|--alphabet NAME) INPUT-FILE (-r|--reading-frame INT) (-u|--universal-code CODE)@@ -444,34 +352,17 @@ -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 Handle phylogenetic trees in Newick format. - tlynx --help+ tlynx --help | head -n -16 - ELynx Suite version 0.2.1.+ ELynx Suite version 0.3.0. Developed by Dominik Schrempf.- Compiled on April 27, 2020, at 12:08 pm, UTC.+ Compiled on July 30, 2020, at 12:40 pm, UTC. Usage: tlynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME] [-f|--force] COMMAND@@ -496,26 +387,10 @@ 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). + 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.- - 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 ## Compare@@ -524,12 +399,12 @@ tlynx compare --help - ELynx Suite version 0.2.1.+ ELynx Suite version 0.3.0. Developed by Dominik Schrempf.- Compiled on April 27, 2020, at 12:08 pm, UTC.+ Compiled on July 30, 2020, at 12:40 pm, UTC. Usage: tlynx compare [-n|--normalize] [-b|--bipartitions] [-t|--intersect] - [-f|--newick-format FORMAT] NAME+ [-f|--newick-format FORMAT] NAMES Compare two phylogenetic trees (compute distances and branch-wise differences). Available options:@@ -540,27 +415,11 @@ -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+ Newick tree format: Standard, IqTree, or RevBayes;+ default: Standard; for detailed help, see 'tlynx+ --help'+ NAMES Tree files -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@@ -569,9 +428,9 @@ tlynx examine --help - ELynx Suite version 0.2.1.+ ELynx Suite version 0.3.0. Developed by Dominik Schrempf.- Compiled on April 27, 2020, at 12:08 pm, UTC.+ Compiled on July 30, 2020, at 12:40 pm, UTC. Usage: tlynx examine INPUT-FILE [-f|--newick-format FORMAT] Compute summary statistics of phylogenetic trees.@@ -581,26 +440,10 @@ -V,--version Show version INPUT-FILE Read trees from INPUT-FILE -f,--newick-format FORMAT- Newick tree format; see 'tlynx- --help' (default: Standard)+ Newick tree format: Standard, IqTree, or RevBayes;+ default: Standard; for detailed help, see 'tlynx+ --help' -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@@ -609,9 +452,9 @@ tlynx simulate --help - ELynx Suite version 0.2.1.+ ELynx Suite version 0.3.0. Developed by Dominik Schrempf.- Compiled on April 27, 2020, at 12:08 pm, UTC.+ Compiled on July 30, 2020, at 12:40 pm, UTC. Usage: tlynx simulate [-t|--nTrees INT] [-n|--nLeaves INT] [-H|--height DOUBLE] [-M|--condition-on-mrca] [-l|--lambda DOUBLE] @@ -640,33 +483,17 @@ 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 Validate and (optionally) redo past ELynx analyses. - elynx --help+ elynx --help | head -n -16 - ELynx Suite version 0.2.1.+ ELynx Suite version 0.3.0. Developed by Dominik Schrempf.- Compiled on April 27, 2020, at 12:08 pm, UTC.+ Compiled on July 30, 2020, at 12:40 pm, UTC. Usage: elynx COMMAND Validate and redo past ELynx analyses@@ -678,37 +505,4 @@ 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-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,20 +1,17 @@-{- |-Module : Main-Description : Work with phylogenetic trees-Copyright : (c) Dominik Schrempf 2020-License : GPL-3.0-or-later--Maintainer : dominik.schrempf@gmail.com-Stability : unstable-Portability : portable--Creation date: Thu Sep 5 21:53:07 2019.---}-+-- |+-- Module : Main+-- Description : Work with phylogenetic trees+-- Copyright : (c) Dominik Schrempf 2020+-- License : GPL-3.0-or-later+--+-- Maintainer : dominik.schrempf@gmail.com+-- Stability : unstable+-- Portability : portable+--+-- Creation date: Thu Sep 5 21:53:07 2019. module Main where -import TLynx.TLynx+import TLynx.TLynx main :: IO () main = rTLynx
src/TLynx/Coalesce/Coalesce.hs view
@@ -1,49 +1,47 @@ {-# 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.+{-# LANGUAGE TemplateHaskell #-}+{-# LANGUAGE TupleSections #-} --}+-- TODO: MERGE WITH SIMULATE. +-- |+-- 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+ ( 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+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.Reader (ask)+import Control.Parallel.Strategies+import qualified Data.ByteString.Builder as L+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 ELynx.Data.Tree+import ELynx.Export.Tree.Newick (toNewick)+import ELynx.Simulate.Coalescent (simulate)+import ELynx.Tools+import System.Random.MWC (initialize)+import TLynx.Coalesce.Options+import TLynx.Simulate.Simulate (nSubSamples) -- | Simulate phylogenetic trees. coalesce :: ELynx CoalesceArguments ()@@ -57,55 +55,91 @@ $(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+ Just _ -> simulateAndSubSampleNTreesConcurrently+ let ls =+ if s+ then parMap rpar (formatNChildSumStat . toNChildSumStat) trs+ else parMap rpar toNewick (map lengthToPhyloTree trs) let res = L.unlines ls out "simulated trees" res ".tree" -simulateNTreesConcurrently :: ELynx CoalesceArguments [Tree (PhyloLabel Int)]+simulateNTreesConcurrently :: ELynx CoalesceArguments (Forest Length Int) simulateNTreesConcurrently = do (CoalesceArguments nT nL _ _ (Fixed s)) <- local <$> ask- c <- liftIO getNumCapabilities+ 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)+ trss <-+ liftIO $+ mapConcurrently+ (\(n, g) -> replicateM n $ simulate nL g)+ (zip chunks gs) return $ concat trss -simulateAndSubSampleNTreesConcurrently- :: ELynx CoalesceArguments [Tree (PhyloLabel Int)]+simulateAndSubSampleNTreesConcurrently ::+ ELynx CoalesceArguments (Forest Length 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 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."+ 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."+ $(logInfo) $ LT.toStrict $ LT.decodeUtf8 $ toNewick $ lengthToPhyloTree 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)+ trss <-+ liftIO $+ mapConcurrently+ (\(nSamples, g) -> nSubSamples nSamples lvs nL tr g)+ (zip chunks gs) let trs = catMaybes $ concat trss return $ map prune trs++-- | Pair of branch length with number of extant children.+type BrLnNChildren = (BranchLength, Int)++-- | Possible summary statistic of phylogenetic trees. A list of tuples+-- (BranchLength, NumberOfExtantChildrenBelowThisBranch).+type NChildSumStat = [BrLnNChildren]++-- | Format the summary statistics in the following form:+-- @+-- nLeaves1 branchLength1+-- nLeaves2 branchLength2+-- ....+formatNChildSumStat :: NChildSumStat -> L.ByteString+formatNChildSumStat s =+ L.toLazyByteString . mconcat $ map formatNChildSumStatLine s++formatNChildSumStatLine :: BrLnNChildren -> L.Builder+formatNChildSumStatLine (l, n) =+ L.intDec n <> L.char8 ' ' <> L.doubleDec l <> L.char8 '\n'++-- | Compute NChilSumStat for a phylogenetic tree.+toNChildSumStat :: Measurable e => Tree e a -> NChildSumStat+toNChildSumStat (Node br _ []) = [(getLen br, 1)]+toNChildSumStat (Node br _ ts) = (getLen br, sumCh) : concat nChSS+ where+ nChSS = map toNChildSumStat ts+ sumCh = sum $ map (snd . head) nChSS++-- TODO: MERGE WITH SIMULATE.
src/TLynx/Coalesce/Options.hs view
@@ -1,39 +1,40 @@ {-# 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+-- 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+ ( CoalesceArguments (..),+ coalesceArguments,+ reportCoalesceArguments,+ coalesceFooter, ) where -import Data.List-import Options.Applicative--import ELynx.Tools+import Data.List+import ELynx.Tools+import Options.Applicative -- | 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'.+ { -- | Number of simulated trees.+ argsNTrees :: Int,+ -- | Number of leaves per tree.+ argsNLeaves :: Int,+ -- | Perform random sub-sampling with given sampling rate.+ argsRho :: Maybe Double,+ -- | Only print summary statistics?+ argsSumStat :: Bool,+ -- | Seed of NRG, random if 'Nothing'.+ argsSeed :: Seed } deriving (Eq, Show, Generic) @@ -41,8 +42,8 @@ inFiles _ = [] outSuffixes _ = [".tree"] getSeed = Just . argsSeed- setSeed a s = a { argsSeed = Fixed s }- parser = coalesceArguments+ 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)."@@ -55,17 +56,18 @@ -- | 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+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@@ -79,43 +81,45 @@ nTreeOpt :: Parser Int nTreeOpt =- option auto- $ long "nTrees"- <> short 't'- <> metavar "INT"- <> value 10- <> showDefault- <> help "Number of trees"+ 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"+ 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)"+ 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"+ 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."+ [ "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
@@ -1,117 +1,112 @@ {-# 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.---}+{-# 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+ ( 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(..)- )+import Control.Monad+import Control.Monad.IO.Class+import Control.Monad.Logger+import Control.Monad.Trans.Reader (ask)+import Data.ByteString.Lazy.Char8 (ByteString)+import qualified Data.ByteString.Lazy.Char8 as L+import Data.List (intercalate)+import qualified Data.Map as M+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 ELynx.Data.Tree+import ELynx.Export.Tree.Newick (toNewick)+import ELynx.Import.Tree.Newick+import ELynx.Tools+ ( Arguments (..),+ ELynx,+ GlobalArguments (..),+ outHandle,+ parseFileWith,+ )+import Graphics.Gnuplot.Simple+import System.IO+import TLynx.Compare.Options+import Text.Printf -treesOneFile- :: FilePath- -> ELynx- CompareArguments- (Tree (PhyloLabel L.ByteString), Tree (PhyloLabel L.ByteString))+treesOneFile ::+ FilePath ->+ ELynx+ CompareArguments+ (Tree PhyloStrict ByteString, Tree PhyloStrict ByteString) treesOneFile tf = do nwF <- argsNewickFormat . local <$> ask $(logInfo) $ T.pack $ "Parse file '" ++ tf ++ "'."- ts <- liftIO $ parseFileWith (manyNewick nwF) tf+ ts <- liftIO $ parseFileWith (someNewick 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)+ EQ ->+ return+ ( either error id $ toStrictTree $ head ts,+ either error id $ toStrictTree $ head . tail $ ts+ ) -treesTwoFiles- :: FilePath- -> FilePath- -> ELynx- CompareArguments- (Tree (PhyloLabel L.ByteString), Tree (PhyloLabel L.ByteString))+treesTwoFiles ::+ FilePath ->+ FilePath ->+ ELynx+ CompareArguments+ (Tree PhyloStrict ByteString, Tree PhyloStrict 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)+ return (either error id $ toStrictTree t1, either error id $ toStrictTree t2) -- | More detailed comparison of two trees. compareCmd :: ELynx CompareArguments () compareCmd = do- l <- local <$> ask+ l <- local <$> ask -- Determine output handle (stdout or file). outH <- outHandle "results" ".out"- -- Read input. let inFiles = argsInFiles l- nFiles = length inFiles+ 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 $ L.hPutStrLn outH $ toNewick $ fromStrictTree tr1 liftIO $ hPutStrLn outH "Tree 2:"- liftIO $ L.hPutStrLn outH $ toNewick tr2+ liftIO $ L.hPutStrLn outH $ toNewick $ fromStrictTree 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)-+ (t1, t2) <-+ if argsIntersect l+ then do+ let [x, y] = either error id $ intersect [tr1, tr2]+ liftIO $ hPutStrLn outH "Intersected trees are:"+ liftIO $ L.hPutStrLn outH $ toNewick $ fromStrictTree x+ liftIO $ L.hPutStrLn outH $ toNewick $ fromStrictTree y+ return (x, y)+ else return (tr1, tr2) -- Check input (moved to library functions). -- let lvs1 = leaves t1 -- lvs2 = leaves t2@@ -127,35 +122,50 @@ -- 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)+ liftIO $+ T.hPutStrLn outH $+ formatD+ "Symmetric"+ (T.pack $ show $ symmetric t1 t2)+ liftIO $+ T.hPutStrLn outH $+ formatD+ "Branch score"+ (T.pack $ show $ branchScore 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'))+ $(logDebug) $ E.decodeUtf8 $ L.toStrict $ toNewick $ fromStrictTree t1'+ $(logDebug) $ E.decodeUtf8 $ L.toStrict $ toNewick $ fromStrictTree 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 (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')@@ -163,112 +173,109 @@ -- 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)"+ ( do+ let bp1 = either error id $ bipartitions (fmap getName t1)+ bp2 = either error id $ 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)+ )+ -- 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)+ )+ -- 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 getLen $ either error id $ bipartitionToBranch t1+ bpToBrLen2 = M.map getLen $ either error id $ bipartitionToBranch t2+ liftIO $+ hPutStrLn+ outH+ "Common bipartitions and their respective differences in branch lengths."+ -- Header.+ liftIO $ hPutStrLn outH header+ forM_+ bpCommon+ ( liftIO+ . hPutStrLn outH+ . getCommonBpStr 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 [%]"]+ 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+getCommonBpStr ::+ (Ord a, Show a, Fractional b, PrintfArg b) =>+ M.Map (Bipartition a) b ->+ M.Map (Bipartition a) b ->+ Bipartition a ->+ String+getCommonBpStr 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 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)+ 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
@@ -1,49 +1,47 @@ {-# 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+-- 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+ ( CompareArguments (..),+ compareArguments, ) where -import Options.Applicative--import ELynx.Tools--import TLynx.Parsers+import ELynx.Tools+import Options.Applicative+import TLynx.Parsers -- | Arguments of compare command. data CompareArguments = CompareArguments- { argsNormalize :: Bool- , argsBipartitions :: Bool- , argsIntersect :: Bool- , argsNewickFormat :: NewickFormat- , argsInFiles :: [FilePath] }+ { 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.++ -- 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 []++ -- outSuffixes a = ".out" : if argsBipartitions a then [".eps", ".pdf"] else [] getSeed _ = Nothing setSeed a _ = a- parser = compareArguments+ parser = compareArguments cmdName = "compare" cmdDsc = [ "Compare two phylogenetic trees (compute distances and branch-wise differences)."@@ -61,23 +59,29 @@ <*> bipartitions <*> intersect <*> newickFormat- <*> some file+ <*> file normalize :: Parser Bool-normalize = switch $ long "normalize" <> short 'n' <> help- "Normalize trees before comparison"+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"+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"+ 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"+file :: Parser [FilePath]+file = some $ strArgument $ metavar "NAMES" <> help "Tree files"
src/TLynx/Connect/Connect.hs view
@@ -1,116 +1,141 @@ {-# 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.---}+{-# 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+ ( 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 Control.Monad.IO.Class+import Control.Monad.Logger+import Control.Monad.Trans.Reader (ask)+import Data.ByteString.Lazy.Char8 (ByteString)+import qualified Data.ByteString.Lazy.Char8 as L+import qualified Data.Set as S+import ELynx.Data.Tree+import ELynx.Export.Tree.Newick (toNewick)+import ELynx.Import.Tree.Newick+ ( oneNewick,+ someNewick,+ )+import ELynx.Tools+ ( Arguments (..),+ ELynx,+ fromBs,+ outHandle,+ parseFileWith,+ tShow,+ )+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 with a branch in all possible ways.+--+-- Introduce a branch between two trees. If the trees have @n>2@, and @m>2@+-- nodes, respectively, there are (n-2)*(m-2) ways to connect them.+--+-- A base node label has to be given which will be used wherever the new node is+-- introduced.+--+-- Return 'Left' if one tree has a non-bifurcating root node.+connect :: (Semigroup e, Splittable e) => e -> a -> Tree e a -> Tree e a -> Either String (Forest e a)+connect br lb l r = do+ ls <- roots l+ rs <- roots r+ return [Node br lb [x, y] | x <- ls, y <- rs] --- | Connect two trees honoring possible constraints. See 'connect'.+-- | Connect two trees honoring possible constraints.+--+-- Introduce a branch between two trees. If the trees have @n>2@, and @m>2@+-- nodes, respectively, there are (n-2)*(m-2) ways to connect them. connectCmd :: ELynx ConnectArguments () connectCmd = do lArgs <- local <$> ask- outH <- outHandle "results" ".out"-+ outH <- outHandle "results" ".out" -- Do we have constraints or not? let cs = constraints lArgs- l = inFileA lArgs- r = inFileB lArgs+ l = inFileA lArgs+ r = inFileB lArgs case cs of Nothing -> connectOnly outH l r- Just c -> connectAndFilter outH c 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)+connectTrees ::+ Tree Length ByteString ->+ Tree Length ByteString ->+ Forest Length ByteString+connectTrees t = either error id . connect 0 "root" t 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+-- Get groups induced by multifurcations. Collect the leaves of all trees+-- induced by multifurcations.+multifurcatingGroups :: Tree e a -> [[a]]+multifurcatingGroups (Node _ _ []) = []+multifurcatingGroups (Node _ _ [x]) = multifurcatingGroups x+multifurcatingGroups (Node _ _ [x, y]) = multifurcatingGroups x ++ multifurcatingGroups y+multifurcatingGroups t = leaves t : concatMap multifurcatingGroups (forest t)++compatibleAll :: (Show a, Ord a) => Tree e a -> [Constraint a] -> Bool+compatibleAll t@(Node _ _ [l, r]) cs =+ all (compatible (either error id $ partition l)) (map getP cs)+ && all (compatible (either error id $ partition r)) (map getP cs)+ where+ lvs = S.fromList $ leaves t+ getP x = either error id $ mp [x, lvs S.\\ x] compatibleAll _ _ = error "Tree is not bifurcating." -compatibleWith- :: (Show b, Ord b) => (a -> b) -> [Constraint a] -> Tree a -> Bool+compatibleWith ::+ (Show b, Ord b) => (a -> b) -> [Constraint a] -> Tree e 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 ::+ FilePath ->+ FilePath ->+ ELynx+ ConnectArguments+ (Tree Length ByteString, Tree Length ByteString) parseTrees l r = do nwF <- nwFormat . local <$> ask- tl <- liftIO $ parseFileWith (oneNewick nwF) l- tr <- liftIO $ parseFileWith (oneNewick nwF) r+ 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)+ return (either error id $ phyloToLengthTree tl, either error id $ phyloToLengthTree 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+ liftIO $ L.hPutStr h $ L.unlines $ map (toNewick . lengthToPhyloTree) ts -connectAndFilter- :: Handle -> FilePath -> FilePath -> FilePath -> ELynx ConnectArguments ()+connectAndFilter ::+ Handle -> FilePath -> FilePath -> FilePath -> ELynx ConnectArguments () connectAndFilter h c l r = do nwF <- nwFormat . local <$> ask- cts <- liftIO $ parseFileWith (manyNewick nwF) c+ cts <- liftIO $ parseFileWith (someNewick 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)]+ let ts = connectTrees tl tr+ cs = map S.fromList $ concatMap multifurcatingGroups cts :: [Constraint 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'+ liftIO $ L.hPutStr h $ L.unlines $ map (toNewick . lengthToPhyloTree) ts'
src/TLynx/Connect/Options.hs view
@@ -1,37 +1,33 @@ {-# 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+-- 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+ ( ConnectArguments (..),+ connectArguments, ) where -import Options.Applicative--import ELynx.Tools--import TLynx.Parsers+import ELynx.Tools+import Options.Applicative+import TLynx.Parsers -- | Arguments of connect command. data ConnectArguments = ConnectArguments- { nwFormat :: NewickFormat- , constraints :: Maybe FilePath- , inFileA :: FilePath- , inFileB :: FilePath }+ { nwFormat :: NewickFormat,+ constraints :: Maybe FilePath,+ inFileA :: FilePath,+ inFileB :: FilePath+ } deriving (Eq, Show, Generic) instance Reproducible ConnectArguments where@@ -39,7 +35,7 @@ outSuffixes _ = [".out"] getSeed _ = Nothing setSeed a _ = a- parser = connectArguments+ parser = connectArguments cmdName = "connect" cmdDsc = [ "Connect two phylogenetic trees in all ways (possibly honoring constraints)."@@ -56,17 +52,23 @@ 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"+ 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"+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"+fileB =+ strArgument $+ metavar "TREE-FILE-B"+ <> help+ "File containing the second Newick tree"
src/TLynx/Distance/Distance.hs view
@@ -1,68 +1,52 @@ {-# 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.---}+{-# 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+ ( 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(..)- )+import Control.Monad+ ( unless,+ when,+ )+import Control.Monad.IO.Class+import Control.Monad.Logger+ ( logDebug,+ logInfo,+ )+import Control.Monad.Trans.Class+import Control.Monad.Trans.Reader hiding (local)+import Data.Bifunctor+import qualified Data.ByteString.Lazy.Char8 as L+import Data.List hiding (intersect)+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 qualified Data.Vector.Unboxed as V+import ELynx.Data.Tree+import ELynx.Export.Tree.Newick+import ELynx.Import.Tree.Newick+import ELynx.Tools+import Statistics.Sample+import System.IO+import TLynx.Distance.Options+import Text.Printf median :: Ord a => [a] -> a median xs = sort xs !! l2 where l2 = length xs `div` 2@@ -72,18 +56,45 @@ 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)+ 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 ::+ (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)+ 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 pairwise distances of a list of input trees. Use given distance+-- measure. Returns a triple, the first two elements are the indices of the+-- compared trees, the third is the distance.+pairwise ::+ -- | Distance function+ (a -> a -> b) ->+ -- | Input trees+ [a] ->+ -- | (index i, index j, distance i j)+ [(Int, Int, b)]+pairwise dist trs =+ [ (i, j, dist x y)+ | (i : is, x : xs) <- zip (tails [0 ..]) (tails trs),+ (j, y) <- zip is xs+ ]++-- -- Compute distances between adjacent pairs of a list of input trees. Use given+-- -- distance measure.+-- adjacent ::+-- -- | Distance function.+-- (a -> a -> b) ->+-- -- | Input values.+-- [a] ->+-- [b]+-- adjacent dist trs = [dist x y | (x, y) <- zip trs (tail trs)]+ -- | Compute distance functions between phylogenetic trees. distance :: ELynx DistanceArguments () distance = do@@ -95,17 +106,17 @@ let mname = argsMasterTreeFile l mtree <- case mname of Nothing -> return Nothing- Just f -> do+ 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."+ [] -> error "No tree input files given." [tf] -> do $(logInfo) "Read trees from single file."- ts <- liftIO $ parseFileWith (manyNewick nwFormat) tf+ ts <- liftIO $ parseFileWith (someNewick nwFormat) tf $(logInfo) $ tShow (length ts) <> " trees found in file." $(logInfo) "Trees are indexed with integers." return (ts, map show [0 .. length ts - 1])@@ -114,99 +125,111 @@ 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 && 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- ++ "."+ $(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+ let distanceMeasure' ::+ (Measurable e1, Measurable e2) =>+ Tree e1 L.ByteString ->+ Tree e2 L.ByteString ->+ Double+ distanceMeasure' t1 t2 = either error id $ case dist of+ Symmetric -> second fromIntegral $ symmetric t1 t2+ IncompatibleSplit _ -> second fromIntegral $ incompatibleSplits t1 t2+ BranchScore -> branchScore t1 t2 -- 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'-+ when (argsIntersect l) $+ $(logInfo) "Intersect trees before calculation of distances."+ let distanceMeasure =+ if argsIntersect l+ then+ ( \t1 t2 ->+ let [t1', t2'] = either error id $ intersect [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-+ when (argsNormalize l) $+ $(logInfo) "Normalize trees before calculation of distances."+ let normalizeF = if argsNormalize l then normalizeBranchLengths 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-+ _ -> id -- The trees can be prepared now.- let trees' = map (collapseF . normalizeF) trees-+ let trees' :: Forest PhyloStrict L.ByteString+ trees' = map (collapseF . normalizeF . either error id . toStrictTree) trees $(logDebug) "The prepared trees are:"- $(logDebug) $ LT.toStrict $ LT.decodeUtf8 $ L.unlines $ map toNewick trees'+ $(logDebug) $ LT.toStrict $ LT.decodeUtf8 $ L.unlines $ map (toNewick . fromStrictTree) 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+ Just pt ->+ let t = (either error id . toStrictTree) pt+ in [(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))+ 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)- )-+ 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
@@ -1,50 +1,50 @@ {-# 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+-- 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+ ( 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+import Data.Scientific (toRealFloat)+import Data.Void+import ELynx.Tools+import Options.Applicative+import TLynx.Parsers+import Text.Megaparsec+ ( Parsec,+ eof,+ try,+ )+import Text.Megaparsec.Char+ ( char,+ string,+ )+import Text.Megaparsec.Char.Lexer (scientific)+import Text.Printf -- | 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.+data DistanceMeasure+ = -- | Symmetric distance.+ Symmetric+ | -- | Incompatible split distance; collapse nodes+ -- with branch support below given value.+ IncompatibleSplit Double+ | -- | Branch score distance.+ BranchScore deriving (Eq, Generic) instance FromJSON DistanceMeasure@@ -52,33 +52,33 @@ instance ToJSON DistanceMeasure instance Show DistanceMeasure where- show Symmetric = "Symmetric"+ show Symmetric = "Symmetric" show (IncompatibleSplit c) = "Incompatible Split (" ++ printf "%.2f" c ++ ")"- show BranchScore = "Branch Score"+ 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]+ { 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+ Just f -> f : argsInFiles a outSuffixes _ = [".out"] getSeed _ = Nothing setSeed = const- parser = distanceArguments+ parser = distanceArguments cmdName = "distance"- cmdDsc = ["Compute distances between many phylogenetic trees."]- cmdFtr = distanceFooter+ cmdDsc = ["Compute distances between many phylogenetic trees."]+ cmdFtr = distanceFooter instance FromJSON DistanceArguments @@ -98,19 +98,19 @@ 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."+ 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"+ 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@@ -143,42 +143,44 @@ 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)"+ 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"+ 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"+ 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"+ 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."+ [ "Distance measures:",+ " symmetric Symmetric distance (Robinson-Foulds distance).",+ " incompatible-split[VAL] Incompatible split distance. Collapse branches with (normalized)",+ " support less than 0.0<=VAL<=1.0 before distance calculation;",+ " if, let's say, VAL>0.7, only well supported differences contribute",+ " to the total distance.",+ " branch-score Branch score distance." ]
src/TLynx/Examine/Examine.hs view
@@ -1,65 +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.---}+{-# 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+ ( 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+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.Containers.ListUtils (nubOrd)+import Data.List ((\\))+import qualified Data.Text as T+import ELynx.Data.Tree+import ELynx.Import.Tree.Newick+import ELynx.Tools+import System.IO+ ( Handle,+ hPutStrLn,+ )+import TLynx.Examine.Options -readTrees :: FilePath -> ELynx ExamineArguments [Tree (PhyloLabel L.ByteString)]+readTrees :: FilePath -> ELynx ExamineArguments (Forest Phylo L.ByteString) readTrees fp = do $(logInfo) $ T.pack $ "Read tree(s) from file " <> fp <> "." nf <- argsNewickFormat . local <$> ask- liftIO $ parseFileWith (manyNewick nf) fp+ liftIO $ parseFileWith (someNewick nf) fp -examineTree :: (Measurable a, Named a) => Handle -> Tree a -> IO ()+examineTree :: Named a => Handle -> Tree Phylo 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+ hPutStrLn h $ "Number of leaves: " ++ show (length lvs)+ let l = phyloToLengthTree t+ case l of+ Left _ -> hPutStrLn h "Branch lengths not available."+ Right t' -> L.hPutStrLn h $ summarizeBranchLengths t'+ unless (null dups) $+ hPutStrLn h ""+ >> hPutStrLn+ h+ ("Duplicate leaves: " ++ show dups)+ where+ lvs = map getName $ leaves t+ dups = lvs \\ nubOrd lvs -- | Examine phylogenetic trees. examine :: ELynx ExamineArguments () examine = do l <- local <$> ask let inFn = argsInFile l- trs <- readTrees inFn+ trs <- readTrees inFn outH <- outHandle "results" ".out" liftIO $ mapM_ (examineTree outH) trs
src/TLynx/Examine/Options.hs view
@@ -1,35 +1,31 @@ {-# 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+-- 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+ ( ExamineArguments (..),+ examineArguments, ) where -import Options.Applicative--import ELynx.Tools--import TLynx.Parsers+import ELynx.Tools+import Options.Applicative+import TLynx.Parsers -- | Arguments needed to examine phylogenetic trees. data ExamineArguments = ExamineArguments- { argsInFile :: FilePath- , argsNewickFormat :: NewickFormat }+ { argsInFile :: FilePath,+ argsNewickFormat :: NewickFormat+ } deriving (Eq, Show, Generic) instance Reproducible ExamineArguments where@@ -37,9 +33,9 @@ outSuffixes _ = [".out"] getSeed _ = Nothing setSeed = const- parser = examineArguments+ parser = examineArguments cmdName = "examine"- cmdDsc = ["Compute summary statistics of phylogenetic trees."]+ cmdDsc = ["Compute summary statistics of phylogenetic trees."] instance FromJSON ExamineArguments
src/TLynx/Options.hs view
@@ -1,42 +1,37 @@ {-# 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+-- 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+ ( 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 )+import ELynx.Tools+import Options.Applicative+import TLynx.Coalesce.Options+import TLynx.Compare.Options+import TLynx.Connect.Options+import TLynx.Distance.Options+import TLynx.Examine.Options+import TLynx.Parsers (newickHelp)+import TLynx.Shuffle.Options+import TLynx.Simulate.Options -- | The different TLynx commands and their arguments.-data CommandArguments =- Coalesce CoalesceArguments+data CommandArguments+ = Coalesce CoalesceArguments | Compare CompareArguments | Connect ConnectArguments | Distance DistanceArguments@@ -47,45 +42,46 @@ instance Reproducible CommandArguments where inFiles (Coalesce a) = inFiles a- inFiles (Compare a) = inFiles a- inFiles (Connect 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 (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 (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 (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 (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 (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 (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 (Examine a) = Examine . setSeed a+ setSeed (Shuffle a) = Shuffle . setSeed a setSeed (Simulate a) = Simulate . setSeed a - parser = commandArguments+ parser = commandArguments cmdName = "tlynx" - cmdDsc = ["Compare, examine, and simulate phylogenetic trees."]+ cmdDsc = ["Compare, examine, and simulate phylogenetic trees."] - cmdFtr = "Available tree file formats:" : indent newickHelp- where indent = map (" " ++)+ cmdFtr = "" : "Available tree file formats:" : indent newickHelp+ where+ indent = map (" " ++) instance FromJSON CommandArguments @@ -114,11 +110,11 @@ commandArguments :: Parser CommandArguments commandArguments =- hsubparser- $ coalesceCommand- <> compareCommand- <> connectCommand- <> distanceCommand- <> examineCommand- <> shuffleCommand- <> simulateCommand+ hsubparser $+ coalesceCommand+ <> compareCommand+ <> connectCommand+ <> distanceCommand+ <> examineCommand+ <> shuffleCommand+ <> simulateCommand
src/TLynx/Parsers.hs view
@@ -1,43 +1,44 @@-{- |-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+-- 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+ ( NewickFormat,+ newickFormat,+ newickHelp, ) where -import Options.Applicative--import ELynx.Import.Tree.Newick ( NewickFormat(..)- , description- )-import ELynx.Tools ( allValues )+import Data.List (intercalate)+import ELynx.Import.Tree.Newick+ ( NewickFormat (..),+ description,+ )+import ELynx.Tools (allValues)+import Options.Applicative -- | 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'"+ option auto $+ long "newick-format"+ <> short 'f'+ <> metavar "FORMAT"+ <> value Standard+ <> help ("Newick tree format: " ++ nwlist ++ "; default: Standard; for detailed help, see 'tlynx --help'")+ where+ nwfs = map show (allValues :: [NewickFormat])+ nwlist = intercalate ", " (init nwfs) <> ", or " <> last nwfs -- | Help for different 'NewickFormat's. newickHelp :: [String] newickHelp = map (toListItem . description) (allValues :: [NewickFormat])- where toListItem = ("- Newick " ++)+ where+ toListItem = ("- Newick " ++)
src/TLynx/Shuffle/Options.hs view
@@ -1,45 +1,41 @@ {-# 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+-- 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+ ( ShuffleArguments (..),+ shuffleArguments, ) where -import Options.Applicative--import ELynx.Tools--import TLynx.Parsers+import ELynx.Tools+import Options.Applicative+import TLynx.Parsers -- | Arguments of shuffle command. data ShuffleArguments = ShuffleArguments- { nwFormat :: NewickFormat- , nReplicates :: Int- , inFile :: FilePath- , argsSeed :: Seed }+ { 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+ setSeed a s = a {argsSeed = Fixed s}+ parser = shuffleArguments cmdName = "shuffle" cmdDsc = [ "Shuffle a phylogenetic tree (keep coalescent times, but shuffle topology and leaves)."@@ -55,8 +51,10 @@ n :: Parser Int n =- option auto $ long "replicates" <> short 'n' <> metavar "N" <> value 1 <> help- "Number of trees to generate"+ option auto $+ long "replicates" <> short 'n' <> metavar "N" <> value 1+ <> help+ "Number of trees to generate" file :: Parser FilePath file =
src/TLynx/Shuffle/Shuffle.hs view
@@ -1,68 +1,42 @@ {-# 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.---}+{-# 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+ ( 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+import qualified Control.Comonad as C+import Control.Monad (when)+import Control.Monad.IO.Class (liftIO)+import Control.Monad.Logger (logDebug, logInfo)+import Control.Monad.Trans.Reader (ask)+import qualified Data.ByteString.Lazy.Char8 as L+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+import System.IO (hClose)+import System.Random.MWC (GenIO, initialize)+import TLynx.Shuffle.Options -- | Shuffle a tree. Get all coalescent times, shuffle them. Get all leaves, -- shuffle them. Connect the shuffled leaves with the shuffled coalescent times.@@ -73,88 +47,42 @@ l <- local <$> ask h <- outHandle "results" ".tree" let nwF = nwFormat l- t <- liftIO $ parseFileWith (oneNewick nwF) (inFile l)+ tPhylo <- 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."- )-+ $(logInfo) $ fromBs $ toNewick tPhylo+ let t = either error id $ phyloToLengthTree tPhylo -- 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 && 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+ let cs = filter (> 0) $ labels $ C.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."+ 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-+ ts <- liftIO $ shuffleT (nReplicates l) (height t) cs ls gen+ liftIO $ L.hPutStr h $ L.unlines $ map (toNewick . lengthToPhyloTree) 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+shuffleT ::+ Int -> -- How many?+ Double -> -- Stem length.+ [Double] -> -- Coalescent times.+ [L.ByteString] -> -- Leave names.+ GenIO ->+ IO (Forest Length L.ByteString)+shuffleT 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+ | (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
@@ -1,45 +1,51 @@ {-# 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+-- 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+ ( SimulateArguments (..),+ simulateArguments,+ reportSimulateArguments,+ simulateFooter, ) where -import Data.List-import Options.Applicative--import ELynx.Tools+import Data.List+import ELynx.Tools+import Options.Applicative -- | 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'.+ { -- | Simulated trees.+ argsNTrees :: Int,+ -- | Number of leaves.+ argsNLeaves :: Int,+ -- | Tree height (time to origin or MRCA).+ argsHeight :: Maybe Double,+ -- | False: condition on origin; True:+ -- condition on MRCA.+ argsConditionMRCA :: Bool,+ -- | Birth rate.+ argsLambda :: Double,+ -- | Death rate.+ argsMu :: Double,+ -- | Smapling rate.+ argsRho :: Double,+ -- | Perform actual sub-sampling.+ argsSubSample :: Bool,+ -- | Only print summary statistics?+ argsSumStat :: Bool,+ -- | Seed of NRG, random if 'Nothing'.+ argsSeed :: Seed } deriving (Eq, Show, Generic) @@ -47,8 +53,8 @@ inFiles _ = [] outSuffixes _ = [".tree"] getSeed = Just . argsSeed- setSeed a s = a { argsSeed = Fixed s }- parser = simulateArguments+ 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)."@@ -61,23 +67,26 @@ -- | 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"+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@@ -96,71 +105,80 @@ nTreeOpt :: Parser Int nTreeOpt =- option auto- $ long "nTrees"- <> short 't'- <> metavar "INT"- <> value 10- <> showDefault- <> help "Number of trees"+ 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"+ 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)"+ 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"+ 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"+ 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"+ 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)"+ 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."+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"+ switch $+ long "summary-statistics" <> short 's' <> showDefault+ <> help+ "Only output number of children for each branch" citation :: String citation =@@ -169,8 +187,8 @@ -- | 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+ [ "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
@@ -1,140 +1,199 @@ {-# 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---}+{-# 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+ ( simulate,+ nSubSamples, ) 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+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.Reader (ask)+import Control.Parallel.Strategies+import qualified Data.ByteString.Builder as L+import qualified Data.ByteString.Lazy.Char8 as L+import Data.Foldable (toList)+import Data.Maybe+import qualified Data.Sequence as Seq+import qualified Data.Set as Set+import qualified Data.Text as T+import qualified Data.Text.Lazy as LT+import qualified Data.Text.Lazy.Encoding as LT+import ELynx.Data.Tree+import ELynx.Export.Tree.Newick (toNewick)+import ELynx.Simulate.PointProcess+ ( TimeSpec,+ simulateNReconstructedTrees,+ simulateReconstructedTree,+ )+import ELynx.Tools+import System.Random.MWC+ ( GenIO,+ initialize,+ )+import TLynx.Simulate.Options -- | 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+ 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)."+ 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 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 $ map lengthToPhyloTree 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 ::+ Int ->+ Double ->+ Double ->+ Double ->+ TimeSpec ->+ [Int] ->+ [GenIO] ->+ ELynx SimulateArguments (Forest Length 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)+ 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 ::+ Int ->+ Double ->+ Double ->+ Double ->+ TimeSpec ->+ [Int] ->+ [GenIO] ->+ ELynx SimulateArguments (Forest Length 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."+ 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."+ $(logInfo) $ LT.toStrict $ LT.decodeUtf8 $ toNewick $ lengthToPhyloTree 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)+ trss <-+ liftIO $+ mapConcurrently+ (\(nSamples, g) -> nSubSamples nSamples lvs nLeaves tr g)+ (zip chunks gs) let trs = catMaybes $ concat trss return $ map prune trs++-- | Extract a random subtree with @N@ leaves of a tree with @M@ leaves, where+-- @M>N@ (otherwise error). The complete list of leaves (names are assumed to be+-- unique) has to be provided as a 'Seq.Seq', and a 'Seq.Set', so that fast+-- sub-sampling as well as lookup are fast and so that these data structures do+-- not have to be recomputed when many sub-samples are requested.+nSubSamples ::+ Ord a =>+ Int ->+ Seq.Seq a ->+ Int ->+ Tree e a ->+ GenIO ->+ IO [Maybe (Tree e a)]+nSubSamples m lvs n tree g+ | Seq.length lvs < n =+ error+ "Given list of leaves is shorter than requested number of leaves."+ | otherwise = do+ lss <- grabble (toList lvs) m n g+ let lsSets = map Set.fromList lss+ return [dropLeavesWith (`Set.notMember` ls) tree | ls <- lsSets]++-- | Pair of branch length with number of extant children.+type BrLnNChildren = (BranchLength, Int)++-- | Possible summary statistic of phylogenetic trees. A list of tuples+-- (BranchLength, NumberOfExtantChildrenBelowThisBranch).+type NChildSumStat = [BrLnNChildren]++-- | Format the summary statistics in the following form:+-- @+-- nLeaves1 branchLength1+-- nLeaves2 branchLength2+-- ....+formatNChildSumStat :: NChildSumStat -> L.ByteString+formatNChildSumStat s =+ L.toLazyByteString . mconcat $ map formatNChildSumStatLine s++formatNChildSumStatLine :: BrLnNChildren -> L.Builder+formatNChildSumStatLine (l, n) =+ L.intDec n <> L.char8 ' ' <> L.doubleDec l <> L.char8 '\n'++-- | Compute NChilSumStat for a phylogenetic tree.+toNChildSumStat :: Measurable e => Tree e a -> NChildSumStat+toNChildSumStat (Node br _ []) = [(getLen br, 1)]+toNChildSumStat (Node br _ ts) = (getLen br, sumCh) : concat nChSS+ where+ nChSS = map toNChildSumStat ts+ sumCh = sum $ map (snd . head) nChSS
src/TLynx/TLynx.hs view
@@ -1,36 +1,32 @@-{- |-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+-- 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+ ( 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+import ELynx.Tools+import TLynx.Coalesce.Coalesce+import TLynx.Compare.Compare+import TLynx.Connect.Connect+import TLynx.Distance.Distance+import TLynx.Examine.Examine+import TLynx.Options+import TLynx.Shuffle.Shuffle+import TLynx.Simulate.Simulate -- 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
tlynx.cabal view
@@ -1,99 +1,91 @@-cabal-version: 1.12-name: tlynx-version: 0.2.2-license: GPL-3-license-file: LICENSE-copyright: Dominik Schrempf (2020)-maintainer: dominik.schrempf@gmail.com-author: Dominik Schrempf-homepage: https://github.com/dschrempf/elynx#readme-bug-reports: https://github.com/dschrempf/elynx/issues-synopsis: Handle phylogenetic trees-description:- Examine, compare, and simulate phylogenetic trees in a reproducible way. Please see the README on GitHub at <https://github.com/dschrempf/elynx>.-category: Bioinformatics-build-type: Simple+cabal-version: 2.2++-- This file has been generated from package.yaml by hpack version 0.34.2.+--+-- see: https://github.com/sol/hpack+--+-- hash: 4c935a10e73a1a6ec515b771ad749be8dcdb27e3daea89833fa5218ceafe0b47++name: tlynx+version: 0.3.0+synopsis: Handle phylogenetic trees+description: Examine, compare, and simulate phylogenetic trees in a reproducible way. Please see the README on GitHub at <https://github.com/dschrempf/elynx>.+category: Bioinformatics+homepage: https://github.com/dschrempf/elynx#readme+bug-reports: https://github.com/dschrempf/elynx/issues+author: Dominik Schrempf+maintainer: dominik.schrempf@gmail.com+copyright: Dominik Schrempf (2020)+license: GPL-3.0-or-later+license-file: LICENSE+build-type: Simple extra-source-files: README.md ChangeLog.md source-repository head- type: git- location: https://github.com/dschrempf/elynx+ type: git+ location: https://github.com/dschrempf/elynx library- exposed-modules:- 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 -fllvm- build-depends:- 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.2 && <0.3,- elynx-tree >=0.2.2 && <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.15.1.0 && <0.16,- parallel >=3.2.2.0 && <3.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.1.2 && <0.13+ exposed-modules:+ 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+ other-modules:+ Paths_tlynx+ autogen-modules:+ Paths_tlynx+ hs-source-dirs:+ src+ ghc-options: -Wall -fllvm+ build-depends:+ aeson+ , array+ , base >=4.7 && <5+ , bytestring+ , comonad+ , containers+ , elynx-tools+ , elynx-tree+ , gnuplot+ , lifted-async+ , megaparsec+ , monad-logger+ , mwc-random+ , optparse-applicative+ , parallel+ , primitive+ , scientific+ , statistics+ , text+ , transformers+ , vector+ default-language: Haskell2010 executable tlynx- main-is: Main.hs- hs-source-dirs: app- other-modules:- Paths_tlynx- default-language: Haskell2010- ghc-options: -Wall -fllvm -threaded -rtsopts -with-rtsopts=-N- build-depends:- 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.2 && <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.15.1.0 && <0.16,- parallel >=3.2.2.0 && <3.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.1.2 && <0.13+ main-is: Main.hs+ other-modules:+ Paths_tlynx+ hs-source-dirs:+ app+ ghc-options: -Wall -fllvm -threaded -rtsopts -with-rtsopts=-N+ build-depends:+ base >=4.7 && <5+ , elynx-tools+ , tlynx+ default-language: Haskell2010