diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -4,6 +4,10 @@ ## Unreleased changes + +## Version 0.5.0.2 + +- Speed up mixture model simulation. - Improve rooting functions. - Improve `Topology` data type (but still a lot to do). - Various additions to the documentation. diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -2,7 +2,7 @@ # The ELynx Suite -Version: 0.5.0.1. +Version: 0.5.0.2. Reproducible evolution made easy.

@@ -73,9 +73,9 @@ slynx --help | head -n -16 - ELynx Suite version 0.5.0.1. + ELynx Suite version 0.5.0.2. Developed by Dominik Schrempf. - Compiled on December 18, 2020, at 10:14 am, UTC. + Compiled on February 19, 2021, at 13:42 pm, UTC. Usage: slynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME] [-f|--force] [--no-elynx-file] COMMAND @@ -90,7 +90,7 @@ Specify base name of output file -f,--force Ignore previous analysis and overwrite existing output files. - --no-elynx-file Do not write files for needed for reproducibility. + --no-elynx-file Do not write data required to reproduce an analysis. Available commands: concatenate Concatenate sequences found in input files. @@ -120,9 +120,9 @@ slynx concatenate --help - ELynx Suite version 0.5.0.1. + ELynx Suite version 0.5.0.2. Developed by Dominik Schrempf. - Compiled on December 18, 2020, at 10:14 am, UTC. + Compiled on February 19, 2021, at 13:42 pm, UTC. Usage: slynx concatenate (-a|--alphabet NAME) INPUT-FILE Concatenate sequences found in input files. @@ -141,9 +141,9 @@ slynx examine --help - ELynx Suite version 0.5.0.1. + ELynx Suite version 0.5.0.2. Developed by Dominik Schrempf. - Compiled on December 18, 2020, at 10:14 am, UTC. + Compiled on February 19, 2021, at 13:42 pm, UTC. Usage: slynx examine (-a|--alphabet NAME) INPUT-FILE [--per-site] Examine sequences. If data is a multi sequence alignment, additionally analyze columns. @@ -163,9 +163,9 @@ slynx filter-rows --help - ELynx Suite version 0.5.0.1. + ELynx Suite version 0.5.0.2. Developed by Dominik Schrempf. - Compiled on December 18, 2020, at 10:14 am, UTC. + Compiled on February 19, 2021, at 13:42 pm, UTC. Usage: slynx filter-rows (-a|--alphabet NAME) INPUT-FILE [--longer-than LENGTH] [--shorter-than LENGTH] [--standard-characters] @@ -186,9 +186,9 @@ slynx filter-columns --help - ELynx Suite version 0.5.0.1. + ELynx Suite version 0.5.0.2. Developed by Dominik Schrempf. - Compiled on December 18, 2020, at 10:14 am, UTC. + Compiled on February 19, 2021, at 13:42 pm, UTC. Usage: slynx filter-columns (-a|--alphabet NAME) INPUT-FILE [--standard-chars DOUBLE] @@ -210,9 +210,9 @@ slynx simulate --help - ELynx Suite version 0.5.0.1. + ELynx Suite version 0.5.0.2. Developed by Dominik Schrempf. - Compiled on December 18, 2020, at 10:14 am, UTC. + Compiled on February 19, 2021, at 13:42 pm, UTC. Usage: slynx simulate (-t|--tree-file Name) [-s|--substitution-model MODEL] [-m|--mixture-model MODEL] [-e|--edm-file NAME] @@ -289,9 +289,9 @@ slynx sub-sample --help - ELynx Suite version 0.5.0.1. + ELynx Suite version 0.5.0.2. Developed by Dominik Schrempf. - Compiled on December 18, 2020, at 10:14 am, UTC. + Compiled on February 19, 2021, at 13:42 pm, UTC. Usage: slynx sub-sample (-a|--alphabet NAME) INPUT-FILE (-n|--number-of-sites INT) @@ -319,9 +319,9 @@ slynx translate --help - ELynx Suite version 0.5.0.1. + ELynx Suite version 0.5.0.2. Developed by Dominik Schrempf. - Compiled on December 18, 2020, at 10:14 am, UTC. + Compiled on February 19, 2021, at 13:42 pm, UTC. Usage: slynx translate (-a|--alphabet NAME) INPUT-FILE (-r|--reading-frame INT) (-u|--universal-code CODE) @@ -344,9 +344,9 @@ tlynx --help | head -n -16 - ELynx Suite version 0.5.0.1. + ELynx Suite version 0.5.0.2. Developed by Dominik Schrempf. - Compiled on December 18, 2020, at 10:14 am, UTC. + Compiled on February 19, 2021, at 13:42 pm, UTC. Usage: tlynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME] [-f|--force] [--no-elynx-file] COMMAND @@ -361,7 +361,7 @@ Specify base name of output file -f,--force Ignore previous analysis and overwrite existing output files. - --no-elynx-file Do not write files for needed for reproducibility. + --no-elynx-file Do not write data required to reproduce an analysis. Available commands: compare Compare two phylogenetic trees (compute distances and branch-wise differences). @@ -384,9 +384,9 @@ tlynx compare --help - ELynx Suite version 0.5.0.1. + ELynx Suite version 0.5.0.2. Developed by Dominik Schrempf. - Compiled on December 18, 2020, at 10:14 am, UTC. + Compiled on February 19, 2021, at 13:42 pm, UTC. Usage: tlynx compare [-n|--normalize] [-b|--bipartitions] [-t|--intersect] [-f|--newick-format FORMAT] NAMES @@ -413,9 +413,9 @@ tlynx examine --help - ELynx Suite version 0.5.0.1. + ELynx Suite version 0.5.0.2. Developed by Dominik Schrempf. - Compiled on December 18, 2020, at 10:14 am, UTC. + Compiled on February 19, 2021, at 13:42 pm, UTC. Usage: tlynx examine INPUT-FILE [-f|--newick-format FORMAT] Compute summary statistics of phylogenetic trees. @@ -437,9 +437,9 @@ tlynx simulate --help - ELynx Suite version 0.5.0.1. + ELynx Suite version 0.5.0.2. Developed by Dominik Schrempf. - Compiled on December 18, 2020, at 10:14 am, UTC. + Compiled on February 19, 2021, at 13:42 pm, UTC. Usage: tlynx simulate (-t|--nTrees INT) (-n|--nLeaves INT) PROCESS [-u|--sub-sample DOUBLE] [-s|--summary-statistics] @@ -474,9 +474,9 @@ elynx --help | head -n -16 - ELynx Suite version 0.5.0.1. + ELynx Suite version 0.5.0.2. Developed by Dominik Schrempf. - Compiled on December 18, 2020, at 10:14 am, UTC. + Compiled on February 19, 2021, at 13:42 pm, UTC. Usage: elynx COMMAND Validate and redo past ELynx analyses diff --git a/slynx.cabal b/slynx.cabal --- a/slynx.cabal +++ b/slynx.cabal @@ -1,6 +1,6 @@ cabal-version: 2.2 name: slynx -version: 0.5.0.1 +version: 0.5.0.2 license: GPL-3.0-or-later license-file: LICENSE copyright: Dominik Schrempf (2020) @@ -47,23 +47,22 @@ default-language: Haskell2010 ghc-options: -Wall -Wunused-packages build-depends: - async >=2.2.2 && <2.3, - attoparsec >=0.13.2.4 && <0.14, + attoparsec >=0.13.2.5, base >=4.7 && <5, - bytestring >=0.10.10.0 && <0.11, - containers >=0.6.2.1 && <0.7, - elynx-markov >=0.5.0.1 && <0.6, - elynx-seq >=0.5.0.1 && <0.6, - elynx-tools >=0.5.0.1 && <0.6, - elynx-tree >=0.5.0.1 && <0.6, - hmatrix >=0.20.0.0 && <0.21, - monad-logger >=0.3.35 && <0.4, - mwc-random >=0.14.0.0 && <0.15, - optparse-applicative >=0.15.1.0 && <0.16, - statistics >=0.15.2.0 && <0.16, - text >=1.2.3.2 && <1.3, - transformers >=0.5.6.2 && <0.6, - vector >=0.12.1.2 && <0.13 + bytestring >=0.10.12.0, + containers >=0.6.2.1, + elynx-markov >=0.5.0.2, + elynx-seq >=0.5.0.2, + elynx-tools >=0.5.0.2, + elynx-tree >=0.5.0.2, + hmatrix >=0.20.1, + monad-logger >=0.3.36, + mwc-random >=0.14.0.0, + optparse-applicative >=0.15.1.0, + statistics >=0.15.2.0, + text >=1.2.4.1, + transformers >=0.5.6.2, + vector >=0.12.2.0 executable slynx main-is: Main.hs diff --git a/src/SLynx/Simulate/PhyloModel.hs b/src/SLynx/Simulate/PhyloModel.hs --- a/src/SLynx/Simulate/PhyloModel.hs +++ b/src/SLynx/Simulate/PhyloModel.hs @@ -21,8 +21,8 @@ import Data.Attoparsec.ByteString.Char8 import qualified Data.ByteString.Char8 as BS import Data.Either (rights) -import Data.List.NonEmpty (fromList) import Data.Maybe +import qualified Data.Vector as V import ELynx.Data.MarkovProcess.AminoAcid import ELynx.Data.MarkovProcess.CXXModels import qualified ELynx.Data.MarkovProcess.MixtureModel as M @@ -170,7 +170,7 @@ then fail $ head errs else return $ - M.fromSubstitutionModels edmName (fromList ws) (fromList $ rights sms) + M.fromSubstitutionModels edmName (V.fromList ws) (V.fromList $ rights sms) cxxModel :: Maybe [M.Weight] -> Parser M.MixtureModel cxxModel mws = do @@ -185,7 +185,7 @@ sms <- parseSubstitutionModel `sepBy1` char separator _ <- char mmEnd -- XXX: The use of `Data.List.NonEmpty.fromList` leads to uninformative error messages. - return $ M.fromSubstitutionModels "MIXTURE" (fromList ws) (fromList sms) + return $ M.fromSubstitutionModels "MIXTURE" (V.fromList ws) (V.fromList sms) mixtureModel :: Maybe [EDMComponent] -> Maybe [M.Weight] -> Parser M.MixtureModel diff --git a/src/SLynx/Simulate/Simulate.hs b/src/SLynx/Simulate/Simulate.hs --- a/src/SLynx/Simulate/Simulate.hs +++ b/src/SLynx/Simulate/Simulate.hs @@ -18,8 +18,6 @@ where import Control.Applicative ((<|>)) -import Control.Concurrent -import Control.Concurrent.Async import Control.Monad import Control.Monad.IO.Class import Control.Monad.Logger @@ -27,19 +25,22 @@ import Control.Monad.Trans.Reader (ask) import qualified Data.ByteString.Builder as BB import qualified Data.ByteString.Lazy.Char8 as BL -import Data.List.NonEmpty (toList) +import Data.List import Data.Maybe import qualified Data.Set as Set import qualified Data.Text as T import qualified Data.Text.Lazy as LT import qualified Data.Text.Lazy.Encoding as LT -import qualified Data.Vector.Unboxed as V +import qualified Data.Vector as V +import qualified Data.Vector.Storable as VS +import qualified Data.Vector.Unboxed as U import ELynx.Data.Alphabet.Alphabet as A +import qualified ELynx.Data.MarkovProcess.AminoAcid as MA import ELynx.Data.MarkovProcess.GammaRateHeterogeneity -import qualified ELynx.Data.MarkovProcess.MixtureModel as M -import qualified ELynx.Data.MarkovProcess.PhyloModel as P -import qualified ELynx.Data.MarkovProcess.SubstitutionModel as SM -import qualified ELynx.Data.Sequence.Alignment as A +import qualified ELynx.Data.MarkovProcess.MixtureModel as MM +import qualified ELynx.Data.MarkovProcess.PhyloModel as MP +import qualified ELynx.Data.MarkovProcess.RateMatrix as MR +import qualified ELynx.Data.MarkovProcess.SubstitutionModel as MS import qualified ELynx.Data.Sequence.Sequence as Seq hiding ( name, ) @@ -49,60 +50,68 @@ import ELynx.Simulate.MarkovProcessAlongTree import ELynx.Tools import ELynx.Tree -import Numeric.LinearAlgebra hiding - ( toList, - (<>), - ) import SLynx.Simulate.Options import SLynx.Simulate.PhyloModel import System.Random.MWC import Text.Printf +getDistLine :: Int -> MR.StationaryDistribution -> BB.Builder +getDistLine i d = + BB.intDec i + <> BB.char8 ' ' + <> s + where + s = mconcat $ intersperse (BB.char8 ' ') $ map BB.doubleDec $ VS.toList d + +writeSiteDists :: [Int] -> V.Vector MR.StationaryDistribution -> ELynx SimulateArguments () +-- writeSiteDists is ds = out "site distributions of distribution mixture model" output ".sitedists" +writeSiteDists componentIs ds = do + mbn <- outFileBaseName . global <$> ask + case mbn of + Nothing -> return () + Just bn -> liftIO $ BL.writeFile (bn <> ".sitedists") output + where + dsPaml = V.map MA.alphaToPamlVec ds + lns = [getDistLine i d | (i, c) <- zip [1 ..] componentIs, let d = dsPaml V.! c] + output = BB.toLazyByteString $ mconcat $ intersperse (BB.char8 '\n') lns + -- Simulate a 'Alignment' for a given phylogenetic model, -- phylogenetic tree, and alignment length. simulateAlignment :: (HasLength e, HasName a) => - P.PhyloModel -> + MP.PhyloModel -> Tree e a -> Int -> GenIO -> - IO A.Alignment + ELynx SimulateArguments () simulateAlignment pm t' n g = do let t = fromLength . getLen <$> toTreeBranchLabels t' - c <- getNumCapabilities - gs <- splitGen c g - let chunks = getChunks c n - leafStatesS <- case pm of - -- TODO @performace: This parallelization is not very intelligent, because - -- the matrix exponentiation is done in all threads. So ten threads will - -- exponentiate the same matrix ten times. - P.SubstitutionModel sm -> - mapConcurrently - (\(num, gen) -> simulateAndFlatten num d e t gen) - (zip chunks gs) + leafStates <- case pm of + MP.SubstitutionModel sm -> liftIO $ simulateAndFlattenPar n d e t g where - d = SM.stationaryDistribution sm - e = SM.exchangeabilityMatrix sm - -- P.MixtureModel mm -> mapConcurrently - -- (\(num, gen) -> simulateAndFlattenNSitesAlongTreeMixtureModel num ws ds es t gen) (zip chunks gs) - P.MixtureModel mm -> simulateAndFlattenMixtureModelPar n ws ds es t g + d = MS.stationaryDistribution sm + e = MS.exchangeabilityMatrix sm + MP.MixtureModel mm -> do + (cs, ss) <- liftIO $ simulateAndFlattenMixtureModelPar n ws ds es t g + -- TODO: Writing site distributions only makes sense for EDM models. + -- Remove this if not needed or improve to be helpful in general. + writeSiteDists cs ds + return ss where - ws = vector . toList $ M.getWeights mm - ds = map SM.stationaryDistribution $ toList $ M.getSubstitutionModels mm - es = map SM.exchangeabilityMatrix $ toList $ M.getSubstitutionModels mm - -- XXX @performace. The horizontal concatenation might be slow. If so, - -- 'concatenateSeqs' or 'concatenateAlignments' can be used, which directly - -- appends vectors. - let leafStates = horizontalConcat leafStatesS - leafNames = map getName $ leaves t' - code = P.getAlphabet pm + ws = MM.getWeights mm + ds = V.map MS.stationaryDistribution $ MM.getSubstitutionModels mm + es = V.map MS.exchangeabilityMatrix $ MM.getSubstitutionModels mm + let leafNames = map getName $ leaves t' + code = MP.getAlphabet pm -- XXX: Probably use type safe stuff here? alph = A.all $ alphabetSpec code sequences = - [ Seq.Sequence (fromName sName) "" code (V.fromList $ map (`Set.elemAt` alph) ss) + [ Seq.Sequence (fromName sName) "" code (U.fromList $ map (`Set.elemAt` alph) ss) | (sName, ss) <- zip leafNames leafStates ] - return $ either error id $ A.fromSequences sequences + output = sequencesToFasta sequences + $(logInfo) "" + out "simulated multi sequence alignment" output ".fasta" -- Summarize EDM components; line to be printed to screen or log. summarizeEDMComponents :: [EDMComponent] -> BL.ByteString @@ -112,17 +121,20 @@ ++ show (length cs) ++ " components." -reportModel :: P.PhyloModel -> ELynx SimulateArguments () +reportModel :: MP.PhyloModel -> ELynx SimulateArguments () reportModel m = do as <- global <$> ask if writeElynxFile as - then (do let bn = outFileBaseName as - case bn of - Nothing -> - $(logInfo) - "No output file provided; omit writing machine-readable phylogenetic model." - Just _ -> - out "model definition (machine readable)" (BL.pack (show m) <> "\n") ".model.gz") + then + ( do + let bn = outFileBaseName as + case bn of + Nothing -> + $(logInfo) + "No output file provided; omit writing machine-readable phylogenetic model." + Just _ -> + out "model definition (machine readable)" (BL.pack (show m) <> "\n") ".model.gz" + ) else $(logInfo) "No elynx file required; omit writing machine-readable phylogenetic model." pretty :: Length -> String @@ -150,60 +162,60 @@ b = totalBranchLength t -- Summarize a substitution model; lines to be printed to screen or log. -summarizeSM :: SM.SubstitutionModel -> [BL.ByteString] +summarizeSM :: MS.SubstitutionModel -> [BL.ByteString] summarizeSM sm = map BL.pack $ - (show (SM.alphabet sm) ++ " substitution model: " ++ SM.name sm ++ ".") : - ["Parameters: " ++ show (SM.params sm) ++ "." | not (null (SM.params sm))] - ++ case SM.alphabet sm of + (show (MS.alphabet sm) ++ " substitution model: " ++ MS.name sm ++ ".") : + ["Parameters: " ++ show (MS.params sm) ++ "." | not (null (MS.params sm))] + ++ case MS.alphabet sm of DNA -> [ "Stationary distribution: " - ++ dispv precision (SM.stationaryDistribution sm) + ++ dispv precision (MS.stationaryDistribution sm) ++ ".", "Exchangeability matrix:\n" - ++ dispmi 2 precision (SM.exchangeabilityMatrix sm) + ++ dispmi 2 precision (MS.exchangeabilityMatrix sm) ++ ".", - "Scale: " ++ show (roundN precision $ SM.totalRate sm) ++ "." + "Scale: " ++ show (roundN precision $ MS.totalRate sm) ++ "." ] Protein -> [ "Stationary distribution: " - ++ dispv precision (SM.stationaryDistribution sm) + ++ dispv precision (MS.stationaryDistribution sm) ++ ".", - "Scale: " ++ show (roundN precision $ SM.totalRate sm) ++ "." + "Scale: " ++ show (roundN precision $ MS.totalRate sm) ++ "." ] _ -> error "Extended character sets are not supported with substitution models." -- Summarize a mixture model component; lines to be printed to screen or log. -summarizeMMComponent :: M.Component -> [BL.ByteString] +summarizeMMComponent :: MM.Component -> [BL.ByteString] summarizeMMComponent c = BL.pack "Weight: " - <> (BB.toLazyByteString . BB.doubleDec $ M.weight c) : - summarizeSM (M.substModel c) + <> (BB.toLazyByteString . BB.doubleDec $ MM.weight c) : + summarizeSM (MM.substModel c) -- Summarize a mixture model; lines to be printed to screen or log. -summarizeMM :: M.MixtureModel -> [BL.ByteString] +summarizeMM :: MM.MixtureModel -> [BL.ByteString] summarizeMM m = - [ BL.pack $ "Mixture model: " ++ M.name m ++ ".", + [ BL.pack $ "Mixture model: " ++ MM.name m ++ ".", BL.pack $ "Number of components: " ++ show n ++ "." ] ++ detail where - n = length $ M.components m + n = length $ MM.components m detail = if n <= 100 then concat [ BL.pack ("Component " ++ show i ++ ":") : summarizeMMComponent c - | (i, c) <- zip [1 :: Int ..] (toList $ M.components m) + | (i, c) <- zip [1 :: Int ..] (V.toList $ MM.components m) ] else [] -- Summarize a phylogenetic model; lines to be printed to screen or log. -summarizePM :: P.PhyloModel -> [BL.ByteString] -summarizePM (P.MixtureModel mm) = summarizeMM mm -summarizePM (P.SubstitutionModel sm) = summarizeSM sm +summarizePM :: MP.PhyloModel -> [BL.ByteString] +summarizePM (MP.MixtureModel mm) = summarizeMM mm +summarizePM (MP.SubstitutionModel sm) = summarizeSM sm -- | Simulate sequences. simulateCmd :: ELynx SimulateArguments () @@ -286,7 +298,4 @@ Random -> error "simulateCmd: seed not available; please contact maintainer." Fixed s -> liftIO $ initialize s - alignment <- liftIO $ simulateAlignment phyloModel t' alignmentLength gen - let output = (sequencesToFasta . A.toSequences) alignment - $(logInfo) "" - out "simulated multi sequence alignment" output ".fasta" + simulateAlignment phyloModel t' alignmentLength gen