elynx-seq 0.0.1 → 0.1.0

raw patch · 28 files changed

+795/−650 lines, 28 filesdep −lensdep ~elynx-toolsdep ~elynx-treedep ~math-functionsPVP ok

version bump matches the API change (PVP)

Dependencies removed: lens

Dependency ranges changed: elynx-tools, elynx-tree, math-functions

API changes (from Hackage documentation)

Files

• ChangeLog.md +2/−2
• README.md +88/−78
• elynx-seq.cabal +14/−13
• src/ELynx/Data/Alphabet/Alphabet.hs +9/−9
• src/ELynx/Data/Alphabet/DistributionDiversity.hs +3/−0
• src/ELynx/Data/Character/BoundaryMutation.hs +4/−4
• src/ELynx/Data/MarkovProcess/GammaRateHeterogeneity.hs +2/−3
• src/ELynx/Data/MarkovProcess/MixtureModel.hs +40/−31
• src/ELynx/Data/MarkovProcess/Nucleotide.hs +5/−0
• src/ELynx/Data/MarkovProcess/PhyloModel.hs +1/−2
• src/ELynx/Data/MarkovProcess/SubstitutionModel.hs +23/−23
• src/ELynx/Data/Sequence/Alignment.hs +244/−0
• src/ELynx/Data/Sequence/Defaults.hs +12/−12
• src/ELynx/Data/Sequence/MultiSequenceAlignment.hs +0/−253
• src/ELynx/Data/Sequence/Sequence.hs +89/−97
• src/ELynx/Data/Sequence/Translate.hs +2/−2
• src/ELynx/Export/Sequence/Fasta.hs +1/−2
• src/ELynx/Import/MarkovProcess/SiteprofilesPhylobayes.hs +77/−0
• src/ELynx/Import/Sequence/Fasta.hs +3/−2
• src/ELynx/Simulate/MarkovProcess.hs +2/−3
• src/ELynx/Simulate/MarkovProcessAlongTree.hs +58/−42
• test/ELynx/Data/MarkovProcess/AminoAcidSpec.hs +6/−7
• test/ELynx/Data/MarkovProcess/NucleotideSpec.hs +1/−1
• test/ELynx/Data/Sequence/AlignmentSpec.hs +40/−0
• test/ELynx/Data/Sequence/MultiSequenceAlignmentSpec.hs +0/−40
• test/ELynx/Import/MarkovProcess/SiteprofilesPhylobayesSpec.hs +46/−0
• test/ELynx/Import/Sequence/FastaSpec.hs +15/−15
• test/ELynx/Simulate/MarkovProcessAlongTreeSpec.hs +8/−9

@@ -1,10 +1,10 @@  # Table of Contents -1.  [Changelog for ELynx](#org9916021)+1.  [Changelog for ELynx](#org263dcac)  -<a id="org9916021"></a>+<a id="org263dcac"></a>  # Changelog for ELynx 

@@ -2,33 +2,33 @@  # The ELynx Suite -Version: 0.0.1.+Version: 0.1.0. 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 code). The command line with all arguments is consistently,-and automatically logged. The work overhead in the beginning usually pays off in-the end.+(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.  The Elynx Suite consists of three library packages and two executables providing a range of sub commands.  The library packages are: --   **elynx-seq:** Handle evolutionary sequences and multi sequence alignments;--   **elynx-tree:** Handle phylogenetic trees;--   **elynx-tools:** Tools for the provided executables;+-   **elynx-seq:** Handle evolutionary sequences and multi sequence alignments.+-   **elynx-tree:** Handle phylogenetic trees.+-   **elynx-tools:** Tools for the provided executables.  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 (FASTA format).+-   **tlynx:** Analyze, modify, and simulate phylogenetic trees (Newick format). -**ELynx is still under development. We happily receive comments, ideas, feature-requests, or pull requests!**+**ELynx is actively developed. We happily receive comments, ideas, feature+requests, and pull requests!**   # Installation@@ -57,8 +57,8 @@              stack install     -    The binaries are installed into `~/.local/bin/` which has to be added [PATH](https://en.wikipedia.org/wiki/PATH_(variable)).-    Then, they can be used directly.+    The binaries are installed into `~/.local/bin/` which has to be added to the+    [PATH](https://en.wikipedia.org/wiki/PATH_(variable)) environment variable. Then, they can be used directly.   # SLynx@@ -67,8 +67,8 @@      slynx --help -    ELynx Suite version 0.5.1. Developed by Dominik Schrempf. Compiled on September-    9, 2019, at 10:48 am, UTC.+    ELynx Suite version 0.1.0. Developed by Dominik Schrempf. Compiled on January+    30, 2020, at 09:18 am, UTC.          Usage: slynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME] COMMAND       Analyze, and simulate multi sequence alignments.@@ -82,17 +82,14 @@                                Specify base name of output file          Available commands:-      concatenate              -      examine                  If data is a multi sequence alignment, additionally-                               analyze columns.-      filter-rows              -      filter-columns           -      simulate                 -      sub-sample               Create a given number of multi sequence alignments,-                               each of which containing a given number of random-                               sites drawn from the original multi sequence-                               alignment.-      translate                +      concatenate              Concatenate sequences found in input files.+      examine                  Examine sequences. If data is a multi sequence+                               alignment, additionally analyze columns.+      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:       - FASTA@@ -106,15 +103,20 @@       - ProteinS (amino acids; including gaps, and translation stops)       - ProteinI (amino acids; including gaps, translation stops, and IUPAC codes)     -    The ELynx Suite.+    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. The work overhead in the beginning usually pays off in the end.+    logged.+         slynx     Analyze, modify, and simulate evolutionary sequences.     tlynx     Analyze, modify, and simulate phylogenetic trees.+    +    Get help for specific commands:+      slynx examine --help   ## Concatenate@@ -123,9 +125,8 @@      slynx concatenate --help -    Concatenate sequences found in input files.-         Usage: slynx concatenate (-a|--alphabet NAME) INPUT-FILE+      Concatenate sequences found in input files.          Available options:       -a,--alphabet NAME       Specify alphabet type NAME@@ -139,10 +140,9 @@      slynx examine --help -    Examine sequences.-         Usage: slynx examine (-a|--alphabet NAME) [INPUT-FILE] [--per-site]-      If data is a multi sequence alignment, additionally analyze columns.+      Examine sequences. If data is a multi sequence alignment, additionally analyze+      columns.          Available options:       -a,--alphabet NAME       Specify alphabet type NAME@@ -157,26 +157,27 @@      slynx filter-rows --help -    Filter rows (or sequences) found in input files.-         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:       -a,--alphabet NAME       Specify alphabet type NAME       INPUT-FILE               Read sequences from INPUT-FILE       --longer-than LENGTH     Only keep sequences longer than LENGTH       --shorter-than LENGTH    Only keep sequences shorter than LENGTH+      --standard-characters    Only keep sequences containing at least one standard+                               (i.e., non-IUPAC) character       -h,--help                Show this help text  Filter columns of multi sequence alignments with `filter-columns`.      slynx filter-columns --help -    Filter columns of multi-sequence alignments.-         Usage: slynx filter-columns (-a|--alphabet NAME) [INPUT-FILE]                                 [--standard-chars DOUBLE]+      Filter columns of multi-sequence alignments.          Available options:       -a,--alphabet NAME       Specify alphabet type NAME@@ -192,13 +193,13 @@      slynx simulate --help -    Simulate multi sequence alignments.-         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@@ -210,6 +211,8 @@                                are shown below (mutually exclusive with -s option)       -e,--edm-file NAME       Empirical distribution model file NAME in Phylobayes                                format+      -p,--siteprofile-files NAMES+                               File names of site profiles in Phylobayes format       -w,--mixture-model-weights "[DOUBLE,DOUBLE,...]"                                Weights of mixture model components       -g,--gamma-rate-heterogeneity "(NCAT,SHAPE)"@@ -257,8 +260,7 @@     Usage: slynx sub-sample (-a|--alphabet NAME) [INPUT-FILE]                             (-n|--number-of-sites INT)                             (-m|--number-of-alignments INT) [-S|--seed [INT]]-      Create a given number of multi sequence alignments, each of which containing a-      given number of random sites drawn from the original multi sequence alignment.+      Sub-sample columns from multi sequence alignments.          Available options:       -a,--alphabet NAME       Specify alphabet type NAME@@ -269,6 +271,9 @@       -S,--seed [INT]          Seed for random number generator; list of 32 bit                                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.   ## Translate@@ -277,10 +282,9 @@      slynx translate --help -    Translate from DNA to Protein or DNAX to ProteinX.-         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:       -a,--alphabet NAME       Specify alphabet type NAME@@ -297,8 +301,8 @@      tlynx --help -    ELynx Suite version 0.5.1. Developed by Dominik Schrempf. Compiled on September-    9, 2019, at 10:48 am, UTC.+    ELynx Suite version 0.1.0. Developed by Dominik Schrempf. Compiled on January+    30, 2020, at 09:18 am, UTC.          Usage: tlynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME] COMMAND       Compare, examine, and simulate phylogenetic trees.@@ -312,26 +316,34 @@                                Specify base name of output file          Available commands:-      compare                  -      examine                  -      simulate                 Simulate reconstructed trees using the point process.-                               See Gernhard, T. (2008). The conditioned-                               reconstructed process. Journal of Theoretical-                               Biology, 253(4), 769–778.-                               http://doi.org/10.1016/j.jtbi.2008.04.005+      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).          File formats:       - Newick     -    The ELynx Suite.+    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. The work overhead in the beginning usually pays off in the end.+    logged.+         slynx     Analyze, modify, and simulate evolutionary sequences.     tlynx     Analyze, modify, and simulate phylogenetic trees.+    +    Get help for specific commands:+      slynx examine --help   ## Compare@@ -340,24 +352,18 @@      tlynx compare --help -    Compute distances between phylogenetic trees.-    -    Usage: tlynx compare (-d|--distance MEASURE) [-s|--summary-statistics]-                         [INPUT-FILES]+    Usage: tlynx compare [-n|--normalize] [-b|--bipartitions] [-i|--newick-iqtree]+                         NAME+      Compare two phylogenetic trees (compute distances and branch-wise+      differences).          Available options:-      -d,--distance MEASURE    Type of distance to calculate (available distance-                               measures are listed below)-      -s,--summary-statistics  Report summary statistics only-      INPUT-FILES              Read tree(s) from INPUT-FILES; if more files are-                               given, one tree is expected per file+      -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)+      NAME                     Tree file       -h,--help                Show this help text-    -    Available distance measures:-      Symmetric distance: -d symmetric-      Incompatible split distance: -d incompatible-split[VAL]-        Collapse branches with support less than VAL before distance calculation;-        in this way, only well supported difference contribute to the distance measure.   ## Examine@@ -366,12 +372,13 @@      tlynx examine --help -    Compute summary statistics of phylogenetic trees.-    -    Usage: tlynx examine [INPUT-FILE]+    Usage: tlynx examine [INPUT-FILE] [-i|--newick-iqtree]+      Compute summary statistics of phylogenetic trees.          Available options:       INPUT-FILE               Read trees from INPUT-FILE+      -i,--newick-iqtree       Use IQ-TREE Newick format (internal node labels are+                               branch support values)       -h,--help                Show this help text  @@ -381,15 +388,11 @@      tlynx simulate --help -    Simulate phylogenetic trees using birth and death processes.-         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 reconstructed trees using the point process. See Gernhard, T. (2008).-      The conditioned reconstructed process. Journal of Theoretical Biology, 253(4),-      769–778. http://doi.org/10.1016/j.jtbi.2008.04.005+      Simulate phylogenetic trees using birth and death processes.          Available options:       -t,--nTrees INT          Number of trees (default: 10)@@ -409,9 +412,16 @@     Height of Trees: if no tree height is given, the heights will be randomly drawn from the expected distribution given the number of leaves, the birth and the death rate.     Summary statistics only: only print (NumberOfExtantChildren BranchLength) pairs for each branch of each tree. The trees are separated by a newline character.     Sub-sampling: simulate one big tree with n'=round(n/rho), n'>=n, leaves, and randomly sample sub-trees with n leaves. Hence, with rho=1.0, the same tree is reported over and over again.+    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 -Documentation of the library can be found on Hackage.+Documentation of the library can be found on [Hackage](https://hackage.haskell.org/):++-   [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)+-   [slynx](https://hackage.haskell.org/package/slynx)+-   [tlynx](https://hackage.haskell.org/package/tlynx) 

@@ -1,6 +1,6 @@ cabal-version: 1.12 name: elynx-seq-version: 0.0.1+version: 0.1.0 license: GPL-3 license-file: LICENSE copyright: Dominik Schrempf (2019)@@ -46,13 +46,14 @@         ELynx.Data.MarkovProcess.PhyloModel         ELynx.Data.MarkovProcess.RateMatrix         ELynx.Data.MarkovProcess.SubstitutionModel+        ELynx.Data.Sequence.Alignment         ELynx.Data.Sequence.Defaults-        ELynx.Data.Sequence.MultiSequenceAlignment         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@@ -60,19 +61,18 @@     other-modules:         Paths_elynx_seq     default-language: Haskell2010-    ghc-options: -Wall+    ghc-options: -Wall -O2     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.0.1 && <0.1,-        elynx-tree >=0.0.1 && <0.1,+        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,-        lens >=4.17.1 && <4.18,-        math-functions >=0.3.1.0 && <0.4,+        math-functions >=0.3.2.1 && <0.4,         matrices >=0.5.0 && <0.6,         megaparsec >=7.0.5 && <7.1,         mwc-random >=0.14.0.0 && <0.15,@@ -93,17 +93,19 @@         ELynx.Data.Alphabet.DistributionDiversitySpec         ELynx.Data.MarkovProcess.AminoAcidSpec         ELynx.Data.MarkovProcess.NucleotideSpec-        ELynx.Data.Sequence.MultiSequenceAlignmentSpec+        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 -eventlog -threaded -rtsopts -with-rtsopts=-N+    ghc-options: -Wall -O2 -eventlog -threaded -rtsopts+                 -with-rtsopts=-N     build-depends:         async >=2.2.2 && <2.3,         base >=4.12.0.0 && <4.13,@@ -111,14 +113,13 @@         containers >=0.6.0.1 && <0.7,         data-memocombinators >=0.5.1 && <0.6,         elynx-seq -any,-        elynx-tools >=0.0.1 && <0.1,-        elynx-tree >=0.0.1 && <0.1,+        elynx-tools >=0.1.0 && <0.2,+        elynx-tree >=0.1.0 && <0.2,         hmatrix >=0.20.0.0 && <0.21,         hspec >=2.7.1 && <2.8,         hspec-megaparsec >=2.0.1 && <2.1,         integration >=0.2.1 && <0.3,-        lens >=4.17.1 && <4.18,-        math-functions >=0.3.1.0 && <0.4,+        math-functions >=0.3.2.1 && <0.4,         matrices >=0.5.0 && <0.6,         megaparsec >=7.0.5 && <7.1,         mwc-random >=0.14.0.0 && <0.15,

@@ -34,7 +34,7 @@     Alphabet (..)   , AlphabetSpec (..)   , alphabetSpec-  , alphabetNameVerbose+  , description   , isStd   , isGap   , isUnknown@@ -53,14 +53,14 @@               deriving (Show, Read, Eq, Ord, Enum, Bounded)  -- | Verbose alphabet name.-alphabetNameVerbose :: Alphabet -> String-alphabetNameVerbose DNA      = "DNA (nucleotides)"-alphabetNameVerbose DNAX     = "DNAX (nucleotides; including gaps)"-alphabetNameVerbose DNAI     = "DNAI (nucleotides; including gaps, and IUPAC codes)"-alphabetNameVerbose Protein  = "Protein (amino acids)"-alphabetNameVerbose ProteinX = "ProteinX (amino acids; including gaps)"-alphabetNameVerbose ProteinS = "ProteinS (amino acids; including gaps, and translation stops)"-alphabetNameVerbose ProteinI = "ProteinI (amino acids; including gaps, translation stops, and IUPAC codes)"+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)"  -- | Alphabet specification. 'S.Set' is used because it provides fast lookups. data AlphabetSpec = AlphabetSpec {

@@ -84,6 +84,9 @@   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 alph d = V.map (`saveDivision` s) counts   where

@@ -12,7 +12,8 @@  * Changelog -TODO: BM states can not be read and written with a single letter, like characters.+TODO: BM states can not be read and written with a single letter, like+characters. Check status of implementation.  -} @@ -34,7 +35,6 @@   , neighbors   ) where -import           Control.Lens import qualified Data.ByteString.Lazy.Char8      as L import           Numeric.SpecFunctions           (choose) @@ -140,12 +140,12 @@   | 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' - p^._1 + 1) (p^._2) (p^._3)+  | 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 ..]]-        p = last $ takeWhile (\e -> e^._1 <= i') l+        (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.

@@ -22,7 +22,6 @@   , expand   ) where -import           Control.Lens import qualified Data.ByteString.Lazy.Char8                 as L import           Numeric.Integration.TanhSinh import           Statistics.Distribution@@ -65,13 +64,13 @@ expandSubstitutionModel :: Int -> Double -> S.SubstitutionModel -> M.MixtureModel expandSubstitutionModel n alpha sm = M.fromSubstitutionModels name (repeat 1.0) sms   where-    name = sm ^. S.name <> getName n alpha+    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 = mm ^. M.name <> getName n alpha+    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 ..]

@@ -1,5 +1,3 @@-{-# LANGUAGE TemplateHaskell #-}- {- | Module      :  ELynx.Data.MarkovProcess.MixtureModel Description :  Mixture models are a set of substitution models with weights@@ -21,7 +19,7 @@     Weight   , Component (Component)   , MixtureModel (MixtureModel)-    -- * Lenses and other accessors+    -- * Getters   , name   , getAlphabet   , getWeights@@ -40,7 +38,6 @@   , summarize   ) where -import           Control.Lens import qualified Data.ByteString.Builder                    as L import qualified Data.ByteString.Lazy.Char8                 as L @@ -53,37 +50,32 @@  -- | A mixture model component has a weight and a substitution model. data Component = Component-  { _weight     :: Weight-  , _substModel :: S.SubstitutionModel+  { weight     :: Weight+  , substModel :: S.SubstitutionModel   }   deriving (Show, Read) -makeLenses ''Component- -- | A mixture model with its components. data MixtureModel = MixtureModel-  { _name       :: S.Name-  , _components :: [Component]+  { name       :: S.Name        -- ^ Name+  , components :: [Component]   }   deriving (Show, Read) -makeLenses ''MixtureModel---- | Get alphabet used with mixture model. Throws error if components use different--- 'Alphabet's.+-- | 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)-            then head $ mm ^.. components . traverse . substModel . S.alphabet+            then S.alphabet . substModel $ head (components mm)             else error "Mixture model is invalid."  -- | Get weights. getWeights :: MixtureModel -> [Weight]-getWeights m = m ^.. components . traverse . weight+getWeights = map weight . components  -- | Get substitution models. getSubstitutionModels :: MixtureModel -> [S.SubstitutionModel]-getSubstitutionModels m = m ^.. components . traverse . substModel+getSubstitutionModels = map substModel . components   -- | Create a mixture model from a list of substitution models.@@ -93,11 +85,17 @@  -- | Concatenate mixture models. concatenate :: S.Name -> [MixtureModel] -> MixtureModel-concatenate n mms = MixtureModel n $ concatMap (view components) mms+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 = over (components . traverse . substModel) (S.scale s)+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.@@ -107,9 +105,15 @@         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 = over (components . traverse . substModel) (S.appendName n)+appendName n m = m {components = cs'}+  where cs  = components m+        cs' = map (appendNameComponent n) cs  -- | Checks if a mixture model is valid. --@@ -118,19 +122,24 @@ -- and providing an algebraic way of creating mixture models (empty and -- addComponent which performs necessary checks). isValid :: MixtureModel -> Bool-isValid mm = not (null $ mm ^. components)-                         && allEqual alphabets-  where alphabets = mm ^.. components . traverse . substModel . S.alphabet+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 mmc =-  L.pack "Weight: " <> (L.toLazyByteString . L.doubleDec $ mmc ^. weight)-  : S.summarize (mmc ^. substModel)+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 mm =-  L.pack ("Mixture model: " ++ mm ^. name ++ ".")-  : concat [ L.pack ("Component " ++ show i ++ ":") : summarizeComponent c-            | (i, c) <- zip [1 :: Int ..] (mm ^. components) ]+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 []

@@ -34,6 +34,11 @@  -- data PhyloModel = MixtureModel | SubstitutionModel +--+-- I think it could simply be:+-- data PhyloModel = [(Weight, SubstitutionModel)]+--+ -- data MixtureModel = [(Weight, SubstitutionModel)]  -- data SubstitutionModel = SMDNA DNASubstitutionModel | SMAA AASubstitutionModel

@@ -24,7 +24,6 @@   , summarize   ) where -import           Control.Lens import qualified Data.ByteString.Lazy.Char8                 as L  import           ELynx.Data.Alphabet.Alphabet@@ -39,7 +38,7 @@ -- | Extract code from phylogenetic model. getAlphabet :: PhyloModel -> Alphabet getAlphabet (MixtureModel mm)      = M.getAlphabet mm-getAlphabet (SubstitutionModel sm) = sm ^. S.alphabet+getAlphabet (SubstitutionModel sm) = S.alphabet sm  -- | Summarize a phylogenetic model; lines to be printed to screen or log. summarize :: PhyloModel -> [L.ByteString]

@@ -1,5 +1,3 @@-{-# LANGUAGE TemplateHaskell #-}- {- | Module      :  ELynx.Data.MarkovProcess.SubstitutionModel Description :  Data type describing substitution model@@ -26,7 +24,7 @@   , name   , stationaryDistribution   , exchangeabilityMatrix-  , getRateMatrix+  , rateMatrix   , totalRate   -- * Building substitution models   , substitutionModel@@ -39,7 +37,6 @@   , summarize   ) where -import           Control.Lens import qualified Data.ByteString.Lazy.Char8          as L import qualified Numeric.LinearAlgebra               as LinAlg @@ -56,25 +53,26 @@ type Params = [Double]  -- | Complete definition of a substitution model. Create instances with--- 'substitutionModel'.+-- '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-  , _name                   :: Name-  , _params                 :: Params-  , _stationaryDistribution :: R.StationaryDistribution-  , _exchangeabilityMatrix  :: R.ExchangeabilityMatrix+  { alphabet               :: Alphabet -- ^ Alphabet+  , name                   :: Name     -- ^ Name+  , params                 :: Params   -- ^ List of parameters+  , stationaryDistribution :: R.StationaryDistribution -- ^ Stationary distribution+  , exchangeabilityMatrix  :: R.ExchangeabilityMatrix  -- ^ Exchangeability matrix   }   deriving (Show, Read) -makeLenses ''SubstitutionModel- -- | Calculate rate matrix from substitution model.-getRateMatrix :: SubstitutionModel -> R.RateMatrix-getRateMatrix sm = R.fromExchangeabilityMatrix (sm^.exchangeabilityMatrix) (sm^.stationaryDistribution)+rateMatrix :: SubstitutionModel -> R.RateMatrix+rateMatrix sm = R.fromExchangeabilityMatrix (exchangeabilityMatrix sm) (stationaryDistribution sm)  -- | Get scale of substitution model. totalRate :: SubstitutionModel -> Double-totalRate sm = R.totalRate (sm^.stationaryDistribution) (getRateMatrix sm)+totalRate sm = R.totalRate (stationaryDistribution sm) (rateMatrix sm)  -- | Create normalized 'SubstitutionModel'. See 'normalize'. substitutionModel :: Alphabet -> Name -> Params@@ -90,7 +88,8 @@  -- | Scale the rate of a substitution model by given factor. scale :: Double -> SubstitutionModel -> SubstitutionModel-scale r = over exchangeabilityMatrix (LinAlg.scale r)+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.@@ -100,19 +99,20 @@  -- | Abbend to name. appendName :: Name -> SubstitutionModel -> SubstitutionModel-appendName n = over name (<> n)+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 (sm^.alphabet) ++ " substitution model: " ++ sm^.name ++ ".") :-  [ "Parameters: " ++ show (sm^.params) ++ "." | not (null (sm^.params))] ++-  case sm^.alphabet of-    DNA -> [ "Stationary distribution: " ++ dispv precision (sm^.stationaryDistribution) ++ "."-           , "Exchangeability matrix:\n" ++ dispmi 2 precision (sm^.exchangeabilityMatrix) ++ "."+  (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 (sm^.stationaryDistribution) ++ "."+    Protein -> [ "Stationary distribution: " ++ dispv precision (stationaryDistribution sm) ++ "."                , "Scale: " ++ show (roundN precision $ totalRate sm) ++ "."                ]     _ -> error "Extended character sets are not supported with substitution models."

@@ -0,0 +1,244 @@+{- |+Module      :  ELynx.Data.Sequence.Alignment+Description :  Multi sequence alignment related types and functions+Copyright   :  (c) Dominik Schrempf 2018+License     :  GPL-3++Maintainer  :  dominik.schrempf@gmail.com+Stability   :  unstable++Portability :  portable++Creation date: Thu Oct  4 18:40:18 2018.++This module is to be imported qualified.++-}+++module ELynx.Data.Sequence.Alignment+  ( Alignment (..)+  , length+  , nSequences+  -- | * Input, output+  , fromSequences+  , toSequences+  , summarize+  -- | * Manipulation+  , join+  , concat+  , concatAlignments+  , filterColsOnlyStd+  , filterColsStd+  , filterColsNoGaps+  -- | * Analysis+  , FrequencyData+  , distribution+  , toFrequencyData+  , kEffEntropy+  , kEffHomoplasy+  , countIUPACChars+  , countGaps+  , countUnknowns+  -- | * Sub sample+  , subSample+  , randomSubSample+  ) where++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           System.Random.MWC++import qualified ELynx.Data.Alphabet.Alphabet              as A+import           ELynx.Data.Alphabet.Character+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++-- | A collection of sequences.+data Alignment = Alignment+                              { names    :: [S.Name]+                              , alphabet :: A.Alphabet+                              , matrix   :: M.Matrix Character+                              }+  deriving (Show, Eq)++-- | Number of sites.+length :: Alignment -> Int+length = M.cols . matrix++-- | Number of sequences.+nSequences :: Alignment -> Int+nSequences = M.rows . matrix++-- | 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++-- | 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++-- -- | Show a 'Alignment', untrimmed.+-- summarizeSeq :: Alignment -> Int -> L.ByteString+-- summarizeSeq m i =+--   L.unwords [ alignLeft nameWidth $ names m !! i+--             , summarizeByteString summaryLength $+--               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 ]++-- | Similar to 'S.summarizeSequenceList' but with different Header.+summarize :: Alignment -> L.ByteString+summarize a = header a <> S.body (toSequences a)++-- | Join two 'Alignment's vertically. That is, add more sequences+-- to an alignment. See also 'concat'.+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++-- | Concatenate two 'Alignment's horizontally. That is, add more+-- sites to an alignment. See also 'join'.+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++-- | Concatenate a list of 'Alignment's horizontally. See+-- 'concat'.+concatAlignments :: [Alignment] -> Alignment+concatAlignments []  = error "concatAlignments: Nothing to concatenate."+concatAlignments [a] = a+concatAlignments as  = foldl' concat (head as) (tail as)++-- Only keep columns from alignment that satisfy given predicate.+filterColsWith :: (V.Vector Character -> Bool) -> Alignment -> Alignment+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+-- characters are removed.+filterColsOnlyStd :: Alignment -> Alignment+filterColsOnlyStd a = filterColsWith (V.all $ A.isStd (alphabet a)) a++-- | 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++-- | Only keep columns without gaps or unknown characters.+filterColsNoGaps :: Alignment -> Alignment+filterColsNoGaps a = filterColsWith (V.all $ not . A.isGap (alphabet a)) a++-- | Frequency data; do not store the actual characters, but their frequencies.+-- The matrix is of size @N x K@, where @N@ is the number of sites, and @K@ is+-- the number of characters.+type FrequencyData = M.Matrix Double++-- | Calculcate frequency of characters at each site of a multi sequence alignment.+toFrequencyData :: Alignment -> FrequencyData+toFrequencyData a = fMapColParChunk 100 (D.frequencyCharacters spec) (matrix a)+  where spec = A.alphabetSpec (alphabet a)++-- | Calculate the distribution of characters.+distribution :: FrequencyData -> [Double]+distribution fd = map (/ fromIntegral nSites) $ V.toList $+  foldl1 (V.zipWith (+)) (M.toColumns fd)+  where nSites = M.cols fd++-- | Diversity analysis. See 'kEffEntropy'.+kEffEntropy :: FrequencyData -> [Double]+kEffEntropy fd = parMapChunk chunksize D.kEffEntropy (M.toColumns fd)++-- | Diversity analysis. See 'kEffEntropy'.+kEffHomoplasy :: FrequencyData -> [Double]+kEffHomoplasy fd = parMapChunk chunksize D.kEffHomoplasy (M.toColumns fd)++-- | Count the number of standard (i.e., not extended IUPAC) characters in the+-- alignment.+countIUPACChars :: Alignment -> Int+countIUPACChars a = V.length . V.filter (A.isIUPAC (alphabet a)) $ allChars+  where allChars = M.flatten $ matrix a++-- | Count the number of gaps in the alignment.+countGaps :: Alignment -> Int+countGaps a = V.length . V.filter (A.isGap (alphabet a)) $ allChars+  where allChars = M.flatten $ matrix a++-- | Count the number of unknown characters in the alignment.+countUnknowns :: Alignment -> Int+countUnknowns a = V.length . V.filter (A.isUnknown (alphabet a)) $ allChars+  where allChars = M.flatten $ matrix a++-- | 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++-- | Randomly sample a given number of sites of the multi sequence 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+  return $ subSample is a

@@ -14,24 +14,24 @@   module ELynx.Data.Sequence.Defaults-  ( defSequenceNameWidth-  , defSequenceSummaryLength-  , defSequenceListSummaryNumber-  , defFieldWidth+  ( nameWidth+  , summaryLength+  , summaryNSequences+  , fieldWidth   ) where  -- | Space reserved for sequence names when printing them.-defSequenceNameWidth :: Int-defSequenceNameWidth = 23+nameWidth :: Int+nameWidth = 23  -- | The length shown when summarizing sequences.-defSequenceSummaryLength :: Int-defSequenceSummaryLength = 60+summaryLength :: Int+summaryLength = 60  -- | How many sequences are shown in summary.-defSequenceListSummaryNumber :: Int-defSequenceListSummaryNumber = 200+summaryNSequences :: Int+summaryNSequences = 200  -- | Field width for tables.-defFieldWidth :: Int-defFieldWidth = 13+fieldWidth :: Int+fieldWidth = 13

@@ -1,253 +0,0 @@-{-# LANGUAGE TemplateHaskell #-}--{- |-Module      :  ELynx.Data.MultiSequenceAlignment-Description :  Multi sequence alignment related types and functions-Copyright   :  (c) Dominik Schrempf 2018-License     :  GPL-3--Maintainer  :  dominik.schrempf@gmail.com-Stability   :  unstable--Portability :  portable--Creation date: Thu Oct  4 18:40:18 2018.---}---module ELynx.Data.Sequence.MultiSequenceAlignment-  ( MultiSequenceAlignment (MultiSequenceAlignment)-  , msaLength-  , msaNSequences-  -- | * Input, output-  , fromSequenceList-  , toSequenceList-  , showMSA-  , summarizeMSA-  -- | * Manipulation-  , msaJoin-  , msaConcatenate-  , msasConcatenate-  , filterColumnsOnlyStd-  , filterColumnsStd-  , filterColumnsNoGaps-  -- | * Analysis-  , FrequencyData-  , toFrequencyData-  , kEffEntropy-  , kEffHomoplasy-  , countIUPACChars-  , countGaps-  , countUnknowns-  -- | * Sub sample-  , subSample-  , randomSubSample-  ) where--import           Control.Lens-import           Control.Monad-import           Control.Monad.Primitive-import qualified Data.ByteString.Lazy.Char8                as L-import           Data.List-import qualified Data.Matrix.Unboxed                       as M-import qualified Data.Vector.Unboxed                       as V-import           System.Random.MWC--import qualified ELynx.Data.Alphabet.Alphabet              as A-import           ELynx.Data.Alphabet.Character-import qualified ELynx.Data.Alphabet.DistributionDiversity as D-import           ELynx.Data.Sequence.Defaults-import qualified ELynx.Data.Sequence.Sequence              as S-import           ELynx.Tools.ByteString-import           ELynx.Tools.Concurrent-import           ELynx.Tools.Equality-import           ELynx.Tools.Matrix---- | A collection of sequences.-data MultiSequenceAlignment = MultiSequenceAlignment-                              { _names    :: [S.SequenceName]-                              , _alphabet :: A.Alphabet-                              , _matrix   :: M.Matrix Character-                              }-  deriving (Read, Show, Eq)--makeLenses ''MultiSequenceAlignment---- | Number of sites.-msaLength :: MultiSequenceAlignment -> Int-msaLength = M.cols . view matrix---- | Number of sequences.-msaNSequences :: MultiSequenceAlignment -> Int-msaNSequences = M.rows . view matrix---- | Create 'MultiSequenceAlignment' from a list of 'S.Sequence's.-fromSequenceList :: [S.Sequence] -> Either String MultiSequenceAlignment-fromSequenceList ss-  | S.equalLength ss && allEqual (map (view S.alphabet) ss) =-      Right $ MultiSequenceAlignment ns a d-  | S.equalLength ss =-      Left "Sequences do not have equal codes."-  | otherwise =-      Left "Sequences do not have equal lengths."-  where-    ns   = map (view S.name) ss-    a    = head ss ^. S.alphabet-    bss  = map (view S.characters) ss-    d    = M.fromRows bss---- | Conversion to list of 'S.Sequence's.-toSequenceList :: MultiSequenceAlignment -> [S.Sequence]-toSequenceList (MultiSequenceAlignment ns a d) = zipWith (\n r -> S.Sequence n a r) ns rows-  where-    rows  = M.toRows d--msaHeader :: L.ByteString-msaHeader = L.unwords [ alignLeft defSequenceNameWidth (L.pack "Name")-                      , L.pack "Sequence" ]---- | Show a 'S.Sequence', untrimmed.-showSequence :: MultiSequenceAlignment -> Int -> L.ByteString-showSequence m i =-  L.unwords [ alignLeft defSequenceNameWidth $ (m ^. names) !! i-            , S.fromCharacters $ M.takeRow (m ^. matrix) i ]---- | Show a 'S.Sequence', untrimmed.-summarizeSequence :: MultiSequenceAlignment -> Int -> L.ByteString-summarizeSequence m i =-  L.unwords [ alignLeft defSequenceNameWidth $ (m ^. names) !! i-            , summarizeByteString defSequenceSummaryLength $-              S.fromCharacters $ M.takeRow (m ^. matrix) i ]---- | Show a 'MultiSequenceAlignment' in text form.-showMSA :: MultiSequenceAlignment -> L.ByteString-showMSA msa = L.unlines $ msaHeader :-  map (showSequence msa) [0 .. (msaNSequences msa - 1)]--summarizeMSAHeader :: MultiSequenceAlignment -> L.ByteString-summarizeMSAHeader msa = L.unlines $-  [ L.pack "Multi sequence alignment."-  , L.pack $ "Code: " ++ A.alphabetNameVerbose (msa^.alphabet) ++ "."-  , L.pack $ "Length: " ++ show (msaLength msa) ++ "." ]-  ++ reportLengthSummary ++ reportNumberSummary-  where reportLengthSummary =-          [ L.pack $ "For each sequence, the "-            ++ show defSequenceSummaryLength ++ " first bases are shown."-          | msaLength msa > defSequenceSummaryLength ]-        reportNumberSummary =-          [ L.pack $ show defSequenceListSummaryNumber ++ " out of " ++-            show (msaNSequences msa) ++ " sequences are shown."-          | msaNSequences msa > defSequenceListSummaryNumber ]---- | Similar to 'S.summarizeSequenceList' but with different Header.-summarizeMSA :: MultiSequenceAlignment -> L.ByteString-summarizeMSA msa = L.unlines $ summarizeMSAHeader msa :-  map (summarizeSequence msa) [0 .. n - 1]-  where n = min (msaNSequences msa) defSequenceListSummaryNumber---- | Join two 'MultiSequenceAlignment's vertically. That is, add more sequences--- to an alignment. See also 'msaConcatenate'.-msaJoin :: MultiSequenceAlignment-        -> MultiSequenceAlignment-        -> MultiSequenceAlignment--- top bottom.-msaJoin t b-  | msaLength t == msaLength b-    && t^.alphabet == b^.alphabet-  = MultiSequenceAlignment ns a (tD === bD)-  | otherwise  = error "msaJoin: Multi sequence alignments do not have equal length."-  where-    ns = t^.names ++ b^.names-    tD = t^.matrix-    bD = t^.matrix-    a  = t^.alphabet---- | Concatenate two 'MultiSequenceAlignment's horizontally. That is, add more--- sites to an alignment. See also 'msaJoin'.-msaConcatenate :: MultiSequenceAlignment-               -> MultiSequenceAlignment-               -> MultiSequenceAlignment--- left right.-msaConcatenate l r-  | msaNSequences l == msaNSequences r-    && l^.alphabet == r^.alphabet-  = MultiSequenceAlignment (l ^. names) a (lD ||| rD)-  | otherwise = error "msaConcatenate: Multi sequence alignments do not have equal length."-  where-    lD = l^.matrix-    rD = r^.matrix-    a  = l^.alphabet---- | Concatenate a list of 'MultiSequenceAlignment's horizontally. See--- 'msaConcatenate'.-msasConcatenate :: [MultiSequenceAlignment] -> MultiSequenceAlignment-msasConcatenate []    = error "msasConcatenate: Nothing to concatenate."-msasConcatenate [msa] = msa-msasConcatenate msas  = foldl' msaConcatenate (head msas) (tail msas)---- Only keep columns from alignment that satisfy given predicate.-filterColumnsWith :: (V.Vector Character -> Bool) -> MultiSequenceAlignment -> MultiSequenceAlignment-filterColumnsWith p = over matrix (M.fromColumns . filter p . M.toColumns)---- | Only keep columns with standard characters. Alignment columns with IUPAC--- characters are removed.-filterColumnsOnlyStd :: MultiSequenceAlignment -> MultiSequenceAlignment-filterColumnsOnlyStd msa = filterColumnsWith (V.all $ A.isStd (msa^.alphabet)) msa---- | Filter columns with proportion of standard character larger than given number.-filterColumnsStd :: Double -> MultiSequenceAlignment -> MultiSequenceAlignment-filterColumnsStd prop msa = filterColumnsWith-  (\col -> prop * nSeqs <= fromIntegral (V.length (V.filter (A.isStd a) col)))-  msa-  where a = msa^.alphabet-        nSeqs = fromIntegral $ msaNSequences msa---- | Only keep columns without gaps or unknown characters.-filterColumnsNoGaps :: MultiSequenceAlignment -> MultiSequenceAlignment-filterColumnsNoGaps msa = filterColumnsWith (V.all $ not . A.isGap (msa^.alphabet)) msa---- | Frequency data; do not store the actual characters, but only their--- frequencies.-type FrequencyData = M.Matrix Double---- | Calculcate frequency of characters in multi sequence alignment.-toFrequencyData :: MultiSequenceAlignment -> FrequencyData-toFrequencyData msa = fMapColParChunk 100 (D.frequencyCharacters spec) (msa^.matrix)-  where spec = A.alphabetSpec (msa^.alphabet)---- | Diversity analysis. See 'kEffEntropy'.-kEffEntropy :: FrequencyData -> [Double]-kEffEntropy fd = parMapChunk 500 D.kEffEntropy (M.toColumns fd)---- | Diversity analysis. See 'kEffEntropy'.-kEffHomoplasy :: FrequencyData -> [Double]-kEffHomoplasy fd = parMapChunk 500 D.kEffHomoplasy (M.toColumns fd)---- | Count the number of standard (i.e., not extended IUPAC) characters in the--- alignment.-countIUPACChars :: MultiSequenceAlignment -> Int-countIUPACChars msa = V.length . V.filter (A.isIUPAC (msa^.alphabet)) $ allChars-  where allChars = M.flatten $ msa^.matrix---- | Count the number of gaps in the alignment.-countGaps :: MultiSequenceAlignment -> Int-countGaps msa = V.length . V.filter (A.isGap (msa^.alphabet)) $ allChars-  where allChars = M.flatten $ msa^.matrix---- | Count the number of unknown characters in the alignment.-countUnknowns :: MultiSequenceAlignment -> Int-countUnknowns msa = V.length . V.filter (A.isUnknown (msa^.alphabet)) $ allChars-  where allChars = M.flatten $ msa^.matrix---- | Sample the given sites from a multi sequence alignment.-subSample :: [Int] -> MultiSequenceAlignment -> MultiSequenceAlignment-subSample is = over matrix (subSampleMatrix is)---- | Randomly sample a given number of sites of the multi sequence alignment.-randomSubSample :: PrimMonad m-          => Int -> MultiSequenceAlignment  -> Gen (PrimState m) -> m MultiSequenceAlignment-randomSubSample n msa g = do-  let l = msaLength msa-  is <- replicateM n $ uniformR (0, l-1) g-  return $ subSample is msa

@@ -1,4 +1,4 @@-{-# LANGUAGE TemplateHaskell #-}+{-# LANGUAGE OverloadedStrings #-}  {- | Module      :  ELynx.Data.Sequence@@ -12,46 +12,44 @@  Creation date: Thu Oct  4 18:54:51 2018. +This module is to be imported qualified.+ -}  module ELynx.Data.Sequence.Sequence   ( -- * Types-    SequenceName-  , SequenceCharacters-  , Sequence (Sequence)-  -- * Lenses-  , name-  , alphabet-  , characters+    Name+  , Characters+  , Sequence (..)   -- * Input-  , toCharacters+  , fromByteString   -- * Output-  , fromCharacters-  , showSequence-  , showSequenceList-  , sequenceListHeader-  , summarizeSequence-  , summarizeSequenceList-  , summarizeSequenceListBody+  , toByteString+  , header+  , summarize+  , summarizeSequences+  , body   -- * Analysis-  , lengthSequence+  , length   , equalLength   , longest   -- * Manipulation-  , trimSequence-  , concatenate-  , concatenateSeqs+  , trim+  , concat+  , concatSequences   -- * Filtering   , filterShorterThan   , filterLongerThan+  , filterStandard   ) where -import           Control.Lens import           Control.Parallel.Strategies import qualified Data.ByteString.Lazy.Char8    as L 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 qualified Text.Printf                   as P  import qualified ELynx.Data.Alphabet.Alphabet  as A@@ -60,122 +58,116 @@ import           ELynx.Tools.ByteString import           ELynx.Tools.Equality --- | For now, 'SequenceName's are just 'L.ByteString's.-type SequenceName = L.ByteString+-- | For now, 'Name's are just 'L.ByteString's.+type Name = L.ByteString  -- | The vector of characters of a sequence.-type SequenceCharacters = V.Vector Character---- | Sequences have a name, a code and hopefully a lot of data.-data Sequence = Sequence { _name       :: SequenceName-                         , _alphabet   :: A.Alphabet-                         , _characters :: SequenceCharacters }-  deriving (Eq)--makeLenses ''Sequence+type Characters = V.Vector Character  -- | Convert byte string to sequence characters.-toCharacters :: L.ByteString -> SequenceCharacters-toCharacters = V.fromList . map fromChar . L.unpack+fromByteString :: L.ByteString -> Characters+fromByteString = V.fromList . map fromChar . L.unpack  -- | Convert sequence characters to byte string.-fromCharacters :: SequenceCharacters -> L.ByteString-fromCharacters = L.pack . map toChar . V.toList+toByteString :: Characters -> L.ByteString+toByteString = L.pack . map toChar . V.toList -showInfo :: Sequence -> L.ByteString-showInfo s = L.unwords [ alignLeft defSequenceNameWidth (s^.name)-                       , alignRight defFieldWidth (L.pack $ show $ s^.alphabet)-                       , alignRight defFieldWidth (L.pack . show $ len)-                       , alignRight defFieldWidth (L.pack $ P.printf "%.3f" pGaps) ]-  where len = lengthSequence s+-- | Sequences have a name, a code and hopefully a lot of data.+data Sequence = Sequence { name       :: Name+                         , 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 -instance Show Sequence where-  show s = L.unpack $ showSequence s---- | Show a 'Sequence', untrimmed.-showSequence :: Sequence -> L.ByteString-showSequence s = L.unwords [showInfo s, fromCharacters $ s^.characters]---- | Show a list of 'Sequence's, untrimmed.-showSequenceList :: [Sequence] -> L.ByteString-showSequenceList = L.unlines . map showSequence---- | Header printed before 'Sequence' list.-sequenceListHeader :: L.ByteString-sequenceListHeader = L.unwords [ alignLeft defSequenceNameWidth (L.pack "Name")-                               , alignRight defFieldWidth (L.pack "Code")-                               , alignRight defFieldWidth (L.pack "Length")-                               , alignRight defFieldWidth (L.pack "Gaps [%]")-                               , L.pack "Sequence" ]- -- | Trim and show a 'Sequence'.-summarizeSequence :: Sequence -> L.ByteString-summarizeSequence s = L.unwords [ showInfo s-                                , summarizeByteString defSequenceSummaryLength-                                  (fromCharacters $ s^.characters) ]+summarize :: Sequence -> L.ByteString+summarize s = L.unwords [ getInfo s , summarizeByteString summaryLength $ toByteString (characters s) ]  -- | Trim and show a list of 'Sequence's.-summarizeSequenceList :: [Sequence] -> L.ByteString-summarizeSequenceList ss = summarizeSequenceListHeader ss <>-                           summarizeSequenceListBody (take defSequenceListSummaryNumber ss)+summarizeSequences :: [Sequence] -> L.ByteString+summarizeSequences ss = header ss <> body (take summaryNSequences ss) -summarizeSequenceListHeader :: [Sequence] -> L.ByteString-summarizeSequenceListHeader ss = L.unlines $+-- | Header printed before 'Sequence' list.+tableHeader :: L.ByteString+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 defSequenceSummaryLength ++ " first bases are shown."-  , L.pack $ "List contains " ++ show (length ss) ++ " sequences."-  , L.pack ""-  , sequenceListHeader ]-  where l = length ss-        s = show defSequenceListSummaryNumber ++ " out of " ++-            show (length ss) ++ " sequences are shown."+  [ 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 > defSequenceListSummaryNumber = [L.pack s]+          | l > summaryNSequences = [L.pack s]           | otherwise = []  -- | Trim and show a list of 'Sequence's.-summarizeSequenceListBody :: [Sequence] -> L.ByteString-summarizeSequenceListBody ss = L.unlines (map summarizeSequence ss `using` parListChunk 5 rdeepseq)+body :: [Sequence] -> L.ByteString+body ss = L.unlines (map summarize ss `using` parListChunk 5 rdeepseq)  -- | Calculate length of 'Sequence'.-lengthSequence :: Sequence -> Int-lengthSequence s = fromIntegral $ V.length $ s ^. characters+length :: Sequence -> Int+length = fromIntegral . V.length . characters  -- | Check if all 'Sequence's have equal length. equalLength :: [Sequence] -> Bool-equalLength = allEqual . map lengthSequence+equalLength = allEqual . map length  -- | Find the longest 'Sequence' in a list. longest :: [Sequence] -> Sequence-longest = maximumBy (comparing lengthSequence)+longest = maximumBy (comparing length)  -- | Count number of gaps or unknown characters in sequence. countGaps :: Sequence -> Int-countGaps s = V.length . V.filter (A.isGap (s^.alphabet)) $ s^.characters+countGaps s = V.length . V.filter (A.isGap $ alphabet s) $ characters s  -- | Trim to given length.-trimSequence :: Int -> Sequence -> Sequence-trimSequence n = over characters (V.take $ fromIntegral n)+trim :: Int -> Sequence -> Sequence+trim n (Sequence nm a cs) = Sequence nm a (V.take (fromIntegral n) cs) --- | Concatenate two sequences. 'SequenceName's have to match.-concatenate :: Sequence -> Sequence -> Sequence-concatenate (Sequence i c cs) (Sequence j k ks)+-- | 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 ++ "." --- | Concatenate a list of sequences, see 'concatenate'.-concatenateSeqs :: [[Sequence]] -> [Sequence]-concatenateSeqs []   = error "concatenateSeqs: Nothing to concatenate."-concatenateSeqs [ss] = ss-concatenateSeqs sss  = foldl1 (zipWith concatenate) sss+-- | Concatenate a list of sequences, see 'concat'.+concatSequences :: [[Sequence]] -> [Sequence]+concatSequences []   = error "concatenateSequences: Nothing to concatenate."+concatSequences [ss] = ss+concatSequences sss  = foldl1 (zipWith concat) sss  -- | Only take 'Sequence's that are shorter than a given number. filterShorterThan :: Int -> [Sequence] -> [Sequence]-filterShorterThan n = filter (\x -> lengthSequence x < n)+filterShorterThan n = filter (\x -> length x < n)  -- | Only take 'Sequence's that are longer than a given number. filterLongerThan :: Int -> [Sequence] -> [Sequence]-filterLongerThan n = filter (\x -> lengthSequence x > n)+filterLongerThan n = filter (\x -> length x > n)++-- | Only take 'Sequence's that contain at least on non-IUPAC character.+filterStandard :: [Sequence] -> [Sequence]+filterStandard = filter (\s -> anyStandard (alphabet s) s)++-- Are all characters IUPAC characters?+anyStandard :: A.Alphabet -> Sequence -> Bool+anyStandard a s = V.any (A.isStd a) cs+  where cs = characters s

@@ -37,8 +37,8 @@ -- 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 "translate: reading frame is larger than 2."-                         | rf < 0    = error "translate: reading frame is negative."+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

@@ -21,7 +21,6 @@   , sequencesToFasta   ) where -import           Control.Lens import qualified Data.ByteString.Lazy.Char8   as L  import           ELynx.Data.Sequence.Sequence@@ -31,7 +30,7 @@  -- | Convert a 'Sequence' to Fasta format. sequenceToFasta :: Sequence -> L.ByteString-sequenceToFasta s = L.unlines [ fastaHeader $ s^.name , fromCharacters $ s^.characters ]+sequenceToFasta s = L.unlines [ fastaHeader $ name s, toByteString $ characters s ]  -- | Convert a list 'Sequence's to Fasta format. A newline is added between any -- two 'Sequence's.

@@ -0,0 +1,77 @@+{- |+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)

@@ -20,7 +20,8 @@   module ELynx.Import.Sequence.Fasta-  ( fastaSequence+  ( Parser+  , fastaSequence   , fasta   ) where @@ -72,7 +73,7 @@                      let !alph  = S.map toWord (A.all . alphabetSpec $ a)                      lns <- some (sequenceLine alph)                      _  <- many eol-                     return $ Sequence hd a (toCharacters $ L.concat lns)+                     return $ Sequence hd a (fromByteString $ L.concat lns)  -- | Parse a Fasta file with given 'Alphabet'. fasta :: Alphabet -> Parser [Sequence]

@@ -37,8 +37,8 @@ probMatrix :: RateMatrix -> Double -> ProbMatrix probMatrix q t | t == 0 = if rows q == cols q                           then ident (rows q)-                          else error "Matrix is not square."-               | t < 0 = error "Time is negative."+                          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@@ -48,7 +48,6 @@ -- 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) 

@@ -1,5 +1,5 @@ {- |-   Description :  Functions to work with transition probability matrices on rooted trees+   Description :  Work with transition probability matrices on rooted trees    Copyright   :  (c) Dominik Schrempf 2017    License     :  GPLv3 @@ -17,27 +17,30 @@  module ELynx.Simulate.MarkovProcessAlongTree   ( -- * Single rate matrix.-    simulateNSitesAlongTree-  , simulateAndFlattenNSitesAlongTree+    simulate+  , simulateAndFlatten     -- * Mixture models.-  , simulateNSitesAlongTreeMixtureModel-  , simulateAndFlattenNSitesAlongTreeMixtureModel+  , 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-import           System.Random.MWC.Distributions+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 -measureableTreeToProbTree :: (Measurable a) => RateMatrix -> Tree a -> Tree ProbMatrix-measureableTreeToProbTree q = fmap (probMatrix q . getLen)+toProbTree :: (Measurable a) => RateMatrix -> Tree a -> Tree ProbMatrix+toProbTree q = fmap (probMatrix q . getLen)  getRootStates :: PrimMonad m   => Int -> StationaryDistribution -> Gen (PrimState m) -> m [State]@@ -46,48 +49,48 @@ -- | 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.-simulateAndFlattenNSitesAlongTree :: (PrimMonad m, Measurable a)+simulateAndFlatten :: (PrimMonad m, Measurable a)   => Int -> StationaryDistribution -> ExchangeabilityMatrix -> Tree a -> Gen (PrimState m) -> m [[State]]-simulateAndFlattenNSitesAlongTree n d e t g = do+simulateAndFlatten n d e t g = do   let q = fromExchangeabilityMatrix e d-      pt = measureableTreeToProbTree q t+      pt = toProbTree q t   is <- getRootStates n d g-  simulateAndFlattenAlongProbTree is pt 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.-simulateAndFlattenAlongProbTree :: (PrimMonad m)+simulateAndFlatten' :: (PrimMonad m)   => [State] -> Tree ProbMatrix -> Gen (PrimState m) -> m [[State]]-simulateAndFlattenAlongProbTree is (Node p f) g = do+simulateAndFlatten' is (Node p f) g = do   is' <- mapM (\i -> jump i p g) is   if null f     then return [is']-    else concat <$> sequence [simulateAndFlattenAlongProbTree is' t g | t <- f]+    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.-simulateNSitesAlongTree :: (PrimMonad m, Measurable a)+simulate :: (PrimMonad m, Measurable a)   => Int -> StationaryDistribution -> ExchangeabilityMatrix -> Tree a -> Gen (PrimState m) -> m (Tree [State])-simulateNSitesAlongTree n d e t g = do+simulate n d e t g = do   let q = fromExchangeabilityMatrix e d-      pt = measureableTreeToProbTree q t+      pt = toProbTree q t   is <- getRootStates n d g-  simulateAlongProbTree is pt 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.-simulateAlongProbTree :: (PrimMonad m)+simulate' :: (PrimMonad m)   => [State] -> Tree ProbMatrix -> Gen (PrimState m) -> m (Tree [State])-simulateAlongProbTree is (Node p f) g = do+simulate' is (Node p f) g = do   is' <- mapM (\i -> jump i p g) is-  f' <- sequence [simulateAlongProbTree is' t g | t <- f]+  f' <- sequence [simulate' is' t g | t <- f]   return $ Node is' f' -measureableTreeToProbTreeMixtureModel :: (Measurable a)+toProbTreeMixtureModel :: (Measurable a)   => [RateMatrix] -> Tree a -> Tree [ProbMatrix]-measureableTreeToProbTreeMixtureModel qs =-  fmap (\a -> [probMatrix q . getLen $ a | q <- qs])+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])@@ -98,41 +101,54 @@  -- | Simulate a number of sites for a given set of substitution models with -- corresponding weights. Forget states at internal nodes. See also--- 'simulateAndFlattenNSitesAlongTree'.-simulateAndFlattenNSitesAlongTreeMixtureModel :: (PrimMonad m, Measurable a)+-- 'simulateAndFlatten'.+simulateAndFlattenMixtureModel :: (PrimMonad m, Measurable a)   => Int -> Vector R -> [StationaryDistribution] -> [ExchangeabilityMatrix] -> Tree a   -> Gen (PrimState m) -> m [[State]]-simulateAndFlattenNSitesAlongTreeMixtureModel n ws ds es t g = do+simulateAndFlattenMixtureModel n ws ds es t g = do   let qs = zipWith fromExchangeabilityMatrix es ds-      pt = measureableTreeToProbTreeMixtureModel qs t+      pt = toProbTreeMixtureModel qs t   (cs, is) <- getComponentsAndRootStates n ws ds g-  simulateAndFlattenAlongProbTreeMixtureModel is cs pt g+  simulateAndFlattenMixtureModel' is cs pt g -simulateAndFlattenAlongProbTreeMixtureModel :: (PrimMonad m)+simulateAndFlattenMixtureModel' :: (PrimMonad m)   => [State] -> [Int] -> Tree [ProbMatrix] -> Gen (PrimState m) -> m [[State]]-simulateAndFlattenAlongProbTreeMixtureModel is cs (Node ps f) g+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 [ simulateAndFlattenAlongProbTreeMixtureModel is' cs t g | t <- f ]+         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--- 'simulateNSitesAlongTree'.-simulateNSitesAlongTreeMixtureModel :: (PrimMonad m, Measurable a)+-- 'simulate'.+simulateMixtureModel :: (PrimMonad m, Measurable a)   => Int -> Vector R -> [StationaryDistribution] -> [ExchangeabilityMatrix] -> Tree a   -> Gen (PrimState m) -> m (Tree [State])-simulateNSitesAlongTreeMixtureModel n ws ds es t g = do+simulateMixtureModel n ws ds es t g = do   let qs = zipWith fromExchangeabilityMatrix es ds-      pt = measureableTreeToProbTreeMixtureModel qs t+      pt = toProbTreeMixtureModel qs t   (cs, is) <- getComponentsAndRootStates n ws ds g-  simulateAlongProbTreeMixtureModel is cs pt 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.-simulateAlongProbTreeMixtureModel :: (PrimMonad m)+simulateMixtureModel' :: (PrimMonad m)   => [State] -> [Int] -> Tree [ProbMatrix] -> Gen (PrimState m) -> m (Tree [State])-simulateAlongProbTreeMixtureModel is cs (Node ps f) g = do+simulateMixtureModel' is cs (Node ps f) g = do   is' <- sequence [ jump i (ps !! c) g | (i, c) <- zip is cs ]-  f'  <- sequence [ simulateAlongProbTreeMixtureModel is' cs t g | t <- f ]+  f'  <- sequence [ simulateMixtureModel' is' cs t g | t <- f ]   return $ Node is' f'

@@ -14,7 +14,6 @@ module ELynx.Data.MarkovProcess.AminoAcidSpec   (spec) where -import           Control.Lens import           Numeric.LinearAlgebra import           Test.Hspec @@ -117,13 +116,13 @@ statDistUniform = vector $ replicate 20 0.05  statDistLG :: R.StationaryDistribution-statDistLG = lg ^. S.stationaryDistribution+statDistLG = S.stationaryDistribution lg  exchLG :: R.ExchangeabilityMatrix-exchLG = lg ^. S.exchangeabilityMatrix+exchLG = S.exchangeabilityMatrix lg  rmLG :: R.RateMatrix-rmLG = S.getRateMatrix lg+rmLG = S.rateMatrix lg  spec :: Spec spec = do@@ -142,14 +141,14 @@    describe "lgCustom" $     it "stationary distribution can be recovered" $ do-    let f = R.getStationaryDistribution $ S.getRateMatrix $ lgCustom Nothing statDistUniform+    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.getRateMatrix poisson) `nearlyEqVec` statDistUniform `shouldBe` True+    R.getStationaryDistribution (S.rateMatrix poisson) `nearlyEqVec` statDistUniform `shouldBe` True    describe "poissonCustom" $     it "stationary distribution can be recovered" $ do-    let f = R.getStationaryDistribution $ S.getRateMatrix $ poissonCustom Nothing statDistLGPython+    let f = R.getStationaryDistribution $ S.rateMatrix $ poissonCustom Nothing statDistLGPython     f `nearlyEqVec` statDistLGPython `shouldBe` True

@@ -32,5 +32,5 @@ spec =   describe "getStationaryDistribution" $   it "extracts the stationary distribution from a rate matrix" $ do-  let sd = getStationaryDistribution (getRateMatrix hkyModel)+  let sd = getStationaryDistribution (rateMatrix hkyModel)   sd `nearlyEqVec` stationaryDist `shouldBe` True

@@ -0,0 +1,40 @@+{- |+Module      :  ELynx.Data.Sequence.AlignmentSpec+Copyright   :  (c) Dominik Schrempf 2018+License     :  GPL-3++Maintainer  :  dominik.schrempf@gmail.com+Stability   :  unstable+Portability :  portable++Creation date: Fri Oct  5 14:25:42 2018.++-}++module ELynx.Data.Sequence.AlignmentSpec+  (spec) where++import qualified Data.ByteString.Lazy.Char8    as L+import qualified Data.Matrix.Unboxed           as M+import           Test.Hspec++import           ELynx.Data.Alphabet.Alphabet+import           ELynx.Data.Alphabet.Character+import           ELynx.Data.Sequence.Alignment+import           ELynx.Import.Sequence.Fasta+import           ELynx.Tools.InputOutput+import           Files++ssData :: M.Matrix Character+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++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

@@ -1,40 +0,0 @@-{- |-Module      :  ELynx.Data.Sequence.MultiSequenceAlignmentSpec-Copyright   :  (c) Dominik Schrempf 2018-License     :  GPL-3--Maintainer  :  dominik.schrempf@gmail.com-Stability   :  unstable-Portability :  portable--Creation date: Fri Oct  5 14:25:42 2018.---}--module ELynx.Data.Sequence.MultiSequenceAlignmentSpec-  (spec) where--import qualified Data.ByteString.Lazy.Char8                 as L-import qualified Data.Matrix.Unboxed                        as M-import           Test.Hspec--import           ELynx.Data.Alphabet.Alphabet-import           ELynx.Data.Alphabet.Character-import           ELynx.Data.Sequence.MultiSequenceAlignment-import           ELynx.Import.Sequence.Fasta-import           ELynx.Tools.InputOutput-import           Files--ssData :: M.Matrix Character-ssData = M.fromLists $ map (reverse . map fromChar) [ "AAA", "GAA", "TAA" ]--ssMSA :: MultiSequenceAlignment-ssMSA = MultiSequenceAlignment (map L.pack ["SEQUENCE_1", "SEQUENCE_2", "SEQUENCE_3"]) DNAI ssData--spec :: Spec-spec =-  describe "subSample" $-  it "correctly sub sample an MSA" $ do-    msa <- either error id . fromSequenceList <$> parseFileWith (fasta DNAI) fastaNucleotideIUPACFN-    let ss = subSample [0,3,5] msa-    ss `shouldBe` ssMSA

@@ -0,0 +1,46 @@+{- |+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

@@ -16,7 +16,7 @@  import           Data.Either import           ELynx.Data.Alphabet.Alphabet-import           ELynx.Data.Sequence.MultiSequenceAlignment+import qualified ELynx.Data.Sequence.Alignment as M import           ELynx.Import.Sequence.Fasta import           ELynx.Tools.InputOutput import           Files@@ -24,27 +24,27 @@  spec :: Spec spec =-  describe "fastaFileMSA" $ do+  describe "fastaFileAlignment" $ do     it "parses a fasta file with nucleotide sequences with equal length" $ do-      msa <- either error id . fromSequenceList <$> parseFileWith (fasta DNA) fastaNucleotideFN-      msaNSequences msa `shouldBe` 3-      msaLength msa `shouldBe` 40+      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-      msa <- either error id . fromSequenceList <$> parseFileWith (fasta DNAI) fastaNucleotideIUPACFN-      msaNSequences msa `shouldBe` 3-      msaLength msa `shouldBe` 40+      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-      emsa <- runParserOnFile (fasta DNAI) fastaErroneousFN-      emsa  `shouldSatisfy` isLeft+      ea <- runParserOnFile (fasta DNAI) fastaErroneousFN+      ea  `shouldSatisfy` isLeft      it "parses a fasta file with amino acid sequences with equal length" $ do-      msa <- either error id . fromSequenceList <$> parseFileWith (fasta Protein) fastaAminoAcidFN-      msaNSequences msa `shouldBe` 2-      msaLength msa `shouldBe` 237+      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-      msa <- runParserOnFile (fasta ProteinI) fastaErroneousFN-      msa  `shouldSatisfy` isLeft+      a <- runParserOnFile (fasta ProteinI) fastaErroneousFN+      a  `shouldSatisfy` isLeft 

@@ -14,7 +14,6 @@ module ELynx.Simulate.MarkovProcessAlongTreeSpec   (spec) where -import           Control.Lens import           Data.Tree import           System.Random.MWC import           Test.Hspec@@ -28,14 +27,14 @@ -- 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 PhyloIntLabel-testTree = Node (PhyloLabel 0 Nothing 0.0) [ Node (PhyloLabel 1 Nothing 1.0) [], Node (PhyloLabel 2 Nothing 1.0) [] ]+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 PhyloIntLabel-nullTree = Node (PhyloLabel 0 Nothing 0.0) []+nullTree :: Tree (PhyloLabel Int)+nullTree = Node (PhyloLabel 0 Nothing (Just 0.0)) []  nullStateTree :: Tree State nullStateTree = fmap (const 0) nullTree@@ -48,13 +47,13 @@   describe "simulateNSitesAlongTree" $ do     it "simulates one site along an easy tree" $ do       gen <- create-      tr <- simulateNSitesAlongTree 1 d e nullTree gen+      tr <- simulate 1 d e nullTree gen       fmap head tr `shouldBe` nullStateTree      it "simulates some sites along a harder tree" $ do       gen <- create-      tr <- simulateNSitesAlongTree 10 d e testTree gen+      tr <- simulate 10 d e testTree gen       (length . rootLabel $ tr) `shouldBe` 10-        where d = jc ^. S.stationaryDistribution-              e = jc ^. S.exchangeabilityMatrix+        where d = S.stationaryDistribution jc+              e = S.exchangeabilityMatrix jc