diff --git a/Bio/Phylogeny/PhyBin.hs b/Bio/Phylogeny/PhyBin.hs --- a/Bio/Phylogeny/PhyBin.hs +++ b/Bio/Phylogeny/PhyBin.hs @@ -5,7 +5,7 @@ {-# OPTIONS_GHC -fwarn-unused-imports #-} -- | This module contains the code that does the tree normalization and binning. --- It's the heart of the prgoram. +-- It's the heart of the program. module Bio.Phylogeny.PhyBin ( driver, binthem, normalize, annotateWLabLists, unitTests, acquireTreeFiles, @@ -15,12 +15,14 @@ import qualified Data.Foldable as F import Data.Function (on) import Data.List (delete, minimumBy, sortBy, foldl1', foldl', intersperse, isPrefixOf) -import Data.Maybe (fromMaybe) +import Data.Maybe (fromMaybe, catMaybes) import Data.Either (partitionEithers) import Data.Time.Clock -import qualified Data.ByteString.Lazy.Char8 as B +import Data.Tuple (swap) +import qualified Data.ByteString.Char8 as B import qualified Data.Map as M import qualified Data.Set as S +import qualified Data.Text.Lazy.IO as T import qualified Data.Vector as V import qualified Data.Vector.Unboxed as U import Control.Monad (forM, forM_, filterM, when, unless) @@ -31,23 +33,27 @@ import System.FilePath (combine) import System.Directory (doesFileExist, doesDirectoryExist, createDirectoryIfMissing, getDirectoryContents, getCurrentDirectory) -import System.IO (openFile, hClose, IOMode(..), stdout) +import System.IO (openFile, hClose, IOMode(..), stdout, withFile) import System.Info (os) import System.Process (system) import System.Exit (ExitCode(..)) import Test.HUnit ((~:),(~=?),Test,test) import qualified Data.Clustering.Hierarchical as C +import qualified Data.GraphViz as Gv -- For vizualization: import Text.PrettyPrint.HughesPJClass hiding (char, Style) +import Data.GraphViz.Printing (renderDot) import Bio.Phylogeny.PhyBin.CoreTypes import Bio.Phylogeny.PhyBin.Parser (parseNewick, parseNewicks) -import Bio.Phylogeny.PhyBin.PreProcessor (collapseBranchLenThresh) +import Bio.Phylogeny.PhyBin.PreProcessor (collapseBranchLenThresh, + collapseBranchBootStrapThresh, pruneTreeLeaves) import Bio.Phylogeny.PhyBin.Visualize (dotToPDF, dotNewickTree, viewNewickTree, dotDendrogram) import Bio.Phylogeny.PhyBin.RFDistance import Bio.Phylogeny.PhyBin.Binning import Bio.Phylogeny.PhyBin.Util + import Debug.Trace ---------------------------------------------------------------------------------------------------- @@ -75,8 +81,9 @@ -- | Driver to put all the pieces together (parse, normalize, bin) driver :: PhyBinConfig -> IO () driver cfg@PBC{ verbose, num_taxa, name_hack, output_dir, inputs=files, - do_graph, branch_collapse_thresh, highlights, - dist_thresh, clust_mode, use_hashrf, print_rfmatrix } = + do_graph, branch_collapse_thresh, bootstrap_collapse_thresh, highlights, + dist_thresh, clust_mode, rfmode, preprune_labels, + print_rfmatrix } = -- Unused: do_draw do -------------------------------------------------------------------------------- @@ -115,48 +122,85 @@ let num_files = length goodFiles bstrs <- mapM B.readFile goodFiles let (labelTab, fulltrees) = parseNewicks name_hack (zip files bstrs) - - highlightTrs <- retrieveHighlights name_hack labelTab highlights + nameToLabel = M.fromList $ map swap $ M.toList labelTab - putStrLn$ "\nTotal unique taxa ("++ show (M.size labelTab) ++"):\n "++ - (unwords$ M.elems labelTab) --- show (nest 2 $ sep $ map text $ M.elems labelTab) + highlightTrs <- retrieveHighlights name_hack labelTab highlights + + let allTaxaLs = M.elems labelTab + totalTaxa = length allTaxaLs + putStrLn$ "\nTotal unique taxa ("++ show (M.size labelTab) ++"):\n "++ (unwords allTaxaLs) + + expected_num_taxa <- + case preprune_labels of + Nothing -> return totalTaxa + Just ls -> return (length ls) + case rfmode of + TolerantNaive -> return () + HashRF -> putStrLn$ "Note: defaulting to expecting ALL "++show expected_num_taxa++" to be present.." -------------------------------------------------------------------------------- + case bootstrap_collapse_thresh of + Just thr -> putStrLn$" !+ Collapsing branches of bootstrap value less than "++show thr + Nothing -> return () + case branch_collapse_thresh of Just thr -> putStrLn$" !+ Collapsing branches of length less than "++show thr Nothing -> return () - let do_one :: FullTree DefDecor -> IO (Int, [FullTree DefDecor], [(Int, String)]) + let do_one :: FullTree DefDecor -> IO (Int, [FullTree DefDecor], Maybe (Int, String)) do_one (FullTree treename lblAcc parsed) = do - let - pruned = case branch_collapse_thresh of - Nothing -> parsed - Just thr -> collapseBranchLenThresh thr parsed - numL = numLeaves pruned + let + pruned0 = parsed + -- Important: MUST collapse bootstrap first if we're doing both: + pruned1 = case bootstrap_collapse_thresh of + Nothing -> pruned0 + Just thr -> collapseBranchBootStrapThresh thr pruned0 + pruned2 = case branch_collapse_thresh of + Nothing -> pruned1 + Just thr -> collapseBranchLenThresh thr pruned1 + pruned3 = case preprune_labels of + Nothing -> Just pruned2 + Just lst -> pruneTreeLeaves (S.fromList (map lkup lst)) pruned2 + + lkup trnm = + case M.lookup trnm nameToLabel of + Nothing -> error$ "Tree name "++show trnm++" did not occur ANYWHERE in the input set." + Just lab -> lab + -- TEMPTOGGLE -- when False $ do putStrLn$ "DRAWING TREE" -- viewNewickTree "Annotated" (FullTree file lblAcc' annot) -- viewNewickTree "Normalized" (FullTree file lblAcc' normal) -- putStrLn$ "WEIGHTS OF NORMALIZED' CHILDREN: "++ -- show (map get_weight$ get_children normal) - - if numL /= num_taxa - then do --putStrLn$ "\n WARNING: file contained an empty or single-node tree: "++ show file - when verbose$ putStrLn$ "\n WARNING: tree contained unexpected number of leaves (" - ++ show numL ++"): "++ treename - return (0, [], [(numL, treename)]) - else do - when verbose$ putStr "." - return$ (numL, [FullTree treename lblAcc pruned], []) - + when verbose$ putStr "." + case pruned3 of + Nothing -> do when verbose$ putStrLn$ "\n WARNING: tree empty after preprocessing, discarding: "++ treename + return (0, [], Just (0, treename)) + Just pruned3' -> + let numL = numLeaves pruned3' + looksOK = return (numL, [FullTree treename lblAcc pruned3'], Nothing) + discardIt = return (0, [], Just (numL, treename)) + in case rfmode of + TolerantNaive -> looksOK + HashRF | numL == expected_num_taxa -> looksOK + | otherwise -> do + when verbose$ + putStrLn$ "\n WARNING: tree contained unexpected number of leaves (" + ++ show numL ++"): "++ treename + discardIt + results <- mapM do_one fulltrees - let (counts::[Int], validtreess, pairs::[[(Int, String)]]) = unzip3 results + let (counts::[Int], validtreess, pairs:: [Maybe (Int, String)]) = unzip3 results let validtrees = concat validtreess - warnings2 = concat pairs - + warnings2 = catMaybes pairs + putStrLn$ "\nNumber of input tree files: " ++ show num_files + mapM_ (print . displayStrippedTree) validtrees -- Debugging. + case preprune_labels of + Nothing -> return () + Just lst -> putStrLn$ "PRUNING trees to just these taxa: "++show lst when (length warnings2 > 0) $ putStrLn$ "Number of bad/unreadable input tree files: " ++ show (length warnings2) putStrLn$ "Number of VALID trees (correct # of leaves/taxa): " ++ show (length validtrees) @@ -179,24 +223,37 @@ M.toList x return (x,binlist,[]) ClusterThem{linkage} -> do - (mat, dendro) <- doCluster use_hashrf num_taxa linkage validtrees + (mat, dendro) <- doCluster (rfmode==HashRF) expected_num_taxa linkage validtrees + -------------------- when print_rfmatrix $ printDistMat stdout mat - hnd <- openFile (combine output_dir ("distance_matrix.txt")) WriteMode - printDistMat hnd mat - hClose hnd + withFile (combine output_dir ("distance_matrix.txt")) WriteMode $ \ hnd -> + printDistMat hnd mat -------------------- writeFile (combine output_dir ("dendrogram.txt")) (show$ fmap treename dendro) putStrLn " [finished] Wrote full dendrogram to file dendrogram.txt" + sanityCheck dendro + -------------------- let plotIt mnameMap = if True -- do_graph - then async (do + then async (do t0 <- getCurrentTime let dot = dotDendrogram cfg "dendrogram" 1.0 dendro mnameMap highlightTrs - _ <- dotToPDF dot (combine output_dir "dendrogram.pdf") - t1 <- getCurrentTime - putStrLn$ " [finished] Wrote dendrogram diagram to file dendrogram.pdf (" - ++show(diffUTCTime t1 t0)++")") + -- Be wary of cycles, something appears to be buggy: + putStrLn$ " [async] writing dendrogram as a graph to dendrogram.dot" + -- writeFile (combine output_dir "dendrogram.dot") + -- (show $ renderDot $ Gv.toDot dot) + mtxt <- safePrintDendro dot + case mtxt of + Nothing -> + putStrLn "WARNING: because we couldn't print it, we're not drawing the dendrogram either." + Just txt -> do + T.writeFile (combine output_dir "dendrogram.dot") txt + putStrLn$ " [async] Next, plot dendrogram.pdf" + _ <- dotToPDF dot (combine output_dir "dendrogram.pdf") + t1 <- getCurrentTime + putStrLn$ " [finished] Writing dendrogram diagram (" + ++show(diffUTCTime t1 t0)++")") else async (return ()) case dist_thresh of Nothing -> do a <- plotIt Nothing @@ -247,12 +304,12 @@ async2 <- case clust_mode of BinThem -> outputBins binlist output_dir do_graph - ClusterThem{} -> outputClusters num_taxa binlist output_dir do_graph + ClusterThem{} -> outputClusters expected_num_taxa binlist output_dir do_graph -- Wait on parallel tasks: - putStrLn$ "Waiting for asynchronous tasks to finish..." + putStrLn$ "Waiting for "++show (length$ async2:asyncs)++" asynchronous tasks to finish..." mapM_ wait (async2:asyncs) - putStrLn$ "Finished." + putStrLn$ "Phybin completed." -------------------------------------------------------------------------------- -- End driver -------------------------------------------------------------------------------- @@ -272,15 +329,15 @@ binthem validtrees return (classes) -doCluster :: Bool -> Int -> C.Linkage -> [FullTree a] -> IO (DistanceMatrix, C.Dendrogram (FullTree a)) -doCluster use_hashrf num_taxa linkage validtrees = do +doCluster :: Bool -> Int -> C.Linkage -> [FullTree DefDecor] -> IO (DistanceMatrix, C.Dendrogram (FullTree DefDecor)) +doCluster use_hashrf expected_num_taxa linkage validtrees = do t0 <- getCurrentTime when use_hashrf$ putStrLn " Using HashRF-style algorithm..." let nwtrees = map nwtree validtrees numtrees = length validtrees mat = if use_hashrf - then hashRF num_taxa nwtrees - else fst (distanceMatrix nwtrees) + then hashRF expected_num_taxa nwtrees + else fst (naiveDistMatrix nwtrees) ixtrees = zip [0..] validtrees dist (i,t1) (j,t2) | j == i = 0 -- | i == numtrees-1 = 0 @@ -383,7 +440,8 @@ async $ do mapM_ wait asyncs putStrLn$ " [finished] Wrote visual representations of consensus trees to "++filePrefix++"_.pdf" - else async (return ()) + else do putStrLn$ "NOT creating processes to build per-cluster .pdf visualizations. (Not asked to.)" + async (return ()) outputBins :: [(Int, OneCluster StandardDecor)] -> String -> Bool -> IO (Async ()) outputBins binlist output_dir do_graph = do @@ -470,7 +528,7 @@ avg ls = sum ls / fromIntegral (length ls) --- | Parse trees in addition to the main inputs (for --highlight). +-- | Parse extra trees in addition to the main inputs (for --highlight). retrieveHighlights :: (String->String) -> LabelTable -> [FilePath] -> IO [[NewickTree ()]] retrieveHighlights name_hack labelTab ls = mapM parseHighlight ls @@ -486,6 +544,10 @@ -- | Create a predicate that tests trees for consistency with the set of --highlight -- (consensus) trees. +-- +-- Note, tree consistency is not the same as an exact match. It's +-- like (<=) rather than (==). All trees are consistent with the +-- "star topology". matchAnyHighlight :: [[NewickTree ()]] -> NewickTree () -> Bool -- matchAnyHighlight :: [[NewickTree ()]] -> NewickTree () -> Maybe Int -- If there is a match, return the index of the highlight that matched. diff --git a/Bio/Phylogeny/PhyBin/Binning.hs b/Bio/Phylogeny/PhyBin/Binning.hs --- a/Bio/Phylogeny/PhyBin/Binning.hs +++ b/Bio/Phylogeny/PhyBin/Binning.hs @@ -15,32 +15,19 @@ ) where -import qualified Data.Foldable as F import Data.Function (on) import Data.List (delete, minimumBy, sortBy, insertBy, intersperse, sort) -import Data.Maybe (fromMaybe, catMaybes) -import qualified Data.ByteString.Lazy.Char8 as B +import qualified Data.ByteString.Char8 as B import qualified Data.Map as M -import qualified Data.Set as S -import Control.Monad (forM, forM_, filterM, when, unless) -import Control.Exception (evaluate) import Control.Applicative ((<$>),(<*>)) import Control.Concurrent (Chan) -import System.FilePath (combine) -import System.Directory (doesFileExist, doesDirectoryExist, - getDirectoryContents, getCurrentDirectory) -import System.IO (openFile, hClose, IOMode(ReadMode)) -import System.Exit (ExitCode(..)) import Test.HUnit ((~:),(~=?),Test,test) -- For vizualization: import Text.PrettyPrint.HughesPJClass hiding (char, Style) import Bio.Phylogeny.PhyBin.CoreTypes import Bio.Phylogeny.PhyBin.Parser (parseNewicks, parseNewick) -import Bio.Phylogeny.PhyBin.PreProcessor (collapseBranches) -import Bio.Phylogeny.PhyBin.Visualize (dotToPDF, dotNewickTree, viewNewickTree) -import Bio.Phylogeny.PhyBin.RFDistance -import Bio.Phylogeny.PhyBin.Util +import Bio.Phylogeny.PhyBin.Visualize (dotNewickTree, viewNewickTree) -- Turn on for extra invariant checking: debug :: Bool @@ -407,8 +394,8 @@ -- Make sure that all of these normalize to the same thing. , "normalize1" ~: "(C, D, E, (A, B))" ~=? show (displayDefaultTree $ deAnnotate$ norm "(A,(C,D,E),B);") - , "normalize2" ~: "(C, D, E, (A, B))" ~=? show (pPrint$ norm "((C,D,E),B,A);") - , "normalize2" ~: "(C, D, E, (A, B))" ~=? show (pPrint$ norm "(D,E,C,(B,A));") + , "normalize2A" ~: "(C, D, E, (A, B))" ~=? show (pPrint$ norm "((C,D,E),B,A);") + , "normalize2B" ~: "(C, D, E, (A, B))" ~=? show (pPrint$ norm "(D,E,C,(B,A));") -- Here's an example from the rhizobia datasetsthat that caused my branch-sorting to fail. , "normalize3" ~: "(((BB, BJ)), (MB, ML), (RE, (SD, SM)))" diff --git a/Bio/Phylogeny/PhyBin/CoreTypes.hs b/Bio/Phylogeny/PhyBin/CoreTypes.hs --- a/Bio/Phylogeny/PhyBin/CoreTypes.hs +++ b/Bio/Phylogeny/PhyBin/CoreTypes.hs @@ -9,7 +9,7 @@ -- * Tree and tree decoration types NewickTree(..), DefDecor, StandardDecor(..), AnnotatedTree, FullTree(..), - ClustMode(..), TreeName, + ClustMode(..), TreeName, NumTaxa(..), -- * Tree operations displayDefaultTree, displayStrippedTree, @@ -21,7 +21,8 @@ avg_branchlen, get_bootstraps, -- * Command line config options - PhyBinConfig(..), default_phybin_config, + PhyBinConfig(..), default_phybin_config, + WhichRFMode(..), -- * General helpers Label, LabelTable, @@ -32,6 +33,7 @@ where import qualified Data.Map as M +import qualified Data.Set as S import Data.Foldable (Foldable(..)) import Data.Maybe (maybeToList) import Data.Monoid (mappend, mconcat) @@ -105,6 +107,7 @@ pPrint mp = pPrint (M.toList mp) -- | Display a tree WITH the bootstrap and branch lengths. +-- This prints in NEWICK format. displayDefaultTree :: FullTree DefDecor -> Doc displayDefaultTree orig = loop tr <> ";" where @@ -117,6 +120,8 @@ Just val -> base <> text ":[" <> text (show val) <> text "]" where base = parens$ sep$ map_but_last (<>text",") $ map loop ls +-- | The same, except with no bootstrap or branch lengths. Any tree annotations +-- ignored. displayStrippedTree :: FullTree a -> Doc displayStrippedTree orig = loop tr <> ";" where @@ -128,10 +133,7 @@ -- Labels ---------------------------------------------------------------------------------------------------- - --- Experimental: toggle this to change the representation of labels: --- Alas I always have problems with the interned string libs (e.g. segfaults)... [2012.11.20] ----------------------------------------- +-- | Labels are inexpensive unique integers. The table is necessary for converting them back. type Label = Int -- | Map labels back onto meaningful names. @@ -208,7 +210,7 @@ -- | Due to the number of configuration options for the driver, we pack them into a record. data PhyBinConfig = PBC { verbose :: Bool - , num_taxa :: Int + , num_taxa :: NumTaxa , name_hack :: String -> String , output_dir :: String , inputs :: [String] @@ -218,17 +220,32 @@ , highlights :: [FilePath] -- [FullTree ()] , show_trees_in_dendro :: Bool , show_interior_consensus :: Bool - , use_hashrf :: Bool + , rfmode :: WhichRFMode + , preprune_labels :: Maybe [String] , print_rfmatrix :: Bool , dist_thresh :: Maybe Int , branch_collapse_thresh :: Maybe Double -- ^ Branches less than this length are collapsed. + , bootstrap_collapse_thresh :: Maybe Int + -- ^ BootStrap values less than this result in the intermediate node being collapsed. } +-- | Supported modes for computing RFDistance. +data WhichRFMode = HashRF | TolerantNaive + deriving (Show, Eq, Ord) + +-- | How many taxa should we expect in the incoming dataset? +data NumTaxa = Expected Int -- ^ Supplied by the user. Committed. + | Unknown -- ^ In the future we may automatically pick a behavior. Now this one is usually an error. + | Variable -- ^ Explicitly ignore this setting in favor of comparing all trees + -- (even if some are missing taxa). This only works with certain modes. + deriving (Show, Read, Eq) + -- | The default phybin configuration. default_phybin_config :: PhyBinConfig default_phybin_config = PBC { verbose = False - , num_taxa = error "must be able to determine the number of taxa expected in the dataset. (Supply it manually with -n.)" +-- , num_taxa = error "must be able to determine the number of taxa expected in the dataset. (Supply it manually with -n.)" + , num_taxa = Unknown , name_hack = id -- Default, no transformation of leaf-labels , output_dir = "./phybin_out/" , inputs = [] @@ -236,17 +253,15 @@ , do_draw = False -- , clust_mode = BinThem , clust_mode = ClusterThem C.UPGMA -#ifdef USE_HASHRF - , use_hashrf = True -#else - , use_hashrf = False -#endif + , rfmode = HashRF + , preprune_labels = Nothing , highlights = [] , show_trees_in_dendro = False , show_interior_consensus = False , print_rfmatrix = False , dist_thresh = Nothing , branch_collapse_thresh = Nothing + , bootstrap_collapse_thresh = Nothing } data ClustMode = BinThem | ClusterThem { linkage :: C.Linkage } diff --git a/Bio/Phylogeny/PhyBin/Parser.hs b/Bio/Phylogeny/PhyBin/Parser.hs --- a/Bio/Phylogeny/PhyBin/Parser.hs +++ b/Bio/Phylogeny/PhyBin/Parser.hs @@ -6,13 +6,12 @@ (newick_parser, parseNewick, parseNewicks, parseNewickFiles, unitTests) where import Control.Exception (evaluate, handle, SomeException) -import qualified Data.ByteString.Lazy.Char8 as B +import qualified Data.ByteString.Char8 as B import Data.Char (isSpace) import Data.Map as M -import Data.Set as S import Data.List as L import Text.Parsec -import Text.Parsec.ByteString.Lazy +import Text.Parsec.ByteString import Test.HUnit ((~:),(~=?),Test,test,assertFailure) import System.FilePath (takeBaseName) @@ -24,11 +23,16 @@ -- | Parse a bytestring into a NewickTree with branch lengths. The -- first argument is file from which the data came and is just for -- better error messages. +-- +-- If the single bytestring contains more than one tree, then a number is appended to +-- the tree names. parseNewick :: LabelTable -> NameHack -> String -> B.ByteString -> (LabelTable, [NewickTree DefDecor]) -parseNewick tbl0 name_hack file input = - extractLabelTable tbl0 $ - runB file (many1$ newick_parser name_hack) $ - B.filter (not . isSpace) input +parseNewick tbl0 name_hack file input = (lbls,trs) + where + (lbls,trs) = + extractLabelTable tbl0 $ + runB file (many1$ newick_parser name_hack) $ + B.filter (not . isSpace) input -- treeFiles <- acquireTreeFiles files -- let fn f = do raw <- B.readFile f @@ -45,6 +49,7 @@ -- | Parse a list of trees, starting with an empty map of labels and accumulating a final map. parseNewickFiles :: NameHack -> [String] -> IO (LabelTable, [FullTree DefDecor]) parseNewickFiles fn nms = do + -- WARNING: Lazy-IO here. We need to be careful about leaking file descriptors. bss <- mapM B.readFile nms return $ parseNewicks fn (zip nms bss) @@ -52,8 +57,17 @@ parseNewicks :: NameHack -> [(String,B.ByteString)] -> (LabelTable, [FullTree DefDecor]) parseNewicks name_hack pairs = (labtbl, fullTrs) where - fullTrs = [ FullTree (takeBaseName file) labtbl tr - | (file,tr) <- trs ] + fullTrs = [ FullTree (tweakName file ind) labtbl tr + | (file,tr) <- trs + | ind <- [(0::Int)..] ] + tweakName file ind = -- Here we do the renaming/numbering business: + if addSuffix + then (takeBaseName file)++"_"++show ind + else (takeBaseName file) + addSuffix = case trs of + [] -> False -- Should this be an error? + [_] -> False + _ -> True (labtbl, trs) = P.foldr fn (M.empty,[]) pairs fn (file,bstr) (!acc,!ls) = let (acc',trs) = parseNewick acc name_hack file bstr diff --git a/Bio/Phylogeny/PhyBin/PreProcessor.hs b/Bio/Phylogeny/PhyBin/PreProcessor.hs new file mode 100644 --- /dev/null +++ b/Bio/Phylogeny/PhyBin/PreProcessor.hs @@ -0,0 +1,95 @@ +{-# LANGUAGE ScopedTypeVariables #-} + +-- | This is a preprocessor which can (optionally) be invoked to +-- collapse short branch lengths. + +module Bio.Phylogeny.PhyBin.PreProcessor + ( collapseBranches, + collapseBranchLenThresh, collapseBranchBootStrapThresh, + pruneTreeLeaves + ) + where + +import qualified Data.Set as S +import Data.Maybe (catMaybes) +import Bio.Phylogeny.PhyBin.CoreTypes + + +-- | Prune the leaves of the tree to only those leaves in the provided set. +-- +-- If ALL leaves are pruned from the set, this function returns nothing. +pruneTreeLeaves :: S.Set Label -> NewickTree DefDecor -> Maybe (NewickTree DefDecor) +pruneTreeLeaves set tr = loop tr + where + loop orig@(NTLeaf _ lab) + | S.member lab set = Just orig + | otherwise = Nothing + loop (NTInterior dec@(_,blen) ls) = + case catMaybes $ map loop ls of + [] -> Nothing + -- Here we ELIMINATE the intermediate node, but we add in its branch length: + [one] -> Just$ fmap (\(bs,bl) -> (bs,blen + bl)) one + ls' -> Just (NTInterior dec ls') + + +-- | Removes branches that do not meet a predicate, leaving a shallower, "bushier" +-- tree. This does NOT change the set of leaves (taxa), it only removes interior +-- nodes. +-- +-- `collapseBranches pred collapser tree` uses `pred` to test the meta-data to see +-- if collapsing the intermediate node below the branch is necessary, and if it is, +-- it uses `collapser` to reduce all the metadata for the collapsed branches into a +-- single piece of metadata. +collapseBranches :: forall a . (a -> Bool) -> (a -> a -> a) -> NewickTree a -> NewickTree a +collapseBranches isCollapsable collapse origtr = final + where + (_,_, final) = loop origtr + + -- This loop returns: + -- (1) a list of leaf "floaters" that can still move upwards, + -- (2) immovable subtrees that can't + -- (3) a final node IF the result is the root. + loop :: NewickTree a -> ([(a,Label)], [NewickTree a], NewickTree a) + loop lf@(NTLeaf dec lb) | isCollapsable dec = ([(dec,lb)], [], lf) + | otherwise = ([], [lf], lf) + loop (NTInterior dec children) = + let (floats, anchors, _) = unzip3 $ map loop children + thenode = NTInterior dec $ concat anchors ++ + map (uncurry NTLeaf) (concat floats) + in + if isCollapsable dec then + -- If we are collapsable we keep floating upwards. + -- We combine our metadata (on the left) into the floatees: + (map (\ (a,l) -> (collapse dec a, l)) + (concat floats), + concat anchors, + thenode) + else + -- Otherwise this is the end of the road for these floaters: + ([], [thenode], thenode) + + +-- | A common configuration. Collapse branches based on a length threshold. +collapseBranchLenThresh :: Double -> NewickTree DefDecor -> NewickTree DefDecor +-- collapseBranchLenThresh :: HasBranchLen a => Double -> NewickTree a -> NewickTree a +collapseBranchLenThresh thr tr = + collapseBranches ((< thr) . getBranchLen) collapser tr + where + -- We REMOVE BootStraps as part of this process, they are not well-defined after this point. + collapser _intermediate@(_bt1, len1) _floatee@(_bt2, len2) = + (Nothing, len1 + len2) + +-- | A common configuration. Collapse branches based on bootstrap values. +collapseBranchBootStrapThresh :: Int -> NewickTree DefDecor -> NewickTree DefDecor +-- collapseBranchLenThresh :: HasBranchLen a => Double -> NewickTree a -> NewickTree a +collapseBranchBootStrapThresh thr tr = + collapseBranches ((< thr) . getBoot) collapser tr + where + getBoot (Nothing,_) = error$"collapseBranchBootStrapThresh: bootstrap value missing on tree node!\n"++ + "They must be present if --minbootstrap is used." + getBoot (Just boot,_) = boot + -- This had better happen BEFORE branch-length based collapsing is done: + collapser (_,len1) (_,len2) = (Nothing, len1+len2) + + + diff --git a/Bio/Phylogeny/PhyBin/RFDistance.hs b/Bio/Phylogeny/PhyBin/RFDistance.hs --- a/Bio/Phylogeny/PhyBin/RFDistance.hs +++ b/Bio/Phylogeny/PhyBin/RFDistance.hs @@ -14,7 +14,7 @@ denseUnions, denseDiff, invertDense, markLabel, -- * Methods for computing distance matrices - distanceMatrix, hashRF, + naiveDistMatrix, hashRF, -- * Output printDistMat) @@ -40,6 +40,7 @@ -- import Data.LVar.NatArray as NA import Bio.Phylogeny.PhyBin.CoreTypes +import Bio.Phylogeny.PhyBin.PreProcessor (pruneTreeLeaves) -- import Data.BitList import qualified Data.Set as S import qualified Data.List as L @@ -155,17 +156,51 @@ -- | Returns a triangular distance matrix encoded as a vector. -- Also return the set-of-BIPs representation for each tree. -distanceMatrix :: [NewickTree a] -> (DistanceMatrix, V.Vector (S.Set DenseLabelSet)) -distanceMatrix lst = +-- +-- This uses a naive method, directly computing the pairwise +-- distance between each pair of trees. +-- +-- This method is TOLERANT of differences in the laba/taxa sets between two trees. +-- It simply prunes to the intersection before doing the distance comparison. +-- Other scoring methods may be added in the future. (For example, penalizing for +-- missing taxa.) +naiveDistMatrix :: [NewickTree DefDecor] -> (DistanceMatrix, V.Vector (S.Set DenseLabelSet)) +naiveDistMatrix lst = let sz = P.length lst - eachbips = V.fromList $ map allBips lst + treeVect = V.fromList lst + labelSets = V.map treeLabels treeVect + eachbips = V.map allBips treeVect mat = V.generate sz $ \ i -> U.generate i $ \ j -> - let diff1 = S.size (S.difference (eachbips V.! i) (eachbips V.! j)) - diff2 = S.size (S.difference (eachbips V.! j) (eachbips V.! i)) - in diff1 + diff2 + let + inI = (labelSets V.! i) + inJ = (labelSets V.! j) + inBoth = S.intersection inI inJ + + -- Match will always succeed due to size==0 test below: + Just prI = pruneTreeLeaves inBoth (treeVect V.! i) + Just prJ = pruneTreeLeaves inBoth (treeVect V.! j) + + -- Memoization: If we are using it at its full size we can use the cached one: + bipsI = if S.size inBoth == S.size inI + then (eachbips V.! i) + else allBips prI + bipsJ = if S.size inBoth == S.size inJ + then (eachbips V.! j) + else allBips prJ + + diff1 = S.size (S.difference bipsI bipsJ) + diff2 = S.size (S.difference bipsJ bipsI) -- symettric difference + in if S.size inBoth == 0 + then 0 -- This is weird, but what other answer could we give? + else diff1 + diff2 in (mat, eachbips) + where + treeLabels :: NewickTree a -> S.Set Label + treeLabels (NTLeaf _ lab) = S.singleton lab + treeLabels (NTInterior _ ls) = S.unions (map treeLabels ls) + -- | The number of bipartitions implied by a tree is one per EDGE in the tree. Thus -- each interior node carries a list of BiPs the same length as its list of children. labelBips :: NewickTree a -> NewickTree (a, [DenseLabelSet]) @@ -305,37 +340,6 @@ T.traverse (U.unsafeFreeze) v1 -#if 0 --- | Returns a (square) distance matrix encoded as a vector. -distanceMatrix :: [AnnotatedTree] -> IO (U.Vector Word) -distanceMatrix lst = do --- IM.IMapSnap (table :: M.Map DenseLabelSet (S.Set TreeID)) <- runParThenFreeze par --- IM.IMapSnap (table :: M.Map DenseLabelSet (Snapshot IS.ISet TreeID)) <- runParThenFreeze par - IM.IMapSnap table <- runParThenFreeze par - let sz = P.length lst - v <- MU.replicate (sz*sz) (0::Word) - let fn set () = - - F.foldrM - undefined - - -- runParThenFreeze -- get bip table - -- followed by ... fill matrix from bip table - where - par = do - table <- IM.newEmptyMap - forM_ lst (insertBips table) - return table - -insertBips :: BipTable s -> AnnotatedTree -> Par d s () -insertBips table tree = do - let bips = allBips tree - fn bip () = do - IM.modify table bip (IS.putInSet tree) - return () - F.foldrM fn () bips -#endif - -------------------------------------------------------------------------------- -- Miscellaneous Helpers -------------------------------------------------------------------------------- @@ -372,13 +376,18 @@ [ tr | tr <- trees , cbips `S.isSubsetOf` allBips tr ] --- | Is a tree compatible with a consensus? +-- | `compatibleWith consensus tree` -- Is a tree compatible with a consensus? -- This is more efficient if partially applied then used repeatedly. +-- +-- Note, tree compatibility is not the same as an exact match. It's +-- like (<=) rather than (==). The "star topology" is consistent with the +-- all trees, because it induces the empty set of bipartitions. compatibleWith :: NewickTree a -> NewickTree b -> Bool -compatibleWith consensus newTr = - S.isSubsetOf (allBips consensus) (allBips newTr) +compatibleWith consensus = + let consBips = allBips consensus in + \ newTr -> S.isSubsetOf consBips (allBips newTr) --- Consensus between two trees, which may even have different label maps. +-- | Consensus between two trees, which may even have different label maps. consensusTreeFull (FullTree n1 l1 t1) (FullTree n2 l2 t2) = error "FINISHME - consensusTreeFull" diff --git a/Bio/Phylogeny/PhyBin/Util.hs b/Bio/Phylogeny/PhyBin/Util.hs --- a/Bio/Phylogeny/PhyBin/Util.hs +++ b/Bio/Phylogeny/PhyBin/Util.hs @@ -9,7 +9,8 @@ module Bio.Phylogeny.PhyBin.Util ( - is_regular_file, acquireTreeFiles + is_regular_file, acquireTreeFiles, + safePrintDendro, sanityCheck ) where @@ -17,9 +18,9 @@ import Data.Function (on) import Data.List (delete, minimumBy, sortBy, insertBy, intersperse, sort) import Data.Maybe (fromMaybe, catMaybes) -import qualified Data.ByteString.Lazy.Char8 as B import qualified Data.Map as M import qualified Data.Set as S +import qualified Data.Text.Lazy as T import Control.Monad (forM, forM_, filterM, when, unless) import Control.Exception (evaluate) import Control.Applicative ((<$>),(<*>)) @@ -30,16 +31,20 @@ import System.IO (openFile, hClose, IOMode(ReadMode), stderr, hPutStr, hPutStrLn) import System.Exit (ExitCode(..)) +import System.Timeout (timeout) import Test.HUnit ((~:),(~=?),Test,test) -- For vizualization: import Text.PrettyPrint.HughesPJClass hiding (char, Style) import Bio.Phylogeny.PhyBin.CoreTypes -import Bio.Phylogeny.PhyBin.Parser (parseNewick) -import Bio.Phylogeny.PhyBin.PreProcessor (collapseBranches) -import Bio.Phylogeny.PhyBin.Visualize (dotToPDF, dotNewickTree, viewNewickTree) -import Bio.Phylogeny.PhyBin.RFDistance +-- import Bio.Phylogeny.PhyBin.Parser (parseNewick) +-- import Bio.Phylogeny.PhyBin.Visualize (dotToPDF, dotNewickTree, viewNewickTree) +import qualified Data.Clustering.Hierarchical as C +import qualified Data.Graph.Inductive as G +import qualified Data.GraphViz as Gv +import Data.GraphViz.Printing (renderDot) +import Data.GraphViz.Types.Canonical (nodeStmts, graphStatements) ---------------------------------------------------------------------------------------------------- -- OS specific bits: @@ -77,7 +82,6 @@ d <- doesDirectoryExist path return (f || d) - -------------------------------------------------------------------------------- -- | Expand out directories to find all the tree files. @@ -109,3 +113,27 @@ return (concat all) +-------------------------------------------------------------------------------- + +-- | Step carefully in case of cycles (argh). +safePrintDendro :: Gv.DotGraph G.Node -> IO (Maybe T.Text) +safePrintDendro dotg= do +-- putStrLn$ "Dendrogram graph size: "++ show (F.foldl' (\a _ -> a+1) 0 dotg) + mx <- timeout (2 * 1000 * 1000) $ do +-- putStrLn$ "Dendrogram graph, is directed?: "++ show (Gv.directedGraph dotg) + putStrLn$ "Dendrogram graph size: "++ show (length $ nodeStmts $ graphStatements dotg) + let str = renderDot $ Gv.toDot dotg + evaluate (T.length str) + return str + case mx of + Nothing -> do putStrLn "WARNING: DotGraph appears to be a cyclic structure. This is probably a bug." + return Nothing + _ -> return mx + +sanityCheck :: C.Dendrogram (FullTree DefDecor) -> IO () +sanityCheck dendro = do + let fn seen elm | S.member (treename elm) seen = + error$"Dendrogram failed sanity check! Tree name occurs multiple times: "++(treename elm) + | otherwise = S.insert (treename elm) seen + sz = S.size $ F.foldl' fn S.empty dendro + putStrLn$ "Sanity checked dendrogram of size: "++show sz diff --git a/Bio/Phylogeny/PhyBin/Visualize.hs b/Bio/Phylogeny/PhyBin/Visualize.hs --- a/Bio/Phylogeny/PhyBin/Visualize.hs +++ b/Bio/Phylogeny/PhyBin/Visualize.hs @@ -23,13 +23,18 @@ import qualified Data.GraphViz.Attributes.Complete as GA import qualified Data.GraphViz.Attributes.Colors as GC import Data.GraphViz.Attributes.Colors (Color(RGB)) + -- import Test.HUnit ((~:),(~=?),Test,test) +import System.Timeout (timeout) + import qualified Data.Clustering.Hierarchical as C import Bio.Phylogeny.PhyBin.CoreTypes import Bio.Phylogeny.PhyBin.RFDistance (filterCompatible, compatibleWith, consensusTree) +import Debug.Trace + ---------------------------------------------------------------------------------------------------- -- Visualization with GraphViz and FGL: ---------------------------------------------------------------------------------------------------- @@ -68,7 +73,7 @@ -- Dendrograms ---------------------------------------------------------------------------------------------------- --- | Some duplicated code with dotNewickTree. +-- | Convert to a dotGraph. Some duplicated code with dotNewickTree. dotDendrogram :: PhyBinConfig -> String -> Double -> C.Dendrogram (FullTree a) -> Maybe (M.Map TreeName Int) -> [[NewickTree ()]] -> Gv.DotGraph G.Node dotDendrogram PBC{show_trees_in_dendro, show_interior_consensus} @@ -235,11 +240,17 @@ -- The channel retuned will carry a single message to signal -- completion of the subprocess. viewNewickTree :: String -> FullTree StandardDecor -> IO (Chan (), FullTree StandardDecor) +-- TODO: UPDATE THIS TO RETURN AN ASYNC!! viewNewickTree title tree@(FullTree{nwtree}) = do chan <- newChan let dot = dotNewickTree title (1.0 / avg_branchlen [nwtree]) tree - runit = do Gv.runGraphvizCanvas default_cmd dot Gv.Xlib + runit = do mx <- -- timeout defaultTimeout $ + -- This one is interactive, so we don't need a timeout. + Gv.runGraphvizCanvas default_cmd dot Gv.Xlib + -- case mx of + -- Nothing -> putStrLn$ "WARNING: call to graphviz TIMED OUT."++file + -- _ -> return () writeChan chan () --str <- prettyPrint d --putStrLn$ "Generating the following graphviz tree:\n " ++ str @@ -262,9 +273,21 @@ then show rnd else printf "%.4f" fl -dotToPDF :: Gv.DotGraph G.Node -> FilePath -> IO FilePath -dotToPDF dot file = - Gv.runGraphvizCommand default_cmd dot Gv.Pdf file +-- | Convert a .dot file to .pdf. +dotToPDF :: Gv.DotGraph G.Node -> FilePath -> IO (Maybe FilePath) +dotToPDF dot file = do + -- If we have any problem with graphviz we want to time that out rather than let + -- our whole run hang: + x <- timeout defaultTimeout $ + Gv.runGraphvizCommand default_cmd dot Gv.Pdf file + case x of + Nothing -> do putStrLn$ "WARNING: call to graphviz TIMED OUT. File not plotted: "++file + return Nothing + _ -> return x + +-- Arbitrary: 15 second timeout. +defaultTimeout :: Int +defaultTimeout = (15 * 1000 * 1000) -- | Convert a NewickTree to a graphviz Dot graph representation. dotNewickTree :: String -> Double -> FullTree StandardDecor -> Gv.DotGraph G.Node diff --git a/Main.hs b/Main.hs --- a/Main.hs +++ b/Main.hs @@ -1,40 +1,39 @@ {-# LANGUAGE RecordWildCards, TupleSections, NamedFieldPuns #-} {-# OPTIONS_GHC -fwarn-unused-imports -fwarn-incomplete-patterns #-} +-- | The MAIN module, of course. This is the script that deals with +-- command line options and calling into the heart of the beast. + module Main where import Data.List (sort, intersperse, foldl') -import qualified Data.ByteString.Lazy.Char8 as B +import Data.List.Split (splitOneOf) +import qualified Data.ByteString.Char8 as B import qualified Data.Map as M import qualified Data.Set as S -import qualified Data.IntSet as IS import Control.Monad -import Control.Concurrent (Chan, readChan, ThreadId, forkIO) +import Control.Concurrent (readChan) import System.Environment (getArgs, withArgs) import System.Console.GetOpt (OptDescr(Option), ArgDescr(..), ArgOrder(..), usageInfo, getOpt) import System.Exit (exitSuccess) import System.IO (stdout) -import Test.HUnit (runTestTT, Test, test, (~:)) -import Control.Applicative ((<$>)) +import Control.Exception (catch, SomeException) import qualified Data.Vector as V import Test.HUnit as HU -import Data.GraphViz (runGraphvizCanvas,GraphvizCommand(Dot),GraphvizCanvas(Xlib)) -import Text.PrettyPrint.HughesPJClass hiding (char, Style) - -import Bio.Phylogeny.PhyBin.CoreTypes -import Bio.Phylogeny.PhyBin (driver, binthem, normalize, annotateWLabLists, +import Bio.Phylogeny.PhyBin.CoreTypes +import Bio.Phylogeny.PhyBin (driver, normalize, annotateWLabLists, unitTests, acquireTreeFiles, deAnnotate, retrieveHighlights, matchAnyHighlight) -import Bio.Phylogeny.PhyBin.Parser (parseNewick, parseNewicks, parseNewickFiles, unitTests) -import Bio.Phylogeny.PhyBin.Visualize (viewNewickTree, dotNewickTree_debug) +import Bio.Phylogeny.PhyBin.Parser (parseNewick, parseNewickFiles, unitTests) +import Bio.Phylogeny.PhyBin.Visualize (viewNewickTree) -- dotNewickTree_debug import Bio.Phylogeny.PhyBin.RFDistance import Bio.Phylogeny.PhyBin.PreProcessor import qualified Data.Clustering.Hierarchical as C - -import Version +import Data.Version (showVersion) +import Paths_phybin (version) ---------------------------------------------------------------------------------------------------- -- MAIN script: Read command line options and call the program. @@ -45,8 +44,10 @@ = Verbose | Version | Output String - | NumTaxa Int - | BranchThresh Double + | SetNumTaxa Int + | BranchThresh Double + | BootStrapThresh Int + | PruneTaxa [String] | NullOpt | Graph | Draw | Force @@ -56,7 +57,7 @@ | Highlight FilePath | ShowTreesInDendro | ShowInterior - | HashRF Bool + | RFMode WhichRFMode | SelfTest | RFMatrix | LineSetDiffMode | PrintNorms | PrintReg | PrintConsensus | PrintMatching @@ -70,8 +71,7 @@ | NameTable String -- Must come after Cutoff/Prefix deriving (Show, Eq, Ord) -- ORD is really used. - - + parseTabDelim :: String -> Flag parseTabDelim _str = TabDelimited 8 9 @@ -102,11 +102,12 @@ "Irrespective of whether this is activated, a hierarchical clustering (dendogram.pdf) is produced." -- , Option [] ["dendogram"] (NoArg DendogramOnly)$ "Report a hierarchical clustering (default)" - , Option [] [] (NoArg$ error "internal problem") " Select Robinson-Foulds (symmetric difference) distance algorithm:" - , Option [] ["simple"] (NoArg$ HashRF False) - ((if hashRF then "" else "(default) ")++ "use direct all-to-all comparison") - , Option [] ["hashrf"] (NoArg$ HashRF True) - ((if hashRF then "(default) " else "")++"use a variant of the HashRF algorithm for the distance matrix") + , Option [] [] (NoArg NullOpt) "" + , Option [] [] (NoArg$ error "internal problem") "--- Select Robinson-Foulds (symmetric difference) distance algorithm: ---" + , Option [] ["hashrf"] (NoArg$ RFMode HashRF) + "(default) use a variant of the HashRF algorithm for the distance matrix" + , Option [] ["tolerant"] (NoArg$ RFMode TolerantNaive) + "use a slower, modified RF metric that tolerates missing taxa" , Option [] [] (NoArg NullOpt) "" , Option [] [] (NoArg$ error "internal problem") "----------------------------- Visualization --------------------------------" @@ -120,19 +121,25 @@ , Option [] ["highlight"] (ReqArg Highlight "FILE") $ "Highlight nodes in the tree-of-trees (dendrogram) consistent with the.\n"++ "given tree file. Multiple highlights are permitted and use different colors." +-- "Highlights may be SUBTREES that use any subset of the taxa." , Option [] ["interior"] (NoArg ShowInterior) "Show the consensus trees for interior nodes in the dendogram, rather than just points." , Option [] [] (NoArg NullOpt) "" , Option [] [] (NoArg$ error "internal problem") "---------------------------- Tree pre-processing -----------------------------" - - , Option ['n'] ["numtaxa"] (ReqArg (NumTaxa . read) "NUM") "expect NUM taxa for this dataset" + , Option [] ["prune"] (ReqArg (PruneTaxa . parseTaxa) "TAXA") + ("Prune trees to only TAXA before doing anything else.\n"++ + "Space and comma separated lists of taxa are allowed. Use quotes.") + , Option ['b'] ["minbranchlen"] (ReqArg (BranchThresh . read) "LEN") "collapse branches less than LEN" - , Option ['b'] ["branchcut"] (ReqArg (BranchThresh . read) "LEN") "collapse branches less than LEN" - + , Option [] ["minbootstrap"] (ReqArg (BootStrapThresh . read) "INT") + "collapse branches with bootstrap values less than INT" + -- , Option ['n'] ["numtaxa"] (ReqArg (SetNumTaxa . read) "NUM") + -- "expect NUM taxa for this dataset, demand trees have all of them" + , Option [] [] (NoArg NullOpt) "" , Option [] [] (NoArg$ error "internal problem") "--------------------------- Extracting taxa names ----------------------------" --- , Option ['n'] ["numtaxa"] (ReqArg (NumTaxa . read) "NUM") "expect NUM taxa for this dataset (otherwise it will guess)" +-- , Option ['n'] ["numtaxa"] (ReqArg (SetNumTaxa . read) "NUM") "expect NUM taxa for this dataset (otherwise it will guess)" -- TODO, FIXME: The "guessing" part doesn't actually work yet -- implement it!! -- What's a good algorithm? Insist they all have the same number? Take the mode? @@ -167,9 +174,13 @@ , Option [] ["consensus"] (NoArg PrintConsensus) "print a strict consensus tree for the inputs, then exit" , Option [] ["matching"] (NoArg PrintMatching) "print a list of tree names that match any --highlight argument" ] - where - hashRF = use_hashrf default_phybin_config + +-- | Parse the list of taxa provided on the command line. +parseTaxa :: String -> [String] +parseTaxa str = filter (not . null) $ + splitOneOf ",; " str + usage :: String usage = "\nUsage: phybin [OPTION...] files or directories...\n\n"++ @@ -190,7 +201,7 @@ "When clustering trees, Phybin computes a complete all-to-all Robinson-Foulds distance matrix.\n"++ "If a threshold distance (tree edit distance) is given, then a flat set of clusters\n"++ - "will be produced in files clusterXX_YY.tr. Otherwise it produces a full dendogram (UNFINISHED).\n"++ + "will be produced in files clusterXX_YY.tr. Otherwise it produces a full dendogram.\n"++ "\n"++ "Binning mode provides an especially quick-and-dirty form of clustering.\n"++ @@ -222,7 +233,7 @@ let process_opt cfg opt = case opt of NullOpt -> return cfg Verbose -> return cfg { verbose= True } - Version -> do putStrLn$ "phybin version "++phybin_version; exitSuccess + Version -> do putStrLn$ "phybin version "++ showVersion version; exitSuccess SelfTest -> do _ <- runTestTT allUnitTests; exitSuccess @@ -248,15 +259,24 @@ PrintMatching -> return cfg Cluster lnk -> return cfg { clust_mode = ClusterThem lnk } - HashRF bl -> return cfg { use_hashrf = bl } + RFMode TolerantNaive + | Expected _ <- num_taxa cfg -> error "should not set --numtaxa AND --tolerant" + | otherwise -> return cfg { rfmode = TolerantNaive } + RFMode HashRF -> return cfg { rfmode = HashRF } BinningMode -> return cfg { clust_mode = BinThem } EditDistThresh n -> return cfg { dist_thresh = Just n } DendogramOnly -> return cfg { dist_thresh = Nothing } Output s -> return cfg { output_dir= s } - NumTaxa n -> return cfg { num_taxa= n } - BranchThresh n -> return cfg { branch_collapse_thresh= Just n } + SetNumTaxa n + | rfmode cfg == TolerantNaive -> error "should not set --tolerant AND --numtaxa" + | otherwise -> return cfg { num_taxa= Expected n } + BranchThresh n -> return cfg { branch_collapse_thresh = Just n } + BootStrapThresh n -> return cfg { bootstrap_collapse_thresh = Just n } + + PruneTaxa ls -> return cfg { preprune_labels = Just ls } + Graph -> return cfg { do_graph= True } Draw -> return cfg { do_draw = True } View -> return cfg -- Handled below @@ -279,7 +299,9 @@ when (null files) $ do defaultErr ["No file arguments!\n"] - ------------------------------------------------------------ + -------------------------------------------------------------------- + ----- Handle SPECIAL modes before going into normal processing ----- + -------------------------------------------------------------------- -- This mode kicks in AFTER config options are processed. when (elem PrintReg opts) $ do (_,fts) <- parseNewickFiles (name_hack config) all_inputs @@ -296,14 +318,17 @@ liftFT (normalize . annotateWLabLists) ft exitSuccess ------------------------------------------------------------ - when (elem PrintConsensus opts) $ do - (_,fts) <- parseNewickFiles (name_hack config) all_inputs - putStrLn $ "Strict Consensus Tree of "++show (length fts)++" trees:" - when (null fts) $ error "No trees provided!" - let ctree = consensusTree (num_taxa config) (map nwtree fts) - FullTree{labelTable} = head fts - print$ displayStrippedTree$ FullTree "" labelTable ctree - exitSuccess + when (elem PrintConsensus opts) $ + case (num_taxa config) of + Expected numtax -> do + (_,fts) <- parseNewickFiles (name_hack config) all_inputs + putStrLn $ "Strict Consensus Tree of "++show (length fts)++" trees:" + when (null fts) $ error "No trees provided!" + let ctree = consensusTree numtax (map nwtree fts) + FullTree{labelTable} = head fts + print$ displayStrippedTree$ FullTree "" labelTable ctree + exitSuccess + _ -> error "--consensus: cannot compute consensus when number of taxa is unknown or variable" ------------------------------------------------------------ when (elem PrintMatching opts) $ do (labs,fts) <- parseNewickFiles (name_hack config) all_inputs @@ -322,9 +347,14 @@ -- Otherwise do the main, normal thing: driver config +-- | Completely optional visualization step. If this fails, we don't +-- sweat it. view_graphs :: PhyBinConfig -> IO () -view_graphs PBC{..} = - do chans <- forM inputs $ \ file -> do +view_graphs PBC{..} = catch act hndl + where + hndl :: SomeException -> IO () + hndl exn = putStrLn$ "Warning: Viewing graph failed, ignoring. Error was: "++show exn + act = do chans <- forM inputs $ \ file -> do putStrLn$ "Drawing "++ file ++"...\n" str <- B.readFile file putStrLn$ "Parsed: " ++ (B.unpack str) @@ -364,7 +394,7 @@ (one:rest) -> [one, unwords rest] temp :: IO () -temp = driver default_phybin_config{ num_taxa=7, inputs=["../datasets/test.tr"] } +temp = driver default_phybin_config{ num_taxa= Expected 7, inputs=["../datasets/test.tr"] } -------------------------------------------------------------------------------- @@ -378,7 +408,9 @@ , Bio.Phylogeny.PhyBin.Parser.unitTests , "norm/Bip1" ~: (testNorm prob1) , "bipTreeConversion" ~: testBipConversion - , "t3" ~: t3_consensusTest, "t4" ~: t4_consensusTest, "t5" ~: t5_consensusTest + , "t3" ~: t3_consensusTest, + "t4" ~: t4_consensusTest, + "t5" ~: t5_consensusTest ] -- [2013.07.23] @@ -412,7 +444,7 @@ -- 112 and 13 rftest :: IO () rftest = do - (mp,[t1,t2]) <- parseNewickFiles (take 2) ["tests/13.tr", "tests/112.tr"] + (_mp,[t1,t2]) <- parseNewickFiles (take 2) ["tests/13.tr", "tests/112.tr"] putStrLn$ "Tree 13 : " ++ show (displayDefaultTree t1) putStrLn$ "Tree 112 : "++ show (displayDefaultTree t2) @@ -425,7 +457,7 @@ putStrLn$ "13 collapsed 0.03: " ++show (disp$ liftFT (collapseBranchLenThresh 0.03) t1) putStrLn$ "112 collapsed 0.03: " ++show (disp$ liftFT (collapseBranchLenThresh 0.03) t2) - let (mat,_) = distanceMatrix [nwtree t1, nwtree t2] + let (mat,_) = naiveDistMatrix [nwtree t1, nwtree t2] printDistMat stdout mat return () where @@ -478,7 +510,7 @@ cbips cbips2 -- (allBips$ fmap (const ()) $ nwtree ctree) putStrLn " Partial distance matrix WITHIN this cluster:" - let (mat,_) = distanceMatrix (map nwtree ftrees) + let (mat,_) = naiveDistMatrix (map nwtree ftrees) printDistMat stdout (V.take 30 mat) HU.assertBool "Consensus should only include bicbips2ps in the members" (S.isSubsetOf cbips totalBips) HU.assertEqual "Consensus tree matches intersected bips" cbips intersectBips diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -1,7 +1,7 @@ -% PhyBin: Binning Trees by Topology +% PhyBin (0.3): Binning/Clustering Newick Trees by Topology -PhyBin is a simple command line tool that classifies (bins) a set of +PhyBin is a simple command line tool that classifies a set of [Newick tree files](http://en.wikipedia.org/wiki/Newick_format) by their topology. The purpose of it is to take a large set of tree files and browse through the most common tree topologies. @@ -10,6 +10,22 @@ phybin run. The file `binXX_YYY`, where `XX` is the rank of the bin and `YYY` is the number of trees having that topology.](trees.jpg) +Change Log +========== + +In version 0.2, PhyBin was extended to do clustering as well as binning: + + * Computee full all-to-all Robinson-Foulds distance matrices (quickly) + * Hierarchical clustering of all trees into a tree-of-trees dendrogram based on + Robinson Foulds symmetric (tree edit) distance. + +In version 0.3, PhyBin gained a number of new features + + * A `--tolerant` mode for computing RF distance matrices even for trees missing taxa. + * A `--prune` option for "zooming in" on a specific set of taxa. + * The `--minboostrap` option was added. + + Invoking PhyBin =============== @@ -18,25 +34,25 @@ [GraphViz program](http://www.graphviz.org/), including the `dot` command, must be installed on the machine. -The following is a typical invocation of PhyBin: +The following is a simple invocation of PhyBin: - phybin *.tree -o output_dir/ + phybin --bin *.tree -o output_dir/ The input trees can be specified directly on the command-line, or, if the name of a directory is provided instead, all contained files are assumed to be trees in Newick format. PhyBin, at minimum, produces files of the form -`output_dir/binXX_YY.tr`, one for each bin. If +`output_dir/clusterXX_YY.tr`, one for each bin. If requested, it will also produce visual representations of each bin in -the form `output_dir/binXX_YY.pdf`. +the form `output_dir/clusterXX_YY.pdf`. Downloading and Installing PhyBin ================================= The source code to PhyBin can be downloaded here: - * [Download Source Tarball](phybin-0.1.2.tar.gz) + * [Download Source from Hackage](http://hackage.haskell.org/package/phybin) PhyBin is written in Haskell and if you have [Haskell Platform](http://hackage.haskell.org/platform/). @@ -44,27 +60,106 @@ cabal install phybin -PhyBin is also available for download as a statically-linked -executable for Mac-OS and Linux: +Otherwise PhyBin is also available for download as a statically-linked +executable for Mac-OS, Linux, and Windows: - * [Download Mac-OS Binary](phybin-0.1.2.mac) - * [Download Linux Binary](phybin-0.1.2.x86_64_linux) + * [Download Mac-OS Binary](phybin-0.3.mac) + * [Download Linux Binary](phybin-0.3.x86_64_linux) + * [Download Windows Binary](phybin-0.3_windows.exe) -It should be possible to build it for Windows, but I haven't done so -yet. - Command-line Options ==================== In addition to input files and directories, `phybin` supports a number of command-line options. Run "phybin --help" to see these options. +Here is a snapshot of the current help output (version 0.2.11): + Usage: phybin [OPTION...] files or directories... + + PhyBin takes Newick tree files as input. Paths of Newick files can + be passed directly on the command line. Or, if directories are provided, + all files in those directories will be read. Taxa are named based on the files + containing them. If a file contains multiple trees, all are read by phybin, and + the taxa name then includes a suffix indicating the position in the file: + e.g. FILENAME_0, FILENAME_1, etc. + + When clustering trees, Phybin computes a complete all-to-all Robinson-Foulds distance matrix. + If a threshold distance (tree edit distance) is given, then a flat set of clusters + will be produced in files clusterXX_YY.tr. Otherwise it produces a full dendogram (UNFINISHED). + + Binning mode provides an especially quick-and-dirty form of clustering. + When running with the --bin option, only exactly equal trees are put in the same cluster. + Tree pre-processing still applies, however: for example collapsing short branches. + + USAGE NOTES: + * Currently phybin ignores input trees with the wrong number of taxa. + * If given a directory as input phybin will assume all contained files are Newick trees. + + + Options include: + + -v --verbose print WARNINGS and other information (recommended at first) + -V --version show version number + -o DIR --output=DIR set directory to contain all output files (default "./phybin_out/") + --selftest run internal unit tests + + ----------------------------- Clustering Options ------------------------------ + --bin Use simple binning, the cheapest form of 'clustering' + --single Use single-linkage clustering (nearest neighbor) + --complete Use complete-linkage clustering (furthest neighbor) + --UPGMA Use Unweighted Pair Group Method (average linkage) - DEFAULT mode + --editdist=DIST Combine all clusters separated by DIST or less. Report a flat list of clusters. + Irrespective of whether this is activated, a hierarchical clustering (dendogram.pdf) is produced. + Select Robinson-Foulds (symmetric difference) distance algorithm: + --simple use direct all-to-all comparison + --hashrf (default) use a variant of the HashRF algorithm for the distance matrix + + ----------------------------- Visualization -------------------------------- + -g --graphbins use graphviz to produce .dot and .pdf output files + -d --drawbins like -g, but open GUI windows to show each bin's tree + -w --view for convenience, "view mode" simply displays input Newick files without binning + --showtrees Print (textual) tree topology inside the nodes of the dendrogram + --highlight=FILE Highlight nodes in the tree-of-trees (dendrogram) consistent with the. + given tree file. Multiple highlights are permitted and use different colors. + --interior Show the consensus trees for interior nodes in the dendogram, rather than just points. + + ---------------------------- Tree pre-processing ----------------------------- + -n NUM --numtaxa=NUM expect NUM taxa for this dataset + -b LEN --branchcut=LEN collapse branches less than LEN + + --------------------------- Extracting taxa names ---------------------------- + + -p NUM --nameprefix=NUM Leaf names in the input Newick trees can be gene names, not taxa. + Then it is typical to extract taxa names from genes. This option extracts + a prefix of NUM characters to serve as the taxa name. + + -s STR --namesep=STR An alternative to --nameprefix, STR provides a set of delimeter characters, + for example '-' or '0123456789'. The taxa name is then a variable-length + prefix of each gene name up to but not including any character in STR. + + -m FILE --namemap=FILE Even once prefixes are extracted it may be necessary to use a lookup table + to compute taxa names, e.g. if multiple genes/plasmids map onto one taxa. + This option specifies a text file with find/replace entries of the form + " ", which are applied AFTER -s and -p. + + --------------------------- Utility Modes ---------------------------- + --rfdist print a Robinson Foulds distance matrix for the input trees + --setdiff for convenience, print the set difference between cluster*.txt files + --print simply print out a concise form of each input tree + --printnorms simply print out a concise and NORMALIZED form of each input tree + --consensus print a strict consensus tree for the inputs, then exit + --matching print a list of tree names that match any --highlight argument + + - - - - - - - - - - - - - - - Authors: Irene and Ryan Newton -Contact email: `irnewton` `indiana` `edu` (with "at" and "dot" inserted). +Contact email: `irnewton` and `rrnewton` at `indiana` `edu` (with "at" and "dot" inserted). + +[Irene's](http://www.bio.indiana.edu/faculty/directory/profile.php?person=irnewton) and +[Ryan](http://www.cs.indiana.edu/~rrnewton/homepage.html) homepages. . diff --git a/TestAll.hs b/TestAll.hs new file mode 100644 --- /dev/null +++ b/TestAll.hs @@ -0,0 +1,41 @@ +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE ScopedTypeVariables, BangPatterns, ParallelListComp, OverloadedStrings #-} +module Main where + +import Data.Vector as V +import Data.Vector.Unboxed as U +import Data.ByteString.Char8 as B +import Bio.Phylogeny.PhyBin +import Bio.Phylogeny.PhyBin.CoreTypes +import Bio.Phylogeny.PhyBin.RFDistance +import Bio.Phylogeny.PhyBin.Parser +import System.IO +import Prelude as P + +import qualified Test.HUnit as HU +-- import Test.Framework +import Test.Framework.Providers.HUnit +import Test.Framework.TH (defaultMainGenerator) + +-- Here we compare a tree against the same tree with 18 removed: +-- Or the same tree with one intermediate node removed. +(_,t30) = parseNewicks id $ + [ ("t30_0", "(2_, (((14, 3_), (19, (5_, 13))), (18, (6_, 7_))), 1_);") + , ("t30_1", "(2_, (((14, 3_), (19, (5_, 13))), (6_, 7_)), 1_);") + , ("t30_2", "(2_, (((14, 3_), (19, (5_, 13))), (18, 6_, 7_)), 1_);") + ] + +t30' = naiveDistMatrix$ P.map nwtree t30 + -- "[[],[1]]" + +case_t30 = HU.assertEqual "3-tree distance matrix" + "[[],[0],[1,0]]" (showMat t30') + +-- m = printDistMat stdout (fst d) + + +-- Simple show for a distance matrix: +showMat (m,_) = show$ V.toList$ V.map U.toList m + +main = $(defaultMainGenerator) + diff --git a/phybin.cabal b/phybin.cabal --- a/phybin.cabal +++ b/phybin.cabal @@ -1,5 +1,5 @@ Name: phybin -Version: 0.2.11 +Version: 0.3 License: BSD3 License-file: LICENSE Stability: Beta @@ -21,6 +21,7 @@ -- 0.2.9 -- Add command line opt --interior -- 0.2.10 -- Add command line opt --matching -- 0.2.11 -- Cleanup and windows compatibility. +-- 0.3 -- significant improvements and new functionality. -- homepage: http://code.haskell.org/phybin -- homepage: http://people.csail.mit.edu/newton/phybin/ @@ -38,7 +39,7 @@ generating GraphViz-based visual representations of the tree topologies). Category: Bioinformatics -Cabal-Version: >=1.8 +Cabal-Version: >= 1.8 Extra-source-files: README.md @@ -62,13 +63,14 @@ Library Exposed-modules: Bio.Phylogeny.PhyBin - Bio.Phylogeny.PhyBin.RFDistance - Bio.Phylogeny.PhyBin.Parser - Bio.Phylogeny.PhyBin.CoreTypes Bio.Phylogeny.PhyBin.Binning + Bio.Phylogeny.PhyBin.CoreTypes + Bio.Phylogeny.PhyBin.Parser + Bio.Phylogeny.PhyBin.PreProcessor + Bio.Phylogeny.PhyBin.RFDistance Bio.Phylogeny.PhyBin.Util Bio.Phylogeny.PhyBin.Visualize --- Other-modules: Data.BitList + -- Data.BitList Build-Depends: base >= 3 && < 5, directory, process, containers, async, time, filepath, @@ -107,7 +109,7 @@ vector >= 0.10, hierarchical-clustering >= 0.4, split >= 0.2 -- lattice-par, - GHC-Options: -O2 -funbox-strict-fields -rtsopts + GHC-Options: -O2 -funbox-strict-fields -rtsopts -threaded if flag(hashrf) CPP-Options: -DUSE_HASHRF @@ -116,3 +118,32 @@ Build-Depends: bitvec >= 0.1 if flag(sequential) CPP-Options: -DSEQUENTIALIZE + + +-- DUPLICATED from executable: +Test-Suite test-phybin + Main-is: TestAll.hs + Type: exitcode-stdio-1.0 +-- Build-Depends: phybin + Build-Depends: base >= 3 && < 5, directory, process, containers, + async, time, + filepath, + graphviz >= 2999.16, + text >= 0.11 && < 0.12, + prettyclass, fgl, + HUnit, bytestring, + parsec >= 3.1.0, + vector >= 0.10, + hierarchical-clustering >= 0.4, split >= 0.2 + -- Additional depends for test: + Build-Depends: HUnit, test-framework, test-framework-hunit, test-framework-th + GHC-Options: -O2 -funbox-strict-fields -rtsopts -threaded + + if flag(hashrf) + CPP-Options: -DUSE_HASHRF + if flag(bitvec) + CPP-Options: -DBITVEC_BIPS + Build-Depends: bitvec >= 0.1 + if flag(sequential) + CPP-Options: -DSEQUENTIALIZE +