phybin 0.1.2.5 → 0.2.2
raw patch · 10 files changed
+2298/−76 lines, 10 filesdep +asyncdep +bitvecdep +hierarchical-clusteringdep −stringtable-atom
Dependencies added: async, bitvec, hierarchical-clustering, phybin, time, vector
Dependencies removed: stringtable-atom
Files
- Bio/Phylogeny/PhyBin.hs +438/−0
- Bio/Phylogeny/PhyBin/Binning.hs +432/−0
- Bio/Phylogeny/PhyBin/CoreTypes.hs +321/−0
- Bio/Phylogeny/PhyBin/Parser.hs +217/−0
- Bio/Phylogeny/PhyBin/RFDistance.hs +198/−0
- Bio/Phylogeny/PhyBin/Util.hs +111/−0
- Bio/Phylogeny/PhyBin/Visualize.hs +264/−0
- Main.hs +198/−46
- README.md +71/−0
- phybin.cabal +48/−30
+ Bio/Phylogeny/PhyBin.hs view
@@ -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<N>_<size>.txt"+ putStrLn$ " Wrote representative trees to cluster<N>_<size>.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++"<N>_<size>.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++"<N>_<binsize>.txt"+ putStrLn$ " Wrote representative trees to "++filePrefix++"<N>_<binsize>.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<N>_<binsize>.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)
+ Bio/Phylogeny/PhyBin/Binning.hs view
@@ -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 + ]+++--------------------------------------------------------------------------------
+ Bio/Phylogeny/PhyBin/CoreTypes.hs view
@@ -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"+
+ Bio/Phylogeny/PhyBin/Parser.hs view
@@ -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 "<unknown>" p (B.pack input)++errortest :: t -> IO ()+errortest x = + --() ~=?+ handle (\ (e::SomeException) -> return ()) $ + do evaluate x+ assertFailure "test was expected to throw an error"+
+ Bio/Phylogeny/PhyBin/RFDistance.hs view
@@ -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 "-----------------------------------------"
+ Bio/Phylogeny/PhyBin/Util.hs view
@@ -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)+
+ Bio/Phylogeny/PhyBin/Visualize.hs view
@@ -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
Main.hs view
@@ -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"++ "\"<string> <taxaname>\", 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);"+-}
+ README.md view
@@ -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.++++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).++.++
phybin.cabal view
@@ -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