diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -1,10 +1,10 @@ # Table of Contents -1. [Changelog for ELynx](#org9916021) +1. [Changelog for ELynx](#org263dcac) - + # Changelog for ELynx diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -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) diff --git a/bench/Bench.hs b/bench/Bench.hs new file mode 100644 --- /dev/null +++ b/bench/Bench.hs @@ -0,0 +1,37 @@ +{- | +Module : Bench +Description : Various benchmarks +Copyright : (c) Dominik Schrempf 2019 +License : GPL-3 + +Maintainer : dominik.schrempf@gmail.com +Stability : unstable +Portability : portable + +Creation date: Mon Dec 16 13:33:27 2019. + +-} + +import Criterion.Main +import qualified Data.ByteString.Lazy.Char8 as L +import Data.Tree + +import ELynx.Data.Tree.Bipartition +import ELynx.Data.Tree.PhyloTree +import ELynx.Import.Tree.Newick +import ELynx.Tools.InputOutput + +treeFileMany :: FilePath +treeFileMany = "data/Many.trees" + +getManyTrees :: IO [Tree (PhyloLabel L.ByteString)] +getManyTrees = parseFileWith manyNewick treeFileMany + +main :: IO () +main = do + ts <- getManyTrees + defaultMain + [ bgroup "bipartition" + [ bench "manyTrees" $ whnf (map bipartitions) ts + ] + ] diff --git a/elynx-tree.cabal b/elynx-tree.cabal --- a/elynx-tree.cabal +++ b/elynx-tree.cabal @@ -1,6 +1,6 @@ cabal-version: 1.12 name: elynx-tree -version: 0.0.1 +version: 0.1.0 license: GPL-3 license-file: LICENSE copyright: Dominik Schrempf (2019) @@ -28,8 +28,10 @@ ELynx.Data.Tree.Distance ELynx.Data.Tree.EvoTree ELynx.Data.Tree.MeasurableTree + ELynx.Data.Tree.Multipartition ELynx.Data.Tree.NamedTree ELynx.Data.Tree.PhyloTree + ELynx.Data.Tree.Subset ELynx.Data.Tree.SumStat ELynx.Data.Tree.Tree ELynx.Distribution.BirthDeath @@ -48,15 +50,15 @@ other-modules: Paths_elynx_tree default-language: Haskell2010 - ghc-options: -Wall + ghc-options: -Wall -O2 build-depends: QuickCheck >=2.13.2 && <2.14, base >=4.12.0.0 && <4.13, bytestring >=0.10.8.2 && <0.11, containers >=0.6.0.1 && <0.7, - elynx-tools >=0.0.1 && <0.1, + elynx-tools >=0.1.0 && <0.2, lifted-async >=0.10.0.4 && <0.11, - math-functions >=0.3.1.0 && <0.4, + math-functions >=0.3.2.1 && <0.4, megaparsec >=7.0.5 && <7.1, mwc-random >=0.14.0.0 && <0.15, optparse-applicative >=0.14.3.0 && <0.15, @@ -74,23 +76,53 @@ other-modules: ELynx.Data.Tree.BipartitionSpec ELynx.Data.Tree.DistanceSpec + ELynx.Data.Tree.MultipartitionSpec ELynx.Data.Tree.TreeSpec ELynx.Export.Tree.NewickSpec ELynx.Import.Tree.NewickSpec Paths_elynx_tree default-language: Haskell2010 - ghc-options: -Wall -eventlog -threaded -rtsopts -with-rtsopts=-N + ghc-options: -Wall -O2 -eventlog -threaded -rtsopts + -with-rtsopts=-N build-depends: QuickCheck >=2.13.2 && <2.14, base >=4.12.0.0 && <4.13, bytestring >=0.10.8.2 && <0.11, containers >=0.6.0.1 && <0.7, - elynx-tools >=0.0.1 && <0.1, + elynx-tools >=0.1.0 && <0.2, elynx-tree -any, hspec >=2.7.1 && <2.8, hspec-megaparsec >=2.0.1 && <2.1, lifted-async >=0.10.0.4 && <0.11, - math-functions >=0.3.1.0 && <0.4, + math-functions >=0.3.2.1 && <0.4, + megaparsec >=7.0.5 && <7.1, + mwc-random >=0.14.0.0 && <0.15, + optparse-applicative >=0.14.3.0 && <0.15, + parallel >=3.2.2.0 && <3.3, + primitive >=0.6.4.0 && <0.7, + quickcheck-instances >=0.3.22 && <0.4, + statistics >=0.15.0.0 && <0.16, + transformers >=0.5.6.2 && <0.6, + vector >=0.12.0.3 && <0.13 + +benchmark tree-bench + type: exitcode-stdio-1.0 + main-is: Bench.hs + hs-source-dirs: bench + other-modules: + Paths_elynx_tree + default-language: Haskell2010 + ghc-options: -Wall -O2 -threaded -rtsopts -with-rtsopts=-N + build-depends: + QuickCheck >=2.13.2 && <2.14, + base >=4.12.0.0 && <4.13, + bytestring >=0.10.8.2 && <0.11, + containers >=0.6.0.1 && <0.7, + criterion >=1.5.6.0 && <1.6, + elynx-tools >=0.1.0 && <0.2, + elynx-tree -any, + lifted-async >=0.10.0.4 && <0.11, + math-functions >=0.3.2.1 && <0.4, megaparsec >=7.0.5 && <7.1, mwc-random >=0.14.0.0 && <0.15, optparse-applicative >=0.14.3.0 && <0.15, diff --git a/src/ELynx/Data/Tree/Bipartition.hs b/src/ELynx/Data/Tree/Bipartition.hs --- a/src/ELynx/Data/Tree/Bipartition.hs +++ b/src/ELynx/Data/Tree/Bipartition.hs @@ -27,41 +27,58 @@ module ELynx.Data.Tree.Bipartition ( -- * The 'Bipartition' data type. Bipartition () + , bps , bp , bpmap + , bphuman -- * Working with 'Bipartition's. + , bipartition , bipartitions - , bipartitionToBranch - , bipartitionsCombined + , bipartitionToBranchLength + , compatible ) where --- import Data.List -import qualified Data.Map as M +import qualified Data.Map as M import Data.Maybe -import qualified Data.Set as S +import qualified Data.Set as S import Data.Tree +import ELynx.Data.Tree.Subset import ELynx.Data.Tree.Tree --- | Bipartitions with 'S.Set's, since order of elements within the leaf sets --- is not important. Also the order of the two leaf sets of the bipartition is --- not important (see 'Eq' instance definition). -newtype Bipartition a = Bipartition (S.Set a, S.Set a) +-- | Each branch of a tree partitions the leaves of the tree into two +-- 'Subset's, or a bipartition. Also the order of the two partitions of the +-- 'Bipartition' is not important (see the 'Eq' instance). +newtype Bipartition a = Bipartition {bps :: (Subset a, Subset a) -- ^ Tuple of partitions + } deriving (Show, Read) --- instance Show a => Show (Bipartition a) where --- show (Bipartition (x, y)) = "(" ++ showSet x ++ "|" ++ showSet y ++ ")" --- where showSet s = intercalate "," $ map show $ S.toList s +-- I decided to NOT provide a human readable show instance because in this case, +-- I need the following identity to hold: +-- +-- > read . show = id +-- +-- This identity is met by the derived instance anyways. A more human readable +-- instance would most likely violate the identity. +-- | Show a bipartition in a human readable form. Use a provided function to +-- extract the valuable information. +bphuman :: (a -> String) -> Bipartition a -> String +bphuman f (Bipartition (x, y)) = "(" ++ sshow f x ++ "|" ++ sshow f y ++ ")" + -- | Create a bipartition from two 'S.Set's. -bp :: Ord a => S.Set a -> S.Set a -> Bipartition a -bp x y = if x >= y - then Bipartition (x, y) - else Bipartition (y, x) +bp :: Ord a => Subset a -> Subset a -> Bipartition a +bp xs ys = if xs >= ys + then Bipartition (xs, ys) + else Bipartition (ys, xs) +-- | Map a function over all elements in the 'Bipartition's. +bpmap :: (Ord a, Ord b) => (a -> b) -> Bipartition a -> Bipartition b +bpmap f (Bipartition (x, y)) = bp (smap f x) (smap f y) + -- | Create a bipartition from two 'S.Set's. -bpWith :: (Ord a, Ord b) => (a -> b) -> S.Set a -> S.Set a -> Bipartition b -bpWith f x y = bpmap f $ bp x y +bpwith :: (Ord a, Ord b) => (a -> b) -> Subset a -> Subset a -> Bipartition b +bpwith f x y = bpmap f $ bp x y instance (Eq a) => Eq (Bipartition a) where Bipartition x == Bipartition y = x == y @@ -69,163 +86,94 @@ instance (Ord a) => Ord (Bipartition a) where Bipartition x `compare` Bipartition y = x `compare` y --- | Map a function over all elements in the 'Bipartition's. -bpmap :: (Ord a, Ord b) => (a -> b) -> Bipartition a -> Bipartition b -bpmap f (Bipartition (x, y)) = bp (S.map f x) (S.map f y) - --- | Each node of a tree is root of a subtree. Get the leaves of the subtree of --- each node. -leavesTree :: (Ord a) => Tree a -> Tree (S.Set a) -leavesTree (Node l []) = Node (S.singleton l) [] -leavesTree (Node _ xs) = Node (S.unions $ map rootLabel xs') xs' - where xs' = map leavesTree xs - --- | Loop through each tree in a forest to report the complementary leaf sets. -subForestGetLeafSets :: (Ord a) - => S.Set a -- ^ Complementary leaf set at the stem - -> Tree (S.Set a) -- ^ Tree with leaf set nodes - -> [S.Set a] -subForestGetLeafSets lvsS t = lvsOthers - where - xs = subForest t - nChildren = length xs - lvsChildren = map rootLabel xs - lvsOtherChildren = [ S.unions $ lvsS - : take i lvsChildren ++ drop (i+1) lvsChildren - | i <- [0 .. (nChildren - 1)] ] - lvsOthers = map (S.union lvsS) lvsOtherChildren +-- Check if a bipartition is valid. For now, only checks if one set is empty. +valid :: Bipartition a -> Bool +valid (Bipartition (xs, ys)) = not $ snull xs || snull ys +-- | For a bifurcating root, get the bipartition induced by the root node. +bipartition :: Ord a => Tree a -> Bipartition a +bipartition (Node _ [x, y]) = bp (sfromlist $ leaves x) (sfromlist $ leaves y) +bipartition _ = error "Root node is not bifurcating." --- | Get all bipartitions. +-- | Get all bipartitions of the tree. bipartitions :: Ord a => Tree a -> S.Set (Bipartition a) -bipartitions t = if S.size (S.fromList lvs) == length lvs - then bipartitionsUnsafe t +bipartitions t = if S.size (S.fromList ls) == length ls + then S.filter valid $ bipartitionsUnsafe sempty pTree else error "bipartitions: The tree contains duplicate leaves." - where lvs = leaves t + where ls = leaves t + pTree = partitionTree t -- | See 'bipartitions', but do not check if leaves are unique. -bipartitionsUnsafe :: Ord a => Tree a -> S.Set (Bipartition a) -bipartitionsUnsafe (Node _ [] ) = S.empty --- If the root stem is split by degree two nodes, just go on since the root stem --- does not induce any bipartitions. -bipartitionsUnsafe (Node _ [x]) = bipartitionsUnsafe x --- We have rose trees, so we need to through the list of children and combine --- each of them with the rest. -bipartitionsUnsafe t = - S.unions [ bipartitions' lvs x - | (lvs, x) <- zip lvsOthers (subForest lvsTree) ] - where - lvsTree = leavesTree t - lvsOthers = subForestGetLeafSets S.empty lvsTree - -bipartitions' :: Ord a => S.Set a -> Tree (S.Set a) -> S.Set (Bipartition a) -bipartitions' lvsStem t@(Node lvs xs) - | S.null lvsStem = error "bipartitions': no complementing leaf set." - -- Leaf; return a singleton map; bipartition with the leaf and the rest of the tree. - | null xs = S.singleton $ bp lvsStem lvs - -- For degree two nodes, pass the creation of the set on. - | length xs == 1 = bipartitions' lvsStem (head xs) - -- We have rose trees, so we need to through the list of children and combine - -- each of them with the rest. Also, we use up the possible branch information - -- 'br' and start afresh with 'mempty'. - | otherwise = S.unions $ S.singleton (bp lvsStem lvs) : zipWith bipartitions' lvsOthers xs +bipartitionsUnsafe :: Ord a => Subset a -> Tree (Subset a) -> S.Set (Bipartition a) +bipartitionsUnsafe p (Node l [] ) = S.singleton $ bp p l +-- Degree two nodes do not induce additional bipartitions. +bipartitionsUnsafe p (Node _ [x]) = bipartitionsUnsafe p x +-- Go through the list of children and combine each of them with the rest. +bipartitionsUnsafe p t@(Node ls xs) = + S.unions $ S.singleton (bp p ls) : [ bipartitionsUnsafe lvs x | (lvs, x) <- zip lsOthers xs ] where - lvsOthers = subForestGetLeafSets lvsStem t + lsOthers = subForestGetSubsets p t --- | Each branch on a 'Tree' defines a unique 'Bipartition' of leaves. Convert a --- tree into a 'M.Map' from each 'Bipartition' to the branch inducing the --- respective 'Bipartition'. The information about the branch is extracted from --- the nodes with a given function. If the tree has degree two nodes, the branch --- values are combined; a unity element is required, and so we need the 'Monoid' --- type class constraint. Checks if leaves are unique. -bipartitionToBranch :: (Ord a, Ord b, Monoid c) - => (a -> b) -- ^ Value to compare on - -> (a -> c) -- ^ Convert node to branch length +-- | For a given rose 'Tree', remove all degree two nodes and reconnect the +-- resulting disconnected pairs of branches and sum their branch lengths. For +-- this operation, a combining binary function and a unity element is required, +-- and so we need the 'Monoid' type class constraint. Now, each branch on the +-- tree defines a unique 'Bipartition' of leaves. Convert a tree into a 'M.Map' +-- from each 'Bipartition' to the length of the branch inducing the respective +-- 'Bipartition'. The relevant information about the leaves is extracted from +-- the (leaf) nodes with a given function. Also check if leaves are unique. +bipartitionToBranchLength :: (Ord a, Ord b, Monoid c) + => (a -> b) -- ^ Convert node labels to leaves (usually + -- leaf names) + -> (a -> c) -- ^ Get length of branch attached to node -> Tree a -- ^ Tree to dissect -> M.Map (Bipartition b) c -bipartitionToBranch f g t = if S.size (S.fromList lvs) == length lvs - then bipartitionToBranchUnsafe f g t - else error "bipartitionToBranch: The tree contains duplicate leaves." - where lvs = leaves t +bipartitionToBranchLength f g t = if S.size (S.fromList ls) == length ls + then M.filterWithKey (const . valid) $ + bipartitionToBranchLengthUnsafe (mempty, sempty) f lAndPTree + else error "bipartitionToBranchLength: The tree contains duplicate leaves." + where ls = leaves t + bTree = fmap g t + pTree = partitionTree t + lAndPTree = fromJust $ merge bTree pTree --- | See 'bipartitionToBranch', but does not check if leaves are unique. -bipartitionToBranchUnsafe :: (Ord a, Ord b, Monoid c) - => (a -> b) -- ^ Value to compare on - -> (a -> c) -- ^ Convert node to branch length - -> Tree a -- ^ Tree to dissect - -> M.Map (Bipartition b) c -bipartitionToBranchUnsafe _ _ (Node _ [] ) = M.empty --- If the root stem is split by degree two nodes, just go on and ignore the --- branch information, because the stem does not induce any bipartition --- anyways.. -bipartitionToBranchUnsafe f g (Node _ [x]) = bipartitionToBranchUnsafe f g x --- We have rose trees, so we need to through the list of children and combine --- each of them with the rest. -bipartitionToBranchUnsafe f g t = - M.unionsWith (<>) [ bipartitionToBranch' lvs mempty f g x - | (lvs, x) <- zip lvsOthers (subForest nodeAndLeavesTrees) ] - where - lvsTree = leavesTree t - nodeAndLeavesTrees = fromJust $ merge t lvsTree - lvsOthers = subForestGetLeafSets S.empty lvsTree +-- | See 'bipartitionToBranchLength'. When calculating the map, branches +-- separated by various degree two nodes have to be combined. Hence, not only +-- the complementary partition towards the stem, but also the node label itself +-- have to be passed along. +type Info c a = (c, Subset a) +-- | See 'bipartitionToBranchLength', but does not check if leaves are unique. -- We need information about the nodes, and also about the leaves of the induced --- sub trees. Hence, we need a somewhat complicated node type @(a, S.Set a)@. -bipartitionToBranch' :: (Ord a, Ord b, Monoid c) - => S.Set a -- ^ Complementary set of leaves - -- towards the stem - -> c -- ^ Maybe we have to pass along some - -- information from above (degree two - -- nodes) - -> (a -> b) -- ^ Extract value to compare on - -> (a -> c) -- ^ Extract information about branch - -- from node - -> Tree (a, S.Set a) -- ^ Tree containing nodes and sub - -- tree leaf sets - -> M.Map (Bipartition b) c -bipartitionToBranch' lvsStem br f g t@(Node l xs ) - | S.null lvsStem = error "bipartitionToBranch': no complementing leaf set." - -- Leaf; return a singleton map; bipartition with the leaf and the rest of the tree. - | null xs = M.singleton (bpWith f lvsStem lvsThisNode) (br <> g label) - -- Pass the creation of the map entry on, but extend the branch. - | length xs == 1 = bipartitionToBranch' lvsStem (br <> g label) f g (head xs) - -- We have rose trees, so we need to through the list of children and combine - -- each of them with the rest. Also, we use up the possible branch information - -- 'br' and start afresh with 'mempty'. - | otherwise = M.insert (bpWith f lvsStem lvsThisNode) (br <> g label) - $ M.unions [ bipartitionToBranch' lvs mempty f g x - | (lvs, x) <- zip lvsOthers xs ] - where - label = fst l - lvsThisNode = snd l - lvsOthers = subForestGetLeafSets lvsStem $ fmap snd t - --- | Get all bipartitions, but combine leaves from multifurcations. This is --- useful to find incompatible splits. See --- 'ELynx.Data.Tree.Distance.incompatibleSplitsDistance'. Assume that a root --- node with three children is actually not a multifurcation (because then we --- would have no induced bypartitions), but rather corresponds to an unrooted --- tree. -bipartitionsCombined :: (Ord a, Show a) => Tree a -> S.Set (Bipartition a) -bipartitionsCombined t@(Node _ xs) - | null xs = S.empty - | length xs == 1 = bipartitionsCombined (head xs) - -- One big multifurcation does not induce any bipartitions. - | length xs > 3 = S.empty - | otherwise = res +-- sub trees. Hence, we need a somewhat complicated node label type +-- +-- > (a, Subset a) +bipartitionToBranchLengthUnsafe :: (Ord a, Ord b, Monoid c) + => Info c a + -> (a -> b) -- ^ Value to compare on + -> Tree (Info c a) -- ^ Tree to dissect + -> M.Map (Bipartition b) c +bipartitionToBranchLengthUnsafe (l, p) f (Node (l', p') [] ) = + M.singleton (bpwith f p p') (l <> l') +-- The branch length has to be added for degree two nodes. +bipartitionToBranchLengthUnsafe (l, p) f (Node (l', _ ) [x]) = + bipartitionToBranchLengthUnsafe (l <> l', p) f x +-- Go through the list of children and combine each of them with the rest. +bipartitionToBranchLengthUnsafe (l, p) f t@(Node (l', p') xs) = + M.unionsWith (<>) $ + M.singleton (bpwith f p p') (l <> l') : + [ bipartitionToBranchLengthUnsafe (mempty, lvs) f x | (lvs, x) <- zip lvsOthers xs ] where - res = S.unions [ bipartitionsCombined' lvs x - | (lvs, x) <- zip lvsOthers (subForest lvsTree) ] - lvsTree = leavesTree t - lvsOthers = subForestGetLeafSets S.empty lvsTree + lvsOthers = subForestGetSubsets p (fmap snd t) -bipartitionsCombined' :: Ord a => S.Set a -> Tree (S.Set a) -> S.Set (Bipartition a) -bipartitionsCombined' lvsStem t@(Node lvs xs) - | S.null lvsStem = error "bipartitionsCombined': no complementing leaf set." - | null xs = S.singleton $ bp lvsStem lvs - | length xs == 1 = bipartitionsCombined' lvsStem (head xs) - | length xs == 2 = S.unions $ - S.singleton (bp lvsStem lvs) : zipWith bipartitionsCombined' lvsOthers xs - | otherwise = S.singleton $ bp lvsStem lvs - where - lvsOthers = subForestGetLeafSets lvsStem t +-- | Determine compatibility between an bipartition and a subset. If both +-- subsets of the bipartition share elements with the given subset, the +-- bipartition is incompatible with this subset. If all elements of the subset +-- are either not in the bipartition or mapping to one of the two subsets of the +-- bipartition, the bipartition and the subset are compatible. See also +-- 'ELynx.Data.Tree.Multipartition.compatible'. +compatible :: (Show a, Ord a) => Bipartition a -> Subset a -> Bool +-- compatible (Bipartition (l, r)) ss = sintersection l ss `sdisjoint` sintersection r ss +compatible (Bipartition (l, r)) ss = snull lOverlap || snull rOverlap + where lOverlap = sintersection l ss + rOverlap = sintersection r ss diff --git a/src/ELynx/Data/Tree/BranchSupportTree.hs b/src/ELynx/Data/Tree/BranchSupportTree.hs --- a/src/ELynx/Data/Tree/BranchSupportTree.hs +++ b/src/ELynx/Data/Tree/BranchSupportTree.hs @@ -14,51 +14,63 @@ module ELynx.Data.Tree.BranchSupportTree ( BranchSupport - , BranchSupportLabel (..) + , BranchSupported (..) , normalize , collapse ) where import Data.List -import Data.Maybe import Data.Tree --- XXX : This is probably the preferred way. --- data BranchSupport = --- BSNothing --- | BSInt Int --- | BSDouble Double --- deriving (Num) +-- TODO: Implementation of a proper data type. When parsing Newick files, we +-- need to check if an Integer can be parsed, otherwise parse a Double. Also, if +-- one branch support value is a Double, all others have to be converted to +-- Double, so that we do not have a mixed data type. In this case, it may be +-- better to use a type family? Or just use Maybe a, and don't restrict the a +-- type? +-- +-- data BranchSupport = BSNothing +-- | BSInt Int +-- | BSDouble Double deriving (Num) -- | At the moment, just use 'Double'. It would be preferable to use a wrapper -- data type that can handle 'Int' or 'Double'. type BranchSupport = Maybe Double --- | A label that supports branch support values. -class BranchSupportLabel a where +-- | A label that supports extraction and setting of branch support values. +class BranchSupported a where -- | For now, branch support is a Double, but one could also think about -- bootstrap values, which are integers. getBranchSupport :: a -> BranchSupport setBranchSupport :: BranchSupport -> a -> a -apply :: BranchSupportLabel a => (Double -> Double) -> a -> a +apply :: BranchSupported a => (Double -> Double) -> a -> a apply f l = setBranchSupport (f <$> s) l where s = getBranchSupport l -- | Normalize branch support values. The maximum branch support value will be -- set to 1.0. -normalize :: BranchSupportLabel a => Tree a -> Tree a -normalize t = if isNothing m then t else fmap (apply (/ fromJust m)) t - where m = maximum $ fmap getBranchSupport t +normalize :: BranchSupported a => Tree a -> Tree a +normalize t = case mm of + Nothing -> t + Just m -> fmap (apply (/ m)) t + where mm = maximum $ fmap getBranchSupport t accept :: Double -> Maybe Double -> Bool accept _ Nothing = True accept thresh (Just s) = s > thresh -- | Collapse branches with support lower than given value. Note, branch length --- is ignored at the moment. -collapse :: BranchSupportLabel a => Double -> Tree a -> Tree a -collapse _ n@(Node _ []) = n -collapse thresh (Node l xs) = Node l $ map (collapse thresh) (highS ++ lowSubForest) - where (highS, lowS) = partition (accept thresh . getBranchSupport . rootLabel) xs +-- is ignored at the moment. Continue collapsing until a fix point is reached. +collapse :: (Show a, Eq a, BranchSupported a) => Double -> Tree a -> Tree a +collapse th tr = if tr == tr' + then tr + else collapse th tr' + where tr' = collapse' th tr + +-- | See 'collapse'. +collapse' :: BranchSupported a => Double -> Tree a -> Tree a +collapse' _ t@(Node _ []) = t +collapse' th (Node l xs) = Node l $ map (collapse' th) (highS ++ lowSubForest) + where (highS, lowS) = partition (accept th . getBranchSupport . rootLabel) xs lowSubForest = concatMap subForest lowS diff --git a/src/ELynx/Data/Tree/Distance.hs b/src/ELynx/Data/Tree/Distance.hs --- a/src/ELynx/Data/Tree/Distance.hs +++ b/src/ELynx/Data/Tree/Distance.hs @@ -16,14 +16,14 @@ -} module ELynx.Data.Tree.Distance - ( symmetricDistance - , symmetricDistanceWith - , incompatibleSplitsDistance - , incompatibleSplitsDistanceWith - , computePairwiseDistances - , computeAdjacentDistances - , branchScoreDistance - , branchScoreDistanceWith + ( symmetric + , symmetricWith + , incompatibleSplits + , incompatibleSplitsWith + , branchScore + , branchScoreWith + , pairwise + , adjacent ) where import Data.List @@ -32,93 +32,125 @@ import qualified Data.Set as S import Data.Tree -import ELynx.Data.Tree.Bipartition -import ELynx.Data.Tree.MeasurableTree +import ELynx.Data.Tree.Bipartition (Bipartition, + bipartitionToBranchLength, + bipartitions) +import ELynx.Data.Tree.MeasurableTree (Measurable, getLen) +import ELynx.Data.Tree.Multipartition (Multipartition, compatible, + fromBipartition, + multipartitions) import ELynx.Data.Tree.NamedTree --- -- Difference between two 'Set's, see 'Set.difference'. Do not compare elements --- -- directly but apply a function beforehand. --- differenceWith :: (Ord a, Ord b) => (a -> b) -> Set.Set a -> Set.Set a -> Set.Set a --- differenceWith f xs ys = Set.filter (\e -> f e `Set.notMember` ys') xs --- where ys' = Set.map f ys - --- -- Symmetric difference between two 'Set's. Do not compare elements directly but --- -- apply a function beforehand. --- symmetricDifferenceWith :: (Ord a, Ord b) => (a -> b) -> Set.Set a -> Set.Set a -> Set.Set a --- symmetricDifferenceWith f xs ys = xsNotInYs `Set.union` ysNotInXs --- where --- xsNotInYs = differenceWith f xs ys --- ysNotInXs = differenceWith f ys xs - -- Symmetric difference between two 'Set's. -symmetricDifferenceS :: Ord a => S.Set a -> S.Set a -> S.Set a -symmetricDifferenceS xs ys = S.difference xs ys `S.union` S.difference ys xs - --- -- Symmetric difference between two 'Map's. --- symmetricDifferenceM :: Ord k => M.Map k a -> M.Map k a -> M.Map k a --- symmetricDifferenceM x y = M.difference x y `M.union` M.difference y x +symmetricDifference :: Ord a => S.Set a -> S.Set a -> S.Set a +symmetricDifference xs ys = S.difference xs ys `S.union` S.difference ys xs -- | Symmetric (Robinson-Foulds) distance between two trees. Before comparing -- the leaf labels, apply a given function. This is useful, for example, to --- compare the labels of 'Named' trees on their names only. The tree is assumed --- to be UNROOTED! +-- compare the labels of 'ELynx.Data.Tree.NamedTree.Named' trees on their names +-- only. The tree is assumed to be UNROOTED! -- -- XXX: Comparing a list of trees with this function recomputes bipartitions. -symmetricDistanceWith :: (Ord b) => (a -> b) -> Tree a -> Tree a -> Int -symmetricDistanceWith f t1 t2 = length $ symmetricDifferenceS (bs t1) (bs t2) +symmetricWith :: (Ord b) => (a -> b) -> Tree a -> Tree a -> Int +symmetricWith f t1 t2 = length $ symmetricDifference (bs t1) (bs t2) where bs t = bipartitions $ fmap f t --- | See 'symmetricDistanceWith', but with 'id' for comparisons. -symmetricDistance :: Ord a => Tree a -> Tree a -> Int -symmetricDistance = symmetricDistanceWith id +-- | See 'symmetricWith', but with 'id' for comparisons. +symmetric :: (Ord a, Named a) => Tree a -> Tree a -> Int +symmetric = symmetricWith getName --- | Number of incompatible splits. Similar to 'symmetricDistanceWith' but --- merges multifurcations. +countIncompatibilities :: (Ord a, Show a) => S.Set (Bipartition a) -> S.Set (Multipartition a) -> Int +countIncompatibilities bs ms = foldl' (\i b -> if any (compatible (fromBipartition b)) ms + then i + else i+1) 0 bs + +-- | Number of incompatible splits. Similar to 'symmetricWith' but all +-- bipartition induced by multifurcations are considered. For a detailed +-- description of how the distance is calculated, see 'compatible'. -- +-- A multifurcation on a tree may (but not necessarily does) represent missing +-- information about the order of bifurcations. In this case, it is interesting +-- to get a set of compatible bifurcations of the tree. For example, the tree +-- +-- > (A,(B,C,D)) +-- +-- induces the following bipartitions: +-- +-- > A|BCD +-- > B|ACD +-- +-- > C|ABD +-- > D|ABC +-- +-- Those are also reported by the function 'bipartitions'. However, the tree is +-- additionally compatible with the following hidden bipartitions: +-- +-- > AB|CD +-- > AC|BD +-- > AD|BC +-- +-- For an explanation of how compatibility of a bipartition with a +-- multipartition is checked, see 'compatible'. Before using 'compatible', +-- bipartitions are simply converted to multipartitions with two subsets. +-- +-- Only if a bipartition is not compatible with all induced multifurcations of +-- the other tree, it is incompatible. +-- -- XXX: Comparing a list of trees with this function recomputes bipartitions. -incompatibleSplitsDistanceWith :: (Ord b, Show b) => (a -> b) -> Tree a -> Tree a -> Int -incompatibleSplitsDistanceWith f t1 t2 = length $ symmetricDifferenceS (ms t1) (ms t2) - where ms t = bipartitionsCombined $ fmap f t +incompatibleSplitsWith :: (Ord b, Show b) => (a -> b) -> Tree a -> Tree a -> Int +incompatibleSplitsWith f t1 t2 = countIncompatibilities putIncBs1 ms2 + + countIncompatibilities putIncBs2 ms1 + where + -- Bipartitions. + bs t = bipartitions $ fmap f t + bs1 = bs t1 + bs2 = bs t2 + -- Putative incompatible bipartitions of trees one and two, respectively. + putIncBs1 = bs1 S.\\ bs2 + putIncBs2 = bs2 S.\\ bs1 + -- Multipartitions. + ms t = multipartitions $ fmap f t + ms1 = ms t1 + ms2 = ms t2 --- | See 'incompatibleSplitsDistanceWith', use 'id' for comparisons. -incompatibleSplitsDistance :: (Ord a, Show a) => Tree a -> Tree a -> Int -incompatibleSplitsDistance = incompatibleSplitsDistanceWith id +-- | See 'incompatibleSplitsWith', use 'id' for comparisons. +incompatibleSplits :: (Ord a, Named a) => Tree a -> Tree a -> Int +incompatibleSplits = incompatibleSplitsWith getName -- | Compute branch score distance between two trees. Before comparing the leaf -- labels, apply a function. This is useful, for example, to compare the labels --- of 'Named' trees on their names only. The branch information which is --- compared to compute the distance is extracted from the nodes with a given --- function. Assumes that the trees are UNROOTED. +-- of 'ELynx.Data.Tree.NamedTree.Named' trees on their names only. The branch +-- information which is compared to compute the distance is extracted from the +-- nodes with a given function. Assumes that the trees are UNROOTED. -- -- XXX: Comparing a list of trees with this function recomputes bipartitions. -branchScoreDistanceWith :: (Ord a, Ord b, Floating c) +branchScoreWith :: (Ord a, Ord b, Floating c) => (a -> b) -- ^ Label to compare on -> (a -> c) -- ^ Branch information (e.g., length) -- associated with a node -> Tree a -> Tree a -> c -branchScoreDistanceWith f g t1 t2 = sqrt dsSquared - where bs = bipartitionToBranch f (Sum . g) +branchScoreWith f g t1 t2 = sqrt dsSquared + where bs = bipartitionToBranchLength f (Sum . g) dBs = M.map getSum $ M.unionWith (-) (bs t1) (bs t2) dsSquared = foldl' (\acc e -> acc + e*e) 0 dBs --- | See 'branchScoreDistanceWith', use 'id' for comparisons. -branchScoreDistance :: (Ord a, Measurable a, Named a) => Tree a -> Tree a -> Double -branchScoreDistance = branchScoreDistanceWith getName getLen +-- | See 'branchScoreWith', use 'id' for comparisons. +branchScore :: (Ord a, Named a, Measurable a) => Tree a -> Tree a -> Double +branchScore = branchScoreWith getName getLen -- | Compute pairwise distances of a list of input trees. Use given distance -- measure. Returns a triple, the first two elements are the indices of the -- compared trees, the third is the distance. -computePairwiseDistances :: (a -> a -> b) -- ^ Distance function - -> [a] -- ^ Input trees - -> [(Int, Int, b)] -- ^ (index i, index j, distance i j) -computePairwiseDistances dist trs = [ (i, j, dist x y) - | (i:is, x:xs) <- zip (tails [0..]) (tails trs) - , (j, y) <- zip is xs ] +pairwise :: (a -> a -> b) -- ^ Distance function + -> [a] -- ^ Input trees + -> [(Int, Int, b)] -- ^ (index i, index j, distance i j) +pairwise dist trs = [ (i, j, dist x y) + | (i:is, x:xs) <- zip (tails [0..]) (tails trs) + , (j, y) <- zip is xs ] -- | Compute distances between adjacent pairs of a list of input trees. Use -- given distance measure. -computeAdjacentDistances :: (Tree a -> Tree a -> b) -- ^ Distance function - -> [Tree a] -- ^ Input trees - -> [b] -computeAdjacentDistances dist trs = [ dist x y | (x, y) <- zip trs (tail trs) ] - +adjacent :: (Tree a -> Tree a -> b) -- ^ Distance function + -> [Tree a] -- ^ Input trees + -> [b] +adjacent dist trs = [ dist x y | (x, y) <- zip trs (tail trs) ] diff --git a/src/ELynx/Data/Tree/EvoTree.hs b/src/ELynx/Data/Tree/EvoTree.hs --- a/src/ELynx/Data/Tree/EvoTree.hs +++ b/src/ELynx/Data/Tree/EvoTree.hs @@ -28,6 +28,7 @@ internal :: n -> Bool internal n = not $ extant n || extinct n + external :: n -> Bool external = not . internal diff --git a/src/ELynx/Data/Tree/MeasurableTree.hs b/src/ELynx/Data/Tree/MeasurableTree.hs --- a/src/ELynx/Data/Tree/MeasurableTree.hs +++ b/src/ELynx/Data/Tree/MeasurableTree.hs @@ -1,3 +1,5 @@ +{-# LANGUAGE OverloadedStrings #-} + {- | Module : ELynx.Data.Tree.MeasurableTree Description : Functions on trees with branch lengths @@ -16,14 +18,18 @@ module ELynx.Data.Tree.MeasurableTree ( Measurable (..) , distancesRootLeaves - , averageDistanceRootLeaves + , distancesOriginLeaves + , averageDistanceOriginLeaves , height - , lengthenRoot - , shortenRoot + , rootHeight + , lengthenStem + , shortenStem , summarize , totalBranchLength , normalize , prune + , removeMultifurcations + , ultrametric ) where import qualified Data.ByteString.Lazy.Char8 as L @@ -31,6 +37,7 @@ import Data.Tree import ELynx.Data.Tree.Tree +import ELynx.Tools.Equality (allNearlyEqual) -- | A 'Node' label with measurable and modifiable branch length to the parent. class Measurable a where @@ -39,40 +46,60 @@ -- | Set attached branch length. setLen :: Double -> a -> a - -- | Elongate branch length. - lengthen :: Double -> a -> a - lengthen dl l = setLen (dl + getLen l) l +-- | Elongate branch length. +lengthen :: Measurable a => Double -> a -> a +lengthen dl l = setLen (dl + getLen l) l - -- | Shorten branch length. - shorten :: Double -> a -> a - shorten dl = lengthen (-dl) +-- -- | Shorten branch length. +-- shorten :: Double -> a -> a +-- shorten dl = lengthen (-dl) --- | Distances from the root of a tree to its leaves. +-- | Distances from the root node of a tree to its leaves (this are not the +-- distances from the origin to the leaves, see 'distancesOriginLeaves'). distancesRootLeaves :: (Measurable a) => Tree a -> [Double] -distancesRootLeaves (Node l []) = [getLen l] -distancesRootLeaves (Node l f ) = concatMap (map (+ getLen l) . distancesRootLeaves) f +distancesRootLeaves (Node _ []) = [0] +distancesRootLeaves (Node _ f ) = concat [map (+ getLen (rootLabel d)) (distancesRootLeaves d) | d <- f] --- | Average distance from the root of a tree to its leaves. -averageDistanceRootLeaves :: (Measurable a) => Tree a -> Double -averageDistanceRootLeaves tr = sum ds / fromIntegral n - where ds = distancesRootLeaves tr +-- -- | Distances from the origin of a tree to its leaves (this is not the distance +-- -- from the root node to the leaves, which would be @distanceOriginLeaves t - +-- -- (getLen $ rootLabel t)@.). +-- distancesOriginLeaves :: (Measurable a) => Tree a -> [Double] +-- distancesOriginLeaves (Node l []) = [getLen l] +-- distancesOriginLeaves (Node l f ) = concatMap (map (+ getLen l) . distancesOriginLeaves) f + +-- | Distances from the origin of a tree to its leaves (this is not the distance +-- from the root node to the leaves, which would be @distanceOriginLeaves t - +-- (getLen $ rootLabel t)@.). See also 'distancesRootLeaves'. +distancesOriginLeaves :: (Measurable a) => Tree a -> [Double] +distancesOriginLeaves t@(Node l _) = map (+ getLen l) (distancesRootLeaves t) + +-- | Average distance from the origin of a tree to its leaves, see +-- 'distancesOriginLeaves'. +averageDistanceOriginLeaves :: (Measurable a) => Tree a -> Double +averageDistanceOriginLeaves tr = sum ds / fromIntegral n + where ds = distancesOriginLeaves tr n = length ds --- | Height of a tree. Returns 0 if the tree is empty. +-- | Height (max distance between origin and leaves) of a tree. Return 0 if the +-- tree is empty. height :: (Measurable a) => Tree a -> Double -height = maximum . distancesRootLeaves +height = maximum . distancesOriginLeaves +-- | Height of root node. Return 0 if the tree is empty. +rootHeight :: (Measurable a) => Tree a -> Double +rootHeight = maximum . distancesRootLeaves + -- | Lengthen the distance between root and origin. -lengthenRoot :: (Measurable a) => Double -> Tree a -> Tree a -lengthenRoot dl (Node lbl chs) = Node (lengthen dl lbl) chs +lengthenStem :: (Measurable a) => Double -> Tree a -> Tree a +lengthenStem dl (Node lbl chs) = Node (lengthen dl lbl) chs -- | Lengthen the distance between root and origin. -shortenRoot :: (Measurable a) => Double -> Tree a -> Tree a -shortenRoot dl = lengthenRoot (-dl) +shortenStem :: (Measurable a) => Double -> Tree a -> Tree a +shortenStem dl = lengthenStem (-dl) -- | Summarize a tree with measureable branch lengths. summarize :: (Measurable a) => Tree a -> L.ByteString -summarize t = L.unlines $ map L.pack +summarize t = L.intercalate "\n" $ map L.pack [ "Leaves: " ++ show n ++ "." , "Height: " ++ show h ++ "." , "Average distance root to leaves: " ++ show h' ++ "." @@ -80,7 +107,7 @@ where n = length . leaves $ t h = height t b = totalBranchLength t - h' = sum (distancesRootLeaves t) / fromIntegral n + h' = sum (distancesOriginLeaves t) / fromIntegral n -- | Total branch length of a tree. totalBranchLength :: (Measurable a) => Tree a -> Double @@ -96,3 +123,16 @@ prune :: (Measurable a) => Tree a -> Tree a prune = pruneWith f where f da pa = lengthen (getLen pa) da + +-- | Remove multifurcations by copying multifurcating nodes and introducing +-- branches of length 0. +removeMultifurcations :: Measurable a => Tree a -> Tree a +removeMultifurcations t@(Node _ [] ) = t +removeMultifurcations (Node l [x]) = Node l [removeMultifurcations x] +removeMultifurcations (Node l [x, y]) = Node l $ map removeMultifurcations [x, y] +removeMultifurcations (Node l (x:xs)) = Node l $ map removeMultifurcations [x, Node l' xs] + where l' = setLen 1.0 l + +-- | Check if a tree is ultrametric. +ultrametric :: Measurable a => Tree a -> Bool +ultrametric = allNearlyEqual . distancesOriginLeaves diff --git a/src/ELynx/Data/Tree/Multipartition.hs b/src/ELynx/Data/Tree/Multipartition.hs new file mode 100644 --- /dev/null +++ b/src/ELynx/Data/Tree/Multipartition.hs @@ -0,0 +1,162 @@ +{- | +Module : ELynx.Data.Tree.Multipartition +Description : Multipartitions on rose trees +Copyright : (c) Dominik Schrempf 2019 +License : GPL-3 + +Maintainer : dominik.schrempf@gmail.com +Stability : unstable +Portability : portable + +Creation date: Thu Dec 12 12:58:49 2019. + +A multifurcation induces a 'Multipartition', similar to branches inducing +'ELynx.Data.Tree.Bipartition's. + +-} + +module ELynx.Data.Tree.Multipartition + ( -- * The 'Multipartition' data type. + Multipartition () + , mps + , mp + , mpmap + , mphuman + , fromBipartition + -- * Working with 'Multipartition's. + , multipartitions + , findSubset + , compatible + ) where + +import Data.List (filter, find, foldl', intercalate) +import Data.Maybe (fromMaybe) +import qualified Data.Set as S +import Data.Tree + +import ELynx.Data.Tree.Bipartition (Bipartition, bps) +import ELynx.Data.Tree.Subset +import ELynx.Data.Tree.Tree + +-- | Each branch of a bifurcating tree partitions the leaves of the tree into +-- three 'Subset's, see 'ELynx.Data.Tree.Bipartition'. In a similar way, each +-- internal node induces a tripartition. Tripartitions are not yet implemented +-- (December 2019) because it is usually sufficient to work with bipartitions. +-- If, however, the tree is multifurcating and a specific node has more than two +-- children, the number of subsets induced by this node is larger than three, +-- and a 'Multipartition'. Multipartitions are interesting in that we can use +-- them for calculating incompatible splits, see 'ELynx.Data.Tree.Distance'. The +-- order of the partitions within a multipartition is unimportant, . +newtype Multipartition a = Multipartition { mps :: S.Set (Subset a) -- ^ Set of partitions + } + deriving (Show, Read) + +-- | Show a multipartition in a human readable form. Use a provided function to +-- extract the valuable information. +mphuman :: (a -> String) -> Multipartition a -> String +mphuman f (Multipartition xs) = "(" ++ intercalate "|" (map (sshow f) (S.toList xs)) ++ ")" + +-- | Create a multipartition. +mp :: Ord a => [Subset a] -> Multipartition a +mp = mp' . filter (not . snull) + +mp' :: Ord a => [Subset a] -> Multipartition a +-- XXX: For now also allow multipartitions with no, one, or two elements. +-- mp' [] = error "mp': Cannot create multipartition from empty list." +-- mp' [_] = error "mp': Cannot create multipartition from list with one element." +-- mp' [_, _] = error "mp': Cannot create multipartition from list with two elements." +mp' xs = Multipartition (S.fromList xs) + +-- | Map a function over all elements in the multipartitions. +mpmap :: (Ord a, Ord b) => (a -> b) -> Multipartition a -> Multipartition b +mpmap f (Multipartition xs) = Multipartition $ S.map (smap f) xs + +instance (Eq a) => Eq (Multipartition a) where + Multipartition xs == Multipartition ys = xs == ys + +instance (Ord a) => Ord (Multipartition a) where + Multipartition xs `compare` Multipartition ys = xs `compare` ys + +-- | Convert bipartition to multipartition. +fromBipartition :: Ord a => Bipartition a -> Multipartition a +fromBipartition bp = mp [l, r] + where (l, r) = bps bp + +-- | Get all multipartitions of a tree. +multipartitions :: Ord a => Tree a -> S.Set (Multipartition a) +multipartitions t = if S.size (S.fromList lvs) == length lvs + then multipartitionsUnsafe sempty (partitionTree t) + else error "multipartitions: The tree contains duplicate leaves." + where lvs = leaves t + +-- | See 'multipartitions', but do not check if leaves are unique. +multipartitionsUnsafe :: Ord a => Subset a -> Tree (Subset a) -> S.Set (Multipartition a) +multipartitionsUnsafe _ (Node _ [] ) = S.empty +multipartitionsUnsafe xs (Node _ [x] ) = multipartitionsUnsafe xs x +multipartitionsUnsafe xs (Node _ [x, y]) = S.union l r + where l = multipartitionsUnsafe (sunion xs (rootLabel x)) y + r = multipartitionsUnsafe (sunion xs (rootLabel y)) x +multipartitionsUnsafe xs t@(Node _ ys ) = S.unions $ + S.singleton (mp (xs : map rootLabel ys)) + : zipWith multipartitionsUnsafe lvsOthers ys + where lvsOthers = subForestGetSubsets xs t + +-- | Find the multipartition containing a given element. +findSubset :: Ord a => a -> Multipartition a -> Subset a +findSubset l m = fromMaybe + -- Return the empty subset if nothing is found. This corresponds + -- to having no information about the leaf in question. + sempty + (find (smember l) ss) + where ss = mps m + +-- Add the subset of a bipartition which contains a given element. +addSubset :: Ord a => Multipartition a -> S.Set (Subset a) -> a -> S.Set (Subset a) +addSubset m ss l = if not $ snull ss' + then ss' `S.insert` ss + else ss + where ss' = findSubset l m + +-- Each subset overlaps with a number of subsets of a bipartition which are +-- returned by this function. +overlap :: Ord a => Multipartition a -> Subset a -> S.Set (Subset a) +overlap m = foldl' (addSubset m) S.empty + +-- | Multipartitions are compatible if they do not contain conflicting +-- information. This function checks if two multipartitions are compatible with +-- each other. Thereby, following algorithm is used: +-- +-- 1. Take each subset of the first multipartition. +-- +-- 2a. Determine the overlap: For each leaf of the chosen subset, add the subset +-- of the second multipartition containing the leaf. The result is a set of +-- subsets, which is the union of the added subsets. +-- +-- The data type "set of subsets" is actually the same data type as a +-- multipartition. However, it is not a partition, because it may and will not +-- span the whole set of leaves, and so, I use @S.Set (Subset a)@. One could +-- define a multiset data type to improve comprehensibility. +-- +-- 2b. Collect the set of subsets from point 1. +-- +-- 3. Each set of subsets needs to be either equal or disjoint with any other +-- set of subsets in the collection. If so, the first multipartition is +-- compatible with the second. +-- +-- 4. Exchange the first with the second multipartition and go through steps 1 +-- to 3. +-- +-- See also 'ELynx.Data.Tree.Bipartition.compatible'. +compatible :: (Ord a, Show a) => Multipartition a -> Multipartition a -> Bool +-- compatible l r = traceShow lOverlaps $ traceShow rOverlaps $ and $ +compatible l r = and $ + [x `S.disjoint` y | x <- lOverlaps, y <- lOverlaps, x /= y] ++ + [x `S.disjoint` y | x <- rOverlaps, y <- rOverlaps, x /= y] + where ls = S.toList $ mps l + rs = S.toList $ mps r + -- The subsets on the left multipartition overlap the subsets of the + -- right multipartition. + lOverlaps = map (overlap r) ls + -- The subsets on the left multipartition overlap the subsets of the + -- right multipartition. + rOverlaps = map (overlap l) rs diff --git a/src/ELynx/Data/Tree/NamedTree.hs b/src/ELynx/Data/Tree/NamedTree.hs --- a/src/ELynx/Data/Tree/NamedTree.hs +++ b/src/ELynx/Data/Tree/NamedTree.hs @@ -16,8 +16,22 @@ ( Named (..) ) where -import qualified Data.ByteString.Lazy.Char8 as L +import qualified Data.ByteString.Lazy.Builder as L +import qualified Data.ByteString.Lazy.Char8 as L -- | Data types with names. class Named a where getName :: a -> L.ByteString + +instance Named Int where + getName = L.toLazyByteString . L.intDec + +instance Named Double where + getName = L.toLazyByteString . L.doubleDec + +instance Named Char where + getName = L.toLazyByteString . L.char8 + +instance Named L.ByteString where + getName = id + diff --git a/src/ELynx/Data/Tree/PhyloTree.hs b/src/ELynx/Data/Tree/PhyloTree.hs --- a/src/ELynx/Data/Tree/PhyloTree.hs +++ b/src/ELynx/Data/Tree/PhyloTree.hs @@ -12,78 +12,66 @@ Creation date: Thu Jan 17 16:08:54 2019. -Phylogenetic nodes have a branch length and a label. - -The easiest label type is 'Int': 'PhyloIntLabel'. - -Also, the 'L.ByteString' label is needed often: 'PhyloByteStringLabel'. - -XXX: This is all too complicated. Maybe I should just define a standard tree object like -> data PhyloTree a = Tree (PhyloLabel a) -and that's it. Forget about type classes like Measurable, Named and so on. +Phylogenetic nodes labels, aka 'PhyloLabel's, have a branch length and an +arbitrary label type, e.g., of type 'Int'. -} module ELynx.Data.Tree.PhyloTree ( PhyloLabel (..) - , PhyloIntLabel - , PhyloByteStringLabel - , removeBrLen + , removeBrInfo ) where -import qualified Data.ByteString.Lazy.Builder as L -import qualified Data.ByteString.Lazy.Char8 as L import Data.Function +import Data.Maybe (fromMaybe) import Data.Tree -import Test.QuickCheck +import Test.QuickCheck hiding (label) import ELynx.Data.Tree.BranchSupportTree import ELynx.Data.Tree.MeasurableTree import ELynx.Data.Tree.NamedTree -- | A primitive label type for phylogenetic trees with a name, possibly a --- branch support value, and a 'Double' branch length. -data PhyloLabel a = PhyloLabel { pLabel :: a - , pBrSup :: Maybe Double - , pBrLen :: Double } - deriving (Read, Show, Eq) +-- branch support value, and possibly a branch length. +data PhyloLabel a = PhyloLabel { label :: a + , brSup :: Maybe Double + , brLen :: Maybe Double } + deriving (Read, Show, Eq) instance Ord a => Ord (PhyloLabel a) where - compare = compare `on` pLabel + compare = compare `on` label +-- | If no branch length is available, 0 is returned. This is probably not the +-- best (and final) behavior. instance Measurable (PhyloLabel a) where - getLen = pBrLen - setLen l (PhyloLabel lbl s _) - | l >= 0 = PhyloLabel lbl s l - | otherwise = error "Branch lengths cannot be negative." + getLen = fromMaybe 0 . brLen + setLen l x + | l >= 0 = x {brLen = Just l} + | otherwise = error $ "Branch lengths cannot be negative: " <> show l -instance BranchSupportLabel (PhyloLabel a) where - getBranchSupport = pBrSup - setBranchSupport Nothing l = l {pBrSup = Nothing} +instance BranchSupported (PhyloLabel a) where + getBranchSupport = brSup + setBranchSupport Nothing l = l {brSup = Nothing} setBranchSupport (Just s) l - | s > 0 = l {pBrSup = Just s} + | s > 0 = l {brSup = Just s} | otherwise = error "Branch support cannot be negative." +-- Of course, the boundaries for branch support and length are chosen pretty +-- arbitrarily :). +-- +-- XXX: This instance does not produce values without branch lengths nor branch +-- supports. instance Arbitrary a => Arbitrary (PhyloLabel a) where arbitrary = PhyloLabel <$> arbitrary <*> (Just <$> choose (0, 100)) - <*> choose (0, 10) - --- | Tree node with 'Int' label. -type PhyloIntLabel = PhyloLabel Int - -instance Named PhyloIntLabel where - getName = L.toLazyByteString . L.intDec . pLabel - --- | Tree node with 'L.ByteString' label. Important for parsing --- 'ELynx.Import.Tree.Newick' files. -type PhyloByteStringLabel = PhyloLabel L.ByteString + <*> (Just <$> choose (0, 10) ) -instance Named PhyloByteStringLabel where - getName = pLabel +instance Named a => Named (PhyloLabel a) where + getName = getName . label --- | Remove branch lengths from tree. -removeBrLen :: Tree (PhyloLabel a) -> Tree a -removeBrLen = fmap pLabel +-- | Remove all branch relevant information from all nodes of the tree; only +-- retain the labels. +removeBrInfo :: Tree (PhyloLabel a) -> Tree a +removeBrInfo = fmap label diff --git a/src/ELynx/Data/Tree/Subset.hs b/src/ELynx/Data/Tree/Subset.hs new file mode 100644 --- /dev/null +++ b/src/ELynx/Data/Tree/Subset.hs @@ -0,0 +1,105 @@ +{-# LANGUAGE GeneralizedNewtypeDeriving #-} + +{- | +Module : ELynx.Data.Tree.Subset +Description : A subset is a set of elements +Copyright : (c) Dominik Schrempf 2019 +License : GPL-3 + +Maintainer : dominik.schrempf@gmail.com +Stability : unstable +Portability : portable + +Creation date: Fri Dec 13 11:02:43 2019. + +-} + +module ELynx.Data.Tree.Subset + ( Subset () + , sfromset + , sfromlist + , smap + , snull + , sempty + , ssingleton + , sunion + , sunions + , sdifference + , sintersection + , sdisjoint + , smember + , sshow + ) where + +import Data.List (intercalate) +import qualified Data.Set as S + +-- | A 'Subset' is a set of elements of type a. For example, on phylogenetic +-- trees, a 'Subset' is a set of leaves. In this case, a 'Subset' is induced, +-- for example, by a node on the (rooted) tree. The 'Subset's of leaf nodes are +-- singletons. The 'Subset's of the root node is the set of all leaves. +-- 'Subset's are the building blocks of partitions. Each branch on the tree +-- induces a bipartition, or a pair of 'Subset's, see +-- 'ELynx.Data.Tree.Bipartition'. Multifurcations induce multipartitions, see +-- 'ELynx.Data.Tree.Multipartition'. +-- +-- Internally, a subset is just an 'S.Set', since the order of elements +-- within the subset is not important, but the uniqueness of elements is. +newtype Subset a = SS {subset :: S.Set a} + deriving (Show, Read, Eq, Ord, Semigroup, Monoid, Foldable) + +-- | Create a subset from a set. +sfromset :: S.Set a -> Subset a +sfromset = SS + +-- | Create a subset from a list. Throws an error if duplicate elements are +-- present in the list. +sfromlist :: Ord a => [a] -> Subset a +sfromlist l = if S.size s == length l + then sfromset s + else error "pfromlist: List contains duplicate elements." + where s = S.fromList l + +-- | Map a function over all elements in a subset. +smap :: Ord b => (a -> b) -> Subset a -> Subset b +smap f = SS . S.map f . subset + +-- | Is the subset empty? +snull :: Subset a -> Bool +snull = S.null . subset + +-- | The empty subset. +sempty :: Subset a +sempty = SS S.empty + +-- | A subset with one element. +ssingleton :: a -> Subset a +ssingleton = SS . S.singleton + +-- | Unite two subsets. +sunion :: Ord a => Subset a -> Subset a -> Subset a +sunion p q = SS $ subset q `S.union` subset p + +-- | Unite a list of subsets. +sunions :: Ord a => [Subset a] -> Subset a +sunions = SS . S.unions . map subset + +-- | Difference of two subsets. +sdifference :: Ord a => Subset a -> Subset a -> Subset a +sdifference p q = SS $ subset p S.\\ subset q + +-- | Intersection of two subsets. +sintersection :: Ord a => Subset a -> Subset a -> Subset a +sintersection p q = SS $ S.intersection (subset p) (subset q) + +-- | Are two subsets disjoint? +sdisjoint :: Ord a => Subset a -> Subset a -> Bool +sdisjoint p q = S.disjoint (subset p) (subset q) + +-- | Check if an element is member of a subset. +smember :: Ord a => a -> Subset a -> Bool +smember x = S.member x . subset + +-- | Show the elements of a subset in a human readable way. +sshow :: (a -> String) -> Subset a -> String +sshow f = intercalate "," . map f . S.toList . subset diff --git a/src/ELynx/Data/Tree/Tree.hs b/src/ELynx/Data/Tree/Tree.hs --- a/src/ELynx/Data/Tree/Tree.hs +++ b/src/ELynx/Data/Tree/Tree.hs @@ -49,6 +49,12 @@ , pruneWith , merge , tZipWith + , partitionTree + , subForestGetSubsets + , bifurcating + , roots + , connect + , clades ) where import Control.Monad @@ -60,6 +66,7 @@ import Data.Tree import System.Random.MWC +import ELynx.Data.Tree.Subset import ELynx.Tools.Random -- | The simplest tree. Usually an extant leaf. @@ -80,7 +87,7 @@ -- rootNodesAgreeWith :: (Ord c) => (a -> c) -> Tree a -> (b -> c) -> Tree b -> Bool -- rootNodesAgreeWith f s g t = -- f (rootLabel s) == g (rootLabel t) && --- S.fromList sDs `S.isSubsetOf` S.fromList tDs +-- S.fromList sDs `S.isSubSetOf` S.fromList tDs -- where sDs = map (f . rootLabel) (subForest s) -- tDs = map (g . rootLabel) (subForest t) @@ -117,10 +124,10 @@ => Int -> Seq.Seq a -> Int -> Tree a -> Gen (PrimState m) -> m [Maybe (Tree a)] nSubSamples nS lvs nL tree g = replicateM nS $ subSample lvs nL tree g --- | Prune degree 2 inner nodes. The information stored in a pruned node can be --- used to change the daughter node. To discard this information, use, --- @pruneWith const tree@, otherwise @pruneWith (\daughter parent -> combined) --- tree@. +-- | Prune or remove degree 2 inner nodes. The information stored in a pruned +-- node can be used to change the daughter node. To discard this information, +-- use, @pruneWith const tree@, otherwise @pruneWith (\daughter parent -> +-- combined) tree@. pruneWith :: (a -> a -> a) -> Tree a -> Tree a pruneWith _ n@(Node _ []) = n pruneWith f (Node paLbl [ch]) = let lbl = f (rootLabel ch) paLbl @@ -141,3 +148,98 @@ tZipWith f xs = sequenceA . snd . mapAccumL pair xs where pair [] _ = ([], Nothing) pair (y:ys) z = (ys, Just (f y z)) + +-- | Each node of a tree is root of a subtree. Get the leaves of the subtree of +-- each node. +partitionTree :: (Ord a) => Tree a -> Tree (Subset a) +partitionTree (Node l []) = Node (ssingleton l) [] +partitionTree (Node _ xs) = Node (sunions $ map rootLabel xs') xs' + where xs' = map partitionTree xs + +-- | Loop through each tree in a forest to report the complementary leaf sets. +subForestGetSubsets :: (Ord a) + => Subset a -- ^ Complementary partition at the stem + -> Tree (Subset a) -- ^ Tree with partition nodes + -> [Subset a] +subForestGetSubsets lvs t = lvsOthers + where + xs = subForest t + nChildren = length xs + lvsChildren = map rootLabel xs + lvsOtherChildren = [ sunions $ lvs + : take i lvsChildren ++ drop (i+1) lvsChildren + | i <- [0 .. (nChildren - 1)] ] + lvsOthers = map (sunion lvs) lvsOtherChildren + +-- | Check if a tree is bifurcating and does not include degree two nodes. I +-- know, one should use a proper data structure to encode bifurcating trees, but +-- I don't have enough time for this now. +bifurcating :: Tree a -> Bool +bifurcating (Node _ [] ) = True +bifurcating (Node _ [_] ) = False +bifurcating (Node _ [x, y]) = bifurcating x && bifurcating y +bifurcating (Node _ _ ) = False + +-- TODO: This bifurcating stuff irritates me. There are two solutions: +-- +-- 1. Use a bifurcating data structure. +-- +-- 2. Do not move down multifurcations. +-- +-- Solution 1 is pretty, but doesn't allow for what we actually want to do, +-- which is solution 2. We want to encode clades that for sure do not contain +-- the root as multifurcations. + +-- | For a rooted, bifurcating tree, get all possible rooted trees. For a tree +-- with @n>2@ leaves, there are @(2n-3)@ rooted trees. Beware, a bifurcating +-- tree without degree two nodes is assumed (see 'bifurcating'). The root node +-- is moved. +roots :: Tree a -> [Tree a] +-- Leaves, and cherries have to be handled separately, because they cannot be +-- rotated. +roots t@(Node _ []) = [t] +roots t@(Node _ [Node _ [], Node _ []]) = [t] +roots t | bifurcating t = t : left t ++ right t + | otherwise = error "roots: Tree is not bifurcating." + +-- Move the root to the left. +left :: Tree a -> [Tree a] +left (Node i [Node j [x] , z]) = let t' = Node i [x , Node j [z] ] + in t' : left t' +left (Node i [Node j [x, y], z]) = let tll = Node i [x , Node j [y,z]] + tlr = Node i [Node j [x,z], y ] + in tll : tlr : left tll ++ right tlr +left (Node _ [Node _ [] , _]) = [] +left (Node _ [] ) = error "left: Encountered a leaf." +left _ = error "left: Tree is not bifurcating." + +-- Move the root to the right. +right :: Tree a -> [Tree a] +right (Node i [x, Node j [z] ]) = let t' = Node i [Node j [x] , z ] + in t' : right t' +right (Node i [x, Node j [y, z]]) = let trl = Node i [y , Node j [x,z]] + trr = Node i [Node j [x,y], z ] + in trl : trr : left trl ++ right trr +right (Node _ [_, Node _ [] ]) = [] +right (Node _ [] ) = error "right: Encountered a leaf." +right (Node _ [_]) = error "right: TODO; this case has to be handled separately." +right _ = error "left: Tree is not bifurcating." + +-- | Connect two trees in all possible ways. +-- +-- Basically, introduce a branch between two trees. If the trees have n, and m +-- branches, respectively, there are n*m ways to connect them. +-- +-- A base node has to be given which will be used wherever the new node is +-- introduced. +connect :: a -> Tree a -> Tree a -> [Tree a] +connect n l r = [ Node n [x, y] | x <- roots l, y <- roots r] + +-- | Get clades induced by multifurcations. +-- +-- XXX: Probably introduce a new module defining a Clade. +clades :: Ord a => Tree a -> [Subset a] +clades (Node _ [] ) = [] +clades (Node _ [x]) = clades x +clades (Node _ [x, y]) = clades x ++ clades y +clades t = sfromlist (leaves t) : concatMap clades (subForest t) diff --git a/src/ELynx/Export/Tree/Newick.hs b/src/ELynx/Export/Tree/Newick.hs --- a/src/ELynx/Export/Tree/Newick.hs +++ b/src/ELynx/Export/Tree/Newick.hs @@ -22,41 +22,36 @@ -- , toNewickPhyloByteStringTree ) where -import qualified Data.ByteString.Lazy.Builder as L -import qualified Data.ByteString.Lazy.Char8 as L -import Data.List (intersperse) +import qualified Data.ByteString.Lazy.Builder as L +import qualified Data.ByteString.Lazy.Char8 as L +import Data.List (intersperse) -- import Data.Maybe import Data.Tree -import ELynx.Data.Tree.BranchSupportTree -import ELynx.Data.Tree.MeasurableTree +-- import ELynx.Data.Tree.BranchSupportTree +-- import ELynx.Data.Tree.MeasurableTree import ELynx.Data.Tree.NamedTree --- import ELynx.Data.Tree.PhyloTree -import ELynx.Tools.ByteString (c2w) +import ELynx.Data.Tree.PhyloTree +import ELynx.Tools.ByteString (c2w) -- | General conversion of a tree into a Newick 'L.Bytestring'. Use provided -- functions to extract node labels and branch lengths builder objects. See also -- Biobase.Newick.Export. -toNewick :: (Named a, Measurable a, BranchSupportLabel a) => Tree a -> L.ByteString +toNewick :: Named a => Tree (PhyloLabel a) -> L.ByteString toNewick t = L.toLazyByteString $ go t <> L.word8 (c2w ';') where go (Node l []) = lbl l go (Node l ts) = L.word8 (c2w '(') <> mconcat (intersperse (L.word8 $ c2w ',') $ map go ts) - <> L.word8 (c2w ')') <> brSup l + <> L.word8 (c2w ')') <> lbl l - lbl l = L.lazyByteString (getName l) - <> L.word8 (c2w ':') - <> L.doubleDec (getLen l) - brSup l = maybe mempty L.doubleDec (getBranchSupport l) - --- -- | Convenience function for exporting trees with 'Int' labels and 'Double' --- -- branch lengths. --- toNewickPhyloIntTree :: Tree PhyloIntLabel -> L.ByteString --- toNewickPhyloIntTree = toNewickWith (L.intDec . pLabel) (L.doubleDec . pBrLen) - --- -- | Convenience function for exporting trees with 'L.ByteString' labels and --- -- 'Double' branch lengths. --- toNewickPhyloByteStringTree :: Tree PhyloByteStringLabel -> L.ByteString --- toNewickPhyloByteStringTree = toNewickWith (L.lazyByteString . pLabel) (L.doubleDec . pBrLen) + brSupStr bs = L.word8 (c2w '[') <> L.doubleDec bs <> L.word8 (c2w ']') + mBrSup l = maybe mempty brSupStr (brSup l) + brLenStr bl = L.word8 (c2w ':') <> L.doubleDec bl + mBrLen l = maybe mempty brLenStr (brLen l) + lbl l = L.lazyByteString (getName l) + <> mBrLen l + -- After reading several discussion, I go for the "more semantical + -- form" with branch support values in square brackets. + <> mBrSup l diff --git a/src/ELynx/Import/Tree/Newick.hs b/src/ELynx/Import/Tree/Newick.hs --- a/src/ELynx/Import/Tree/Newick.hs +++ b/src/ELynx/Import/Tree/Newick.hs @@ -21,13 +21,19 @@ module ELynx.Import.Tree.Newick ( Parser + -- * Newick tree format , newick + , oneNewick , manyNewick , forest , leaf , node , name , branchLength + -- * Newick tree format with branch support as node names (e.g., used by IQ-TREE) + , newickIqTree + , oneNewickIqTree + , manyNewickIqTree ) where import qualified Data.ByteString.Lazy as L @@ -44,28 +50,30 @@ -- | Shortcut. type Parser = Parsec Void L.ByteString --- | Parse many Newick trees. -manyNewick :: Parser [Tree PhyloByteStringLabel] -manyNewick = some (newick <* space) <* eof "manyNewick" - --- | Parse a Newick tree. -newick :: Parser (Tree PhyloByteStringLabel) +-- | Parse a single Newick tree. Also succeeds when more trees follow. +newick :: Parser (Tree (PhyloLabel L.ByteString)) newick = tree <* char (c2w ';') "newick" -tree :: Parser (Tree PhyloByteStringLabel) +-- | Parse a single Newick tree. Fails when end of file is not reached. +oneNewick :: Parser (Tree (PhyloLabel L.ByteString)) +oneNewick = newick <* space <* eof "oneNewick" + +-- | Parse many Newick trees until end of file. +manyNewick :: Parser [Tree (PhyloLabel L.ByteString)] +manyNewick = some (newick <* space) <* eof "manyNewick" + +tree :: Parser (Tree (PhyloLabel L.ByteString)) tree = space *> (branched <|> leaf) "tree" -branched :: Parser (Tree PhyloByteStringLabel) +branched :: Parser (Tree (PhyloLabel L.ByteString)) branched = do f <- forest - s <- branchSupport n <- node "branched" - let n' = n {pBrSup = s} - return $ Node n' f + return $ Node n f -- | A 'forest' is a set of trees separated by @,@ and enclosed by parentheses. -forest :: Parser [Tree PhyloByteStringLabel] +forest :: Parser [Tree (PhyloLabel L.ByteString)] forest = do _ <- char (c2w '(') f <- tree `sepBy1` char (c2w ',') @@ -74,22 +82,27 @@ return f branchSupport :: Parser (Maybe Double) -branchSupport = optional $ try float <|> try decimalAsDouble +branchSupport = optional $ do + _ <- try $ char (c2w '[') + s <- try float <|> try decimalAsDouble + _ <- try $ char (c2w ']') + return s -- | A 'leaf' is a 'node' without children. -leaf :: Parser (Tree PhyloByteStringLabel) +leaf :: Parser (Tree (PhyloLabel L.ByteString)) leaf = do n <- node "leaf" return $ Node n [] -- | A 'node' has a name and a 'branchLength'. -node :: Parser PhyloByteStringLabel +node :: Parser (PhyloLabel L.ByteString) node = do n <- name b <- branchLength + s <- branchSupport "node" - return $ PhyloLabel n Nothing b + return $ PhyloLabel n s b checkNameCharacter :: Word8 -> Bool checkNameCharacter c = c `notElem` map c2w " :;()[]," @@ -99,12 +112,62 @@ name :: Parser L.ByteString name = L.pack <$> many (satisfy checkNameCharacter) "name" --- | Branch lengths default to 0. -branchLength :: Parser Double -branchLength = char (c2w ':') *> branchLengthGiven <|> pure 0 "branchLength" +-- | Branch length. +branchLength :: Parser (Maybe Double) +branchLength = (optional $ char (c2w ':') *> branchLengthGiven) "branchLength" branchLengthGiven :: Parser Double branchLengthGiven = try float <|> decimalAsDouble decimalAsDouble :: Parser Double decimalAsDouble = fromIntegral <$> (decimal :: Parser Int) + +-------------------------------------------------------------------------------- +-- IQ-TREE STUFF. + +-- | IQ-TREE stores the branch support as node names after the closing bracket of a forest. +newickIqTree :: Parser (Tree (PhyloLabel L.ByteString)) +newickIqTree = treeIqTree <* char (c2w ';') "newickIqTree" + +-- | IQ-TREE stores the branch support as node names after the closing bracket of a forest. +oneNewickIqTree :: Parser (Tree (PhyloLabel L.ByteString)) +oneNewickIqTree = newickIqTree <* space <* eof "oneNewickIqTree" + +-- | IQ-TREE stores the branch support as node names after the closing bracket of a forest. +manyNewickIqTree :: Parser [Tree (PhyloLabel L.ByteString)] +manyNewickIqTree = some (newickIqTree <* space) <* eof "manyNewickIqTree" + +-- IQ-TREE stores the branch support as node names after the closing bracket of a forest. +treeIqTree :: Parser (Tree (PhyloLabel L.ByteString)) +treeIqTree = space *> (branchedIqTree <|> leaf) "treeIqTree" + +-- IQ-TREE stores the branch support as node names after the closing bracket of a forest. +forestIqTree :: Parser [Tree (PhyloLabel L.ByteString)] +forestIqTree = do + _ <- char (c2w '(') + f <- treeIqTree `sepBy1` char (c2w ',') + _ <- char (c2w ')') + "forestIqTree" + return f + +-- IQ-TREE stores the branch support as node names after the closing bracket of a forest. +branchedIqTree :: Parser (Tree (PhyloLabel L.ByteString)) +branchedIqTree = do + f <- forestIqTree + s <- branchSupportIqTree + n <- nodeIqTree + "branchedIqTree" + let n' = n {brSup = s} + return $ Node n' f + +-- IQ-TREE stores the branch support as node names after the closing bracket of a forest. +branchSupportIqTree :: Parser (Maybe Double) +branchSupportIqTree = optional $ try float <|> try decimalAsDouble + +-- IQ-TREE stores the branch support as node names after the closing bracket of a forest. +nodeIqTree :: Parser (PhyloLabel L.ByteString) +nodeIqTree = do + n <- name + b <- branchLength + "nodeIqTree" + return $ PhyloLabel n Nothing b diff --git a/src/ELynx/Simulate/Coalescent.hs b/src/ELynx/Simulate/Coalescent.hs --- a/src/ELynx/Simulate/Coalescent.hs +++ b/src/ELynx/Simulate/Coalescent.hs @@ -10,6 +10,8 @@ Creation date: Wed May 16 13:13:11 2018. +TODO: This is unused, is it? + -} @@ -32,16 +34,16 @@ simulate :: (PrimMonad m) => Int -- ^ Number of leaves. -> Gen (PrimState m) - -> m (Tree PhyloIntLabel) + -> m (Tree (PhyloLabel Int)) simulate n = simulate' n 0 trs - where trs = [ singleton (PhyloLabel i Nothing 0.0) | i <- [0..n-1] ] + where trs = [ singleton (PhyloLabel i Nothing Nothing) | i <- [0..n-1] ] simulate' :: (PrimMonad m) => Int -> Int - -> [Tree PhyloIntLabel] + -> [Tree (PhyloLabel Int)] -> Gen (PrimState m) - -> m (Tree PhyloIntLabel) + -> m (Tree (PhyloLabel Int)) simulate' n a trs g | n <= 0 = error "Cannot construct trees without leaves." | n == 1 && length trs /= 1 = error "Too many trees provided." @@ -52,11 +54,11 @@ i <- uniformR (1, n-1) g -- The time of the coalescent event. t <- genContVar (coalescentDistributionCont n) g - let trs' = map (lengthenRoot t) trs -- Move time 't' up on the tree. + let trs' = map (lengthenStem t) trs -- Move time 't' up on the tree. tl = trs' !! (i-1) tr = trs' !! i -- Join the two chosen trees. - tm = Node (PhyloLabel a Nothing 0.0) [tl, tr] + tm = Node (PhyloLabel a Nothing Nothing) [tl, tr] -- Take the trees on the left, the merged tree, and the trees on the right. trs'' = take (i-1) trs' ++ [tm] ++ drop (i+1) trs' simulate' (n-1) a trs'' g diff --git a/src/ELynx/Simulate/PointProcess.hs b/src/ELynx/Simulate/PointProcess.hs --- a/src/ELynx/Simulate/PointProcess.hs +++ b/src/ELynx/Simulate/PointProcess.hs @@ -163,7 +163,7 @@ -> Rate -- ^ Birth rate -> Rate -- ^ Death rate -> Gen (PrimState m) -- ^ Generator (see 'System.Random.MWC') - -> m [Tree PhyloIntLabel] + -> m [Tree (PhyloLabel Int)] simulateNReconstructedTrees nT nP t l m g | nT <= 0 = return [] | otherwise = replicateM nT $ simulateReconstructedTree nP t l m g @@ -178,8 +178,8 @@ -> Rate -- ^ Birth rate -> Rate -- ^ Death rate -> Gen (PrimState m) -- ^ Generator (see 'System.Random.MWC') - -> m (Tree PhyloIntLabel) -simulateReconstructedTree n t l m g = toReconstructedTree <$> simulate n t l m g + -> m (Tree (PhyloLabel Int)) +simulateReconstructedTree n t l m g = toReconstructedTree 0 <$> simulate n t l m g -- | Convert a point process to a reconstructed tree. See Lemma 2.2. @@ -191,29 +191,31 @@ -- contain extinct leaves. I wanted to use a Monoid constraint to get the unit -- element, but this fails for classical 'Int's. So, I rather have another -- (useless) argument. -toReconstructedTree :: PointProcess Int Double - -> Tree PhyloIntLabel -toReconstructedTree pp@(PointProcess ps vs o) +toReconstructedTree :: a -- Default node label. + -> PointProcess a Double + -> Tree (PhyloLabel a) +toReconstructedTree l pp@(PointProcess ps vs o) | length ps /= length vs + 1 = error "Too few or too many points." | length vs <= 1 = error "Too few values." -- -- XXX: Test is deactivated. -- -- | otherwise = if isReconstructed treeOrigin then treeOrigin else error "Error in algorithm." | otherwise = treeOrigin where (vsSorted, isSorted) = sort pp - !lvs = [ singleton (PhyloLabel p Nothing 0) | p <- ps ] + !lvs = [ singleton (PhyloLabel p Nothing Nothing) | p <- ps ] !heights = replicate (length ps) 0 - !treeRoot = toReconstructedTree' isSorted vsSorted lvs heights + !treeRoot = toReconstructedTree' isSorted vsSorted l lvs heights !h = last vsSorted - !treeOrigin = lengthenRoot (o-h) treeRoot + !treeOrigin = lengthenStem (o-h) treeRoot -- Move up the tree, connect nodes when they join according to the point process. -toReconstructedTree' :: [Int] -- Sorted indices, see 'sort'. - -> [Double] -- Sorted merge values. - -> [Tree PhyloIntLabel] -- Leaves with accumulated root branch lengths. - -> [Double] -- Accumulated heights of the leaves. - -> Tree PhyloIntLabel -toReconstructedTree' [] [] trs _ = head trs -toReconstructedTree' is vs trs hs = toReconstructedTree' is' vs' trs'' hs' +toReconstructedTree' :: [Int] -- Sorted indices, see 'sort'. + -> [Double] -- Sorted merge values. + -> a -- Default node label. + -> [Tree (PhyloLabel a)] -- Leaves with accumulated root branch lengths. + -> [Double] -- Accumulated heights of the leaves. + -> Tree (PhyloLabel a) +toReconstructedTree' [] [] _ trs _ = head trs +toReconstructedTree' is vs l trs hs = toReconstructedTree' is' vs' l trs'' hs' -- For the algorithm, see 'ELynx.Coalescent.simulate', but index starts -- at zero. where !i = head is @@ -225,9 +227,9 @@ !hr = hs !! (i+1) !dvl = v - hl !dvr = v - hr - !tl = lengthenRoot dvl $ trs !! i - !tr = lengthenRoot dvr $ trs !! (i+1) + !tl = lengthenStem dvl $ trs !! i + !tr = lengthenStem dvr $ trs !! (i+1) !h' = hl + dvl -- Should be the same as 'hr + dvr'. - !tm = Node (PhyloLabel 0 Nothing 0) [tl, tr] + !tm = Node (PhyloLabel l Nothing Nothing) [tl, tr] !trs'' = take i trs ++ [tm] ++ drop (i+2) trs !hs' = take i hs ++ [h'] ++ drop (i+2) hs diff --git a/test/ELynx/Data/Tree/BipartitionSpec.hs b/test/ELynx/Data/Tree/BipartitionSpec.hs --- a/test/ELynx/Data/Tree/BipartitionSpec.hs +++ b/test/ELynx/Data/Tree/BipartitionSpec.hs @@ -17,57 +17,62 @@ (spec ) where -import qualified Data.ByteString.Lazy as L -import qualified Data.Map as M +import qualified Data.ByteString.Lazy as L +import qualified Data.Map as M import Data.Monoid -import qualified Data.Set as S +import qualified Data.Set as S import Data.Tree import Test.Hspec import ELynx.Data.Tree.Bipartition +import ELynx.Data.Tree.MeasurableTree import ELynx.Data.Tree.PhyloTree +import ELynx.Data.Tree.Subset import ELynx.Import.Tree.Newick import ELynx.Tools.InputOutput +sfrom :: [L.ByteString] -> Subset L.ByteString +sfrom = sfromset . S.fromList + treeFileSimple :: FilePath treeFileSimple = "data/TreeDist.trees" -getSimpleTrees :: IO [Tree PhyloByteStringLabel] +getSimpleTrees :: IO [Tree (PhyloLabel L.ByteString)] getSimpleTrees = parseFileWith manyNewick treeFileSimple bipartitionToBranchAnswer :: M.Map (Bipartition L.ByteString) (Sum Double) bipartitionToBranchAnswer = - M.fromList [ (bp (S.fromList ["B"]) (S.fromList ["A","C","D","E"]), Sum {getSum = 0.3}) - , (bp (S.fromList ["B","C","D","E"]) (S.fromList ["A"]), Sum {getSum = 0.1}) - , (bp (S.fromList ["B","C","E"]) (S.fromList ["A","D"]), Sum {getSum = 5.0e-2}) - , (bp (S.fromList ["B","E"]) (S.fromList ["A","C","D"]), Sum {getSum = 0.4}) - , (bp (S.fromList ["C"]) (S.fromList ["A","B","D","E"]), Sum {getSum = 1.0e-2}) - , (bp (S.fromList ["D"]) (S.fromList ["A","B","C","E"]), Sum {getSum = 0.25}) - , (bp (S.fromList ["E"]) (S.fromList ["A","B","C","D"]), Sum {getSum = 0.8}) ] + M.fromList [ (bp (sfrom ["B"]) (sfrom ["A","C","D","E"]), Sum {getSum = 0.3}) + , (bp (sfrom ["B","C","D","E"]) (sfrom ["A"]), Sum {getSum = 0.1}) + , (bp (sfrom ["B","C","E"]) (sfrom ["A","D"]), Sum {getSum = 5.0e-2}) + , (bp (sfrom ["B","E"]) (sfrom ["A","C","D"]), Sum {getSum = 0.4}) + , (bp (sfrom ["C"]) (sfrom ["A","B","D","E"]), Sum {getSum = 1.0e-2}) + , (bp (sfrom ["D"]) (sfrom ["A","B","C","E"]), Sum {getSum = 0.25}) + , (bp (sfrom ["E"]) (sfrom ["A","B","C","D"]), Sum {getSum = 0.8}) ] bipartitionsFirstTree :: S.Set (Bipartition L.ByteString) -bipartitionsFirstTree = S.fromList [ bp (S.fromList ["B"]) (S.fromList ["A","C","D","E"]) - , bp (S.fromList ["B","C","D","E"]) (S.fromList ["A"]) - , bp (S.fromList ["B","D","E"]) (S.fromList ["A","C"]) - , bp (S.fromList ["B","E"]) (S.fromList ["A","C","D"]) - , bp (S.fromList ["C"]) (S.fromList ["A","B","D","E"]) - , bp (S.fromList ["D"]) (S.fromList ["A","B","C","E"]) - , bp (S.fromList ["E"]) (S.fromList ["A","B","C","D"]) ] +bipartitionsFirstTree = S.fromList [ bp (sfrom ["B"]) (sfrom ["A","C","D","E"]) + , bp (sfrom ["B","C","D","E"]) (sfrom ["A"]) + , bp (sfrom ["B","D","E"]) (sfrom ["A","C"]) + , bp (sfrom ["B","E"]) (sfrom ["A","C","D"]) + , bp (sfrom ["C"]) (sfrom ["A","B","D","E"]) + , bp (sfrom ["D"]) (sfrom ["A","B","C","E"]) + , bp (sfrom ["E"]) (sfrom ["A","B","C","D"]) ] bipartitionsSecondTree :: S.Set (Bipartition L.ByteString) -bipartitionsSecondTree = S.fromList [ bp (S.fromList ["B"]) (S.fromList ["A","C","D","E"]) - , bp (S.fromList ["B","C","D","E"]) (S.fromList ["A"]) - , bp (S.fromList ["B","C","E"]) (S.fromList ["A","D"]) - , bp (S.fromList ["B","E"]) (S.fromList ["A","C","D"]) - , bp (S.fromList ["C"]) (S.fromList ["A","B","D","E"]) - , bp (S.fromList ["D"]) (S.fromList ["A","B","C","E"]) - , bp (S.fromList ["E"]) (S.fromList ["A","B","C","D"])] +bipartitionsSecondTree = S.fromList [ bp (sfrom ["B"]) (sfrom ["A","C","D","E"]) + , bp (sfrom ["B","C","D","E"]) (sfrom ["A"]) + , bp (sfrom ["B","C","E"]) (sfrom ["A","D"]) + , bp (sfrom ["B","E"]) (sfrom ["A","C","D"]) + , bp (sfrom ["C"]) (sfrom ["A","B","D","E"]) + , bp (sfrom ["D"]) (sfrom ["A","B","C","E"]) + , bp (sfrom ["E"]) (sfrom ["A","B","C","D"])] spec :: Spec spec = do describe "bipartitions" $ it "calculates correct bipartitions for sample trees" $ do - simpleTrees <- map removeBrLen <$> getSimpleTrees + simpleTrees <- map removeBrInfo <$> getSimpleTrees let t1 = head simpleTrees t2 = simpleTrees !! 1 bipartitions t1 `shouldBe` bipartitionsFirstTree @@ -75,5 +80,5 @@ describe "bipartitionToBranch" $ it "creates a map from bipartitions to branch lengths" $ do simpleTrees <- getSimpleTrees - bipartitionToBranch pLabel (Sum . pBrLen) (simpleTrees !! 2) + bipartitionToBranchLength label (Sum . getLen) (simpleTrees !! 2) `shouldBe` bipartitionToBranchAnswer diff --git a/test/ELynx/Data/Tree/DistanceSpec.hs b/test/ELynx/Data/Tree/DistanceSpec.hs --- a/test/ELynx/Data/Tree/DistanceSpec.hs +++ b/test/ELynx/Data/Tree/DistanceSpec.hs @@ -1,4 +1,5 @@ {-# LANGUAGE FlexibleInstances #-} + {- | Module : ELynx.Data.Tree.DistanceSpec Copyright : (c) Dominik Schrempf 2019 @@ -16,6 +17,7 @@ (spec ) where +import qualified Data.ByteString.Lazy.Char8 as L import Data.Tree import Test.Hspec import Test.QuickCheck @@ -31,25 +33,25 @@ treeFileSimple :: FilePath treeFileSimple = "data/TreeDist.trees" -getSimpleTrees :: IO [Tree PhyloByteStringLabel] +getSimpleTrees :: IO [Tree (PhyloLabel L.ByteString)] getSimpleTrees = parseFileWith manyNewick treeFileSimple treeFileMany :: FilePath treeFileMany = "data/Many.trees" -getManyTrees :: IO [Tree PhyloByteStringLabel] +getManyTrees :: IO [Tree (PhyloLabel L.ByteString)] getManyTrees = parseFileWith manyNewick treeFileMany -- I used treedist from Phylip to get the correct results. -- See http://evolution.genetics.washington.edu/phylip/doc/treedist.html. -symmetricDistanceAnswers :: [Int] -symmetricDistanceAnswers = +symmetricAnswers :: [Int] +symmetricAnswers = [ 6, 8, 0, 0, 12, 20, 18, 20, 10, 2, 10, 4, 4, 4, 4, 4, 10, 16, 8, 2, 4, 0, 0, 0, 10, 4, 0, 0, 2, 2, 0, 0, 4, 0, 2, 0, 8, 6, 2, 6, 4, 4, 8, 0, 0, 4, 2, 0, 10, 0, 0, 10 ] -branchScoreDistanceAnswers :: [Double] -branchScoreDistanceAnswers = +branchScoreAnswers :: [Double] +branchScoreAnswers = [ 8.567916e-02, 9.570577e-02, 1.704571e-02, 7.603990e-03, 6.149761e-01, 3.557070e-01, 2.329811e-01, 3.820208e-01, 1.895421e-02, 6.302364e-03, 2.083286e-02, 1.023777e-03, 2.138244e-02, 1.444380e-02, 1.958628e-02, @@ -65,35 +67,53 @@ prop_dist_same_tree :: (Num b, Eq b) => (Tree a -> Tree a -> b) -> Tree a -> Bool prop_dist_same_tree distanceMeasure t = distanceMeasure t t == 0 --- TODO: Microsporidia trees with branch support values. --- getMicrospoPoissonTree :: IO (Tree PhyloByteStringLabel) --- getMicrospoPoissonTree = parseFileWith newick "data/MicrospoPoisson.tree" --- getMicrospoUDM32Tree = parseFileWith newick "data/MicrospoEDM32.tree" --- getMicrospoUDM64Tree = parseFileWith newick "data/MicrospoEDM64.tree" - each :: Int -> [a] -> [a] each n = map head . takeWhile (not . null) . iterate (drop n) +multifurcating :: Tree Char +multifurcating = Node ' ' [Node 'A' [], Node 'B' [], Node ' ' [Node 'C' [], Node 'D' [], Node 'E' []]] + +bifurcatingComp :: Tree Char +bifurcatingComp = Node ' ' [ Node ' ' [Node 'A' [], Node 'B'[]] + , Node ' ' [Node 'C' [], Node ' ' [Node 'D' [], Node 'E' []]]] + +bifurcatingIncomp :: Tree Char +bifurcatingIncomp = Node ' ' [ Node ' ' [Node 'A' [], Node 'C'[]] + , Node ' ' [Node 'B' [], Node ' ' [Node 'D' [], Node 'E' []]]] + spec :: Spec spec = do - describe "symmetricDistance" $ + describe "symmetric" $ do it "calculates correct distances for sample trees" $ do simpleTrees <- getSimpleTrees - symmetricDistance (head simpleTrees) (simpleTrees !! 1) `shouldBe` 2 + symmetric (head simpleTrees) (simpleTrees !! 1) `shouldBe` 2 manyTrees <- getManyTrees -- Since treedist computes the distance between adjacent pairs, in the -- following manner: [tr0, tr1, tr2, tr3] -> [dist tr0 tr1, dist tr2 tr3], -- we have to skip some distances. - each 2 (computeAdjacentDistances (symmetricDistanceWith getName) manyTrees) - `shouldBe` symmetricDistanceAnswers + each 2 (adjacent (symmetricWith getName) manyTrees) + `shouldBe` symmetricAnswers + it "is zero for a collection of random trees" $ + property $ prop_dist_same_tree + (symmetric :: Tree (PhyloLabel Double) -> Tree (PhyloLabel Double) -> Int) - describe "incompatibleSplitDistance" $ - it "calculates correct distances for completely collapsed trees" $ - property $ prop_dist_same_tree (incompatibleSplitsDistance :: Tree PhyloIntLabel -> Tree PhyloIntLabel -> Int) + describe "incompatibleSplit" $ do + it "calculates correct distances for sample trees" $ do + incompatibleSplits multifurcating bifurcatingComp `shouldBe` 0 + incompatibleSplits bifurcatingComp multifurcating `shouldBe` 0 + incompatibleSplits bifurcatingIncomp multifurcating `shouldBe` 2 + incompatibleSplits multifurcating bifurcatingIncomp `shouldBe` 2 + it "is zero for a collection of random trees" $ + property $ prop_dist_same_tree + (incompatibleSplits :: Tree (PhyloLabel Double) -> Tree (PhyloLabel Double) -> Int) - describe "branchScoreDistance" $ + + describe "branchScore" $ do it "calculates correct distances for sample trees" $ do manyTrees <- getManyTrees - print branchScoreDistanceAnswers - each 2 (computeAdjacentDistances branchScoreDistance manyTrees) - `shouldSatisfy` nearlyEqListWith 1e-5 branchScoreDistanceAnswers + -- print branchScoreAnswers + each 2 (adjacent branchScore manyTrees) + `shouldSatisfy` nearlyEqListWith 1e-5 branchScoreAnswers + it "is zero for a collection of random trees" $ + property $ prop_dist_same_tree + (branchScore :: Tree (PhyloLabel Double) -> Tree (PhyloLabel Double) -> Double) diff --git a/test/ELynx/Data/Tree/MultipartitionSpec.hs b/test/ELynx/Data/Tree/MultipartitionSpec.hs new file mode 100644 --- /dev/null +++ b/test/ELynx/Data/Tree/MultipartitionSpec.hs @@ -0,0 +1,59 @@ +{-# LANGUAGE OverloadedStrings #-} + +{- | +Module : ELynx.Data.Tree.MultipartitionSpec +Copyright : (c) Dominik Schrempf 2019 +License : GPL-3 + +Maintainer : dominik.schrempf@gmail.com +Stability : unstable +Portability : portable + +Creation date: Fri Aug 30 09:38:50 2019. + +-} + +module ELynx.Data.Tree.MultipartitionSpec + (spec + ) where + +import Data.Set (Set, empty, fromList) +import Data.Tree +import Test.Hspec +import Test.QuickCheck +import Test.QuickCheck.Instances.Containers () + +import ELynx.Data.Tree.BranchSupportTree +import ELynx.Data.Tree.MeasurableTree +import ELynx.Data.Tree.Multipartition (Multipartition, mp, + multipartitions) +import ELynx.Data.Tree.NamedTree +import ELynx.Data.Tree.PhyloTree +import ELynx.Data.Tree.Subset (sfromlist) + +ex1 :: Tree Int +ex1 = Node 0 [Node 1 [], Node 2 [Node 4 [], Node 5 [], Node 6 []], Node 3 []] + +sol1 :: Set (Multipartition Int) +sol1 = fromList [ mp [sfromlist [1], sfromlist [3], sfromlist [4,5,6]] + , mp [sfromlist [1,3], sfromlist [4], sfromlist [5], sfromlist [6]] ] + +ex2 :: Tree Int +ex2 = Node 0 [Node 1 [], Node 2 [], Node 0 [Node 3 [], Node 4 []], Node 5 []] + +sol2 :: Set (Multipartition Int) +sol2 = fromList [ mp [sfromlist [1], sfromlist [2], sfromlist [3,4], sfromlist [5]] ] + +prop_bifurcating_tree :: (Ord a, Measurable a, Named a, BranchSupported a) => Tree a -> Bool +prop_bifurcating_tree t = multipartitions (removeMultifurcations t) == empty + +spec :: Spec +spec = + describe "multipartitions" $ do + it "calculates correct multipartitions for a sample tree" $ do + -- t <- removeBrInfo <$> getMultifurcatingTree + -- print $ multipartitions ex1 + multipartitions ex1 `shouldBe` sol1 + multipartitions ex2 `shouldBe` sol2 + it "is empty for a collection of random bifurcating trees" $ + property (prop_bifurcating_tree :: Tree (PhyloLabel Double) -> Bool) diff --git a/test/ELynx/Data/Tree/TreeSpec.hs b/test/ELynx/Data/Tree/TreeSpec.hs --- a/test/ELynx/Data/Tree/TreeSpec.hs +++ b/test/ELynx/Data/Tree/TreeSpec.hs @@ -1,3 +1,5 @@ +{-# LANGUAGE OverloadedStrings #-} + {- | Module : ELynx.Data.Tree.TreeSpec Copyright : (c) Dominik Schrempf 2019 @@ -15,16 +17,29 @@ (spec ) where -import qualified Data.ByteString.Lazy.Char8 as L +import qualified Data.ByteString.Lazy.Char8 as L import Data.Maybe import Data.Tree -import Test.Hspec +import Test.Hspec (Spec, describe, it, + shouldBe) +import Test.Hspec.QuickCheck (modifyMaxSize) +import Test.QuickCheck hiding (label) +import Test.QuickCheck.Instances.Containers () -import ELynx.Data.Tree.MeasurableTree -import ELynx.Data.Tree.PhyloTree -import ELynx.Data.Tree.Tree -import ELynx.Import.Tree.Newick hiding (node) -import ELynx.Tools.InputOutput (parseByteStringWith) +import ELynx.Data.Tree.Bipartition (bipartition, compatible) +import ELynx.Data.Tree.MeasurableTree (height, prune) +import ELynx.Data.Tree.NamedTree (getName) +import ELynx.Data.Tree.PhyloTree (PhyloLabel (PhyloLabel), + label) +import ELynx.Data.Tree.Subset (Subset, smap) +import ELynx.Data.Tree.Tree (bifurcating, clades, + connect, leaves, + pruneWith, roots, + subTree) +import ELynx.Import.Tree.Newick (manyNewick, newick, + oneNewick) +import ELynx.Tools.InputOutput (parseByteStringWith, + parseFileWith) node :: Int -> Tree Int node n = Node n [] @@ -41,13 +56,36 @@ sampleTreeBS :: L.ByteString sampleTreeBS = L.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);" - -largeTree :: Tree PhyloByteStringLabel +largeTree :: Tree (PhyloLabel L.ByteString) largeTree = parseByteStringWith "Sample newick byte string" newick sampleTreeBS -subSampleLargeTree :: Tree PhyloByteStringLabel -subSampleLargeTree = fromJust $ subTree ((== 'P') . L.head . pLabel) largeTree +subSampleLargeTree :: Tree (PhyloLabel L.ByteString) +subSampleLargeTree = fromJust $ subTree ((== 'P') . L.head . label) largeTree +prop_roots :: Tree a -> Bool +prop_roots t + -- XXX: Skip not bifurcating trees. This is ugly, I know. + | not $ bifurcating t = True + | length (leaves t) < 3 = length (roots t) == 1 + | otherwise = length (roots t) == 2 * length (leaves t) - 3 + +prop_connect :: a -> Tree a -> Tree a -> Bool +prop_connect n l r + -- XXX: Skip not bifurcating trees. This is ugly, I know. + | not (bifurcating l) || not (bifurcating r) = True + | length (leaves l) < 3 || length (leaves r) < 3 = length (connect n l r) == 1 + | otherwise = length (connect n l r) == length (leaves l) * length (leaves r) + +type Constraint a = Subset a + +compatibleAll :: (Show a, Ord a) => Tree a -> [Constraint a] -> Bool +compatibleAll (Node _ [l, r]) cs = all (compatible (bipartition l)) cs && + all (compatible (bipartition r)) cs +compatibleAll _ _ = error "Tree is not bifurcating." + +compatibleWith :: (Show b, Ord b) => (a -> b) -> [Constraint a] -> Tree a -> Bool +compatibleWith f cs t = compatibleAll (fmap f t) (map (smap f) cs) + spec :: Spec spec = do describe "subTree" $ do @@ -64,3 +102,36 @@ it "leaves height constant for Measurable trees" $ height (prune subSampleLargeTree) `shouldBe` height subSampleLargeTree + describe "roots" $ do + it "correctly handles leaves and cherries" $ do + let tleaf = Node 0 [] :: Tree Int + tcherry = Node 0 [Node 1 [], Node 2 []] :: Tree Int + roots tleaf `shouldBe` [tleaf] + roots tcherry `shouldBe` [tcherry] + it "correctly handles simple trees" $ do + let simpleTre = Node "i" [ Node "j" [Node "x" [], Node "y" []], Node "z" [] ] :: Tree String + simpleSol = + [ Node "i" [ Node "j" [ Node "x" [] , Node "y" [] ] + , Node "z" [] ] + , Node "i" [ Node "x" [] + , Node "j" [ Node "y" [] , Node "z" [] ] ] + , Node "i" [ Node "j" [ Node "x" [] , Node "z" [] ] + , Node "y" [] ] ] + roots simpleTre `shouldBe` simpleSol + modifyMaxSize (*100) $ + it "returns the correct number of rooted trees for arbitrary trees" $ + property (prop_roots :: (Tree Int -> Bool)) + + describe "connect" $ + modifyMaxSize (*100) $ do + it "returns the correct number of rooted trees for arbitrary trees" $ + property (prop_connect :: Int -> Tree Int -> Tree Int -> Bool) + it "correctly connects sample trees without and with constraints" $ do + a <- parseFileWith oneNewick "data/ConnectA.tree" + b <- parseFileWith oneNewick "data/ConnectB.tree" + c <- parseFileWith manyNewick "data/ConnectConstraints.tree" + let ts = connect (PhyloLabel "" Nothing (Just 1.0)) a b + cs = concatMap clades c :: [Constraint (PhyloLabel L.ByteString)] + ts' = filter (compatibleWith getName cs) ts + length ts `shouldBe` 63 + length ts' `shouldBe` 15 diff --git a/test/ELynx/Export/Tree/NewickSpec.hs b/test/ELynx/Export/Tree/NewickSpec.hs --- a/test/ELynx/Export/Tree/NewickSpec.hs +++ b/test/ELynx/Export/Tree/NewickSpec.hs @@ -15,25 +15,38 @@ (spec ) where -import qualified Data.ByteString.Lazy.Char8 as L +import qualified Data.ByteString.Lazy.Char8 as L import Data.Tree import Test.Hspec +import Test.QuickCheck hiding (label) +import Test.QuickCheck.Instances.Containers () +import ELynx.Data.Tree.NamedTree import ELynx.Data.Tree.PhyloTree import ELynx.Export.Tree.Newick import ELynx.Import.Tree.Newick import ELynx.Tools.InputOutput -samplePhyloByteStringTree :: Tree PhyloByteStringLabel +samplePhyloByteStringTree :: Tree (PhyloLabel L.ByteString) samplePhyloByteStringTree = parseByteStringWith "Sample newick byte string" newick sampleNewickByteString1 sampleNewickByteString1 :: L.ByteString sampleNewickByteString1 = L.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);" +prop_invariant :: Tree (PhyloLabel Int) -> Bool +prop_invariant t = parseByteStringWith "Sample newick byte string" newick (toNewick t) == + -- We have to convert Int to ByteString, but we cannot use + -- arbitrary instances with ByteStrings because many + -- characters are disallowed and break the Newick string. + fmap (\l -> l{label = getName l}) t + spec :: Spec -spec = describe "parseByteStringWith newick $ toNewickPhyloByteString" $ - it "should be an invariant" $ do - let bs = toNewick samplePhyloByteStringTree - parseByteStringWith "Newick string converted from tree object" - newick bs `shouldBe` samplePhyloByteStringTree +spec = describe "parseByteStringWith newick $ toNewickPhyloByteString" $ do + it "should be an invariant for a sample tree" $ do + let bs = toNewick samplePhyloByteStringTree + parseByteStringWith "Newick string converted from tree object" + newick bs `shouldBe` samplePhyloByteStringTree + + it "should be an invariant for a sample tree" $ + property prop_invariant diff --git a/test/ELynx/Import/Tree/NewickSpec.hs b/test/ELynx/Import/Tree/NewickSpec.hs --- a/test/ELynx/Import/Tree/NewickSpec.hs +++ b/test/ELynx/Import/Tree/NewickSpec.hs @@ -25,24 +25,24 @@ sampleLabelByteString :: L.ByteString sampleLabelByteString = L.pack "name:0.3" -sampleLabel :: PhyloByteStringLabel -sampleLabel = PhyloLabel (L.pack "name") Nothing 0.3 +sampleLabel :: (PhyloLabel L.ByteString) +sampleLabel = PhyloLabel (L.pack "name") Nothing (Just 0.3) sampleForestByteString :: L.ByteString sampleForestByteString = L.pack "(l,l,(a,b))" -sampleForest :: [Tree PhyloByteStringLabel] +sampleForest :: [Tree (PhyloLabel L.ByteString)] sampleForest = - [ Node { rootLabel = PhyloLabel (L.pack "l") Nothing 0 + [ Node { rootLabel = PhyloLabel (L.pack "l") Nothing Nothing , subForest = [] } - , Node { rootLabel = PhyloLabel (L.pack "l") Nothing 0 + , Node { rootLabel = PhyloLabel (L.pack "l") Nothing Nothing , subForest = [] } - , Node { rootLabel = PhyloLabel (L.pack "") Nothing 0 + , Node { rootLabel = PhyloLabel (L.pack "") Nothing Nothing , subForest = - [ Node {rootLabel = PhyloLabel (L.pack "a") Nothing 0, subForest = []} - , Node {rootLabel = PhyloLabel (L.pack "b") Nothing 0, subForest = []} + [ Node {rootLabel = PhyloLabel (L.pack "a") Nothing Nothing, subForest = []} + , Node {rootLabel = PhyloLabel (L.pack "b") Nothing Nothing, subForest = []} ] } ] @@ -57,18 +57,18 @@ sampleNewickEmptyByteString :: L.ByteString sampleNewickEmptyByteString = L.pack "(,(,,),);" -sampleNewickEmpty :: Tree PhyloByteStringLabel +sampleNewickEmpty :: Tree (PhyloLabel L.ByteString) sampleNewickEmpty = - Node { rootLabel = PhyloLabel (L.pack "") Nothing 0 + Node { rootLabel = PhyloLabel (L.pack "") Nothing Nothing , subForest = - [ Node {rootLabel = PhyloLabel (L.pack "") Nothing 0, subForest = []} - , Node {rootLabel = PhyloLabel (L.pack "") Nothing 0, subForest = - [ Node {rootLabel = PhyloLabel (L.pack "") Nothing 0, subForest = []} - , Node {rootLabel = PhyloLabel (L.pack "") Nothing 0, subForest = []} - , Node {rootLabel = PhyloLabel (L.pack "") Nothing 0, subForest = []} + [ Node {rootLabel = PhyloLabel (L.pack "") Nothing Nothing, subForest = []} + , Node {rootLabel = PhyloLabel (L.pack "") Nothing Nothing, subForest = + [ Node {rootLabel = PhyloLabel (L.pack "") Nothing Nothing, subForest = []} + , Node {rootLabel = PhyloLabel (L.pack "") Nothing Nothing, subForest = []} + , Node {rootLabel = PhyloLabel (L.pack "") Nothing Nothing, subForest = []} ] } - , Node {rootLabel = PhyloLabel (L.pack "") Nothing 0, subForest = []} + , Node {rootLabel = PhyloLabel (L.pack "") Nothing Nothing, subForest = []} ] } @@ -76,10 +76,10 @@ spec = do describe "branchLength" $ do it "parses a colon and a branch length" $ - parse branchLength "" (L.pack ":13.2") `shouldParse` 13.2 + parse branchLength "" (L.pack ":13.2") `shouldParse` Just 13.2 - it "returns 0 if no branch length is given" $ - parse branchLength "" (L.pack "") `shouldParse` 0 + it "returns Nothing if no branch length is given" $ + parse branchLength "" (L.pack "") `shouldParse` Nothing describe "name" $ do it "parses a string of printable characters" $ @@ -96,7 +96,7 @@ parse node "" sampleLabelByteString `shouldParse` sampleLabel it "parses tree nodes with empty names and branch lengths" $ - parse node "" (L.pack "") `shouldParse` PhyloLabel (L.pack "") Nothing 0 + parse node "" (L.pack "") `shouldParse` PhyloLabel (L.pack "") Nothing Nothing describe "leaf" $ it "parses a leaf of a tree" $