diff --git a/Bio/Phylogeny/PhyBin.hs b/Bio/Phylogeny/PhyBin.hs new file mode 100644 --- /dev/null +++ b/Bio/Phylogeny/PhyBin.hs @@ -0,0 +1,438 @@ +{-# LANGUAGE ScopedTypeVariables, RecordWildCards, TypeSynonymInstances, CPP #-} +{-# LANGUAGE NamedFieldPuns #-} +{-# LANGUAGE OverloadedStrings #-} +{-# OPTIONS_GHC -fwarn-incomplete-patterns #-} +{-# OPTIONS_GHC -fwarn-unused-imports #-} + +-- | This module contains the code that does the tree normalization and binning. +-- It's the heart of the prgoram. + +module Bio.Phylogeny.PhyBin + ( driver, binthem, normalize, annotateWLabLists, unitTests, acquireTreeFiles, + deAnnotate ) + where + +import qualified Data.Foldable as F +import Data.Function (on) +import Data.List (delete, minimumBy, sortBy, foldl1') +import Data.Maybe (fromMaybe) +import Data.Either (partitionEithers) +import Data.Time.Clock +import qualified Data.ByteString.Lazy.Char8 as B +import qualified Data.Map as M +import qualified Data.Set as S +import qualified Data.Vector as V +import qualified Data.Vector.Unboxed as U +import Control.Monad (forM, forM_, filterM, when, unless) +import Control.Concurrent.Async +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(..), stdout) +import System.Process (system) +import System.Exit (ExitCode(..)) +import Test.HUnit ((~:),(~=?),Test,test) +import qualified Data.Clustering.Hierarchical as C + +-- For vizualization: +import Text.PrettyPrint.HughesPJClass hiding (char, Style) +import Bio.Phylogeny.PhyBin.CoreTypes +import Bio.Phylogeny.PhyBin.Parser (parseNewick, parseNewicks) +import Bio.Phylogeny.PhyBin.PreProcessor (collapseBranchLenThresh) +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 + +-- Turn on for extra invariant checking: +debug :: Bool +debug = True + +---------------------------------------------------------------------------------------------------- + +-- | Driver to put all the pieces together (parse, normalize, bin) +driver :: PhyBinConfig -> IO () +driver PBC{ verbose, num_taxa, name_hack, output_dir, inputs=files, + do_graph, branch_collapse_thresh, + dist_thresh, clust_mode, print_rfmatrix } = + -- Unused: do_draw + do + -------------------------------------------------------------------------------- + -- First, find out where we are and open the files: + -------------------------------------------------------------------------------- + cd <- getCurrentDirectory + --putStrLn$ "PHYBIN RUNNING IN DIRECTORY: "++ cd + + bl <- doesDirectoryExist output_dir + unless bl $ do + c <- system$ "mkdir -p "++output_dir + case c of + ExitSuccess -> return () + ExitFailure cde -> error$"Could not create output directory. 'mkdir' command failed with: "++show cde + + putStrLn$ "Cleaning away previous phybin outputs..." + system$ "rm -f "++output_dir++"/dendrogram.*" + system$ "rm -f "++output_dir++"/cluster*" + system$ "rm -f "++output_dir++"/distance_matrix.txt" + system$ "rm -f "++output_dir++"/WARNINGS.txt" + + putStrLn$ "Parsing "++show (length files)++" Newick tree files." + --putStrLn$ "\nFirst ten \n"++ concat (map (++"\n") $ map show $ take 10 files) + + -------------------------------------------------------------------------------- + -- Next, parse the files and do error checking and annotation. + -------------------------------------------------------------------------------- + + (goodFiles,warnings1) <- fmap partitionEithers $ + forM files $ \ file -> do + reg <- is_regular_file file + if reg + then return$ Left file + else return$ Right file + let num_files = length goodFiles + bstrs <- mapM B.readFile goodFiles + let (labelTab, fulltrees) = parseNewicks name_hack (zip files bstrs) + + -------------------------------------------------------------------------------- + + 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)]) + do_one (FullTree treename lblAcc parsed) = do + let + pruned = case branch_collapse_thresh of + Nothing -> parsed + Just thr -> collapseBranchLenThresh thr parsed + numL = numLeaves pruned + + -- 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], []) + + results <- mapM do_one fulltrees + let (counts::[Int], validtreess, pairs::[[(Int, String)]]) = unzip3 results + let validtrees = concat validtreess + warnings2 = concat pairs + + putStrLn$ "\nNumber of input tree files: " ++ show num_files + 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) + putStrLn$ "Total tree nodes contained in valid trees: "++ show (sum counts) + let all_branches = concatMap (F.foldr' (\x ls -> getBranchLen x:ls) []) validtrees + putStrLn$ "Average branch len over valid trees: "++ show (avg all_branches) + putStrLn$ "Max/Min branch lengths: "++ show (foldl1' max all_branches, + foldl1' min all_branches) + + -------------------------------------------------------------------------------- + -- Next, dispatch on the mode and do the actual clustering or binning. + -------------------------------------------------------------------------------- + + (classes,binlist,asyncs) <- case clust_mode of + BinThem -> do x <- doBins validtrees + -- A list of bins, sorted by size: + let binlist = reverse $ sortBy (compare `on` fst3) $ + map (\ (tr,OneCluster ls) -> (length ls, tr, OneCluster ls)) $ + M.toList x + return (x,binlist,[]) + ClusterThem lnk -> do + (mat, dendro) <- doCluster lnk validtrees + case print_rfmatrix of + False -> return () + True -> do -- treeFiles <- acquireTreeFiles files + -- let fn f = do raw <- B.readFile f + -- let ls = map (`B.append` (B.pack ";")) $ + -- B.splitWith (== ';') raw + -- return (map (f,) ls) + -- trees0 <- concat <$> mapM fn treeFiles + -- FIXME: no name_hack here: + -- let (lbls, trees) = parseNewicks id trees0 + -- putStrLn$ "Read trees! "++show (length trees) + -- putStrLn$ "Taxa: "++show (pPrint lbls) + -- putStrLn$ "First tree: "++show (displayDefaultTree (head trees)) + printDistMat stdout mat + writeFile (combine output_dir ("dendrogram.txt")) + (show$ fmap treename dendro) + putStrLn "Wrote full dendrogram to file dendrogram.txt" + + gvizAsync <- async $ do + t0 <- getCurrentTime + let dot = dotDendrogram "dendrogram" 1.0 dendro + _ <- dotToPDF dot (combine output_dir "dendrogram.pdf") + t1 <- getCurrentTime + putStrLn$ "Wrote dendrogram diagram to file dendrogram.pdf ("++show(diffUTCTime t1 t0)++")" + + hnd <- openFile (combine output_dir ("distance_matrix.txt")) WriteMode + printDistMat hnd mat + hClose hnd + + case dist_thresh of + Nothing -> error "Fully hierarchical cluster output is not finished! Use --editdist." + Just dstThresh -> do + putStrLn$ "Combining all clusters at distance less than or equal to "++show dstThresh + let clusts = sliceDendro (fromIntegral dstThresh) dendro + wlens = map (\ (OneCluster l) -> (length l, error "need consensus tree", OneCluster l)) clusts + sorted0 = reverse$ sortBy (compare `on` fst3) wlens +-- sorted = map thd3 sorted0 + putStrLn$ "After flattening, cluster sizes are: "++show (map fst3 sorted0) + -- Flatten out the dendogram: + return (clustsToMap clusts, sorted0, [gvizAsync]) + + reportClusts clust_mode binlist + + ---------------------------------------- + -- TEST, TEMPTOGGLE: print out edge weights : + -- forM_ (map snd3 results) $ \parsed -> do + -- let weights = all_edge_weights (head$ S.toList taxa) parsed + -- trace ("weights of "++ show parsed ++" "++ show weights) $ + -- return () + -- exitSuccess + ---------------------------------------- + + -------------------------------------------------------------------------------- + -- Finally, produce all the required outputs. + -------------------------------------------------------------------------------- + + putStrLn$ "\nTotal unique taxa ("++ show (M.size labelTab) ++"):\n"++ + show (nest 2 $ sep $ map text $ M.elems labelTab) + + putStrLn$ "Final number of tree bins: "++ show (M.size classes) + + unless (null warnings1 && null warnings2) $ + writeFile (combine output_dir "WARNINGS.txt") + ("This file was generated to record all of the files which WERE NOT incorporated successfully into the results.\n" ++ + "Each of these files had some kind of problem, likely one of the following:\n"++ + " (1) a mismatched number of taxa (leaves) in the tree relative to the rest of the dataset\n"++ + " (2) a file that could not be read.\n"++ + " (3) a file that could not be parsed.\n\n"++ + concat (map (\ file -> "Not a regular/readable file: "++ file++"\n") + warnings1) ++ + concat (map (\ (n,file) -> + "Wrong number of taxa ("++ show n ++"): "++ file++"\n") + warnings2)) + + case clust_mode of + BinThem -> outputBins binlist output_dir do_graph + ClusterThem lnk -> outputClusters binlist output_dir do_graph + + -- Wait on parallel tasks: + putStrLn$ "Waiting for asynchronous tasks to finish..." + mapM_ wait asyncs + putStrLn$ "Finished." + -------------------------------------------------------------------------------- + -- End driver + -------------------------------------------------------------------------------- + +-- filePrefix = "bin" +filePrefix = "cluster" + +-------------------------------------------------------------------------------- +-- Driver helpers: +-------------------------------------------------------------------------------- + +-- doBins :: [FullTree DefDecor] -> t1 -> t2 -> IO BinResults +doBins :: [FullTree DefDecor] -> IO (BinResults StandardDecor) +doBins validtrees = do + putStrLn$ "Creating equivalence classes (bins)..." + let classes = --binthem_normed$ zip files $ concat$ map snd3 results + binthem validtrees + return (classes) + +doCluster :: C.Linkage -> [FullTree a] -> IO (DistanceMatrix, C.Dendrogram (FullTree a)) +doCluster linkage validtrees = do + putStrLn$ "Clustering using method "++show linkage + let nwtrees = map nwtree validtrees + numtrees = length validtrees + mat = distanceMatrix nwtrees + ixtrees = zip [0..] validtrees + dist (i,t1) (j,t2) | j == i = 0 +-- | i == numtrees-1 = 0 + | j < i = fromIntegral ((mat V.! i) U.! j) + | otherwise = fromIntegral ((mat V.! j) U.! i) + dist1 a b = trace ("Taking distance between "++show (fst a, fst b)) $ dist a b + dendro = fmap snd $ C.dendrogram linkage ixtrees dist + -- putStrLn$ "Got numtrees ... "++show numtrees + -- putStrLn$ "Got the distance matrix ... "++show (V.length mat) + return (mat,dendro) + +-- reportClusts :: ClustMode -> BinResults StandardDecor -> IO [(Int, StrippedTree, OneCluster StandardDecor)] +reportClusts mode binlist = do + let + taxa :: S.Set Int + taxa = S.unions$ map (S.fromList . all_labels . snd3) binlist + binsizes = map fst3 binlist + + putStrLn$ " Outcome: "++show (length binlist)++" clusters found, "++show (length$ takeWhile (>1) binsizes) + ++" non-singleton, top bin sizes: "++show(take 10 binsizes) + putStrLn$" First 50 bin sizes, excluding singletons:" + forM_ (zip [1..50] binlist) $ \ (ind, (len, tr, OneCluster ftrees)) -> do + when (len > 1) $ -- Omit that long tail of single element classes... + putStrLn$show$ + hcat [text (" * cluster#"++show ind++", members "++ show len ++", "), + case mode of + BinThem -> vcat [text ("avg bootstraps "++show (get_bootstraps$ avg_trees$ map nwtree ftrees)++", "), + text "all: " <> pPrint (filter (not . null) $ map (get_bootstraps . nwtree) ftrees)] + ClusterThem _ -> hcat [] ] + return binlist + +-- | Convert a flat list of clusters into a map from individual trees to clusters. +clustsToMap :: [OneCluster StandardDecor] -> BinResults StandardDecor +clustsToMap clusts = F.foldl' fn M.empty clusts + where + fn acc theclust@(OneCluster ftrs) = + F.foldl' (fn2 theclust) acc ftrs + fn2 theclust acc (FullTree{nwtree}) = + M.insert (anonymize_annotated nwtree) theclust acc + +flattenDendro :: (C.Dendrogram (FullTree DefDecor)) -> OneCluster StandardDecor +flattenDendro dendro = + case dendro of + C.Leaf (FullTree{treename,labelTable,nwtree}) -> + OneCluster [FullTree treename labelTable (annotateWLabLists nwtree)] + C.Branch _ left right -> + -- TODO: fix quadratic append + flattenDendro left `appendClusts` flattenDendro right + where + appendClusts (OneCluster ls1) (OneCluster ls2) = OneCluster (ls1++ls2) + +-- | Turn a hierarchical clustering into a flat clustering. +sliceDendro :: Double -> (C.Dendrogram (FullTree DefDecor)) -> [OneCluster StandardDecor] +sliceDendro dstThresh den = loop den + where + loop br@(C.Branch dist left right) + -- Too far apart to combine: + | dist > dstThresh = loop left ++ loop right + | otherwise = [flattenDendro br] + loop br@(C.Leaf _) = [flattenDendro br] + + +-------------------------------------------------------------------------------- + +outputClusters binlist output_dir do_graph = do + let numbins = length binlist + let base i size = combine output_dir (filePrefix ++ show i ++"_"++ show size) + + forM_ (zip [1::Int ..] binlist) $ \ (i, (size, _tr, OneCluster ftrees)) -> do + writeFile (base i size ++".txt") (concat$ map ((++"\n") . treename) ftrees) + -- TODO: CONSENSUS TREE: + -- writeFile (base i size ++".tr") (show (displayDefaultTree$ deAnnotate fullAvgTr) ++ ";\n") -- FIXME + + putStrLn$ "[finished] Wrote contents of each cluster to cluster_.txt" + putStrLn$ " Wrote representative trees to cluster_.tr" + when (do_graph) $ do + putStrLn$ "Next do the time consuming operation of writing out graphviz visualizations:" + forM_ (zip [1::Int ..] binlist) $ \ (i, (size, _tr, OneCluster membs)) -> do + when (size > 1 || numbins < 100) $ do + -- TODO CONSENSUS TREE +-- let dot = dotNewickTree ("cluster #"++ show i) (1.0 / avg_branchlen (map nwtree membs)) fullAvgTr +-- _ <- dotToPDF dot (base i size ++ ".pdf") + return () + putStrLn$ "[finished] Wrote visual representations of trees to "++filePrefix++"_.pdf" + + return () + + +outputBins binlist output_dir do_graph = do + let numbins = length binlist + let base i size = combine output_dir (filePrefix ++ show i ++"_"++ show size) + let avgs = map (avg_trees . map nwtree . clustMembers . thd3) binlist + forM_ (zip3 [1::Int ..] binlist avgs) $ \ (i, (size, _tr, OneCluster ftrees), avgTree) -> do + let FullTree fstName labs _ = head ftrees + fullAvgTr = FullTree fstName labs avgTree + + --putStrLn$ (" WRITING " ++ combine output_dir (filePrefix ++ show i ++"_"++ show size ++".txt")) + writeFile (base i size ++".txt") (concat$ map ((++"\n") . treename) ftrees) + -- writeFile (base i size ++".tr") (show (pPrint tr) ++ ";\n") + -- Printing the average tree instead of the stripped cannonical one: + when debug$ do + writeFile (base i size ++".dbg") (show (pPrint avgTree) ++ "\n") + writeFile (base i size ++".tr") (show (displayDefaultTree$ deAnnotate fullAvgTr) ++ ";\n") -- FIXME + + putStrLn$ "[finished] Wrote contents of each bin to "++filePrefix++"_.txt" + putStrLn$ " Wrote representative trees to "++filePrefix++"_.tr" + when (do_graph) $ do + putStrLn$ "Next do the time consuming operation of writing out graphviz visualizations:" + forM_ (zip3 [1::Int ..] binlist avgs) $ \ (i, (size, _tr, OneCluster membs), avgTree) -> do + let FullTree fstName labs _ = head membs + fullAvgTr = FullTree fstName labs avgTree + when (size > 1 || numbins < 100) $ do + let dot = dotNewickTree ("cluster #"++ show i) (1.0 / avg_branchlen (map nwtree membs)) + --(annotateWLabLists$ fmap (const 0) tr) + -- TEMP FIXME -- using just ONE representative tree: + ( --trace ("WEIGHTED: "++ show (head$ trees bentry)) $ + --(head$ trees bentry) + fullAvgTr) + _ <- dotToPDF dot (base i size ++ ".pdf") + return () + putStrLn$ "[finished] Wrote visual representations of trees to bin_.pdf" + + +-------------------------------------------------------------------------------- + +-- Monadic mapAccum +mapAccumM :: Monad m => (acc -> x -> m (acc,y)) -> acc -> [x] -> m (acc,[y]) +mapAccumM fn acc ls = F.foldrM fn' (acc,[]) ls + where + fn' x (acc,ls) = do (acc',y) <- fn acc x + return (acc',y:ls) + +fst3 :: (t, t1, t2) -> t +snd3 :: (t, t1, t2) -> t1 +thd3 :: (t, t1, t2) -> t2 +fst3 (a,_,_) = a +snd3 (_,b,_) = b +thd3 (_,_,c) = c + + +-- Come up with an average tree from a list of isomorphic trees. +-- This comes up with some blending of edge lengths. +avg_trees :: [AnnotatedTree] -> AnnotatedTree +avg_trees origls = + fmap unlift $ + foldIsomorphicTrees (foldl1 sumthem) $ + map (fmap lift) origls + where + totalCount = fromIntegral$ length origls + + -- Here we do the actual averaging: + unlift :: TempDecor -> StandardDecor + unlift (bl, (bs, bootcount), wt, ls) = + (StandardDecor (bl / totalCount) + (case bootcount of + 0 -> Nothing + _ -> Just$ round (fromIntegral bs / fromIntegral bootcount)) + wt ls) + + lift :: StandardDecor -> TempDecor + lift (StandardDecor bl bs wt ls) = (bl, (fromMaybe 0 bs, countMayb bs), wt, ls) + + sumthem :: TempDecor -> TempDecor -> TempDecor + sumthem (bl1, (bs1, cnt1), wt, ls) + (bl2, (bs2, cnt2), _, _) = + (bl1+bl2, (bs1 + bs2, cnt1+cnt2), wt, ls) + countMayb Nothing = 0 + countMayb (Just _) = 1 + +-- Used only by avg_trees above... +type TempDecor = (Double, (Int, Int), Int, [Label]) + +avg ls = sum ls / fromIntegral (length ls) diff --git a/Bio/Phylogeny/PhyBin/Binning.hs b/Bio/Phylogeny/PhyBin/Binning.hs new file mode 100644 --- /dev/null +++ b/Bio/Phylogeny/PhyBin/Binning.hs @@ -0,0 +1,432 @@ +{-# LANGUAGE ScopedTypeVariables, RecordWildCards, TypeSynonymInstances, CPP #-} +{-# LANGUAGE NamedFieldPuns #-} +{-# LANGUAGE OverloadedStrings #-} +{-# OPTIONS_GHC -fwarn-incomplete-patterns #-} +{-# OPTIONS_GHC -fwarn-unused-imports #-} + +-- | This module contains the code that does the tree normalization and binning. + +module Bio.Phylogeny.PhyBin.Binning + ( -- * Binning and normalization + binthem, normalize, normalizeFT, annotateWLabLists, + deAnnotate, OneCluster(..), BinResults, StrippedTree, + -- * Utilities and unit tests + get_weight, unitTests, anonymize_annotated + ) + 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.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.Process (system) +import System.Exit (ExitCode(..)) +import Test.HUnit ((~:),(~=?),Test,test) +import qualified HSH + +-- 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 + +-- Turn on for extra invariant checking: +debug :: Bool +debug = True + +---------------------------------------------------------------------------------------------------- +-- Normal form for unordered, unrooted trees +---------------------------------------------------------------------------------------------------- + +-- The basic idea is that what we *want* is the following, +-- ROOT: most balanced point +-- ORDER: sorted in increasing subtree weight + +-- But that's not quite good enough. There are ties to break. To do +-- that we fall back on the (totally ordered) leaf labels. + +-------------------------------------------------------------------------------- + + +-- Our sorting criteria for the children of interior nodes: +compare_childtrees :: AnnotatedTree -> AnnotatedTree -> Ordering +compare_childtrees node1 node2 = + case (subtreeWeight $ get_dec node1) `compare` (subtreeWeight $ get_dec node2) of + -- Comparisons on atoms cause problems WRT to determinism between runs if parallelism is introduced. + -- Can consider it an optimization for the serial case perhaps: +-- EQ -> case map deAtom (get_label_list node1) `compare` +-- map deAtom (get_label_list node2) of + EQ -> case (sortedLabels$ get_dec node1) `compare` (sortedLabels$ get_dec node2) of + EQ -> error$ "Internal invariant broken. These two children have equal ordering priority:\n" + ++ "Pretty printing:\n " + ++ show (pPrint node1) ++ "\n " ++ show (pPrint node2) + ++ "\nFull data structures:\n " + ++ show (node1) ++ "\n " ++ show (node2) + x -> x + x -> x + +-- | A version lifted to operate over full trees. +normalizeFT :: FullTree StandardDecor -> FullTree StandardDecor +normalizeFT (FullTree nm labs tr) = FullTree nm labs (normalize tr) + +-- | This is it, here's the routine that transforms a tree into normal form. +-- This relies HEAVILY on lazy evaluation. +normalize :: AnnotatedTree -> AnnotatedTree +normalize tree = snd$ loop tree Nothing + where + + add_context dec Nothing = dec + add_context dec (Just c) = add_weight dec c + + -- loop: Walk over the tree, turning it inside-out in the process. + -- Inputs: + -- 1. node: the NewickTree node to process ("us") + -- 2. context: all nodes connected through the parent, "flipped" as though *we* were root + -- The "flipped" part has ALREADY been normalized. + -- Outputs: + -- 1. new node + -- 3. the best candidate root anywhere under this subtree + loop :: AnnotatedTree -> Maybe (AnnotatedTree) -> (AnnotatedTree, AnnotatedTree) + loop node ctxt = case node of + NTLeaf (StandardDecor _ _ w sorted) _name -> + (node, + -- If the leaf becomes the root... we could introduce another node: + NTInterior (add_context (StandardDecor 0 Nothing w sorted) ctxt) $ + (verify_sorted "1" id$ maybeInsert compare_childtrees ctxt [node]) + + -- It may be reasonable to not support leaves becoming root.. that changes the number of nodes! + --error "normalize: leaf becoming root not currently supported." + ) + + NTInterior dec ls -> + let + -- If this node becomes the root, the parent becomes one of our children: + inverted = NTInterior inverted_dec inverted_children + inverted_dec = add_context dec ctxt + inverted_children = verify_sorted "2" id$ maybeInsert compare_childtrees ctxt newchildren + + newchildren = --trace ("SORTED "++ show (map (get_label_list . fst) sorted)) $ + map fst sorted + sorted = sortBy (compare_childtrees `on` fst) possibs + + possibs = + flip map ls $ \ child -> + let + + -- Will this diverge??? Probably depends on how equality (for delete) is defined... + + -- Reconstruct the current node missing one child (because it became a parent): + -- Update its metadata appropriately: + newinverted = NTInterior (subtract_weight inverted_dec child) + (verify_sorted "3" id$ delete newnode inverted_children) + (newnode, _) = result + + result = loop child (Just newinverted) + in + result + + -- Either us or a candidate suggested by one of the children: + rootcandidates = inverted : map snd sorted + + -- Who wins? The "most balanced". Minimize max subtree weight. + -- The compare operator is NOT allowed to return EQ here. Therefore there will be a unique minima. + winner = --trace ("Candidates: \n"++ show (nest 6$ vcat (map pPrint (zip (map max_subtree_weight rootcandidates) rootcandidates )))) $ + minimumBy cmpr_subtree_weight rootcandidates + + max_subtree_weight = maximum . map get_weight . get_children + fat_id = map get_label_list . get_children + + cmpr_subtree_weight tr1 tr2 = + case max_subtree_weight tr1 `compare` max_subtree_weight tr2 of + EQ -> -- As a fallback we compare the alphabetic order of the "bignames" of the children: + case fat_id tr1 `compare` fat_id tr2 of + EQ -> error$ "\nInternal invariant broken. These two were equally good roots:\n" + ++ show (pPrint tr1) ++ "\n" ++ show (pPrint tr2) + x -> x + x -> x + + in (NTInterior dec newchildren, winner) + + +-- Verify that our invariants are met: +verify_sorted :: (Show a, Pretty a) => String -> (a -> AnnotatedTree) -> [a] -> [a] +verify_sorted msg = + if debug + then \ project nodes -> + let weights = map (get_weight . project) nodes in + if sort weights == weights + then nodes +-- else error$ "Child list failed verification: "++ show (pPrint nodes) + else error$ msg ++ ": Child list failed verification, not sorted: "++ show (weights) + ++"\n "++ show (sep $ map pPrint nodes) ++ + "\n\nFull output:\n " ++ (concat$ intersperse "\n " $ map show nodes) + else \ _ nodes -> nodes + + +-- TODO: Salvage any of these tests that are worthwhile and get them into the unit tests: +tt :: AnnotatedTree +tt = normalize $ annotateWLabLists $ snd$ parseNewick M.empty id "" "(A,(C,D,E),B);" + +norm4 :: FullTree StandardDecor +norm4 = norm "((C,D,E),B,A);" + +norm5 :: AnnotatedTree +norm5 = normalize$ annotateWLabLists$ snd$ parseNewick M.empty id "" "(D,E,C,(B,A));" + + +---------------------------------------------------------------------------------------------------- +-- Equivalence classes on Trees: +---------------------------------------------------------------------------------------------------- + +-- | When binning, the members of a OneCluster are isomorphic trees. When clustering +-- based on robinson-foulds distance they are merely similar trees. +newtype OneCluster a = OneCluster { clustMembers :: [FullTree a] } + deriving Show + +-- | Ignore metadata (but keep weights) for the purpose of binning +type StrippedTree = NewickTree Int + +-- | Index the results of binning by topology-only stripped trees +-- that have their decorations removed. +type BinResults a = M.Map StrippedTree (OneCluster a) + +-- | The binning function. +-- Takes labeled trees, classifies labels into equivalence classes. +binthem :: [FullTree DefDecor] -> BinResults StandardDecor +binthem ls = binthem_normed normalized + where + normalized = map (\ (FullTree n lab tree) -> + FullTree n lab (normalize $ annotateWLabLists tree)) ls + +-- | This version accepts trees that are already normalized: +binthem_normed :: [FullTree StandardDecor] -> BinResults StandardDecor +binthem_normed normalized = +-- foldl (\ acc (lab,tree) -> M.insertWith update tree (OneCluster{ members=[lab] }) acc) + foldl (\ acc ft@(FullTree treename _ tree) -> + M.insertWith update (anonymize_annotated tree) (OneCluster [ft]) acc) + M.empty normalized + --(map (mapSnd$ fmap (const ())) normalized) -- still need to STRIP them + where +-- update new old = OneCluster{ members= (members new ++ members old) } + update (OneCluster new) (OneCluster old) = OneCluster (new++old) + --strip = fmap (const ()) + +-- | For binning. Remove branch lengths and labels but leave weights. +anonymize_annotated :: AnnotatedTree -> StrippedTree +anonymize_annotated = fmap (\ (StandardDecor bl bs w labs) -> w) + + +---------------------------------------------------------------------------------------------------- +-- Other tools and algorithms. +---------------------------------------------------------------------------------------------------- + +-- Extract all edges connected to a particular node in every tree. Return branch lengths. +all_edge_weights lab trees = + concat$ map (loop []) trees + where + loop acc (NTLeaf len name) | lab == name = len:acc + loop acc (NTLeaf _ _) = acc + loop acc (NTInterior _ ls) = foldl loop acc ls + + +---------------------------------------------------------------------------------------------------- +-- Bitvector based normalization. +---------------------------------------------------------------------------------------------------- + +-- TODO: This approach is probably faster. Give it a try. + +{- +int NumberOfSetBits(int i) +{ + i = i - ((i >> 1) & 0x55555555); + i = (i & 0x33333333) + ((i >> 2) & 0x33333333); + return ((i + (i >> 4) & 0xF0F0F0F) * 0x1010101) >> 24; +} + +int __builtin_popcount (unsigned int x); +-} + + +---------------------------------------------------------------------------------------------------- + + + +{- + ---------------------------------------- + PARSING TIMING TEST: + ---------------------------------------- + + Compiling this with GHC 6.12 on my laptop -O2... + It takes 0.043s startup to parse ten files. + And 0.316 seconds to parse 2648.. so we can say that's almost all time spent parsing/building/traversing. + (All nodes summed to 14966) + (The tested version uses Strings for labels... not Atoms) + + Comparing against the original mzscheme version (with Racket 5.0) + with default optimization (there's no obvious -O2), well the + generated .exe has a ~0.5 second startup time overhead... + 0.881 seconds total to do the parsing, or about 380ms just for parsing. + But that doesn't do the counting! + Ok, this mzscheme version is in a messed up state at this point, but hacking + it to do a count (and it gets a different one.. 12319), I get 0.882 seconds real time, + that is neglibly more. + + If anything parsec should be at a disadvantage because of the lack of + a preprocessing phase to generate the FSM... + + Btw, switching node labels over to Atoms made no difference. (But + didn't slow down at least.) We wouldn't expect this to save anything + on the construction side... parsec still allocates/builds the strings + before we intern them. + + -} + + + + +---------------------------------------------------------------------------------------------------- +-- General helper/utility functions: +---------------------------------------------------------------------------------------------------- + +merge :: Ord a => [a] -> [a] -> [a] +merge [] ls = ls +merge ls [] = ls +merge l@(a:b) r@(x:y) = + if a < x + then a : merge b r + else x : merge y l + +-- Set subtraction for sorted lists: +demerge :: (Ord a, Show a) => [a] -> [a] -> [a] +demerge ls [] = ls +demerge [] ls = error$ "demerge: first list did not contain all of second, remaining: " ++ show ls +demerge (a:b) r@(x:y) = + case a `compare` x of + EQ -> demerge b y + LT -> a : demerge b r + GT -> error$ "demerge: element was missing from first list: "++ show x + +-- maybeCons :: Maybe a -> [a] -> [a] +-- maybeCons Nothing ls = ls +-- maybeCons (Just x) ls = x : ls + +maybeInsert :: (a -> a -> Ordering) -> Maybe a -> [a] -> [a] +maybeInsert _ Nothing ls = ls +maybeInsert fn (Just x) ls = insertBy fn x ls + +-- | Add the metadata that is used for binning +annotateWLabLists :: NewickTree DefDecor -> AnnotatedTree +annotateWLabLists tr = case tr of + NTLeaf (bs,bl) n -> NTLeaf (StandardDecor bl bs 1 [n]) n + NTInterior (bs,bl) ls -> + let children = map annotateWLabLists ls in + NTInterior (StandardDecor bl bs + (sum $ map (subtreeWeight . get_dec) children) + (foldl1 merge $ map (sortedLabels . get_dec) children)) + children + +---------------------------------------------------------------------------------------------------- +-- Simple Helper Functions +---------------------------------------------------------------------------------------------------- + +-- | Take the extra annotations away. Inverse of `annotateWLabLists`. +deAnnotate :: FullTree StandardDecor -> FullTree DefDecor +deAnnotate (FullTree a b tr) = FullTree a b (fmap (\ (StandardDecor bl bs _ _) -> (bs,bl)) tr) + +-- Number of LEAVES contained in subtree: +get_weight :: AnnotatedTree -> Int +get_weight = subtreeWeight . get_dec + +-- Sorted list of leaf labels contained in subtree +get_label_list :: AnnotatedTree -> [Label] +get_label_list = sortedLabels . get_dec + +add_weight :: StandardDecor -> AnnotatedTree -> StandardDecor +add_weight (StandardDecor l1 bs1 w1 sorted1) node = + let (StandardDecor _ bs2 w2 sorted2) = get_dec node in + StandardDecor l1 ((+) <$> bs1 <*> bs2) (w1+w2) (merge sorted1 sorted2) + +-- Remove the influence of one subtree from the metadata of another. +subtract_weight :: StandardDecor -> AnnotatedTree -> StandardDecor +subtract_weight (StandardDecor l1 bs1 w1 sorted1) node = + let (StandardDecor _ bs2 w2 sorted2) = get_dec node in + StandardDecor l1 ((-) <$> bs1 <*> bs2) (w1-w2) (demerge sorted1 sorted2) + + +---------------------------------------------------------------------------------------------------- +-- UNIT TESTING +---------------------------------------------------------------------------------------------------- + +tre1 :: (LabelTable, NewickTree DefDecor) +tre1 = parseNewick M.empty id "" "(A:0.1,B:0.2,(C:0.3,D:0.4):0.5);" + +tre1draw :: IO (Chan (), FullTree StandardDecor) +tre1draw = viewNewickTree "tre1"$ (FullTree "" (fst tre1) (annotateWLabLists (snd tre1))) + +tre1dot :: IO () +tre1dot = print $ dotNewickTree "" 1.0 $ (FullTree "" (fst tre1) (annotateWLabLists$ snd tre1)) + +norm :: String -> FullTree StandardDecor +norm = norm2 . B.pack + +norm2 :: B.ByteString -> FullTree StandardDecor +norm2 bstr = FullTree "" tbl (normalize $ annotateWLabLists tr) + where + (tbl,tr) = parseNewick M.empty id "test" bstr + +unitTests :: Test +unitTests = + let ntl s = NTLeaf (Nothing,0.0) s + (tbl,tre1_) = tre1 + in + test + [ + "merge" ~: [1,2,3,4,5,6] ~=? merge [1,3,5] [2,4,6::Int] + , "demerge" ~: [2,4,6] ~=? demerge [1,2,3,4,5,6] [1,3,5::Int] + , "demerge" ~: [1,3,5] ~=? demerge [1,2,3,4,5,6] [2,4,6::Int] + , "annotateWLabLists" ~: + ["A","B","C","D"] -- ORD on atoms is expensive... it must use the whole string. + ~=? map (tbl M.!) + (sortedLabels (get_dec (annotateWLabLists tre1_))) + + -- 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));") + + -- Here's an example from the rhizobia datasetsthat that caused my branch-sorting to fail. + , "normalize3" ~: "(((BB, BJ)), (MB, ML), (RE, (SD, SM)))" + ~=? show (pPrint$ norm2 (B.pack "(((ML,MB),(RE,(SD,SM))),(BB,BJ));")) + +-- "((BB: 2.691831, BJ: 1.179707): 0.000000, ((ML: 0.952401, MB: 1.020319): 0.000000, (RE: 2.031345, (SD: 0.180786, SM: 0.059988): 0.861187): 0.717913): 0.000000);" + + , "phbin: these 3 trees should fall in the same category" ~: + 1 ~=? (length $ M.toList $ + binthem $ snd $ + parseNewicks id [("one", "(A,(C,D,E),B);"), + ("two", "((C,D,E),B,A);"), + ("three","(D,E,C,(B,A));")] + -- [("one", snd$parseNewick M.empty id "" "(A,(C,D,E),B);"), + -- ("two", snd$parseNewick M.empty id "" "((C,D,E),B,A);"), + -- ("three", snd$parseNewick M.empty id "" "(D,E,C,(B,A));")] + ) + + , "dotConversion" ~: True ~=? 100 < length (show $ dotNewickTree "" 1.0$ norm "(D,E,C,(B,A));") -- 444 + ] + + +-------------------------------------------------------------------------------- diff --git a/Bio/Phylogeny/PhyBin/CoreTypes.hs b/Bio/Phylogeny/PhyBin/CoreTypes.hs new file mode 100644 --- /dev/null +++ b/Bio/Phylogeny/PhyBin/CoreTypes.hs @@ -0,0 +1,321 @@ +{-# LANGUAGE NamedFieldPuns, BangPatterns #-} +{-# LANGUAGE OverloadedStrings, TypeSynonymInstances, FlexibleInstances #-} +{-# LANGUAGE CPP #-} +{-# OPTIONS_GHC -fwarn-incomplete-patterns #-} +{-# OPTIONS_GHC -fwarn-unused-imports #-} + +module Bio.Phylogeny.PhyBin.CoreTypes + ( + -- * Tree and tree decoration types + NewickTree(..), + DefDecor, StandardDecor(..), AnnotatedTree, FullTree(..), + ClustMode(..), TreeName, + + -- * Tree operations + displayDefaultTree, + treeSize, numLeaves, liftFT, + get_dec, set_dec, get_children, + map_labels, all_labels, foldIsomorphicTrees, + + -- * Utilities specific to StandardDecor: + avg_branchlen, get_bootstraps, + + -- * Command line config options + PhyBinConfig(..), default_phybin_config, + + -- * General helpers + Label, LabelTable, + + -- * Experimenting with abstracting decoration operations + HasBranchLen(..) + ) + where + +import qualified Data.Map as M +import Data.Foldable (Foldable(..)) +import Data.Maybe (maybeToList) +import Data.Monoid (mappend, mconcat) +import Text.PrettyPrint.HughesPJClass hiding (char, Style) + +import qualified Data.Clustering.Hierarchical as C + +#define NO_ATOMS +#ifndef NO_ATOMS +import StringTable.Atom +#endif + +---------------------------------------------------------------------------------------------------- +-- Type definitions +---------------------------------------------------------------------------------------------------- + +type BranchLen = Double + +-- | Even though the Newick format allows it, here we ignore interior node +-- labels. (They are not commonly used.) +-- +-- Note that these trees are rooted. The `normalize` function ensures that a +-- single, canonical rooted representation is chosen. +data NewickTree a = + NTLeaf a {-# UNPACK #-} !Label + | NTInterior a [NewickTree a] + deriving (Show, Eq, Ord) + +-- TODO: Ordering maybe shouldn't need to touch the metadata. At least on the fast +-- path. + +{- +-- [2010.09.22] Disabling: +instance NFData Atom where + rnf a = rnf (fromAtom a :: Int) + +instance NFData a => NFData (NewickTree a) where + rnf (NTLeaf l n) = rnf (l,n) + rnf (NTInterior l ls) = rnf (l,ls) +-} + +instance Functor NewickTree where + fmap fn (NTLeaf dec x) = NTLeaf (fn dec) x + fmap fn (NTInterior dec ls) = NTInterior (fn dec) (map (fmap fn) ls) + +instance Foldable NewickTree where + foldMap f (NTLeaf dec x) = f dec + foldMap f (NTInterior dec ls) = mappend (f dec) $ + mconcat (map (foldMap f) ls) + +instance Foldable FullTree where + foldMap f (FullTree _ _ tr) = foldMap f tr + +instance Pretty dec => Pretty (NewickTree dec) where + -- I'm using displayDefaultTree for the "prettiest" printing and + -- replacing pPrint with a whitespace-improved version of show: + pPrint (NTLeaf dec name) = "NTLeaf" <+> pPrint dec <+> text (show name) + pPrint (NTInterior dec ls) = "NTInterior" <+> pPrint dec <+> pPrint ls + +instance Pretty a => Pretty (FullTree a) where + pPrint (FullTree name mp tr) = + "FullTree " <+> text name <+> loop tr + where + loop (NTLeaf dec ind) = "NTLeaf" <+> pPrint dec <+> text (mp M.! ind) + loop (NTInterior dec ls) = "NTInterior" <+> pPrint dec <+> pPrint ls + +instance (Pretty k, Pretty v) => Pretty (M.Map k v) where + pPrint mp = pPrint (M.toList mp) + +-- | Display a tree WITH the bootstrap and branch lengths. +displayDefaultTree :: FullTree DefDecor -> Doc +displayDefaultTree orig = loop tr <> ";" + where + (FullTree _ mp tr) = orig -- normalize orig + loop (NTLeaf (Nothing,_) name) = text (mp M.! name) + loop (NTLeaf _ _) = error "WEIRD -- why did a leaf node have a bootstrap value?" + loop (NTInterior (bootstrap,_) ls) = + case bootstrap of + Nothing -> base + Just val -> base <> text ":[" <> text (show val) <> text "]" + where base = parens$ sep$ map_but_last (<>text",") $ map loop ls + +---------------------------------------------------------------------------------------------------- +-- 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] +---------------------------------------- +type Label = Int + +-- | Map labels back onto meaningful names. +type LabelTable = M.Map Label String + +---------------------------------------------------------------------------------------------------- +-- Tree metadata (decorators) +---------------------------------------------------------------------------------------------------- + + +-- | The barebones default decorator for NewickTrees contains BOOTSTRAP and +-- BRANCHLENGTH. The bootstrap values, if present, will range in [0..100] +type DefDecor = (Maybe Int, BranchLen) + +-- | Additionally includes some scratch data that is used by the binning algorithm. +type AnnotatedTree = NewickTree StandardDecor + +-- | The standard decoration includes everything in `DefDecor` plus +-- some extra cached data: +-- +-- (1) branch length from parent to "this" node +-- (2) bootstrap values for the node +-- +-- (3) subtree weights for future use +-- (defined as number of LEAVES, not counting intermediate nodes) +-- (4) sorted lists of labels for symmetry breaking +data StandardDecor = StandardDecor { + branchLen :: BranchLen, + bootStrap :: Maybe Int, + + -- The rest of these are used by the computations below. These are + -- cached (memoized) values that could be recomputed: + ---------------------------------------- + subtreeWeight :: Int, + sortedLabels :: [Label] + } + deriving (Show,Read,Eq,Ord) + +class HasBranchLen a where + getBranchLen :: a -> BranchLen + +instance HasBranchLen StandardDecor where + getBranchLen = branchLen + +-- This one is kind of sloppy: +instance HasBranchLen DefDecor where + getBranchLen = snd + +-- | A common type of tree contains the standard decorator and also a table for +-- restoring the human-readable node names. +data FullTree a = + FullTree { treename :: TreeName + , labelTable :: LabelTable + , nwtree :: NewickTree a + } + deriving (Show) + +liftFT :: (NewickTree t -> NewickTree a) -> FullTree t -> FullTree a +liftFT fn (FullTree nm labs x) = FullTree nm labs (fn x) + +type TreeName = String + +instance Pretty StandardDecor where + pPrint (StandardDecor bl bs wt ls) = parens$ + "StandardDecor" <+> hsep [pPrint bl, pPrint bs +-- , pPrint wt, pPrint ls + ] + + +---------------------------------------------------------------------------------------------------- +-- * Configuring and running the command line tool. +---------------------------------------------------------------------------------------------------- + +-- | Due to the number of configuration options for the driver, we pack them into a record. +data PhyBinConfig = + PBC { verbose :: Bool + , num_taxa :: Int + , name_hack :: String -> String + , output_dir :: String + , inputs :: [String] + , do_graph :: Bool + , do_draw :: Bool + , clust_mode :: ClustMode + , print_rfmatrix :: Bool + , dist_thresh :: Maybe Int + , branch_collapse_thresh :: Maybe Double -- ^ Branches less than this length are collapsed. + } + +-- | 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.)" + , name_hack = id -- Default, no transformation of leaf-labels + , output_dir = "./phybin_out/" + , inputs = [] + , do_graph = False + , do_draw = False + , clust_mode = BinThem + , print_rfmatrix = False + , dist_thresh = Nothing + , branch_collapse_thresh = Nothing + } + +data ClustMode = BinThem | ClusterThem C.Linkage + +---------------------------------------------------------------------------------------------------- +-- * Simple utility functions for the core types: +---------------------------------------------------------------------------------------------------- + +-- | How many nodes (leaves and interior) are contained in a NewickTree? +treeSize :: NewickTree a -> Int +treeSize (NTLeaf _ _) = 1 +treeSize (NTInterior _ ls) = 1 + sum (map treeSize ls) + +-- | This counts only leaf nodes, which should include all taxa. +numLeaves :: NewickTree a -> Int +numLeaves (NTLeaf _ _) = 1 +numLeaves (NTInterior _ ls) = sum (map numLeaves ls) + + +map_but_last :: (a -> a) -> [a] -> [a] +map_but_last _ [] = [] +map_but_last _ [h] = [h] +map_but_last fn (h:t) = fn h : map_but_last fn t + + + +get_dec :: NewickTree t -> t +get_dec (NTLeaf dec _) = dec +get_dec (NTInterior dec _) = dec + +-- Set all the decorations to a constant: +set_dec :: b -> NewickTree a -> NewickTree b +set_dec d = fmap (const d) +--set_dec d (NTLeaf _ x) = NTLeaf d x +--set_dec d (NTInterior _ ls) = NTInterior d $ map (set_dec d) ls + +get_children :: NewickTree t -> [NewickTree t] +get_children (NTLeaf _ _) = [] +get_children (NTInterior _ ls) = ls + + +-- | Average branch length across all branches in all all trees. +avg_branchlen :: HasBranchLen a => [NewickTree a] -> Double +avg_branchlen origls = fst total / snd total + where + total = sum_ls $ map sum_tree origls + sum_ls ls = (sum$ map fst ls, sum$ map snd ls) + sum_tree (NTLeaf dec _) | getBranchLen dec == 0 = (0,0) + | otherwise = (abs (getBranchLen dec),1) + sum_tree (NTInterior dec ls) = + let branchLen = getBranchLen dec + (x,y) = sum_ls$ map sum_tree ls in + if branchLen == 0 then (x, y) else ((abs branchLen) + x, 1+y) + +-- | Retrieve all the bootstraps values actually present in a tree. +get_bootstraps :: NewickTree StandardDecor -> [Int] +get_bootstraps (NTLeaf (StandardDecor{bootStrap}) _) = maybeToList bootStrap +get_bootstraps (NTInterior (StandardDecor{bootStrap}) ls) = + maybeToList bootStrap ++ concatMap get_bootstraps ls + +-- | Apply a function to all the *labels* (leaf names) in a tree. +map_labels :: (Label -> Label) -> NewickTree a -> NewickTree a +map_labels fn (NTLeaf dec lbl) = NTLeaf dec $ fn lbl +map_labels fn (NTInterior dec ls) = NTInterior dec$ map (map_labels fn) ls + +-- | Return all the labels contained in the tree. +all_labels :: NewickTree t -> [Label] +all_labels (NTLeaf _ lbl) = [lbl] +all_labels (NTInterior _ ls) = concat$ map all_labels ls + + + +-- | This function allows one to collapse multiple trees while looking +-- only at the "horizontal slice" of all the annotations *at a given +-- position* in the tree. +-- +-- "Isomorphic" must apply both to the shape and the name labels or it +-- is an error to apply this function. +foldIsomorphicTrees :: ([a] -> b) -> [NewickTree a] -> NewickTree b +foldIsomorphicTrees _ [] = error "foldIsomorphicTrees: empty list of input trees" +foldIsomorphicTrees fn ls@(hd:_) = fmap fn horiztrees + where + -- Preserve the input order: + horiztrees = Prelude.foldr consTrees (fmap (const []) hd) ls + -- We use the tree datatype itself as the intermediate data + -- structure. This is VERY allocation-expensive, it would be + -- possible to trade compute for allocation here: + consTrees a b = case (a,b) of + (NTLeaf dec nm1, NTLeaf decls nm2) | nm1 /= nm2 -> error$"foldIsomorphicTrees: mismatched names: "++show (nm1,nm2) + | otherwise -> + NTLeaf (dec : decls) nm1 + (NTInterior dec ls1, NTInterior decls ls2) -> + NTInterior (dec:decls) $ zipWith consTrees ls1 ls2 + _ -> error "foldIsomorphicTrees: difference in tree shapes" + diff --git a/Bio/Phylogeny/PhyBin/Parser.hs b/Bio/Phylogeny/PhyBin/Parser.hs new file mode 100644 --- /dev/null +++ b/Bio/Phylogeny/PhyBin/Parser.hs @@ -0,0 +1,217 @@ +{-# LANGUAGE ScopedTypeVariables, BangPatterns, ParallelListComp #-} +{-# OPTIONS_GHC -fwarn-incomplete-patterns #-} +{-# OPTIONS_GHC -fwarn-unused-imports #-} + +module Bio.Phylogeny.PhyBin.Parser + (newick_parser, parseNewick, parseNewicks, parseNewickFiles, unitTests) + where +import Control.Exception (evaluate, handle, SomeException) +import qualified Data.ByteString.Lazy.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 Test.HUnit ((~:),(~=?),Test,test,assertFailure) +import System.FilePath (takeBaseName) + +import Bio.Phylogeny.PhyBin.CoreTypes +import Prelude as P + +type NameHack = (String->String) + +-- | 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. +parseNewick :: LabelTable -> NameHack -> String -> B.ByteString -> (LabelTable, NewickTree DefDecor) +parseNewick tbl0 name_hack file input = + extractLabelTable tbl0 $ + runB file (newick_parser name_hack) $ + B.filter (not . isSpace) input + +-- | 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 + bss <- mapM B.readFile nms + return $ parseNewicks fn (zip nms bss) + +-- | A version which takes in-memory trees as ByteStrings. +parseNewicks :: NameHack -> [(String,B.ByteString)] -> (LabelTable, [FullTree DefDecor]) +parseNewicks name_hack pairs = (labtbl, fullTrs) + where + fullTrs = [ FullTree (takeBaseName file) labtbl tr + | (file,_) <- pairs + | tr <- trs ] + (labtbl, trs) = P.foldr fn (M.empty,[]) pairs + fn (file,bstr) (!acc,!ls) = + let (acc',tr) = parseNewick acc name_hack file bstr + in (acc', tr:ls) + +runB :: Show a => String -> Parser a -> B.ByteString -> a +runB file p input = case (parse p "" input) of + Left err -> error ("parse error in file "++ show file ++" at "++ show err) + Right x -> x + +extractLabelTable :: LabelTable -> TempTree -> (LabelTable, NewickTree DefDecor) +extractLabelTable tbl0 tr = (finMap,finTree) + where + flipped = M.fromList $ L.map (\(x,y)->(y,x)) $ M.toList tbl0 + -- (_,finMap,finTree) = loop (S.fromList (M.elems tbl0)) tbl0 tr + (_,finMap,finTree) = loop flipped tbl0 tr + + loop seen acc (NTLeaf (d,Just nm) _) + | M.member nm seen = (seen, acc, NTLeaf d (seen M.! nm)) + | otherwise = let nxt = M.size acc in + (M.insert nm nxt seen, + M.insert nxt nm acc, NTLeaf d nxt) + loop seen1 acc1 (NTInterior (d,Nothing) chlds) = + let (seen',acc',ls') = + P.foldr (\ x (seen2,acc2,ls) -> + let (seen3,acc3,x') = loop seen2 acc2 x in + (seen3, acc3, x':ls)) + (seen1,acc1,[]) + chlds + in (seen',acc', NTInterior d ls') + +---------------------------------------------------------------------------------------------------- +-- Newick file format parser definitions: +---------------------------------------------------------------------------------------------------- + +-- | Hack: we store the names in the leaves. +type TempTree = NewickTree (DefDecor,Maybe String) + +-- | This is used to post-facto splice metadata into the data structure. +tag :: DefDecor -> TempTree -> TempTree +tag l s = + case s of + NTLeaf (_,m) n -> NTLeaf (l,m) n + NTInterior (_,Nothing) ls -> NTInterior (l,Nothing) ls + +-- | This parser ASSUMES that whitespace has been prefiltered from the input. +newick_parser :: NameHack -> Parser TempTree +newick_parser name_hack = + do x <- subtree name_hack + -- Get the top-level metadata: + l <- branchMetadat name_hack + _ <- char ';' + return$ tag l x + +subtree :: NameHack -> Parser TempTree +subtree name_hack = internal name_hack <|> leaf name_hack + +defaultMeta :: (Maybe Int, Double) +defaultMeta = (Nothing,0.0) + +leaf :: NameHack -> Parser TempTree +leaf name_hack = + do nm <- name + let nm' = name_hack nm + return$ NTLeaf (defaultMeta,Just nm') 0 + +internal :: NameHack -> Parser TempTree +internal name_hack = + do _ <- char '(' + bs <- branchset name_hack + _ <- char ')' + _nm <- name -- IGNORED, internal names. + return$ NTInterior (defaultMeta,Nothing) bs + +branchset :: NameHack -> Parser [TempTree] +branchset name_hack = + do b <- branch name_hack "at least one branch" + rest <- option [] $ try$ do char ','; branchset name_hack + return (b:rest) + +branch :: NameHack -> Parser TempTree +branch name_hack = + do s <- subtree name_hack + l <- branchMetadat name_hack + return$ tag l s + +-- If the length is omitted, it is implicitly zero. +branchMetadat :: NameHack -> Parser DefDecor +branchMetadat name_hack = option defaultMeta $ do + char ':' + n <- (try sciNotation <|> number) + -- IF the branch length succeeds then we go for the bracketed bootstrap value also: + bootstrap <- option Nothing $ do + char '[' + s <- many1 digit + char ']' + return (Just (read s)) + return (bootstrap,n) + +-- | Parse a normal, decimal number. +number :: Parser Double +number = + do sign <- option "" $ string "-" + first <- many1 digit + second <- option "0" $ try$ do char '.'; many1 digit + return (read (sign ++ first++"."++second) :: Double) + +-- | Parse a number in scientific notation. +sciNotation :: Parser Double +sciNotation = + do coeff <- do first <- many1 digit + second <- option "0" $ try$ do char '.'; many1 digit + return $ first++"."++second + char 'e' + sign <- option "" $ string "-" + expon <- many1 digit + return (read (coeff++"e"++sign++expon)) + +-- | Names are a mess... they contain all kinds of garbage sometimes it seems. +-- Thus we are very permissive. We allow anything that is not something we specifically need to reserve. +name :: Parser String +name = option "" $ many1 (noneOf "()[]:;, \t\n\r\f\v") +-- name = option "" $ many1 (letter <|> digit <|> oneOf "_.-") + + +-------------------------------------------------------------------------------- +-- Unit Tests +-------------------------------------------------------------------------------- + +tre1 :: NewickTree DefDecor +tre1 = snd $ parseNewick M.empty id "" $ B.pack "(A:0.1,B:0.2,(C:0.3,D:0.4):0.5);" + +unitTests :: Test +unitTests = test + [ "test name" ~: "foo" ~=? run name "foo" + , "test number" ~: 3.3 ~=? run number "3.3" + , "test number" ~: 3.0 ~=? run number "3" + , "test number" ~: -3.0 ~=? run number "-3" + + , "leaf" ~: ntl "A" ~=? run (leaf id) "A" + , "subtree" ~: ntl "A" ~=? run (subtree id) "A" + + -- These are not allowed: + , "null branchset" ~: errortest$ run (branchset id) "" + + , "internal" ~: NTInterior nada [ntl "A"] ~=? run (internal id) "(A);" + , "example: no nodes are named" ~: NTInterior nada + [ntl "", ntl "", NTInterior nada [ntl "", ntl ""]] + ~=? run (newick_parser id) "(,,(,));" + , "example: leaf nodes are named" ~: 6 ~=? treeSize (run (newick_parser id) "(A,B,(C,D));") + , "example: all nodes are named" ~: 6 ~=? treeSize (run (newick_parser id) "(A,B,(C,D)E)F;") + + , "example: all but root node have a distance to parent" ~: 6 ~=? treeSize (run (newick_parser id) "(:0.1,:0.2,(:0.3,:0.4):0.5);") + , "example: all have a distance to parent" ~: 6 ~=? treeSize (run (newick_parser id) "(:0.1,:0.2,(:0.3,:0.4):0.5):0.6;") + , "example: distances and leaf names (popular)" ~: 6 ~=? treeSize tre1 + , "example: distances and all names" ~: 6 ~=? treeSize (run (newick_parser id) "(A:0.1,B:0.2,(C:0.3,D:0.4)E:0.5)F;") + , "example: a tree rooted on a leaf node (rare)" ~: 6 ~=? treeSize (run (newick_parser id) "((B:0.2,(C:0.3,D:0.4)E:0.5)F:0.1)A;") + ] + where + ntl s = NTLeaf ((Nothing,0.0),Just s) 0 + nada = ((Nothing,0),Nothing) + +run :: Show a => Parser a -> String -> a +run p input = runB "" p (B.pack input) + +errortest :: t -> IO () +errortest x = + --() ~=? + handle (\ (e::SomeException) -> return ()) $ + do evaluate x + assertFailure "test was expected to throw an error" + diff --git a/Bio/Phylogeny/PhyBin/RFDistance.hs b/Bio/Phylogeny/PhyBin/RFDistance.hs new file mode 100644 --- /dev/null +++ b/Bio/Phylogeny/PhyBin/RFDistance.hs @@ -0,0 +1,198 @@ +{-# LANGUAGE ScopedTypeVariables, CPP #-} + +module Bio.Phylogeny.PhyBin.RFDistance + (DenseLabelSet, DistanceMatrix, + allBips, foldBips, + distanceMatrix, printDistMat) + where + +import Control.Monad +import Data.Word +import qualified Data.Vector as V +import qualified Data.Vector.Unboxed.Mutable as MV +import qualified Data.Vector.Unboxed as U +import qualified Data.Vector.Unboxed.Bit as UB +import qualified Data.Bit as B +import Text.PrettyPrint.HughesPJClass hiding (char, Style) +import System.IO (hPutStrLn, hPutStr, Handle) + +-- import Control.LVish +-- import qualified Data.LVar.Set as IS +-- import qualified Data.LVar.SLSet as SL + +-- import Data.LVar.Map as IM +-- import Data.LVar.NatArray as NA + +import Bio.Phylogeny.PhyBin.CoreTypes +-- import Data.BitList +import qualified Data.Set as S +import qualified Data.IntSet as SI +import qualified Data.Map.Strict as M +import qualified Data.Foldable as F +import Data.Monoid +import Prelude as P +import Debug.Trace + +-------------------------------------------------------------------------------- +-- A data structure choice +-------------------------------------------------------------------------------- + +-- | Dense sets of taxa, aka Bipartitions or BiPs +-- We assume that taxa labels have been mapped onto a dense, contiguous range of integers [0,N). +-- +-- NORMALIZATION Rule: Bipartitions are really two disjoint sets. But as long as +-- the parent set (the union of the partitions, aka "all taxa") then a bipartition +-- can be represented just by *one* subset. Yet we must choose WHICH subset for +-- consistency. We use the rule that we always choose the SMALLER. Thus the +-- DenseLabelSet should always be half the size or less, compared to the total +-- number of taxa. +-- +-- A set that is more than a majority of the taxa can be normalized by "flipping", +-- i.e. taking the taxa that are NOT in that set. +type DenseLabelSet = SI.IntSet +-- type DenseLabelSet s = BitList +-- type DenseLabelSet = UB.Vector B.Bit + +-- M.write vec lab (B.fromBool True) +-- mkEmptyDense size = U.replicate size (B.fromBool False) + +-- markLabel lab set = IS.putInSet lab set +-- mkEmptyDense _size = IS.newEmptySet + +markLabel :: Label -> DenseLabelSet -> DenseLabelSet +mkEmptyDense :: Int -> DenseLabelSet +denseUnions :: Int -> [DenseLabelSet] -> DenseLabelSet +bipSize :: DenseLabelSet -> Int + +markLabel lab set = SI.insert lab set +mkEmptyDense _size = SI.empty +denseUnions _size = SI.unions +bipSize = SI.size + + +-------------------------------------------------------------------------------- +-- Dirt-simple reference implementation +-------------------------------------------------------------------------------- + +type DistanceMatrix = V.Vector (U.Vector Int) + +-- | Returns a triangular distance matrix encoded as a vector. +distanceMatrix :: [NewickTree a] -> DistanceMatrix +distanceMatrix lst = + let sz = P.length lst + eachbips = V.fromList $ map allBips lst +-- in V.generate (sz-1) $ \ i -> + in V.generate sz $ \ i -> + U.generate i $ \ j -> + S.size (S.difference (eachbips V.! i) (eachbips V.! j)) + +-- | 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]) +labelBips tr = + trace ("labelbips "++show allLeaves++" "++show size) $ + loop tr + where + size = numLeaves tr + zero = mkEmptyDense size + loop (NTLeaf dec lab) = NTLeaf (dec, [markLabel lab zero]) lab + loop (NTInterior dec chlds) = + let chlds' = map loop chlds + sets = map (normBip . denseUnions size . snd . get_dec) chlds' in + NTInterior (dec, sets) chlds' + + halfSize = size `quot` 2 + normBip bip = + let flipped = SI.difference allLeaves bip in + case compare (SI.size bip) halfSize of + LT -> bip + GT -> flipped -- Flip it + EQ -> -- This is a painful case, we need a tie-breaker + min bip flipped + + allLeaves = leafSet tr + leafSet (NTLeaf _ lab) = SI.singleton lab + leafSet (NTInterior _ ls) = denseUnions size $ map leafSet ls + +foldBips :: Monoid m => (DenseLabelSet -> m) -> NewickTree a -> m +foldBips f tr = F.foldMap f' (labelBips tr) + where f' (_,bips) = F.foldMap f bips + +-- | Get all non-singleton BiPs implied by a tree. +allBips :: NewickTree a -> S.Set DenseLabelSet +allBips tr = S.filter ((> 1) . bipSize) $ foldBips S.insert tr S.empty + +-------------------------------------------------------------------------------- +-- Optimized, LVish version +-------------------------------------------------------------------------------- +-- First, necessary types: + +#if 0 +-- | A collection of all observed bipartitons (bips) with a mapping of which trees +-- contain which Bips. +type BipTable s = IMap DenseLabelSet s (SparseTreeSet s) +-- type BipTable = IMap BitList (U.Vector Bool) +-- type BipTable s = IMap BitList s (NA.NatArray s Word8) + +-- | Sets of taxa (BiPs) that are expected to be sparse. +type SparseTreeSet s = IS.ISet s TreeID +-- TODO: make this a set of numeric tree IDs... +-- NA.NatArray s Word8 + +type TreeID = AnnotatedTree +-- | Tree's are identified simply by their order within the list of input trees. +-- type TreeID = Int +#endif +-------------------------------------------------------------------------------- + +-- The distance matrix is an atomically-bumped matrix of numbers. +-- type DistanceMat s = NA.NatArray s Word32 +-- Except... bump isn't supported by our idempotent impl. + +#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 <- MV.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 + +-------------------------------------------------------------------------------- + +instance Pretty a => Pretty (S.Set a) where + pPrint s = pPrint (S.toList s) + + +printDistMat :: Handle -> V.Vector (U.Vector Int) -> IO () +printDistMat h mat = do + hPutStrLn h "Robinson-Foulds distance (matrix format):" + hPutStrLn h "-----------------------------------------" + V.forM_ mat $ \row -> do + U.forM_ row $ \elem -> do + hPutStr h (show elem) + hPutStr h " " + hPutStr h "0\n" + hPutStrLn h "-----------------------------------------" diff --git a/Bio/Phylogeny/PhyBin/Util.hs b/Bio/Phylogeny/PhyBin/Util.hs new file mode 100644 --- /dev/null +++ b/Bio/Phylogeny/PhyBin/Util.hs @@ -0,0 +1,111 @@ +{-# LANGUAGE ScopedTypeVariables #-} +-- RecordWildCards, TypeSynonymInstances, CPP +-- {-# LANGUAGE NamedFieldPuns #-} +-- {-# LANGUAGE OverloadedStrings #-} +-- {-# OPTIONS_GHC -fwarn-incomplete-patterns #-} +-- {-# OPTIONS_GHC -fwarn-unused-imports #-} + +-- | This module contains misc bits used by (multiple) other modules. + +module Bio.Phylogeny.PhyBin.Util + ( + is_regular_file, acquireTreeFiles + ) + 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.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.Process (system) +import System.Exit (ExitCode(..)) +import Test.HUnit ((~:),(~=?),Test,test) +import qualified HSH + +-- 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 + + +---------------------------------------------------------------------------------------------------- +-- OS specific bits: +---------------------------------------------------------------------------------------------------- +-- #ifdef WIN32 +-- is_regular_file = undefined +-- is_directory path = +-- getFileAttributes +-- --getFileInformationByHandle +-- -- bhfiFileAttributes +-- file_exists = undefined +-- #else +-- is_regular_file :: FilePath -> IO Bool +-- is_regular_file file = +-- do stat <- getFileStatus file; +-- -- Hmm, this is probably bad practice... hard to know its exhaustive: +-- return$ isRegularFile stat || isNamedPipe stat || isSymbolicLink stat +-- is_directory :: FilePath -> IO Bool +-- is_directory path = +-- do stat <- getFileStatus path +-- return (isDirectory stat) +-- file_exists = fileExist +-- #endif + +-- Here we ASSUME it exists, then these functions are good enough: +is_regular_file :: FilePath -> IO Bool +is_regular_file = doesFileExist + +is_directory :: FilePath -> IO Bool +is_directory = doesDirectoryExist + +file_exists :: FilePath -> IO Bool +file_exists path = + do f <- doesFileExist path + d <- doesDirectoryExist path + return (f || d) + + +-------------------------------------------------------------------------------- + +-- | Expand out directories to find all the tree files. +acquireTreeFiles :: [String] -> IO [String] +acquireTreeFiles inputs = do + all :: [[String]] <- forM inputs $ \ path -> do + exists <- file_exists path + + --stat <- if exists then getFileStatus path else return (error "internal invariant") + -- [2010.09.23] This is no longer really necessary: + if not exists then do + putStr$ "Input not a file/directory, assuming wildcard, using 'find' for expansion" + entries <- HSH.run$ "find " ++ path + putStrLn$ "("++show (length entries)++" files found): "++ show path + return entries + else do + isdir <- is_directory path + reg <- is_regular_file path + if isdir then do + putStr$ "Input is a directory, reading all regular files contained " + children <- getDirectoryContents path + filtered <- filterM is_regular_file $ map (combine path) children + putStrLn$ "("++show (length filtered)++" regular files found): "++ show path + return$ filtered + else if reg then do + return [path] + else error$ "phybin: Unhandled input path: " ++ path + + return (concat all) + diff --git a/Bio/Phylogeny/PhyBin/Visualize.hs b/Bio/Phylogeny/PhyBin/Visualize.hs new file mode 100644 --- /dev/null +++ b/Bio/Phylogeny/PhyBin/Visualize.hs @@ -0,0 +1,264 @@ +{-# LANGUAGE NamedFieldPuns #-} + +module Bio.Phylogeny.PhyBin.Visualize + (dotNewickTree, dotToPDF, viewNewickTree, + dotNewickTree_debug, + + -- * Dendrogram visualization + dendrogramToGraph, dotDendrogram + ) + where +import Text.Printf (printf) +import Data.List (elemIndex, isPrefixOf) +import Data.Maybe (fromJust) +import Data.Map ((!)) +import Data.Text.Lazy (pack) +import Control.Monad (void) +import Control.Concurrent (Chan, newChan, writeChan, forkIO) +import qualified Data.Graph.Inductive as G hiding (run) +import qualified Data.GraphViz as Gv hiding (parse, toLabel) +import qualified Data.GraphViz.Attributes.Complete as GA +import qualified Data.GraphViz.Attributes.Colors as GC +-- import Test.HUnit ((~:),(~=?),Test,test) + +import qualified Data.Clustering.Hierarchical as C + +import Bio.Phylogeny.PhyBin.CoreTypes + +---------------------------------------------------------------------------------------------------- +-- Visualization with GraphViz and FGL: +---------------------------------------------------------------------------------------------------- + +-- First we need to be able to convert our trees to FGL graphs: +toGraph :: FullTree StandardDecor -> G.Gr String Double +toGraph (FullTree _ tbl tree) = G.run_ G.empty $ loop tree + where + fromLabel ix = tbl ! ix + loop (NTLeaf _ name) = + do let str = fromLabel name + _ <- G.insMapNodeM str + return str + loop (NTInterior (StandardDecor{sortedLabels}) ls) = + do let bigname = concatMap fromLabel sortedLabels + names <- mapM loop ls + _ <- G.insMapNodeM bigname + mapM_ (\x -> G.insMapEdgeM (bigname, x, 0.0)) names + return bigname + +-- This version uses the tree nodes themselves as graph labels. +toGraph2 :: FullTree StandardDecor -> G.Gr (NewickTree StandardDecor) Double +toGraph2 (FullTree _ tbl tree) = G.run_ G.empty $ loop tree + where + loop node@(NTLeaf _ _) = + do _ <- G.insMapNodeM node + return () + loop node@(NTInterior _ ls) = + do mapM_ loop ls + _ <- G.insMapNodeM node + -- Edge weights as just branchLen (not bootstrap): + mapM_ (\x -> G.insMapEdgeM (node, x, branchLen$ get_dec x)) ls + return () + +-- dendrogramToGraph :: C.Dendrogram (FullTree a) -> G.Gr (Either String String) Double +dendrogramToGraph :: C.Dendrogram (FullTree a) -> G.Gr String Double +dendrogramToGraph x = G.run_ G.empty $ void$ loop x + where + -- deEither (Left s) = s + -- deEither (Right s) = s + loop node@(C.Leaf FullTree{treename}) = G.insMapNodeM (treename) + loop node@(C.Branch dist left right) = + do (_,l) <- loop left + (_,r) <- loop right + let ndname = "DUMMY_"++(l++"_"++r) + -- Right (deEither l ++"_"++ deEither r) + (midN,mid) <- G.insMapNodeM ndname + G.insMapEdgeM (l, mid, dist) + G.insMapEdgeM (r, mid, dist) + return (midN,mid) + + +-- | Open a GUI window to displaya tree. +-- +-- Fork a thread that then runs graphviz. +-- The channel retuned will carry a single message to signal +-- completion of the subprocess. +viewNewickTree :: String -> FullTree StandardDecor -> IO (Chan (), FullTree StandardDecor) +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 + writeChan chan () + --str <- prettyPrint d + --putStrLn$ "Generating the following graphviz tree:\n " ++ str + _ <- forkIO runit + --do runGraphvizCanvas Dot dot Xlib; return () + + return (chan, tree) + + +--default_cmd = TwoPi -- Totally ignores edge lengths. +default_cmd :: Gv.GraphvizCommand +default_cmd = Gv.Neato + +-- Show a float without scientific notation: +myShowFloat :: Double -> String +-- showFloat weight = showEFloat (Just 2) weight "" +myShowFloat fl = + let rnd = round fl in + if fl == fromIntegral rnd + 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 NewickTree to a graphviz Dot graph representation. +dotNewickTree :: String -> Double -> FullTree StandardDecor -> Gv.DotGraph G.Node +dotNewickTree title edge_scale atree@(FullTree _ tbl tree) = + --trace ("EDGE SCALE: " ++ show edge_scale) $ + Gv.graphToDot myparams graph + where + graph = toGraph2 atree + fromLabel ix = tbl ! ix + myparams :: Gv.GraphvizParams G.Node (NewickTree StandardDecor) Double () (NewickTree StandardDecor) + myparams = Gv.defaultParams { Gv.globalAttributes= [Gv.GraphAttrs [GA.Label$ GA.StrLabel$ pack title]], + Gv.fmtNode= nodeAttrs, Gv.fmtEdge= edgeAttrs } + nodeAttrs :: (Int,NewickTree StandardDecor) -> [GA.Attribute] + nodeAttrs (_num, node) = + let children = get_children node in + [ GA.Label$ GA.StrLabel$ pack$ + concatMap fromLabel $ + sortedLabels $ get_dec node + , GA.Shape (if null children then {-PlainText-} GA.Ellipse else GA.PointShape) + , GA.Style [GA.SItem GA.Filled []] + ] + edgeAttrs = getEdgeAttrs edge_scale + + +getEdgeAttrs :: Double -> (t, t1, Double) -> [GA.Attribute] +getEdgeAttrs edge_scale = edgeAttrs + where + -- TOGGLE: + -- edgeAttrs (_,_,weight) = [ArrowHead noArrow, Len (weight * edge_scale + bump), GA.Label (StrLabel$ show (weight))] + edgeAttrs (_,_,weight) = + let draw_weight = compute_draw_weight weight edge_scale in + --trace ("EDGE WEIGHT "++ show weight ++ " drawn at "++ show draw_weight) $ + [GA.ArrowHead Gv.noArrow, + GA.Label$ GA.StrLabel$ pack$ myShowFloat weight] ++ -- TEMPTOGGLE + --[ArrowHead noArrow, GA.Label (StrLabel$ show draw_weight)] ++ -- TEMPTOGGLE + if weight == 0.0 + then [GA.Color [weighted$ GA.X11Color Gv.Red], + GA.LWidth 3.0, GA.Len minlen] + else [GA.Len draw_weight] + + weighted c = GC.WC {GC.wColor=c, GC.weighting=Nothing} + minlen = 0.7 + maxlen = 3.0 + compute_draw_weight w scale = + let scaled = (abs w) * scale + minlen in + -- Don't draw them too big or it gets annoying: + (min scaled maxlen) + +-- | Some duplicated code with dotNewickTree. +dotDendrogram :: String -> Double -> C.Dendrogram (FullTree a) -> Gv.DotGraph G.Node +dotDendrogram title edge_scale dendro = + Gv.graphToDot myparams graph + where +-- graph :: G.Gr (Either String String) Double + graph :: G.Gr String Double + graph = dendrogramToGraph dendro + myparams :: Gv.GraphvizParams G.Node String Double () String -- (Either String String) + myparams = Gv.defaultParams { Gv.globalAttributes= [Gv.GraphAttrs [GA.Label$ GA.StrLabel$ pack title]], + Gv.fmtNode= nodeAttrs, + Gv.fmtEdge= edgeAttrs + } +-- nodeAttrs :: (Int, C.Dendrogram(FullTree StandardDecor)) -> [GA.Attribute] + nodeAttrs :: (Int, String) -> [GA.Attribute] + nodeAttrs (_num, eith) = + let (tag,shp) = -- case eith of + -- Left treename -> (take 60 treename, GA.Ellipse) + -- Right _ -> ("", GA.PointShape) + if isPrefixOf "DUMMY_" eith + then ("", GA.PointShape) + else (take 60 eith, GA.Ellipse) + in + [ GA.Label$ GA.StrLabel$ pack tag + , GA.Shape shp + , GA.Style [GA.SItem GA.Filled []] + ] + edgeAttrs = getEdgeAttrs edge_scale + + +-- | This version shows the ordered/rooted structure of the normalized tree. +dotNewickTree_debug :: String -> FullTree StandardDecor -> Gv.DotGraph G.Node +dotNewickTree_debug title atree@(FullTree _ tbl tree) = Gv.graphToDot myparams graph + where + graph = toGraph2 atree + fromLabel ix = tbl ! ix + myparams :: Gv.GraphvizParams G.Node (NewickTree StandardDecor) Double () (NewickTree StandardDecor) + myparams = Gv.defaultParams { Gv.globalAttributes= [Gv.GraphAttrs [GA.Label$ GA.StrLabel$ pack title]], + Gv.fmtNode= nodeAttrs, Gv.fmtEdge= edgeAttrs } + nodeAttrs :: (Int,(NewickTree StandardDecor)) -> [GA.Attribute] + nodeAttrs (num,node) = + let children = get_children node in + [ GA.Label (if null children + then GA.StrLabel$ pack$ concatMap fromLabel $ sortedLabels $ get_dec node + else GA.RecordLabel$ take (length children) $ + -- This will leave interior nodes unlabeled: + map (GA.PortName . GA.PN . pack) $ map show [1..] + -- This version gives some kind of name to interior nodes: +-- map (\ (i,ls) -> LabelledTarget (PN$ show i) (fromLabel$ head ls)) $ +-- zip [1..] (map (thd3 . get_dec) children) + ) + , GA.Shape GA.Record + , GA.Style [GA.SItem GA.Filled []] + ] + + edgeAttrs (num1, num2, _weight) = + let node1 = fromJust$ G.lab graph num1 + node2 = fromJust$ G.lab graph num2 + ind = fromJust$ elemIndex node2 (get_children node1) + in [GA.TailPort$ GA.LabelledPort (GA.PN$ pack$ show$ 1+ind) (Just GA.South)] + + + + +-------------------------------------------------------------------------------- +-- Unit Testing +-------------------------------------------------------------------------------- + + +-- tt0 :: IO (Chan (), FullTree) +-- tt0 = drawNewickTree "tt0" $ annotateWLabLists $ parseNewick "" "(A,(C,D,E),B);" + +-- tt3 :: IO (Chan (), FullTree) +-- tt3 = drawNewickTree "tt3" tt + +-- tt4 :: IO (Chan (), FullTree) +-- tt4 = drawNewickTree "tt4"$ trace ("FINAL: "++ show (pPrint norm4)) $ norm4 + +-- tt5 :: IO (Chan (), FullTree) +-- tt5 = drawNewickTree "tt5"$ norm5 + + +-- tt5' :: String +-- tt5' = prettyPrint' $ dotNewickTree "norm5" 1.0 norm5 + +-- ttall :: IO (Chan (), FullTree) +-- ttall = do tt3; tt4; tt5 + +-- tt2 :: G.Gr String Double +-- tt2 = toGraph tt + + +-- unitTests :: Test +-- unitTests = test +-- [ +-- ] + + +-- TEMP / HACK: +prettyPrint' :: Show a => a -> String +prettyPrint' = show diff --git a/Main.hs b/Main.hs --- a/Main.hs +++ b/Main.hs @@ -1,26 +1,35 @@ -{-# LANGUAGE RecordWildCards #-} -{-# OPTIONS_GHC -fwarn-unused-imports #-} +{-# LANGUAGE RecordWildCards, TupleSections #-} +{-# OPTIONS_GHC -fwarn-unused-imports -fwarn-incomplete-patterns #-} module Main where import Data.List (sort) import qualified Data.ByteString.Lazy.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 System.Environment (getArgs, withArgs) import System.Console.GetOpt (OptDescr(Option), ArgDescr(..), ArgOrder(..), usageInfo, getOpt) import System.Exit (exitSuccess) -import Test.HUnit (runTestTT, Test, test) +import System.IO (stdout) +import Test.HUnit (runTestTT, Test, test, (~:)) +import Control.Applicative ((<$>)) import Data.GraphViz (runGraphvizCanvas,GraphvizCommand(Dot),GraphvizCanvas(Xlib)) -import Bio.Phylogeny.PhyBin.CoreTypes - ( NewickTree(..), PhyBinConfig(..), default_phybin_config, DefDecor, StandardDecor(..), - toLabel, fromLabel, Label, set_dec, map_labels ) -import Bio.Phylogeny.PhyBin (driver, binthem, normalize, annotateWLabLists, unitTests) -import Bio.Phylogeny.PhyBin.Parser (parseNewick, unitTests) +import Text.PrettyPrint.HughesPJClass hiding (char, Style) + +import Bio.Phylogeny.PhyBin.CoreTypes +import Bio.Phylogeny.PhyBin (driver, binthem, normalize, annotateWLabLists, + unitTests, acquireTreeFiles, deAnnotate) +import Bio.Phylogeny.PhyBin.Parser (parseNewick, parseNewicks, parseNewickFiles, unitTests) import Bio.Phylogeny.PhyBin.Visualize (viewNewickTree, dotNewickTree_debug) +import Bio.Phylogeny.PhyBin.RFDistance (distanceMatrix, printDistMat, allBips) +import Bio.Phylogeny.PhyBin.PreProcessor +import qualified Data.Clustering.Hierarchical as C + + import Version ---------------------------------------------------------------------------------------------------- @@ -41,6 +50,11 @@ | TabDelimited Int Int | SelfTest + | RFMatrix | LineSetDiffMode | PrintNorms | PrintReg + | Cluster C.Linkage + | BinningMode + | EditDistThresh Int + | DendogramOnly | NameCutoff String | NamePrefix Int @@ -48,6 +62,7 @@ deriving (Show, Eq, Ord) -- ORD is really used. + parseTabDelim :: String -> Flag parseTabDelim _str = TabDelimited 8 9 @@ -57,43 +72,55 @@ [ Option ['v'] ["verbose"] (NoArg Verbose) "print WARNINGS and other information (recommended at first)" , Option ['V'] ["version"] (NoArg Version) "show version number" - , Option ['o'] ["output"] (ReqArg Output "DIR") "set directory to contain all output files (default \"./\")" - + , Option ['o'] ["output"] (ReqArg Output "DIR") "set directory to contain all output files (default \"./phybin_out/\")" , Option [] ["selftest"] (NoArg SelfTest) "run internal unit tests" - - + {- -- TODO: FIXME: IMPLEMENT THIS: , Option [] [] (NoArg NullOpt) "" , Option ['t'] ["tabbed"] (ReqArg parseTabDelim "NUM1:NUM2")$ "assume the input is a tab-delimited file with gene names \n"++ "in column NUM1 and Newick trees in NUM2" --} +-} + , Option [] [] (NoArg NullOpt) "" + , Option [] [] (NoArg$ error "internal problem") "----------------------------- Clustering Options ------------------------------" + , Option [] ["bin"] (NoArg BinningMode)$ "Use simple binning, the cheapest form of 'clustering'" + , Option [] ["single"] (NoArg$ Cluster C.SingleLinkage) $ "Use single-linkage clustering (nearest neighbor)" + , Option [] ["complete"] (NoArg$ Cluster C.CompleteLinkage)$ "Use complete-linkage clustering (furthest neighbor)" + , Option [] ["UPGMA"] (NoArg$ Cluster C.UPGMA) $ "Use Unweighted Pair Group Method (average linkage)" + + , Option [] ["editdist"] (ReqArg (EditDistThresh . read) "DIST")$ + "Combine all clusters separated by DIST or less. Report a flat list of clusters." + , Option [] ["dendogram"] (NoArg DendogramOnly)$ "Report a hierarchical clustering (default)" + , Option [] [] (NoArg NullOpt) "" , Option [] [] (NoArg$ error "internal problem") "----------------------------- Visualization --------------------------------" - , Option ['g'] ["graphbins"] (NoArg Graph) "use graphviz to produce .dot and .pdf output files named bin1.*, bin2.*, etc" + , Option ['g'] ["graphbins"] (NoArg Graph) "use graphviz to produce .dot and .pdf output files named cluster1.*, cluster2.*, etc" +-- TODO: Produce the consensus tree as well as the individual trees. , Option ['d'] ["drawbins"] (NoArg Draw) "like -g, but open GUI windows to show a tree for each bin" , Option ['w'] ["view"] (NoArg View)$ "for convenience, \"view mode\" simply displays input Newick files without binning" , Option [] [] (NoArg NullOpt) "" - , Option [] [] (NoArg$ error "internal problem") "--------------------------- Handling taxa names ----------------------------" --- , Option ['n'] ["numtaxa"] (ReqArg (NumTaxa . 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? - + , Option [] [] (NoArg$ error "internal problem") "---------------------------- Tree pre-processing -----------------------------" + , Option ['n'] ["numtaxa"] (ReqArg (NumTaxa . read) "NUM") "expect NUM taxa for this dataset" , Option ['b'] ["branchcut"] (ReqArg (BranchThresh . read) "LEN") "collapse branches less than LEN" + + , 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)" + -- 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? - {- -- TODO: FIXME: IMPLEMENT THIS: , Option ['f'] ["force"] (NoArg Force) "force phybin to consume and bin trees with different numbers of taxa" -} , Option [] [] (NoArg NullOpt) "" , Option ['p'] ["nameprefix"] (ReqArg (NamePrefix . read) "NUM") $ - "Leaf names in the input Newick trees are usually gene names, not taxa.\n"++ - "It is typical to extract taxa names from genes. This option extracts a\n"++ - "prefix of NUM characters to serve as the taxa name." + "Leaf names in the input Newick trees can be gene names, not taxa.\n"++ + "Then it is typical to extract taxa names from genes. This option extracts\n"++ + "a prefix of NUM characters to serve as the taxa name." , Option [] [] (NoArg NullOpt) "" , Option ['s'] ["namesep"] (ReqArg NameCutoff "STR") $ --"\n"++ @@ -107,31 +134,94 @@ "to compute taxa names, e.g. if multiple genes/plasmids map onto one taxa.\n"++ "This option specifies a text file with find/replace entries of the form\n"++ "\" \", which are applied AFTER -s and -p." + + , Option [] [] (NoArg NullOpt) "" + , Option [] [] (NoArg$ error "internal problem") "--------------------------- Utility Modes ----------------------------" + , Option [] ["rfdist"] (NoArg RFMatrix) "print a Robinson Foulds distance matrix for the input trees" + , Option [] ["setdiff"] (NoArg LineSetDiffMode) "for convenience, print the set difference between cluster*.txt files" + , Option [] ["print"] (NoArg PrintReg) "simply print out a concise form of each input tree" + , Option [] ["printnorms"] (NoArg PrintNorms) "simply print out a concise and NORMALIZED form of each input tree" ] usage :: String usage = "\nUsage: phybin [OPTION...] files or directories...\n\n"++ - "PhyBin takes Newick tree files as input and produces, at minimum, files of the form\n"++ - "binXX_YY.tr, each containing a list of input file paths that fall into that bin.\n\n"++ +-- "MODE must be one of 'bin' or 'cluster'.\n\n" ++ - "USAGE NOTES: Currently phybin ignores input trees with the wrong number of taxa.\n"++ - "If given a directory as input phybin will assume all contained files are Newick trees.\n\n"++ + "PhyBin takes Newick tree files as input. Paths of Newick files can\n"++ + "be passed directly on the command line. Or, if directories are provided,\n"++ + "all files in those directories will be read. Taxa are named based on the files\n"++ + "containing them. If a file contains multiple trees, all are read by phybin, and\n"++ + "the taxa name then includes a suffix indicating the position in the file:\n"++ + " e.g. FILENAME_0, FILENAME_1, etc.\n"++ + "\n"++ + -- "In binning mode, Phybin output contains, at minimum, files of the form binXX_YY.tr,\n"++ + -- "each containing a list of input file paths that fall into that bin. XX signifies\n"++ + -- "the rank of the bin (biggest first), and YY is how many trees it contains.\n"++ + -- "\n"++ + + "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"++ + "\n"++ + + "Binning mode provides an especially quick-and-dirty form of clustering.\n"++ + "When running with the --bin option, only exactly equal trees are put in the same cluster.\n"++ + "Tree pre-processing still applies, however: for example collapsing short branches.\n"++ + "\n"++ + + "USAGE NOTES: \n"++ + " * Currently phybin ignores input trees with the wrong number of taxa.\n"++ + " * If given a directory as input phybin will assume all contained files are Newick trees.\n\n"++ + "\nOptions include:\n" defaultErr :: [String] -> t defaultErr errs = error $ "ERROR!\n" ++ (concat errs ++ usageInfo usage options) +-------------------------------------------------------------------------------- +-- Aggregated Unit Tests +-------------------------------------------------------------------------------- allUnitTests :: Test -- allUnitTests = unitTests ++ allUnitTests = test [ Bio.Phylogeny.PhyBin.unitTests , Bio.Phylogeny.PhyBin.Parser.unitTests + , "norm/Bip1" ~: (testNorm prob1) ] -- Bio.Phylogeny.PhyBin.Parser.unitTests +-- "(5_, (19, ((3_, 14), ((2_, 1_), (7_, (6_, 18))))), 13)" +-- [2013.07.23] +-- This was INCORRECTLY normalizing to: +-- ((1_, 2_), (7_, (18, 6_)), ((14, 3_), (19, (13, 5_)))) +prob1 = "(5_, (19, ((3_, 14), ((2_, 1_), (7_, (6_, 18))))), 13);" + +-- | Make sure that the normalized version of a tree yields the same bipartitions as +-- the unnormalized one. +testNorm :: String -> IO () +testNorm str = do + let (labs,parsed) = parseNewick M.empty id "test" (B.pack str) + normed = normalize $ annotateWLabLists parsed + bips1 = allBips parsed + bips2 = allBips normed + added = S.difference bips2 bips1 + removed = S.difference bips1 bips2 + dispBips bip = show$ + map (map (labs M.!)) $ + map IS.toList$ S.toList bip + unless (bips1 == bips2) $ do + putStrLn$ "Normalized this: "++show (displayDefaultTree $ FullTree "" labs parsed) + putStrLn$ "To this : "++show (displayDefaultTree $ deAnnotate $ FullTree "" labs normed) + error$ "Normalization added and removed these bipartitions, respectively:\n " + ++dispBips added ++"\n " + ++dispBips removed + + +-------------------------------------------------------------------------------- + main :: IO () main = do argv <- getArgs @@ -141,13 +231,34 @@ (o,n,[] ) -> return (o,n) (_,_,errs) -> defaultErr errs + all_inputs <- acquireTreeFiles files let process_opt cfg opt = case opt of NullOpt -> return cfg Verbose -> return cfg { verbose= True } Version -> do putStrLn$ "phybin version "++phybin_version; exitSuccess SelfTest -> do _ <- runTestTT allUnitTests; exitSuccess - + + RFMatrix -> return cfg { print_rfmatrix= True } + + LineSetDiffMode -> do + bss <- mapM B.readFile files + case map (S.fromList . B.lines) bss of + [set1,set2] -> do let [f1,f2] = files + let diff = S.difference set1 set2 + putStrLn$" Found "++show(S.size diff)++" lines occuring in "++f1++" but not "++f2 + mapM_ B.putStrLn $ S.toList diff + oth -> error $"Line set difference mode expects two files as input, got "++show(length oth) + exitSuccess + + PrintNorms -> return cfg + PrintReg -> return cfg + + Cluster lnk -> return cfg { clust_mode = ClusterThem lnk } + 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 } @@ -160,27 +271,43 @@ Force -> error "force option not handled yet" - NameCutoff str -> let set = S.fromList str - new = toLabel . takeWhile (not . flip S.member set) . fromLabel - in return cfg { name_hack = new . name_hack cfg } - NamePrefix n -> let new = toLabel . (take n) . fromLabel - in return cfg { name_hack = new . name_hack cfg } + NameCutoff str -> let chopper = takeWhile (not . flip S.member (S.fromList str)) + in return cfg { name_hack = chopper . name_hack cfg } + NamePrefix n -> return cfg { name_hack = (take n) . name_hack cfg } -- This should always be after cutoff/prefix: NameTable file -> do reader <- name_table_reader file return cfg { name_hack = reader . name_hack cfg } - - config <- foldM process_opt default_phybin_config{ inputs=files } + config <- foldM process_opt default_phybin_config{ inputs= all_inputs } (sort opts) -- NOTE: options processed in sorted order. when (null files) $ do defaultErr ["No file arguments!\n"] - if View `elem` opts - then view_graphs config - --else driver config{ name_hack= name_hack_legionella } - else driver config + ------------------------------------------------------------ + -- This mode kicks in AFTER config options are processed. + when (elem PrintReg opts) $ do + (_,fts) <- parseNewickFiles (name_hack config) all_inputs + forM_ fts $ \ ft@(FullTree name _ _) -> do + putStrLn $ "Tree "++show name + putStrLn$ show$ displayDefaultTree ft + exitSuccess + ------------------------------------------------------------ + when (elem PrintNorms opts) $ do + (_,fts) <- parseNewickFiles (name_hack config) all_inputs + forM_ fts $ \ ft@(FullTree name _ _) -> do + putStrLn $ "Tree , NORMALIZED:"++show name + putStrLn$ show$ displayDefaultTree$ deAnnotate $ + liftFT (normalize . annotateWLabLists) ft + exitSuccess + ------------------------------------------------------------ + when (View `elem` opts) $ do + view_graphs config + exitSuccess + ------------------------------------------------------------ + -- Otherwise do the main, normal thing: + driver config view_graphs :: PhyBinConfig -> IO () view_graphs PBC{..} = @@ -188,10 +315,8 @@ putStrLn$ "Drawing "++ file ++"...\n" str <- B.readFile file putStrLn$ "Parsed: " ++ (B.unpack str) - (chan, _tr) <- viewNewickTree file $ - annotateWLabLists$ - map_labels name_hack $ - parseNewick file str + let (tbl,tr) = parseNewick M.empty name_hack file str + (chan, _tr) <- viewNewickTree file (FullTree "" tbl (annotateWLabLists tr)) return chan forM_ chans readChan return () @@ -212,7 +337,8 @@ lines contents return$ \ name_to_hack -> - case M.lookup name_to_hack mp of -- Could use a trie + + case M.lookup name_to_hack mp of -- Could use a trie Just x -> x Nothing -> name_to_hack where @@ -227,6 +353,31 @@ temp :: IO () temp = driver default_phybin_config{ num_taxa=7, inputs=["../datasets/test.tr"] } + +-- 112 and 13 +rftest = do + (mp,[t1,t2]) <- parseNewickFiles (take 2) ["tests/13.tr", "tests/112.tr"] + putStrLn$ "Tree 13 : " ++ show (displayDefaultTree t1) + putStrLn$ "Tree 112 : "++ show (displayDefaultTree t2) + + putStrLn$ "Tree 13 normed : "++ show (disp t1) + putStrLn$ "Tree 112 normed : "++ show (disp t2) + + putStrLn$ "13 collapsed 0.02: " ++show (disp$ liftFT (collapseBranchLenThresh 0.02) t1) + putStrLn$ "112 collapsed 0.02: " ++show (disp$ liftFT (collapseBranchLenThresh 0.02) t2) + + 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] + printDistMat stdout mat + return () + where + disp (FullTree nm labs tr) = + let collapsed :: AnnotatedTree + collapsed = normalize$ annotateWLabLists tr + in displayDefaultTree$ deAnnotate $ FullTree nm labs collapsed + ---------------------------------------------------------------------------------------------------- -- TODO: expose a command line argument for testing. -- The below test exposed my normalization bug relating to "deriving Ord". @@ -237,6 +388,8 @@ withArgs ["-w","~/newton_and_newton_local/datasets/yersinia/yersinia_trees/111.dnd","-m","../datasets/yersinia/name_table_hack_yersinia.txt"] main +-- [2013.07.22] Disabling for the new Label representation: +{- pa :: NewickTree DefDecor pa = set_dec (Nothing,1) $ NTInterior () [NTInterior () [NTLeaf () "RE",NTInterior () [NTLeaf () "SD",NTLeaf () "SM"]],NTInterior () [NTLeaf () "BB",NTLeaf () "BJ"],NTInterior () [NTLeaf () "MB",NTLeaf () "ML"]] @@ -270,8 +423,6 @@ a_norm :: NewickTree StandardDecor a_norm = normalize (annotateWLabLists$ snd a_) ---a_norm = NTInterior (0.13899,7,["BB","BJ","MB","ML","RE","SD","SM"]) [NTInterior (0.0,3,["RE","SD","SM"]) [NTLeaf (0.13143,1,["RE"]) "RE",NTInterior (5.697e-2,2,["SD","SM"]) [NTLeaf (5.977e-2,1,["SD"]) "SD",NTLeaf (3.95e-2,1,["SM"]) "SM"]],NTInterior (9.019e-2,2,["BB","BJ"]) [NTLeaf (0.11856,1,["BB"]) "BB",NTLeaf (0.13592,1,["BJ"]) "BJ"],NTInterior (0.13194,2,["MB","ML"]) [NTLeaf (0.19456,1,["MB"]) "MB",NTLeaf (0.16603,1,["ML"]) "ML"]] - b_norm_ :: NewickTree (Double, Int, [Label]) b_norm_ = NTInterior (0.0,7,["BB","BJ","MB","ML","RE","SD","SM"]) [NTInterior (0.23143,2,["BB","BJ"]) [NTLeaf (9.192e-2,1,["BB"]) "BB",NTLeaf (0.10125,1,["BJ"]) "BJ"],NTInterior (6.621e-2,2,["MB","ML"]) [NTLeaf (0.16184,1,["MB"]) "MB",NTLeaf (0.15233,1,["ML"]) "ML"],NTInterior (6.527e-2,3,["RE","SD","SM"]) [NTLeaf (0.18443,1,["RE"]) "RE",NTInterior (0.13734,2,["SD","SM"]) [NTLeaf (3.002e-2,1,["SD"]) "SD",NTLeaf (2.975e-2,1,["SM"]) "SM"]]] @@ -298,3 +449,4 @@ withBootstrap2 :: String withBootstrap2 = "((((A8F330_:0.01131438136322714984,(G0GWK2_:0.00568050636963043226,(Q92FV4_:0.00284163304504484121,((B0BVQ5_:0.00319487112504297311,A8GU65_:0.00000122123005994819):0.00279881991324161267[74],(C3PM27_:0.00560787769333294297,C4K2Z0_:0.00559642713265556899):0.00000122123005994819[15]):0.00000122123005994819[4]):0.00276851661606284868[56]):0.00283144414216590342[60]):0.00886304965525876697[76],(A8GQC0_:0.05449879836105625541,(A8F0B2_:0.04736199885985507840,Q4UJN9_:0.02648399728559588939):0.00905997055810744446[64]):0.00323255855543533657[28]):0.02237505187863457132[29],(Q1RGK5_:0.00000122123005994819,A8GYD7_:0.00000122123005994819):0.28299884298270094884[100]):0.05776841634437222123[100],(Q9ZC84_:0.00000122123005994819,D5AYH5_:0.00000122123005994819):0.00951976341375833368[99],Q68VM9_:0.04408933524904214141);" +-} diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,71 @@ +% PhyBin: Binning Trees by Topology + + +PhyBin is a simple command line tool that classifies (bins) 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. + +![(Above figure) Trees corresponding to the three largest bins resulting from a + 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) + +Invoking PhyBin +=============== + +PhyBin is a command-line program that produces output in the form of +text files and pdfs, but to produce pdfs (to visualize trees) the + [GraphViz program](http://www.graphviz.org/), +including the `dot` command, must be installed on the machine. + +The following is a typical invocation of PhyBin: + + phybin *.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 +requested, it will also produce visual representations of each bin in +the form `output_dir/binXX_YY.pdf`. + +Downloading and Installing PhyBin +================================= + +The source code to PhyBin can be downloaded here: + + * [Download Source Tarball](phybin-0.1.2.tar.gz) + +PhyBin is written in Haskell and if you have + [Haskell Platform](http://hackage.haskell.org/platform/). +installed you can install phybin with this one-liner: + + cabal install phybin + +PhyBin is also available for download as a statically-linked +executable for Mac-OS and Linux: + + * [Download Mac-OS Binary](phybin-0.1.2.mac) + * [Download Linux Binary](phybin-0.1.2.x86_64_linux) + +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. + +- - - - - - - - - - - - - - - +Authors: Irene and Ryan Newton + +Contact email: `irnewton` `indiana` `edu` (with "at" and "dot" inserted). + +. + + diff --git a/phybin.cabal b/phybin.cabal --- a/phybin.cabal +++ b/phybin.cabal @@ -1,5 +1,5 @@ Name: phybin -Version: 0.1.2.5 +Version: 0.2.2 License: BSD3 License-file: LICENSE Stability: Beta @@ -12,53 +12,71 @@ -- 0.1.2.1 -- -- 0.1.2.4 -- new release, several new features -- 0.1.2.5 -- bump for graphviz API changes +-- 0.2 -- Add Robinson-Foulds distance, use Int labels. +-- 0.2.2 -- misc changes and expose library -- homepage: http://code.haskell.org/phybin homepage: http://people.csail.mit.edu/newton/phybin/ Copyright: Copyright (c) 2010 Ryan Newton -Synopsis: Utility for binning phylogenetic trees in Newick format. +Synopsis: Utility for clustering phylogenetic trees in Newick format based on Robinson-Foulds distance. Description: - Classifies (bins) input Newick trees by topology, creating output files that - characterize the size and contents of each bin (including + This package provides a libary and executable for dealing with Newick tree files. + . + It can do simple binning of identical trees or more complex clustering based on + an all-to-all Robinson-Foulds distance matrix. + . + Output files that characterize the size and contents of each bin or cluster (including generating GraphViz-based visual representations of the tree topologies). Category: Bioinformatics -Cabal-Version: >=1.6 +Cabal-Version: >=1.8 -build-type: Simple +Extra-source-files: README.md -source-repository head +Build-type: Simple + +Source-repository head type: git location: git://github.com/rrnewton/PhyBin.git --- location: https://github.com/rrnewton/PhyBin --- - -- type: darcs - -- location: http://code.haskell.org/phybin/repo/ +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.Util + Bio.Phylogeny.PhyBin.Visualize +-- Other-modules: Data.BitList + Build-Depends: base >= 3 && < 5, directory, process, containers, unix, + async, time, + filepath, + graphviz >= 2999.16, + text >= 0.11 && < 0.12, + prettyclass, fgl, + HSH, HUnit, bytestring, + -- For bytestring.lazy support: + parsec >= 3.1.0, + bitvec >= 0.1, vector >= 0.10, + hierarchical-clustering >= 0.4 +-- lattice-par, + GHC-Options: -O2 -funbox-strict-fields -rtsopts + Executable phybin Main-is: Main.hs --- Main-is: Bio/Phylogeny/PhyBin/Main.hs + Build-Depends: phybin + -- DUPLICATE: Build-Depends: base >= 3 && < 5, directory, process, containers, unix, - stringtable-atom, filepath, - -- graphviz == 2999.11.0.0, - -- graphviz >= 2999.12, - -- More API Changes: [2013.07.19] - graphviz >= 2999.16, + async, time, + filepath, + graphviz >= 2999.16, text >= 0.11 && < 0.12, prettyclass, fgl, HSH, HUnit, bytestring, -- For bytestring.lazy support: - parsec >= 3.1.0 --- extensions: CPP --- extensions: ScopedTypeVariables, RecordWildCards, TypeSynonymInstances - GHC-Options: -O2 --- GHC-Options: -O2 -threaded - --- Doesn't work for an executable? --- extra-source-files: README.txt - --- Oh, well this doesn't seem right but it achieves the effect for now: --- c-sources: README.txt --- includes: README.txt - + parsec >= 3.1.0, + bitvec >= 0.1, vector >= 0.10, + hierarchical-clustering >= 0.4 +-- lattice-par, + GHC-Options: -O2 -funbox-strict-fields -rtsopts