diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -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. diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -2,7 +2,7 @@ # The ELynx Suite -Version: 0.3.3. +Version: 0.3.4. Reproducible evolution made easy.

@@ -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 diff --git a/elynx-tree.cabal b/elynx-tree.cabal --- a/elynx-tree.cabal +++ b/elynx-tree.cabal @@ -1,6 +1,6 @@ cabal-version: 2.2 name: elynx-tree -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 . category: Bioinformatics @@ -101,6 +101,7 @@ ELynx.Data.Tree.PartitionSpec ELynx.Data.Tree.PhylogenySpec ELynx.Data.Tree.RootedSpec + ELynx.Data.Tree.SupportedSpec ELynx.Export.Tree.NewickSpec ELynx.Export.Tree.NexusSpec ELynx.Import.Tree.NewickSpec diff --git a/src/ELynx/Data/Tree/Measurable.hs b/src/ELynx/Data/Tree/Measurable.hs --- a/src/ELynx/Data/Tree/Measurable.hs +++ b/src/ELynx/Data/Tree/Measurable.hs @@ -9,6 +9,10 @@ -- Portability : portable -- -- Creation date: Thu Jan 17 14:16:34 2019. +-- +-- Non-negativity of branch lengths is not (yet) ensured. To ensure +-- non-negativity, a newtype wrapper could be used, but this would be a major +-- refactor. module ELynx.Data.Tree.Measurable ( BranchLength, Measurable (..), @@ -73,12 +77,6 @@ rootHeight :: Measurable e => Tree e a -> BranchLength rootHeight (Node _ _ []) = 0 rootHeight t = maximum $ concatMap distancesOriginLeaves (forest t) - --- -- | The maximum distance between root and leaves. --- -- --- -- Does not include the length of the stem, see 'height'. --- rootHeight :: Measurable e => Tree e a -> BranchLength --- rootHeight t = subtract (getStem t) $ height t -- | Distances from the origin of a tree to the leaves. -- diff --git a/src/ELynx/Data/Tree/Phylogeny.hs b/src/ELynx/Data/Tree/Phylogeny.hs --- a/src/ELynx/Data/Tree/Phylogeny.hs +++ b/src/ELynx/Data/Tree/Phylogeny.hs @@ -31,12 +31,13 @@ -- Note: 'Tree's encoded in Newick format correspond to rooted trees. By -- convention only, a tree parsed from Newick format is usually thought to be -- unrooted, when the root node is multifurcating and has three children. This --- convention is not enforced here. Newick trees are just parsed as they are, --- and a rooted tree is returned. +-- convention is not used here. Newick trees are just parsed as they are, and a +-- rooted tree is returned. -- --- The bifurcating root of a tree can be changed with 'roots', or 'rootAt'. +-- The bifurcating root of a tree can be changed with 'rootAt' or 'midpoint'; a +-- list of all rooted trees is returned by 'roots'. -- --- Trees with multifurcating root nodes can be properly rooted using 'outgroup'. +-- Trees with multifurcating root nodes can be rooted using 'outgroup'. module ELynx.Data.Tree.Phylogeny ( -- * Functions equal, @@ -49,15 +50,16 @@ -- * Branch labels Phylo (..), + toPhyloTree, + measurableToPhyloTree, + supportedToPhyloTree, Length (..), phyloToLengthTree, - lengthToPhyloTree, Support (..), phyloToSupportTree, phyloToSupportTreeUnsafe, - PhyloStrict (..), - toStrictTree, - fromStrictTree, + PhyloExplicit (..), + toExplicitTree, ) where @@ -122,20 +124,22 @@ bifurcating (Node _ _ [x, y]) = bifurcating x && bifurcating y bifurcating _ = False --- TODO. Adapt for trees with branches. +-- I believe that manual treatment with 'outgroup' is preferable. + -- -- | Remove multifurcations. -- -- --- -- A caterpillar like bifurcating tree is used to resolve all multifurcations on --- -- a tree. The multifurcating nodes are copied. +-- -- A caterpillar like bifurcating structure is used to resolve all +-- -- multifurcations on a tree. -- -- --- -- Branch labels are not handled. --- resolve :: Tree e a -> Tree e a +-- -- Multifurcating nodes are copied and branches are 'split'. +-- resolve :: Splittable e => Tree e a -> Tree e a -- resolve t@(Node _ _ []) = t --- resolve (Node _ l [x]) = Node () l [resolve x] --- resolve (Node _ l [x, y]) = Node () l $ map resolve [x, y] --- resolve (Node _ l (x : xs)) = Node () l $ map resolve [x, Node () l xs] +-- resolve (Node br lb [x]) = Node br lb [resolve x] +-- resolve (Node br lb [x, y]) = Node br lb $ map resolve [x, y] +-- resolve (Node br lb (Node brL lbL xsL : xs)) = Node br lb [Node brL' lbL (map resolve xsL), Node brL' lb (map resolve xs)] +-- where brL' = split brL --- | Resolve a multifurcation at the root using an outgroup. +-- | Resolve a multifurcating root using an outgroup. -- -- A bifurcating root node with the provided label is introduced. The affected -- branch is 'split'. @@ -163,11 +167,11 @@ -- Introduce a bifurcating root node. t' = Node b r [Node (split brO) lbO tsO, Node (split brO) l (tail ts)] --- XXX: The 'midpoint' algorithm is pretty stupid because it calculates all --- rooted trees and then finds the one minimizing the difference between the --- heights of the left and right sub tree. Actually, one just needs to move left --- or right, with the aim to minimize the height difference between the left and --- right sub tree. +-- The 'midpoint' algorithm is pretty stupid because it calculates all rooted +-- trees and then finds the one minimizing the difference between the heights of +-- the left and right sub tree. Actually, one just needs to move left or right, +-- with the aim to minimize the height difference between the left and right sub +-- tree. -- | Root tree at the midpoint. -- @@ -323,6 +327,33 @@ instance FromJSON Phylo +-- | Set all branch length and support values to 'Just' the value. +-- +-- Useful to export a tree with branch lengths in Newick format. +toPhyloTree :: (Measurable e, Supported e) => Tree e a -> Tree Phylo a +toPhyloTree = first toPhyloLabel + +toPhyloLabel :: (Measurable e, Supported e) => e -> Phylo +toPhyloLabel x = Phylo (Just $ getLen x) (Just $ getSup x) + +-- | Set all branch support values to 'Nothing'. +-- +-- Useful to export a tree with branch lengths to Newick format. +measurableToPhyloTree :: Measurable e => Tree e a -> Tree Phylo a +measurableToPhyloTree = first measurableToPhyloLabel + +measurableToPhyloLabel :: Measurable e => e -> Phylo +measurableToPhyloLabel x = Phylo (Just $ getLen x) Nothing + +-- | Set all branch lengths to 'Nothing'. +-- +-- Useful to export a tree with branch support to Newick format. +supportedToPhyloTree :: Supported e => Tree e a -> Tree Phylo a +supportedToPhyloTree = first supportedToPhyloLabel + +supportedToPhyloLabel :: Supported e => e -> Phylo +supportedToPhyloLabel x = Phylo Nothing (Just $ getSup x) + -- | Branch length label. -- -- For conversion, see 'phyloToLengthTree' and 'lengthToPhyloTree'. @@ -358,15 +389,6 @@ toLength :: Phylo -> Maybe Length toLength p = Length <$> brLen p --- | Set all branch support values to 'Nothing'. --- --- Useful, for example, to export a tree with branch lengths in Newick format. -lengthToPhyloTree :: Tree Length a -> Tree Phylo a -lengthToPhyloTree = first fromLengthLabel - -fromLengthLabel :: Length -> Phylo -fromLengthLabel (Length b) = Phylo (Just b) Nothing - -- | Branch support label. -- -- For conversion, see 'phyloToSupportTree'. @@ -426,51 +448,40 @@ cleanSupport :: BranchSupport -> Tree Phylo a -> Tree Support a cleanSupport maxSup (Node (Phylo _ s) l xs) = Node (Support $ fromMaybe maxSup s) l $ map (cleanSupport maxSup) xs --- XXX: Probably change name. Strict is reserved for "not lazy". - --- | Strict branch label for phylogenetic trees. -data PhyloStrict = PhyloStrict +-- | Explicit branch label for phylogenetic trees. +data PhyloExplicit = PhyloExplicit { sBrLen :: BranchLength, sBrSup :: BranchSupport } deriving (Read, Show, Eq, Ord, Generic) -instance Semigroup PhyloStrict where - PhyloStrict bL sL <> PhyloStrict bR sR = PhyloStrict (bL + bR) (min sL sR) +instance Semigroup PhyloExplicit where + PhyloExplicit bL sL <> PhyloExplicit bR sR = PhyloExplicit (bL + bR) (min sL sR) -instance Measurable PhyloStrict where +instance Measurable PhyloExplicit where getLen = sBrLen setLen b l = l {sBrLen = b} -instance Splittable PhyloStrict where +instance Splittable PhyloExplicit where split l = l {sBrLen = b'} where b' = sBrLen l / 2.0 -instance Supported PhyloStrict where +instance Supported PhyloExplicit where getSup = sBrSup setSup s l = l {sBrSup = s} -instance ToJSON PhyloStrict +instance ToJSON PhyloExplicit -instance FromJSON PhyloStrict +instance FromJSON PhyloExplicit --- | Conversion to a 'PhyloStrict' tree. +-- | Conversion to a 'PhyloExplicit' tree. -- -- See 'phyloToLengthTree' and 'phyloToSupportTree'. -toStrictTree :: Tree Phylo a -> Either String (Tree PhyloStrict a) -toStrictTree t = do +toExplicitTree :: Tree Phylo a -> Either String (Tree PhyloExplicit a) +toExplicitTree t = do lt <- first fromLength <$> phyloToLengthTree t st <- first fromSupport <$> phyloToSupportTree t - case zipTreesWith PhyloStrict const lt st of - Nothing -> error "toStrictTree: This is a bug. Can not zip two trees with the same topology." + case zipTreesWith PhyloExplicit const lt st of + Nothing -> error "toExplicitTree: This is a bug. Can not zip two trees with the same topology." Just zt -> return zt - --- | Set all branch length and support values to 'Just' the value. --- --- Useful, for example, to export a tree with branch lengths in Newick format. -fromStrictTree :: Tree PhyloStrict a -> Tree Phylo a -fromStrictTree = first fromStrictLabel - -fromStrictLabel :: PhyloStrict -> Phylo -fromStrictLabel (PhyloStrict b s) = Phylo (Just b) (Just s) diff --git a/src/ELynx/Data/Tree/Supported.hs b/src/ELynx/Data/Tree/Supported.hs --- a/src/ELynx/Data/Tree/Supported.hs +++ b/src/ELynx/Data/Tree/Supported.hs @@ -9,6 +9,10 @@ -- Portability : portable -- -- Creation date: Thu Jun 13 14:06:45 2019. +-- +-- Non-negativity of branch support values is not (yet) ensured. To ensure +-- non-negativity, a newtype wrapper could be used, but this would be a major +-- refactor. module ELynx.Data.Tree.Supported ( BranchSupport, Supported (..), @@ -41,9 +45,6 @@ where m = bimaximum $ bimap getSup (const 0) t --- TODO: Something was wrong here. @collapse 1.0 t@ should be a star tree but it --- was a leaf. Is this still so? - -- | Collapse branches with support lower than given value. -- -- The branch and node labels of the collapsed branches are discarded. @@ -52,10 +53,14 @@ let tr' = collapse' th tr in if tr == tr' then tr else collapse th tr' +-- A leaf has full support. +highP :: Supported e => Double -> Tree e a -> Bool +highP _ (Node _ _ []) = True +highP th (Node br _ _) = getSup br >= th + -- See 'collapse'. collapse' :: Supported e => BranchSupport -> Tree e a -> Tree e a -collapse' _ t@(Node _ _ []) = t collapse' th (Node br lb ts) = Node br lb $ map (collapse' th) (highSupport ++ lowSupportForest) where - (highSupport, lowSupport) = partition ((>= th) . getSup . branch) ts + (highSupport, lowSupport) = partition (highP th) ts lowSupportForest = concatMap forest lowSupport diff --git a/test/ELynx/Data/Tree/SupportedSpec.hs b/test/ELynx/Data/Tree/SupportedSpec.hs new file mode 100644 --- /dev/null +++ b/test/ELynx/Data/Tree/SupportedSpec.hs @@ -0,0 +1,42 @@ +{-# LANGUAGE OverloadedStrings #-} + +-- | +-- Module : ELynx.Data.Tree.SupportedSpec +-- Description : Unit tests for ELynx.Data.Tree.SupportedSpec +-- Copyright : (c) Dominik Schrempf, 2020 +-- License : GPL-3.0-or-later +-- +-- Maintainer : dominik.schrempf@gmail.com +-- Stability : unstable +-- Portability : portable +-- +-- Creation date: Fri Aug 21 14:20:09 2020. +module ELynx.Data.Tree.SupportedSpec + ( spec, + ) +where + +import qualified Data.ByteString.Char8 as BS +import ELynx.Data.Tree +import ELynx.Import.Tree.Newick +-- import ELynx.Export.Tree.Newick +import ELynx.Tools +import Test.Hspec + +collapseTree :: Tree Phylo BS.ByteString +collapseTree = parseByteStringWith (oneNewick IqTree) "((a,b),(c,d));" + +collapseStarTree :: Tree Phylo BS.ByteString +collapseStarTree = parseByteStringWith (oneNewick Standard) "(a[1.0],b[1.0],c[1.0],d[1.0])[1.0];" + +spec :: Spec +spec = do + describe "collapse" $ do + it "creates a star tree for 1.0" $ do + let t = phyloToSupportTreeUnsafe collapseTree + s = phyloToSupportTreeUnsafe collapseStarTree + collapse 0 t `shouldBe` t + collapse 0.01 t `shouldBe` t + collapse 0.99 t `shouldBe` t + collapse 1.0 t `shouldBe` t + collapse 1.1 t `shouldBe` s