tlynx 0.0.1 → 0.1.0
raw patch · 17 files changed
+1096/−339 lines, 17 filesdep +MonadRandomdep +arraydep ~elynx-seqdep ~elynx-toolsdep ~elynx-tree
Dependencies added: MonadRandom, array
Dependency ranges changed: elynx-seq, elynx-tools, elynx-tree, math-functions
Files
- ChangeLog.md +2/−2
- README.md +88/−78
- app/Main.hs +9/−5
- app/Options.hs +58/−26
- src/Compare/Compare.hs +158/−103
- src/Compare/Options.hs +33/−94
- src/Connect/Connect.hs +107/−0
- src/Connect/Options.hs +64/−0
- src/Distance/Distance.hs +157/−0
- src/Distance/Options.hs +147/−0
- src/Examine/Examine.hs +21/−7
- src/Examine/Options.hs +14/−5
- src/Shuffle/Options.hs +68/−0
- src/Shuffle/Shuffle.hs +142/−0
- src/Simulate/Options.hs +3/−5
- src/Simulate/Simulate.hs +4/−4
- tlynx.cabal +21/−10
ChangeLog.md view
@@ -1,10 +1,10 @@ # Table of Contents -1. [Changelog for ELynx](#org9916021)+1. [Changelog for ELynx](#org263dcac) -<a id="org9916021"></a>+<a id="org263dcac"></a> # Changelog for ELynx
README.md view
@@ -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)
app/Main.hs view
@@ -1,5 +1,3 @@-{-# LANGUAGE TemplateHaskell #-}- {- | Module : Main Description : Work with phylogenetic trees@@ -21,7 +19,10 @@ import Options import Compare.Compare+import Connect.Connect+import Distance.Distance import Examine.Examine+import Shuffle.Shuffle import Simulate.Simulate import ELynx.Tools.Logger@@ -34,6 +35,9 @@ lvl = verbosity g lf = (++ ".log") <$> fn case c of- Compare a -> runReaderT (eLynxWrapper lvl lf compareHeader $ compareTrees fn) a- Examine a -> runReaderT (eLynxWrapper lvl lf examineHeader $ examine fn) a- Simulate a -> runReaderT (eLynxWrapper lvl lf simulateHeader $ simulate fn) a+ Distance a -> runReaderT (eLynxWrapper lvl lf distanceDescription $ distance fn) a+ Examine a -> runReaderT (eLynxWrapper lvl lf examineDescription $ examine fn) a+ Simulate a -> runReaderT (eLynxWrapper lvl lf simulateDescription $ simulate fn) a+ Compare a -> runReaderT (eLynxWrapper lvl lf compareDescription $ compareCmd fn) a+ Connect a -> runReaderT (eLynxWrapper lvl lf connectDescription $ connectCmd fn) a+ Shuffle a -> runReaderT (eLynxWrapper lvl lf shuffleDescription $ shuffleCmd fn) a
app/Options.hs view
@@ -16,65 +16,97 @@ ( Arguments (..) , CommandArguments (..) , parseArguments- , compareHeader- , examineHeader- , simulateHeader+ , distanceDescription+ , examineDescription+ , simulateDescription+ , compareDescription+ , connectDescription+ , shuffleDescription ) where import Options.Applicative import Options.Applicative.Help.Pretty import Compare.Options+import Connect.Options+import Distance.Options import Examine.Options+import Shuffle.Options import Simulate.Options import ELynx.Tools.Options data CommandArguments =- Compare CompareArguments+ Distance DistanceArguments | Examine ExamineArguments | Simulate SimulateArguments+ | Compare CompareArguments+ | Connect ConnectArguments+ | Shuffle ShuffleArguments -compareHeader, examineHeader, simulateHeader :: String-compareHeader = "Compute distances between phylogenetic trees."-examineHeader = "Compute summary statistics of phylogenetic trees."-simulateHeader = "Simulate phylogenetic trees using birth and death processes."+distanceDescription :: String+distanceDescription = "Compute distances between many phylogenetic trees." -compareCommand :: Mod CommandFields CommandArguments-compareCommand = command "compare" $- info (Compare <$> compareArguments)+examineDescription :: String+examineDescription = "Compute summary statistics of phylogenetic trees."++simulateDescription :: String+simulateDescription = "Simulate phylogenetic trees using birth and death processes."++compareDescription :: String+compareDescription = "Compare two phylogenetic trees (compute distances and branch-wise differences)."++connectDescription :: String+connectDescription = "Connect two phylogenetic trees in all ways (possibly honoring constraints)."++shuffleDescription :: String+shuffleDescription = "Shuffle a phylogenetic tree (keep coalescent times, but shuffle topology and leaves)."++distanceCommand :: Mod CommandFields CommandArguments+distanceCommand = command "distance" $+ info (Distance <$> distanceArguments) ( fullDesc- <> header compareHeader- <> footerDoc (Just $ pretty compareFooter) )+ <> progDesc distanceDescription+ <> footerDoc (Just $ pretty distanceFooter) ) examineCommand :: Mod CommandFields CommandArguments examineCommand = command "examine" $ info (Examine <$> examineArguments)- $ header examineHeader+ $ progDesc examineDescription simulateCommand :: Mod CommandFields CommandArguments simulateCommand = command "simulate" $ info (Simulate <$> simulateArguments) ( fullDesc- <> header simulateHeader- <> progDesc simulateDesc+ <> progDesc simulateDescription <> footerDoc (Just $ pretty simulateFooter) ) +compareCommand :: Mod CommandFields CommandArguments+compareCommand = command "compare" $+ info (Compare <$> compareArguments)+ $ progDesc compareDescription++connectCommand :: Mod CommandFields CommandArguments+connectCommand = command "connect" $+ info (Connect <$> connectArguments)+ $ progDesc connectDescription++shuffleCommand :: Mod CommandFields CommandArguments+shuffleCommand = command "shuffle" $+ info (Shuffle <$> shuffleArguments)+ $ progDesc shuffleDescription+ commandArguments :: Parser CommandArguments commandArguments = hsubparser $- compareCommand+ distanceCommand <> examineCommand <> simulateCommand--data Arguments = Arguments- { globalArgs :: GlobalArguments- , cmdArgs :: CommandArguments }+ <> compareCommand+ <> connectCommand+ <> shuffleCommand -parseArguments :: IO Arguments-parseArguments = parseArgumentsWith desc ftr $- Arguments- <$> globalArguments- <*> commandArguments+parseArguments :: IO (Arguments CommandArguments)+parseArguments = parseArgumentsWith desc ftr commandArguments desc :: [String] desc = [ "Compare, examine, and simulate phylogenetic trees." ]
src/Compare/Compare.hs view
@@ -2,7 +2,8 @@ {-# LANGUAGE TemplateHaskell #-} {- |-Description : Compute distances between trees+Module : Compare.Compare+Description : Compare two phylogenies Copyright : (c) Dominik Schrempf 2019 License : GPL-3 @@ -10,17 +11,13 @@ Stability : unstable Portability : portable -Creation date: Wed May 29 18:09:39 2019.--- Symmetric (Robinson-Foulds) distance.-- Incompatible splits distance.+Creation date: Thu Sep 19 15:01:52 2019. -} module Compare.Compare- ( compareTrees- )-where+ ( compareCmd+ ) where import Control.Monad import Control.Monad.IO.Class@@ -28,117 +25,175 @@ import Control.Monad.Trans.Class import Control.Monad.Trans.Reader import qualified Data.ByteString.Lazy.Char8 as L+import Data.List+import qualified Data.Map as M+import Data.Monoid+import qualified Data.Set as S import qualified Data.Text as T+import qualified Data.Text.Encoding as E import qualified Data.Text.IO as T-import qualified Data.Text.Lazy as LT-import qualified Data.Text.Lazy.Encoding as LT+ import Data.Tree-import qualified Data.Vector.Unboxed as V-import Statistics.Sample import System.IO import Text.Printf import Compare.Options -import ELynx.Data.Tree.BranchSupportTree as B+import ELynx.Data.Tree.Bipartition+import ELynx.Data.Tree.BranchSupportTree as BS import ELynx.Data.Tree.Distance-import ELynx.Data.Tree.MeasurableTree as M+import ELynx.Data.Tree.MeasurableTree import ELynx.Data.Tree.NamedTree import ELynx.Data.Tree.PhyloTree-import ELynx.Export.Tree.Newick+import ELynx.Data.Tree.Tree+import ELynx.Export.Tree.Newick (toNewick) import ELynx.Import.Tree.Newick-import ELynx.Tools.ByteString (alignLeft, alignRight) import ELynx.Tools.InputOutput-import ELynx.Tools.Logger -pf :: String-pf = "%.3f"--header :: Int -> Distance -> L.ByteString-header n d = alignLeft (n+2) "Tree 1"- <> alignLeft (n+2) "Tree 2"- <> alignRight 20 (L.pack $ show d)+treesOneFile :: FilePath+ -> Compare (Tree (PhyloLabel L.ByteString), Tree (PhyloLabel L.ByteString))+treesOneFile tf = do+ a <- lift ask+ let nw = if argsNewickIqTree a then manyNewickIqTree else manyNewick+ $(logInfo) $ T.pack $ "Parse file '" ++ tf ++ "'."+ ts <- liftIO $ parseFileWith nw tf+ let n = length ts+ case compare n 2 of+ LT -> error "Not enough trees in file."+ GT -> error "Too many trees in file."+ EQ -> return (head ts, head . tail $ ts) -showTriplet :: (PrintfArg a) => Int -> [String] -> (Int, Int, a) -> L.ByteString-showTriplet n args (i, j, d) = i' <> j' <> d'- where i' = alignLeft (n+2) $ L.pack (args !! i)- j' = alignLeft (n+2) $ L.pack (args !! j)- d' = alignRight 20 $ L.pack (printf pf d)+treesTwoFiles :: FilePath -> FilePath+ -> Compare (Tree (PhyloLabel L.ByteString), Tree (PhyloLabel L.ByteString))+treesTwoFiles tf1 tf2 = do+ a <- lift ask+ let nw = if argsNewickIqTree a then oneNewickIqTree else oneNewick+ $(logInfo) $ T.pack $ "Parse first tree file '" ++ tf1 ++ "'."+ t1 <- liftIO $ parseFileWith nw tf1+ $(logInfo) $ T.pack $ "Parse second tree file '" ++ tf2 ++ "'."+ t2 <- liftIO $ parseFileWith nw tf2+ return (t1, t2) --- | Compute distance functions between phylogenetic trees.-compareTrees :: Maybe FilePath -> Compare ()-compareTrees outFileBN = do+-- | More detailed comparison of two trees.+compareCmd :: Maybe FilePath -> Compare ()+compareCmd outFile = do a <- lift ask -- Determine output handle (stdout or file).- let outFile = (++ ".out") <$> outFileBN- outH <- liftIO $ maybe (pure stdout) (`openFile` WriteMode) outFile- let tfps = argsInFiles a- (trees, names) <-- if length tfps <= 1- then- do ts <- if null tfps- then do $(logInfo) "Read trees from standard input."- liftIO $ parseIOWith manyNewick- else do let f = head tfps- $(logInfo) $ T.pack $ "Read trees from file: " <> f <> "."- liftIO $ parseFileWith manyNewick f- let n = length ts- when (n <= 1) (error "Not enough trees found in file.")- $(logInfo) "Compute pairwise distances between trees in the same file."- $(logInfo) $ T.pack $ "Trees are numbered from 0 to " ++ show (n-1) ++ "."- return (ts, take n (map show [0 :: Int ..]))- else- do $(logInfo) "Read trees from files."- ts <- liftIO $ mapM (parseFileWith newick) tfps- when (length ts <= 1) (error "Not enough trees found in files.")- $(logInfo) "Compute pairwise distances between trees from different files."- $(logInfo) "Trees are named according to their file names."- return (ts, tfps)- $(logDebug) "The trees are:"- $(logDebug) $ LT.toStrict $ LT.decodeUtf8 $ L.unlines $ map toNewick trees- case outFile of- Nothing -> logNewSection "Write results to standard output."- Just f -> logNewSection $ T.pack $ "Write results to file " <> f <> "."- let n = maximum $ 6 : map length names- distance = argsDistance a- case distance of- Symmetric -> $(logInfo) "Use symmetric (Robinson-Foulds) distance."- IncompatibleSplit val -> do- $(logInfo) "Use incompatible split distance."- $(logInfo) $ T.pack $ "Collapse nodes with support less than " ++ show val ++ "."- BranchScore False -> $(logInfo) "Use branch score distance without normalization."- BranchScore True -> $(logInfo) "Use branch score distance with normalization."- let distanceMeasure :: Tree PhyloByteStringLabel -> Tree PhyloByteStringLabel -> Double- distanceMeasure = case distance of- Symmetric -> \t1 t2 -> fromIntegral $ symmetricDistanceWith getName t1 t2- IncompatibleSplit _ -> \t1 t2 -> fromIntegral $ incompatibleSplitsDistanceWith getName t1 t2- BranchScore _ -> branchScoreDistance- normalizeF = case distance of- BranchScore True -> M.normalize- _ -> id- collapseF = case distance of- -- For the incompatible split distance we have to collapse branches with- -- support lower than the given value. Before doing so, we normalize the- -- branch support values.- IncompatibleSplit val -> collapse val . B.normalize- _ -> id- trees' = map (collapseF . normalizeF) trees- $(logDebug) "The prepared trees are:"- $(logDebug) $ LT.toStrict $ LT.decodeUtf8 $ L.unlines $ map toNewick trees'- let dsTriplets = computePairwiseDistances distanceMeasure trees'- ds = map (\(_, _, x) -> x) dsTriplets- dsVec = V.fromList ds- -- XXX: It may be good to use the common 'io' function also here.- liftIO $ hPutStrLn outH $ "Summary statistics of " ++ show distance ++ " Distance:"- liftIO $ T.hPutStrLn outH $ T.justifyLeft 10 ' ' "Mean: " <> T.pack (printf pf (mean dsVec))- liftIO $ T.hPutStrLn outH $ T.justifyLeft 10 ' ' "Variance: " <> T.pack (printf pf (variance dsVec))- -- L.putStrLn $ L.unlines $ map toNewick ts- -- L.putStrLn $ L.unlines $ map toNewick tsN- -- L.putStrLn $ L.unlines $ map toNewick tsC- lift $ unless (argsSummaryStatistics a) (- do- lift $ hPutStrLn outH ""- lift $ L.hPutStrLn outH $ header n distance- lift $ L.hPutStr outH $ L.unlines (map (showTriplet n names) dsTriplets)- )+ let outFn = (++ ".out") <$> outFile+ outH <- outHandle "results" outFn+ -- liftIO $ hPutStrLn outH ""++ -- Read input.+ let inFiles = argsInFiles a+ nFiles = length inFiles+ (t1, t2) <- case nFiles of+ 1 -> treesOneFile (head inFiles)+ 2 -> treesTwoFiles (head inFiles) (head . tail $ inFiles)+ _ -> error "Need two input files with one tree each or one input file with two trees."+ -- liftIO $ hPutStrLn outH ""++ liftIO $ hPutStrLn outH "Tree 1:"+ liftIO $ L.hPutStrLn outH $ toNewick t1+ -- liftIO $ hPrint outH t1+ liftIO $ hPutStrLn outH "Tree 2:"+ liftIO $ L.hPutStrLn outH $ toNewick t2+ -- liftIO $ hPrint outH t2+ liftIO $ hPutStrLn outH ""++ -- Check input.+ let lvs1 = leaves t1+ lvs2 = leaves t2+ lfns1 = map getName lvs1+ lfns2 = map getName lvs2+ s1 = S.fromList lfns1+ s2 = S.fromList lfns2+ if s1 == s2+ then liftIO $ hPutStrLn outH "Trees have the same set of leaf names."+ else liftIO $ hPutStrLn outH "Trees do not have the same set of leaf names."+ liftIO $ hPutStrLn outH ""++ -- Distances.+ let formatD str val = T.justifyLeft 25 ' ' str <> " " <> val+ liftIO $ hPutStrLn outH "Distances."+ liftIO $ T.hPutStrLn outH $ formatD "Symmetric"+ (T.pack $ show $ symmetric t1 t2)+ liftIO $ T.hPutStrLn outH $ formatD "Branch score"+ (T.pack $ show $ branchScoreWith getName getLen t1 t2)+ let t1' = BS.normalize t1+ t2' = BS.normalize t2+ $(logDebug) "Trees with normalized branch support values:"+ $(logDebug) $ E.decodeUtf8 $ L.toStrict $ toNewick t1'+ $(logDebug) $ E.decodeUtf8 $ L.toStrict $ toNewick t2'+ liftIO $ T.hPutStrLn outH $ formatD "Incompatible split (0.10)"+ (T.pack $ show $ incompatibleSplits (collapse 0.1 t1') (collapse 0.1 t2'))+ liftIO $ T.hPutStrLn outH $ formatD "Incompatible split (0.50)"+ (T.pack $ show $ incompatibleSplits (collapse 0.5 t1') (collapse 0.5 t2'))+ -- liftIO $ T.hPutStrLn outH $ formatD "Incompatible split (0.60)"+ -- (T.pack $ show $ incompatibleSplits (collapse 0.6 t1') (collapse 0.6 t2'))+ -- liftIO $ T.hPutStrLn outH $ formatD "Incompatible split (0.70)"+ -- (T.pack $ show $ incompatibleSplits (collapse 0.7 t1') (collapse 0.7 t2'))+ liftIO $ T.hPutStrLn outH $ formatD "Incompatible split (0.80)"+ (T.pack $ show $ incompatibleSplits (collapse 0.8 t1') (collapse 0.8 t2'))+ liftIO $ T.hPutStrLn outH $ formatD "Incompatible split (0.90)"+ (T.pack $ show $ incompatibleSplits (collapse 0.9 t1') (collapse 0.9 t2'))+ -- liftIO $ T.hPutStrLn outH $ formatD "Incompatible split (1.01)"+ -- (T.pack $ show $ incompatibleSplits (collapse 1.01 t1') (collapse 1.01 t2'))+ -- liftIO $ L.hPutStrLn outH $ toNewick (collapse 1.01 t1')++ -- Bipartitions.+ when (argsBipartitions a) (do+ let bp1 = bipartitions (fmap getName t1)+ bp2 = bipartitions (fmap getName t2)+ bp1Only = bp1 S.\\ bp2+ bp2Only = bp2 S.\\ bp1+ unless (S.null bp1Only)+ (do+ liftIO $ hPutStrLn outH ""+ liftIO $ hPutStrLn outH "Bipartitions in Tree 1 that are not in Tree 2."+ -- let bp1Strs = map (bphuman L.unpack . bpmap getName) (S.toList bp1Only)+ forM_ bp1Only (liftIO . hPutStrLn outH . bphuman L.unpack))+ -- let bp1Strs = map (bphuman L.unpack) (S.toList bp1Only)+ -- liftIO $ hPutStrLn outH $ intercalate "\n" bp1Strs)+ unless (S.null bp2Only)+ (do+ liftIO $ hPutStrLn outH ""+ liftIO $ hPutStrLn outH "Bipartitions in Tree 2 that are not in Tree 1."+ forM_ bp2Only (liftIO . hPutStrLn outH . bphuman L.unpack))++ -- Common bipartitions and their respective differences in branch lengths.+ liftIO $ hPutStrLn outH ""+ let bpCommon = bp1 `S.intersection` bp2+ if S.null bpCommon+ then liftIO $ hPutStrLn outH "There are no common bipartitions."+ else do+ let bpToBrLen1 = M.map getSum $ bipartitionToBranchLength getName (Sum . getLen) t1+ bpToBrLen2 = M.map getSum $ bipartitionToBranchLength getName (Sum . getLen) t2+ liftIO $ hPutStrLn outH "Common bipartitions and their respective differences in branch lengths."+ -- Header.+ liftIO $ hPutStrLn outH header+ forM_ bpCommon (liftIO . hPutStrLn outH . getCommonBpStr L.unpack bpToBrLen1 bpToBrLen2))+ liftIO $ hClose outH+ -- TODO. Plot using haskell-chart.+ -- See https://github.com/timbod7/haskell-chart/wiki.++header :: String+header = intercalate " " $ cols ++ ["Bipartition"]+ where cols = map (T.unpack . T.justifyRight 12 ' ')+ [ "Length 1", "Length 2", "Delta", "Relative [%]"]++getCommonBpStr :: (Ord a, Fractional b, PrintfArg b)+ => (a -> String)+ -> M.Map (Bipartition a) b -> M.Map (Bipartition a) b+ -> Bipartition a -> String+getCommonBpStr f m1 m2 p = intercalate " "+ [ printf "% 12.7f" l1+ , printf "% 12.7f" l2+ , printf "% 12.7f" d+ , printf "% 12.7f" rd+ , s ]+ where l1 = m1 M.! p+ l2 = m2 M.! p+ d = l1 - l2+ rd = 2 * d / (l1 + l2)+ s = bphuman f p
src/Compare/Options.hs view
@@ -1,6 +1,6 @@ {- | Module : Compare.Options-Description : Options of tree-dist+Description : Options for the compare subcommand Copyright : (c) Dominik Schrempf 2019 License : GPL-3 @@ -8,119 +8,58 @@ Stability : unstable Portability : portable -Creation date: Thu Aug 29 13:02:22 2019.+Creation date: Thu Sep 19 15:02:21 2019. -} module Compare.Options ( CompareArguments (..) , Compare- , Distance (..) , compareArguments- , compareFooter ) where import Control.Monad.Logger import Control.Monad.Trans.Reader-import Data.List-import Data.Void import Options.Applicative-import Text.Megaparsec (Parsec, try)-import Text.Megaparsec.Char (char, string)-import Text.Megaparsec.Char.Lexer (float)-import Text.Printf -import ELynx.Tools.Options---- | Supported distance measures.-data Distance =- Symmetric -- ^ Symmetric distance.- | IncompatibleSplit Double -- ^ Incompatible split distance; collapse nodes- -- with branch support below given value.- | BranchScore Bool -- ^ Branch score distance. If given, normalize- -- the trees before comparison.--instance Show Distance where- show Symmetric = "Symmetric"- show (IncompatibleSplit c) = "Incompatible Split (" ++ printf "%.1f" c ++ ")"- show (BranchScore False) = "Branch Score"- show (BranchScore True ) = "Branch Score (normalized)"---- | Arguments needed to compute distance measures.+-- | Arguments of compare command. data CompareArguments = CompareArguments- { argsDistance :: Distance- , argsSummaryStatistics :: Bool- , argsInFiles :: [FilePath]- }+ { argsNormalize :: Bool+ , argsBipartitions :: Bool+ , argsNewickIqTree :: Bool+ , argsInFiles :: [FilePath] } --- | Lgger and reader data type.++-- | Logger and reader data type. type Compare = LoggingT (ReaderT CompareArguments IO) --- | COmmand line parser.+-- | Parse arguments of compare command. compareArguments :: Parser CompareArguments compareArguments = CompareArguments <$>- distanceOpt- <*> summaryStatisticsSwitch- <*> many inFilesArg--inFilesArg :: Parser FilePath-inFilesArg = strArgument $- metavar "INPUT-FILES" <>- help "Read tree(s) from INPUT-FILES; if more files are given, one tree is expected per file"--symmetric :: Parsec Void String Distance-symmetric = do- _ <- string "symmetric"- pure Symmetric--incompatibleSplit :: Parsec Void String Distance-incompatibleSplit = do- _ <- string "incompatible-split"- _ <- char '['- f <- float- _ <- char ']'- if (0 < f) && (f < 1)- then pure $ IncompatibleSplit f- else error "Branch support has to be between 0 and 1."--branchScore :: Parsec Void String Distance-branchScore = do- _ <- string "branch-score"- pure (BranchScore False)--branchScoreNormalized :: Parsec Void String Distance-branchScoreNormalized = do- _ <- string "branch-score-normalized"- pure (BranchScore True)+ normalize+ <*> bipartitions+ <*> newickIqTree+ <*> some file -distanceParser :: Parsec Void String Distance-distanceParser = try symmetric- <|> try incompatibleSplit- -- Try first the normalized one, since the normal branch score- -- parser also succeeds in this case.- <|> try branchScoreNormalized- <|> try branchScore+normalize :: Parser Bool+normalize = switch $+ long "normalize"+ <> short 'n'+ <> help "Normalize trees before comparison" -distanceOpt :: Parser Distance-distanceOpt = option (megaReadM distanceParser) $- long "distance" <>- short 'd' <>- metavar "MEASURE" <>- help "Type of distance to calculate (available distance measures are listed below)"+bipartitions :: Parser Bool+bipartitions = switch $+ long "bipartitions"+ <> short 'b'+ <> help "Print common and missing bipartitions" -summaryStatisticsSwitch :: Parser Bool-summaryStatisticsSwitch = switch $- long "summary-statistics" <>- short 's' <>- help "Report summary statistics only"+file :: Parser FilePath+file = strArgument $+ metavar "NAME"+ <> help "Tree file" --- | Information about provided distance types.-compareFooter :: String-compareFooter = intercalate "\n"- [ "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."- , " Branch score distance : -d branch-score"- , " -d branch-score-normalized (normalize trees before comparison)" ]+newickIqTree :: Parser Bool+newickIqTree = switch $+ long "newick-iqtree"+ <> short 'i'+ <> help "Use IQ-TREE Newick format (internal node labels are branch support values)"
+ src/Connect/Connect.hs view
@@ -0,0 +1,107 @@+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE TemplateHaskell #-}++{- |+Module : Connect.Connect+Description : Connect two phylogenies+Copyright : (c) Dominik Schrempf 2019+License : GPL-3++Maintainer : dominik.schrempf@gmail.com+Stability : unstable+Portability : portable++Creation date: Thu Sep 19 15:01:52 2019.++-}++module Connect.Connect+ ( connectCmd+ ) where++import Control.Monad.IO.Class+import Control.Monad.Logger+import Control.Monad.Trans.Class+import Control.Monad.Trans.Reader+import qualified Data.ByteString.Lazy.Char8 as L+import Data.Tree+import System.IO++import Connect.Options++import ELynx.Data.Tree.Bipartition (bipartition, compatible)+import ELynx.Data.Tree.NamedTree (getName)+import ELynx.Data.Tree.PhyloTree (PhyloLabel (PhyloLabel))+import ELynx.Data.Tree.Subset (Subset, smap)+import ELynx.Data.Tree.Tree (clades, connect)+import ELynx.Export.Tree.Newick (toNewick)+import ELynx.Import.Tree.Newick (manyNewick, oneNewick)+import ELynx.Tools.InputOutput (outHandle, parseFileWith)+import ELynx.Tools.Text (fromBs, tShow)++-- TODO: Write a proper documentation.++-- | Connect two trees honoring possible constraints.+connectCmd :: Maybe FilePath -> Connect ()+connectCmd outFile = do+ -- Determine output handle (stdout or file).+ a <- lift ask+ let outFn = (++ ".out") <$> outFile+ outH <- outHandle "results" outFn++ -- Do we have constraints or not?+ let cs = constraints a+ l = inFileA a+ r = inFileB a+ case cs of+ Nothing -> connectOnly outH l r+ Just c -> connectAndFilter outH c l r++ liftIO $ hClose outH++connectTrees :: Tree (PhyloLabel L.ByteString)+ -> Tree (PhyloLabel L.ByteString)+ -> [Tree (PhyloLabel L.ByteString)]+connectTrees = connect (PhyloLabel "" Nothing Nothing)++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)++parseTrees :: FilePath -> FilePath+ -> Connect (Tree (PhyloLabel L.ByteString), Tree (PhyloLabel L.ByteString))+parseTrees l r = do+ tl <- liftIO $ parseFileWith oneNewick l+ tr <- liftIO $ parseFileWith oneNewick r+ $(logInfo) "Tree 1:"+ $(logInfo) $ fromBs $ toNewick tl+ $(logInfo) "Tree 2:"+ $(logInfo) $ fromBs $ toNewick tr+ return (tl, tr)++connectOnly :: Handle -> FilePath -> FilePath -> Connect ()+connectOnly h l r = do+ (tl, tr) <- parseTrees l r+ let ts = connectTrees tl tr+ $(logInfo) $ "Connected trees: " <> tShow (length ts)+ liftIO $ L.hPutStr h $ L.unlines $ map toNewick ts++connectAndFilter :: Handle -> FilePath -> FilePath -> FilePath -> Connect ()+connectAndFilter h c l r = do+ cts <- liftIO $ parseFileWith manyNewick c+ $(logInfo) "Constraints:"+ $(logInfo) $ fromBs $ L.intercalate "\n" $ map toNewick cts+ (tl, tr) <- parseTrees l r+ let ts = connectTrees tl tr+ cs = concatMap clades cts :: [Constraint (PhyloLabel L.ByteString)]+ -- Only collect trees that are compatible with the constraints.+ ts' = filter (compatibleWith getName cs) ts+ $(logInfo) $ "Connected trees: " <> tShow (length ts)+ $(logInfo) $ "Compatible trees: " <> tShow (length ts')+ liftIO $ L.hPutStr h $ L.unlines $ map toNewick ts'
+ src/Connect/Options.hs view
@@ -0,0 +1,64 @@+{- |+Module : Connect.Options+Description : Options for the connect subcommand+Copyright : (c) Dominik Schrempf 2019+License : GPL-3++Maintainer : dominik.schrempf@gmail.com+Stability : unstable+Portability : portable++Creation date: Thu Sep 19 15:02:21 2019.++-}++module Connect.Options+ ( ConnectArguments (..)+ , Connect+ , connectArguments+ ) where++import Control.Monad.Logger+import Control.Monad.Trans.Reader+import Options.Applicative++-- | Arguments of connect command.+data ConnectArguments = ConnectArguments+ { newickIqTreeFlag :: Bool+ , constraints :: Maybe FilePath+ , inFileA :: FilePath+ , inFileB :: FilePath }++-- | Logger and reader data type.+type Connect = LoggingT (ReaderT ConnectArguments IO)++-- | Parse arguments of connect command.+connectArguments :: Parser ConnectArguments+connectArguments = ConnectArguments+ <$> newickIqTree+ <*> constraintsFile+ <*> fileA+ <*> fileB++newickIqTree :: Parser Bool+newickIqTree = switch $+ long "newick-iqtree"+ <> short 'i'+ <> help "Use IQ-TREE Newick format (internal node labels are branch support values)"++constraintsFile :: Parser (Maybe FilePath)+constraintsFile = optional $ strOption $+ metavar "CONSTRAINTS"+ <> short 'c'+ <> long "contraints"+ <> help "File containing one or more Newick trees to be used as constraints"++fileA :: Parser FilePath+fileA = strArgument $+ metavar "TREE-FILE-A"+ <> help "File containing the first Newick tree"++fileB :: Parser FilePath+fileB = strArgument $+ metavar "TREE-FILE-B"+ <> help "File containing the second Newick tree"
+ src/Distance/Distance.hs view
@@ -0,0 +1,157 @@+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE TemplateHaskell #-}++{- |+Description : Compute distances between trees+Copyright : (c) Dominik Schrempf 2019+License : GPL-3++Maintainer : dominik.schrempf@gmail.com+Stability : unstable+Portability : portable++Creation date: Wed May 29 18:09:39 2019.++- Symmetric (Robinson-Foulds) distance.+- Incompatible splits distance.++-}++module Distance.Distance+ ( distance+ )+where++import Control.Monad+import Control.Monad.IO.Class+import Control.Monad.Logger+import Control.Monad.Trans.Class+import Control.Monad.Trans.Reader+import qualified Data.ByteString.Lazy.Char8 as L+import Data.Maybe+import qualified Data.Text as T+import qualified Data.Text.IO as T+import qualified Data.Text.Lazy as LT+import qualified Data.Text.Lazy.Encoding as LT+import Data.Tree+import qualified Data.Vector.Unboxed as V+import Statistics.Sample+import System.IO+import Text.Printf++import Distance.Options++import ELynx.Data.Tree.BranchSupportTree as B+import ELynx.Data.Tree.Distance+import ELynx.Data.Tree.MeasurableTree as M+import ELynx.Data.Tree.PhyloTree+import ELynx.Export.Tree.Newick+import ELynx.Import.Tree.Newick+import ELynx.Tools.ByteString (alignLeft, alignRight)+import ELynx.Tools.InputOutput+import ELynx.Tools.Logger++pf :: String+pf = "%.3f"++header :: Int -> DistanceMeasure -> L.ByteString+header n d = alignLeft (n+2) "Tree 1"+ <> alignLeft (n+2) "Tree 2"+ <> alignRight 20 (L.pack $ show d)++showTriplet :: (PrintfArg a) => Int -> [String] -> (Int, Int, a) -> L.ByteString+showTriplet n args (i, j, d) = i' <> j' <> d'+ where i' = alignLeft (n+2) $ L.pack (args !! i)+ j' = alignLeft (n+2) $ L.pack (args !! j)+ d' = alignRight 20 $ L.pack (printf pf d)++-- | Compute distance functions between phylogenetic trees.+distance :: Maybe FilePath -> Distance ()+distance outFileBN = do+ a <- lift ask+ let nw = if argsNewickIqTree a then manyNewickIqTree else manyNewick+ -- Determine output handle (stdout or file).+ let outFile = (++ ".out") <$> outFileBN+ outH <- outHandle "results" outFile+ -- Master tree (in case it is given).+ let mname = argsMasterTreeFile a+ mtree <- case mname of+ Nothing -> return Nothing+ Just f -> do $(logInfo) $ T.pack $ "Read master tree from file: " <> f <> "."+ ts <- liftIO $ parseFileWith nw f+ let n = length ts+ when (n > 1) (error "More than one tree found in master file.")+ $(logInfo) "Compute distances between all trees and master tree."+ $(logInfo) $ T.pack $ "Trees are numbered from 0 to " ++ show (n-1) ++ "."+ return $ Just (head ts)+ let tfps = argsInFiles a+ (trees, names) <-+ if length tfps <= 1+ then+ do ts <- if null tfps+ then do $(logInfo) "Read trees from standard input."+ liftIO $ parseIOWith nw+ else do let f = head tfps+ $(logInfo) $ T.pack $ "Read trees from file: " <> f <> "."+ liftIO $ parseFileWith nw f+ let n = length ts+ when (n < 1) (error "Not enough trees found in file.")+ when (isNothing mtree) $ $(logInfo) "Compute pairwise distances between trees in the same file."+ $(logInfo) $ T.pack $ "Trees are numbered from 0 to " ++ show (n-1) ++ "."+ return (ts, take n (map show [0 :: Int ..]))+ else+ do $(logInfo) "Read trees from files."+ ts <- liftIO $ mapM (parseFileWith newick) tfps+ when (length ts <= 1) (error "Not enough trees found in files.")+ when (isNothing mtree) $ $(logInfo) "Compute pairwise distances between trees from different files."+ $(logInfo) "Trees are named according to their file names."+ return (ts, tfps)+ $(logDebug) "The trees are:"+ $(logDebug) $ LT.toStrict $ LT.decodeUtf8 $ L.unlines $ map toNewick trees+ case outFile of+ Nothing -> logNewSection "Write results to standard output."+ Just f -> logNewSection $ T.pack $ "Write results to file " <> f <> "."+ let n = maximum $ 6 : map length names+ dist = argsDistance a+ case dist of+ Symmetric -> $(logInfo) "Use symmetric (Robinson-Foulds) distance."+ IncompatibleSplit val -> do+ $(logInfo) "Use incompatible split distance."+ $(logInfo) $ T.pack $ "Collapse nodes with support less than " ++ show val ++ "."+ BranchScore -> $(logInfo) "Use branch score distance."+ when (argsNormalize a) $ $(logInfo) "Normalize trees before calculation of distances."+ let distanceMeasure :: Tree (PhyloLabel L.ByteString) -> Tree (PhyloLabel L.ByteString) -> Double+ distanceMeasure = case dist of+ Symmetric -> \t1 t2 -> fromIntegral $ symmetric t1 t2+ IncompatibleSplit _ -> \t1 t2 -> fromIntegral $ incompatibleSplits t1 t2+ BranchScore -> branchScore+ normalizeF = if argsNormalize a then M.normalize else id+ collapseF = case dist of+ -- For the incompatible split distance we have to collapse branches with+ -- support lower than the given value. Before doing so, we normalize the+ -- branch support values.+ IncompatibleSplit val -> collapse val . B.normalize+ _ -> id+ trees' = map (collapseF . normalizeF) trees+ $(logDebug) "The prepared trees are:"+ $(logDebug) $ LT.toStrict $ LT.decodeUtf8 $ L.unlines $ map toNewick trees'+ let dsTriplets = case mtree of+ Nothing -> pairwise distanceMeasure trees'+ Just t -> [ (0, i, distanceMeasure t t') | ( i, t' ) <- zip [1..] trees' ]+ ds = map (\(_, _, x) -> x) dsTriplets+ dsVec = V.fromList ds+ liftIO $ hPutStrLn outH $ "Summary statistics of " ++ show dist ++ " Distance:"+ liftIO $ T.hPutStrLn outH $ T.justifyLeft 10 ' ' "Mean: " <> T.pack (printf pf (mean dsVec))+ liftIO $ T.hPutStrLn outH $ T.justifyLeft 10 ' ' "Variance: " <> T.pack (printf pf (variance dsVec))+ -- L.putStrLn $ L.unlines $ map toNewick ts+ -- L.putStrLn $ L.unlines $ map toNewick tsN+ -- L.putStrLn $ L.unlines $ map toNewick tsC+ lift $ unless (argsSummaryStatistics a) (+ do+ lift $ hPutStrLn outH ""+ lift $ L.hPutStrLn outH $ header n dist+ case mname of+ Nothing -> lift $ L.hPutStr outH $ L.unlines (map (showTriplet n names) dsTriplets)+ Just mn -> lift $ L.hPutStr outH $ L.unlines (map (showTriplet n (mn : names)) dsTriplets)+ )+ liftIO $ hClose outH
+ src/Distance/Options.hs view
@@ -0,0 +1,147 @@+{- |+Module : Distance.Options+Description : Options of tree-dist+Copyright : (c) Dominik Schrempf 2019+License : GPL-3++Maintainer : dominik.schrempf@gmail.com+Stability : unstable+Portability : portable++Creation date: Thu Aug 29 13:02:22 2019.++-}++module Distance.Options+ ( DistanceArguments (..)+ , DistanceMeasure (..)+ , Distance+ , distanceArguments+ , distanceFooter+ ) where++import Control.Monad.Logger+import Control.Monad.Trans.Reader+import Data.List+import Data.Void+import Options.Applicative+import Text.Megaparsec (Parsec, eof, try)+import Text.Megaparsec.Char (char, string)+import Text.Megaparsec.Char.Lexer (float)+import Text.Printf++import ELynx.Tools.Options++-- | Supported distance measures.+data DistanceMeasure =+ Symmetric -- ^ Symmetric distance.+ | IncompatibleSplit Double -- ^ Incompatible split distance; collapse nodes+ -- with branch support below given value.+ | BranchScore -- ^ Branch score distance.++instance Show DistanceMeasure where+ show Symmetric = "Symmetric"+ show (IncompatibleSplit c) = "Incompatible Split (" ++ printf "%.1f" c ++ ")"+ show BranchScore = "Branch Score"++-- | Arguments needed to compute distance measures.+data DistanceArguments = DistanceArguments+ { argsDistance :: DistanceMeasure+ , argsNormalize :: Bool+ , argsSummaryStatistics :: Bool+ , argsMasterTreeFile :: Maybe FilePath+ , argsNewickIqTree :: Bool+ , argsInFiles :: [FilePath]+ }++-- | Logger and reader data type.+type Distance = LoggingT (ReaderT DistanceArguments IO)++-- | COmmand line parser.+distanceArguments :: Parser DistanceArguments+distanceArguments = DistanceArguments <$>+ distanceOpt+ <*> normalizeSwitch+ <*> summaryStatisticsSwitch+ <*> masterTreeFile+ <*> newickIqTree+ <*> many inFilesArg++masterTreeFile :: Parser (Maybe FilePath)+masterTreeFile = optional $ strOption $+ long "master-tree-file" <>+ short 'm' <>+ metavar "MASTER-TREE-File" <>+ help "Compare all trees to the tree in the master tree file."++inFilesArg :: Parser FilePath+inFilesArg = strArgument $+ metavar "INPUT-FILES" <>+ help "Read tree(s) from INPUT-FILES; if more files are given, one tree is expected per file"++symmetric :: Parsec Void String DistanceMeasure+symmetric = do+ _ <- string "symmetric"+ _ <- eof+ pure Symmetric++incompatibleSplit :: Parsec Void String DistanceMeasure+incompatibleSplit = do+ _ <- string "incompatible-split"+ _ <- char '['+ f <- float+ _ <- char ']'+ _ <- eof+ if (0 < f) && (f < 1)+ then pure $ IncompatibleSplit f+ else error "Branch support has to be between 0 and 1."++branchScore :: Parsec Void String DistanceMeasure+branchScore = do+ _ <- string "branch-score"+ _ <- eof+ pure BranchScore++distanceParser :: Parsec Void String DistanceMeasure+distanceParser = try symmetric+ <|> try incompatibleSplit+ -- Try first the normalized one, since the normal branch score+ -- parser also succeeds in this case.+ <|> branchScore++distanceOpt :: Parser DistanceMeasure+distanceOpt = option (megaReadM distanceParser) $+ long "distance" <>+ short 'd' <>+ metavar "MEASURE" <>+ help "Type of distance to calculate (available distance measures are listed below)"++summaryStatisticsSwitch :: Parser Bool+summaryStatisticsSwitch = switch $+ long "summary-statistics" <>+ short 's' <>+ help "Report summary statistics only"++normalizeSwitch :: Parser Bool+normalizeSwitch = switch $+ long "normalize" <>+ short 'n' <>+ help "Normalize trees before distance calculation; only affect distances depending on branch lengths"++newickIqTree :: Parser Bool+newickIqTree = switch $+ long "newick-iqtree"+ <> short 'i'+ <> help "Use IQ-TREE Newick format (internal node labels are branch support values)"++-- | Information about provided distance types.+distanceFooter :: String+distanceFooter = intercalate "\n"+ [ "Available distance measures:"+ , " symmetric Symmetric distance (Robinson-Foulds distance)."+ , " incompatible-split[VAL] Incompatible split distance. Collapse branches"+ , " with support less than VAL before distance calculation;"+ , " in this way, only well supported difference contribute"+ , " to the distance measure."+ , " branch-score Branch score distance."+ ]
src/Examine/Examine.hs view
@@ -19,35 +19,49 @@ ) where +import Control.Monad (unless) import Control.Monad.IO.Class import Control.Monad.Logger import Control.Monad.Trans.Class import Control.Monad.Trans.Reader import qualified Data.ByteString.Lazy.Char8 as L+import Data.List (nub, (\\)) import qualified Data.Text as T import Data.Tree+import System.IO (Handle, hPutStrLn) import Examine.Options import ELynx.Data.Tree.MeasurableTree+import ELynx.Data.Tree.NamedTree import ELynx.Data.Tree.PhyloTree+import ELynx.Data.Tree.Tree import ELynx.Import.Tree.Newick import ELynx.Tools.InputOutput-import ELynx.Tools.Logger -readTrees :: Maybe FilePath -> Examine [Tree PhyloByteStringLabel]+readTrees :: Maybe FilePath -> Examine [Tree (PhyloLabel L.ByteString)] readTrees mfp = do case mfp of Nothing -> $(logInfo) "Read tree(s) from standard input." Just fp -> $(logInfo) $ T.pack $ "Read tree(s) from file " <> fp <> "."- liftIO $ parseFileOrIOWith manyNewick mfp+ a <- lift ask+ let nw = if argsNewickIqTree a then manyNewickIqTree else manyNewick+ liftIO $ parseFileOrIOWith nw mfp +examineTree :: (Measurable a, Named a)+ => Handle -> Tree a -> IO ()+examineTree h t = do+ L.hPutStrLn h $ summarize t+ unless (null dups) (hPutStrLn h $ "Duplicate leaves: " ++ show dups)+ where lvs = map getName $ leaves t+ dups = lvs \\ nub lvs+ -- | Examine phylogenetic trees. examine :: Maybe FilePath -> Examine () examine outFn = do- ExamineArguments inFn <- lift ask+ a <- lift ask+ let inFn = argsInFile a trs <- readTrees inFn- let lsStrs = map summarize trs let outFilePath = (++ ".out") <$> outFn- logNewSection "Results."- io "results of tree analysis" (L.intercalate (L.pack "\n") lsStrs) outFilePath+ outH <- outHandle "results" outFilePath+ liftIO $ mapM_ (examineTree outH) trs
src/Examine/Options.hs view
@@ -23,17 +23,26 @@ import Options.Applicative -- | Arguments needed to examine phylogenetic trees.-newtype ExamineArguments = ExamineArguments- { inFile :: Maybe FilePath }+data ExamineArguments = ExamineArguments+ { argsInFile :: Maybe FilePath+ , argsNewickIqTree :: Bool } -- | Logger and reader data type. type Examine = LoggingT (ReaderT ExamineArguments IO) -- | Command line parser. examineArguments :: Parser ExamineArguments-examineArguments = ExamineArguments <$> optional inFileArg+examineArguments = ExamineArguments <$>+ optional inFile+ <*> newickIqTree -inFileArg :: Parser FilePath-inFileArg = strArgument $+inFile :: Parser FilePath+inFile = strArgument $ metavar "INPUT-FILE" <> help "Read trees from INPUT-FILE"++newickIqTree :: Parser Bool+newickIqTree = switch $+ long "newick-iqtree"+ <> short 'i'+ <> help "Use IQ-TREE Newick format (internal node labels are branch support values)"
+ src/Shuffle/Options.hs view
@@ -0,0 +1,68 @@+{- |+Module : Shuffle.Options+Description : Options for the connect subcommand+Copyright : (c) Dominik Schrempf 2019+License : GPL-3++Maintainer : dominik.schrempf@gmail.com+Stability : unstable+Portability : portable++Creation date: Thu Sep 19 15:02:21 2019.++-}++module Shuffle.Options+ ( ShuffleArguments (..)+ , Shuffle+ , shuffleArguments+ ) where++import Control.Monad.Logger+import Control.Monad.Trans.Reader+-- import Data.Word+import Options.Applicative++-- import ELynx.Tools.Options++-- | Arguments of shuffle command.+data ShuffleArguments = ShuffleArguments+ { newickIqTreeFlag :: Bool+ , nReplicates :: Int+ , inFile :: FilePath }+ -- -- TODO. Too complicated at the moment with MWC, use MonadRandom, but then+ -- -- seed option would not match other commands, so leave the option out for+ -- -- now. Maybe, when working on this, try using splitmix consistently.+ --+ -- -- | Seed of NRG, random if 'Nothing'.+ -- , argsSeed :: Maybe [Word32] }++-- | Logger and reader data type.+type Shuffle = LoggingT (ReaderT ShuffleArguments IO)++-- | Parse arguments of shuffle command.+shuffleArguments :: Parser ShuffleArguments+shuffleArguments = ShuffleArguments+ <$> newickIqTree+ <*> n+ <*> file+ -- <*> seedOpt++newickIqTree :: Parser Bool+newickIqTree = switch $+ long "newick-iqtree"+ <> short 'i'+ <> help "Use IQ-TREE Newick format (internal node labels are branch support values)"++n :: Parser Int+n = option auto $+ long "replicates"+ <> short 'n'+ <> metavar "N"+ <> value 1+ <> help "Number of trees to generate"++file :: Parser FilePath+file = strArgument $+ metavar "TREE-FILE"+ <> help "File containing a Newick tree"
+ src/Shuffle/Shuffle.hs view
@@ -0,0 +1,142 @@+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE TemplateHaskell #-}++{- |+Module : Shuffle.Shuffle+Description : Shuffle a phylogeny.+Copyright : (c) Dominik Schrempf 2019+License : GPL-3++Maintainer : dominik.schrempf@gmail.com+Stability : unstable+Portability : portable++Creation date: Thu Sep 19 15:01:52 2019.++The coalescent times are unaffected. The topology and the leaf order is+shuffled. Branch support values are ignored and lost.++-}++module Shuffle.Shuffle+ ( shuffleCmd+ ) where++import Control.Monad.IO.Class+import Control.Monad.Logger+import Control.Monad.Random+import Control.Monad.Trans.Reader+import Data.Array (elems)+import Data.Array.ST (newListArray, readArray,+ runSTArray, writeArray)+import qualified Data.ByteString.Lazy.Char8 as L+import Data.List (filter)+import Data.Maybe (isNothing)+import Data.Tree (Tree (Node), flatten,+ rootLabel)+import System.IO++import Shuffle.Options++import ELynx.Data.Tree.MeasurableTree (distancesOriginLeaves, height,+ rootHeight)+import ELynx.Data.Tree.NamedTree (getName)+import ELynx.Data.Tree.PhyloTree (PhyloLabel, brLen)+import ELynx.Data.Tree.Tree (leaves)+import ELynx.Export.Tree.Newick (toNewick)+import ELynx.Import.Tree.Newick (oneNewick)+import ELynx.Simulate.PointProcess (PointProcess (PointProcess),+ toReconstructedTree)+import ELynx.Tools.Definitions (eps)+import ELynx.Tools.InputOutput (outHandle, parseFileWith)+import ELynx.Tools.Text (fromBs, tShow)++-- | Shuffle a tree. Get all coalescent times, shuffle them. Get all leaves,+-- shuffle them. Connect the shuffled leaves with the shuffled coalescent times.+-- The shuffled tree has a new topology while keeping the same set of coalescent+-- times and leaves.+shuffleCmd :: Maybe FilePath -> Shuffle ()+shuffleCmd outFile = do+ a <- lift ask++ let outFn = (++ ".tree") <$> outFile+ h <- outHandle "results" outFn++ t <- liftIO $ parseFileWith oneNewick (inFile a)+ $(logInfo) "Input tree:"+ $(logInfo) $ fromBs $ toNewick t++ -- Check if all branches have a given length. However, the length of the stem is not important.+ let r = rootLabel t+ r' = r {brLen = Just 0}+ t' = t {rootLabel = r'}+ when (isNothing $ traverse brLen t')+ (do+ $(logDebug) $ tShow t'+ error "Not all branches have a given length.")++ -- Check if tree is ultrametric enough.+ let dh = sum $ map (height t -) (distancesOriginLeaves t)+ $(logDebug) $ "Distance in branch length to being ultrametric: " <> tShow dh+ when (dh > 2e-4)+ (error "Tree is not ultrametric.")+ when (dh > eps && dh < 2e-4) $+ $(logInfo) "Tree is nearly ultrametric, ignore branch length differences smaller than 2e-4."+ when (dh < eps) $+ $(logInfo) "Tree is ultrametric."++ let cs = filter (>0) $ flatten $ mapTree rootHeight t+ ls = map getName $ leaves t+ $(logDebug) $ "Number of coalescent times: " <> tShow (length cs)+ $(logDebug) $ "Number of leaves: " <> tShow (length ls)+ $(logDebug) "The coalescent times are: "+ $(logDebug) $ tShow cs++ ts <- liftIO $ shuffle (nReplicates a) (height t) cs ls+ liftIO $ L.hPutStr h $ L.unlines $ map toNewick ts++ liftIO $ hClose h++shuffle :: MonadRandom m+ => Int -- How many?+ -> Double -- Stem length.+ -> [Double] -- Coalescent times.+ -> [L.ByteString] -- Leave names.+ -> m [Tree (PhyloLabel L.ByteString)]+shuffle n o cs ls = do+ css <- grabble cs n (length cs)+ lss <- grabble ls n (length ls)+ return [toReconstructedTree "" (PointProcess names times o) | (times, names) <- zip css lss]++-- TODO: This seems to be 'extend' from comonad!+mapTree :: (Tree a -> b) -> Tree a -> Tree b+mapTree f t@(Node _ cs) = Node (f t) (map (mapTree f) cs)++-- | From https://wiki.haskell.org/Random_shuffle.+--+-- @grabble xs m n@ is /O(m*n')/, where @n' = min n (length xs)@. Choose @n@+-- elements from @xs@, without replacement, and that @m@ times.+grabble :: MonadRandom m => [a] -> Int -> Int -> m [[a]]+grabble xs m n = do+ swapss <- replicateM m $ forM [0 .. min (maxIx - 1) n] $ \i -> do+ j <- getRandomR (i, maxIx)+ return (i, j)+ return $ map (take n . swapElems xs) swapss+ where+ maxIx = length xs - 1++-- grabbleOnce :: MonadRandom m => [a] -> Int -> m [a]+-- grabbleOnce xs n = head `liftM` grabble xs 1 n++swapElems :: [a] -> [(Int, Int)] -> [a]+swapElems xs swaps = elems $ runSTArray (do+ arr <- newListArray (0, maxIx) xs+ mapM_ (swap arr) swaps+ return arr)+ where+ maxIx = length xs - 1+ swap arr (i,j) = do+ vi <- readArray arr i+ vj <- readArray arr j+ writeArray arr i vj+ writeArray arr j vi
src/Simulate/Options.hs view
@@ -17,7 +17,6 @@ , Simulate , simulateArguments , reportSimulateArguments- , simulateDesc , simulateFooter ) where @@ -151,9 +150,8 @@ <> showDefault <> help "Only output number of children for each branch" --- | Simulation is a little more involved, so we provide a description.-simulateDesc :: String-simulateDesc = "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"+citation :: String+citation = "Gernhard, T. (2008). The conditioned reconstructed process. Journal of Theoretical Biology, 253(4), 769–778. http://doi.org/10.1016/j.jtbi.2008.04.005" -- | And a footer. simulateFooter :: String@@ -161,4 +159,4 @@ [ "Height of Trees: if no tree height is given, the heights will be randomly drawn from the expected distribution given the number of leaves, the birth and the death rate." , "Summary statistics only: only print (NumberOfExtantChildren BranchLength) pairs for each branch of each tree. The trees are separated by a newline character." , "Sub-sampling: simulate one big tree with n'=round(n/rho), n'>=n, leaves, and randomly sample sub-trees with n leaves. Hence, with rho=1.0, the same tree is reported over and over again."- ]+ , citation ]
src/Simulate/Simulate.hs view
@@ -47,7 +47,7 @@ import Simulate.Options import ELynx.Data.Tree.MeasurableTree-import ELynx.Data.Tree.PhyloTree (PhyloIntLabel)+import ELynx.Data.Tree.PhyloTree (PhyloLabel) import ELynx.Data.Tree.SumStat (formatNChildSumStat, toNChildSumStat) import ELynx.Data.Tree.Tree@@ -78,9 +78,9 @@ else parMap rpar toNewick trs let outFile' = (++ ".tree") <$> outFile let res = L.unlines ls- io "simulated trees" res outFile'+ out "simulated trees" res outFile' -simulateNTreesConcurrently :: Int -> Simulate [Tree PhyloIntLabel]+simulateNTreesConcurrently :: Int -> Simulate [Tree (PhyloLabel Int)] simulateNTreesConcurrently c = do (SimulateArguments nT nL h cM l m r _ _ s) <- lift ask let l' = l * r@@ -93,7 +93,7 @@ (zip chunks gs) return $ concat trss -simulateAndSubSampleNTreesConcurrently :: Int -> Simulate [Tree PhyloIntLabel]+simulateAndSubSampleNTreesConcurrently :: Int -> Simulate [Tree (PhyloLabel Int)] simulateAndSubSampleNTreesConcurrently c = do (SimulateArguments nT nL h cM l m r _ _ s) <- lift ask let nLeavesBigTree = (round $ fromIntegral nL / r) :: Int
tlynx.cabal view
@@ -1,6 +1,6 @@ cabal-version: 1.12 name: tlynx-version: 0.0.1+version: 0.1.0 license: GPL-3 license-file: LICENSE copyright: Dominik Schrempf (2019)@@ -25,25 +25,33 @@ exposed-modules: Compare.Compare Compare.Options+ Connect.Connect+ Connect.Options+ Distance.Distance+ Distance.Options Examine.Examine Examine.Options+ Shuffle.Options+ Shuffle.Shuffle Simulate.Options Simulate.Simulate hs-source-dirs: src other-modules: Paths_tlynx default-language: Haskell2010- ghc-options: -Wall+ ghc-options: -Wall -O2 build-depends:+ MonadRandom >=0.5.1.1 && <0.6, QuickCheck >=2.13.2 && <2.14,+ array >=0.5.3.0 && <0.6, base >=4.12.0.0 && <4.13, bytestring >=0.10.8.2 && <0.11, containers >=0.6.0.1 && <0.7,- elynx-seq >=0.0.1 && <0.1,- elynx-tools >=0.0.1 && <0.1,- elynx-tree >=0.0.1 && <0.1,+ elynx-seq >=0.1.0 && <0.2,+ elynx-tools >=0.1.0 && <0.2,+ elynx-tree >=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, monad-logger >=0.3.30 && <0.4, mwc-random >=0.14.0.0 && <0.15,@@ -63,16 +71,19 @@ Options Paths_tlynx default-language: Haskell2010- ghc-options: -Wall -eventlog -threaded -rtsopts -with-rtsopts=-N+ ghc-options: -Wall -O2 -eventlog -threaded -rtsopts+ -with-rtsopts=-N build-depends:+ MonadRandom >=0.5.1.1 && <0.6, QuickCheck >=2.13.2 && <2.14,+ array >=0.5.3.0 && <0.6, 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-tree >=0.0.1 && <0.1,+ elynx-tools >=0.1.0 && <0.2,+ elynx-tree >=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, monad-logger >=0.3.30 && <0.4, mwc-random >=0.14.0.0 && <0.15,