elynx-seq 0.1.0 → 0.2.1
raw patch · 49 files changed
+1267/−3189 lines, 49 filesdep +aesondep −asyncdep −data-memocombinatorsdep −elynx-treedep ~basedep ~bytestringdep ~containers
Dependencies added: aeson
Dependencies removed: async, data-memocombinators, elynx-tree, hmatrix, integration, math-functions, optparse-applicative, statistics, transformers
Dependency ranges changed: base, bytestring, containers, elynx-tools, hspec-megaparsec, megaparsec, primitive, vector, vector-th-unbox
Files
- ChangeLog.md +8/−5
- README.md +375/−88
- elynx-seq.cabal +23/−62
- src/ELynx/Data/Alphabet/Alphabet.hs +44/−25
- src/ELynx/Data/Alphabet/Character.hs +9/−8
- src/ELynx/Data/Alphabet/DistributionDiversity.hs +38/−35
- src/ELynx/Data/Character/AminoAcid.hs +28/−26
- src/ELynx/Data/Character/AminoAcidI.hs +35/−33
- src/ELynx/Data/Character/AminoAcidS.hs +31/−29
- src/ELynx/Data/Character/AminoAcidX.hs +31/−29
- src/ELynx/Data/Character/BoundaryMutation.hs +0/−281
- src/ELynx/Data/Character/Character.hs +10/−11
- src/ELynx/Data/Character/Codon.hs +186/−54
- src/ELynx/Data/Character/Nucleotide.hs +12/−10
- src/ELynx/Data/Character/NucleotideI.hs +27/−25
- src/ELynx/Data/Character/NucleotideX.hs +16/−14
- src/ELynx/Data/MarkovProcess/AminoAcid.hs +0/−250
- src/ELynx/Data/MarkovProcess/CXXModels.hs +0/−128
- src/ELynx/Data/MarkovProcess/CXXModelsData.hs +0/−294
- src/ELynx/Data/MarkovProcess/EDMModel.hs +0/−38
- src/ELynx/Data/MarkovProcess/GammaRateHeterogeneity.hs +0/−109
- src/ELynx/Data/MarkovProcess/MixtureModel.hs +0/−145
- src/ELynx/Data/MarkovProcess/Nucleotide.hs +0/−81
- src/ELynx/Data/MarkovProcess/PhyloModel.hs +0/−46
- src/ELynx/Data/MarkovProcess/RateMatrix.hs +0/−95
- src/ELynx/Data/MarkovProcess/SubstitutionModel.hs +0/−119
- src/ELynx/Data/Sequence/Alignment.hs +102/−82
- src/ELynx/Data/Sequence/Defaults.hs +3/−2
- src/ELynx/Data/Sequence/Sequence.hs +87/−45
- src/ELynx/Data/Sequence/Translate.hs +21/−15
- src/ELynx/Export/Sequence/CountsFile.hs +0/−109
- src/ELynx/Export/Sequence/Fasta.hs +9/−6
- src/ELynx/Import/MarkovProcess/EDMModelPhylobayes.hs +0/−70
- src/ELynx/Import/MarkovProcess/SiteprofilesPhylobayes.hs +0/−77
- src/ELynx/Import/Sequence/Fasta.hs +24/−16
- src/ELynx/Simulate/MarkovProcess.hs +0/−76
- src/ELynx/Simulate/MarkovProcessAlongTree.hs +0/−154
- test/ELynx/Data/Alphabet/DistributionDiversitySpec.hs +21/−20
- test/ELynx/Data/MarkovProcess/AminoAcidSpec.hs +0/−154
- test/ELynx/Data/MarkovProcess/NucleotideSpec.hs +0/−36
- test/ELynx/Data/Sequence/AlignmentSpec.hs +20/−12
- test/ELynx/Data/Sequence/SequenceSpec.hs +24/−15
- test/ELynx/Data/Sequence/TranslateSpec.hs +12/−6
- test/ELynx/Export/Sequence/FastaSpec.hs +22/−13
- test/ELynx/Import/MarkovProcess/EDMModelPhylobayesSpec.hs +0/−55
- test/ELynx/Import/MarkovProcess/SiteprofilesPhylobayesSpec.hs +0/−46
- test/ELynx/Import/Sequence/FastaSpec.hs +49/−26
- test/ELynx/Simulate/MarkovProcessAlongTreeSpec.hs +0/−59
- test/Files.hs +0/−55
ChangeLog.md view
@@ -1,10 +1,13 @@ -# Table of Contents--1. [Changelog for ELynx](#org263dcac)+# Changelog for ELynx -<a id="org263dcac"></a>+## Version 0.2.1 -# Changelog for ELynx+- Validation and repetition of previous analyses is finally possible with the+ new `elynx` binary.+- A library `elynx-markov` for running Markov processes along phylogenetic trees+ has been split off `elynx-seq`. This library performs the computations when+ executing `slynx simulate ...`.+- Many other small improvements.
README.md view
@@ -2,30 +2,31 @@ # The ELynx Suite -Version: 0.1.0.+Version: 0.2.1. Reproducible evolution made easy. -The ELynx Suite is a Haskell library and a tool set for computational biology.-The goal of the ELynx Suite is reproducible research. Evolutionary sequences and-phylogenetic trees can be read, viewed, modified and simulated. Exact-specification of all options is necessary, and nothing is assumed about the data-(e.g., the type of the genetic code). The command line with all arguments is-consistently, and automatically logged. The work overhead in the beginning-usually pays off in the end.+A Haskell library and tool set for computational biology. The goal of ELynx is+reproducible research. Evolutionary sequences and phylogenetic trees can be+read, viewed, modified and simulated. The command line with all arguments is+logged consistently, and automatically. Data integrity is verified using SHA256+sums so that validation of past analyses is possible without the need to+recompute the result. -The Elynx Suite consists of three library packages and two executables providing+The Elynx Suite consists of four library packages and three executables providing a range of sub commands. The library packages are: +- **elynx-markov:** Simulate multi sequence alignments along phylogenetic trees. - **elynx-seq:** Handle evolutionary sequences and multi sequence alignments.-- **elynx-tree:** Handle phylogenetic trees. - **elynx-tools:** Tools for the provided executables.+- **elynx-tree:** Handle phylogenetic trees. The executables are: -- **slynx:** Analyze, modify, and simulate evolutionary sequences (FASTA format).-- **tlynx:** Analyze, modify, and simulate phylogenetic trees (Newick format).+- **slynx:** Analyze, modify, and simulate evolutionary sequences.+- **tlynx:** Analyze, modify, and simulate phylogenetic trees.+- **elynx:** Validate and redo past analyses. **ELynx is actively developed. We happily receive comments, ideas, feature requests, and pull requests!**@@ -67,10 +68,12 @@ slynx --help - ELynx Suite version 0.1.0. Developed by Dominik Schrempf. Compiled on January- 30, 2020, at 09:18 am, UTC.+ ELynx Suite version 0.2.1.+ Developed by Dominik Schrempf.+ Compiled on April 27, 2020, at 06:31 am, UTC. - Usage: slynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME] COMMAND+ Usage: slynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME] + [-f|--force] COMMAND Analyze, and simulate multi sequence alignments. Available options:@@ -80,21 +83,22 @@ Debug (default: Info) -o,--output-file-basename NAME Specify base name of output file+ -f,--force Ignore previous analysis and overwrite existing+ output files. Available commands: concatenate Concatenate sequences found in input files.- examine Examine sequences. If data is a multi sequence- alignment, additionally analyze columns.+ examine Examine sequences. If data is a multi sequence alignment, additionally analyze columns.+ filter-columns Filter columns of multi sequence alignments. filter-rows Filter rows (or sequences) found in input files.- filter-columns Filter columns of multi-sequence alignments. simulate Simulate multi sequence alignments. sub-sample Sub-sample columns from multi sequence alignments. translate Translate from DNA to Protein or DNAX to ProteinX. - File formats:+ Available sequence file formats: - FASTA - Alphabet types:+ Available alphabets: - DNA (nucleotides) - DNAX (nucleotides; including gaps) - DNAI (nucleotides; including gaps, and IUPAC codes)@@ -103,19 +107,20 @@ - ProteinS (amino acids; including gaps, and translation stops) - ProteinI (amino acids; including gaps, translation stops, and IUPAC codes) - The ELynx Suite- ---------------- A Haskell library and a tool set for computational biology. The goal of the- ELynx Suite is reproducible research. Evolutionary sequences and phylogenetic- trees can be read, viewed, modified and simulated. Exact specification of all- options is necessary, and nothing is assumed about the data (e.g., the type of- code). The command line with all arguments is consistently, and automatically- logged.+ ELynx+ -----+ A Haskell library and tool set for computational biology. The goal of ELynx is+ reproducible research. Evolutionary sequences and phylogenetic trees can be+ read, viewed, modified and simulated. The command line with all arguments is+ logged consistently, and automatically. Data integrity is verified using SHA256+ sums so that validation of past analyses is possible without the need to+ recompute the result. slynx Analyze, modify, and simulate evolutionary sequences. tlynx Analyze, modify, and simulate phylogenetic trees.+ elynx Validate and redo past analyses. - Get help for specific commands:+ Get help for sub commands: slynx examine --help @@ -125,13 +130,36 @@ slynx concatenate --help + ELynx Suite version 0.2.1.+ Developed by Dominik Schrempf.+ Compiled on April 27, 2020, at 06:31 am, UTC.+ Usage: slynx concatenate (-a|--alphabet NAME) INPUT-FILE Concatenate sequences found in input files. Available options:+ -h,--help Show this help text+ -V,--version Show version -a,--alphabet NAME Specify alphabet type NAME INPUT-FILE Read sequences from INPUT-FILE -h,--help Show this help text+ + + ELynx+ -----+ A Haskell library and tool set for computational biology. The goal of ELynx is+ reproducible research. Evolutionary sequences and phylogenetic trees can be+ read, viewed, modified and simulated. The command line with all arguments is+ logged consistently, and automatically. Data integrity is verified using SHA256+ sums so that validation of past analyses is possible without the need to+ recompute the result.+ + slynx Analyze, modify, and simulate evolutionary sequences.+ tlynx Analyze, modify, and simulate phylogenetic trees.+ elynx Validate and redo past analyses.+ + Get help for sub commands:+ slynx examine --help ## Examine@@ -140,15 +168,37 @@ slynx examine --help - Usage: slynx examine (-a|--alphabet NAME) [INPUT-FILE] [--per-site]- Examine sequences. If data is a multi sequence alignment, additionally analyze- columns.+ ELynx Suite version 0.2.1.+ Developed by Dominik Schrempf.+ Compiled on April 27, 2020, at 06:31 am, UTC. + Usage: slynx examine (-a|--alphabet NAME) INPUT-FILE [--per-site]+ Examine sequences. If data is a multi sequence alignment, additionally analyze columns.+ Available options:+ -h,--help Show this help text+ -V,--version Show version -a,--alphabet NAME Specify alphabet type NAME INPUT-FILE Read sequences from INPUT-FILE --per-site Report per site summary statistics -h,--help Show this help text+ + + ELynx+ -----+ A Haskell library and tool set for computational biology. The goal of ELynx is+ reproducible research. Evolutionary sequences and phylogenetic trees can be+ read, viewed, modified and simulated. The command line with all arguments is+ logged consistently, and automatically. Data integrity is verified using SHA256+ sums so that validation of past analyses is possible without the need to+ recompute the result.+ + slynx Analyze, modify, and simulate evolutionary sequences.+ tlynx Analyze, modify, and simulate phylogenetic trees.+ elynx Validate and redo past analyses.+ + Get help for sub commands:+ slynx examine --help ## Filter@@ -157,12 +207,17 @@ slynx filter-rows --help - Usage: slynx filter-rows (-a|--alphabet NAME) [INPUT-FILE]- [--longer-than LENGTH] [--shorter-than LENGTH]- [--standard-characters]+ ELynx Suite version 0.2.1.+ Developed by Dominik Schrempf.+ Compiled on April 27, 2020, at 06:31 am, UTC.+ + Usage: slynx filter-rows (-a|--alphabet NAME) INPUT-FILE [--longer-than LENGTH] + [--shorter-than LENGTH] [--standard-characters] Filter rows (or sequences) found in input files. Available options:+ -h,--help Show this help text+ -V,--version Show version -a,--alphabet NAME Specify alphabet type NAME INPUT-FILE Read sequences from INPUT-FILE --longer-than LENGTH Only keep sequences longer than LENGTH@@ -170,21 +225,61 @@ --standard-characters Only keep sequences containing at least one standard (i.e., non-IUPAC) character -h,--help Show this help text+ + + ELynx+ -----+ A Haskell library and tool set for computational biology. The goal of ELynx is+ reproducible research. Evolutionary sequences and phylogenetic trees can be+ read, viewed, modified and simulated. The command line with all arguments is+ logged consistently, and automatically. Data integrity is verified using SHA256+ sums so that validation of past analyses is possible without the need to+ recompute the result.+ + slynx Analyze, modify, and simulate evolutionary sequences.+ tlynx Analyze, modify, and simulate phylogenetic trees.+ elynx Validate and redo past analyses.+ + Get help for sub commands:+ slynx examine --help Filter columns of multi sequence alignments with `filter-columns`. slynx filter-columns --help - Usage: slynx filter-columns (-a|--alphabet NAME) [INPUT-FILE]+ ELynx Suite version 0.2.1.+ Developed by Dominik Schrempf.+ Compiled on April 27, 2020, at 06:31 am, UTC.+ + Usage: slynx filter-columns (-a|--alphabet NAME) INPUT-FILE [--standard-chars DOUBLE]- Filter columns of multi-sequence alignments.+ Filter columns of multi sequence alignments. Available options:+ -h,--help Show this help text+ -V,--version Show version -a,--alphabet NAME Specify alphabet type NAME INPUT-FILE Read sequences from INPUT-FILE --standard-chars DOUBLE Keep columns with a proportion standard (non-IUPAC) characters larger than DOUBLE in [0,1] -h,--help Show this help text+ + + ELynx+ -----+ A Haskell library and tool set for computational biology. The goal of ELynx is+ reproducible research. Evolutionary sequences and phylogenetic trees can be+ read, viewed, modified and simulated. The command line with all arguments is+ logged consistently, and automatically. Data integrity is verified using SHA256+ sums so that validation of past analyses is possible without the need to+ recompute the result.+ + slynx Analyze, modify, and simulate evolutionary sequences.+ tlynx Analyze, modify, and simulate phylogenetic trees.+ elynx Validate and redo past analyses.+ + Get help for sub commands:+ slynx examine --help ## Simulate@@ -193,16 +288,22 @@ slynx simulate --help - Usage: slynx simulate (-t|--tree-file Name) [-s|--substitution-model MODEL]- [-m|--mixture-model MODEL] [-e|--edm-file NAME]- [-p|--siteprofile-files NAMES]- [-w|--mixture-model-weights "[DOUBLE,DOUBLE,...]"]+ ELynx Suite version 0.2.1.+ Developed by Dominik Schrempf.+ Compiled on April 27, 2020, at 06:31 am, UTC.+ + Usage: slynx simulate (-t|--tree-file Name) [-s|--substitution-model MODEL] + [-m|--mixture-model MODEL] [-e|--edm-file NAME] + [-p|--siteprofile-files NAMES] + [-w|--mixture-model-weights "[DOUBLE,DOUBLE,...]"] [-g|--gamma-rate-heterogeneity "(NCAT,SHAPE)"] (-l|--length NUMBER) [-S|--seed [INT]] Simulate multi sequence alignments. Available options:- -t,--tree-file Name Read trees from file NAME+ -h,--help Show this help text+ -V,--version Show version+ -t,--tree-file Name Read tree from Newick file NAME -s,--substitution-model MODEL Set the phylogenetic substitution model; available models are shown below (mutually exclusive with -m@@ -224,17 +325,24 @@ Substitution models: -s "MODEL[PARAMETER,PARAMETER,...]{STATIONARY_DISTRIBUTION}"- Supported DNA models: JC, HKY.+ Supported DNA models: JC, F81, HKY, GTR4. For example, -s HKY[KAPPA]{DOUBLE,DOUBLE,DOUBLE,DOUBLE}- Supported Protein models: Poisson, Poisson-Custom, LG, LG-Custom, WAG, WAG-Custom.+ -s GTR4[e_AC,e_AG,e_AT,e_CG,e_CT,e_GT]{DOUBLE,DOUBLE,DOUBLE,DOUBLE}+ where the 'e_XY' are the exchangeabilities from nucleotide X to Y.+ Supported Protein models: Poisson, Poisson-Custom, LG, LG-Custom, WAG, WAG-Custom, GTR20. MODEL-Custom means that only the exchangeabilities of MODEL are used, and a custom stationary distribution is provided. For example,+ -s LG -s LG-Custom{...}+ -s GTR20[e_AR,e_AN,...]{...}+ the 'e_XY' are the exchangeabilities from amino acid X to Y (alphabetical order).+ Notes: The F81 model for DNA is equivalent to the Poisson-Custom for proteins.+ The GTR4 model for DNA is equivalent to the GTR20 for proteins. Mixture models:- -m "MIXTURE(SUBSTITUTION_MODEL_1,SUBSTITUTION_MODEL_2)"+ -m "MIXTURE(SUBSTITUTION_MODEL_1,SUBSTITUTION_MODEL_2[PARAMETERS]{STATIONARY_DISTRIBUTION},...)" For example, -m "MIXTURE(JC,HKY[6.0]{0.3,0.2,0.2,0.3})" Mixture weights have to be provided with the -w option.@@ -249,6 +357,22 @@ LG exchangeabilities with stationary distributions given in FILE. -m "EDM(LG-Custom)" -e FILE For special mixture models, mixture weights are optional.+ + ELynx+ -----+ A Haskell library and tool set for computational biology. The goal of ELynx is+ reproducible research. Evolutionary sequences and phylogenetic trees can be+ read, viewed, modified and simulated. The command line with all arguments is+ logged consistently, and automatically. Data integrity is verified using SHA256+ sums so that validation of past analyses is possible without the need to+ recompute the result.+ + slynx Analyze, modify, and simulate evolutionary sequences.+ tlynx Analyze, modify, and simulate phylogenetic trees.+ elynx Validate and redo past analyses.+ + Get help for sub commands:+ slynx examine --help ## Sub-sample@@ -257,12 +381,18 @@ slynx sub-sample --help - Usage: slynx sub-sample (-a|--alphabet NAME) [INPUT-FILE]+ ELynx Suite version 0.2.1.+ Developed by Dominik Schrempf.+ Compiled on April 27, 2020, at 06:31 am, UTC.+ + Usage: slynx sub-sample (-a|--alphabet NAME) INPUT-FILE (-n|--number-of-sites INT) (-m|--number-of-alignments INT) [-S|--seed [INT]] Sub-sample columns from multi sequence alignments. Available options:+ -h,--help Show this help text+ -V,--version Show version -a,--alphabet NAME Specify alphabet type NAME INPUT-FILE Read sequences from INPUT-FILE -n,--number-of-sites INT Number of sites randomly drawn with replacement@@ -272,8 +402,23 @@ integers with up to 256 elements (default: random) -h,--help Show this help text - Create a given number of multi sequence alignments, each of which contains a- given number of random sites drawn from the original multi sequence alignment.+ Create a given number of multi sequence alignments, each of which contains a given number of random sites drawn from the original multi sequence alignment.+ + ELynx+ -----+ A Haskell library and tool set for computational biology. The goal of ELynx is+ reproducible research. Evolutionary sequences and phylogenetic trees can be+ read, viewed, modified and simulated. The command line with all arguments is+ logged consistently, and automatically. Data integrity is verified using SHA256+ sums so that validation of past analyses is possible without the need to+ recompute the result.+ + slynx Analyze, modify, and simulate evolutionary sequences.+ tlynx Analyze, modify, and simulate phylogenetic trees.+ elynx Validate and redo past analyses.+ + Get help for sub commands:+ slynx examine --help ## Translate@@ -282,17 +427,40 @@ slynx translate --help - Usage: slynx translate (-a|--alphabet NAME) [INPUT-FILE]- (-r|--reading-frame INT) (-u|--universal-code CODE)+ ELynx Suite version 0.2.1.+ Developed by Dominik Schrempf.+ Compiled on April 27, 2020, at 06:31 am, UTC.+ + Usage: slynx translate (-a|--alphabet NAME) INPUT-FILE (-r|--reading-frame INT)+ (-u|--universal-code CODE) Translate from DNA to Protein or DNAX to ProteinX. Available options:+ -h,--help Show this help text+ -V,--version Show version -a,--alphabet NAME Specify alphabet type NAME INPUT-FILE Read sequences from INPUT-FILE -r,--reading-frame INT Reading frame [0|1|2]. -u,--universal-code CODE universal code; one of: Standard, VertebrateMitochondrial. -h,--help Show this help text+ + + ELynx+ -----+ A Haskell library and tool set for computational biology. The goal of ELynx is+ reproducible research. Evolutionary sequences and phylogenetic trees can be+ read, viewed, modified and simulated. The command line with all arguments is+ logged consistently, and automatically. Data integrity is verified using SHA256+ sums so that validation of past analyses is possible without the need to+ recompute the result.+ + slynx Analyze, modify, and simulate evolutionary sequences.+ tlynx Analyze, modify, and simulate phylogenetic trees.+ elynx Validate and redo past analyses.+ + Get help for sub commands:+ slynx examine --help # TLynx@@ -301,10 +469,12 @@ tlynx --help - ELynx Suite version 0.1.0. Developed by Dominik Schrempf. Compiled on January- 30, 2020, at 09:18 am, UTC.+ ELynx Suite version 0.2.1.+ Developed by Dominik Schrempf.+ Compiled on April 27, 2020, at 06:31 am, UTC. - Usage: tlynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME] COMMAND+ Usage: tlynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME] + [-f|--force] COMMAND Compare, examine, and simulate phylogenetic trees. Available options:@@ -314,35 +484,37 @@ Debug (default: Info) -o,--output-file-basename NAME Specify base name of output file+ -f,--force Ignore previous analysis and overwrite existing+ output files. Available commands:+ coalesce Simulate phylogenetic trees using the coalescent processes (see also the 'simulate' command for simulations using the birth and death process).+ compare Compare two phylogenetic trees (compute distances and branch-wise differences).+ connect Connect two phylogenetic trees in all ways (possibly honoring constraints). distance Compute distances between many phylogenetic trees. examine Compute summary statistics of phylogenetic trees.- simulate Simulate phylogenetic trees using birth and death- processes.- compare Compare two phylogenetic trees (compute distances and- branch-wise differences).- connect Connect two phylogenetic trees in all ways (possibly- honoring constraints).- shuffle Shuffle a phylogenetic tree (keep coalescent times,- but shuffle topology and leaves).+ shuffle Shuffle a phylogenetic tree (keep coalescent times, but shuffle topology and leaves).+ simulate Simulate phylogenetic trees using birth and death processes (see also the 'coalesce' command for simulations using the coalescent process). - File formats:- - Newick+ Available tree file formats:+ - Newick Standard: Branch support values are stored in square brackets after branch lengths.+ - Newick IqTree: Branch support values are stored as node names after the closing bracket of forests.+ - Newick RevBayes Key-value pairs is provided in square brackets after node names as well as branch lengths. XXX: Key value pairs are IGNORED at the moment. - The ELynx Suite- ---------------- A Haskell library and a tool set for computational biology. The goal of the- ELynx Suite is reproducible research. Evolutionary sequences and phylogenetic- trees can be read, viewed, modified and simulated. Exact specification of all- options is necessary, and nothing is assumed about the data (e.g., the type of- code). The command line with all arguments is consistently, and automatically- logged.+ ELynx+ -----+ A Haskell library and tool set for computational biology. The goal of ELynx is+ reproducible research. Evolutionary sequences and phylogenetic trees can be+ read, viewed, modified and simulated. The command line with all arguments is+ logged consistently, and automatically. Data integrity is verified using SHA256+ sums so that validation of past analyses is possible without the need to+ recompute the result. slynx Analyze, modify, and simulate evolutionary sequences. tlynx Analyze, modify, and simulate phylogenetic trees.+ elynx Validate and redo past analyses. - Get help for specific commands:+ Get help for sub commands: slynx examine --help @@ -352,18 +524,43 @@ tlynx compare --help - Usage: tlynx compare [-n|--normalize] [-b|--bipartitions] [-i|--newick-iqtree]- NAME- Compare two phylogenetic trees (compute distances and branch-wise- differences).+ ELynx Suite version 0.2.1.+ Developed by Dominik Schrempf.+ Compiled on April 27, 2020, at 06:31 am, UTC. + Usage: tlynx compare [-n|--normalize] [-b|--bipartitions] [-t|--intersect] + [-f|--newick-format FORMAT] NAME+ Compare two phylogenetic trees (compute distances and branch-wise differences).+ Available options:+ -h,--help Show this help text+ -V,--version Show version -n,--normalize Normalize trees before comparison- -b,--bipartitions Print common and missing bipartitions- -i,--newick-iqtree Use IQ-TREE Newick format (internal node labels are- branch support values)+ -b,--bipartitions Print and plot common and missing bipartitions+ -t,--intersect Compare intersections; i.e., before comparison, drop+ leaves that are not present in the other tree+ -f,--newick-format FORMAT+ Newick tree format; see 'tlynx+ --help' (default: Standard) NAME Tree file -h,--help Show this help text+ + + ELynx+ -----+ A Haskell library and tool set for computational biology. The goal of ELynx is+ reproducible research. Evolutionary sequences and phylogenetic trees can be+ read, viewed, modified and simulated. The command line with all arguments is+ logged consistently, and automatically. Data integrity is verified using SHA256+ sums so that validation of past analyses is possible without the need to+ recompute the result.+ + slynx Analyze, modify, and simulate evolutionary sequences.+ tlynx Analyze, modify, and simulate phylogenetic trees.+ elynx Validate and redo past analyses.+ + Get help for sub commands:+ slynx examine --help ## Examine@@ -372,14 +569,38 @@ tlynx examine --help - Usage: tlynx examine [INPUT-FILE] [-i|--newick-iqtree]+ ELynx Suite version 0.2.1.+ Developed by Dominik Schrempf.+ Compiled on April 27, 2020, at 06:31 am, UTC.+ + Usage: tlynx examine INPUT-FILE [-f|--newick-format FORMAT] Compute summary statistics of phylogenetic trees. Available options:+ -h,--help Show this help text+ -V,--version Show version INPUT-FILE Read trees from INPUT-FILE- -i,--newick-iqtree Use IQ-TREE Newick format (internal node labels are- branch support values)+ -f,--newick-format FORMAT+ Newick tree format; see 'tlynx+ --help' (default: Standard) -h,--help Show this help text+ + + ELynx+ -----+ A Haskell library and tool set for computational biology. The goal of ELynx is+ reproducible research. Evolutionary sequences and phylogenetic trees can be+ read, viewed, modified and simulated. The command line with all arguments is+ logged consistently, and automatically. Data integrity is verified using SHA256+ sums so that validation of past analyses is possible without the need to+ recompute the result.+ + slynx Analyze, modify, and simulate evolutionary sequences.+ tlynx Analyze, modify, and simulate phylogenetic trees.+ elynx Validate and redo past analyses.+ + Get help for sub commands:+ slynx examine --help ## Simulate@@ -388,13 +609,19 @@ tlynx simulate --help - Usage: tlynx simulate [-t|--nTrees INT] [-n|--nLeaves INT] [-H|--height DOUBLE]- [-M|--condition-on-mrca] [-l|--lambda DOUBLE]- [-m|--mu DOUBLE] [-r|--rho DOUBLE] [-u|--sub-sample]+ ELynx Suite version 0.2.1.+ Developed by Dominik Schrempf.+ Compiled on April 27, 2020, at 06:31 am, UTC.+ + Usage: tlynx simulate [-t|--nTrees INT] [-n|--nLeaves INT] [-H|--height DOUBLE] + [-M|--condition-on-mrca] [-l|--lambda DOUBLE] + [-m|--mu DOUBLE] [-r|--rho DOUBLE] [-u|--sub-sample] [-s|--summary-statistics] [-S|--seed [INT]]- Simulate phylogenetic trees using birth and death processes.+ Simulate phylogenetic trees using birth and death processes (see also the 'coalesce' command for simulations using the coalescent process). Available options:+ -h,--help Show this help text+ -V,--version Show version -t,--nTrees INT Number of trees (default: 10) -n,--nLeaves INT Number of leaves per tree (default: 5) -H,--height DOUBLE Fix tree height (no default)@@ -413,15 +640,75 @@ Summary statistics only: only print (NumberOfExtantChildren BranchLength) pairs for each branch of each tree. The trees are separated by a newline character. Sub-sampling: simulate one big tree with n'=round(n/rho), n'>=n, leaves, and randomly sample sub-trees with n leaves. Hence, with rho=1.0, the same tree is reported over and over again. Gernhard, T. (2008). The conditioned reconstructed process. Journal of Theoretical Biology, 253(4), 769–778. http://doi.org/10.1016/j.jtbi.2008.04.005+ + ELynx+ -----+ A Haskell library and tool set for computational biology. The goal of ELynx is+ reproducible research. Evolutionary sequences and phylogenetic trees can be+ read, viewed, modified and simulated. The command line with all arguments is+ logged consistently, and automatically. Data integrity is verified using SHA256+ sums so that validation of past analyses is possible without the need to+ recompute the result.+ + slynx Analyze, modify, and simulate evolutionary sequences.+ tlynx Analyze, modify, and simulate phylogenetic trees.+ elynx Validate and redo past analyses.+ + Get help for sub commands:+ slynx examine --help # ELynx -Documentation of the library can be found on [Hackage](https://hackage.haskell.org/):+Validate and (optionally) redo past ELynx analyses. + elynx --help++ ELynx Suite version 0.2.1.+ Developed by Dominik Schrempf.+ Compiled on April 27, 2020, at 06:31 am, UTC.+ + Usage: elynx COMMAND+ Validate and redo past ELynx analyses+ + Available options:+ -h,--help Show this help text+ -V,--version Show version+ + Available commands:+ validate Validate an ELynx analysis+ redo Redo an ELynx analysis+ + + ELynx+ -----+ A Haskell library and tool set for computational biology. The goal of ELynx is+ reproducible research. Evolutionary sequences and phylogenetic trees can be+ read, viewed, modified and simulated. The command line with all arguments is+ logged consistently, and automatically. Data integrity is verified using SHA256+ sums so that validation of past analyses is possible without the need to+ recompute the result.+ + slynx Analyze, modify, and simulate evolutionary sequences.+ tlynx Analyze, modify, and simulate phylogenetic trees.+ elynx Validate and redo past analyses.+ + Get help for sub commands:+ slynx examine --help+++# Library documentation++Documentation of the libraries can be found on [Hackage](https://hackage.haskell.org/):++- [elynx-markov](https://hackage.haskell.org/package/elynx-markov) - [elynx-seq](https://hackage.haskell.org/package/elynx-seq)-- [elynx-tree](https://hackage.haskell.org/package/elynx-tree) - [elynx-tools](https://hackage.haskell.org/package/elynx-tools)+- [elynx-tree](https://hackage.haskell.org/package/elynx-tree)++Documentation of the executables is also available:++- [elynx](https://hackage.haskell.org/package/elynx) - [slynx](https://hackage.haskell.org/package/slynx) - [tlynx](https://hackage.haskell.org/package/tlynx)
elynx-seq.cabal view
@@ -1,9 +1,9 @@ cabal-version: 1.12 name: elynx-seq-version: 0.1.0+version: 0.2.1 license: GPL-3 license-file: LICENSE-copyright: Dominik Schrempf (2019)+copyright: Dominik Schrempf (2020) maintainer: dominik.schrempf@gmail.com author: Dominik Schrempf homepage: https://github.com/dschrempf/elynx#readme@@ -30,59 +30,35 @@ ELynx.Data.Character.AminoAcidI ELynx.Data.Character.AminoAcidS ELynx.Data.Character.AminoAcidX- ELynx.Data.Character.BoundaryMutation ELynx.Data.Character.Character ELynx.Data.Character.Codon ELynx.Data.Character.Nucleotide ELynx.Data.Character.NucleotideI ELynx.Data.Character.NucleotideX- ELynx.Data.MarkovProcess.AminoAcid- ELynx.Data.MarkovProcess.CXXModels- ELynx.Data.MarkovProcess.CXXModelsData- ELynx.Data.MarkovProcess.EDMModel- ELynx.Data.MarkovProcess.GammaRateHeterogeneity- ELynx.Data.MarkovProcess.MixtureModel- ELynx.Data.MarkovProcess.Nucleotide- ELynx.Data.MarkovProcess.PhyloModel- ELynx.Data.MarkovProcess.RateMatrix- ELynx.Data.MarkovProcess.SubstitutionModel ELynx.Data.Sequence.Alignment ELynx.Data.Sequence.Defaults ELynx.Data.Sequence.Sequence ELynx.Data.Sequence.Translate- ELynx.Export.Sequence.CountsFile ELynx.Export.Sequence.Fasta- ELynx.Import.MarkovProcess.EDMModelPhylobayes- ELynx.Import.MarkovProcess.SiteprofilesPhylobayes ELynx.Import.Sequence.Fasta- ELynx.Simulate.MarkovProcess- ELynx.Simulate.MarkovProcessAlongTree hs-source-dirs: src other-modules: Paths_elynx_seq default-language: Haskell2010- ghc-options: -Wall -O2+ ghc-options: -Wall -fllvm build-depends:- async >=2.2.2 && <2.3,- base >=4.12.0.0 && <4.13,- bytestring >=0.10.8.2 && <0.11,- containers >=0.6.0.1 && <0.7,- data-memocombinators >=0.5.1 && <0.6,- elynx-tools >=0.1.0 && <0.2,- elynx-tree >=0.1.0 && <0.2,- hmatrix >=0.20.0.0 && <0.21,- integration >=0.2.1 && <0.3,- math-functions >=0.3.2.1 && <0.4,+ aeson >=1.4.7.1 && <1.5,+ base >=4.13.0.0 && <4.14,+ bytestring >=0.10.10.0 && <0.11,+ containers >=0.6.2.1 && <0.7,+ elynx-tools >=0.2.1 && <0.3, matrices >=0.5.0 && <0.6,- megaparsec >=7.0.5 && <7.1,+ megaparsec >=8.0.0 && <8.1, mwc-random >=0.14.0.0 && <0.15,- optparse-applicative >=0.14.3.0 && <0.15, parallel >=3.2.2.0 && <3.3,- primitive >=0.6.4.0 && <0.7,- statistics >=0.15.0.0 && <0.16,- transformers >=0.5.6.2 && <0.6,- vector >=0.12.0.3 && <0.13,- vector-th-unbox >=0.2.1.6 && <0.3,+ primitive >=0.7.0.1 && <0.8,+ vector >=0.12.1.2 && <0.13,+ vector-th-unbox >=0.2.1.7 && <0.3, word8 >=0.1.3 && <0.2 test-suite seq-test@@ -91,43 +67,28 @@ hs-source-dirs: test other-modules: ELynx.Data.Alphabet.DistributionDiversitySpec- ELynx.Data.MarkovProcess.AminoAcidSpec- ELynx.Data.MarkovProcess.NucleotideSpec ELynx.Data.Sequence.AlignmentSpec ELynx.Data.Sequence.SequenceSpec ELynx.Data.Sequence.TranslateSpec ELynx.Export.Sequence.FastaSpec- ELynx.Import.MarkovProcess.EDMModelPhylobayesSpec- ELynx.Import.MarkovProcess.SiteprofilesPhylobayesSpec ELynx.Import.Sequence.FastaSpec- ELynx.Simulate.MarkovProcessAlongTreeSpec- Files Paths_elynx_seq default-language: Haskell2010- ghc-options: -Wall -O2 -eventlog -threaded -rtsopts- -with-rtsopts=-N+ ghc-options: -Wall -fllvm -threaded -rtsopts -with-rtsopts=-N build-depends:- async >=2.2.2 && <2.3,- base >=4.12.0.0 && <4.13,- bytestring >=0.10.8.2 && <0.11,- containers >=0.6.0.1 && <0.7,- data-memocombinators >=0.5.1 && <0.6,+ aeson >=1.4.7.1 && <1.5,+ base >=4.13.0.0 && <4.14,+ bytestring >=0.10.10.0 && <0.11,+ containers >=0.6.2.1 && <0.7, elynx-seq -any,- elynx-tools >=0.1.0 && <0.2,- elynx-tree >=0.1.0 && <0.2,- hmatrix >=0.20.0.0 && <0.21,+ elynx-tools >=0.2.1 && <0.3, hspec >=2.7.1 && <2.8,- hspec-megaparsec >=2.0.1 && <2.1,- integration >=0.2.1 && <0.3,- math-functions >=0.3.2.1 && <0.4,+ hspec-megaparsec >=2.1.0 && <2.2, matrices >=0.5.0 && <0.6,- megaparsec >=7.0.5 && <7.1,+ megaparsec >=8.0.0 && <8.1, mwc-random >=0.14.0.0 && <0.15,- optparse-applicative >=0.14.3.0 && <0.15, parallel >=3.2.2.0 && <3.3,- primitive >=0.6.4.0 && <0.7,- statistics >=0.15.0.0 && <0.16,- transformers >=0.5.6.2 && <0.6,- vector >=0.12.0.3 && <0.13,- vector-th-unbox >=0.2.1.6 && <0.3,+ primitive >=0.7.0.1 && <0.8,+ vector >=0.12.1.2 && <0.13,+ vector-th-unbox >=0.2.1.7 && <0.3, word8 >=0.1.3 && <0.2
src/ELynx/Data/Alphabet/Alphabet.hs view
@@ -1,8 +1,10 @@+{-# LANGUAGE DeriveGeneric #-}+ {- | Module : ELynx.Data.Alphabet.Alphabet Description : Alphabets store hereditary information-Copyright : (c) Dominik Schrempf 2019-License : GPL-3+Copyright : (c) Dominik Schrempf 2020+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -30,37 +32,48 @@ -} module ELynx.Data.Alphabet.Alphabet- (- Alphabet (..)- , AlphabetSpec (..)+ ( Alphabet(..)+ , AlphabetSpec(..) , alphabetSpec- , description+ , alphabetDescription , isStd , isGap , isUnknown , isIUPAC , isMember- ) where+ )+where import qualified Data.Set as S-import Prelude hiding (all)+import Prelude hiding ( all )+import Data.Aeson ( FromJSON+ , ToJSON+ )+import GHC.Generics ( Generic ) import ELynx.Data.Alphabet.Character -- | Available alphabets; for details see 'alphabetSpec'. data Alphabet = DNA | DNAX | DNAI | Protein | ProteinX | ProteinS | ProteinI- deriving (Show, Read, Eq, Ord, Enum, Bounded)+ deriving (Show, Read, Eq, Ord, Enum, Bounded, Generic) +instance FromJSON Alphabet++instance ToJSON Alphabet+ -- | Verbose alphabet name.-description :: Alphabet -> String-description DNA = "DNA (nucleotides)"-description DNAX = "DNAX (nucleotides; including gaps)"-description DNAI = "DNAI (nucleotides; including gaps, and IUPAC codes)"-description Protein = "Protein (amino acids)"-description ProteinX = "ProteinX (amino acids; including gaps)"-description ProteinS = "ProteinS (amino acids; including gaps, and translation stops)"-description ProteinI = "ProteinI (amino acids; including gaps, translation stops, and IUPAC codes)"+alphabetDescription :: Alphabet -> String+alphabetDescription DNA = "DNA (nucleotides)"+alphabetDescription DNAX = "DNAX (nucleotides; including gaps)"+alphabetDescription DNAI =+ "DNAI (nucleotides; including gaps, and IUPAC codes)"+alphabetDescription Protein = "Protein (amino acids)"+alphabetDescription ProteinX = "ProteinX (amino acids; including gaps)"+alphabetDescription ProteinS =+ "ProteinS (amino acids; including gaps, and translation stops)"+alphabetDescription ProteinI =+ "ProteinI (amino acids; including gaps, translation stops, and IUPAC codes)" -- | Alphabet specification. 'S.Set' is used because it provides fast lookups. data AlphabetSpec = AlphabetSpec {@@ -111,14 +124,20 @@ isMember :: Alphabet -> Character -> Bool isMember = isWith all -fromChars :: String -> String -> String -> String -> (Char -> String) -> AlphabetSpec-fromChars st ga un iu to = AlphabetSpec st' ga' un' iu' al (fromString . to . toChar)- where- st' = S.fromList $ fromString st- ga' = S.fromList $ fromString ga- un' = S.fromList $ fromString un- iu' = S.fromList $ fromString iu- al = S.unions [st', ga', un', iu']+fromChars+ :: String -> String -> String -> String -> (Char -> String) -> AlphabetSpec+fromChars st ga un iu to = AlphabetSpec st'+ ga'+ un'+ iu'+ al+ (fromString . to . toChar)+ where+ st' = S.fromList $ fromString st+ ga' = S.fromList $ fromString ga+ un' = S.fromList $ fromString un+ iu' = S.fromList $ fromString iu+ al = S.unions [st', ga', un', iu'] dna :: AlphabetSpec dna = fromChars "ACGT" [] [] [] toStdDNA
src/ELynx/Data/Alphabet/Character.hs view
@@ -5,8 +5,8 @@ {- | Module : ELynx.Data.Alphabet.Character Description : Alphabet characters-Copyright : (c) Dominik Schrempf 2019-License : GPL-3+Copyright : (c) Dominik Schrempf 2020+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -17,8 +17,7 @@ -} module ELynx.Data.Alphabet.Character- (- Character+ ( Character , toWord , fromWord , toChar@@ -27,14 +26,16 @@ , fromString , toCVec , fromCVec- ) where+ )+where -import qualified Data.Vector.Unboxed as V+import qualified Data.Vector.Unboxed as V import Data.Vector.Unboxed.Deriving import Data.Word8 -import qualified ELynx.Data.Character.Character as C-import ELynx.Tools.ByteString (c2w, w2c)+import qualified ELynx.Data.Character.Character+ as C+import ELynx.Tools -- | Alphabet characters; abstracted so that representation can be changed at -- some point.
src/ELynx/Data/Alphabet/DistributionDiversity.hs view
@@ -1,8 +1,10 @@+{-# LANGUAGE FlexibleContexts #-}+ {- | Module : ELynx.Data.Alphabet.DistributionDiversity Description : Summarize statistics for alphabets-Copyright : (c) Dominik Schrempf 2019-License : GPL-3+Copyright : (c) Dominik Schrempf 2020+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -21,74 +23,75 @@ , kEffHomoplasy -- * Count characters , frequencyCharacters- ) where+ )+where import qualified Data.Set as S-import qualified Data.Vector.Unboxed as V+import qualified Data.Vector.Generic as V+import Data.Vector.Generic ( Vector+ , toList+ ) import ELynx.Data.Alphabet.Alphabet import ELynx.Data.Alphabet.Character-import ELynx.Tools.Definitions-import ELynx.Tools.Numeric-import ELynx.Tools.Vector+import ELynx.Tools -- | Entropy of vector.-entropy :: V.Vector Double -> Double+entropy :: (Vector v Double) => v Double -> Double entropy v = if isNaN res- then error ("entropy: Sesult of following vector is NaN: " ++ show v ++ ".")+ then error+ ("entropy: Sesult of following vector is NaN: " ++ show (toList v) ++ ".") else res where res = negate $ sumVec $ V.map xLogX v -- | Effective number of used characters measured using 'entropy'. The result -- only makes sense when the sum of the array is 1.0.-kEffEntropy :: V.Vector Double -> Double-kEffEntropy v = if e < eps- then 1.0- else exp e- where e = entropy v+kEffEntropy :: Vector v Double => v Double -> Double+kEffEntropy v = if e < eps then 1.0 else exp e where e = entropy v -- | Probability of homoplasy of vector. The result is the probability of -- binomially sampling the same character twice and only makes sense when the -- sum of the array is 1.0.-homoplasy :: V.Vector Double -> Double-homoplasy v = sumVec $ V.map (\x -> x*x) v+homoplasy :: Vector v Double => v Double -> Double+homoplasy v = sumVec $ V.map (\x -> x * x) v -- | Effective number of used characters measured using 'homoplasy'. The result -- only makes sense when the sum of the array is 1.0.-kEffHomoplasy :: V.Vector Double -> Double+kEffHomoplasy :: Vector v Double => v Double -> Double kEffHomoplasy v = 1.0 / homoplasy v -- XXX: Use mutable vector; then V.// is much faster. -- Increment element at index in vector by one.-incrementElemIndexByOne :: [Int] -> V.Vector Int -> V.Vector Int+incrementElemIndexByOne :: Vector v Int => [Int] -> v Int -> v Int incrementElemIndexByOne is v = v V.// zip is es'- where es' = [v V.! i + 1 | i <- is]+ where es' = [ v V.! i + 1 | i <- is ] -- For a given code and counts vector, increment the count of the given character.-acc :: AlphabetSpec -> V.Vector Int -> Character -> V.Vector Int+acc :: Vector v Int => AlphabetSpec -> v Int -> Character -> v Int acc alph vec char = incrementElemIndexByOne is vec- where- is = [ S.findIndex c (std alph) | c <- toStd alph char ]+ where is = [ S.findIndex c (std alph) | c <- toStd alph char ] -countCharacters :: AlphabetSpec -> V.Vector Character -> V.Vector Int-countCharacters alph =- V.foldl' (acc alph) zeroCounts- where- nChars = length (std alph)- zeroCounts = V.replicate nChars (0 :: Int)+countCharacters+ :: (Vector v Character, Vector v Int) => AlphabetSpec -> v Character -> v Int+countCharacters alph = V.foldl' (acc alph) zeroCounts+ where+ nChars = length (std alph)+ zeroCounts = V.replicate nChars (0 :: Int) saveDivision :: Int -> Int -> Double saveDivision value divisor =- if divisor == 0- then 0.0- else fromIntegral value / fromIntegral divisor+ if divisor == 0 then 0.0 else fromIntegral value / fromIntegral divisor -- | For a given code vector of characters, calculate frequency of characters. -- The input vector has arbitrary length (most often the number of sequences in -- an alignment), the length of the output vector is the number of characters in -- the alphabet.-frequencyCharacters :: AlphabetSpec -> V.Vector Character -> V.Vector Double+frequencyCharacters+ :: (Vector v Character, Vector v Int, Vector v Double)+ => AlphabetSpec+ -> v Character+ -> v Double frequencyCharacters alph d = V.map (`saveDivision` s) counts- where- counts = countCharacters alph d- s = sumVec counts+ where+ counts = countCharacters alph d+ s = sumVec counts
src/ELynx/Data/Character/AminoAcid.hs view
@@ -6,7 +6,7 @@ Module : ELynx.Data.AminoAcid Description : Amino acid related types and functions Copyright : (c) Dominik Schrempf 2018-License : GPL-3+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -46,14 +46,16 @@ -} module ELynx.Data.Character.AminoAcid- ( AminoAcid (..)- ) where+ ( AminoAcid(..)+ )+where import Data.Vector.Unboxed.Deriving import Data.Word8 -import qualified ELynx.Data.Character.Character as C-import ELynx.Tools.ByteString (c2w, w2c)+import qualified ELynx.Data.Character.Character+ as C+import ELynx.Tools -- | Amino acids. data AminoAcid = A | C | D | E | F | G | H | I | K | L | M | N | P | Q | R | S | T | V | W | Y@@ -83,27 +85,27 @@ fromWord :: Word8 -> AminoAcid fromWord w = case w2c w of- 'A' -> A- 'C' -> C- 'D' -> D- 'E' -> E- 'F' -> F- 'G' -> G- 'H' -> H- 'I' -> I- 'K' -> K- 'L' -> L- 'M' -> M- 'N' -> N- 'P' -> P- 'Q' -> Q- 'R' -> R- 'S' -> S- 'T' -> T- 'V' -> V- 'W' -> W- 'Y' -> Y- _ -> error "fromWord: Cannot convert to AminoAcid."+ 'A' -> A+ 'C' -> C+ 'D' -> D+ 'E' -> E+ 'F' -> F+ 'G' -> G+ 'H' -> H+ 'I' -> I+ 'K' -> K+ 'L' -> L+ 'M' -> M+ 'N' -> N+ 'P' -> P+ 'Q' -> Q+ 'R' -> R+ 'S' -> S+ 'T' -> T+ 'V' -> V+ 'W' -> W+ 'Y' -> Y+ _ -> error "fromWord: Cannot convert to AminoAcid." derivingUnbox "AminoAcid" [t| AminoAcid -> Word8 |]
src/ELynx/Data/Character/AminoAcidI.hs view
@@ -6,7 +6,7 @@ Module : ELynx.Data.AminoAcid Description : Amino acid related types and functions Copyright : (c) Dominik Schrempf 2018-License : GPL-3+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -58,14 +58,16 @@ -} module ELynx.Data.Character.AminoAcidI- ( AminoAcidI (..)- ) where+ ( AminoAcidI(..)+ )+where import Data.Vector.Unboxed.Deriving import Data.Word8 -import qualified ELynx.Data.Character.Character as C-import ELynx.Tools.ByteString (c2w, w2c)+import qualified ELynx.Data.Character.Character+ as C+import ELynx.Tools -- | Amino acids. data AminoAcidI = A | C | D | E | F | G | H | I | K | L | M | N | P | Q | R | S | T | V | W | Y@@ -105,34 +107,34 @@ fromWord :: Word8 -> AminoAcidI fromWord w = case w2c w of- 'A' -> A- 'C' -> C- 'D' -> D- 'E' -> E- 'F' -> F- 'G' -> G- 'H' -> H- 'I' -> I- 'K' -> K- 'L' -> L- 'M' -> M- 'N' -> N- 'P' -> P- 'Q' -> Q- 'R' -> R- 'S' -> S- 'T' -> T- 'V' -> V- 'W' -> W- 'Y' -> Y- 'J' -> J- 'B' -> B- 'Z' -> Z- 'X' -> X- '*' -> Stop- '-' -> Gap- '.' -> Gap- _ -> error "fromWord: Cannot convert Word8 to AminoAcidI"+ 'A' -> A+ 'C' -> C+ 'D' -> D+ 'E' -> E+ 'F' -> F+ 'G' -> G+ 'H' -> H+ 'I' -> I+ 'K' -> K+ 'L' -> L+ 'M' -> M+ 'N' -> N+ 'P' -> P+ 'Q' -> Q+ 'R' -> R+ 'S' -> S+ 'T' -> T+ 'V' -> V+ 'W' -> W+ 'Y' -> Y+ 'J' -> J+ 'B' -> B+ 'Z' -> Z+ 'X' -> X+ '*' -> Stop+ '-' -> Gap+ '.' -> Gap+ _ -> error "fromWord: Cannot convert Word8 to AminoAcidI" derivingUnbox "AminoAcidI" [t| AminoAcidI -> Word8 |]
src/ELynx/Data/Character/AminoAcidS.hs view
@@ -6,7 +6,7 @@ Module : ELynx.Data.AminoAcid Description : Amino acid related types and functions Copyright : (c) Dominik Schrempf 2018-License : GPL-3+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -51,14 +51,16 @@ -} module ELynx.Data.Character.AminoAcidS- ( AminoAcidS (..)- ) where+ ( AminoAcidS(..)+ )+where import Data.Vector.Unboxed.Deriving import Data.Word8 -import qualified ELynx.Data.Character.Character as C-import ELynx.Tools.ByteString (c2w, w2c)+import qualified ELynx.Data.Character.Character+ as C+import ELynx.Tools -- | Amino acids. data AminoAcidS = A | C | D | E | F | G | H | I | K | L | M | N | P | Q | R | S | T | V | W | Y@@ -92,30 +94,30 @@ fromWord :: Word8 -> AminoAcidS fromWord w = case w2c w of- 'A' -> A- 'C' -> C- 'D' -> D- 'E' -> E- 'F' -> F- 'G' -> G- 'H' -> H- 'I' -> I- 'K' -> K- 'L' -> L- 'M' -> M- 'N' -> N- 'P' -> P- 'Q' -> Q- 'R' -> R- 'S' -> S- 'T' -> T- 'V' -> V- 'W' -> W- 'Y' -> Y- '*' -> Stop- '-' -> Gap- '.' -> Gap- _ -> error "fromWord: Cannot convert to AminoAcidS."+ 'A' -> A+ 'C' -> C+ 'D' -> D+ 'E' -> E+ 'F' -> F+ 'G' -> G+ 'H' -> H+ 'I' -> I+ 'K' -> K+ 'L' -> L+ 'M' -> M+ 'N' -> N+ 'P' -> P+ 'Q' -> Q+ 'R' -> R+ 'S' -> S+ 'T' -> T+ 'V' -> V+ 'W' -> W+ 'Y' -> Y+ '*' -> Stop+ '-' -> Gap+ '.' -> Gap+ _ -> error "fromWord: Cannot convert to AminoAcidS." derivingUnbox "AminoAcidS" [t| AminoAcidS -> Word8 |]
src/ELynx/Data/Character/AminoAcidX.hs view
@@ -6,7 +6,7 @@ Module : ELynx.Data.AminoAcid Description : Amino acid related types and functions Copyright : (c) Dominik Schrempf 2018-License : GPL-3+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -51,14 +51,16 @@ -} module ELynx.Data.Character.AminoAcidX- ( AminoAcidX (..)- ) where+ ( AminoAcidX(..)+ )+where import Data.Vector.Unboxed.Deriving import Data.Word8 -import qualified ELynx.Data.Character.Character as C-import ELynx.Tools.ByteString (c2w, w2c)+import qualified ELynx.Data.Character.Character+ as C+import ELynx.Tools -- | Amino acids. data AminoAcidX = A | C | D | E | F | G | H | I | K | L | M | N | P | Q | R | S | T | V | W | Y@@ -90,29 +92,29 @@ fromWord :: Word8 -> AminoAcidX fromWord w = case w2c w of- 'A' -> A- 'C' -> C- 'D' -> D- 'E' -> E- 'F' -> F- 'G' -> G- 'H' -> H- 'I' -> I- 'K' -> K- 'L' -> L- 'M' -> M- 'N' -> N- 'P' -> P- 'Q' -> Q- 'R' -> R- 'S' -> S- 'T' -> T- 'V' -> V- 'W' -> W- 'Y' -> Y- '-' -> Gap- '.' -> Gap- _ -> error "fromWord: Cannot convert to AminoAcidX."+ 'A' -> A+ 'C' -> C+ 'D' -> D+ 'E' -> E+ 'F' -> F+ 'G' -> G+ 'H' -> H+ 'I' -> I+ 'K' -> K+ 'L' -> L+ 'M' -> M+ 'N' -> N+ 'P' -> P+ 'Q' -> Q+ 'R' -> R+ 'S' -> S+ 'T' -> T+ 'V' -> V+ 'W' -> W+ 'Y' -> Y+ '-' -> Gap+ '.' -> Gap+ _ -> error "fromWord: Cannot convert to AminoAcidX." derivingUnbox "AminoAcidX" [t| AminoAcidX -> Word8 |]@@ -124,4 +126,4 @@ fromWord = fromWord instance C.CharacterX AminoAcidX where- gap = Gap+ gap = Gap
− src/ELynx/Data/Character/BoundaryMutation.hs
@@ -1,281 +0,0 @@-{- |-Description : State space of the boundary mutation model-Copyright : (c) Dominik Schrempf 2017-License : GPLv3--Maintainer : dominik.schrempf@gmail.com-Stability : unstable-Portability : non-portable (not tested)--The boundary mutation model is a discrete-state, continuous-time Markov process-that allows mutations only when the population is monomorphic.--* Changelog--TODO: BM states can not be read and written with a single letter, like-characters. Check status of implementation.---}--module ELynx.Data.Character.BoundaryMutation- ( -- * Types- Nucleotide- , Allele- , PopulationSize- , AlleleCount- , State(..)- , showCounts- , nFixed- -- * Functions- , setPopulationSize- , fromIndexWith- , toIndex- , stateSpace- , stateSpaceSize- , neighbors- ) where--import qualified Data.ByteString.Lazy.Char8 as L-import Numeric.SpecFunctions (choose)--import ELynx.Data.Character.Nucleotide-import ELynx.Tools.Misc---- | Alleles are just nucleotides at the moment. However, I want to keep the--- code such that it can be extended easily to codons or amino acids.-type Allele = Nucleotide--- | The population size has to be larger than one otherwise there be dragons.-type PopulationSize = Int--- | The absolute frequency of an allele.-type AlleleCount = Int---- | The number of alleles.-nAlleles :: Int-nAlleles = 1 + fromEnum (maxBound :: Allele)---- | A boundary mutation model state is either a boundary state or a polymorphic--- state. The population size has to be larger than one; the allele count has to--- be larger than one and lower than the population size, otherwise there be--- dragons.------ Another possibility would be:--- @--- data State = Bnd Allele | Ply AlleleCount Allele Allele--- data StateComplete = StateComplete PopulationSize State--- @--- But then, I think it is more important that the information is kept in one,--- at the cost of some overhead.-data State = Bnd { bndN :: PopulationSize -- | Population size.- , bndA :: Allele }- | Ply { plyN :: PopulationSize -- | Population size.- , plyI :: AlleleCount -- | Allele count.- , plyA :: Allele- , plyB :: Allele }- deriving (Read, Eq)---- | L.ByteString representation of 'State'; without surrounding brackets.-showCounts :: State -> L.ByteString-showCounts (Bnd n a) = L.intersperse ',' $ L.concat $ map (L.pack . toCounts) allValues- where toCounts b- | a == b = show n- | otherwise = "0"-showCounts (Ply n i a b) = L.intersperse ',' $ L.concat $ map (L.pack . toCounts) allValues- where toCounts c- | c == a = show i- | c == b = show (n-i)- | otherwise = "0"--showState :: State -> L.ByteString-showState s = L.singleton '(' <> showCounts s <> L.singleton ')'---- instance Show State where--- show s = "(" ++ showCounts s ++ ")"---- | A total order on the boundary mutation model states. In general, Bnd < Ply.--- Then, sorting happens according to the order population size, first allele,--- second allele, allele count. It may be beneficial to reverse the allele count--- order (i.e., make a polymorphic state with higher allele count show up before--- a polymorphic state with lower allele count, this would move some polymorphic--- states closer to their respective boundaries),-instance Ord State where- Bnd {} <= Ply {} = True- Ply {} <= Bnd {} = False- s@(Bnd n a) <= t@(Bnd m b)- | s == t = True- | n /= m = n <= m- | otherwise = a <= b- s@(Ply n i a b) <= t@(Ply m j c d)- | s == t = True- | n /= m = n <= m- | a < c = True- | a > c = False- -- We can be sure that a = c now.- | b < d = True- | b > d = False- -- Now we can be sure that both nucleotides are the same.- | otherwise = i <= j---- | Fixed population size when converting a 'State' to or from a number. In--- this case, a fixed population size is necessary so that @toEnum . fromEnum ==--- id@. When converting from a number to 'State', the population size has to be--- given or assumed (see 'fromIndexWith') anyways. Especially when performing IO,--- the same number should always correspond to the same 'State' (bijection).--- 'nFixed' has been set such that the size of the state space is 256.-nFixed :: Int-nFixed = 43---- | Set the population size of a 'State'; validity of resulting 'State' is checked.-setPopulationSize :: PopulationSize -> State -> Maybe State-setPopulationSize n s = if valid s' then Just s' else Nothing- where s' = unsafeSetPopulationSize n s---- | See 'setPopulationSize'. Does not check if resulting 'State' is valid.-unsafeSetPopulationSize :: Int -> State -> State-unsafeSetPopulationSize n (Bnd _ s) = Bnd n s-unsafeSetPopulationSize n (Ply _ i a b) = Ply n i a b---- | For a given population size 'PopulationSize', convert a number 'Int' to 'State'.-fromIndexWith :: PopulationSize -> Int -> State-fromIndexWith n i- | i >= stateSpaceSize n = error $- "Index " ++ show i ++ "out of bounds when population size is " ++ show n ++ "."- | i < nAlleles = Bnd n (toEnum i)- | otherwise = Ply n (i' - j + 1) x k- where i' = i - nAlleles- l = [ (enumCombination a b * (n-1), a, b)- | a <- [minBound .. pred maxBound]- , b <- [succ a ..]]- (j, x, k) = last $ takeWhile (\(e, _, _) -> e <= i') l---- | Convert 'State' to a number 'Int' for the given population size 'PopulationSize'.--- Back conversion can be done with 'fromIndexWith', with the same population size.-toIndex :: State -> Int-toIndex (Bnd _ a) = fromEnum a--- We also have to shift the enumeration value by the number of boundary--- states, which is 'nAlleles'.-toIndex (Ply n i a b) = nAlleles + enumCombination a b * (n-1) + i-1---- | Enumeration only works when the population size is 'nFixed'. Only then,--- @toEnum . fromEnum == id@ can be guaranteed. This is because @toEnum ::--- State@ is only defined if the population size is known. See also--- 'fromIndexWith', and 'toIndex', as well as, 'setPopulationSize'.-instance Enum State where- fromEnum s = if getPopulationSize s /= nFixed- then error $ "State is not enumerable: " ++ (L.unpack . showState) s ++ "."- else toIndex s- toEnum = fromIndexWith nFixed---- The formula is a little complicated. Sketch of derivation: Order the states--- in the following way:--- @--- AC CG GT--- AG CT--- AT--- @--- The edge length of the triangle is @'nAlleles' - 1@. Use Gauss's triangle--- equation @area=binom(length+1, 2)@ twice to count the number of combinations--- up to a certain allele. E.g., up to, but excluding G:--- @--- AC CG--- AG CT--- AT--- @-countCombinationsUpToAllele :: Allele -> Int-countCombinationsUpToAllele a = round $ nAlleles `choose` 2 - (nAlleles - fromEnum a) `choose` 2---- See 'countCombinationsUpToAllele'. The @-1@ pops up because we start counting--- from 0. For example, the enumeration value of @GT@ (with @fromEnum G = k = 2@--- and @fromEnum T = 3@) is then @enumCombinationsUpToK 2 + (3-2)@.-enumCombination :: Allele -> Allele -> Int-enumCombination a b = countCombinationsUpToAllele a - 1 + (fromEnum b - fromEnum a)---- | A fixed population size 'nFixed' is assumed.-instance Bounded State where- minBound = Bnd nFixed minBound- maxBound = Ply nFixed (nFixed-1) (pred maxBound) maxBound---- -- I am not sure if I should remove the 'Character' instance because writing--- -- Fasta files with boundary mutation model states is not really promising--- -- anyways. However, the 'toIndex' and 'fromIndexWith' function provide a--- -- convenient way to map states to integers. This functionality is needed when--- -- working with matrices.--- -- | A fixed population size 'nFixed' is assumed.--- instance Character State where--- fromWord = toEnum . fromEnum--- toWord = toEnum . fromEnum--- -- FIXME: This requires more thought. Are polymorphic characters standard?--- isStandard _ = error "Requires more thought."--- -- FIXME: This requires more work. (0,0,0,0) should be a gap!--- isGapOrUnknown _ = error "Not implemented."--valid :: State -> Bool-valid (Bnd n _)- | n <= 1 = False- | otherwise = True-valid (Ply n i a b)- | n <= 1 = False- | a >= b = False- | i <= 0 = False- | i >= n = False- | otherwise = True--filterValidStates :: [State] -> [State]-filterValidStates = filter valid--getPopulationSize :: State -> PopulationSize-getPopulationSize (Bnd n _) = n-getPopulationSize (Ply n _ _ _) = n---- CCC: This function is a not very efficient. A better would be something like:--- @--- | Sorted list of all possible PoMo states for a specific population size.-stateSpace :: PopulationSize -> [State]-stateSpace n = map (fromIndexWith n) [0 .. stateSpaceSize n - 1]---- An easier, but slower implementation.--- stateSpace n--- | n <= 1 = error "The population size has to be larger than one."--- | otherwise = sort $ filterValidStates ( allBndStates ++ allPlyStates )--- where allBndStates = [ Bnd n a |--- a <- [minBound .. maxBound] :: [Allele] ]--- allPlyStates = [ Ply n i a b |--- i <- [0..n],--- a <- [minBound .. maxBound] :: [Allele],--- b <- [minBound .. maxBound] :: [Allele] ]---- | The state space of the boundary mutation model for four alleles and a--- population size N is 4 + 6*(N-1).-stateSpaceSize :: PopulationSize -> Int-stateSpaceSize n = k + k*(k-1) `div` 2 * (n-1)- where k = nAlleles---- -- This is a very convenient version of toIndex, but can we always guarantee--- -- that the state space is sorted the same way?--- -- | Convert a boundary state to its ID (integer). See also 'idToState'.--- stateToId :: State -> Maybe Int--- stateToId s = elemIndex s (stateSpace $ getPopulationSize s)---- -- Same here.--- -- | Convert an ID to a boundary state. See also 'stateID'.--- idToState :: PopulationSize -> Int -> State--- idToState n i = stateSpace n !! i---- | Check if two states are connected. By definition, states are NOT connected--- with themselves.-neighbors :: State -> State -> Bool-neighbors s t = s `elem` getNeighbors t--getNeighbors :: State -> [State]-getNeighbors (Bnd n a) = filterValidStates allNeighbors- where allNeighbors = [ Ply n (n-1) a b |- b <- [minBound .. maxBound] :: [Allele] ]- ++- [ Ply n 1 b a |- b <- [minBound .. maxBound] :: [Allele] ]-getNeighbors (Ply n i a b)- -- Careful when the population size is two, because then each polymorphic- -- states has two boundary states as neighbors.- | i == 1 && n == 2 = Bnd n a : [Bnd n b]- | i == 1 = Bnd n b : [Ply n 2 a b]- | i == (n-1) = Bnd n a : [Ply n (n-2) a b]- | otherwise = Ply n (i+1) a b : [Ply n (i-1) a b]
src/ELynx/Data/Character/Character.hs view
@@ -2,7 +2,7 @@ Module : Character Description : Character interface Copyright : (c) Dominik Schrempf 2018-License : GPL-3+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -15,27 +15,26 @@ -} module ELynx.Data.Character.Character- (- Character (..)+ ( Character(..) , fromChar , toChar , fromString , toString- , CharacterX (..)+ , CharacterX(..) , isGap- , CharacterI (..)+ , CharacterI(..) , isUnknown , isIUPAC , isStandard , convert- ) where--import qualified Data.Set as S-import Data.Vector.Unboxed.Base (Unbox)-import Data.Word8 (Word8)+ )+where -import ELynx.Tools.ByteString (c2w, w2c)+import qualified Data.Set as S+import Data.Vector.Unboxed.Base ( Unbox )+import Data.Word8 ( Word8 ) +import ELynx.Tools -- XXX: Remove name clash with ELynx.Data.Alphabet.Alphabet.Character? -- | A set of characters forms an 'ELynx.Data.Alphabet.Alphabet'. At the
src/ELynx/Data/Character/Codon.hs view
@@ -1,8 +1,10 @@+{-# LANGUAGE DeriveGeneric #-}+ {- | Module : ELynx.Data.Character.Codon Description : Codons are triplets of nucleotides-Copyright : (c) Dominik Schrempf 2019-License : GPL-3+Copyright : (c) Dominik Schrempf 2020+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -18,25 +20,34 @@ -} module ELynx.Data.Character.Codon- ( Codon (Codon)- , unsafeFromList+ ( Codon(Codon) , unsafeFromVec- , UniversalCode (..)+ , UniversalCode(..) , translate , translateX , translateI- ) where+ )+where +import Data.Aeson ( FromJSON+ , ToJSON+ ) import Data.List-import qualified Data.Map as M-import qualified Data.Vector.Generic as V+import GHC.Generics ( Generic )+import qualified Data.Map as M+import qualified Data.Vector.Generic as V -import qualified ELynx.Data.Character.AminoAcidI as AI+import qualified ELynx.Data.Character.AminoAcidI+ as AI import ELynx.Data.Character.AminoAcidS-import qualified ELynx.Data.Character.Character as C-import qualified ELynx.Data.Character.Nucleotide as N-import qualified ELynx.Data.Character.NucleotideI as NI-import qualified ELynx.Data.Character.NucleotideX as NX+import qualified ELynx.Data.Character.Character+ as C+import qualified ELynx.Data.Character.Nucleotide+ as N+import qualified ELynx.Data.Character.NucleotideI+ as NI+import qualified ELynx.Data.Character.NucleotideX+ as NX -- | Codons are triplets of characters. newtype Codon a = Codon (a, a, a)@@ -45,56 +56,175 @@ convert :: (C.Character a, C.Character b) => Codon a -> Codon b convert (Codon (x, y, z)) = Codon (C.convert x, C.convert y, C.convert z) --- | Unsafe conversion from list with three elements.-unsafeFromList :: [a] -> Codon a-unsafeFromList xs = Codon (head xs, head . tail $ xs, head . tail . tail $ xs)---- | Unsafe conversion from vector with three elements.+-- | Unsafe conversion from vector with at least three elements; only the first+-- three elements are used, the rest is discarded. unsafeFromVec :: V.Vector v a => v a -> Codon a-unsafeFromVec xs = Codon (V.head xs, V.head . V.tail $ xs, V.head . V.tail . V.tail $ xs)+unsafeFromVec xs =+ Codon (V.head xs, V.head . V.tail $ xs, V.head . V.tail . V.tail $ xs) -- | Universal codes. data UniversalCode = Standard | VertebrateMitochondrial- deriving (Show, Read, Eq, Ord, Enum, Bounded)+ deriving (Show, Read, Eq, Ord, Enum, Bounded, Generic) +instance FromJSON UniversalCode++instance ToJSON UniversalCode+ -- It is important that the map is lazy, because some keys have errors as values.-mapFromLists :: Ord a => [a] -> [a] -> [a]- -> [b] -> M.Map (Codon a) b-mapFromLists xs ys zs as = M.fromList $- zipWith4 (\f s t a -> (Codon (f, s, t), a)) xs ys zs as+mapFromLists :: Ord a => [a] -> [a] -> [a] -> [b] -> M.Map (Codon a) b+mapFromLists xs ys zs as =+ M.fromList $ zipWith4 (\f s t a -> (Codon (f, s, t), a)) xs ys zs as nucs :: Enum a => [a]-nucs = map toEnum [3,1,0,2] -- Order T, C, A , G.+nucs = map toEnum [3, 1, 0, 2] -- Order T, C, A , G. -- Permutation of the triplets PLUS GAPS! I avoid 'Z' because I do not want to -- translate DNAI. base1, base2, base3 :: Enum a => [a]-base1 = [n | n <- nucs- , _ <- [0..3 :: Int]- , _ <- [0..3 :: Int]]+base1 = [ n | n <- nucs, _ <- [0 .. 3 :: Int], _ <- [0 .. 3 :: Int] ] -- base1 = "TTTTTTTTTTTTTTTTCCCCCCCCCCCCCCCCAAAAAAAAAAAAAAAAGGGGGGGGGGGGGGGG" ++ "-."-base2 = [n | _ <- [0..3 :: Int]- , n <- nucs- , _ <- [0..3 :: Int]]+base2 = [ n | _ <- [0 .. 3 :: Int], n <- nucs, _ <- [0 .. 3 :: Int] ] -- base2 = "TTTTCCCCAAAAGGGGTTTTCCCCAAAAGGGGTTTTCCCCAAAAGGGGTTTTCCCCAAAAGGGG" ++ "-."-base3 = [n | _ <- [0..3 :: Int]- , _ <- [0..3 :: Int]- , n <- nucs]+base3 = [ n | _ <- [0 .. 3 :: Int], _ <- [0 .. 3 :: Int], n <- nucs ] -- base3 = "TCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAG" ++ "-." -- The actual codes. standard :: [AminoAcidS]-standard = [ F, F, L, L, S, S, S, S, Y, Y, Stop, Stop, C, C, Stop, W, L, L, L,- L, P, P, P, P, H, H, Q, Q, R, R, R, R, I, I, I, M, T, T, T, T, N,- N, K, K, S, S, R, R, V, V, V, V, A, A, A, A, D, D, E, E, G, G, G,- G]+standard =+ [ F+ , F+ , L+ , L+ , S+ , S+ , S+ , S+ , Y+ , Y+ , Stop+ , Stop+ , C+ , C+ , Stop+ , W+ , L+ , L+ , L+ , L+ , P+ , P+ , P+ , P+ , H+ , H+ , Q+ , Q+ , R+ , R+ , R+ , R+ , I+ , I+ , I+ , M+ , T+ , T+ , T+ , T+ , N+ , N+ , K+ , K+ , S+ , S+ , R+ , R+ , V+ , V+ , V+ , V+ , A+ , A+ , A+ , A+ , D+ , D+ , E+ , E+ , G+ , G+ , G+ , G+ ] -- "FFLLSSSSYY**CC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG" ++ "--" vertebrateMitochondrial :: [AminoAcidS]-vertebrateMitochondrial = [F, F, L, L, S, S, S, S, Y, Y, Stop, Stop, C, C, W, W,- L, L, L, L, P, P, P, P, H, H, Q, Q, R, R, R, R, I, I,- M, M, T, T, T, T, N, N, K, K, S, S, Stop, Stop, V, V,- V, V, A, A, A, A, D, D, E, E, G, G, G, G]+vertebrateMitochondrial =+ [ F+ , F+ , L+ , L+ , S+ , S+ , S+ , S+ , Y+ , Y+ , Stop+ , Stop+ , C+ , C+ , W+ , W+ , L+ , L+ , L+ , L+ , P+ , P+ , P+ , P+ , H+ , H+ , Q+ , Q+ , R+ , R+ , R+ , R+ , I+ , I+ , M+ , M+ , T+ , T+ , T+ , T+ , N+ , N+ , K+ , K+ , S+ , S+ , Stop+ , Stop+ , V+ , V+ , V+ , V+ , A+ , A+ , A+ , A+ , D+ , D+ , E+ , E+ , G+ , G+ , G+ , G+ ] -- "FFLLSSSSYY**CCWWLLLLPPPPHHQQRRRRIIMMTTTTNNKKSS**VVVVAAAADDEEGGGG" ++ "--" -- | Translate a codon to amino acids including translation stops.@@ -102,29 +232,31 @@ translate code = (M.!) (universalCode code) -- | Translate a codon to amino acids including translation stops. Translate--- codons including gaps to amino acid gaps. XXX: Be careful, single or two--- character gaps could have led to a reading frame shift and hence, the--- translated sequence may be bogus.+-- codons including gaps to amino acid gaps. Be careful, single or two character+-- gaps can lead to a reading frame shift and hence, the translated sequence may+-- be bogus. translateX :: UniversalCode -> Codon NX.NucleotideX -> AminoAcidS -- translateX _ (Codon (NX.Gap, NX.Gap, NX.Gap)) = Gap -- translateX code codon = C.convert . translate code . convert $ codon-translateX code codon@(Codon (x,y,z)) | C.isGap x || C.isGap y || C.isGap z = Gap- | otherwise = C.convert . translate code . convert $ codon+translateX code codon@(Codon (x, y, z))+ | C.isGap x || C.isGap y || C.isGap z = Gap+ | otherwise = C.convert . translate code . convert $ codon -- | Translate a codon to amino acids including translation stops. Translate gap -- triplets to amino acid gaps, and triplets including unknowns to amino acid--- unknowns. XXX: Be careful, also translates other IUPAC characters at the--- moment (to amino acid Xs)!+-- unknowns. Be careful, also translates other IUPAC characters to amino acid Xs!+translateI :: UniversalCode -> Codon NI.NucleotideI -> AI.AminoAcidI+translateI code codon@(Codon (x, y, z))+ | C.isIUPAC x || C.isIUPAC y || C.isIUPAC z = AI.X+ | otherwise = C.convert . translateX code . convert $ codon -- translateI :: UniversalCode -> Codon NI.NucleotideI -> AI.AminoAcidI -- translateI _ (Codon (NI.N, _, _ )) = AI.X -- translateI _ (Codon (_ , NI.N, _ )) = AI.X -- translateI _ (Codon (_, _, NI.N)) = AI.X -- translateI code codon = C.convert . translateX code . convert $ codon-translateI :: UniversalCode -> Codon NI.NucleotideI -> AI.AminoAcidI-translateI code codon@(Codon (x,y,z)) | C.isIUPAC x || C.isIUPAC y || C.isIUPAC z = AI.X- | otherwise = C.convert . translateX code . convert $ codon -- Map from 'Codon' to amino acid character. universalCode :: UniversalCode -> M.Map (Codon N.Nucleotide) AminoAcidS-universalCode Standard = mapFromLists base1 base2 base3 standard-universalCode VertebrateMitochondrial = mapFromLists base1 base2 base3 vertebrateMitochondrial+universalCode Standard = mapFromLists base1 base2 base3 standard+universalCode VertebrateMitochondrial =+ mapFromLists base1 base2 base3 vertebrateMitochondrial
src/ELynx/Data/Character/Nucleotide.hs view
@@ -6,7 +6,7 @@ Module : ELynx.Data.Nucleotide Description : Nucleotides Copyright : (c) Dominik Schrempf 2018-License : GPL-3+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -28,14 +28,16 @@ -} module ELynx.Data.Character.Nucleotide- ( Nucleotide (..)- ) where+ ( Nucleotide(..)+ )+where import Data.Vector.Unboxed.Deriving import Data.Word8 -import qualified ELynx.Data.Character.Character as C-import ELynx.Tools.ByteString (c2w, w2c)+import qualified ELynx.Data.Character.Character+ as C+import ELynx.Tools -- | Nucleotides. data Nucleotide = A | C | G | T@@ -53,11 +55,11 @@ fromWord :: Word8 -> Nucleotide fromWord w = case w2c w of- 'A' -> A- 'C' -> C- 'G' -> G- 'T' -> T- c -> error $ "fromWord: Cannot convert " ++ show c ++ " to Nucleotide."+ 'A' -> A+ 'C' -> C+ 'G' -> G+ 'T' -> T+ c -> error $ "fromWord: Cannot convert " ++ show c ++ " to Nucleotide." derivingUnbox "Nucleotide" [t| Nucleotide -> Word8 |]
src/ELynx/Data/Character/NucleotideI.hs view
@@ -6,7 +6,7 @@ Module : ELynx.Data.NucleotideI Description : Nucleotides with IUPAC characters Copyright : (c) Dominik Schrempf 2018-License : GPL-3+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -48,14 +48,16 @@ -} module ELynx.Data.Character.NucleotideI- ( NucleotideI (..)- ) where+ ( NucleotideI(..)+ )+where import Data.Vector.Unboxed.Deriving import Data.Word8 -import qualified ELynx.Data.Character.Character as C-import ELynx.Tools.ByteString (c2w, w2c)+import qualified ELynx.Data.Character.Character+ as C+import ELynx.Tools -- | NucleotideIs. data NucleotideI = A | C | G | T@@ -89,25 +91,25 @@ fromWord :: Word8 -> NucleotideI fromWord w = case w2c w of- 'A' -> A- 'C' -> C- 'G' -> G- 'T' -> T- 'U' -> U- 'W' -> W- 'S' -> S- 'M' -> M- 'K' -> K- 'R' -> R- 'Y' -> Y- 'B' -> B- 'D' -> D- 'H' -> H- 'V' -> V- 'N' -> N- '-' -> Gap- '.' -> Gap- _ -> error "fromWord: Cannot convert to NucleotideI."+ 'A' -> A+ 'C' -> C+ 'G' -> G+ 'T' -> T+ 'U' -> U+ 'W' -> W+ 'S' -> S+ 'M' -> M+ 'K' -> K+ 'R' -> R+ 'Y' -> Y+ 'B' -> B+ 'D' -> D+ 'H' -> H+ 'V' -> V+ 'N' -> N+ '-' -> Gap+ '.' -> Gap+ _ -> error "fromWord: Cannot convert to NucleotideI." derivingUnbox "NucleotideI" [t| NucleotideI -> Word8 |]@@ -138,7 +140,7 @@ toStandard Gap = [] instance C.CharacterX NucleotideI where- gap = Gap+ gap = Gap instance C.CharacterI NucleotideI where unknown = N
src/ELynx/Data/Character/NucleotideX.hs view
@@ -7,7 +7,7 @@ Description : Extended nucleotides including gaps and unknowns Copyright : (c) Dominik Schrempf 2018 -License : GPL-3+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -34,14 +34,16 @@ -} module ELynx.Data.Character.NucleotideX- ( NucleotideX (..)- ) where+ ( NucleotideX(..)+ )+where import Data.Vector.Unboxed.Deriving import Data.Word8 -import qualified ELynx.Data.Character.Character as C-import ELynx.Tools.ByteString (c2w, w2c)+import qualified ELynx.Data.Character.Character+ as C+import ELynx.Tools -- | Extended nucleotides. data NucleotideX = A | C | G | T@@ -53,17 +55,17 @@ toWord C = c2w 'C' toWord G = c2w 'G' toWord T = c2w 'T'-toWord Gap = c2w '-'+toWord Gap = c2w '-' fromWord :: Word8 -> NucleotideX fromWord w = case w2c w of- 'A' -> A- 'C' -> C- 'G' -> G- 'T' -> T- '-' -> Gap- '.' -> Gap- c -> error $ "fromWord: Cannot convert " ++ show c ++ " to NucleotideX."+ 'A' -> A+ 'C' -> C+ 'G' -> G+ 'T' -> T+ '-' -> Gap+ '.' -> Gap+ c -> error $ "fromWord: Cannot convert " ++ show c ++ " to NucleotideX." derivingUnbox "NucleotideX" [t| NucleotideX -> Word8 |]@@ -75,4 +77,4 @@ fromWord = fromWord instance C.CharacterX NucleotideX where- gap = Gap+ gap = Gap
− src/ELynx/Data/MarkovProcess/AminoAcid.hs
@@ -1,250 +0,0 @@-{- |-Module : ELynx.Data.MarkovProcess.AminoAcid-Description : Amino acid rate matrices such as LG-Copyright : (c) Dominik Schrempf 2019-License : GPL-3--Maintainer : dominik.schrempf@gmail.com-Stability : unstable-Portability : portable--Creation date: Tue Jan 29 09:29:19 2019.--The order of amino acids is alphabetic.---}--module ELynx.Data.MarkovProcess.AminoAcid- ( lg- , lgCustom- , lgCustomUnnormalized- , wag- , wagCustom- , wagCustomUnnormalized- , poisson- , poissonCustom- ) where--import Data.List (elemIndex)-import Data.Maybe (fromMaybe)-import Data.Word (Word8)-import Numeric.LinearAlgebra-import Numeric.SpecFunctions--import ELynx.Data.Alphabet.Alphabet-import ELynx.Data.MarkovProcess.RateMatrix-import qualified ELynx.Data.MarkovProcess.SubstitutionModel as S-import ELynx.Tools.ByteString (c2w)-import ELynx.Tools.LinearAlgebra-import ELynx.Tools.Vector--n :: Int-n = 20---- Some matrices have to be converted from PAML order to alphabetical order. See--- 'pamlToAlphaVec' and 'pamlToAlphaMat'.---- Amno acids in alphabetical order.-aaAlphaOrder :: [Word8]-aaAlphaOrder = map c2w [ 'A', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'K', 'L', 'M'- , 'N', 'P', 'Q', 'R', 'S', 'T', 'V', 'W', 'Y' ]---- Amino acids in PAML oder.-aaPamlOrder :: [Word8]-aaPamlOrder = map c2w [ 'A', 'R', 'N', 'D', 'C', 'Q', 'E', 'G', 'H', 'I', 'L', 'K'- , 'M', 'F', 'P', 'S', 'T', 'W', 'Y', 'V' ]---- -- This is a very slow implementation; since I only convert matrices once it--- -- should not be a problem. A map would be better if performance is an issue.--- pamlIndexToAlphaIndex :: Int -> Int--- pamlIndexToAlphaIndex i = fromMaybe--- (error $ "Could not convert index " ++ show i ++ ".")--- (elemIndex aa aaAlphaOrder)--- where aa = aaPamlOrder !! i---- This is a very slow implementation; since I only convert matrices once it--- should not be a problem. A map would be better if performance is an issue.-alphaIndexToPamlIndex :: Int -> Int-alphaIndexToPamlIndex i = fromMaybe- (error $ "Could not convert index " ++ show i ++ ".")- (elemIndex aa aaPamlOrder)- where aa = aaAlphaOrder !! i---- Convert an amino acid vector in PAML order to a vector in alphabetical order.-pamlToAlphaVec :: Vector R -> Vector R-pamlToAlphaVec v = build n (\i -> v ! alphaIndexToPamlIndex (round i))----- Convert an amino acid matrix in PAML order to a matrix in alphabetical order.-pamlToAlphaMat :: Matrix R -> Matrix R-pamlToAlphaMat m = build (n,n) (\i j -> m- ! alphaIndexToPamlIndex (round i)- ! alphaIndexToPamlIndex (round j))---- The next functions tackle the somewhat stupid, but not easy solvable problem--- of converting a lower triangular matrix (excluding the diagonal) given as a--- list into a symmetric matrix. The diagonal entries are set to zero. This is--- how the exchangeabilities are specified in PAML.---- Conversion from matrix indices (i,j) to list index k.--- (i,j) k------ (0,0) ---- (1,0) 0 (1,1) ---- (2,0) 1 (2,1) 2 (2,2) ---- (3,0) 3 (3,1) 4 (3,2) 5 (3,3) ---- (4,0) 6 (4,1) 7 (4,2) 8 (4,3) 9 (4,4) ------- k = (i choose 2) + j.-ijToK :: Int -> Int -> Int-ijToK i j = round (i `choose` 2) + j---- The function is a little weird because HMatrix uses Double indices for Matrix--- Double builders.-fromListBuilder :: [Double] -> Double -> Double -> Double-fromListBuilder es i j | i > j = es !! ijToK iI jI- | i == j = 0.0- | i < j = es !! ijToK jI iI- | otherwise = error "Float indices could not be compared during matrix creation."- where iI = round i :: Int- jI = round j :: Int---- Exchangeability matrix from list denoting lower triangular matrix, and--- excluding diagonal. This is how the exchangeabilities are specified in PAML.-exchFromList :: [Double] -> ExchangeabilityMatrix-exchFromList es = build (n,n) (fromListBuilder es)---- Lower triangular matrix of LG exchangeabilities in PAML order and in form of--- a list.-lgExchRawPaml :: [Double]-lgExchRawPaml = [0.425093, 0.276818, 0.751878, 0.395144, 0.123954, 5.076149,- 2.489084, 0.534551, 0.528768, 0.062556, 0.969894, 2.807908,- 1.695752, 0.523386, 0.084808, 1.038545, 0.363970, 0.541712,- 5.243870, 0.003499, 4.128591, 2.066040, 0.390192, 1.437645,- 0.844926, 0.569265, 0.267959, 0.348847, 0.358858, 2.426601,- 4.509238, 0.927114, 0.640543, 4.813505, 0.423881, 0.311484,- 0.149830, 0.126991, 0.191503, 0.010690, 0.320627, 0.072854,- 0.044265, 0.008705, 0.108882, 0.395337, 0.301848, 0.068427,- 0.015076, 0.594007, 0.582457, 0.069673, 0.044261, 0.366317,- 4.145067, 0.536518, 6.326067, 2.145078, 0.282959, 0.013266,- 3.234294, 1.807177, 0.296636, 0.697264, 0.159069, 0.137500,- 1.124035, 0.484133, 0.371004, 0.025548, 0.893680, 1.672569,- 0.173735, 0.139538, 0.442472, 4.273607, 6.312358, 0.656604,- 0.253701, 0.052722, 0.089525, 0.017416, 1.105251, 0.035855,- 0.018811, 0.089586, 0.682139, 1.112727, 2.592692, 0.023918,- 1.798853, 1.177651, 0.332533, 0.161787, 0.394456, 0.075382,- 0.624294, 0.419409, 0.196961, 0.508851, 0.078281, 0.249060,- 0.390322, 0.099849, 0.094464, 4.727182, 0.858151, 4.008358,- 1.240275, 2.784478, 1.223828, 0.611973, 1.739990, 0.990012,- 0.064105, 0.182287, 0.748683, 0.346960, 0.361819, 1.338132,- 2.139501, 0.578987, 2.000679, 0.425860, 1.143480, 1.080136,- 0.604545, 0.129836, 0.584262, 1.033739, 0.302936, 1.136863,- 2.020366, 0.165001, 0.571468, 6.472279, 0.180717, 0.593607,- 0.045376, 0.029890, 0.670128, 0.236199, 0.077852, 0.268491,- 0.597054, 0.111660, 0.619632, 0.049906, 0.696175, 2.457121,- 0.095131, 0.248862, 0.140825, 0.218959, 0.314440, 0.612025,- 0.135107, 1.165532, 0.257336, 0.120037, 0.054679, 5.306834,- 0.232523, 0.299648, 0.131932, 0.481306, 7.803902, 0.089613,- 0.400547, 0.245841, 3.151815, 2.547870, 0.170887, 0.083688,- 0.037967, 1.959291, 0.210332, 0.245034, 0.076701, 0.119013,- 10.649107, 1.702745, 0.185202, 1.898718, 0.654683, 0.296501,- 0.098369, 2.188158, 0.189510, 0.249313]----- Exchangeabilities of LG model in alphabetical order.-lgExch :: ExchangeabilityMatrix-lgExch = pamlToAlphaMat $ exchFromList lgExchRawPaml---- Stationary distribution in PAML order.-lgStatDistPaml :: StationaryDistribution-lgStatDistPaml = normalizeSumVec 1.0 $- fromList [ 0.079066, 0.055941, 0.041977, 0.053052, 0.012937, 0.040767- , 0.071586, 0.057337, 0.022355, 0.062157, 0.099081, 0.064600- , 0.022951, 0.042302, 0.044040, 0.061197, 0.053287, 0.012066- , 0.034155, 0.069147 ]---- Stationary distribution of LG model in alphabetical order.-lgStatDist :: StationaryDistribution-lgStatDist = pamlToAlphaVec lgStatDistPaml---- | LG substitution model.-lg :: S.SubstitutionModel-lg = S.substitutionModel Protein "LG" [] lgStatDist lgExch---- | LG substitution model with maybe a name and a custom stationary distribution.-lgCustom :: Maybe String -> StationaryDistribution -> S.SubstitutionModel-lgCustom mn d = S.substitutionModel Protein name [] d lgExch- where name = fromMaybe "LG-Custom" mn---- | LG substitution model with maybe a name and a custom stationary distribution.-lgCustomUnnormalized :: Maybe String -> StationaryDistribution -> S.SubstitutionModel-lgCustomUnnormalized mn d = S.unnormalized Protein name [] d lgExch- where name = fromMaybe "LG-Custom-Unnormalized" mn---- WAG exchangeability list in PAML order.-wagExchRawPaml :: [Double]-wagExchRawPaml =- [ 55.15710- , 50.98480, 63.53460- , 73.89980, 14.73040, 542.94200- , 102.70400, 52.81910, 26.52560, 3.02949- , 90.85980, 303.55000, 154.36400, 61.67830, 9.88179- , 158.28500, 43.91570, 94.71980, 617.41600, 2.13520, 546.94700- , 141.67200, 58.46650, 112.55600, 86.55840, 30.66740, 33.00520, 56.77170- , 31.69540, 213.71500, 395.62900, 93.06760, 24.89720, 429.41100, 57.00250, 24.94100- , 19.33350, 18.69790, 55.42360, 3.94370, 17.01350, 11.39170, 12.73950, 3.04501, 13.81900- , 39.79150, 49.76710, 13.15280, 8.48047, 38.42870, 86.94890, 15.42630, 6.13037, 49.94620, 317.09700- , 90.62650, 535.14200, 301.20100, 47.98550, 7.40339, 389.49000, 258.44300, 37.35580, 89.04320, 32.38320, 25.75550- , 89.34960, 68.31620, 19.82210, 10.37540, 39.04820, 154.52600, 31.51240, 17.41000, 40.41410, 425.74600, 485.40200, 93.42760- , 21.04940, 10.27110, 9.61621, 4.67304, 39.80200, 9.99208, 8.11339, 4.99310, 67.93710, 105.94700, 211.51700, 8.88360, 119.06300- , 143.85500, 67.94890, 19.50810, 42.39840, 10.94040, 93.33720, 68.23550, 24.35700, 69.61980, 9.99288, 41.58440, 55.68960, 17.13290, 16.14440- , 337.07900, 122.41900, 397.42300, 107.17600, 140.76600, 102.88700, 70.49390, 134.18200, 74.01690, 31.94400, 34.47390, 96.71300, 49.39050, 54.59310, 161.32800- , 212.11100, 55.44130, 203.00600, 37.48660, 51.29840, 85.79280, 82.27650, 22.58330, 47.33070, 145.81600, 32.66220, 138.69800, 151.61200, 17.19030, 79.53840, 437.80200- , 11.31330, 116.39200, 7.19167, 12.97670, 71.70700, 21.57370, 15.65570, 33.69830, 26.25690, 21.24830, 66.53090, 13.75050, 51.57060, 152.96400, 13.94050, 52.37420, 11.08640- , 24.07350, 38.15330, 108.60000, 32.57110, 54.38330, 22.77100, 19.63030, 10.36040, 387.34400, 42.01700, 39.86180, 13.32640, 42.84370, 645.42800, 21.60460, 78.69930, 29.11480, 248.53900- , 200.60100, 25.18490, 19.62460, 15.23350, 100.21400, 30.12810, 58.87310, 18.72470, 11.83580, 782.13000, 180.03400, 30.54340, 205.84500, 64.98920, 31.48870, 23.27390, 138.82300, 36.53690, 31.47300 ]---- WAG exchangeability matrix n alphabetical order.-wagExch :: ExchangeabilityMatrix-wagExch = pamlToAlphaMat $ exchFromList wagExchRawPaml---- WAG stationary distribution in PAML order.-wagStatDistPaml :: StationaryDistribution-wagStatDistPaml = normalizeSumVec 1.0 $- fromList [ 0.0866279, 0.043972, 0.0390894, 0.0570451, 0.0193078, 0.0367281- , 0.0580589, 0.0832518, 0.0244313, 0.048466, 0.086209, 0.0620286- , 0.0195027, 0.0384319, 0.0457631, 0.0695179, 0.0610127, 0.0143859- , 0.0352742, 0.0708957 ]---- WAG stationary distribution in alphabetical order.-wagStatDist :: StationaryDistribution-wagStatDist = pamlToAlphaVec wagStatDistPaml---- | LG substitution model.-wag :: S.SubstitutionModel-wag = S.substitutionModel Protein "WAG" [] wagStatDist wagExch---- | LG substitution model with maybe a name and a custom stationary distribution.-wagCustom :: Maybe String -> StationaryDistribution -> S.SubstitutionModel-wagCustom mn d = S.substitutionModel Protein name [] d wagExch- where name = fromMaybe "WAG-Custom" mn---- | LG substitution model with maybe a name and a custom stationary distribution.-wagCustomUnnormalized :: Maybe String -> StationaryDistribution -> S.SubstitutionModel-wagCustomUnnormalized mn d = S.unnormalized Protein name [] d wagExch- where name = fromMaybe "WAG-Custom-Unnormalized" mn--uniformExch :: ExchangeabilityMatrix-uniformExch = matrixSetDiagToZero $ matrix n $ replicate (n*n) 1.0--poissonExch :: ExchangeabilityMatrix-poissonExch = uniformExch---- | Poisson substitution model.-poisson :: S.SubstitutionModel-poisson = S.substitutionModel Protein "Poisson" [] (uniformVec n) poissonExch---- | Poisson substitution model with maybe a name and a custom stationary distribution.-poissonCustom :: Maybe String -> StationaryDistribution -> S.SubstitutionModel-poissonCustom mn d = S.substitutionModel Protein name [] d poissonExch- where name = fromMaybe "Poisson-Custom" mn
− src/ELynx/Data/MarkovProcess/CXXModels.hs
@@ -1,128 +0,0 @@-{- |-Module : ELynx.Data.MarkovProcess.CXXModels-Description : C10 to C60 models-Copyright : (c) Dominik Schrempf 2019-License : GPL-3--Maintainer : dominik.schrempf@gmail.com-Stability : unstable-Portability : portable--Creation date: Tue Feb 26 16:44:33 2019.--Quang, L. S., Gascuel, O. & Lartillot, N. Empirical profile mixture models for-phylogenetic reconstruction. Bioinformatics 24, 2317–2323 (2008).--XXX: For now, I only provide Poisson exchangeabilities.---}--module ELynx.Data.MarkovProcess.CXXModels- (- cxx- ) where--import ELynx.Data.MarkovProcess.AminoAcid-import ELynx.Data.MarkovProcess.CXXModelsData-import qualified ELynx.Data.MarkovProcess.MixtureModel as M-import ELynx.Data.MarkovProcess.RateMatrix-import qualified ELynx.Data.MarkovProcess.SubstitutionModel as S---- | Create CXX model with given number of components and probably with custom--- weights.-cxx :: Int -> Maybe [M.Weight] -> Maybe M.MixtureModel-cxx 10 (Just ws) = Just $ c10CustomWeights ws-cxx 20 (Just ws) = Just $ c20CustomWeights ws-cxx 30 (Just ws) = Just $ c30CustomWeights ws-cxx 40 (Just ws) = Just $ c40CustomWeights ws-cxx 50 (Just ws) = Just $ c50CustomWeights ws-cxx 60 (Just ws) = Just $ c60CustomWeights ws-cxx 10 Nothing = Just c10-cxx 20 Nothing = Just c20-cxx 30 Nothing = Just c30-cxx 40 Nothing = Just c40-cxx 50 Nothing = Just c50-cxx 60 Nothing = Just c60-cxx _ _ = Nothing---- | C10 model.-c10 :: M.MixtureModel-c10 = cxxFromStatDistsAndWeights c10Weights c10StatDists---- | C20 model.-c20 :: M.MixtureModel-c20 = cxxFromStatDistsAndWeights c20Weights c20StatDists---- | C30 model.-c30 :: M.MixtureModel-c30 = cxxFromStatDistsAndWeights c30Weights c30StatDists---- | C40 model.-c40 :: M.MixtureModel-c40 = cxxFromStatDistsAndWeights c40Weights c40StatDists---- | C50 model.-c50 :: M.MixtureModel-c50 = cxxFromStatDistsAndWeights c50Weights c50StatDists---- | C60 model.-c60 :: M.MixtureModel-c60 = cxxFromStatDistsAndWeights c60Weights c60StatDists---- | C10 model with custom weights.-c10CustomWeights :: [M.Weight] -> M.MixtureModel-c10CustomWeights ws- | length ws == 10 = cxxFromStatDistsAndWeights ws c10StatDists- | otherwise = error "Number of weights does not match C10 model."---- | C20 model with custom weights.-c20CustomWeights :: [M.Weight] -> M.MixtureModel-c20CustomWeights ws- | length ws == 20 = cxxFromStatDistsAndWeights ws c20StatDists- | otherwise = error "Number of weights does not match C20 model."---- | C30 model with custom weights.-c30CustomWeights :: [M.Weight] -> M.MixtureModel-c30CustomWeights ws- | length ws == 30 = cxxFromStatDistsAndWeights ws c30StatDists- | otherwise = error "Number of weights does not match C30 model."---- | C40 model with custom weights.-c40CustomWeights :: [M.Weight] -> M.MixtureModel-c40CustomWeights ws- | length ws == 40 = cxxFromStatDistsAndWeights ws c40StatDists- | otherwise = error "Number of weights does not match C40 model."---- | C50 model with custom weights.-c50CustomWeights :: [M.Weight] -> M.MixtureModel-c50CustomWeights ws- | length ws == 50 = cxxFromStatDistsAndWeights ws c50StatDists- | otherwise = error "Number of weights does not match C50 model."---- | C60 model with custom weights.-c60CustomWeights :: [M.Weight] -> M.MixtureModel-c60CustomWeights ws- | length ws == 60 = cxxFromStatDistsAndWeights ws c60StatDists- | otherwise = error "Number of weights does not match C60 model."--cxxName :: Int -> String-cxxName nComps = 'C' : show nComps--componentName :: Int -> Int -> String-componentName nComps comp = cxxName nComps ++ "; component " ++ show comp---- Keep in mind, that when using different exchangeabilities, I have to decide--- about global or local normalization.-cxxSubstitutionModelFromStatDist :: Int -> Int -> StationaryDistribution -> S.SubstitutionModel-cxxSubstitutionModelFromStatDist nComps comp = poissonCustom (Just name)- where name = componentName nComps comp--cxxSubstitutionModelsFromStatDists :: [StationaryDistribution] -> [S.SubstitutionModel]-cxxSubstitutionModelsFromStatDists ds = zipWith (cxxSubstitutionModelFromStatDist nComp) [1..] ds- where nComp = length ds--cxxFromStatDistsAndWeights :: [M.Weight] -> [StationaryDistribution] -> M.MixtureModel-cxxFromStatDistsAndWeights ws ds = M.MixtureModel (cxxName nComps) comps- where- nComps = length ds- comps = zipWith M.Component ws (cxxSubstitutionModelsFromStatDists ds)
− src/ELynx/Data/MarkovProcess/CXXModelsData.hs
@@ -1,294 +0,0 @@-{- |-Module : ELynx.Data.MarkovProcess.CXXModelsData-Description : Stationary distributions and weights-Copyright : (c) Dominik Schrempf 2019-License : GPL-3--Maintainer : dominik.schrempf@gmail.com-Stability : unstable-Portability : portable--Creation date: Tue Feb 26 17:17:35 2019.--Quang, L. S., Gascuel, O. & Lartillot, N. Empirical profile mixture models for-phylogenetic reconstruction. Bioinformatics 24, 2317–2323 (2008).---}--module ELynx.Data.MarkovProcess.CXXModelsData- (- c10StatDists- , c10Weights- , c20StatDists- , c20Weights- , c30StatDists- , c30Weights- , c40StatDists- , c40Weights- , c50StatDists- , c50Weights- , c60StatDists- , c60Weights- ) where--import qualified Data.Vector.Storable as V--import ELynx.Data.MarkovProcess.RateMatrix---- | Stationary distribution of C10 model.-c10StatDists :: [StationaryDistribution]-c10StatDists = map V.fromList- [[0.408257, 0.0349388, 0.00698709, 0.00978467, 0.00616043, 0.122161, 0.00391518, 0.0125784, 0.00596702, 0.0158339, 0.00813132, 0.00962854, 0.0394156, 0.00752797, 0.0081783, 0.168245, 0.0658133, 0.0604427, 0.00187516, 0.00415797],- [0.102776, 0.0149663, 0.0155944, 0.0419667, 0.0180729, 0.0138806, 0.0158865, 0.106608, 0.0436344, 0.113194, 0.04378, 0.0213272, 0.0223251, 0.0440685, 0.0418664, 0.0529608, 0.108174, 0.160665, 0.00451472, 0.0137374],- [0.0351766, 0.00787065, 0.000676874, 0.00196868, 0.0126221, 0.00224206, 0.00128783, 0.351582, 0.00188565, 0.127818, 0.0242632, 0.00165915, 0.00297716, 0.00165596, 0.00196786, 0.00499981, 0.0255378, 0.388864, 0.00119078, 0.00375393],- [0.0408514, 0.00376029, 0.233381, 0.0901239, 0.00251082, 0.115833, 0.0373197, 0.00255236, 0.0485017, 0.00521646, 0.00225718, 0.218565, 0.0108334, 0.0380451, 0.0269887, 0.0804527, 0.030288, 0.00444811, 0.00108153, 0.00698909],- [0.0185493, 0.00704165, 0.000977506, 0.00248916, 0.073333, 0.00289529, 0.0040104, 0.163242, 0.00435709, 0.444308, 0.120282, 0.00248957, 0.00488276, 0.00835394, 0.00623624, 0.00516424, 0.0131807, 0.0968581, 0.00687598, 0.0144734],- [0.110675, 0.00148349, 0.163644, 0.263846, 0.00232568, 0.0325228, 0.0163804, 0.00683349, 0.0677158, 0.014068, 0.00489881, 0.0405186, 0.0298982, 0.0877962, 0.035219, 0.0562888, 0.0426922, 0.0181079, 0.0010339, 0.00405223],- [0.0522658, 0.0143325, 0.0297745, 0.0388387, 0.0624033, 0.0228101, 0.155164, 0.0187406, 0.0439469, 0.065378, 0.0207189, 0.0714837, 0.0145475, 0.073654, 0.0668295, 0.0549018, 0.037014, 0.0267512, 0.0193757, 0.111069],- [0.0116587, 0.0105341, 0.00217425, 0.00242511, 0.365099, 0.00347091, 0.0366787, 0.0187185, 0.00266947, 0.067649, 0.0143535, 0.00640111, 0.00311599, 0.00402037, 0.00509901, 0.00948485, 0.00737139, 0.0206341, 0.0509565, 0.357486],- [0.0627196, 0.00526629, 0.0236193, 0.0686285, 0.00391818, 0.0256175, 0.0332612, 0.0128968, 0.227084, 0.0305628, 0.0124037, 0.0428629, 0.0140441, 0.109811, 0.203878, 0.0483152, 0.0463378, 0.0197063, 0.00251435, 0.00655211],- [0.114552, 0.00985495, 0.0416192, 0.0364908, 0.0046606, 0.0503818, 0.0165233, 0.00929495, 0.0423027, 0.0139154, 0.00822408, 0.0750615, 0.0379222, 0.0339625, 0.0324009, 0.261065, 0.184583, 0.0195769, 0.0017549, 0.00585383]]---- | Weights of C10 model.-c10Weights :: [Double]-c10Weights = [0.119134, 0.0874372, 0.103711, 0.0922585, 0.107049, 0.132995, 0.0538028, 0.0691986, 0.131994, 0.10242]---- | Stationary distribution of C20 model.-c20StatDists :: [StationaryDistribution]-c20StatDists = map V.fromList- [[0.0862413, 0.0130505, 0.0329909, 0.0184527, 0.00441553, 0.0366905, 0.0108013, 0.00979071, 0.0220195, 0.0112826, 0.00878215, 0.0791293, 0.0189273, 0.0169047, 0.0171944, 0.317815, 0.27117, 0.0179753, 0.00153173, 0.00483429],- [0.203558, 0.0348667, 0.00316561, 0.00708594, 0.0112429, 0.0195236, 0.0024392, 0.115257, 0.00423808, 0.0789777, 0.0309187, 0.00770524, 0.0164189, 0.00640441, 0.00509808, 0.0496777, 0.111895, 0.284906, 0.00177626, 0.00484482],- [0.0211547, 0.00481886, 0.000549287, 0.00145396, 0.0128252, 0.00114309, 0.00113464, 0.392846, 0.00135799, 0.125064, 0.0209789, 0.0012755, 0.00202472, 0.00123288, 0.00149462, 0.00262407, 0.0171914, 0.386068, 0.00115911, 0.0036028],- [0.0376904, 0.00640738, 0.0109469, 0.0358365, 0.00363498, 0.0191107, 0.0329514, 0.0101712, 0.289763, 0.0237496, 0.00965289, 0.0365411, 0.0105337, 0.0893564, 0.28852, 0.0356314, 0.0355927, 0.0144622, 0.00279252, 0.00665572],- [0.00845978, 0.0084909, 0.00244879, 0.00250555, 0.342046, 0.00242771, 0.0433214, 0.0097713, 0.0026741, 0.0380507, 0.00807248, 0.00725259, 0.00214187, 0.00427815, 0.00535899, 0.00804189, 0.00553221, 0.012141, 0.049484, 0.4375],- [0.17599, 0.00175587, 0.130126, 0.218217, 0.0025277, 0.0409535, 0.0130708, 0.00856221, 0.0542946, 0.0159531, 0.00540458, 0.0332846, 0.037102, 0.0707184, 0.0290429, 0.0793481, 0.0540083, 0.0249553, 0.00105921, 0.00362591],- [0.16344, 0.00886599, 0.0374273, 0.0220612, 0.00306413, 0.529672, 0.00900061, 0.00175694, 0.0167118, 0.00611563, 0.00293908, 0.0438702, 0.0126458, 0.0137555, 0.0195541, 0.0829343, 0.0142836, 0.00579857, 0.00286407, 0.00323983],- [0.0917469, 0.0284015, 0.0133819, 0.0196876, 0.0998479, 0.0249899, 0.0449766, 0.0583556, 0.0164916, 0.115501, 0.0395995, 0.0290699, 0.0209916, 0.0255085, 0.0265853, 0.0736483, 0.0661518, 0.0831856, 0.0246464, 0.0972327],- [0.0646701, 0.00771176, 0.0168734, 0.0544978, 0.0219148, 0.0148894, 0.0313852, 0.0505983, 0.0907931, 0.184428, 0.077484, 0.0228907, 0.0105004, 0.0996415, 0.0988016, 0.0321196, 0.0411766, 0.0505824, 0.0084303, 0.0206106],- [0.0135994, 0.010009, 0.00079517, 0.00180118, 0.264097, 0.00267946, 0.00724019, 0.0814027, 0.00251581, 0.366142, 0.0734965, 0.00184694, 0.00389941, 0.00464208, 0.00434084, 0.00436688, 0.00752485, 0.0573473, 0.0261565, 0.0660971],- [0.147804, 0.00488258, 0.0534743, 0.0727246, 0.00299039, 0.0907726, 0.0262289, 0.00357811, 0.105166, 0.0126777, 0.00596218, 0.072663, 0.0156558, 0.0757166, 0.0842845, 0.14599, 0.0634877, 0.00927198, 0.00159285, 0.00507607],- [0.0186377, 0.00549689, 0.00083297, 0.00202485, 0.0385383, 0.00217135, 0.0023666, 0.202081, 0.00291207, 0.437038, 0.124186, 0.00198652, 0.00406723, 0.00658901, 0.00420552, 0.00461774, 0.0149904, 0.118938, 0.00268717, 0.00563241],- [0.0477624, 0.00757917, 0.0141349, 0.0462688, 0.0130691, 0.00523279, 0.0165352, 0.17415, 0.0577575, 0.112125, 0.0330288, 0.0209574, 0.0124375, 0.0429297, 0.0505743, 0.0264989, 0.0951755, 0.20937, 0.00316605, 0.0112466],- [0.416419, 0.0406938, 0.00451317, 0.00632298, 0.00484384, 0.0946185, 0.00310574, 0.00764432, 0.00389418, 0.00998854, 0.00693232, 0.00917014, 0.0187841, 0.00613205, 0.00561008, 0.236077, 0.0746275, 0.0459225, 0.00121726, 0.00348258],- [0.0402296, 0.0124783, 0.0365524, 0.0372197, 0.0459095, 0.0233618, 0.210831, 0.00934787, 0.0482411, 0.0360561, 0.010029, 0.103665, 0.0098504, 0.0826558, 0.0735203, 0.0533383, 0.0310209, 0.015248, 0.0140077, 0.106438],- [0.0323453, 0.00359763, 0.24315, 0.0710274, 0.00244293, 0.101607, 0.0366225, 0.00314108, 0.0470129, 0.00519805, 0.00240287, 0.252045, 0.00948378, 0.0330831, 0.0236283, 0.0848355, 0.0359083, 0.00487046, 0.000873093, 0.00672477],- [0.147626, 0.00323272, 0.0403052, 0.0576893, 0.00471772, 0.0330851, 0.0146393, 0.0108267, 0.0451351, 0.0256201, 0.00586514, 0.0211973, 0.347371, 0.0371554, 0.0334507, 0.0892065, 0.0485899, 0.0282336, 0.00163587, 0.00441772],- [0.103145, 0.00617625, 0.0386402, 0.0923369, 0.00676664, 0.0202338, 0.0246762, 0.0376904, 0.0921699, 0.0376284, 0.0161883, 0.0435172, 0.0128302, 0.0786603, 0.0717748, 0.095145, 0.137857, 0.0740454, 0.00221447, 0.00830416],- [0.0837543, 0.00207351, 0.0804871, 0.194776, 0.00230634, 0.022903, 0.0268459, 0.00740798, 0.145929, 0.019025, 0.00673952, 0.0518811, 0.0085616, 0.14565, 0.0899383, 0.045574, 0.0451081, 0.0150303, 0.00107713, 0.00493253],- [0.0578736, 0.00111308, 0.294674, 0.34021, 0.00170349, 0.0293911, 0.0139817, 0.00305257, 0.0363365, 0.00626119, 0.0027296, 0.0491422, 0.0156106, 0.059825, 0.0138314, 0.0358045, 0.0249942, 0.00876742, 0.000866434, 0.0038313]]---- | Weights of C20 model.-c20Weights :: [Double]-c20Weights = [0.0559911, 0.0514825, 0.0812922, 0.0721977, 0.0556719, 0.0331003, 0.0589502, 0.0263757, 0.0307584, 0.0376701, 0.0303058, 0.0808776, 0.0263349, 0.0579101, 0.0371248, 0.0586868, 0.0561479, 0.0349811, 0.0544937, 0.0596472]---- | Stationary distribution of C30 model.-c30StatDists :: [StationaryDistribution]-c30StatDists = map V.fromList- [[0.110045, 0.00190472, 0.159541, 0.109896, 0.00166295, 0.0684302, 0.0137951, 0.00262831, 0.0358554, 0.00733965, 0.00247064, 0.0640338, 0.166936, 0.0310187, 0.0171295, 0.138179, 0.0568343, 0.00823656, 0.000466112, 0.00359702],- [0.0874125, 0.00498264, 0.032612, 0.0951701, 0.00489966, 0.0144043, 0.0210627, 0.0399884, 0.11472, 0.0301585, 0.0126489, 0.0382152, 0.0137397, 0.0798169, 0.080632, 0.087377, 0.155862, 0.0793881, 0.00151228, 0.00539745],- [0.0225477, 0.00500182, 0.000595928, 0.00150305, 0.0089216, 0.0011571, 0.000937432, 0.394469, 0.00136009, 0.0889573, 0.0189103, 0.00130346, 0.0018312, 0.00114366, 0.00149005, 0.00283364, 0.0189813, 0.425056, 0.000669375, 0.00233037],- [0.0602158, 0.000952546, 0.290008, 0.361087, 0.00146256, 0.0281926, 0.0130501, 0.00305162, 0.0352705, 0.00604019, 0.00274606, 0.0414988, 0.0127175, 0.0621611, 0.0136833, 0.0318109, 0.022528, 0.00932584, 0.000794803, 0.00340246],- [0.0101224, 0.00859894, 0.000637919, 0.00112496, 0.278538, 0.00240852, 0.00477534, 0.0701153, 0.00167485, 0.413591, 0.0744863, 0.00129289, 0.00404666, 0.00350286, 0.00283449, 0.00370872, 0.00523793, 0.040886, 0.0200223, 0.0523939],- [0.133583, 0.00105145, 0.112578, 0.209957, 0.00167936, 0.0207552, 0.012133, 0.00735265, 0.0771772, 0.0133278, 0.00305717, 0.0213892, 0.18902, 0.0565844, 0.028479, 0.0484054, 0.0373318, 0.0225174, 0.000926699, 0.00269464],- [0.0408277, 0.0153918, 0.00306349, 0.00660109, 0.157505, 0.00581131, 0.0245212, 0.148751, 0.00759232, 0.16378, 0.0385527, 0.00804649, 0.00583522, 0.0102922, 0.0124492, 0.0151579, 0.033222, 0.154771, 0.0264937, 0.121334],- [0.246906, 0.103941, 0.00274183, 0.00549448, 0.0251776, 0.0373263, 0.00857523, 0.0292404, 0.00561231, 0.0535091, 0.0302246, 0.016893, 0.00780989, 0.0103988, 0.0106279, 0.164235, 0.123989, 0.0955868, 0.00531559, 0.0163954],- [0.0549429, 0.0099281, 0.00929153, 0.0417085, 0.024386, 0.0105564, 0.0363512, 0.0569585, 0.115252, 0.168183, 0.0592328, 0.0202958, 0.00830554, 0.0906036, 0.130543, 0.0283779, 0.0412594, 0.0592101, 0.00963554, 0.024978],- [0.0462773, 0.0172727, 0.0182504, 0.0224266, 0.133632, 0.0160971, 0.135785, 0.0164967, 0.0239396, 0.0598936, 0.0164507, 0.0412365, 0.0117413, 0.0348991, 0.0362984, 0.0454156, 0.0304388, 0.0253979, 0.0330338, 0.235016],- [0.047438, 0.00823324, 0.0112117, 0.0388101, 0.0116644, 0.00559986, 0.0149157, 0.183217, 0.0467851, 0.110069, 0.0356444, 0.0222454, 0.0100245, 0.0374051, 0.041018, 0.0317171, 0.111435, 0.219931, 0.00266856, 0.00996601],- [0.0213608, 0.00610249, 0.00129412, 0.00362973, 0.0371808, 0.0030017, 0.00384259, 0.130947, 0.00545678, 0.456699, 0.194784, 0.0039879, 0.00407473, 0.0139566, 0.00699762, 0.00769915, 0.0198019, 0.0693208, 0.00340864, 0.00645457],- [0.0919632, 0.0117031, 0.0306717, 0.0171908, 0.00415894, 0.0370685, 0.0100793, 0.00931237, 0.0205386, 0.0097241, 0.00757673, 0.0764682, 0.0179686, 0.016006, 0.0160005, 0.325447, 0.274438, 0.0178218, 0.00138874, 0.00447397],- [0.464925, 0.0233329, 0.00507317, 0.00579942, 0.00255882, 0.149524, 0.00232984, 0.00433612, 0.00285254, 0.00559955, 0.00393132, 0.00753048, 0.0186467, 0.00435713, 0.00430132, 0.215019, 0.047703, 0.0292668, 0.00090381, 0.00200872],- [0.205133, 0.00209159, 0.107098, 0.198973, 0.00182351, 0.0487574, 0.0127143, 0.00581247, 0.0667787, 0.0133472, 0.00437834, 0.0339418, 0.0110998, 0.0822742, 0.0439661, 0.0873962, 0.0519782, 0.0193174, 0.00073616, 0.00238214],- [0.026369, 0.00394624, 0.311772, 0.0694354, 0.00233519, 0.0799842, 0.0309615, 0.00275212, 0.0288631, 0.00385797, 0.00183637, 0.272716, 0.00624576, 0.0218979, 0.0133614, 0.0798729, 0.0324143, 0.00427053, 0.000722966, 0.00638577],- [0.15265, 0.00402153, 0.0237392, 0.0377949, 0.00585844, 0.0306214, 0.014193, 0.0123318, 0.0375543, 0.0290062, 0.00644732, 0.016823, 0.38645, 0.0341734, 0.0332784, 0.0880336, 0.0489543, 0.0313553, 0.00182526, 0.00488778],- [0.00802476, 0.00234121, 0.000338596, 0.000781895, 0.0311968, 0.000540383, 0.00100264, 0.317706, 0.00125983, 0.373789, 0.0488212, 0.00063274, 0.00206875, 0.00155079, 0.00174086, 0.00120951, 0.00656968, 0.19446, 0.00163092, 0.00433436],- [0.059995, 0.0016885, 0.0889731, 0.0962701, 0.00118007, 0.0588796, 0.0327277, 0.00214673, 0.182586, 0.00708764, 0.00339794, 0.133492, 0.0141408, 0.0864856, 0.100054, 0.0779002, 0.0448817, 0.00422414, 0.000624903, 0.00326411],- [0.0393521, 0.00990719, 0.0431198, 0.039464, 0.0188681, 0.0282866, 0.209572, 0.00761014, 0.0596434, 0.025834, 0.00845867, 0.142548, 0.00968405, 0.101979, 0.0838171, 0.0624999, 0.0347088, 0.0127686, 0.00546458, 0.0564145],- [0.00727155, 0.00678525, 0.00276038, 0.00250694, 0.127148, 0.00808347, 0.0113218, 0.00566096, 0.00177351, 0.03942, 0.00798661, 0.00197568, 0.00410981, 0.00433393, 0.0140999, 0.00522444, 0.00430223, 0.00942264, 0.627357, 0.108456],- [0.090707, 0.00636999, 0.0745873, 0.0386802, 0.00288592, 0.475005, 0.016809, 0.00147211, 0.0343856, 0.00551498, 0.00246921, 0.0860678, 0.0112151, 0.0259377, 0.0290062, 0.073111, 0.0153706, 0.00392819, 0.00229148, 0.00418607],- [0.00552919, 0.00721059, 0.00114134, 0.000968304, 0.374409, 0.0016071, 0.0325831, 0.00829184, 0.00138494, 0.0353678, 0.00744868, 0.00460866, 0.00133746, 0.00226922, 0.00246263, 0.00574027, 0.00372796, 0.00988026, 0.0330334, 0.460998],- [0.244326, 0.027388, 0.00315909, 0.00767156, 0.00606194, 0.019609, 0.00201894, 0.101744, 0.00454923, 0.0468424, 0.0201286, 0.00624227, 0.0185219, 0.0053594, 0.00453866, 0.0497754, 0.11708, 0.310493, 0.000957725, 0.00353337],- [0.0863111, 0.00244198, 0.0600903, 0.20361, 0.00293931, 0.0175221, 0.0245245, 0.0105994, 0.148579, 0.027121, 0.00958244, 0.0313963, 0.00682768, 0.167235, 0.0984812, 0.03478, 0.0408211, 0.0202271, 0.00140013, 0.00551054],- [0.0643926, 0.00231659, 0.162821, 0.17627, 0.00430667, 0.0297139, 0.0303504, 0.0088163, 0.072714, 0.0148017, 0.00567484, 0.103121, 0.00992703, 0.0752535, 0.0369049, 0.0926434, 0.083313, 0.0161551, 0.00112371, 0.00938015],- [0.173668, 0.00835966, 0.0285985, 0.0483894, 0.0058274, 0.0781901, 0.0266135, 0.00686419, 0.0964012, 0.0219499, 0.0112303, 0.0520405, 0.0169661, 0.0722447, 0.0943629, 0.15478, 0.0751702, 0.0172524, 0.00287745, 0.00821304],- [0.0347857, 0.00627549, 0.00923552, 0.0323568, 0.00343035, 0.0170288, 0.0306439, 0.00919323, 0.302085, 0.0224429, 0.00937208, 0.0314157, 0.0104447, 0.0861073, 0.307598, 0.0326883, 0.0328713, 0.0130832, 0.00252449, 0.00641713],- [0.108774, 0.0247877, 0.00158232, 0.00417699, 0.0316523, 0.0134703, 0.00247658, 0.164346, 0.00270004, 0.233715, 0.0539213, 0.00326798, 0.0154887, 0.0057932, 0.0051781, 0.0188188, 0.0474912, 0.251283, 0.00376565, 0.00731064],- [0.110101, 0.00726998, 0.0579269, 0.0828181, 0.0269154, 0.0314463, 0.0308557, 0.0530866, 0.029386, 0.109679, 0.0458729, 0.0435099, 0.0296431, 0.0615197, 0.0324325, 0.0715888, 0.0685882, 0.0754042, 0.00623241, 0.0257238]]---- | Weights of C30 model.-c30Weights :: [Double]-c30Weights = [0.00957833, 0.0248476, 0.0636309, 0.0537939, 0.0295886, 0.0117588, 0.0132013, 0.0236869, 0.0261688, 0.0239822, 0.0257101, 0.0465072, 0.0546795, 0.0536085, 0.0270623, 0.0403914, 0.0474213, 0.0458816, 0.0214037, 0.0290386, 0.0123392, 0.056935, 0.0419688, 0.0339027, 0.0388777, 0.0196344, 0.0233086, 0.0622723, 0.0184803, 0.0203395]---- | Stationary distribution of C40 model.-c40StatDists :: [StationaryDistribution]-c40StatDists = map V.fromList- [[0.066026, 0.00565867, 0.105447, 0.0440361, 0.00131048, 0.0711239, 0.0168195, 0.00390887, 0.036669, 0.0055316, 0.00374124, 0.159982, 0.0176359, 0.0273928, 0.0231862, 0.249769, 0.150708, 0.0065529, 0.000672321, 0.00382902],- [0.0232377, 0.00379875, 0.353209, 0.0739378, 0.00240321, 0.0576668, 0.0315867, 0.00310928, 0.0259363, 0.00387116, 0.00173556, 0.275965, 0.00631169, 0.0197339, 0.0122683, 0.0657068, 0.0270484, 0.00475317, 0.000760289, 0.00696025],- [0.0166487, 0.00366657, 0.000565145, 0.00133563, 0.00827757, 0.000889475, 0.000823185, 0.412937, 0.00119041, 0.0884689, 0.0186055, 0.00126222, 0.001403, 0.00106698, 0.00125948, 0.00213394, 0.0162167, 0.420686, 0.000608205, 0.00195532],- [0.239474, 0.0283812, 0.00447417, 0.010553, 0.00559911, 0.013511, 0.00389298, 0.0765957, 0.0071093, 0.0358495, 0.0199496, 0.0120537, 0.0114266, 0.00865589, 0.00729013, 0.0847799, 0.179728, 0.245468, 0.0009838, 0.00422407],- [0.119461, 0.0150527, 0.0134273, 0.0192173, 0.0550467, 0.0337676, 0.0214746, 0.0579002, 0.0147261, 0.144631, 0.0561243, 0.0294552, 0.0631355, 0.0301538, 0.0233256, 0.0925267, 0.083123, 0.0811758, 0.0131636, 0.0331118],- [0.0567044, 0.00089248, 0.29555, 0.379515, 0.00129723, 0.023047, 0.0118361, 0.0031182, 0.0314206, 0.00601375, 0.00285841, 0.0364734, 0.0124746, 0.0609517, 0.0117359, 0.0300335, 0.0227051, 0.00946396, 0.000773876, 0.00313438],- [0.0179027, 0.016076, 0.000887041, 0.00231821, 0.334486, 0.00398298, 0.0127293, 0.0404651, 0.00279947, 0.167614, 0.0424172, 0.00356977, 0.00201151, 0.00453955, 0.00409671, 0.00758416, 0.00682273, 0.0326045, 0.0518381, 0.245254],- [0.271217, 0.200383, 0.0021017, 0.002323, 0.020299, 0.0502501, 0.0053728, 0.0150685, 0.00206463, 0.0330003, 0.0154811, 0.0141045, 0.0045351, 0.00482641, 0.00564808, 0.17642, 0.0839578, 0.0741934, 0.00462652, 0.0141271],- [0.0894737, 0.00455383, 0.0272183, 0.127508, 0.00565902, 0.0115686, 0.0215746, 0.0469424, 0.138205, 0.0512035, 0.0147657, 0.0190192, 0.00955465, 0.116809, 0.104003, 0.0383954, 0.0836653, 0.0819556, 0.00170794, 0.00621813],- [0.0495441, 0.0182506, 0.0143641, 0.0215379, 0.141805, 0.01402, 0.110854, 0.0247066, 0.0258142, 0.0700288, 0.0188272, 0.0315864, 0.0112101, 0.0316504, 0.0375346, 0.0456094, 0.0361428, 0.0369178, 0.0371985, 0.222397],- [0.170431, 0.000974733, 0.109856, 0.253646, 0.00133213, 0.0249846, 0.010139, 0.00587494, 0.0903324, 0.0116526, 0.00365127, 0.0271109, 0.0293614, 0.09173, 0.0415784, 0.0561766, 0.0479046, 0.02033, 0.000669682, 0.00226373],- [0.0162725, 0.00506141, 0.00101821, 0.00251413, 0.0376246, 0.00219354, 0.00299143, 0.132817, 0.00401204, 0.490444, 0.192993, 0.00218762, 0.00343332, 0.0104414, 0.00548261, 0.00401221, 0.0127074, 0.064772, 0.00321076, 0.00581006],- [0.0823766, 0.00656943, 0.0311745, 0.0675531, 0.00647179, 0.0178962, 0.0251144, 0.0291162, 0.0982302, 0.0287904, 0.0168023, 0.059839, 0.0114045, 0.0686451, 0.0734226, 0.1303, 0.182037, 0.0540271, 0.00227246, 0.00795733],- [0.359497, 0.0251417, 0.00314844, 0.00649627, 0.00920205, 0.119468, 0.00229704, 0.0458767, 0.00501688, 0.0468054, 0.0215569, 0.00334215, 0.0443916, 0.00490143, 0.00724072, 0.0465271, 0.0477755, 0.194216, 0.00245402, 0.00464504],- [0.201558, 0.00323653, 0.095415, 0.153491, 0.00256433, 0.0667292, 0.0155219, 0.00677408, 0.0547323, 0.0165114, 0.0060163, 0.0425386, 0.00919706, 0.0772011, 0.0430162, 0.118598, 0.0625473, 0.0202798, 0.000956551, 0.003115],- [0.104273, 0.00224501, 0.242407, 0.177482, 0.00125703, 0.169782, 0.0132649, 0.00189295, 0.0220652, 0.00425426, 0.00164412, 0.0621646, 0.0317042, 0.0356499, 0.0147062, 0.0778636, 0.0288516, 0.00602502, 0.00069309, 0.00177419],- [0.0781183, 0.0194449, 0.00415417, 0.0116634, 0.0262794, 0.0111524, 0.00635894, 0.135453, 0.00937298, 0.245757, 0.108778, 0.015927, 0.0055294, 0.0240152, 0.0111498, 0.0408519, 0.0860514, 0.148276, 0.00315476, 0.00851085],- [0.0856592, 0.0136073, 0.0135062, 0.00786026, 0.0047153, 0.0245401, 0.00553791, 0.0100592, 0.0127319, 0.0103344, 0.00806758, 0.0441923, 0.0175274, 0.00925906, 0.0101233, 0.340648, 0.357329, 0.019367, 0.00142431, 0.00350998],- [0.0674595, 0.00216342, 0.0662588, 0.0865501, 0.00182127, 0.0368557, 0.0381149, 0.00332388, 0.189974, 0.009384, 0.00394874, 0.116311, 0.0151208, 0.093936, 0.116173, 0.0842204, 0.0565954, 0.00645142, 0.00071873, 0.00461894],- [0.0572262, 0.00153015, 0.179393, 0.199226, 0.00137018, 0.0316472, 0.0291392, 0.00458046, 0.101562, 0.010074, 0.00402046, 0.108388, 0.00636741, 0.0903669, 0.0494724, 0.0621143, 0.0496102, 0.00859413, 0.000666929, 0.00464976],- [0.00360202, 0.00454848, 0.00208716, 0.00178577, 0.0855715, 0.00563916, 0.00649688, 0.00292929, 0.00104198, 0.0232635, 0.00445923, 0.00134555, 0.0024992, 0.00327181, 0.0102713, 0.00306718, 0.00259003, 0.00586684, 0.761782, 0.067881],- [0.203202, 0.00981316, 0.0135012, 0.00838182, 0.00196196, 0.618489, 0.00277479, 0.00118285, 0.00445989, 0.00398268, 0.00206318, 0.0143744, 0.00858704, 0.00445146, 0.00838957, 0.073992, 0.0108922, 0.00607691, 0.00186061, 0.00156387],- [0.00508988, 0.00617765, 0.00161262, 0.00120404, 0.356313, 0.00163342, 0.0393461, 0.00590888, 0.00137137, 0.0249344, 0.00497952, 0.0057093, 0.00141364, 0.00246931, 0.00287408, 0.00595277, 0.00365368, 0.00819341, 0.0357987, 0.485365],- [0.0403336, 0.00815495, 0.00982186, 0.0375407, 0.0119141, 0.00479344, 0.0176736, 0.189342, 0.0607377, 0.105186, 0.03056, 0.0216052, 0.00775506, 0.0383639, 0.0540186, 0.025711, 0.100991, 0.221091, 0.002878, 0.0115277],- [0.0790086, 0.00249479, 0.0546012, 0.199788, 0.00237734, 0.0192656, 0.02707, 0.00756675, 0.155311, 0.0254542, 0.00980244, 0.0309384, 0.00566407, 0.184338, 0.106544, 0.0332371, 0.0359575, 0.0145306, 0.00116828, 0.00488208],- [0.0722241, 0.00647553, 0.119488, 0.134589, 0.0348233, 0.0287815, 0.0699011, 0.0173589, 0.0490342, 0.051987, 0.0154411, 0.067893, 0.0145597, 0.070897, 0.0489728, 0.058958, 0.0425973, 0.0317884, 0.00879138, 0.0554387],- [0.108584, 0.00125026, 0.152967, 0.166485, 0.00145535, 0.0336098, 0.0134902, 0.00388218, 0.0576227, 0.00898614, 0.0024339, 0.0441956, 0.19901, 0.0405398, 0.020645, 0.084675, 0.0454715, 0.0113416, 0.000590283, 0.00276502],- [0.0309526, 0.0065594, 0.00823521, 0.0291974, 0.00368916, 0.0154206, 0.0310385, 0.00982516, 0.306263, 0.02379, 0.00970717, 0.0301337, 0.00950291, 0.0832608, 0.319589, 0.0295285, 0.0303052, 0.0133037, 0.00281253, 0.00688506],- [0.00989537, 0.00282767, 0.000374823, 0.00091821, 0.0298607, 0.000699707, 0.00104195, 0.311504, 0.00139605, 0.375039, 0.0474451, 0.000730793, 0.00252963, 0.0017337, 0.00196045, 0.0014628, 0.0075739, 0.1973, 0.00167998, 0.00402599],- [0.116321, 0.00347923, 0.0731918, 0.138088, 0.00941177, 0.0193193, 0.0160241, 0.0712243, 0.035512, 0.0771474, 0.0242841, 0.0250164, 0.0508927, 0.0586677, 0.0273321, 0.047556, 0.0726552, 0.123571, 0.00268927, 0.0076166],- [0.128522, 0.0172929, 0.040275, 0.0250692, 0.0199118, 0.112703, 0.0606981, 0.010935, 0.028875, 0.0258416, 0.0167593, 0.117984, 0.0180675, 0.0439706, 0.0373073, 0.174149, 0.0648968, 0.0182067, 0.0063575, 0.0321772],- [0.0372287, 0.014494, 0.00237032, 0.00485851, 0.146377, 0.00464339, 0.0186795, 0.182046, 0.00581985, 0.17801, 0.0371334, 0.00533773, 0.00485386, 0.00790971, 0.0094528, 0.0103571, 0.0284162, 0.184992, 0.0211294, 0.0958905],- [0.0535644, 0.00962562, 0.0113537, 0.0391699, 0.0264214, 0.0120279, 0.0384888, 0.0522748, 0.0996038, 0.189239, 0.0712219, 0.0239173, 0.00837206, 0.0928585, 0.11598, 0.0299114, 0.0389485, 0.0500948, 0.0104232, 0.026503],- [0.133242, 0.0246514, 0.00127392, 0.00404615, 0.0156995, 0.00891392, 0.00158647, 0.197128, 0.00237132, 0.125129, 0.0286947, 0.0022705, 0.0118846, 0.00308435, 0.00331477, 0.0171462, 0.0563298, 0.356621, 0.00173418, 0.00487784],- [0.149866, 0.00354374, 0.0280355, 0.0435381, 0.00475339, 0.0311113, 0.0140626, 0.0101953, 0.0393125, 0.0251434, 0.00515483, 0.0176453, 0.39238, 0.0348151, 0.0326607, 0.0874497, 0.0473307, 0.0271597, 0.00152152, 0.00432083],- [0.421437, 0.018761, 0.00733051, 0.00868378, 0.00271561, 0.0902333, 0.0030262, 0.00393628, 0.00515087, 0.00471933, 0.00383066, 0.012159, 0.020894, 0.00727486, 0.0061426, 0.290119, 0.0651922, 0.0252211, 0.000810824, 0.00236228],- [0.177071, 0.00783489, 0.02265, 0.0509767, 0.00405142, 0.089739, 0.0220667, 0.00595198, 0.125769, 0.020537, 0.00929825, 0.0311657, 0.0264088, 0.0752471, 0.133278, 0.116959, 0.0565567, 0.0165087, 0.00299471, 0.00493467],- [0.0293984, 0.00317291, 0.109971, 0.046427, 0.00150396, 0.422242, 0.0272495, 0.000799733, 0.0622314, 0.00376343, 0.00166571, 0.148362, 0.00564818, 0.0388688, 0.0370902, 0.0472252, 0.0086569, 0.00203639, 0.000917602, 0.00276931],- [0.0265779, 0.0101369, 0.0280314, 0.0269057, 0.0276961, 0.0173377, 0.281513, 0.0064647, 0.0474749, 0.026821, 0.00723753, 0.13186, 0.0083015, 0.0989711, 0.0791105, 0.0426277, 0.0259043, 0.0100147, 0.00785289, 0.0891598],- [0.00960965, 0.00773017, 0.000633186, 0.00104719, 0.263017, 0.00202274, 0.00390014, 0.0733098, 0.00149315, 0.445169, 0.0732575, 0.00131044, 0.00427681, 0.00338994, 0.00271362, 0.00361174, 0.00579284, 0.0425173, 0.0181276, 0.0370698]]---- | Weights of C40 model.-c40Weights :: [Double]-c40Weights = [0.0223854, 0.0338892, 0.0577169, 0.0252416, 0.0108608, 0.0462374, 0.0102293, 0.0147524, 0.0143161, 0.0182303, 0.0204025, 0.0425505, 0.0248627, 0.0105893, 0.0188239, 0.00866634, 0.0148496, 0.0343037, 0.0225335, 0.0174069, 0.0112208, 0.0443532, 0.0392573, 0.0196757, 0.028769, 0.0114441, 0.0112339, 0.0582694, 0.0444272, 0.0112011, 0.0145176, 0.0114629, 0.0239628, 0.0266266, 0.0481201, 0.0371147, 0.0160477, 0.0237249, 0.0235226, 0.0261998]---- | Stationary distribution of C50 model.-c50StatDists :: [StationaryDistribution]-c50StatDists = map V.fromList- [[0.115682, 0.00212412, 0.125699, 0.0919063, 0.0016175, 0.0391977, 0.0108596, 0.00344834, 0.0344439, 0.0125587, 0.00467132, 0.0530851, 0.206661, 0.0264814, 0.0147638, 0.154493, 0.0841319, 0.0134754, 0.000392551, 0.00430762],- [0.0983768, 0.00160659, 0.0456633, 0.0935028, 0.00112009, 0.0251388, 0.0191777, 0.00367483, 0.151407, 0.0115555, 0.00258422, 0.0266318, 0.229718, 0.0694385, 0.0855866, 0.0696565, 0.0498854, 0.0123983, 0.00069296, 0.00218381],- [0.0214598, 0.00446298, 0.000555089, 0.00135258, 0.00855351, 0.00112246, 0.000813505, 0.399631, 0.00128242, 0.0868331, 0.0170323, 0.00117488, 0.00183519, 0.00101872, 0.00142423, 0.00233511, 0.0161891, 0.430205, 0.00061736, 0.00210206],- [0.117155, 0.032654, 0.00620066, 0.00998032, 0.0657132, 0.0333659, 0.0134621, 0.0557335, 0.00794829, 0.199356, 0.0938018, 0.0221831, 0.0115415, 0.0204903, 0.0146873, 0.0842555, 0.0920314, 0.0798031, 0.0127314, 0.0269048],- [0.0371284, 0.0131607, 0.00255462, 0.00465054, 0.156922, 0.00416133, 0.019977, 0.178808, 0.00498614, 0.165159, 0.0352649, 0.0057612, 0.00512406, 0.00679862, 0.00826841, 0.0102949, 0.0275498, 0.181807, 0.0221255, 0.109498],- [0.0540087, 0.000836323, 0.306245, 0.38187, 0.00111469, 0.0245756, 0.0116353, 0.00254509, 0.030937, 0.00536416, 0.00256903, 0.0359123, 0.0122515, 0.0586791, 0.0115134, 0.0277199, 0.020314, 0.00828799, 0.000720693, 0.00290032],- [0.0404572, 0.028652, 0.00495693, 0.0068813, 0.242849, 0.0103569, 0.0318377, 0.0329142, 0.0105561, 0.143431, 0.0485779, 0.0105988, 0.00489299, 0.0118557, 0.0126873, 0.0229136, 0.0212457, 0.0349069, 0.0593874, 0.220042],- [0.319667, 0.138968, 0.00144847, 0.00200138, 0.00997903, 0.0503464, 0.00296588, 0.0182219, 0.00202126, 0.0260477, 0.0136422, 0.00936868, 0.00591365, 0.00358722, 0.00435407, 0.1656, 0.118413, 0.0996664, 0.00164228, 0.00614589],- [0.0977328, 0.00251083, 0.0724816, 0.197027, 0.00556062, 0.0117997, 0.0154288, 0.0563751, 0.0613631, 0.047368, 0.0159721, 0.0241106, 0.0149, 0.0909736, 0.034432, 0.0408092, 0.0929289, 0.111482, 0.00122744, 0.00551621],- [0.0195453, 0.0107627, 0.0127492, 0.011794, 0.14032, 0.0061666, 0.215761, 0.00826718, 0.0118897, 0.0328258, 0.00664433, 0.0423533, 0.00528472, 0.0258824, 0.024497, 0.0203711, 0.0124449, 0.0133889, 0.0330587, 0.345993],- [0.185338, 0.00185224, 0.161345, 0.234878, 0.00192286, 0.054278, 0.0109591, 0.00550122, 0.0382027, 0.00925008, 0.00343293, 0.033406, 0.0233824, 0.0551124, 0.0201268, 0.0884267, 0.051789, 0.0177494, 0.000835994, 0.00221142],- [0.0143603, 0.0036047, 0.000411253, 0.00108887, 0.0196479, 0.00187385, 0.000991664, 0.285685, 0.00213672, 0.317288, 0.150556, 0.00155703, 0.00122034, 0.0028111, 0.00194519, 0.00344669, 0.0196686, 0.167294, 0.000843478, 0.00356936],- [0.150577, 0.00957632, 0.0227565, 0.0358641, 0.0104881, 0.0285062, 0.0171219, 0.0577109, 0.0276149, 0.0531781, 0.0298859, 0.0514957, 0.0243637, 0.0354202, 0.0272841, 0.132633, 0.169554, 0.104387, 0.00219421, 0.00938724],- [0.277615, 0.0280749, 0.00194439, 0.00537302, 0.0070772, 0.0234861, 0.0012257, 0.0949678, 0.00308357, 0.0474515, 0.0189616, 0.0034483, 0.0134815, 0.00374076, 0.00354151, 0.0423818, 0.0932708, 0.326207, 0.000963513, 0.00370386],- [0.205879, 0.00249091, 0.075426, 0.161651, 0.00210732, 0.0463391, 0.0159993, 0.00515537, 0.0913036, 0.0166911, 0.00536596, 0.0366042, 0.00849072, 0.0973331, 0.0696431, 0.088244, 0.0501877, 0.0172747, 0.000776174, 0.0030381],- [0.0704473, 0.00149732, 0.187469, 0.208434, 0.0013217, 0.0409359, 0.0247326, 0.00339884, 0.0739754, 0.00761491, 0.00270559, 0.0965448, 0.00673196, 0.0820493, 0.0355811, 0.0832004, 0.0610002, 0.00788473, 0.000646918, 0.00382803],- [0.108767, 0.00735776, 0.0100122, 0.0249269, 0.015562, 0.0114116, 0.00730171, 0.0923916, 0.0138538, 0.112308, 0.0251594, 0.00831075, 0.259979, 0.0187962, 0.0142911, 0.0329858, 0.0608007, 0.162244, 0.0049709, 0.00857003],- [0.109878, 0.0295589, 0.000867374, 0.00267084, 0.028604, 0.0101436, 0.00159952, 0.171472, 0.0015037, 0.241889, 0.0488425, 0.00201796, 0.00631863, 0.00393786, 0.00336164, 0.0148413, 0.0393773, 0.273667, 0.00359303, 0.00585493],- [0.058736, 0.00692528, 0.0144651, 0.0412403, 0.0154694, 0.015436, 0.0243546, 0.0388984, 0.0996366, 0.207192, 0.103832, 0.0289715, 0.00833815, 0.108218, 0.0976056, 0.0350037, 0.0432645, 0.0338957, 0.00652556, 0.0119916],- [0.0642736, 0.00525071, 0.144305, 0.153274, 0.0286093, 0.0265312, 0.0545105, 0.0231701, 0.0452879, 0.058794, 0.0192273, 0.071063, 0.0146855, 0.066931, 0.0405954, 0.054046, 0.043602, 0.0350962, 0.00655203, 0.0441953],- [0.00487716, 0.00503139, 0.00209095, 0.00192293, 0.0983777, 0.00615619, 0.00725617, 0.00400789, 0.000993075, 0.0277171, 0.00537714, 0.00145717, 0.0029039, 0.00350461, 0.0112195, 0.00329284, 0.00295187, 0.00708726, 0.724652, 0.0791231],- [0.0552756, 0.0115894, 0.00846155, 0.0432312, 0.00961795, 0.0074182, 0.0218432, 0.0843826, 0.163327, 0.0631707, 0.0228359, 0.0221225, 0.00652522, 0.0639538, 0.146077, 0.0387464, 0.100713, 0.118118, 0.00298298, 0.00960875],- [0.00449551, 0.00619902, 0.000890609, 0.000674178, 0.400884, 0.0013356, 0.0223834, 0.00805107, 0.00105813, 0.0341606, 0.00701178, 0.00318917, 0.00130624, 0.00157564, 0.00181707, 0.00502515, 0.00337803, 0.00946203, 0.0345278, 0.452575],- [0.0505469, 0.00768183, 0.0101983, 0.0470721, 0.0160596, 0.00517598, 0.0176387, 0.169463, 0.0497501, 0.209444, 0.0474166, 0.0132404, 0.00844474, 0.053753, 0.051076, 0.0149657, 0.03846, 0.172926, 0.00437246, 0.0123139],- [0.0749936, 0.0124353, 0.0203869, 0.0100161, 0.0038092, 0.0270453, 0.00581932, 0.00743392, 0.014732, 0.00838354, 0.00640172, 0.0530316, 0.0176075, 0.0096947, 0.0117746, 0.3482, 0.348235, 0.0157071, 0.00137596, 0.00291601],- [0.0569785, 0.00176591, 0.102737, 0.106157, 0.00132179, 0.0360979, 0.034465, 0.00311879, 0.185459, 0.00858711, 0.00431374, 0.139374, 0.00855798, 0.0898676, 0.101042, 0.0669467, 0.0422839, 0.00552428, 0.000686497, 0.00471478],- [0.0746973, 0.00394557, 0.0423454, 0.0906599, 0.00363883, 0.0183504, 0.0229406, 0.0156166, 0.12297, 0.0167059, 0.00901115, 0.0672502, 0.0127747, 0.0811334, 0.0755385, 0.132886, 0.170553, 0.0324157, 0.00152093, 0.00504616],- [0.17313, 0.00938612, 0.0203504, 0.0387504, 0.00448721, 0.0867075, 0.0227177, 0.00492238, 0.125544, 0.0156727, 0.00787132, 0.0400476, 0.0180405, 0.0680443, 0.133335, 0.142136, 0.0670614, 0.0136107, 0.00287027, 0.00531475],- [0.00751784, 0.00212117, 0.000341917, 0.000667876, 0.0310677, 0.000466913, 0.000898937, 0.320302, 0.00103903, 0.384419, 0.0335748, 0.000602706, 0.00243131, 0.00137426, 0.00152778, 0.00115498, 0.00552329, 0.199295, 0.00172716, 0.00394622],- [0.113459, 0.0281671, 0.0268294, 0.0283293, 0.0928796, 0.0494616, 0.0733971, 0.0152335, 0.0176739, 0.037425, 0.0140126, 0.0539481, 0.0259048, 0.0303861, 0.0259068, 0.121621, 0.0621544, 0.0350518, 0.0240853, 0.124074],- [0.0810974, 0.00899714, 0.0930492, 0.0372909, 0.002713, 0.0907986, 0.0275326, 0.00533957, 0.0362449, 0.00831441, 0.00615441, 0.177093, 0.0120848, 0.0316411, 0.025247, 0.233557, 0.106302, 0.00814021, 0.000950692, 0.00745168],- [0.00757098, 0.00716117, 0.00049456, 0.000891326, 0.297029, 0.00155605, 0.00342729, 0.0711141, 0.00107067, 0.416918, 0.0635143, 0.000864035, 0.00349893, 0.00251506, 0.00186372, 0.00271598, 0.00424124, 0.041488, 0.0202374, 0.0518279],- [0.086899, 0.00542266, 0.0274343, 0.117345, 0.00568342, 0.0162895, 0.0386176, 0.0182492, 0.149526, 0.048245, 0.0150062, 0.0265186, 0.00577334, 0.159013, 0.149227, 0.0416278, 0.0492401, 0.0283951, 0.00272523, 0.00876225],- [0.0501736, 0.0113806, 0.00768113, 0.0191841, 0.0119036, 0.00521331, 0.00988204, 0.203482, 0.0206446, 0.0870652, 0.036001, 0.0197961, 0.0061624, 0.0180415, 0.0178215, 0.0399964, 0.162255, 0.262373, 0.00172457, 0.00921833],- [0.177098, 0.00357587, 0.0740854, 0.0713923, 0.00385864, 0.211235, 0.0161904, 0.0040472, 0.0370753, 0.00826202, 0.00369757, 0.0494737, 0.0543871, 0.0396583, 0.0340392, 0.140615, 0.0532238, 0.0108131, 0.00212448, 0.00514764],- [0.242384, 0.0115188, 0.00720908, 0.00544261, 0.00229719, 0.593273, 0.00205723, 0.00208173, 0.00287156, 0.00491, 0.00271204, 0.00899551, 0.00763059, 0.00285962, 0.00600369, 0.0738385, 0.0108329, 0.00968775, 0.00178961, 0.00160391],- [0.032912, 0.00319267, 0.143899, 0.0532755, 0.00166462, 0.463509, 0.0175283, 0.000876108, 0.0317056, 0.00297349, 0.0012156, 0.131603, 0.00688601, 0.0251373, 0.0184844, 0.0491059, 0.00984294, 0.00239178, 0.0008917, 0.00290456],- [0.0388306, 0.00489943, 0.0491655, 0.0413874, 0.00170354, 0.345044, 0.0421499, 0.00103798, 0.123725, 0.00726049, 0.00335164, 0.119727, 0.00562115, 0.0664068, 0.0881106, 0.0460192, 0.00791813, 0.0020819, 0.00195492, 0.00360516],- [0.0373301, 0.00916334, 0.0396965, 0.0362168, 0.0144959, 0.0269081, 0.220512, 0.00638711, 0.0634517, 0.0225147, 0.00718608, 0.152708, 0.00995802, 0.108449, 0.0882553, 0.0609295, 0.0332848, 0.0105312, 0.00429315, 0.047729],- [0.0171508, 0.00565843, 0.00111564, 0.00268175, 0.0490847, 0.00219974, 0.00370226, 0.116938, 0.00407925, 0.500942, 0.1832, 0.00280387, 0.00400697, 0.0115255, 0.00613355, 0.00539456, 0.0140434, 0.0579479, 0.00419195, 0.00719895],- [0.0238117, 0.00329378, 0.370777, 0.0831676, 0.0022275, 0.053433, 0.0290274, 0.00288609, 0.0258098, 0.00370314, 0.00159506, 0.259554, 0.00668902, 0.0201877, 0.0117854, 0.0626213, 0.0278627, 0.00464381, 0.000688992, 0.0062349],- [0.209736, 0.0105937, 0.0311812, 0.0613839, 0.00991105, 0.0855144, 0.0226071, 0.0296058, 0.0455281, 0.0792591, 0.0306389, 0.019447, 0.0213982, 0.0610672, 0.0588191, 0.0820738, 0.0527883, 0.0722389, 0.00513, 0.0110785],- [0.0295833, 0.00588894, 0.00852018, 0.0296344, 0.00323424, 0.0154495, 0.02992, 0.00752729, 0.317117, 0.0212988, 0.00874154, 0.0312424, 0.00843308, 0.0831654, 0.323587, 0.0294177, 0.0283085, 0.0100345, 0.00257807, 0.00631797],- [0.0937098, 0.00289961, 0.0255803, 0.0577481, 0.0140839, 0.0182984, 0.0300827, 0.0312422, 0.0867966, 0.067914, 0.01767, 0.0206426, 0.213942, 0.0611154, 0.0902988, 0.0526724, 0.0464314, 0.0509513, 0.00365078, 0.01427],- [0.419507, 0.0176715, 0.00719436, 0.00880545, 0.00285112, 0.0865515, 0.00330323, 0.00419601, 0.00503365, 0.0048203, 0.00393323, 0.012443, 0.0192167, 0.00750218, 0.00609116, 0.295973, 0.0660748, 0.0253174, 0.000942289, 0.00257126],- [0.0368106, 0.0124259, 0.0102859, 0.0343205, 0.0759532, 0.00772969, 0.0938663, 0.0458547, 0.0701169, 0.133403, 0.031916, 0.0265402, 0.00749485, 0.0673302, 0.10388, 0.0271625, 0.0353789, 0.0485031, 0.0254607, 0.105568],- [0.0665479, 0.00185956, 0.183575, 0.149215, 0.00127224, 0.0634659, 0.019684, 0.00262142, 0.0655694, 0.00750609, 0.00151282, 0.0699101, 0.182105, 0.0446839, 0.0292035, 0.075452, 0.0278833, 0.00484561, 0.000723776, 0.00236318],- [0.16002, 0.00397537, 0.019682, 0.0286161, 0.00458768, 0.0326662, 0.0114325, 0.00722291, 0.0232558, 0.0202542, 0.00379358, 0.0142832, 0.460797, 0.0260458, 0.0217875, 0.0913633, 0.0427326, 0.0227118, 0.0012566, 0.0035159],- [0.0843584, 0.00129715, 0.0847585, 0.260552, 0.00143564, 0.0191268, 0.0174708, 0.00556223, 0.144172, 0.014456, 0.00550323, 0.0336163, 0.00724868, 0.160704, 0.0742581, 0.0323411, 0.0355032, 0.0136662, 0.000760206, 0.00320994],- [0.152149, 0.000925399, 0.116471, 0.220438, 0.00150071, 0.0172039, 0.0095692, 0.00747151, 0.0529323, 0.0148756, 0.00356173, 0.0157411, 0.217335, 0.0489135, 0.0171707, 0.044102, 0.0323715, 0.023915, 0.000908902, 0.00244332]]---- | Weights of C50 model.-c50Weights :: [Double]-c50Weights = [0.00509464, 0.00753314, 0.0595801, 0.0100458, 0.0109109, 0.0427412, 0.00632106, 0.0192362, 0.0119168, 0.0154209, 0.0148037, 0.014293, 0.0164064, 0.0272538, 0.0178737, 0.0125183, 0.00698661, 0.0140027, 0.0140399, 0.0109043, 0.011579, 0.0138661, 0.0377329, 0.0118586, 0.0356688, 0.0161327, 0.0200601, 0.0133087, 0.0357641, 0.0119353, 0.0241183, 0.0243098, 0.0178039, 0.0200988, 0.00579569, 0.0397441, 0.0240799, 0.00868054, 0.0239822, 0.042803, 0.0335193, 0.00636191, 0.0547546, 0.00544598, 0.0354219, 0.0120158, 0.0047506, 0.036392, 0.0270346, 0.00709747]---- | Stationary distribution of C60 model.-c60StatDists :: [StationaryDistribution]-c60StatDists = map V.fromList- [[0.153436, 0.00473066, 0.0546757, 0.0529324, 0.00355543, 0.110378, 0.017448, 0.00503439, 0.0603929, 0.0130294, 0.00755506, 0.0796727, 0.0249524, 0.0514333, 0.0444389, 0.202963, 0.0957668, 0.0101809, 0.00144445, 0.00598003],- [0.0281985, 0.00195035, 0.00915493, 0.038853, 0.00159193, 0.00840283, 0.0155385, 0.0107873, 0.347604, 0.0217787, 0.0109905, 0.0312955, 0.00715399, 0.0939884, 0.303106, 0.0197479, 0.0328352, 0.0135303, 0.000920999, 0.0025714],- [0.00836807, 0.00204339, 0.000222837, 0.000473154, 0.0123984, 0.000482575, 0.000518999, 0.376845, 0.000629617, 0.260833, 0.0315701, 0.000612345, 0.000959592, 0.000949869, 0.000731977, 0.000974638, 0.00499908, 0.293808, 0.000865776, 0.00171323],- [0.222723, 0.164029, 0.000799759, 0.00272808, 0.0228667, 0.0202703, 0.00371839, 0.0455406, 0.00228329, 0.088335, 0.0348095, 0.00612065, 0.00354716, 0.00510519, 0.00648461, 0.085004, 0.101285, 0.169858, 0.00484248, 0.00965],- [0.041214, 0.0204716, 0.00174344, 0.00372634, 0.2239, 0.00691075, 0.0180294, 0.115428, 0.00429003, 0.169356, 0.0414067, 0.00510677, 0.00584164, 0.00575385, 0.00676271, 0.0149106, 0.0239548, 0.120034, 0.0332237, 0.137935],- [0.048055, 0.00378018, 0.208461, 0.185518, 0.0196001, 0.0191402, 0.0872162, 0.00946854, 0.0375741, 0.0277341, 0.00883084, 0.0940629, 0.00812678, 0.0747677, 0.0308438, 0.0439681, 0.0324589, 0.0181769, 0.00346657, 0.03875],- [0.00628487, 0.00730276, 0.000532415, 0.000974149, 0.397455, 0.00195782, 0.0102225, 0.0180592, 0.00162059, 0.117906, 0.0234722, 0.00303425, 0.00201656, 0.00344091, 0.00262469, 0.00569033, 0.00370918, 0.0152754, 0.0598639, 0.318557],- [0.181501, 0.420563, 0.00251455, 0.000708814, 0.00689764, 0.0461854, 0.00143746, 0.00417455, 0.000647425, 0.00983105, 0.00416114, 0.0148485, 0.00387672, 0.0014097, 0.00268454, 0.186454, 0.068719, 0.0345742, 0.00270835, 0.0061033],- [0.0600741, 0.00707198, 0.0170572, 0.0127487, 0.114091, 0.0343564, 0.0305859, 0.0204571, 0.0212367, 0.0994551, 0.0318166, 0.0134869, 0.0297628, 0.014247, 0.0367643, 0.0505793, 0.0339368, 0.0217422, 0.231281, 0.119249],- [0.0708395, 0.0131641, 0.0324483, 0.0508264, 0.059436, 0.0183309, 0.0567273, 0.0650369, 0.0343618, 0.128226, 0.0390363, 0.0416822, 0.0135608, 0.0494875, 0.0474098, 0.0551343, 0.064226, 0.078961, 0.0137118, 0.0673934],- [0.0617689, 0.000719984, 0.343023, 0.379251, 0.000981419, 0.0284659, 0.00795921, 0.00169996, 0.0216077, 0.00399453, 0.00197343, 0.0303082, 0.0174791, 0.0307856, 0.00763329, 0.0337832, 0.0203427, 0.00589923, 0.000613027, 0.00171028],- [0.042156, 0.00444503, 0.0164385, 0.0531949, 0.00867283, 0.0134711, 0.0177765, 0.0267728, 0.132374, 0.196724, 0.118283, 0.0286981, 0.00518379, 0.139369, 0.104207, 0.0255853, 0.0333292, 0.0217067, 0.00458523, 0.00702815],- [0.281481, 0.0258338, 0.00643857, 0.0115046, 0.00504817, 0.0270055, 0.00546297, 0.0188216, 0.00987128, 0.0190993, 0.015872, 0.0172868, 0.012951, 0.0133102, 0.0100367, 0.188681, 0.24271, 0.0837524, 0.00122746, 0.00360529],- [0.276919, 0.0318293, 0.00116725, 0.00247525, 0.00708439, 0.0199647, 0.000517086, 0.0983109, 0.00162322, 0.0489264, 0.0173415, 0.00213153, 0.00701797, 0.00182163, 0.0017227, 0.0336349, 0.0814141, 0.362092, 0.000711814, 0.00329423],- [0.15778, 0.00533181, 0.064829, 0.0913587, 0.00250478, 0.090621, 0.053381, 0.00291716, 0.0782149, 0.0156226, 0.0045759, 0.0570403, 0.00670778, 0.106537, 0.111214, 0.092931, 0.0393123, 0.00850409, 0.00285758, 0.00775903],- [0.0593735, 0.00153322, 0.218948, 0.181942, 0.00108383, 0.0813707, 0.0220478, 0.00209936, 0.0750172, 0.00561913, 0.00218717, 0.115118, 0.0109737, 0.0688752, 0.0354741, 0.0726449, 0.0380238, 0.00426697, 0.000734646, 0.00266649],- [0.0978066, 0.00641384, 0.0120781, 0.0291283, 0.0198964, 0.0128913, 0.0128037, 0.0723904, 0.0245631, 0.127944, 0.0303267, 0.0101844, 0.272352, 0.0307876, 0.0265576, 0.0350549, 0.0484557, 0.115203, 0.00468425, 0.0104774],- [0.0124023, 0.00418857, 0.000603732, 0.0012524, 0.0610615, 0.00113068, 0.00135421, 0.234417, 0.00167182, 0.455056, 0.0667308, 0.00092391, 0.00370762, 0.00329574, 0.00306804, 0.00194209, 0.00676129, 0.129048, 0.00389372, 0.00749118],- [0.0794231, 0.00199195, 0.0587236, 0.12469, 0.001287, 0.0325432, 0.0238605, 0.00362772, 0.21476, 0.0097988, 0.00418462, 0.0662793, 0.014241, 0.114388, 0.129474, 0.0615807, 0.0477334, 0.00781872, 0.000652537, 0.00294202],- [0.0248149, 0.00278397, 0.427883, 0.143239, 0.00146804, 0.0643968, 0.0185737, 0.00225069, 0.0179274, 0.00345586, 0.00157145, 0.188892, 0.0073768, 0.0210778, 0.00835529, 0.0377003, 0.0187768, 0.00460727, 0.000589186, 0.00426027],- [0.00170034, 0.00290594, 0.00107115, 0.00117317, 0.0513377, 0.00355796, 0.00276305, 0.00122912, 0.00042738, 0.0127421, 0.00256713, 0.00042229, 0.00135367, 0.00164242, 0.00606743, 0.00118717, 0.00140331, 0.00331373, 0.864044, 0.0390913],- [0.0468361, 0.00599541, 0.0185616, 0.0705436, 0.00801772, 0.00454353, 0.0152063, 0.155061, 0.0866248, 0.0824253, 0.0245854, 0.0205604, 0.00814856, 0.0557031, 0.0639797, 0.0237026, 0.0962054, 0.204749, 0.00183687, 0.00671317],- [0.0258765, 0.0142504, 0.0107037, 0.010565, 0.187643, 0.0073065, 0.141108, 0.0114778, 0.0119179, 0.0407992, 0.0098799, 0.0298384, 0.00512288, 0.0158529, 0.0201097, 0.02757, 0.0170765, 0.0187502, 0.0405125, 0.353639],- [0.0296285, 0.00900465, 0.000885102, 0.00330465, 0.0408798, 0.00279695, 0.00507015, 0.25834, 0.00466909, 0.266809, 0.0770825, 0.00349444, 0.00269182, 0.00554811, 0.00464003, 0.00685381, 0.0322266, 0.22719, 0.00355061, 0.0153353],- [0.0555726, 0.00974186, 0.014039, 0.00694432, 0.00344967, 0.0157957, 0.00416313, 0.00698265, 0.0139225, 0.00752712, 0.00587627, 0.0409064, 0.0119733, 0.00687277, 0.00984479, 0.348247, 0.421366, 0.0131119, 0.00100615, 0.00265768],- [0.0682671, 0.0141433, 0.0360279, 0.0497415, 0.0803376, 0.0268697, 0.112767, 0.0132647, 0.0482609, 0.0544494, 0.0170034, 0.0530661, 0.019195, 0.0612274, 0.0615207, 0.067184, 0.0443996, 0.0267714, 0.0199958, 0.125507],- [0.0792619, 0.00495033, 0.0436646, 0.0829639, 0.00491311, 0.0183429, 0.0233169, 0.0249427, 0.0932578, 0.0221483, 0.0120893, 0.0635289, 0.012636, 0.0666365, 0.0638377, 0.133485, 0.191675, 0.0503102, 0.00180401, 0.00623531],- [0.0731759, 0.0123459, 0.0110007, 0.0349943, 0.00931488, 0.0224509, 0.0386903, 0.0246226, 0.178334, 0.0508307, 0.0185741, 0.0324201, 0.0148723, 0.0858952, 0.210534, 0.0603181, 0.0649575, 0.038504, 0.00510464, 0.0130597],- [0.0878403, 0.0171147, 0.00328039, 0.0101466, 0.0191954, 0.00870909, 0.00549022, 0.198776, 0.00820969, 0.175646, 0.0417233, 0.00608012, 0.0111284, 0.0109832, 0.0110332, 0.0209863, 0.0697719, 0.28252, 0.0031744, 0.00819055],- [0.0990216, 0.0108996, 0.0118798, 0.0279, 0.0214442, 0.024025, 0.0123445, 0.0776565, 0.0322818, 0.235451, 0.120767, 0.021115, 0.00756555, 0.0557711, 0.0349352, 0.052417, 0.0649785, 0.0771067, 0.00470758, 0.00773287],- [0.0601641, 0.0315067, 0.0337189, 0.00598395, 0.00695259, 0.0543081, 0.0531524, 0.00706508, 0.0139598, 0.00706981, 0.00882987, 0.278352, 0.00758343, 0.0210646, 0.0161995, 0.28298, 0.0860317, 0.0100954, 0.00149666, 0.0134849],- [0.00497817, 0.00529787, 0.000531537, 0.000763055, 0.335112, 0.00150513, 0.00415752, 0.0443829, 0.00116151, 0.441742, 0.0602807, 0.000701221, 0.00278477, 0.00240899, 0.00184123, 0.00257958, 0.00302885, 0.0237677, 0.0171303, 0.0458456],- [0.0251997, 0.00970995, 0.00410123, 0.0262572, 0.0429092, 0.00380673, 0.0431939, 0.0974043, 0.0824313, 0.24472, 0.0539323, 0.0142032, 0.00526585, 0.0620071, 0.111447, 0.00960931, 0.0251183, 0.0783303, 0.0146572, 0.0456965],- [0.0230362, 0.00485801, 0.000626444, 0.00157626, 0.0072127, 0.00115653, 0.000889924, 0.39529, 0.00146635, 0.0576537, 0.0140987, 0.00135344, 0.00201779, 0.000987005, 0.00147487, 0.00287702, 0.0205581, 0.460366, 0.000547769, 0.00195391],- [0.140878, 0.00322137, 0.0285077, 0.0423839, 0.00412311, 0.0299558, 0.0131321, 0.00669145, 0.0383782, 0.019512, 0.00362769, 0.0171298, 0.438347, 0.0320632, 0.0297808, 0.08514, 0.0422766, 0.0198195, 0.00132369, 0.00370876],- [0.444249, 0.0235405, 0.0057193, 0.00648899, 0.00244132, 0.0822688, 0.00255057, 0.00336151, 0.00380971, 0.00409901, 0.00286831, 0.0102305, 0.0162891, 0.00559976, 0.00502166, 0.299997, 0.0559665, 0.0226608, 0.000773543, 0.00206398],- [0.0898718, 0.00441315, 0.0129167, 0.0058957, 0.000807357, 0.808756, 0.00165179, 0.000433928, 0.0033939, 0.00155645, 0.000425342, 0.0130068, 0.00341281, 0.00238065, 0.00709583, 0.0362877, 0.00328875, 0.00202896, 0.00152239, 0.000853745],- [0.055084, 0.00358249, 0.12738, 0.0655885, 0.0040375, 0.0637608, 0.0404883, 0.00755742, 0.0867683, 0.0136305, 0.00816842, 0.187783, 0.0110682, 0.0527969, 0.0472254, 0.126338, 0.0752544, 0.0102728, 0.00135637, 0.011859],- [0.0117681, 0.00703883, 0.0144712, 0.0110901, 0.0338266, 0.00496266, 0.433719, 0.0061337, 0.0137929, 0.0298795, 0.00762376, 0.0844145, 0.00813461, 0.103834, 0.0442559, 0.0140571, 0.0108277, 0.00837578, 0.00806831, 0.143725],- [0.0159286, 0.00805199, 0.00109669, 0.00213242, 0.0743651, 0.00250864, 0.00491921, 0.0501544, 0.00355994, 0.530763, 0.216009, 0.00326926, 0.00452474, 0.0134552, 0.00480987, 0.00669222, 0.0119092, 0.0281012, 0.00709281, 0.0106565],- [0.0195973, 0.0034539, 0.309985, 0.0250778, 0.00197065, 0.0627529, 0.0295961, 0.00326504, 0.0240964, 0.00282469, 0.000842506, 0.32891, 0.00490628, 0.0116197, 0.0105143, 0.106498, 0.0438054, 0.0040302, 0.000633396, 0.0056198],- [0.0833804, 0.00819811, 0.0686821, 0.0227415, 0.00241326, 0.0982292, 0.00739549, 0.00174712, 0.0129342, 0.00396531, 0.0019558, 0.0969824, 0.0355924, 0.0121521, 0.0125871, 0.311561, 0.211337, 0.00471386, 0.0016329, 0.00179911],- [0.0181409, 0.00851436, 0.00333335, 0.00965319, 0.00354269, 0.0224553, 0.0375239, 0.00350905, 0.320807, 0.0149147, 0.00460989, 0.0233205, 0.00871975, 0.0526694, 0.412966, 0.0262309, 0.0131791, 0.00447465, 0.0034767, 0.00795882],- [0.249423, 0.0122863, 0.0169234, 0.0370491, 0.00773609, 0.0347468, 0.00870696, 0.0595718, 0.0177065, 0.0451065, 0.0204556, 0.016412, 0.0686404, 0.0228502, 0.0185482, 0.0889296, 0.0986017, 0.167252, 0.00286039, 0.00619385],- [0.141974, 0.00688567, 0.0537452, 0.040754, 0.00402861, 0.398811, 0.0152896, 0.00166276, 0.0314274, 0.00923483, 0.00554255, 0.0576297, 0.0180329, 0.028624, 0.0373946, 0.112373, 0.0242478, 0.00484629, 0.00259091, 0.00490542],- [0.0178903, 0.0290304, 0.00310543, 0.00405034, 0.0238913, 0.010048, 0.12527, 0.0016708, 0.023334, 0.0722357, 0.0116142, 0.0155854, 0.000993842, 0.105182, 0.195884, 0.0181676, 0.018626, 0.0113865, 0.226055, 0.0859787],- [0.145476, 0.00141863, 0.129425, 0.247135, 0.00175521, 0.0319651, 0.0130426, 0.00585254, 0.0818154, 0.0123593, 0.00441789, 0.040042, 0.0151136, 0.0906469, 0.0420979, 0.0656688, 0.0511289, 0.0169401, 0.000773144, 0.00292584],- [0.0169799, 0.00220872, 0.104366, 0.037638, 0.000678059, 0.540447, 0.0181137, 0.000745968, 0.0508618, 0.0021147, 0.000947377, 0.131805, 0.00387544, 0.026935, 0.0242347, 0.029703, 0.00458362, 0.00158917, 0.000603189, 0.00157041],- [0.0402646, 0.00393887, 0.0293968, 0.160352, 0.00411534, 0.012926, 0.0617968, 0.00984913, 0.14488, 0.0354919, 0.0124819, 0.0323829, 0.00433742, 0.249701, 0.115202, 0.0243247, 0.0305645, 0.0156027, 0.00266766, 0.00972278],- [0.225691, 0.00433901, 0.0637125, 0.115483, 0.00274401, 0.0867335, 0.0131067, 0.00850862, 0.0660469, 0.0193314, 0.00648772, 0.0244309, 0.0611149, 0.0578159, 0.0523417, 0.107088, 0.0507677, 0.0299883, 0.00136959, 0.00289829],- [0.00331642, 0.00518907, 0.000826647, 0.000513413, 0.380176, 0.00104328, 0.0278451, 0.00418953, 0.000714992, 0.0111212, 0.00236218, 0.00308302, 0.000836508, 0.00110249, 0.00153108, 0.00358767, 0.00236089, 0.00607663, 0.0333347, 0.510789],- [0.199501, 0.058459, 0.00949551, 0.0125535, 0.0634434, 0.0939787, 0.0341857, 0.014021, 0.0118723, 0.0449388, 0.0246991, 0.0392255, 0.0145385, 0.0254265, 0.0236079, 0.166392, 0.0533159, 0.0367303, 0.0129803, 0.0606346],- [0.0319449, 0.00660403, 0.0378075, 0.0279488, 0.0108916, 0.0365541, 0.208864, 0.00475423, 0.0868665, 0.0156546, 0.00432533, 0.208471, 0.00608996, 0.0766373, 0.101167, 0.0577657, 0.0302051, 0.00605852, 0.00260019, 0.0387897],- [0.07768, 0.0140093, 0.00806517, 0.0112622, 0.0173585, 0.0136868, 0.013373, 0.123964, 0.0129145, 0.0724671, 0.0420745, 0.0403217, 0.0117084, 0.0179996, 0.0142519, 0.0922724, 0.23169, 0.167914, 0.00281536, 0.0141727],- [0.118366, 0.00356248, 0.0495595, 0.153759, 0.00686657, 0.0194993, 0.0229373, 0.0302661, 0.0982304, 0.0571251, 0.0171727, 0.0259525, 0.0175153, 0.120492, 0.0805193, 0.0486588, 0.0635796, 0.0553336, 0.00230083, 0.00830274],- [0.052856, 0.000862588, 0.209276, 0.402419, 0.00126507, 0.0155838, 0.0124149, 0.00548648, 0.0497017, 0.00902565, 0.00423571, 0.0264744, 0.00631856, 0.121241, 0.0193569, 0.0197263, 0.0235464, 0.0159764, 0.000838161, 0.00339487],- [0.0344366, 0.0020986, 0.113901, 0.054178, 0.00315576, 0.336164, 0.0461777, 0.000346342, 0.0667553, 0.00483557, 0.00197045, 0.163672, 0.00403698, 0.0605414, 0.0426222, 0.0481173, 0.00891481, 0.00201106, 0.00065101, 0.00541456],- [0.115309, 0.00149624, 0.175552, 0.174941, 0.00157002, 0.0394181, 0.0132402, 0.0056913, 0.0433118, 0.010141, 0.00303321, 0.0458477, 0.16658, 0.0366731, 0.0151279, 0.0871536, 0.0468261, 0.0141932, 0.00075157, 0.00314327],- [0.386515, 0.0218928, 0.00223668, 0.00313878, 0.00392814, 0.369435, 0.00146729, 0.00853761, 0.00188405, 0.0127257, 0.00805817, 0.00304205, 0.0158688, 0.00214647, 0.00375793, 0.0808877, 0.0305196, 0.0410635, 0.000992288, 0.00190203],- [0.0146571, 0.00243179, 0.000521057, 0.00148744, 0.00916375, 0.00209533, 0.00101819, 0.19139, 0.00228984, 0.44328, 0.221743, 0.00129983, 0.000768515, 0.00493628, 0.00288413, 0.00272515, 0.0170997, 0.0778557, 0.000877938, 0.0014756]]---- | Weights of C60 model.-c60Weights :: [Double]-c60Weights = [0.0169699, 0.0211683, 0.0276589, 0.0065676, 0.0141221, 0.00687748, 0.014691, 0.00672258, 0.00183967, 0.0102547, 0.0230896, 0.0057941, 0.0125395, 0.0204526, 0.00706296, 0.0117983, 0.00683347, 0.0433776, 0.0318279, 0.0222546, 0.0102265, 0.0150546, 0.013416, 0.0148552, 0.0239112, 0.0128776, 0.0222319, 0.0247445, 0.0214275, 0.0115002, 0.00760174, 0.0130259, 0.0093702, 0.0467194, 0.044194, 0.0322263, 0.0403, 0.0150234, 0.010459, 0.0214742, 0.0154958, 0.010179, 0.022798, 0.0123205, 0.00667776, 0.000415008, 0.0344385, 0.0113663, 0.0127143, 0.0124324, 0.0262124, 0.0064995, 0.0103203, 0.0142464, 0.02156, 0.0199151, 0.00389642, 0.0113449, 0.0128596, 0.0117657]
− src/ELynx/Data/MarkovProcess/EDMModel.hs
@@ -1,38 +0,0 @@-{- |-Module : ELynx.Data.MarkovProcess.EDMModel-Description : Empiricial distribution mixture models-Copyright : (c) Dominik Schrempf 2019-License : GPL-3--Maintainer : dominik.schrempf@gmail.com-Stability : unstable-Portability : portable--Creation date: Tue Mar 19 13:25:55 2019.---}--module ELynx.Data.MarkovProcess.EDMModel- (- EDMComponent- , edmModel- ) where--import qualified Data.Vector.Storable as V--import ELynx.Data.MarkovProcess.MixtureModel-import ELynx.Data.MarkovProcess.RateMatrix-import ELynx.Data.MarkovProcess.SubstitutionModel---- | An empirical mixture model component has a weight and a stationary--- distribution.-type EDMComponent = (Weight, V.Vector Double)---- | Create an EDM model from components and a substitution module construction--- manual from stationary distributions.-edmModel :: [EDMComponent] -> (StationaryDistribution -> SubstitutionModel)- -> MixtureModel-edmModel cs f =- MixtureModel n [ Component w (f d)- | (w, d) <- cs ]- where n = "EDM" ++ show (length cs)
− src/ELynx/Data/MarkovProcess/GammaRateHeterogeneity.hs
@@ -1,109 +0,0 @@-{- |-Module : ELynx.Data.MarkovProcess.GammaRateHeterogeneity-Description : Discrete gamma rate heterogeneity-Copyright : (c) Dominik Schrempf 2019-License : GPL-3--Maintainer : dominik.schrempf@gmail.com-Stability : unstable-Portability : portable--Creation date: Thu Feb 28 14:09:11 2019.--At the moment, a mixture model is used to emulate gamma rate heterogeneity. This-does not come with huge run time increases when simulating data. For inference-however, it would make a lot of sense to reuse the Eigendecomposition for all-rate heterogeneity components though.---}--module ELynx.Data.MarkovProcess.GammaRateHeterogeneity- ( summarizeGammaRateHeterogeneity- , expand- ) where--import qualified Data.ByteString.Lazy.Char8 as L-import Numeric.Integration.TanhSinh-import Statistics.Distribution-import Statistics.Distribution.Gamma--import qualified ELynx.Data.MarkovProcess.MixtureModel as M-import qualified ELynx.Data.MarkovProcess.PhyloModel as P-import qualified ELynx.Data.MarkovProcess.SubstitutionModel as S---- | Short summary of gamma rate heterogeneity parameters.-summarizeGammaRateHeterogeneity :: Int -> Double -> [L.ByteString]-summarizeGammaRateHeterogeneity n alpha = map L.pack- [ "Discrete gamma rate heterogeneity."- , "Number of categories: " ++ show n- , "Shape parameter of gamma distribution: "++ show alpha- , "Rates: " ++ show (getMeans n alpha) ]---- | For a given number of rate categories, a gamma shape parameter alpha and a--- substitution model, compute the scaled substitution models corresponding to--- the gamma rates.-expand :: Int -> Double -> P.PhyloModel -> P.PhyloModel-expand n alpha (P.SubstitutionModel sm)- = P.MixtureModel $ expandSubstitutionModel n alpha sm-expand n alpha (P.MixtureModel mm)- = P.MixtureModel $ expandMixtureModel n alpha mm--getName :: Int -> Double -> String-getName n alpha = " with discrete gamma rate heterogeneity; "- ++ show n ++ " categories; "- ++ "shape parameter " ++ show alpha--splitSubstitutionModel :: Int -> Double -> S.SubstitutionModel -> [S.SubstitutionModel]-splitSubstitutionModel n alpha sm = renamedSMs- where- means = getMeans n alpha- scaledSMs = map (`S.scale` sm) means- names = map (("; gamma rate category " ++) . show) [1 :: Int ..]- renamedSMs = zipWith S.appendName names scaledSMs--expandSubstitutionModel :: Int -> Double -> S.SubstitutionModel -> M.MixtureModel-expandSubstitutionModel n alpha sm = M.fromSubstitutionModels name (repeat 1.0) sms- where- name = S.name sm <> getName n alpha- sms = splitSubstitutionModel n alpha sm--expandMixtureModel :: Int -> Double -> M.MixtureModel -> M.MixtureModel-expandMixtureModel n alpha mm = M.concatenate name renamedMMs- where- name = M.name mm <> getName n alpha- means = getMeans n alpha- scaledMMs = map (`M.scale` mm) means- names = map (("; gamma rate category " ++) . show) [1 :: Int ..]- renamedMMs = zipWith M.appendName names scaledMMs---- For a given number of rate categories 'n' and a shape parameter 'alpha' (the--- rate or scale is set such that the mean is 1.0), return a list of rates that--- represent the respective categories. Use the mean rate for each category.-getMeans :: Int -> Double -> [Double]-getMeans n alpha = means ++ lastMean- where gamma = gammaDistr alpha (1.0/alpha)- quantiles = [ quantile gamma (fromIntegral i / fromIntegral n) | i <- [0..n] ]- -- Calculate the mean rate. Multiplication with the number of rate- -- categories 'n' is necessary because in each n-quantile the- -- probability mass is 1/n.- meanFunc x = fromIntegral n * x * density gamma x- -- Only calculate the first (n-1) categories with normal integration.- means = [ integralAToB meanFunc (quantiles !! i) (quantiles !! (i+1)) | i <- [0..n-2] ]- -- The last category has to be calculated with an improper integration.- lastMean = [integralAToInf meanFunc (quantiles !! (n-1))]---- The error of integration.-eps :: Double-eps = 1e-6---- The integration method to use-method :: (Double -> Double ) -> Double -> Double -> [Result]-method = parSimpson---- Helper function for a normal integral from 'a' to 'b'.-integralAToB :: (Double -> Double) -> Double -> Double -> Double-integralAToB f a b = result . absolute eps $ method f a b---- Helper function for an improper integral from 'a' to infinity.-integralAToInf :: (Double -> Double) -> Double -> Double-integralAToInf f a = (result . absolute eps $ nonNegative method f) - integralAToB f 0 a
− src/ELynx/Data/MarkovProcess/MixtureModel.hs
@@ -1,145 +0,0 @@-{- |-Module : ELynx.Data.MarkovProcess.MixtureModel-Description : Mixture models are a set of substitution models with weights-Copyright : (c) Dominik Schrempf 2019-License : GPL-3--Maintainer : dominik.schrempf@gmail.com-Stability : unstable-Portability : portable--Creation date: Tue Jan 29 19:17:40 2019.--To be imported qualified.---}--module ELynx.Data.MarkovProcess.MixtureModel- ( -- * Types- Weight- , Component (Component)- , MixtureModel (MixtureModel)- -- * Getters- , name- , getAlphabet- , getWeights- , getSubstitutionModels- -- * Building mixture models- , fromSubstitutionModels- -- * Transformations- , concatenate- , scale- , normalize- , appendName- -- * Tests- , isValid- -- * Output- , summarizeComponent- , summarize- ) where--import qualified Data.ByteString.Builder as L-import qualified Data.ByteString.Lazy.Char8 as L--import ELynx.Data.Alphabet.Alphabet-import qualified ELynx.Data.MarkovProcess.SubstitutionModel as S-import ELynx.Tools.Equality (allEqual)---- | Mixture model component weight.-type Weight = Double---- | A mixture model component has a weight and a substitution model.-data Component = Component- { weight :: Weight- , substModel :: S.SubstitutionModel- }- deriving (Show, Read)---- | A mixture model with its components.-data MixtureModel = MixtureModel- { name :: S.Name -- ^ Name- , components :: [Component]- }- deriving (Show, Read)---- | Get alphabet used with mixture model. Throws error if components use--- different 'Alphabet's.-getAlphabet :: MixtureModel -> Alphabet-getAlphabet mm = if isValid mm- then S.alphabet . substModel $ head (components mm)- else error "Mixture model is invalid."---- | Get weights.-getWeights :: MixtureModel -> [Weight]-getWeights = map weight . components---- | Get substitution models.-getSubstitutionModels :: MixtureModel -> [S.SubstitutionModel]-getSubstitutionModels = map substModel . components----- | Create a mixture model from a list of substitution models.-fromSubstitutionModels :: S.Name -> [Weight] -> [S.SubstitutionModel] -> MixtureModel-fromSubstitutionModels n ws sms = MixtureModel n comps- where comps = zipWith Component ws sms---- | Concatenate mixture models.-concatenate :: S.Name -> [MixtureModel] -> MixtureModel-concatenate n mms = MixtureModel n $ concatMap components mms--scaleComponent :: Double -> Component -> Component-scaleComponent s c = c {substModel = s'}- where s' = S.scale s $ substModel c---- | Scale all substitution models of the mixture model.-scale :: Double -> MixtureModel -> MixtureModel-scale s m = m {components = cs'}- where cs = components m- cs' = map (scaleComponent s) cs---- | Globally normalize a mixture model so that on average one event happens per--- unit time.-normalize :: MixtureModel -> MixtureModel-normalize mm = scale (1/c) mm- where c = sum $ zipWith (*) weights scales- weights = getWeights mm- scales = map S.totalRate $ getSubstitutionModels mm--appendNameComponent :: S.Name -> Component -> Component-appendNameComponent n c = c {substModel = s'}- where s' = S.appendName n $ substModel c---- | Append byte string to all substitution models of mixture model.-appendName :: S.Name -> MixtureModel -> MixtureModel-appendName n m = m {components = cs'}- where cs = components m- cs' = map (appendNameComponent n) cs---- | Checks if a mixture model is valid.------ XXX: For the future, a proper way of creating mixture models might be of--- interest. For example, not exporting the constructor nor the record fields--- and providing an algebraic way of creating mixture models (empty and--- addComponent which performs necessary checks).-isValid :: MixtureModel -> Bool-isValid m = not (null $ components m) && allEqual alphabets- where alphabets = map (S.alphabet . substModel) $ components m---- | Summarize a mixture model component; lines to be printed to screen or log.-summarizeComponent :: Component -> [L.ByteString]-summarizeComponent c =- L.pack "Weight: " <> (L.toLazyByteString . L.doubleDec $ weight c)- : S.summarize (substModel c)---- | Summarize a mixture model; lines to be printed to screen or log.-summarize :: MixtureModel -> [L.ByteString]-summarize m =- [ L.pack $ "Mixture model: " ++ name m ++ "."- , L.pack $ "Number of components: " ++ show n ++ "." ]- ++ detail- where- n = length $ components m- detail = if n <= 100- then concat [ L.pack ("Component " ++ show i ++ ":") : summarizeComponent c- | (i, c) <- zip [1 :: Int ..] (components m) ]- else []
− src/ELynx/Data/MarkovProcess/Nucleotide.hs
@@ -1,81 +0,0 @@-{- |-Module : ELynx.Data.MarkovProcess.Nucleotide-Description : Substitution models using nucleotides-Copyright : (c) Dominik Schrempf 2019-License : GPL-3--Maintainer : dominik.schrempf@gmail.com-Stability : unstable-Portability : portable--Creation date: Thu Jan 24 08:33:26 2019.--XXX: Maybe rename to something like /DNA substitution models/. Nucleotide ~-Alphabet; DNA ~ Character.--The order of nucleotides is A, C, G, T; see 'ELynx.Data.Character.Nucleotide'.---}--module ELynx.Data.MarkovProcess.Nucleotide- ( jc- , hky- ) where--import Numeric.LinearAlgebra hiding (normalize)--import ELynx.Data.Alphabet.Alphabet-import ELynx.Data.MarkovProcess.RateMatrix-import ELynx.Data.MarkovProcess.SubstitutionModel-import ELynx.Tools.Vector---- XXX: Another idea of structuring the code. This would probably be cleaner in--- the long run.---- data PhyloModel = MixtureModel | SubstitutionModel------- I think it could simply be:--- data PhyloModel = [(Weight, SubstitutionModel)]------- data MixtureModel = [(Weight, SubstitutionModel)]---- data SubstitutionModel = SMDNA DNASubstitutionModel | SMAA AASubstitutionModel---- data DNASubstitutionModel = JC | HKY Double StationaryDistribution---- data AASubstitutionModel = LG | ...--n :: Int--- n = length (alphabet :: [Nucleotide])--- Hard code this here. Reduces model dependencies, and number of nucleotides--- will not change.-n = 4---- | JC model matrix.-jcExch :: ExchangeabilityMatrix-jcExch =- (n><n)- [ 0.0, 1.0, 1.0, 1.0- , 1.0, 0.0, 1.0, 1.0- , 1.0, 1.0, 0.0, 1.0- , 1.0, 1.0, 1.0, 0.0 ]---- | JC substitution model.-jc :: SubstitutionModel-jc = substitutionModel DNA "JC" [] f jcExch- where f = uniformVec n--hkyExch :: Double -> ExchangeabilityMatrix-hkyExch k =- (n><n)- [ 0.0, 1.0, k, 1.0- , 1.0, 0.0, 1.0, k- , k, 1.0, 0.0, 1.0- , 1.0, k, 1.0, 0.0 ]---- | HKY substitution model.-hky :: Double -> StationaryDistribution -> SubstitutionModel-hky k f = substitutionModel DNA "HKY" [k] f e- where e = hkyExch k
− src/ELynx/Data/MarkovProcess/PhyloModel.hs
@@ -1,46 +0,0 @@-{- |-Module : ELynx.Data.MarkovProcess.PhyloModel-Description : Phylogenetic model-Copyright : (c) Dominik Schrempf 2019-License : GPL-3--Maintainer : dominik.schrempf@gmail.com-Stability : unstable-Portability : portable--Creation date: Fri Feb 1 12:43:06 2019.--A phylogenetic model is a complete description of the evolutionary process. At-the moment, it is either a mixture model or a plain substitution model, but more-complicated models may be added in the future.--To be imported qualified.---}--module ELynx.Data.MarkovProcess.PhyloModel- ( PhyloModel (..)- , getAlphabet- , summarize- ) where--import qualified Data.ByteString.Lazy.Char8 as L--import ELynx.Data.Alphabet.Alphabet-import qualified ELynx.Data.MarkovProcess.MixtureModel as M-import qualified ELynx.Data.MarkovProcess.SubstitutionModel as S---- | A phylogenetic model is a mixture model or a substitution model. More--- complicated models may be added.-data PhyloModel = MixtureModel M.MixtureModel | SubstitutionModel S.SubstitutionModel- deriving (Show, Read)---- | Extract code from phylogenetic model.-getAlphabet :: PhyloModel -> Alphabet-getAlphabet (MixtureModel mm) = M.getAlphabet mm-getAlphabet (SubstitutionModel sm) = S.alphabet sm---- | Summarize a phylogenetic model; lines to be printed to screen or log.-summarize :: PhyloModel -> [L.ByteString]-summarize (MixtureModel mm) = M.summarize mm-summarize (SubstitutionModel sm) = S.summarize sm
− src/ELynx/Data/MarkovProcess/RateMatrix.hs
@@ -1,95 +0,0 @@-{- |-Description : Rate matrix helper functions-Copyright : (c) Dominik Schrempf 2017-License : GPLv3--Maintainer : dominik.schrempf@gmail.com-Stability : unstable-Portability : non-portable (not tested)--Some helper functions that come handy when working with rate matrices of-continuous-time discrete-state Markov processes.--* Changelog--To be imported qualified.---}--module ELynx.Data.MarkovProcess.RateMatrix- ( RateMatrix- , ExchangeabilityMatrix- , StationaryDistribution- , totalRate- , normalize- , normalizeWith- , setDiagonal- , toExchangeabilityMatrix- , fromExchangeabilityMatrix- , getStationaryDistribution- ) where--import Numeric.LinearAlgebra hiding (normalize)-import Prelude hiding ((<>))--import ELynx.Tools.Equality-import ELynx.Tools.LinearAlgebra-import ELynx.Tools.Vector---- | A rate matrix is just a real matrix.-type RateMatrix = Matrix R---- | A matrix of exchangeabilities, we have q = e * pi, where q is a rate--- matrix, e is the exchangeability matrix and pi is the diagonal matrix--- containing the stationary frequency distribution.-type ExchangeabilityMatrix = Matrix R---- | Stationary distribution of a rate matrix.-type StationaryDistribution = Vector R---- | Get average number of substitutions per unit time.-totalRate :: StationaryDistribution -> RateMatrix -> Double-totalRate d m = norm_1 $ d <# matrixSetDiagToZero m---- | Normalizes a Markov process generator such that one event happens per unit--- time. Calculates stationary distribution from rate matrix.-normalize :: RateMatrix -> RateMatrix-normalize m = normalizeWith (getStationaryDistribution m) m---- | Normalizes a Markov process generator such that one event happens per unit--- time. Stationary distribution has to be given.-normalizeWith :: StationaryDistribution -> RateMatrix -> RateMatrix-normalizeWith d m = scale (1.0 / totalRate d m) m---- | Set the diagonal entries of a matrix such that the rows sum to 0.-setDiagonal :: RateMatrix -> RateMatrix-setDiagonal m = diagZeroes - diag (fromList rowSums)- where diagZeroes = matrixSetDiagToZero m- rowSums = map norm_1 $ toRows diagZeroes---- | Extract the exchangeability matrix from a rate matrix.-toExchangeabilityMatrix :: RateMatrix -> StationaryDistribution -> ExchangeabilityMatrix-toExchangeabilityMatrix m f = m <> diag oneOverF- where oneOverF = cmap (1.0/) f---- | Convert exchangeability matrix to rate matrix.-fromExchangeabilityMatrix :: ExchangeabilityMatrix -> StationaryDistribution -> RateMatrix-fromExchangeabilityMatrix em d = setDiagonal $ em <> diag d---- | Get stationary distribution from 'RateMatrix'. Involves eigendecomposition.--- If the given matrix does not satisfy the required properties of transition--- rate matrices and no eigenvector with an eigenvalue nearly equal to 0 is--- found, an error is thrown. Is there an easier way to calculate the stationary--- distribution or a better way to handle errors (of course I could use the--- Maybe monad, but then the error report is just delayed to the calling--- function)?-getStationaryDistribution :: RateMatrix -> StationaryDistribution-getStationaryDistribution m =- if magnitude (eVals ! i) `nearlyEq` 0- then normalizeSumVec 1.0 distReal- else error "Could not retrieve stationary distribution."- where- (eVals, eVecs) = eig (tr m)- i = minIndex eVals- distComplex = toColumns eVecs !! i- distReal = cmap realPart distComplex
− src/ELynx/Data/MarkovProcess/SubstitutionModel.hs
@@ -1,119 +0,0 @@-{- |-Module : ELynx.Data.MarkovProcess.SubstitutionModel-Description : Data type describing substitution model-Copyright : (c) Dominik Schrempf 2019-License : GPL-3--Maintainer : dominik.schrempf@gmail.com-Stability : unstable-Portability : portable--Creation date: Tue Jan 29 19:10:46 2019.--To be imported qualified.---}--module ELynx.Data.MarkovProcess.SubstitutionModel- ( -- * Types- Name- , Params- , SubstitutionModel- -- * Lenses and other accessors- , alphabet- , name- , stationaryDistribution- , exchangeabilityMatrix- , rateMatrix- , totalRate- -- * Building substitution models- , substitutionModel- , unnormalized- -- * Transformations- , scale- , normalize- , appendName- -- * Output- , summarize- ) where--import qualified Data.ByteString.Lazy.Char8 as L-import qualified Numeric.LinearAlgebra as LinAlg--import ELynx.Data.Alphabet.Alphabet-import qualified ELynx.Data.MarkovProcess.RateMatrix as R-import ELynx.Tools.Definitions-import ELynx.Tools.LinearAlgebra-import ELynx.Tools.Numeric---- | Name of substitution model; abstracted and subject to change.-type Name = String---- | Parameters of substitution model. May be the empty list.-type Params = [Double]---- | Complete definition of a substitution model. Create instances with--- 'substitutionModel'. A substitution model has an alphabet, a name, and a list--- of parameters (e.g., the kappa value for the HKY model). Further, the--- transition rate matrix is defined by a stationary distribution and a set of--- exchangeabilities.-data SubstitutionModel = SubstitutionModel- { alphabet :: Alphabet -- ^ Alphabet- , name :: Name -- ^ Name- , params :: Params -- ^ List of parameters- , stationaryDistribution :: R.StationaryDistribution -- ^ Stationary distribution- , exchangeabilityMatrix :: R.ExchangeabilityMatrix -- ^ Exchangeability matrix- }- deriving (Show, Read)---- | Calculate rate matrix from substitution model.-rateMatrix :: SubstitutionModel -> R.RateMatrix-rateMatrix sm = R.fromExchangeabilityMatrix (exchangeabilityMatrix sm) (stationaryDistribution sm)---- | Get scale of substitution model.-totalRate :: SubstitutionModel -> Double-totalRate sm = R.totalRate (stationaryDistribution sm) (rateMatrix sm)---- | Create normalized 'SubstitutionModel'. See 'normalize'.-substitutionModel :: Alphabet -> Name -> Params- -> R.StationaryDistribution -> R.ExchangeabilityMatrix- -> SubstitutionModel-substitutionModel c n ps d e = normalize $ SubstitutionModel c n ps d e---- | Create UNNORMALIZED 'SubstitutionModel'. See 'substitutionModel'.-unnormalized :: Alphabet -> Name -> Params- -> R.StationaryDistribution -> R.ExchangeabilityMatrix- -> SubstitutionModel-unnormalized = SubstitutionModel---- | Scale the rate of a substitution model by given factor.-scale :: Double -> SubstitutionModel -> SubstitutionModel-scale r sm = sm {exchangeabilityMatrix = em'}- where em' = LinAlg.scale r $ exchangeabilityMatrix sm---- | Normalize a substitution model, so that, on average, one substitution--- happens per unit time.-normalize :: SubstitutionModel -> SubstitutionModel-normalize sm = scale (1.0/r) sm- where r = totalRate sm---- | Abbend to name.-appendName :: Name -> SubstitutionModel -> SubstitutionModel-appendName n sm = sm {name = n'}- where n' = name sm <> n---- | Summarize a substitution model; lines to be printed to screen or log.-summarize :: SubstitutionModel -> [L.ByteString]-summarize sm = map L.pack $- (show (alphabet sm) ++ " substitution model: " ++ name sm ++ ".") :- [ "Parameters: " ++ show (params sm) ++ "." | not (null (params sm))] ++- case alphabet sm of- DNA -> [ "Stationary distribution: " ++ dispv precision (stationaryDistribution sm) ++ "."- , "Exchangeability matrix:\n" ++ dispmi 2 precision (exchangeabilityMatrix sm) ++ "."- , "Scale: " ++ show (roundN precision $ totalRate sm) ++ "."- ]- Protein -> [ "Stationary distribution: " ++ dispv precision (stationaryDistribution sm) ++ "."- , "Scale: " ++ show (roundN precision $ totalRate sm) ++ "."- ]- _ -> error "Extended character sets are not supported with substitution models."-
src/ELynx/Data/Sequence/Alignment.hs view
@@ -2,7 +2,7 @@ Module : ELynx.Data.Sequence.Alignment Description : Multi sequence alignment related types and functions Copyright : (c) Dominik Schrempf 2018-License : GPL-3+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -17,7 +17,7 @@ module ELynx.Data.Sequence.Alignment- ( Alignment (..)+ ( Alignment(..) , length , nSequences -- | * Input, output@@ -43,34 +43,36 @@ -- | * Sub sample , subSample , randomSubSample- ) where+ )+where -import Control.Monad hiding (join)+import Control.Monad hiding ( join ) import Control.Monad.Primitive-import qualified Data.ByteString.Lazy.Char8 as L-import Data.List hiding (concat,- length)-import qualified Data.Matrix.Unboxed as M-import qualified Data.Vector.Unboxed as V-import Prelude hiding (concat,- length)+import qualified Data.ByteString.Lazy.Char8 as L+import Data.List hiding ( concat+ , length+ )+import qualified Data.Matrix.Unboxed as M+import qualified Data.Vector.Unboxed as V+import Prelude hiding ( concat+ , length+ ) import System.Random.MWC -import qualified ELynx.Data.Alphabet.Alphabet as A+import qualified ELynx.Data.Alphabet.Alphabet as A import ELynx.Data.Alphabet.Character-import qualified ELynx.Data.Alphabet.DistributionDiversity as D+import qualified ELynx.Data.Alphabet.DistributionDiversity+ as D import ELynx.Data.Sequence.Defaults-import qualified ELynx.Data.Sequence.Sequence as S-import ELynx.Tools.Concurrent-import ELynx.Tools.Definitions-import ELynx.Tools.Equality-import ELynx.Tools.Matrix+import qualified ELynx.Data.Sequence.Sequence as S+import ELynx.Tools -- | A collection of sequences. data Alignment = Alignment- { names :: [S.Name]- , alphabet :: A.Alphabet- , matrix :: M.Matrix Character+ { names :: [S.Name]+ , descriptions :: [S.Description]+ , alphabet :: A.Alphabet+ , matrix :: M.Matrix Character } deriving (Show, Eq) @@ -85,23 +87,24 @@ -- | Create 'Alignment' from a list of 'S.Sequence's. fromSequences :: [S.Sequence] -> Either String Alignment fromSequences ss- | S.equalLength ss && allEqual (map S.alphabet ss) =- Right $ Alignment ns a d- | S.equalLength ss =- Left "Sequences do not have equal codes."- | otherwise =- Left "Sequences do not have equal lengths."- where- ns = map S.name ss- a = S.alphabet $ head ss- bss = map S.characters ss- d = M.fromRows bss+ | S.equalLength ss && allEqual (map S.alphabet ss) = Right+ $ Alignment ns ds a d+ | S.equalLength ss = Left "Sequences do not have equal codes."+ | otherwise = Left "Sequences do not have equal lengths."+ where+ ns = map S.name ss+ ds = map S.description ss+ a = S.alphabet $ head ss+ bss = map S.characters ss+ d = M.fromRows bss -- | Conversion to list of 'S.Sequence's. toSequences :: Alignment -> [S.Sequence]-toSequences (Alignment ns a d) = zipWith (\n r -> S.Sequence n a r) ns rows- where- rows = M.toRows d+toSequences (Alignment ns ds a da) = zipWith3 (\n d r -> S.Sequence n d a r)+ ns+ ds+ rows+ where rows = M.toRows da -- -- | Show a 'Alignment', untrimmed. -- summarizeSeq :: Alignment -> Int -> L.ByteString@@ -111,19 +114,30 @@ -- L.pack $ V.toList $ V.map toChar $ M.takeRow (matrix m) i ] header :: Alignment -> L.ByteString-header a = L.unlines $- [ L.pack "Multi sequence alignment."- , L.pack $ "Code: " ++ A.description (alphabet a) ++ "."- , L.pack $ "Length: " ++ show (length a) ++ "." ]- ++ reportLengthSummary ++ reportNumberSummary- where reportLengthSummary =- [ L.pack $ "For each sequence, the "- ++ show summaryLength ++ " first bases are shown."- | length a > summaryLength ]- reportNumberSummary =- [ L.pack $ show summaryNSequences ++ " out of " ++- show (nSequences a) ++ " sequences are shown."- | nSequences a > summaryNSequences ]+header a =+ L.unlines+ $ [ L.pack "Multi sequence alignment."+ , L.pack $ "Code: " ++ A.alphabetDescription (alphabet a) ++ "."+ , L.pack $ "Length: " ++ show (length a) ++ "."+ ]+ ++ reportLengthSummary+ ++ reportNumberSummary+ where+ reportLengthSummary =+ [ L.pack+ $ "For each sequence, the "+ ++ show summaryLength+ ++ " first bases are shown."+ | length a > summaryLength+ ]+ reportNumberSummary =+ [ L.pack+ $ show summaryNSequences+ ++ " out of "+ ++ show (nSequences a)+ ++ " sequences are shown."+ | nSequences a > summaryNSequences+ ] -- | Similar to 'S.summarizeSequenceList' but with different Header. summarize :: Alignment -> L.ByteString@@ -131,36 +145,40 @@ -- | Join two 'Alignment's vertically. That is, add more sequences -- to an alignment. See also 'concat'.-join :: Alignment- -> Alignment- -> Alignment+join :: Alignment -> Alignment -> Alignment -- top bottom. join t b- | length t == length b &&- al == alphabet b- = Alignment ns al (tD === bD)- | otherwise = error "join: Multi sequence alignments do not have equal length."- where- ns = names t ++ names b- tD = matrix t- bD = matrix b- al = alphabet t+ | length t /= length b = error+ "join: Multi sequence alignments do not have equal lengths."+ | alphabet t /= alphabet b = error+ "join: Multi sequence alignments do not have equal alphabets."+ | otherwise = Alignment ns ds al (tD === bD)+ where+ ns = names t ++ names b+ ds = descriptions t ++ descriptions b+ tD = matrix t+ bD = matrix b+ al = alphabet t -- | Concatenate two 'Alignment's horizontally. That is, add more -- sites to an alignment. See also 'join'.-concat :: Alignment- -> Alignment- -> Alignment+concat :: Alignment -> Alignment -> Alignment -- left right. concat l r- | nSequences l == nSequences r &&- al == alphabet r- = Alignment (names l) al (lD ||| rD)- | otherwise = error "concat: Multi sequence alignments do not have equal length."- where- lD = matrix l- rD = matrix r- al = alphabet l+ | nSequences l /= nSequences r+ = error+ "concat: Multi sequence alignments do not have an equal number of sequences."+ | alphabet l /= alphabet r+ = error "concat: Multi sequence alignments do not have an equal alphabets."+ | names l /= names r+ = error "concat: Multi sequence alignments do not have an equal names."+ | descriptions l /= descriptions r+ = error "concat: Multi sequence alignments do not have an equal descriptions."+ | otherwise+ = Alignment (names l) (descriptions l) (alphabet l) (lD ||| rD)+ where+ lD = matrix l+ rD = matrix r -- | Concatenate a list of 'Alignment's horizontally. See -- 'concat'.@@ -171,7 +189,7 @@ -- Only keep columns from alignment that satisfy given predicate. filterColsWith :: (V.Vector Character -> Bool) -> Alignment -> Alignment-filterColsWith p a = a {matrix = m'}+filterColsWith p a = a { matrix = m' } where m' = M.fromColumns . filter p . M.toColumns $ matrix a -- | Only keep columns with standard characters. Alignment columns with IUPAC@@ -182,9 +200,11 @@ -- | Filter columns with proportion of standard character larger than given number. filterColsStd :: Double -> Alignment -> Alignment filterColsStd prop a = filterColsWith- (\col -> prop * n <= fromIntegral (V.length (V.filter (A.isStd al) col))) a- where al = alphabet a- n = fromIntegral $ nSequences a+ (\col -> prop * n <= fromIntegral (V.length (V.filter (A.isStd al) col)))+ a+ where+ al = alphabet a+ n = fromIntegral $ nSequences a -- | Only keep columns without gaps or unknown characters. filterColsNoGaps :: Alignment -> Alignment@@ -202,8 +222,9 @@ -- | Calculate the distribution of characters. distribution :: FrequencyData -> [Double]-distribution fd = map (/ fromIntegral nSites) $ V.toList $- foldl1 (V.zipWith (+)) (M.toColumns fd)+distribution fd = map (/ fromIntegral nSites) $ V.toList $ foldl1+ (V.zipWith (+))+ (M.toColumns fd) where nSites = M.cols fd -- | Diversity analysis. See 'kEffEntropy'.@@ -232,13 +253,12 @@ -- | Sample the given sites from a multi sequence alignment. subSample :: [Int] -> Alignment -> Alignment-subSample is a = a {matrix = m'}- where m' = subSampleMatrix is $ matrix a+subSample is a = a { matrix = m' } where m' = subSampleMatrix is $ matrix a -- | Randomly sample a given number of sites of the multi sequence alignment.-randomSubSample :: PrimMonad m- => Int -> Alignment -> Gen (PrimState m) -> m Alignment+randomSubSample+ :: PrimMonad m => Int -> Alignment -> Gen (PrimState m) -> m Alignment randomSubSample n a g = do let l = length a- is <- replicateM n $ uniformR (0, l-1) g+ is <- replicateM n $ uniformR (0, l - 1) g return $ subSample is a
src/ELynx/Data/Sequence/Defaults.hs view
@@ -2,7 +2,7 @@ Module : ELynx.Defaults Description : Various default values Copyright : (c) Dominik Schrempf 2018-License : GPL-3+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -18,7 +18,8 @@ , summaryLength , summaryNSequences , fieldWidth- ) where+ )+where -- | Space reserved for sequence names when printing them. nameWidth :: Int
src/ELynx/Data/Sequence/Sequence.hs view
@@ -4,7 +4,7 @@ Module : ELynx.Data.Sequence Description : Hereditary sequences Copyright : (c) Dominik Schrempf 2018-License : GPL-3+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -19,8 +19,9 @@ module ELynx.Data.Sequence.Sequence ( -- * Types Name+ , Description , Characters- , Sequence (..)+ , Sequence(..) -- * Input , fromByteString -- * Output@@ -41,26 +42,32 @@ , filterShorterThan , filterLongerThan , filterStandard- ) where+ )+where import Control.Parallel.Strategies import qualified Data.ByteString.Lazy.Char8 as L-import Data.List (maximumBy)-import Data.Ord (comparing)+import Data.List ( maximumBy )+import Data.Ord ( comparing ) import qualified Data.Vector.Unboxed as V-import Prelude hiding (concat, length)-import qualified Prelude as Pr (length)+import Prelude hiding ( concat+ , length+ )+import qualified Prelude as Pr+ ( length ) import qualified Text.Printf as P import qualified ELynx.Data.Alphabet.Alphabet as A import ELynx.Data.Alphabet.Character import ELynx.Data.Sequence.Defaults-import ELynx.Tools.ByteString-import ELynx.Tools.Equality+import ELynx.Tools -- | For now, 'Name's are just 'L.ByteString's. type Name = L.ByteString +-- | The description of a sequence.+type Description = L.ByteString+ -- | The vector of characters of a sequence. type Characters = V.Vector Character @@ -72,24 +79,30 @@ toByteString :: Characters -> L.ByteString toByteString = L.pack . map toChar . V.toList --- | Sequences have a name, a code and hopefully a lot of data.-data Sequence = Sequence { name :: Name- , alphabet :: A.Alphabet- , characters :: Characters }+-- | Sequences have a name, a possibly empty description, a code and hopefully a+-- lot of data.+data Sequence = Sequence { name :: Name+ , description :: Description+ , alphabet :: A.Alphabet+ , characters :: Characters } deriving (Show, Eq) getInfo :: Sequence -> L.ByteString-getInfo s = L.unwords [ alignLeft nameWidth (name s)- , alignRight fieldWidth (L.pack $ show $ alphabet s)- , alignRight fieldWidth (L.pack . show $ len)- , alignRight fieldWidth (L.pack $ P.printf "%.3f" pGaps) ]- where len = length s- nGaps = countGaps s- pGaps = fromIntegral nGaps / fromIntegral len :: Double+getInfo s = L.unwords+ [ alignLeft nameWidth (name s)+ , alignRight fieldWidth (L.pack $ show $ alphabet s)+ , alignRight fieldWidth (L.pack . show $ len)+ , alignRight fieldWidth (L.pack $ P.printf "%2.2f" pGaps)+ ]+ where+ len = length s+ nGaps = countGaps s+ pGaps = 100 * fromIntegral nGaps / fromIntegral len :: Double -- | Trim and show a 'Sequence'. summarize :: Sequence -> L.ByteString-summarize s = L.unwords [ getInfo s , summarizeByteString summaryLength $ toByteString (characters s) ]+summarize s = L.unwords+ [getInfo s, summarizeByteString summaryLength $ toByteString (characters s)] -- | Trim and show a list of 'Sequence's. summarizeSequences :: [Sequence] -> L.ByteString@@ -97,26 +110,36 @@ -- | Header printed before 'Sequence' list. tableHeader :: L.ByteString-tableHeader = L.unwords [ alignLeft nameWidth "Name"- , alignRight fieldWidth "Code"- , alignRight fieldWidth "Length"- , alignRight fieldWidth "Gaps [%]"- , "Sequence" ]+tableHeader = L.unwords+ [ alignLeft nameWidth "Name"+ , alignRight fieldWidth "Code"+ , alignRight fieldWidth "Length"+ , alignRight fieldWidth "Gaps [%]"+ , "Sequence"+ ] -- | A short description of the sequence. header :: [Sequence] -> L.ByteString-header ss = L.unlines $- reportIfSubsetIsShown ++- [ L.pack $ "For each sequence, the " ++ show summaryLength ++ " first bases are shown."- , L.pack $ "List contains " ++ show (Pr.length ss) ++ " sequences."- , ""- , tableHeader ]- where l = Pr.length ss- s = show summaryNSequences ++ " out of " ++- show (Pr.length ss) ++ " sequences are shown."- reportIfSubsetIsShown- | l > summaryNSequences = [L.pack s]- | otherwise = []+header ss =+ L.unlines+ $ reportIfSubsetIsShown+ ++ [ L.pack+ $ "For each sequence, the "+ ++ show summaryLength+ ++ " first bases are shown."+ , L.pack $ "List contains " ++ show (Pr.length ss) ++ " sequences."+ , ""+ , tableHeader+ ]+ where+ l = Pr.length ss+ s =+ show summaryNSequences+ ++ " out of "+ ++ show (Pr.length ss)+ ++ " sequences are shown."+ reportIfSubsetIsShown | l > summaryNSequences = [L.pack s]+ | otherwise = [] -- | Trim and show a list of 'Sequence's. body :: [Sequence] -> L.ByteString@@ -140,14 +163,34 @@ -- | Trim to given length. trim :: Int -> Sequence -> Sequence-trim n (Sequence nm a cs) = Sequence nm a (V.take (fromIntegral n) cs)+trim n (Sequence nm d a cs) = Sequence nm d a (V.take (fromIntegral n) cs) -- | Concatenate two sequences. 'Name's have to match. concat :: Sequence -> Sequence -> Sequence-concat (Sequence i c cs) (Sequence j k ks)- | i == j && c == k = Sequence i c (cs <> ks)- | otherwise = error $ "concatenate: Sequences do not have equal names: "- ++ L.unpack i ++ ", " ++ L.unpack j ++ "."+concat (Sequence i d c cs) (Sequence j f k ks)+ | i /= j+ = error+ $ "concatenate: Sequences do not have equal names: "+ ++ L.unpack i+ ++ ", "+ ++ L.unpack j+ ++ "."+ | d /= f+ = error+ $ "concatenate: Sequences do not have equal descriptions: "+ ++ L.unpack d+ ++ ", "+ ++ L.unpack f+ ++ "."+ | c /= k+ = error+ $ "concatenate: Sequences do not have equal alphabets: "+ ++ show c+ ++ ", "+ ++ show k+ ++ "."+ | otherwise+ = Sequence i d c (cs <> ks) -- | Concatenate a list of sequences, see 'concat'. concatSequences :: [[Sequence]] -> [Sequence]@@ -169,5 +212,4 @@ -- Are all characters IUPAC characters? anyStandard :: A.Alphabet -> Sequence -> Bool-anyStandard a s = V.any (A.isStd a) cs- where cs = characters s+anyStandard a s = V.any (A.isStd a) cs where cs = characters s
src/ELynx/Data/Sequence/Translate.hs view
@@ -1,8 +1,8 @@ {- | Module : ELynx.Data.Sequence.Translate Description : Translate sequences-Copyright : (c) Dominik Schrempf 2019-License : GPL-3+Copyright : (c) Dominik Schrempf 2020+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -23,22 +23,28 @@ import qualified ELynx.Data.Alphabet.Character as C import ELynx.Data.Character.Codon import ELynx.Data.Sequence.Sequence-import ELynx.Tools.Vector+import ELynx.Tools -- | Translate a sequence from 'DNA' or 'DNAX' to 'ProteinS'. translateSeq :: UniversalCode -> Int -> Sequence -> Sequence-translateSeq uc rf (Sequence n a cs) = case a of- DNA -> Sequence n ProteinS (cs' $ translate uc)- DNAX -> Sequence n ProteinS (cs' $ translateX uc)- DNAI -> Sequence n ProteinI (cs' $ translateI uc)- _ -> error "translate: can only translate DNA, DNAX, and DNAI."+translateSeq uc rf (Sequence n d a cs) = case a of+ DNA -> Sequence n d ProteinS (cs' $ translate uc)+ DNAX -> Sequence n d ProteinS (cs' $ translateX uc)+ DNAI -> Sequence n d ProteinI (cs' $ translateI uc)+ _ -> error "translate: can only translate DNA, DNAX, and DNAI." where cs' f = C.fromCVec $ translateVecWith f rf (C.toCVec cs) -- Translate from DNA to Protein with given reading frame (0, 1, 2).-translateVecWith :: (V.Unbox a, Ord a, V.Unbox b)- => (Codon a -> b) -> Int -> V.Vector a -> V.Vector b-translateVecWith f rf cs | rf > 2 = error "translateVecWith: reading frame is larger than 2."- | rf < 0 = error "translateVecWith: reading frame is negative."- | otherwise = aas- where codons = map unsafeFromVec $ chop 3 $ V.drop rf cs- aas = V.fromList $ map f codons+translateVecWith+ :: (V.Unbox a, Ord a, V.Unbox b)+ => (Codon a -> b)+ -> Int+ -> V.Vector a+ -> V.Vector b+translateVecWith f rf cs+ | rf > 2 = error "translateVecWith: reading frame is larger than 2."+ | rf < 0 = error "translateVecWith: reading frame is negative."+ | otherwise = aas+ where+ codons = map unsafeFromVec $ chopVec 3 $ V.drop rf cs+ aas = V.fromList $ map f codons
− src/ELynx/Export/Sequence/CountsFile.hs
@@ -1,109 +0,0 @@-{- |-Description : Write a counts file-Copyright : (c) Dominik Schrempf 2017-License : GPLv3--Maintainer : dominik.schrempf@gmail.com-Stability : unstable-Portability : non-portable (not tested)--TODO: Import.--* The Counts Format--The input of PoMo is allele frequency data. Especially, when-populations have many individuals it is preferable to count the-number of bases at each position. This decreases file size and speeds-up the parser.--Counts files contain:--- One headerline that specifies the file as counts file and states the- number of populations as well as the number of sites (separated by- white space).--- A second headerline with white space separated headers: CRHOM- (chromosome), POS (position) and sequence names.--- Many lines with counts of A, C, G and T bases and their respective- positions.--Comments:--- Lines starting with # before the first headerline are treated as- comments.--A toy example:--@- COUNTSFILE NPOP 5 NSITES N- CHROM POS Sheep BlackSheep RedSheep Wolf RedWolf- 1 1 0,0,1,0 0,0,1,0 0,0,1,0 0,0,5,0 0,0,0,1- 1 2 0,0,0,1 0,0,0,1 0,0,0,1 0,0,0,5 0,0,0,1- .- .- .- 9 8373 0,0,0,1 1,0,0,0 0,1,0,0 0,1,4,0 0,0,1,0- .- .- .- Y 9999 0,0,0,1 0,1,0,0 0,1,0,0 0,5,0,0 0,0,1,0-@---}--module ELynx.Export.Sequence.CountsFile- ( Chrom- , Pos- , DataOneSite- , PopulationNames- , toCountsFile- ) where--import qualified Data.ByteString.Lazy.Char8 as L-import Data.Maybe (fromMaybe)--import ELynx.Data.Character.BoundaryMutation-import ELynx.Tools.ByteString (alignLeft, alignRight)---- | The number of sites that will be printed.-type NSites = Int---- | The names of the populations.-type PopulationNames = [L.ByteString]---- Desired column width of the counts file.-colW :: Int-colW = 11---- | Compose the header using the number of sites and the population names.-header :: NSites -> PopulationNames -> L.ByteString-header nSites popNames = L.unlines [lineOne, lineTwo]- where nPop = length popNames- lineOne = L.pack $ "COUNTSFILE NPOP " ++ show nPop ++ " NSITES " ++ show nSites- lineTwo = L.unwords $- [ alignLeft colW $ L.pack "CHROM"- , alignRight colW $ L.pack "POS" ]- ++ map (alignLeft colW) popNames---- | The chromosome name.-type Chrom = L.ByteString---- | The position on the chromosome.-type Pos = Int---- | The set of boundary states for one site.-type DataOneSite = [State]---- | Get a data line in the counts file.-dataLine :: Maybe Chrom -> Maybe Pos -> DataOneSite -> L.ByteString-dataLine chrom mPos bstates = L.unwords $- [ alignLeft colW (fromMaybe (L.pack "NA") chrom)- , alignRight colW (L.pack (maybe "NaN" show mPos)) ]- ++ map (alignRight colW . showCounts) bstates---- | Convert data to a counts file.-toCountsFile :: PopulationNames -> [(Maybe Chrom, Maybe Pos, DataOneSite)] -> L.ByteString-toCountsFile ns d = L.unlines $ header l ns : zipWith3 dataLine cs ps ds- where l = length d- (cs, ps, ds) = unzip3 d
src/ELynx/Export/Sequence/Fasta.hs view
@@ -2,7 +2,7 @@ Module : ELynx.Export.Sequence.Fasta Description : Export Fasta sequences Copyright : (c) Dominik Schrempf 2018-License : GPL-3+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com@@ -19,18 +19,21 @@ module ELynx.Export.Sequence.Fasta ( sequenceToFasta , sequencesToFasta- ) where+ )+where -import qualified Data.ByteString.Lazy.Char8 as L+import qualified Data.ByteString.Lazy.Char8 as L import ELynx.Data.Sequence.Sequence -fastaHeader :: L.ByteString -> L.ByteString-fastaHeader i = L.singleton '>' <> i+fastaHeader :: L.ByteString -> L.ByteString -> L.ByteString+fastaHeader n d =+ L.singleton '>' <> n <> if L.null d then L.empty else L.pack " " <> d -- | Convert a 'Sequence' to Fasta format. sequenceToFasta :: Sequence -> L.ByteString-sequenceToFasta s = L.unlines [ fastaHeader $ name s, toByteString $ characters s ]+sequenceToFasta s =+ L.unlines [fastaHeader (name s) (description s), toByteString $ characters s] -- | Convert a list 'Sequence's to Fasta format. A newline is added between any -- two 'Sequence's.
− src/ELynx/Import/MarkovProcess/EDMModelPhylobayes.hs
@@ -1,70 +0,0 @@-{- |-Module : ELynx.Import.MarkovProcess.EDMModelPhylobayes-Description : Import stationary distributions from Phylobayes format-Copyright : (c) Dominik Schrempf 2019-License : GPL-3--Maintainer : dominik.schrempf@gmail.com-Stability : unstable-Portability : portable--Creation date: Tue Jan 29 12:12:55 2019.---}--module ELynx.Import.MarkovProcess.EDMModelPhylobayes- ( Parser- , EDMComponent- , phylobayes- ) where--import Control.Monad-import qualified Data.ByteString.Lazy.Char8 as L-import qualified Data.Vector.Storable as V-import Data.Void-import Text.Megaparsec-import Text.Megaparsec.Byte-import Text.Megaparsec.Byte.Lexer--import ELynx.Data.MarkovProcess.EDMModel-import ELynx.Tools.ByteString (c2w)---- | Shortcut.-type Parser = Parsec Void L.ByteString---- | Parse stationary distributions from Phylobayes format.-phylobayes :: Parser [EDMComponent]-phylobayes = do- n <- headerLine- k <- kComponentsLine- cs <- count k $ dataLine n- _ <- many newline *> eof- <?> "phylobayes"- return cs--horizontalSpace :: Parser ()-horizontalSpace = skipMany $ char (c2w ' ') <|> tab--headerLine :: Parser Int-headerLine = do- n <- decimal- _ <- horizontalSpace- -- FIXME: This should be more general, but then we also want to ensure that- -- the order of states is correct.- _ <- chunk (L.pack "A C D E F G H I K L M N P Q R S T V W Y") <|> chunk (L.pack "A C G T")- _ <- many newline- <?> "headerLine"- return n--kComponentsLine :: Parser Int-kComponentsLine = decimal <* newline <?> "kComponentsLine"--dataLine :: Int -> Parser EDMComponent-dataLine n = do- weight <- float- _ <- horizontalSpace- vals <- float `sepBy` horizontalSpace- when (length vals /= n) (error "Did not find correct number of entries.")- _ <- many newline- <?> "dataLine"- return (weight, V.fromList vals)
− src/ELynx/Import/MarkovProcess/SiteprofilesPhylobayes.hs
@@ -1,77 +0,0 @@-{- |-Module : ELynx.Import.MarkovProcess.SiteprofilesPhylobayes-Description : Import site profiles in Phylobayes format-Copyright : (c) Dominik Schrempf 2019-License : GPL-3--Maintainer : dominik.schrempf@gmail.com-Stability : unstable-Portability : portable--Creation date: Tue Jan 29 12:12:55 2019.--For now I just try to go with a huge empirical distribution mixture model. Let's-see if performance is good enough.---}--module ELynx.Import.MarkovProcess.SiteprofilesPhylobayes- ( Parser- , EDMComponent- , siteprofiles- ) where--import Control.Monad-import qualified Data.ByteString.Lazy.Char8 as L-import Data.List (nub)-import qualified Data.Vector.Storable as V-import Data.Void-import Text.Megaparsec-import Text.Megaparsec.Byte-import Text.Megaparsec.Byte.Lexer--import ELynx.Data.MarkovProcess.EDMModel-import ELynx.Tools.ByteString (c2w)---- | Shortcut.-type Parser = Parsec Void L.ByteString---- | Parse stationary distributions from Phylobayes format.-siteprofiles :: Parser [EDMComponent]-siteprofiles = do- _ <- headerLines- cs <- many dataLine- _ <- many newline *> eof- <?> "phylobayes siteprofiles"- let ls = map length cs- nLs = length $ nub ls- when (nLs /= 1) (error "The site profiles have a different number of entries.")- return cs--horizontalSpace :: Parser ()-horizontalSpace = skipMany $ char (c2w ' ') <|> tab--line :: Parser ()-line = do- _ <- many $ noneOf [c2w '\n']- pure ()---- For now, just ignore the header.-headerLines :: Parser ()-headerLines = do- _ <- line- _ <- many newline- <?> "headerLine"- pure ()--dataLine :: Parser EDMComponent-dataLine = do- -- Ignore site number.- _ <- decimal :: Parser Integer- _ <- horizontalSpace- -- Also ignore additional white space on line.- vals <- float `sepEndBy1` horizontalSpace- _ <- many newline- <?> "dataLine"- -- Set the weight to 1.0 for all sites.- return (1.0, V.fromList vals)
src/ELynx/Import/Sequence/Fasta.hs view
@@ -4,7 +4,7 @@ Module : ELynx.Import.Sequence.Fasta Description : Import Fasta sequences Copyright : (c) Dominik Schrempf 2018-License : GPL-3+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com@@ -23,20 +23,23 @@ ( Parser , fastaSequence , fasta- ) where+ )+where import Control.Monad import qualified Data.ByteString.Lazy.Char8 as L import qualified Data.Set as S import Data.Void-import Data.Word8+import Data.Word8 ( Word8+ , isAlphaNum+ ) import Text.Megaparsec import Text.Megaparsec.Byte import ELynx.Data.Alphabet.Alphabet as A import ELynx.Data.Alphabet.Character import ELynx.Data.Sequence.Sequence-import ELynx.Tools.ByteString (c2w)+import ELynx.Tools -- | Shortcut. type Parser = Parsec Void L.ByteString@@ -44,16 +47,20 @@ isSpecial :: Word8 -> Bool isSpecial w = w `elem` map c2w ['_', '|', '.', '-'] -isHeaderChar :: Word8 -> Bool-isHeaderChar w = isAlphaNum w || isSpecial w+isHeader :: Word8 -> Bool+isHeader w = isAlphaNum w || isSpecial w -sequenceHeader :: Parser L.ByteString+isHorizontalSpace :: Word8 -> Bool+isHorizontalSpace w = (w == c2w ' ') || (w == c2w '\t')++sequenceHeader :: Parser (L.ByteString, L.ByteString) sequenceHeader = do _ <- char (c2w '>')- h <- takeWhile1P (Just "Header character") isHeaderChar- -- XXX: Allow description.+ n <- takeWhile1P (Just "Name") isHeader+ _ <- takeWhileP (Just "Horizontal space") isHorizontalSpace+ d <- takeWhileP (Just "Description") isHeader _ <- eol- return h+ return (n, d) -- It is a little faster to directly pass the set of allowed characters. Then, -- this set only has to be calculcated once per sequence in 'fastaSequence'.@@ -61,7 +68,7 @@ sequenceLine s = do -- FIXME: Will fail for non-capital letters. !xs <- takeWhile1P (Just "Alphabet character") (`S.member` s)- _ <- void eol <|> eof+ _ <- void eol <|> eof return xs -- XXX: If sequences are parsed line by line, the lines have to be copied when@@ -69,11 +76,12 @@ -- | Parse a sequence of characters. fastaSequence :: Alphabet -> Parser Sequence-fastaSequence a = do hd <- sequenceHeader- let !alph = S.map toWord (A.all . alphabetSpec $ a)- lns <- some (sequenceLine alph)- _ <- many eol- return $ Sequence hd a (fromByteString $ L.concat lns)+fastaSequence a = do+ (n, d) <- sequenceHeader+ let !alph = S.map toWord (A.all . alphabetSpec $ a)+ lns <- some (sequenceLine alph)+ _ <- many eol+ return $ Sequence n d a (fromByteString $ L.concat lns) -- | Parse a Fasta file with given 'Alphabet'. fasta :: Alphabet -> Parser [Sequence]
− src/ELynx/Simulate/MarkovProcess.hs
@@ -1,76 +0,0 @@-{- |-Module : ELynx.Simulate.MarkovProcess-Description : Markov process helpers-Copyright : (c) Dominik Schrempf 2019-License : GPL-3--Maintainer : dominik.schrempf@gmail.com-Stability : unstable-Portability : portable--Creation date: Thu Jan 24 09:02:25 2019.---}--module ELynx.Simulate.MarkovProcess- ( ProbMatrix- , State- , probMatrix- , jump- ) where--import Control.Monad.Primitive-import Numeric.LinearAlgebra-import System.Random.MWC-import System.Random.MWC.Distributions--import ELynx.Data.MarkovProcess.RateMatrix---- | A probability matrix, P_ij(t) = Pr (X_t = j | X_0 = i).-type ProbMatrix = Matrix R---- | Make type signatures a little clearer.-type State = Int---- | The important matrix that gives the probabilities to move from one state to--- another in a specific time (branch length).-probMatrix :: RateMatrix -> Double -> ProbMatrix-probMatrix q t | t == 0 = if rows q == cols q- then ident (rows q)- else error "probMatrix: Matrix is not square."- | t < 0 = error "probMatrix: Time is negative."- | otherwise = expm $ scale t q---- | Move from a given state to a new one according to a transition probability--- matrix (for performance reasons this probability matrix needs to be given as--- a list of generators, see--- https://hackage.haskell.org/package/distribution-1.1.0.0/docs/Data-Distribution-Sample.html).--- This function is the bottleneck of the simulator and takes up most of the--- computation time. However, I was not able to find a faster implementation--- than the one from Data.Distribution.-jump :: (PrimMonad m) => State -> ProbMatrix -> Gen (PrimState m) -> m State-jump i p = categorical (p ! i)---- XXX: Maybe for later, use condensed tables.------ Write storable instance, compilation is really slow otherwise. instance--- Storable (Int, R) where sizeOf (x, y) = sizeOf x + sizeOf y------ Do not generate table for each jump.------ jump :: (PrimMonad m) => State -> ProbMatrix -> Gen (PrimState m) -> m State--- jump i p = genFromTable table--- where--- ws = toList $ p ! i--- vsAndWs = fromList [ (v, w) | (v, w) <- zip [(0 :: Int) ..] ws--- , w > 0 ]--- table = tableFromProbabilities vsAndWs---- -- | Perform N jumps from a given state and according to a transition--- -- probability matrix transformed to a list of generators. This implementation--- -- uses 'foldM' and I am not sure how to access or store the actual chain. This--- -- could be done by an equivalent of 'scanl' for general monads, which I was--- -- unable to find. This function is neat, but will most likely not be needed.--- -- However, it is instructive and is left in place.--- jumpN :: (MonadRandom m) => State -> [Generator State] -> Int -> m State--- jumpN s p n = foldM jump s (replicate n p)
− src/ELynx/Simulate/MarkovProcessAlongTree.hs
@@ -1,154 +0,0 @@-{- |- Description : Work with transition probability matrices on rooted trees- Copyright : (c) Dominik Schrempf 2017- License : GPLv3-- Maintainer : dominik.schrempf@gmail.com- Stability : unstable- Portability : non-portable (not tested)--Calculate transition probability matrices, map rate matrices on trees, populate-a tree with states according to a stationary distribution, etc.--The implementation of the Markov process is more than basic and can be improved-in a lot of ways.---}--module ELynx.Simulate.MarkovProcessAlongTree- ( -- * Single rate matrix.- simulate- , simulateAndFlatten- -- * Mixture models.- , simulateMixtureModel- , simulateAndFlattenMixtureModel- , simulateAndFlattenMixtureModelPar- )- where--import Control.Monad-import Control.Monad.Primitive-import Control.Parallel.Strategies-import Data.Tree-import Numeric.LinearAlgebra-import System.Random.MWC (Gen, GenIO)-import System.Random.MWC.Distributions (categorical)--import ELynx.Data.MarkovProcess.RateMatrix-import ELynx.Data.Tree.MeasurableTree-import ELynx.Simulate.MarkovProcess-import ELynx.Tools.Concurrent--toProbTree :: (Measurable a) => RateMatrix -> Tree a -> Tree ProbMatrix-toProbTree q = fmap (probMatrix q . getLen)--getRootStates :: PrimMonad m- => Int -> StationaryDistribution -> Gen (PrimState m) -> m [State]-getRootStates n d g = replicateM n $ categorical d g---- | Simulate a number of sites for a given substitution model. Only the states--- at the leafs are retained. The states at internal nodes are removed. This has--- a lower memory footprint.-simulateAndFlatten :: (PrimMonad m, Measurable a)- => Int -> StationaryDistribution -> ExchangeabilityMatrix -> Tree a -> Gen (PrimState m) -> m [[State]]-simulateAndFlatten n d e t g = do- let q = fromExchangeabilityMatrix e d- pt = toProbTree q t- is <- getRootStates n d g- simulateAndFlatten' is pt g---- This is the heart of the simulation. Take a tree and a list of root states.--- Recursively jump down the branches to the leafs. Forget states at internal.-simulateAndFlatten' :: (PrimMonad m)- => [State] -> Tree ProbMatrix -> Gen (PrimState m) -> m [[State]]-simulateAndFlatten' is (Node p f) g = do- is' <- mapM (\i -> jump i p g) is- if null f- then return [is']- else concat <$> sequence [simulateAndFlatten' is' t g | t <- f]---- | Simulate a number of sites for a given substitution model. Keep states at--- internal nodes. The result is a tree with the list of simulated states as--- node labels.-simulate :: (PrimMonad m, Measurable a)- => Int -> StationaryDistribution -> ExchangeabilityMatrix -> Tree a -> Gen (PrimState m) -> m (Tree [State])-simulate n d e t g = do- let q = fromExchangeabilityMatrix e d- pt = toProbTree q t- is <- getRootStates n d g- simulate' is pt g---- This is the heart of the simulation. Take a tree and a list of root states.--- Recursively jump down the branches to the leafs.-simulate' :: (PrimMonad m)- => [State] -> Tree ProbMatrix -> Gen (PrimState m) -> m (Tree [State])-simulate' is (Node p f) g = do- is' <- mapM (\i -> jump i p g) is- f' <- sequence [simulate' is' t g | t <- f]- return $ Node is' f'--toProbTreeMixtureModel :: (Measurable a)- => [RateMatrix] -> Tree a -> Tree [ProbMatrix]-toProbTreeMixtureModel qs =- fmap (\a -> [probMatrix q . getLen $ a | q <- qs] `using` parList rpar)--getComponentsAndRootStates :: PrimMonad m- => Int -> Vector R -> [StationaryDistribution] -> Gen (PrimState m) -> m ([Int], [State])-getComponentsAndRootStates n ws ds g = do- cs <- replicateM n $ categorical ws g- is <- sequence [ categorical (ds !! c) g | c <- cs ]- return (cs, is)---- | Simulate a number of sites for a given set of substitution models with--- corresponding weights. Forget states at internal nodes. See also--- 'simulateAndFlatten'.-simulateAndFlattenMixtureModel :: (PrimMonad m, Measurable a)- => Int -> Vector R -> [StationaryDistribution] -> [ExchangeabilityMatrix] -> Tree a- -> Gen (PrimState m) -> m [[State]]-simulateAndFlattenMixtureModel n ws ds es t g = do- let qs = zipWith fromExchangeabilityMatrix es ds- pt = toProbTreeMixtureModel qs t- (cs, is) <- getComponentsAndRootStates n ws ds g- simulateAndFlattenMixtureModel' is cs pt g--simulateAndFlattenMixtureModel' :: (PrimMonad m)- => [State] -> [Int] -> Tree [ProbMatrix] -> Gen (PrimState m) -> m [[State]]-simulateAndFlattenMixtureModel' is cs (Node ps f) g- = do is' <- sequence [ jump i (ps !! c) g | (i, c) <- zip is cs ]- if null f- then return [is']- else concat <$> sequence [ simulateAndFlattenMixtureModel' is' cs t g | t <- f ]---- | See 'simulateAndFlattenMixtureModel', parallel version;--- needs to be run in IO monad.-simulateAndFlattenMixtureModelPar- :: Measurable a- => Int -> Vector R -> [StationaryDistribution] -> [ExchangeabilityMatrix] -> Tree a- -> GenIO -> IO [[[State]]]-simulateAndFlattenMixtureModelPar n ws ds es t g = do- let qs = zipWith fromExchangeabilityMatrix es ds- pt = toProbTreeMixtureModel qs t- parComp n (\n' g' -> do- (cs, is) <- getComponentsAndRootStates n' ws ds g'- simulateAndFlattenMixtureModel' is cs pt g') g---- | Simulate a number of sites for a given set of substitution models with--- corresponding weights. Keep states at internal nodes. See also--- 'simulate'.-simulateMixtureModel :: (PrimMonad m, Measurable a)- => Int -> Vector R -> [StationaryDistribution] -> [ExchangeabilityMatrix] -> Tree a- -> Gen (PrimState m) -> m (Tree [State])-simulateMixtureModel n ws ds es t g = do- let qs = zipWith fromExchangeabilityMatrix es ds- pt = toProbTreeMixtureModel qs t- (cs, is) <- getComponentsAndRootStates n ws ds g- simulateMixtureModel' is cs pt g---- See 'simulateAlongProbTree', only we have a number of mixture components. The--- starting states and the components for each site have to be provided.-simulateMixtureModel' :: (PrimMonad m)- => [State] -> [Int] -> Tree [ProbMatrix] -> Gen (PrimState m) -> m (Tree [State])-simulateMixtureModel' is cs (Node ps f) g = do- is' <- sequence [ jump i (ps !! c) g | (i, c) <- zip is cs ]- f' <- sequence [ simulateMixtureModel' is' cs t g | t <- f ]- return $ Node is' f'
test/ELynx/Data/Alphabet/DistributionDiversitySpec.hs view
@@ -1,7 +1,7 @@ {- | Module : ELynx.Data.Alphabet.DistributionDiversitySpec-Copyright : (c) Dominik Schrempf 2019-License : GPL-3+Copyright : (c) Dominik Schrempf 2020+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -12,14 +12,15 @@ -} module ELynx.Data.Alphabet.DistributionDiversitySpec- (spec- ) where+ ( spec+ )+where -import qualified Data.Vector.Unboxed as V+import qualified Data.Vector.Unboxed as V import Test.Hspec import ELynx.Data.Alphabet.DistributionDiversity-import ELynx.Tools.Equality+import ELynx.Tools testArr1 :: V.Vector Double testArr1 = V.replicate 20 0.0@@ -33,28 +34,28 @@ spec :: Spec spec = do- describe "entropy" $- it "calculates entropy of vectors" $ do+ describe "entropy" $ it "calculates entropy of vectors" $ do entropy testArr1 `shouldBe` 0.0 entropy testArr2 `shouldBe` 0.0 entropy testArr3 `shouldSatisfy` nearlyEq 0.9773805948045555 - describe "kEffEntropy" $- it "calculates the effective number of used states using entropy" $ do- kEffEntropy testArr1 `shouldBe` 1.0- kEffEntropy testArr2 `shouldBe` 1.0- kEffEntropy testArr3 `shouldSatisfy` nearlyEq 2.6574860842252765+ describe "kEffEntropy"+ $ it "calculates the effective number of used states using entropy"+ $ do+ kEffEntropy testArr1 `shouldBe` 1.0+ kEffEntropy testArr2 `shouldBe` 1.0+ kEffEntropy testArr3 `shouldSatisfy` nearlyEq 2.6574860842252765 - describe "homoplasy" $- it "calculates homoplasy of vectors" $ do+ describe "homoplasy" $ it "calculates homoplasy of vectors" $ do homoplasy testArr1 `shouldBe` 0.0 homoplasy testArr2 `shouldBe` 1.0 homoplasy testArr3 `shouldSatisfy` nearlyEq 0.74 - describe "kEffHomoplasy" $- it "calculates the effective number of used states using homoplasy" $ do- kEffHomoplasy testArr1 `shouldSatisfy` isInfinite- kEffHomoplasy testArr2 `shouldBe` 1.0- kEffHomoplasy testArr3 `shouldSatisfy` nearlyEq 1.3513513513513513+ describe "kEffHomoplasy"+ $ it "calculates the effective number of used states using homoplasy"+ $ do+ kEffHomoplasy testArr1 `shouldSatisfy` isInfinite+ kEffHomoplasy testArr2 `shouldBe` 1.0+ kEffHomoplasy testArr3 `shouldSatisfy` nearlyEq 1.3513513513513513
− test/ELynx/Data/MarkovProcess/AminoAcidSpec.hs
@@ -1,154 +0,0 @@-{- |-Module : ELynx.Data.MarkovProcess.AminoAcidSpec-Copyright : (c) Dominik Schrempf 2019-License : GPL-3--Maintainer : dominik.schrempf@gmail.com-Stability : unstable-Portability : portable--Creation date: Tue Jan 29 10:47:40 2019.---}--module ELynx.Data.MarkovProcess.AminoAcidSpec- (spec) where--import Numeric.LinearAlgebra-import Test.Hspec--import ELynx.Data.MarkovProcess.AminoAcid-import qualified ELynx.Data.MarkovProcess.RateMatrix as R-import qualified ELynx.Data.MarkovProcess.SubstitutionModel as S-import ELynx.Tools.Equality-import ELynx.Tools.Vector--statDistLGPython :: R.StationaryDistribution-statDistLGPython = normalizeSumVec 1.0 $- fromList [ 0.079066, 0.012937, 0.053052, 0.071586, 0.042302, 0.057337, 0.022355- , 0.062157, 0.064600, 0.099081, 0.022951, 0.041977, 0.044040, 0.040767- , 0.055941, 0.061197, 0.053287, 0.069147, 0.012066, 0.034155]--exchLGPython :: R.ExchangeabilityMatrix-exchLGPython = fromLists- [ [ 0.0000000e+00, 2.4890840e+00, 3.9514400e-01, 1.0385450e+00, 2.5370100e-01- , 2.0660400e+00, 3.5885800e-01, 1.4983000e-01, 5.3651800e-01, 3.9533700e-01- , 1.1240350e+00, 2.7681800e-01, 1.1776510e+00, 9.6989400e-01, 4.2509300e-01- , 4.7271820e+00, 2.1395010e+00, 2.5478700e+00, 1.8071700e-01, 2.1895900e-01]- , [ 2.4890840e+00, 0.0000000e+00, 6.2556000e-02, 3.4990000e-03, 1.1052510e+00- , 5.6926500e-01, 6.4054300e-01, 3.2062700e-01, 1.3266000e-02, 5.9400700e-01- , 8.9368000e-01, 5.2876800e-01, 7.5382000e-02, 8.4808000e-02, 5.3455100e-01- , 2.7844780e+00, 1.1434800e+00, 1.9592910e+00, 6.7012800e-01, 1.1655320e+00]- , [ 3.9514400e-01, 6.2556000e-02, 0.0000000e+00, 5.2438700e+00, 1.7416000e-02- , 8.4492600e-01, 9.2711400e-01, 1.0690000e-02, 2.8295900e-01, 1.5076000e-02- , 2.5548000e-02, 5.0761490e+00, 3.9445600e-01, 5.2338600e-01, 1.2395400e-01- , 1.2402750e+00, 4.2586000e-01, 3.7967000e-02, 2.9890000e-02, 1.3510700e-01]- , [ 1.0385450e+00, 3.4990000e-03, 5.2438700e+00, 0.0000000e+00, 1.8811000e-02- , 3.4884700e-01, 4.2388100e-01, 4.4265000e-02, 1.8071770e+00, 6.9673000e-02- , 1.7373500e-01, 5.4171200e-01, 4.1940900e-01, 4.1285910e+00, 3.6397000e-01- , 6.1197300e-01, 6.0454500e-01, 2.4503400e-01, 7.7852000e-02, 1.2003700e-01]- , [ 2.5370100e-01, 1.1052510e+00, 1.7416000e-02, 1.8811000e-02, 0.0000000e+00- , 8.9586000e-02, 6.8213900e-01, 1.1127270e+00, 2.3918000e-02, 2.5926920e+00- , 1.7988530e+00, 8.9525000e-02, 9.4464000e-02, 3.5855000e-02, 5.2722000e-02- , 3.6181900e-01, 1.6500100e-01, 6.5468300e-01, 2.4571210e+00, 7.8039020e+00]- , [ 2.0660400e+00, 5.6926500e-01, 8.4492600e-01, 3.4884700e-01, 8.9586000e-02- , 0.0000000e+00, 3.1148400e-01, 8.7050000e-03, 2.9663600e-01, 4.4261000e-02- , 1.3953800e-01, 1.4376450e+00, 1.9696100e-01, 2.6795900e-01, 3.9019200e-01- , 1.7399900e+00, 1.2983600e-01, 7.6701000e-02, 2.6849100e-01, 5.4679000e-02]- , [ 3.5885800e-01, 6.4054300e-01, 9.2711400e-01, 4.2388100e-01, 6.8213900e-01- , 3.1148400e-01, 0.0000000e+00, 1.0888200e-01, 6.9726400e-01, 3.6631700e-01- , 4.4247200e-01, 4.5092380e+00, 5.0885100e-01, 4.8135050e+00, 2.4266010e+00- , 9.9001200e-01, 5.8426200e-01, 1.1901300e-01, 5.9705400e-01, 5.3068340e+00]- , [ 1.4983000e-01, 3.2062700e-01, 1.0690000e-02, 4.4265000e-02, 1.1127270e+00- , 8.7050000e-03, 1.0888200e-01, 0.0000000e+00, 1.5906900e-01, 4.1450670e+00- , 4.2736070e+00, 1.9150300e-01, 7.8281000e-02, 7.2854000e-02, 1.2699100e-01- , 6.4105000e-02, 1.0337390e+00, 1.0649107e+01, 1.1166000e-01, 2.3252300e-01]- , [ 5.3651800e-01, 1.3266000e-02, 2.8295900e-01, 1.8071770e+00, 2.3918000e-02- , 2.9663600e-01, 6.9726400e-01, 1.5906900e-01, 0.0000000e+00, 1.3750000e-01- , 6.5660400e-01, 2.1450780e+00, 3.9032200e-01, 3.2342940e+00, 6.3260670e+00- , 7.4868300e-01, 1.1368630e+00, 1.8520200e-01, 4.9906000e-02, 1.3193200e-01]- , [ 3.9533700e-01, 5.9400700e-01, 1.5076000e-02, 6.9673000e-02, 2.5926920e+00- , 4.4261000e-02, 3.6631700e-01, 4.1450670e+00, 1.3750000e-01, 0.0000000e+00- , 6.3123580e+00, 6.8427000e-02, 2.4906000e-01, 5.8245700e-01, 3.0184800e-01- , 1.8228700e-01, 3.0293600e-01, 1.7027450e+00, 6.1963200e-01, 2.9964800e-01]- , [ 1.1240350e+00, 8.9368000e-01, 2.5548000e-02, 1.7373500e-01, 1.7988530e+00- , 1.3953800e-01, 4.4247200e-01, 4.2736070e+00, 6.5660400e-01, 6.3123580e+00- , 0.0000000e+00, 3.7100400e-01, 9.9849000e-02, 1.6725690e+00, 4.8413300e-01- , 3.4696000e-01, 2.0203660e+00, 1.8987180e+00, 6.9617500e-01, 4.8130600e-01]- , [ 2.7681800e-01, 5.2876800e-01, 5.0761490e+00, 5.4171200e-01, 8.9525000e-02- , 1.4376450e+00, 4.5092380e+00, 1.9150300e-01, 2.1450780e+00, 6.8427000e-02- , 3.7100400e-01, 0.0000000e+00, 1.6178700e-01, 1.6957520e+00, 7.5187800e-01- , 4.0083580e+00, 2.0006790e+00, 8.3688000e-02, 4.5376000e-02, 6.1202500e-01]- , [ 1.1776510e+00, 7.5382000e-02, 3.9445600e-01, 4.1940900e-01, 9.4464000e-02- , 1.9696100e-01, 5.0885100e-01, 7.8281000e-02, 3.9032200e-01, 2.4906000e-01- , 9.9849000e-02, 1.6178700e-01, 0.0000000e+00, 6.2429400e-01, 3.3253300e-01- , 1.3381320e+00, 5.7146800e-01, 2.9650100e-01, 9.5131000e-02, 8.9613000e-02]- , [ 9.6989400e-01, 8.4808000e-02, 5.2338600e-01, 4.1285910e+00, 3.5855000e-02- , 2.6795900e-01, 4.8135050e+00, 7.2854000e-02, 3.2342940e+00, 5.8245700e-01- , 1.6725690e+00, 1.6957520e+00, 6.2429400e-01, 0.0000000e+00, 2.8079080e+00- , 1.2238280e+00, 1.0801360e+00, 2.1033200e-01, 2.3619900e-01, 2.5733600e-01]- , [ 4.2509300e-01, 5.3455100e-01, 1.2395400e-01, 3.6397000e-01, 5.2722000e-02- , 3.9019200e-01, 2.4266010e+00, 1.2699100e-01, 6.3260670e+00, 3.0184800e-01- , 4.8413300e-01, 7.5187800e-01, 3.3253300e-01, 2.8079080e+00, 0.0000000e+00- , 8.5815100e-01, 5.7898700e-01, 1.7088700e-01, 5.9360700e-01, 3.1444000e-01]- , [ 4.7271820e+00, 2.7844780e+00, 1.2402750e+00, 6.1197300e-01, 3.6181900e-01- , 1.7399900e+00, 9.9001200e-01, 6.4105000e-02, 7.4868300e-01, 1.8228700e-01- , 3.4696000e-01, 4.0083580e+00, 1.3381320e+00, 1.2238280e+00, 8.5815100e-01- , 0.0000000e+00, 6.4722790e+00, 9.8369000e-02, 2.4886200e-01, 4.0054700e-01]- , [ 2.1395010e+00, 1.1434800e+00, 4.2586000e-01, 6.0454500e-01, 1.6500100e-01- , 1.2983600e-01, 5.8426200e-01, 1.0337390e+00, 1.1368630e+00, 3.0293600e-01- , 2.0203660e+00, 2.0006790e+00, 5.7146800e-01, 1.0801360e+00, 5.7898700e-01- , 6.4722790e+00, 0.0000000e+00, 2.1881580e+00, 1.4082500e-01, 2.4584100e-01]- , [ 2.5478700e+00, 1.9592910e+00, 3.7967000e-02, 2.4503400e-01, 6.5468300e-01- , 7.6701000e-02, 1.1901300e-01, 1.0649107e+01, 1.8520200e-01, 1.7027450e+00- , 1.8987180e+00, 8.3688000e-02, 2.9650100e-01, 2.1033200e-01, 1.7088700e-01- , 9.8369000e-02, 2.1881580e+00, 0.0000000e+00, 1.8951000e-01, 2.4931300e-01]- , [ 1.8071700e-01, 6.7012800e-01, 2.9890000e-02, 7.7852000e-02, 2.4571210e+00- , 2.6849100e-01, 5.9705400e-01, 1.1166000e-01, 4.9906000e-02, 6.1963200e-01- , 6.9617500e-01, 4.5376000e-02, 9.5131000e-02, 2.3619900e-01, 5.9360700e-01- , 2.4886200e-01, 1.4082500e-01, 1.8951000e-01, 0.0000000e+00, 3.1518150e+00]- , [ 2.1895900e-01, 1.1655320e+00, 1.3510700e-01, 1.2003700e-01, 7.8039020e+00- , 5.4679000e-02, 5.3068340e+00, 2.3252300e-01, 1.3193200e-01, 2.9964800e-01- , 4.8130600e-01, 6.1202500e-01, 8.9613000e-02, 2.5733600e-01, 3.1444000e-01- , 4.0054700e-01, 2.4584100e-01, 2.4931300e-01, 3.1518150e+00, 0.0000000e+00]]--statDistUniform :: R.StationaryDistribution-statDistUniform = vector $ replicate 20 0.05--statDistLG :: R.StationaryDistribution-statDistLG = S.stationaryDistribution lg--exchLG :: R.ExchangeabilityMatrix-exchLG = S.exchangeabilityMatrix lg--rmLG :: R.RateMatrix-rmLG = S.rateMatrix lg--spec :: Spec-spec = do- describe "statDistLG" $- it "matches distribution from python library" $- statDistLG `nearlyEqVec` statDistLGPython `shouldBe` True-- describe "exchLG" $- it "matches exchangeability matrix from python library" $ do- exchLG `shouldSatisfy` nearlyEqMatWith 1e-4 exchLGPython- exchLG `nearlyEqMat` rmLG `shouldBe` False-- describe "lg" $- it "stationary distribution can be extracted" $- nearlyEqVecWith 1e-4 (R.getStationaryDistribution rmLG) statDistLG `shouldBe` True-- describe "lgCustom" $- it "stationary distribution can be recovered" $ do- let f = R.getStationaryDistribution $ S.rateMatrix $ lgCustom Nothing statDistUniform- f `nearlyEqVec` statDistUniform `shouldBe` True-- describe "poisson" $- it "stationary distribution is uniform 1/20" $- R.getStationaryDistribution (S.rateMatrix poisson) `nearlyEqVec` statDistUniform `shouldBe` True-- describe "poissonCustom" $- it "stationary distribution can be recovered" $ do- let f = R.getStationaryDistribution $ S.rateMatrix $ poissonCustom Nothing statDistLGPython- f `nearlyEqVec` statDistLGPython `shouldBe` True
− test/ELynx/Data/MarkovProcess/NucleotideSpec.hs
@@ -1,36 +0,0 @@-{- |-Module : ELynx.Data.MarkovProcess.NucleotideSpec-Copyright : (c) Dominik Schrempf 2019-License : GPL-3--Maintainer : dominik.schrempf@gmail.com-Stability : unstable-Portability : portable--Creation date: Fri Jan 25 16:47:28 2019.---}--module ELynx.Data.MarkovProcess.NucleotideSpec- (spec) where--import Data.Vector.Generic-import Test.Hspec--import ELynx.Data.MarkovProcess.Nucleotide-import ELynx.Data.MarkovProcess.RateMatrix-import ELynx.Data.MarkovProcess.SubstitutionModel-import ELynx.Tools.Equality (nearlyEqVec)--stationaryDist :: StationaryDistribution-stationaryDist = fromList [0.2, 0.3, 0.3, 0.2]--hkyModel :: SubstitutionModel-hkyModel = hky 6.0 stationaryDist--spec :: Spec-spec =- describe "getStationaryDistribution" $- it "extracts the stationary distribution from a rate matrix" $ do- let sd = getStationaryDistribution (rateMatrix hkyModel)- sd `nearlyEqVec` stationaryDist `shouldBe` True
test/ELynx/Data/Sequence/AlignmentSpec.hs view
@@ -1,7 +1,7 @@ {- | Module : ELynx.Data.Sequence.AlignmentSpec Copyright : (c) Dominik Schrempf 2018-License : GPL-3+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -12,7 +12,9 @@ -} module ELynx.Data.Sequence.AlignmentSpec- (spec) where+ ( spec+ )+where import qualified Data.ByteString.Lazy.Char8 as L import qualified Data.Matrix.Unboxed as M@@ -22,19 +24,25 @@ import ELynx.Data.Alphabet.Character import ELynx.Data.Sequence.Alignment import ELynx.Import.Sequence.Fasta-import ELynx.Tools.InputOutput-import Files+import ELynx.Tools +fastaNucleotideIUPACFN :: FilePath+fastaNucleotideIUPACFN = "data/NucleotideIUPAC.fasta"+ ssData :: M.Matrix Character-ssData = M.fromLists $ map (reverse . map fromChar) [ "AAA", "GAA", "TAA" ]+ssData = M.fromLists $ map (reverse . map fromChar) ["AAA", "GAA", "TAA"] ssA :: Alignment-ssA = Alignment (map L.pack ["SEQUENCE_1", "SEQUENCE_2", "SEQUENCE_3"]) DNAI ssData+ssA = Alignment (map L.pack ["SEQUENCE_1", "SEQUENCE_2", "SEQUENCE_3"])+ (replicate 3 L.empty)+ DNAI+ ssData spec :: Spec-spec =- describe "subSample" $- it "correctly sub sample an Alignment" $ do- a <- either error id . fromSequences <$> parseFileWith (fasta DNAI) fastaNucleotideIUPACFN- let ss = subSample [0,3,5] a- ss `shouldBe` ssA+spec = describe "subSample" $ it "correctly sub sample an Alignment" $ do+ a <-+ either error id+ . fromSequences+ <$> parseFileWith (fasta DNAI) fastaNucleotideIUPACFN+ let ss = subSample [0, 3, 5] a+ ss `shouldBe` ssA
test/ELynx/Data/Sequence/SequenceSpec.hs view
@@ -1,7 +1,7 @@ {- | Module : ELynx.Data.Sequence.SequenceSpec Copyright : (c) Dominik Schrempf 2018-License : GPL-3+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -12,33 +12,42 @@ -} module ELynx.Data.Sequence.SequenceSpec- (spec) where+ ( spec+ )+where -import qualified Data.ByteString.Lazy.Char8 as L+import qualified Data.ByteString.Lazy.Char8 as L import Test.Hspec import ELynx.Data.Alphabet.Alphabet import ELynx.Data.Sequence.Sequence import ELynx.Import.Sequence.Fasta-import ELynx.Tools.InputOutput-import Files+import ELynx.Tools +fastaDifferentLengthFN :: FilePath+fastaDifferentLengthFN = "data/NucleotideDifferentLength.fasta"++fastaDifferentLengthTrimmedFN :: FilePath+fastaDifferentLengthTrimmedFN = "data/NucleotideDifferentLengthTrimmed.fasta"+ longestSequenceInFileBS :: L.ByteString-longestSequenceInFileBS = L.unlines $ map L.pack [ ">SEQUENCE_3"- , "ATTTAAAAAAACCCAAAACCCGGGCCCCGGGTTTTTTTA" ]+longestSequenceInFileBS = L.unlines+ $ map L.pack [">SEQUENCE_3", "ATTTAAAAAAACCCAAAACCCGGGCCCCGGGTTTTTTTA"] longestSequenceInFile :: Sequence-longestSequenceInFile = parseByteStringWith "Fasta byte string" (fastaSequence DNA) longestSequenceInFileBS+longestSequenceInFile = parseByteStringWith "Fasta byte string"+ (fastaSequence DNA)+ longestSequenceInFileBS spec :: Spec spec = do- describe "longest" $- it "finds the longest sequence"$ do+ describe "longest" $ it "finds the longest sequence" $ do ss <- parseFileWith (fasta DNA) fastaDifferentLengthFN longest ss `shouldBe` longestSequenceInFile - describe "filterLongerThan" $- it "filters sequences that are longer than a specified length" $ do- ss <- parseFileWith (fasta DNA) fastaDifferentLengthFN- ss' <- parseFileWith (fasta DNA) fastaDifferentLengthTrimmedFN- filterLongerThan 10 ss `shouldBe` ss'+ describe "filterLongerThan"+ $ it "filters sequences that are longer than a specified length"+ $ do+ ss <- parseFileWith (fasta DNA) fastaDifferentLengthFN+ ss' <- parseFileWith (fasta DNA) fastaDifferentLengthTrimmedFN+ filterLongerThan 10 ss `shouldBe` ss'
test/ELynx/Data/Sequence/TranslateSpec.hs view
@@ -1,7 +1,7 @@ {- | Module : ELynx.Data.Sequence.TranslateSpec Copyright : (c) Dominik Schrempf 2018-License : GPL-3+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -12,7 +12,9 @@ -} module ELynx.Data.Sequence.TranslateSpec- (spec) where+ ( spec+ )+where import Test.Hspec @@ -20,13 +22,17 @@ import ELynx.Data.Character.Codon import ELynx.Data.Sequence.Translate import ELynx.Import.Sequence.Fasta-import ELynx.Tools.InputOutput-import Files+import ELynx.Tools +fastaTranslateDNAFN :: FilePath+fastaTranslateDNAFN = "data/TranslateMitochondrialVertebrateDNA.fasta"++fastaTranslateProteinFN :: FilePath+fastaTranslateProteinFN = "data/TranslateMitochondrialVertebrateProtein.fasta"+ spec :: Spec spec =- describe "translateDNAX" $- it "correctly translates a test sequence" $ do+ describe "translateDNAX" $ it "correctly translates a test sequence" $ do ss <- parseFileWith (fasta DNAX) fastaTranslateDNAFN ss' <- parseFileWith (fasta ProteinS) fastaTranslateProteinFN map (translateSeq VertebrateMitochondrial 0) ss `shouldBe` ss'
test/ELynx/Export/Sequence/FastaSpec.hs view
@@ -1,7 +1,7 @@ {- | Module : ELynx.Export.Sequence.FastaSpec-Copyright : (c) Dominik Schrempf 2019-License : GPL-3+Copyright : (c) Dominik Schrempf 2020+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -12,21 +12,30 @@ -} module ELynx.Export.Sequence.FastaSpec- (spec) where+ ( spec+ )+where +import Test.Hspec+ import ELynx.Data.Alphabet.Alphabet import ELynx.Export.Sequence.Fasta import ELynx.Import.Sequence.Fasta-import ELynx.Tools.InputOutput-import Files-import Test.Hspec+import ELynx.Tools +fastaNucleotideIUPACFN :: FilePath+fastaNucleotideIUPACFN = "data/NucleotideIUPAC.fasta"+ spec :: Spec spec =- describe "sequencesToFasta" $- it "should create a fasta bytestring that, when parsed again, is the original sequence" $ do- ss <- parseFileWith (fasta DNAI) fastaNucleotideIUPACFN- let f = sequencesToFasta ss- ss' = parseByteStringWith- "Fasta byte string created from sequence object" (fasta DNAI) f- ss `shouldBe` ss'+ describe "sequencesToFasta"+ $ it+ "should create a fasta bytestring that, when parsed again, is the original sequence"+ $ do+ ss <- parseFileWith (fasta DNAI) fastaNucleotideIUPACFN+ let f = sequencesToFasta ss+ ss' = parseByteStringWith+ "Fasta byte string created from sequence object"+ (fasta DNAI)+ f+ ss `shouldBe` ss'
− test/ELynx/Import/MarkovProcess/EDMModelPhylobayesSpec.hs
@@ -1,55 +0,0 @@-{- |-Module : ELynx.Import.MarkovProcess.EDMModelPhylobayesSpec-Copyright : (c) Dominik Schrempf 2019-License : GPL-3--Maintainer : dominik.schrempf@gmail.com-Stability : unstable-Portability : portable--Creation date: Tue Jan 29 12:26:49 2019.---}--module ELynx.Import.MarkovProcess.EDMModelPhylobayesSpec- (spec) where--import Numeric.LinearAlgebra (fromList)-import Test.Hspec--import ELynx.Import.MarkovProcess.EDMModelPhylobayes-import ELynx.Tools.InputOutput (parseFileWith)-import Files (componentsFilePhylobayes)--testComponents :: [EDMComponent]-testComponents =- [ (0.3690726059430788 , fromList- [ 0.0746116859730418, 0.0069967050701822, 0.0792378063875672, 0.1103113280153337- , 0.0084753541579630, 0.0420976864270153, 0.0359995156499732, 0.0176053658394712- , 0.1140305648643406, 0.0306656386818212, 0.0143923751236750, 0.0805907035564168- , 0.0136556223619922, 0.0854829959418090, 0.0907856579607629, 0.0955855590305961- , 0.0585231905811694, 0.0276404860443685, 0.0024024774767569, 0.0109092808557433 ])- , (0.28019163862430846, fromList- [ 0.0379912838025543, 0.0085917039746034, 0.0016977963056655, 0.0029780993213938- , 0.0369086689176780, 0.0051768704164023, 0.0027496797757250, 0.2650202435999584- , 0.0035114356263963, 0.2705514660371282, 0.0634852940497682, 0.0036512748214048- , 0.0030666766077233, 0.0046946399741140, 0.0039700770877914, 0.0137110414055372- , 0.0286159990065166, 0.2351259204392804, 0.0018121721806795, 0.0066896566496796 ])- , (0.23699880225859807, fromList- [ 0.2699288544670601, 0.0262695884363299, 0.0119030286238640, 0.0143694121662917- , 0.0113341089668980, 0.1028437949687162, 0.0097376980401260, 0.0250662077996463- , 0.0137176088824149, 0.0302792626622677, 0.0151697996735645, 0.0255041240872322- , 0.0194279535358810, 0.0155053062932145, 0.0146597721555602, 0.2261838872893110- , 0.1040919090386891, 0.0534187313840882, 0.0024899475832214, 0.0080990039456234 ])- , (0.11373695317401472, fromList- [ 0.0325183213724537, 0.0153142424565790, 0.0072297704620506, 0.0101526783552517- , 0.2682049848838039, 0.0107298905683658, 0.0420960604309369, 0.0414751553493385- , 0.0131985876449809, 0.1152120909184748, 0.0300072241461393, 0.0173787316286536- , 0.0062421326435569, 0.0173124967325860, 0.0181776303531701, 0.0328778350287987- , 0.0250382111285986, 0.0437257936268565, 0.0258103804738129, 0.2272977817955918 ]) ]--spec :: Spec-spec = describe "phylobayes" $- it "parses a text file with stationary distributions in phylobayes format" $ do- cs <- parseFileWith phylobayes componentsFilePhylobayes- cs `shouldBe` testComponents
− test/ELynx/Import/MarkovProcess/SiteprofilesPhylobayesSpec.hs
@@ -1,46 +0,0 @@-{- |-Module : ELynx.Import.MarkovProcess.SiteprofilesPhylobayesSpec-Copyright : (c) Dominik Schrempf 2019-License : GPL-3--Maintainer : dominik.schrempf@gmail.com-Stability : unstable-Portability : portable--Creation date: Wed Sep 18 09:51:35 2019.---}--module ELynx.Import.MarkovProcess.SiteprofilesPhylobayesSpec- ( spec- ) where--import qualified Data.Vector.Storable as V-import Test.Hspec--import ELynx.Data.MarkovProcess.EDMModel-import ELynx.Import.MarkovProcess.SiteprofilesPhylobayes-import ELynx.Tools.Equality-import ELynx.Tools.InputOutput--fn :: FilePath-fn = "data/HSSPMany.siteprofiles"--getProfiles :: IO [EDMComponent]-getProfiles = parseFileWith siteprofiles fn--firstProfile :: V.Vector Double-firstProfile = V.fromList- [ 0.0267009, 0.013874, 0.022432, 0.0440502, 0.0485174, 0.0407515- , 0.0170806, 0.348043, 0.0371379, 0.0715536, 0.0454168, 0.0342213- , 0.0146872, 0.0239681, 0.0768379, 0.0210387, 0.0123336- , 0.0512678, 0.0149093, 0.0351776 ]--spec :: Spec-spec = describe "import phylobayes siteprofiles" $- it "parses a text file with siteprofiles in phylobayes format" $ do- profiles <- getProfiles- length profiles `shouldBe` 701- map fst profiles `shouldBe` replicate 701 1.0- map (V.sum . snd) profiles `shouldSatisfy` nearlyEqListWith 1e-5 (replicate 701 1.0)- snd (head profiles) `shouldBe` firstProfile
test/ELynx/Import/Sequence/FastaSpec.hs view
@@ -1,7 +1,7 @@ {- | Module : ELynx.Import.Sequence.FastaSpec-Copyright : (c) Dominik Schrempf 2019-License : GPL-3+Copyright : (c) Dominik Schrempf 2020+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -12,39 +12,62 @@ -} module ELynx.Import.Sequence.FastaSpec- (spec) where+ ( spec+ )+where import Data.Either+import Test.Hspec+ import ELynx.Data.Alphabet.Alphabet import qualified ELynx.Data.Sequence.Alignment as M import ELynx.Import.Sequence.Fasta-import ELynx.Tools.InputOutput-import Files-import Test.Hspec+import ELynx.Tools +fastaNucleotideFN :: FilePath+fastaNucleotideFN = "data/Nucleotide.fasta"++fastaNucleotideIUPACFN :: FilePath+fastaNucleotideIUPACFN = "data/NucleotideIUPAC.fasta"++fastaErroneousFN :: FilePath+fastaErroneousFN = "data/Erroneous.fasta"++fastaAminoAcidFN :: FilePath+fastaAminoAcidFN = "data/AminoAcid.fasta"+ spec :: Spec-spec =- describe "fastaFileAlignment" $ do- it "parses a fasta file with nucleotide sequences with equal length" $ do- a <- either error id . M.fromSequences <$> parseFileWith (fasta DNA) fastaNucleotideFN- M.nSequences a `shouldBe` 3- M.length a `shouldBe` 40+spec = describe "fastaFileAlignment" $ do+ it "parses a fasta file with nucleotide sequences with equal length" $ do+ a <-+ either error id+ . M.fromSequences+ <$> parseFileWith (fasta DNA) fastaNucleotideFN+ M.nSequences a `shouldBe` 3+ M.length a `shouldBe` 40 - it "parses a fasta file with nucleotide IUPAC sequences with equal length" $ do- a <- either error id . M.fromSequences <$> parseFileWith (fasta DNAI) fastaNucleotideIUPACFN- M.nSequences a `shouldBe` 3- M.length a `shouldBe` 40+ it "parses a fasta file with nucleotide IUPAC sequences with equal length"+ $ do+ a <-+ either error id+ . M.fromSequences+ <$> parseFileWith (fasta DNAI) fastaNucleotideIUPACFN+ M.nSequences a `shouldBe` 3+ M.length a `shouldBe` 40 - it "should not parse erroneous files" $ do- ea <- runParserOnFile (fasta DNAI) fastaErroneousFN- ea `shouldSatisfy` isLeft+ it "should not parse erroneous files" $ do+ ea <- runParserOnFile (fasta DNAI) fastaErroneousFN+ ea `shouldSatisfy` isLeft - it "parses a fasta file with amino acid sequences with equal length" $ do- a <- either error id . M.fromSequences <$> parseFileWith (fasta Protein) fastaAminoAcidFN- M.nSequences a `shouldBe` 2- M.length a `shouldBe` 237+ it "parses a fasta file with amino acid sequences with equal length" $ do+ a <-+ either error id+ . M.fromSequences+ <$> parseFileWith (fasta Protein) fastaAminoAcidFN+ M.nSequences a `shouldBe` 2+ M.length a `shouldBe` 237 - it "should not parse erroneous files" $ do- a <- runParserOnFile (fasta ProteinI) fastaErroneousFN- a `shouldSatisfy` isLeft+ it "should not parse erroneous files" $ do+ a <- runParserOnFile (fasta ProteinI) fastaErroneousFN+ a `shouldSatisfy` isLeft
− test/ELynx/Simulate/MarkovProcessAlongTreeSpec.hs
@@ -1,59 +0,0 @@-{- |-Module : ELynx.Simulate.MarkovProcessAlongTreeSpec-Copyright : (c) Dominik Schrempf 2019-License : GPL-3--Maintainer : dominik.schrempf@gmail.com-Stability : unstable-Portability : portable--Creation date: Thu Jan 24 15:25:49 2019.---}--module ELynx.Simulate.MarkovProcessAlongTreeSpec- (spec) where--import Data.Tree-import System.Random.MWC-import Test.Hspec--import ELynx.Data.MarkovProcess.Nucleotide-import qualified ELynx.Data.MarkovProcess.SubstitutionModel as S-import ELynx.Data.Tree.PhyloTree-import ELynx.Simulate.MarkovProcess-import ELynx.Simulate.MarkovProcessAlongTree---- sampleNewickText :: T.Text--- sampleNewickText = T.pack "(Aeropyrum0:0.5478645225,(((((((((Arabidopsi:0.0701001024,Oryza_sati:0.0765988261):0.0309636193,Gymnosperm:0.0520325624):0.0338982245,Physcomitr:0.0768008916):0.0895714685,(Chlamydomo:0.1136227755,Dunaliella:0.1406347323):0.1117340620):0.0818876186,Rhodophyta:0.3405656487):0.0363527066,((((((Babesia_bo:0.1646969208,Theileria0:0.1519889486):0.1908081096,Plasmodium:0.3250696762):0.0637865908,(Toxoplasma:0.1153570425,Eimeria000:0.1671916078):0.0980136930):0.0518956330,Cryptospor:0.3175062809):0.1607708388,Ciliophora:0.5687502950):0.0624078848,(Phytophtho:0.2016424948,((Thalassios:0.1202730781,Phaeodacty:0.1290341329):0.1772775509,Phaeophyce:0.1989260715):0.0312359673):0.1154768302):0.0311952864):0.0149160316,(((((((((Candida_al:0.1027755272,Saccharomy:0.1190206560):0.1333487870,Neurospora:0.1977309079):0.0522926266,Schizosacc:0.2019603227):0.0567441011,(Cryptococc:0.1948614959,Ustilago_m:0.1564451295):0.0775729694):0.0323959951,Glomus_int:0.1573670796):0.0194701292,Chytridiom:0.2228415254):0.0384370601,Encephalit:1.4622174644):0.0416231688,(((Drosophila:0.2160627753,(Mammalians:0.1080484094,Tunicates0:0.1739253014):0.0289624371):0.0346633757,Hydrozoa00:0.2058137032):0.0480963050,Monosiga_b:0.3020637584):0.0654894239):0.0380915725,(Dictyostel:0.3453588998,Mastigamoe:0.3844779231):0.0478795653):0.0129578395):1.7592083381,((Archaeoglo:0.5402784445,Methanococ:0.4088567459):0.0993669265,Pyrococcus:0.4058713829):0.1734405968):0.2193511807,Pyrobaculu:0.7507718047):0.1646616482,Sulfolobus:0.5404967897);"--testTree :: Tree (PhyloLabel Int)-testTree = Node (PhyloLabel 0 Nothing (Just 0.0)) [ Node (PhyloLabel 1 Nothing (Just 1.0)) [], Node (PhyloLabel 2 Nothing (Just 1.0)) [] ]---- testStateTree :: Tree State--- testStateTree = Node 0 [Node 3 [], Node 2 []]--nullTree :: Tree (PhyloLabel Int)-nullTree = Node (PhyloLabel 0 Nothing (Just 0.0)) []--nullStateTree :: Tree State-nullStateTree = fmap (const 0) nullTree---- uniformStationaryDistribution :: Int -> Vector R--- uniformStationaryDistribution n = fromList $ replicate n (1.0 / fromIntegral n)--spec :: Spec-spec =- describe "simulateNSitesAlongTree" $ do- it "simulates one site along an easy tree" $ do- gen <- create- tr <- simulate 1 d e nullTree gen- fmap head tr `shouldBe` nullStateTree-- it "simulates some sites along a harder tree" $ do- gen <- create- tr <- simulate 10 d e testTree gen- (length . rootLabel $ tr) `shouldBe` 10- where d = S.stationaryDistribution jc- e = S.exchangeabilityMatrix jc-
− test/Files.hs
@@ -1,55 +0,0 @@-{- |-Module : Files-Description : File names of test data-Copyright : (c) Dominik Schrempf 2019-License : GPL-3--Maintainer : dominik.schrempf@gmail.com-Stability : unstable-Portability : portable--Creation date: Fri Jan 18 10:03:02 2019.---}--module Files- ( fastaNucleotideFN- , fastaNucleotideIUPACFN- , fastaAminoAcidFN- , fastaDifferentLengthFN- , fastaDifferentLengthTrimmedFN- , fastaErroneousFN- , componentsFilePhylobayes- , fastaTranslateDNAFN- , fastaTranslateProteinFN- ) where--dataDir :: FilePath-dataDir = "data/"--fastaNucleotideFN :: FilePath-fastaNucleotideFN = dataDir ++ "Nucleotide.fasta"--fastaNucleotideIUPACFN :: FilePath-fastaNucleotideIUPACFN = dataDir ++ "NucleotideIUPAC.fasta"--fastaAminoAcidFN :: FilePath-fastaAminoAcidFN = dataDir ++ "AminoAcid.fasta"--fastaDifferentLengthFN :: FilePath-fastaDifferentLengthFN = dataDir ++ "NucleotideDifferentLength.fasta"--fastaDifferentLengthTrimmedFN :: FilePath-fastaDifferentLengthTrimmedFN = dataDir ++ "NucleotideDifferentLengthTrimmed.fasta"--fastaErroneousFN :: FilePath-fastaErroneousFN = dataDir ++ "Erroneous.fasta"--componentsFilePhylobayes :: FilePath-componentsFilePhylobayes = dataDir ++ "EDMDistsPhylobayes.txt"--fastaTranslateDNAFN :: FilePath-fastaTranslateDNAFN = dataDir ++ "TranslateMitochondrialVertebrateDNA.fasta"--fastaTranslateProteinFN :: FilePath-fastaTranslateProteinFN = dataDir ++ "TranslateMitochondrialVertebrateProtein.fasta"