phybin 0.2.11 → 0.3
raw patch · 12 files changed
+649/−217 lines, 12 filesdep +test-frameworkdep +test-framework-hunitdep +test-framework-th
Dependencies added: test-framework, test-framework-hunit, test-framework-th
Files
- Bio/Phylogeny/PhyBin.hs +111/−49
- Bio/Phylogeny/PhyBin/Binning.hs +4/−17
- Bio/Phylogeny/PhyBin/CoreTypes.hs +29/−14
- Bio/Phylogeny/PhyBin/Parser.hs +23/−9
- Bio/Phylogeny/PhyBin/PreProcessor.hs +95/−0
- Bio/Phylogeny/PhyBin/RFDistance.hs +51/−42
- Bio/Phylogeny/PhyBin/Util.hs +35/−7
- Bio/Phylogeny/PhyBin/Visualize.hs +28/−5
- Main.hs +84/−52
- README.md +110/−15
- TestAll.hs +41/−0
- phybin.cabal +38/−7
Bio/Phylogeny/PhyBin.hs view
@@ -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++"<N>_<size>.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.
Bio/Phylogeny/PhyBin/Binning.hs view
@@ -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)))"
Bio/Phylogeny/PhyBin/CoreTypes.hs view
@@ -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 }
Bio/Phylogeny/PhyBin/Parser.hs view
@@ -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
+ Bio/Phylogeny/PhyBin/PreProcessor.hs view
@@ -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)+++
Bio/Phylogeny/PhyBin/RFDistance.hs view
@@ -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"
Bio/Phylogeny/PhyBin/Util.hs view
@@ -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
Bio/Phylogeny/PhyBin/Visualize.hs view
@@ -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
Main.hs view
@@ -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
README.md view
@@ -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+ "<string> <taxaname>", 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. .
+ TestAll.hs view
@@ -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)+
phybin.cabal view
@@ -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+