packages feed

tlynx 0.3.3 → 0.3.4

raw patch · 13 files changed

+294/−482 lines, 13 files

Files

ChangeLog.md view
@@ -5,6 +5,14 @@ ## Unreleased changes  +## Version 0.3.4++-   Improve `slynx examine`; show hamming distance; show constant sites.+-   PhyloStrict -> PhyloExplicit; and some conversion functions were changed.+-   `tlynx coalesce` was merged into `tlynx simulate`, the syntax has changed; see+    `tlynx simulate --help`.++ ## Version 0.3.3  -   Fix test suites.
README.md view
@@ -2,7 +2,7 @@  # The ELynx Suite -Version: 0.3.3.+Version: 0.3.4. Reproducible evolution made easy.  <p align="center"><img src="https://travis-ci.org/dschrempf/elynx.svg?branch=master"/></p>@@ -90,9 +90,9 @@      slynx --help | head -n -16 -    ELynx Suite version 0.3.3.+    ELynx Suite version 0.3.4.     Developed by Dominik Schrempf.-    Compiled on August 18, 2020, at 08:07 am, UTC.+    Compiled on August 21, 2020, at 09:40 am, UTC.          Usage: slynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME]                   [-f|--force] COMMAND@@ -136,9 +136,9 @@      slynx concatenate --help -    ELynx Suite version 0.3.3.+    ELynx Suite version 0.3.4.     Developed by Dominik Schrempf.-    Compiled on August 18, 2020, at 08:07 am, UTC.+    Compiled on August 21, 2020, at 09:40 am, UTC.          Usage: slynx concatenate (-a|--alphabet NAME) INPUT-FILE       Concatenate sequences found in input files.@@ -157,9 +157,9 @@      slynx examine --help -    ELynx Suite version 0.3.3.+    ELynx Suite version 0.3.4.     Developed by Dominik Schrempf.-    Compiled on August 18, 2020, at 08:07 am, UTC.+    Compiled on August 21, 2020, at 09:40 am, UTC.          Usage: slynx examine (-a|--alphabet NAME) INPUT-FILE [--per-site]       Examine sequences. If data is a multi sequence alignment, additionally analyze columns.@@ -179,9 +179,9 @@      slynx filter-rows --help -    ELynx Suite version 0.3.3.+    ELynx Suite version 0.3.4.     Developed by Dominik Schrempf.-    Compiled on August 18, 2020, at 08:07 am, UTC.+    Compiled on August 21, 2020, at 09:40 am, UTC.          Usage: slynx filter-rows (-a|--alphabet NAME) INPUT-FILE [--longer-than LENGTH]                               [--shorter-than LENGTH] [--standard-characters]@@ -202,9 +202,9 @@      slynx filter-columns --help -    ELynx Suite version 0.3.3.+    ELynx Suite version 0.3.4.     Developed by Dominik Schrempf.-    Compiled on August 18, 2020, at 08:07 am, UTC.+    Compiled on August 21, 2020, at 09:40 am, UTC.          Usage: slynx filter-columns (-a|--alphabet NAME) INPUT-FILE                                  [--standard-chars DOUBLE]@@ -226,9 +226,9 @@      slynx simulate --help -    ELynx Suite version 0.3.3.+    ELynx Suite version 0.3.4.     Developed by Dominik Schrempf.-    Compiled on August 18, 2020, at 08:07 am, UTC.+    Compiled on August 21, 2020, at 09:40 am, UTC.          Usage: slynx simulate (-t|--tree-file Name) [-s|--substitution-model MODEL]                            [-m|--mixture-model MODEL] [-e|--edm-file NAME] @@ -305,9 +305,9 @@      slynx sub-sample --help -    ELynx Suite version 0.3.3.+    ELynx Suite version 0.3.4.     Developed by Dominik Schrempf.-    Compiled on August 18, 2020, at 08:07 am, UTC.+    Compiled on August 21, 2020, at 09:40 am, UTC.          Usage: slynx sub-sample (-a|--alphabet NAME) INPUT-FILE                             (-n|--number-of-sites INT)@@ -335,9 +335,9 @@      slynx translate --help -    ELynx Suite version 0.3.3.+    ELynx Suite version 0.3.4.     Developed by Dominik Schrempf.-    Compiled on August 18, 2020, at 08:07 am, UTC.+    Compiled on August 21, 2020, at 09:40 am, UTC.          Usage: slynx translate (-a|--alphabet NAME) INPUT-FILE (-r|--reading-frame INT)                            (-u|--universal-code CODE)@@ -360,9 +360,9 @@      tlynx --help | head -n -16 -    ELynx Suite version 0.3.3.+    ELynx Suite version 0.3.4.     Developed by Dominik Schrempf.-    Compiled on August 18, 2020, at 08:07 am, UTC.+    Compiled on August 21, 2020, at 09:40 am, UTC.          Usage: tlynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME]                   [-f|--force] COMMAND@@ -379,13 +379,12 @@                                output files.          Available commands:-      coalesce                 Simulate phylogenetic trees using the coalescent processes (see also the 'simulate' command for simulations using the birth and death process).       compare                  Compare two phylogenetic trees (compute distances and branch-wise differences).       connect                  Connect two phylogenetic trees in all ways (possibly honoring constraints).       distance                 Compute distances between many phylogenetic trees.       examine                  Compute summary statistics of phylogenetic trees.       shuffle                  Shuffle a phylogenetic tree (keep coalescent times, but shuffle topology and leaves).-      simulate                 Simulate phylogenetic trees using birth and death processes (see also the 'coalesce' command for simulations using the coalescent process).+      simulate                 Simulate phylogenetic trees using a birth and death or coalescent process.               Available tree file formats:@@ -399,9 +398,9 @@      tlynx compare --help -    ELynx Suite version 0.3.3.+    ELynx Suite version 0.3.4.     Developed by Dominik Schrempf.-    Compiled on August 18, 2020, at 08:07 am, UTC.+    Compiled on August 21, 2020, at 09:40 am, UTC.          Usage: tlynx compare [-n|--normalize] [-b|--bipartitions] [-t|--intersect]                           [-f|--newick-format FORMAT] NAMES@@ -428,9 +427,9 @@      tlynx examine --help -    ELynx Suite version 0.3.3.+    ELynx Suite version 0.3.4.     Developed by Dominik Schrempf.-    Compiled on August 18, 2020, at 08:07 am, UTC.+    Compiled on August 21, 2020, at 09:40 am, UTC.          Usage: tlynx examine INPUT-FILE [-f|--newick-format FORMAT]       Compute summary statistics of phylogenetic trees.@@ -452,37 +451,35 @@      tlynx simulate --help -    ELynx Suite version 0.3.3.+    ELynx Suite version 0.3.4.     Developed by Dominik Schrempf.-    Compiled on August 18, 2020, at 08:07 am, UTC.+    Compiled on August 21, 2020, at 09:40 am, UTC.     -    Usage: tlynx simulate [-t|--nTrees INT] [-n|--nLeaves INT] [-H|--height DOUBLE] -                          [-M|--condition-on-mrca] [-l|--lambda DOUBLE] -                          [-m|--mu DOUBLE] [-r|--rho DOUBLE] [-u|--sub-sample] -                          [-s|--summary-statistics] [-S|--seed [INT]]-      Simulate phylogenetic trees using birth and death processes (see also the 'coalesce' command for simulations using the coalescent process).+    Usage: tlynx simulate (-t|--nTrees INT) (-n|--nLeaves INT) PROCESS +                          [-u|--sub-sample DOUBLE] [-s|--summary-statistics] +                          [-S|--seed [INT]]+      Simulate phylogenetic trees using a birth and death or coalescent process.          Available options:       -h,--help                Show this help text       -V,--version             Show version-      -t,--nTrees INT          Number of trees (default: 10)-      -n,--nLeaves INT         Number of leaves per tree (default: 5)-      -H,--height DOUBLE       Fix tree height (no default)-      -M,--condition-on-mrca   Do not condition on height of origin but on height of-                               MRCA-      -l,--lambda DOUBLE       Birth rate lambda (default: 1.0)-      -m,--mu DOUBLE           Death rate mu (default: 0.9)-      -r,--rho DOUBLE          Sampling probability rho (default: 1.0)-      -u,--sub-sample          Perform sub-sampling; see below.-      -s,--summary-statistics  Only output number of children for each branch+      -t,--nTrees INT          Number of trees+      -n,--nLeaves INT         Number of leaves per tree+      -u,--sub-sample DOUBLE   Perform sub-sampling; see below.+      -s,--summary-statistics  For each branch, print length and number of children       -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     -    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+    Available processes:+      birthdeath               Birth and death process+      coalescent               Coalescent process+    +    See, for example, 'tlynx simulate birthdeath --help'.+    Sub-sample with probability p:+      1. Simulate one big tree with n'=round(n/p), n'>=n, leaves;+      2. Randomly sample sub-trees with n leaves.+      (With p=1.0, the same tree is reported over and over again.)   # ELynx@@ -491,9 +488,9 @@      elynx --help | head -n -16 -    ELynx Suite version 0.3.3.+    ELynx Suite version 0.3.4.     Developed by Dominik Schrempf.-    Compiled on August 18, 2020, at 08:07 am, UTC.+    Compiled on August 21, 2020, at 09:40 am, UTC.          Usage: elynx COMMAND       Validate and redo past ELynx analyses
− src/TLynx/Coalesce/Coalesce.hs
@@ -1,145 +0,0 @@-{-# LANGUAGE OverloadedStrings #-}-{-# LANGUAGE TemplateHaskell #-}-{-# LANGUAGE TupleSections #-}---- TODO: MERGE WITH SIMULATE.---- |---   Description :  Simulate reconstructed trees using the coalescent process---   Copyright   :  (c) Dominik Schrempf 2020---   License     :  GPL-3.0-or-later------   Maintainer  :  dominik.schrempf@gmail.com---   Stability   :  unstable---   Portability :  portable------ Creation date: Tue Feb 3 17:00:00 2020.-module TLynx.Coalesce.Coalesce-  ( coalesce,-  )-where--import Control.Concurrent (getNumCapabilities)-import Control.Concurrent.Async.Lifted.Safe-  ( mapConcurrently,-  )-import Control.Monad-import Control.Monad.IO.Class-import Control.Monad.Logger-import Control.Monad.Trans.Reader (ask)-import Control.Parallel.Strategies-import qualified Data.ByteString.Builder as BB-import qualified Data.ByteString.Lazy.Char8 as BL-import Data.Maybe-import qualified Data.Sequence as Seq-import qualified Data.Text as T-import qualified Data.Text.Lazy as LT-import qualified Data.Text.Lazy.Encoding as LT-import ELynx.Data.Tree-import ELynx.Export.Tree.Newick (toNewick)-import ELynx.Simulate.Coalescent (simulate)-import ELynx.Tools-import System.Random.MWC (initialize)-import TLynx.Coalesce.Options-import TLynx.Simulate.Simulate (nSubSamples)---- | Simulate phylogenetic trees.-coalesce :: ELynx CoalesceArguments ()-coalesce = do-  l <- local <$> ask-  let s = argsSumStat l-  logNewSection "Arguments"-  $(logInfo) $ T.pack $ reportCoalesceArguments l-  logNewSection "Simulation"-  c <- liftIO getNumCapabilities-  $(logInfo) $ T.pack $ "Number of used cores: " <> show c-  trs <- case argsRho l of-    Nothing -> simulateNTreesConcurrently-    Just _ -> simulateAndSubSampleNTreesConcurrently-  let ls =-        if s-          then parMap rpar (formatNChildSumStat . toNChildSumStat) trs-          else parMap rpar toNewick (map lengthToPhyloTree trs)-  let res = BL.unlines ls-  out "simulated trees" res ".tree"--simulateNTreesConcurrently :: ELynx CoalesceArguments (Forest Length Int)-simulateNTreesConcurrently = do-  (CoalesceArguments nT nL _ _ (Fixed s)) <- local <$> ask-  c <- liftIO getNumCapabilities-  gs <- liftIO $ initialize s >>= \g -> splitGen c g-  let chunks = getChunks c nT-  trss <--    liftIO $-      mapConcurrently-        (\(n, g) -> replicateM n $ simulate nL g)-        (zip chunks gs)-  return $ concat trss--simulateAndSubSampleNTreesConcurrently ::-  ELynx CoalesceArguments (Forest Length Int)-simulateAndSubSampleNTreesConcurrently = do-  (CoalesceArguments nT nL mR _ (Fixed s)) <- local <$> ask-  c <- liftIO getNumCapabilities-  let r =-        fromMaybe-          ( error-              "cimulateAndSubSampleNTreesConcurrently: no sampling probability given."-          )-          mR-  let nLeavesBigTree = (round $ fromIntegral nL / r) :: Int-  gs <- liftIO $ initialize s >>= \g -> splitGen c g-  let chunks = getChunks c nT-  tr <- liftIO $ simulate nLeavesBigTree (head gs)-  logNewSection $-    T.pack $-      "Simulate one big tree with "-        <> show nLeavesBigTree-        <> " leaves."-  -- Log the base tree.-  $(logInfo) $ LT.toStrict $ LT.decodeUtf8 $ toNewick $ lengthToPhyloTree tr-  logNewSection $-    T.pack $-      "Sub sample "-        <> show nT-        <> " trees with "-        <> show nL-        <> " leaves."-  let lvs = Seq.fromList $ leaves tr-  trss <--    liftIO $-      mapConcurrently-        (\(nSamples, g) -> nSubSamples nSamples lvs nL tr g)-        (zip chunks gs)-  let trs = catMaybes $ concat trss-  return $ map prune trs---- | Pair of branch length with number of extant children.-type BrLnNChildren = (BranchLength, Int)---- | Possible summary statistic of phylogenetic trees. A list of tuples--- (BranchLength, NumberOfExtantChildrenBelowThisBranch).-type NChildSumStat = [BrLnNChildren]---- | Format the summary statistics in the following form:--- @---    nLeaves1 branchLength1---    nLeaves2 branchLength2---    ....-formatNChildSumStat :: NChildSumStat -> BL.ByteString-formatNChildSumStat s =-  BB.toLazyByteString . mconcat $ map formatNChildSumStatLine s--formatNChildSumStatLine :: BrLnNChildren -> BB.Builder-formatNChildSumStatLine (l, n) =-  BB.intDec n <> BB.char8 ' ' <> BB.doubleDec l <> BB.char8 '\n'---- | Compute NChilSumStat for a phylogenetic tree.-toNChildSumStat :: Measurable e => Tree e a -> NChildSumStat-toNChildSumStat (Node br _ []) = [(getLen br, 1)]-toNChildSumStat (Node br _ ts) = (getLen br, sumCh) : concat nChSS-  where-    nChSS = map toNChildSumStat ts-    sumCh = sum $ map (snd . head) nChSS---- TODO: MERGE WITH SIMULATE.
− src/TLynx/Coalesce/Options.hs
@@ -1,125 +0,0 @@-{-# LANGUAGE DeriveGeneric #-}---- |--- Module      :  TLynx.Coalesce.Options--- Description :  Argument parser for @slynx coalesce@--- Copyright   :  (c) Dominik Schrempf 2020--- License     :  GPL-3.0-or-later------ Maintainer  :  dominik.schrempf@gmail.com--- Stability   :  unstable--- Portability :  portable------ Creation date: Fri May  3 11:51:07 2019.-module TLynx.Coalesce.Options-  ( CoalesceArguments (..),-    coalesceArguments,-    reportCoalesceArguments,-    coalesceFooter,-  )-where--import Data.List-import ELynx.Tools-import Options.Applicative---- | Arguments need to simulate phylogenetic trees using the coalescent process.-data CoalesceArguments = CoalesceArguments-  { -- | Number of simulated trees.-    argsNTrees :: Int,-    -- | Number of leaves per tree.-    argsNLeaves :: Int,-    -- | Perform random sub-sampling with given sampling rate.-    argsRho :: Maybe Double,-    -- | Only print summary statistics?-    argsSumStat :: Bool,-    -- | Seed of NRG, random if 'Nothing'.-    argsSeed :: Seed-  }-  deriving (Eq, Show, Generic)--instance Reproducible CoalesceArguments where-  inFiles _ = []-  outSuffixes _ = [".tree"]-  getSeed = Just . argsSeed-  setSeed a s = a {argsSeed = Fixed s}-  parser = coalesceArguments-  cmdName = "coalesce"-  cmdDsc =-    [ "Simulate phylogenetic trees using the coalescent processes (see also the 'simulate' command for simulations using the birth and death process)."-    ]-  cmdFtr = coalesceFooter--instance FromJSON CoalesceArguments--instance ToJSON CoalesceArguments---- | Print useful information about the provided arguments.-reportCoalesceArguments :: CoalesceArguments -> String-reportCoalesceArguments a =-  intercalate-    "\n"-    [ "Number of simulated trees: " ++ show (argsNTrees a),-      "Number of leaves per tree: " ++ show (argsNLeaves a),-      "Sub-sampling: " ++ ssStr,-      "Summary statistics only: " ++ show (argsSumStat a)-    ]-  where-    ssStr = case argsRho a of-      Nothing -> "No"-      Just r -> "Yes, with probability " <> show r---- | Command line parser.-coalesceArguments :: Parser CoalesceArguments-coalesceArguments =-  CoalesceArguments-    <$> nTreeOpt-    <*> nLeavesOpt-    <*> rhoOpt-    <*> sumStatOpt-    <*> seedOpt--nTreeOpt :: Parser Int-nTreeOpt =-  option auto $-    long "nTrees"-      <> short 't'-      <> metavar "INT"-      <> value 10-      <> showDefault-      <> help "Number of trees"--nLeavesOpt :: Parser Int-nLeavesOpt =-  option auto $-    long "nLeaves"-      <> short 'n'-      <> metavar "INT"-      <> value 5-      <> showDefault-      <> help "Number of leaves per tree"--rhoOpt :: Parser (Maybe Double)-rhoOpt =-  optional $-    option auto $-      long "sub-sample"-        <> short 'r'-        <> metavar "DOUBLE"-        <> showDefault-        <> help-          "Perform random sub-sampling with sampling probability rho (see below)"--sumStatOpt :: Parser Bool-sumStatOpt =-  switch $-    long "summary-statistics" <> short 's' <> showDefault-      <> help-        "Only output number of children for each branch"---- | And a footer.-coalesceFooter :: [String]-coalesceFooter =-  [ "Sub-sampling: simulate one big tree with n'=round(n/rho), n'>=n, leaves, and randomly sample sub-trees with n leaves. Hence, with rho=1.0, the same tree is reported over and over again.",-    "Summary statistics: only print (NumberOfExtantChildren BranchLength) pairs for each branch of each tree. The trees are separated by a newline character."-  ]
src/TLynx/Compare/Compare.hs view
@@ -48,7 +48,7 @@   FilePath ->   ELynx     CompareArguments-    (Tree PhyloStrict BS.ByteString, Tree PhyloStrict BS.ByteString)+    (Tree PhyloExplicit BS.ByteString, Tree PhyloExplicit BS.ByteString) treesOneFile tf = do   nwF <- argsNewickFormat . local <$> ask   $(logInfo) $ T.pack $ "Parse file '" ++ tf ++ "'."@@ -59,8 +59,8 @@     GT -> error "Too many trees in file."     EQ ->       return-        ( either error id $ toStrictTree $ head ts,-          either error id $ toStrictTree $ head . tail $ ts+        ( either error id $ toExplicitTree $ head ts,+          either error id $ toExplicitTree $ head . tail $ ts         )  treesTwoFiles ::@@ -68,14 +68,14 @@   FilePath ->   ELynx     CompareArguments-    (Tree PhyloStrict BS.ByteString, Tree PhyloStrict BS.ByteString)+    (Tree PhyloExplicit BS.ByteString, Tree PhyloExplicit BS.ByteString) treesTwoFiles tf1 tf2 = do   nwF <- argsNewickFormat . local <$> ask   $(logInfo) $ T.pack $ "Parse first tree file '" ++ tf1 ++ "'."   t1 <- liftIO $ parseFileWith (oneNewick nwF) tf1   $(logInfo) $ T.pack $ "Parse second tree file '" ++ tf2 ++ "'."   t2 <- liftIO $ parseFileWith (oneNewick nwF) tf2-  return (either error id $ toStrictTree t1, either error id $ toStrictTree t2)+  return (either error id $ toExplicitTree t1, either error id $ toExplicitTree t2)  -- | More detailed comparison of two trees. compareCmd :: ELynx CompareArguments ()@@ -93,9 +93,9 @@       error         "Need two input files with one tree each or one input file with two trees."   liftIO $ hPutStrLn outH "Tree 1:"-  liftIO $ BL.hPutStrLn outH $ toNewick $ fromStrictTree tr1+  liftIO $ BL.hPutStrLn outH $ toNewick $ toPhyloTree tr1   liftIO $ hPutStrLn outH "Tree 2:"-  liftIO $ BL.hPutStrLn outH $ toNewick $ fromStrictTree tr2+  liftIO $ BL.hPutStrLn outH $ toNewick $ toPhyloTree tr2   liftIO $ hPutStrLn outH ""   -- Intersect trees.   (t1, t2) <-@@ -103,8 +103,8 @@       then do         let [x, y] = either error id $ intersect [tr1, tr2]         liftIO $ hPutStrLn outH "Intersected trees are:"-        liftIO $ BL.hPutStrLn outH $ toNewick $ fromStrictTree x-        liftIO $ BL.hPutStrLn outH $ toNewick $ fromStrictTree y+        liftIO $ BL.hPutStrLn outH $ toNewick $ toPhyloTree x+        liftIO $ BL.hPutStrLn outH $ toNewick $ toPhyloTree y         return (x, y)       else return (tr1, tr2)   -- Check input (moved to library functions).@@ -135,8 +135,8 @@   let t1' = normalizeBranchSupport t1       t2' = normalizeBranchSupport t2   $(logDebug) "Trees with normalized branch support values:"-  $(logDebug) $ E.decodeUtf8 $ BL.toStrict $ toNewick $ fromStrictTree t1'-  $(logDebug) $ E.decodeUtf8 $ BL.toStrict $ toNewick $ fromStrictTree t2'+  $(logDebug) $ E.decodeUtf8 $ BL.toStrict $ toNewick $ toPhyloTree t1'+  $(logDebug) $ E.decodeUtf8 $ BL.toStrict $ toNewick $ toPhyloTree t2'   liftIO $     T.hPutStrLn outH $       formatD
src/TLynx/Connect/Connect.hs view
@@ -122,7 +122,7 @@   (tl, tr) <- parseTrees l r   let ts = connectTrees tl tr   $(logInfo) $ "Connected trees: " <> tShow (length ts)-  liftIO $ BL.hPutStr h $ BL.unlines $ map (toNewick . lengthToPhyloTree) ts+  liftIO $ BL.hPutStr h $ BL.unlines $ map (toNewick . measurableToPhyloTree) ts  connectAndFilter ::   Handle -> FilePath -> FilePath -> FilePath -> ELynx ConnectArguments ()@@ -138,4 +138,4 @@       ts' = filter (compatibleWith getName cs) ts   $(logInfo) $ "Connected  trees: " <> tShow (length ts)   $(logInfo) $ "Compatible trees: " <> tShow (length ts')-  liftIO $ BL.hPutStr h $ BL.unlines $ map (toNewick . lengthToPhyloTree) ts'+  liftIO $ BL.hPutStr h $ BL.unlines $ map (toNewick . measurableToPhyloTree) ts'
src/TLynx/Distance/Distance.hs view
@@ -35,7 +35,6 @@ import qualified Data.ByteString.Lazy.Char8 as BL import Data.List hiding (intersect) import Data.Maybe--- TODO: Remove text dependency. import qualified Data.Text as T import qualified Data.Text.IO as T import qualified Data.Text.Lazy as LT@@ -179,14 +178,14 @@         IncompatibleSplit val -> collapse val . normalizeBranchSupport         _ -> id   -- The trees can be prepared now.-  let trees' :: Forest PhyloStrict BS.ByteString-      trees' = map (collapseF . normalizeF . either error id . toStrictTree) trees+  let trees' :: Forest PhyloExplicit BS.ByteString+      trees' = map (collapseF . normalizeF . either error id . toExplicitTree) trees   $(logDebug) "The prepared trees are:"-  $(logDebug) $ LT.toStrict $ LT.decodeUtf8 $ BL.unlines $ map (toNewick . fromStrictTree) trees'+  $(logDebug) $ LT.toStrict $ LT.decodeUtf8 $ BL.unlines $ map (toNewick . toPhyloTree) trees'   let dsTriplets = case mtree of         Nothing -> pairwise distanceMeasure trees'         Just pt ->-          let t = (either error id . toStrictTree) pt+          let t = (either error id . toExplicitTree) pt            in [(0, i, distanceMeasure t t') | (i, t') <- zip [1 ..] trees']       ds = map (\(_, _, x) -> x) dsTriplets       dsVec = V.fromList ds
src/TLynx/Options.hs view
@@ -20,7 +20,6 @@  import ELynx.Tools import Options.Applicative-import TLynx.Coalesce.Options import TLynx.Compare.Options import TLynx.Connect.Options import TLynx.Distance.Options@@ -31,8 +30,7 @@  -- | The different TLynx commands and their arguments. data CommandArguments-  = Coalesce CoalesceArguments-  | Compare CompareArguments+  = Compare CompareArguments   | Connect ConnectArguments   | Distance DistanceArguments   | Examine ExamineArguments@@ -41,7 +39,6 @@   deriving (Eq, Show, Generic)  instance Reproducible CommandArguments where-  inFiles (Coalesce a) = inFiles a   inFiles (Compare a) = inFiles a   inFiles (Connect a) = inFiles a   inFiles (Distance a) = inFiles a@@ -49,7 +46,6 @@   inFiles (Shuffle a) = inFiles a   inFiles (Simulate a) = inFiles a -  outSuffixes (Coalesce a) = outSuffixes a   outSuffixes (Compare a) = outSuffixes a   outSuffixes (Connect a) = outSuffixes a   outSuffixes (Distance a) = outSuffixes a@@ -57,7 +53,6 @@   outSuffixes (Shuffle a) = outSuffixes a   outSuffixes (Simulate a) = outSuffixes a -  getSeed (Coalesce a) = getSeed a   getSeed (Compare a) = getSeed a   getSeed (Connect a) = getSeed a   getSeed (Distance a) = getSeed a@@ -65,7 +60,6 @@   getSeed (Shuffle a) = getSeed a   getSeed (Simulate a) = getSeed a -  setSeed (Coalesce a) = Coalesce . setSeed a   setSeed (Compare a) = Compare . setSeed a   setSeed (Connect a) = Connect . setSeed a   setSeed (Distance a) = Distance . setSeed a@@ -87,9 +81,6 @@  instance ToJSON CommandArguments -coalesceCommand :: Mod CommandFields CommandArguments-coalesceCommand = createCommandReproducible Coalesce- compareCommand :: Mod CommandFields CommandArguments compareCommand = createCommandReproducible Compare @@ -111,8 +102,7 @@ commandArguments :: Parser CommandArguments commandArguments =   hsubparser $-    coalesceCommand-      <> compareCommand+    compareCommand       <> connectCommand       <> distanceCommand       <> examineCommand
src/TLynx/Shuffle/Shuffle.hs view
@@ -69,7 +69,7 @@     Random -> error "Seed not available; please contact maintainer."     Fixed s -> liftIO $ initialize s   ts <- liftIO $ shuffleT (nReplicates l) (height t) cs ls gen-  liftIO $ BL.hPutStr h $ BL.unlines $ map (toNewick . lengthToPhyloTree) ts+  liftIO $ BL.hPutStr h $ BL.unlines $ map (toNewick . measurableToPhyloTree) ts   liftIO $ hClose h  shuffleT ::
src/TLynx/Simulate/Options.hs view
@@ -1,4 +1,5 @@ {-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE OverloadedStrings #-}  -- | -- Module      :  TLynx.Simulate.Options@@ -12,10 +13,11 @@ -- -- Creation date: Fri May  3 11:51:07 2019. module TLynx.Simulate.Options-  ( SimulateArguments (..),+  ( Height (..),+    Process (..),+    SimulateArguments (..),     simulateArguments,     reportSimulateArguments,-    simulateFooter,   ) where @@ -23,25 +25,59 @@ import ELynx.Tools import Options.Applicative +-- | Condition on tree height (origin or most recent common ancestor).+data Height = Origin Double | Mrca Double+  deriving (Eq, Generic)++instance Show Height where+  show (Origin o) = "Condition on height of origin: " ++ show o+  show (Mrca m) = "Condition on height of MRCA: " ++ show m++instance FromJSON Height++instance ToJSON Height++-- | Process to be used for simulation.+data Process+  = BirthDeath+      { -- | Birth rate.+        bdLambda :: Double,+        -- | Death rate.+        bdMu :: Double,+        -- | Sampling rate.+        bdRho :: Maybe Double,+        -- | Condition on height?+        bdHeight :: Maybe Height+      }+  | Coalescent+  deriving (Eq, Show, Generic)++instance FromJSON Process++instance ToJSON Process++reportProcess :: Process -> String+reportProcess (BirthDeath l m mr mh) =+  intercalate+    "\n"+    [ "Model: Birth and death process",+      "  Birth rate: " ++ show l,+      "  Death rate: " ++ show m,+      "  Sampling probability: " ++ maybe "1.0" show mr,+      "  Height specification: " ++ maybe "Random" show mh+    ]+reportProcess Coalescent = "Model: Coalescent process"+ -- | Arguments need to simulate phylogenetic trees using birth and death processes. data SimulateArguments = SimulateArguments   { -- | Simulated trees.     argsNTrees :: Int,     -- | Number of leaves.     argsNLeaves :: Int,-    -- | Tree height (time to origin or MRCA).-    argsHeight :: Maybe Double,-    -- | False: condition on origin; True:-    -- condition on MRCA.-    argsConditionMRCA :: Bool,-    -- | Birth rate.-    argsLambda :: Double,-    -- | Death rate.-    argsMu :: Double,-    -- | Smapling rate.-    argsRho :: Double,-    -- | Perform actual sub-sampling.-    argsSubSample :: Bool,+    -- | Process.+    argsProcess :: Process,+    -- | Perform sub-sampling with given probability.+    argsSubSample :: Maybe Double,     -- | Only print summary statistics?     argsSumStat :: Bool,     -- | Seed of NRG, random if 'Nothing'.@@ -56,9 +92,7 @@   setSeed a s = a {argsSeed = Fixed s}   parser = simulateArguments   cmdName = "simulate"-  cmdDsc =-    [ "Simulate phylogenetic trees using birth and death processes (see also the 'coalesce' command for simulations using the coalescent process)."-    ]+  cmdDsc = ["Simulate phylogenetic trees using a birth and death or coalescent process."]   cmdFtr = simulateFooter  instance FromJSON SimulateArguments@@ -72,21 +106,14 @@     "\n"     [ "Number of simulated trees: " ++ show (argsNTrees a),       "Number of leaves per tree: " ++ show (argsNLeaves a),-      "Height of trees: " ++ hStr,-      "Birth rate: " ++ show (argsLambda a),-      "Death rate: " ++ show (argsMu a),-      "Sampling probability: " ++ show (argsRho a),-      "Perform sub-sampling: " ++ show (argsSubSample a),+      reportProcess (argsProcess a),+      "Perform sub-sampling: " ++ ssStr,       "Summary statistics only: " ++ show (argsSumStat a)     ]   where-    hStr = case argsHeight a of-      Nothing -> "Random height of origin"-      Just h ->-        show h ++ ", conditioned on "-          ++ if argsConditionMRCA a-            then "MRCA"-            else "origin"+    ssStr = case argsSubSample a of+      Nothing -> "No"+      Just p -> "Yes, with probability " ++ show p  -- | Command line parser. simulateArguments :: Parser SimulateArguments@@ -94,11 +121,7 @@   SimulateArguments     <$> nTreeOpt     <*> nLeavesOpt-    <*> treeHeightOpt-    <*> conditionMRCAOpt-    <*> lambdaOpt-    <*> muOpt-    <*> rhoOpt+    <*> process     <*> subSampleOpt     <*> sumStatOpt     <*> seedOpt@@ -109,8 +132,6 @@     long "nTrees"       <> short 't'       <> metavar "INT"-      <> value 10-      <> showDefault       <> help "Number of trees"  nLeavesOpt :: Parser Int@@ -119,34 +140,14 @@     long "nLeaves"       <> short 'n'       <> metavar "INT"-      <> value 5-      <> showDefault       <> help "Number of leaves per tree" -treeHeightOpt :: Parser (Maybe Double)-treeHeightOpt =-  optional $-    option auto $-      long "height"-        <> short 'H'-        <> metavar "DOUBLE"-        <> help "Fix tree height (no default)"--conditionMRCAOpt :: Parser Bool-conditionMRCAOpt =-  switch $-    long "condition-on-mrca" <> short 'M' <> showDefault-      <> help-        "Do not condition on height of origin but on height of MRCA"- lambdaOpt :: Parser Double lambdaOpt =   option auto $     long "lambda"       <> short 'l'       <> metavar "DOUBLE"-      <> value 1.0-      <> showDefault       <> help "Birth rate lambda"  muOpt :: Parser Double@@ -155,40 +156,85 @@     long "mu"       <> short 'm'       <> metavar "DOUBLE"-      <> value 0.9-      <> showDefault       <> help "Death rate mu"  rhoOpt :: Parser Double rhoOpt =   option auto $-    long "rho" <> short 'r' <> metavar "DOUBLE" <> value 1.0-      <> help-        "Sampling probability rho (default: 1.0)"+    long "rho"+      <> short 'r'+      <> metavar "DOUBLE"+      <> help "Sampling probability rho" -subSampleOpt :: Parser Bool+mrca :: Parser Height+mrca =+  Mrca+    <$> option+      auto+      ( long "mrca"+          <> metavar "DOUBLE"+          <> help "Condition on height of most recent common ancestor"+      )++origin :: Parser Height+origin =+  Origin+    <$> option+      auto+      ( long "origin"+          <> metavar "DOUBLE"+          <> help "Condition on height of origin"+      )++birthDeath :: Parser Process+birthDeath = BirthDeath <$> lambdaOpt <*> muOpt <*> optional rhoOpt <*> optional (mrca <|> origin)++coalescent :: Parser Process+coalescent = pure Coalescent++process :: Parser Process+process =+  hsubparser $+    ( command+        "birthdeath"+        ( info+            birthDeath+            ( progDesc "Birth and death process"+                <> footer "Height: 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 assuming a uniform prior."+            )+        )+        <> command+          "coalescent"+          (info coalescent (progDesc "Coalescent process"))+    )+      <> metavar "PROCESS"+      <> commandGroup "Available processes:"++subSampleOpt :: Parser (Maybe Double) subSampleOpt =-  switch $-    long "sub-sample" <> short 'u' <> showDefault-      <> help-        "Perform sub-sampling; see below."+  optional $+    option auto $+      long+        "sub-sample"+        <> short 'u'+        <> metavar "DOUBLE"+        <> showDefault+        <> help "Perform sub-sampling; see below."  sumStatOpt :: Parser Bool sumStatOpt =   switch $     long "summary-statistics" <> short 's' <> showDefault       <> help-        "Only output number of children for each branch"+        "For each branch, print length and number of children" -citation :: String-citation =-  "Gernhard, T. (2008). The conditioned reconstructed process. Journal of Theoretical Biology, 253(4), 769–778. http://doi.org/10.1016/j.jtbi.2008.04.005"+-- citation :: String+-- citation =+--   "Gernhard, T. (2008). The conditioned reconstructed process. Journal of Theoretical Biology, 253(4), 769–778. http://doi.org/10.1016/j.jtbi.2008.04.005"  -- | And a footer. simulateFooter :: [String] simulateFooter =-  [ "Height of Trees: if no tree height is given, the heights will be randomly drawn from the expected distribution given the number of leaves, the birth and the death rate.",-    "Summary statistics only: only print (NumberOfExtantChildren BranchLength) pairs for each branch of each tree. The trees are separated by a newline character.",-    "Sub-sampling: simulate one big tree with n'=round(n/rho), n'>=n, leaves, and randomly sample sub-trees with n leaves. Hence, with rho=1.0, the same tree is reported over and over again.",-    citation+  [ "See, for example, 'tlynx simulate birthdeath --help'.",+    "Sub-sample with probability p:\n  1. Simulate one big tree with n'=round(n/p), n'>=n, leaves;\n  2. Randomly sample sub-trees with n leaves.\n  (With p=1.0, the same tree is reported over and over again.)"   ]
src/TLynx/Simulate/Simulate.hs view
@@ -33,11 +33,11 @@ import Control.Monad import Control.Monad.IO.Class import Control.Monad.Logger-import Control.Monad.Trans.Reader (ask)+import Control.Monad.Trans.Reader hiding (local) import Control.Parallel.Strategies import qualified Data.ByteString.Builder as BB import qualified Data.ByteString.Lazy.Char8 as BL-import Data.Foldable (toList)+import Data.Foldable import Data.Maybe import qualified Data.Sequence as Seq import qualified Data.Set as Set@@ -45,28 +45,17 @@ import qualified Data.Text.Lazy as LT import qualified Data.Text.Lazy.Encoding as LT import ELynx.Data.Tree-import ELynx.Export.Tree.Newick (toNewick)-import ELynx.Simulate.PointProcess-  ( TimeSpec,-    simulateNReconstructedTrees,-    simulateReconstructedTree,-  )+import ELynx.Export.Tree.Newick+import qualified ELynx.Simulate.Coalescent as CS+import qualified ELynx.Simulate.PointProcess as PP import ELynx.Tools import System.Random.MWC-  ( GenIO,-    initialize,-  ) import TLynx.Simulate.Options --- | Simulate phylogenetic trees.+-- | Simulate phylogenetic trees using birth and death process. simulate :: ELynx SimulateArguments () simulate = do-  l <- local <$> ask-  let SimulateArguments nTrees nLeaves tHeight mrca lambda mu rho subS sumS (Fixed s) =-        l-  -- error "simulate: seed not available; please contact maintainer."-  when (isNothing tHeight && mrca) $-    error "Cannot condition on MRCA (-M) when height is not given (-H)."+  l@(SimulateArguments nTrees nLeaves pr subS sumS (Fixed s)) <- local <$> ask   c <- liftIO getNumCapabilities   logNewSection "Arguments"   $(logInfo) $ T.pack $ reportSimulateArguments l@@ -74,64 +63,88 @@   $(logInfo) $ T.pack $ "Number of used cores: " <> show c   gs <- liftIO $ initialize s >>= \gen -> splitGen c gen   let chunks = getChunks c nTrees-      timeSpec = fmap (,mrca) tHeight-  trs <--    if subS-      then-        simulateAndSubSampleNTreesConcurrently-          nLeaves-          lambda-          mu-          rho-          timeSpec-          chunks-          gs-      else simulateNTreesConcurrently nLeaves lambda mu rho timeSpec chunks gs+  trs <- case pr of+    (BirthDeath lambda mu mRho mHeight) -> do+      let rho = fromMaybe 1.0 mRho+          -- This is bad code, but I don't want to change the definition of 'TimeSpec'.+          timeSpec = case mHeight of+            Nothing -> Nothing+            Just (Mrca h) -> Just (h, True)+            Just (Origin h) -> Just (h, False)+      case subS of+        Nothing -> liftIO $ bdSimulateNTreesConcurrently nLeaves lambda mu rho timeSpec chunks gs+        Just p ->+          bdSimulateAndSubSampleNTreesConcurrently+            nLeaves+            lambda+            mu+            rho+            p+            timeSpec+            chunks+            gs+    Coalescent -> case subS of+      Nothing -> liftIO $ coalSimulateNTreesConcurrently nLeaves chunks gs+      Just p -> coalSimulateAndSubSampleNTreesConcurrently nLeaves p chunks gs   let ls =         if sumS           then parMap rpar (formatNChildSumStat . toNChildSumStat) trs-          else parMap rpar toNewick $ map lengthToPhyloTree trs+          else parMap rpar toNewick $ map measurableToPhyloTree trs   let res = BL.unlines ls   out "simulated trees" res ".tree" -simulateNTreesConcurrently ::+bdSimulateNTreesConcurrently ::   Int ->   Double ->   Double ->   Double ->-  TimeSpec ->+  PP.TimeSpec ->   [Int] ->   [GenIO] ->-  ELynx SimulateArguments (Forest Length Int)-simulateNTreesConcurrently nLeaves l m r timeSpec chunks gs = do+  IO (Forest Length Int)+bdSimulateNTreesConcurrently nLeaves l m r timeSpec chunks gs = do   let l' = l * r       m' = m - l * (1.0 - r)   trss <--    liftIO $+    mapConcurrently+      (\(n, g) -> PP.simulateNReconstructedTrees n nLeaves timeSpec l' m' g)+      (zip chunks gs)+  return $ concat trss++coalSimulateNTreesConcurrently ::+  Int ->+  [Int] ->+  [GenIO] ->+  IO (Forest Length Int)+coalSimulateNTreesConcurrently nL chunks gs = do+  trss <-       mapConcurrently-        (\(n, g) -> simulateNReconstructedTrees n nLeaves timeSpec l' m' g)+        (\(n, g) -> replicateM n $ CS.simulate nL g)         (zip chunks gs)   return $ concat trss -simulateAndSubSampleNTreesConcurrently ::+bdSimulateAndSubSampleNTreesConcurrently ::   Int ->   Double ->   Double ->   Double ->-  TimeSpec ->+  Double ->+  PP.TimeSpec ->   [Int] ->   [GenIO] ->   ELynx SimulateArguments (Forest Length Int)-simulateAndSubSampleNTreesConcurrently nLeaves l m r timeSpec chunks gs = do-  let nLeavesBigTree = (round $ fromIntegral nLeaves / r) :: Int+bdSimulateAndSubSampleNTreesConcurrently nLeaves l m r p timeSpec chunks gs = do+  let nLeavesBigTree = (round $ fromIntegral nLeaves / p) :: Int+      l' = l * r+      m' = m - l * (1.0 - r)   logNewSection $     T.pack $       "Simulate one big tree with "         <> show nLeavesBigTree         <> " leaves."-  tr <- liftIO $ simulateReconstructedTree nLeavesBigTree timeSpec l m (head gs)+  tr <- liftIO $ PP.simulateReconstructedTree nLeavesBigTree timeSpec l' m' (head gs)   -- Log the base tree.-  $(logInfo) $ LT.toStrict $ LT.decodeUtf8 $ toNewick $ lengthToPhyloTree tr+  $(logInfo) $ LT.toStrict $ LT.decodeUtf8 $ toNewick $ measurableToPhyloTree tr   logNewSection $     T.pack $       "Sub sample "@@ -148,7 +161,39 @@   let trs = catMaybes $ concat trss   return $ map prune trs --- | Extract a random subtree with @N@ leaves of a tree with @M@ leaves, where+coalSimulateAndSubSampleNTreesConcurrently ::+  Int ->+  Double ->+  [Int] ->+  [GenIO] ->+  ELynx SimulateArguments (Forest Length Int)+coalSimulateAndSubSampleNTreesConcurrently nL p chunks gs = do+  let nLeavesBigTree = (round $ fromIntegral nL / p) :: Int+  logNewSection $+    T.pack $+      "Simulate one big tree with "+        <> show nLeavesBigTree+        <> " leaves."+  tr <- liftIO $ CS.simulate nLeavesBigTree (head gs)+  -- Log the base tree.+  $(logInfo) $ LT.toStrict $ LT.decodeUtf8 $ toNewick $ measurableToPhyloTree tr+  logNewSection $+    T.pack $+      "Sub sample "+        <> show (sum chunks)+        <> " trees with "+        <> show nL+        <> " leaves."+  let lvs = Seq.fromList $ leaves tr+  trss <-+    liftIO $+      mapConcurrently+        (\(nSamples, g) -> nSubSamples nSamples lvs nL tr g)+        (zip chunks gs)+  let trs = catMaybes $ concat trss+  return $ map prune trs++-- Extract a random subtree with @N@ leaves of a tree with @M@ leaves, where -- @M>N@ (otherwise error). The complete list of leaves (names are assumed to be -- unique) has to be provided as a 'Seq.Seq', and a 'Seq.Set', so that fast -- sub-sampling as well as lookup are fast and so that these data structures do@@ -170,18 +215,20 @@     let lsSets = map Set.fromList lss     return [dropLeavesWith (`Set.notMember` ls) tree | ls <- lsSets] --- | Pair of branch length with number of extant children.+-- Pair of branch length with number of extant children. type BrLnNChildren = (BranchLength, Int) --- | Possible summary statistic of phylogenetic trees. A list of tuples+-- Possible summary statistic of phylogenetic trees. A list of tuples -- (BranchLength, NumberOfExtantChildrenBelowThisBranch). type NChildSumStat = [BrLnNChildren] --- | Format the summary statistics in the following form:+-- Format the summary statistics in the following form:+-- -- @---    nLeaves1 branchLength1---    nLeaves2 branchLength2---    ....+--   nLeaves1 branchLength1+--   nLeaves2 branchLength2+--   ....+-- @ formatNChildSumStat :: NChildSumStat -> BL.ByteString formatNChildSumStat s =   BB.toLazyByteString . mconcat $ map formatNChildSumStatLine s@@ -190,7 +237,7 @@ formatNChildSumStatLine (l, n) =   BB.intDec n <> BB.char8 ' ' <> BB.doubleDec l <> BB.char8 '\n' --- | Compute NChilSumStat for a phylogenetic tree.+-- Compute NChilSumStat for a phylogenetic tree. toNChildSumStat :: Measurable e => Tree e a -> NChildSumStat toNChildSumStat (Node br _ []) = [(getLen br, 1)] toNChildSumStat (Node br _ ts) = (getLen br, sumCh) : concat nChSS
src/TLynx/TLynx.hs view
@@ -16,7 +16,6 @@ where  import ELynx.Tools-import TLynx.Coalesce.Coalesce import TLynx.Compare.Compare import TLynx.Connect.Connect import TLynx.Distance.Distance@@ -30,8 +29,6 @@ -- | Run TLynx with given arguments. tlynx :: Arguments CommandArguments -> IO () tlynx c = case local c of-  Coalesce _ ->-    eLynxWrapper c (\(Arguments g (Coalesce l)) -> Arguments g l) coalesce   Compare _ ->     eLynxWrapper c (\(Arguments g (Compare l)) -> Arguments g l) compareCmd   Connect _ ->
tlynx.cabal view
@@ -1,6 +1,6 @@ cabal-version:  2.2 name:           tlynx-version:        0.3.3+version:        0.3.4 synopsis:       Handle phylogenetic trees description:    Examine, compare, and simulate phylogenetic trees in a reproducible way. Please see the README on GitHub at <https://github.com/dschrempf/elynx>. category:       Bioinformatics@@ -23,8 +23,6 @@  library   exposed-modules:-      TLynx.Coalesce.Coalesce-      TLynx.Coalesce.Options       TLynx.Compare.Compare       TLynx.Compare.Options       TLynx.Connect.Connect