phybin 0.2.2 → 0.2.11
raw patch · 9 files changed
+994/−320 lines, 9 filesdep +splitdep −HSHdep −unix
Dependencies added: split
Dependencies removed: HSH, unix
Files
- Bio/Phylogeny/PhyBin.hs +188/−128
- Bio/Phylogeny/PhyBin/Binning.hs +8/−6
- Bio/Phylogeny/PhyBin/CoreTypes.hs +29/−6
- Bio/Phylogeny/PhyBin/Parser.hs +36/−15
- Bio/Phylogeny/PhyBin/RFDistance.hs +274/−37
- Bio/Phylogeny/PhyBin/Util.hs +10/−10
- Bio/Phylogeny/PhyBin/Visualize.hs +211/−49
- Main.hs +191/−58
- phybin.cabal +47/−11
Bio/Phylogeny/PhyBin.hs view
@@ -9,12 +9,12 @@ module Bio.Phylogeny.PhyBin ( driver, binthem, normalize, annotateWLabLists, unitTests, acquireTreeFiles,- deAnnotate )+ deAnnotate, retrieveHighlights, matchAnyHighlight ) where import qualified Data.Foldable as F import Data.Function (on)-import Data.List (delete, minimumBy, sortBy, foldl1')+import Data.List (delete, minimumBy, sortBy, foldl1', foldl', intersperse, isPrefixOf) import Data.Maybe (fromMaybe) import Data.Either (partitionEithers) import Data.Time.Clock@@ -24,14 +24,15 @@ 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 qualified Control.Concurrent.Async as Async import Control.Exception (evaluate) import Control.Applicative ((<$>),(<*>)) import Control.Concurrent (Chan) import System.FilePath (combine)-import System.Directory (doesFileExist, doesDirectoryExist,+import System.Directory (doesFileExist, doesDirectoryExist, createDirectoryIfMissing, getDirectoryContents, getCurrentDirectory) import System.IO (openFile, hClose, IOMode(..), stdout)+import System.Info (os) import System.Process (system) import System.Exit (ExitCode(..)) import Test.HUnit ((~:),(~=?),Test,test)@@ -48,18 +49,34 @@ import Bio.Phylogeny.PhyBin.Util import Debug.Trace+---------------------------------------------------------------------------------------------------- -- Turn on for extra invariant checking: debug :: Bool debug = True +#ifdef SEQUENTIALIZE+#warning "SEQUENTIALIING execution. Disabling all parallelism."+type Async t = t+async a = a+wait x = return x +#else+type Async t = Async.Async t+async = Async.async+wait = Async.wait+#endif++-- | A dendrogram PLUS consensus trees at the intermediate points.+data DendroPlus a = DPLeaf (FullTree a)+ | DPBranch !Double (FullTree ()) (DendroPlus a) (DendroPlus a)+ ---------------------------------------------------------------------------------------------------- -- | 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 } =+driver cfg@PBC{ verbose, num_taxa, name_hack, output_dir, inputs=files,+ do_graph, branch_collapse_thresh, highlights, + dist_thresh, clust_mode, use_hashrf, print_rfmatrix } = -- Unused: do_draw do --------------------------------------------------------------------------------@@ -69,17 +86,18 @@ --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+ unless bl $ createDirectoryIfMissing True output_dir - 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"+ if isPrefixOf "mingw" os then+ -- TODO: Need a portable version of this. 'filemanip' would do:+ putStrLn$ "Not cleaning away previous phybin outputs (TODO: port this to Windows)."+ else do + 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"+ return () putStrLn$ "Parsing "++show (length files)++" Newick tree files." --putStrLn$ "\nFirst ten \n"++ concat (map (++"\n") $ map show $ take 10 files)@@ -97,7 +115,12 @@ let num_files = length goodFiles bstrs <- mapM B.readFile goodFiles let (labelTab, fulltrees) = parseNewicks name_hack (zip files bstrs)-+ + highlightTrs <- retrieveHighlights name_hack labelTab highlights+ + putStrLn$ "\nTotal unique taxa ("++ show (M.size labelTab) ++"):\n "+++ (unwords$ M.elems labelTab)+-- show (nest 2 $ sep $ map text $ M.elems labelTab) -------------------------------------------------------------------------------- case branch_collapse_thresh of @@ -138,10 +161,11 @@ 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)+ let all_branches = concatMap (F.foldr' (\x ls -> getBranchLen x:ls) []) validtrees + unless (null all_branches) $ do+ 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.@@ -150,54 +174,48 @@ (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)) $+ let binlist = reverse $ sortBy (compare `on` fst) $+ map (\ (tr,OneCluster ls) -> (length ls, 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)++")"-+ ClusterThem{linkage} -> do+ (mat, dendro) <- doCluster use_hashrf num_taxa linkage validtrees + when print_rfmatrix $ printDistMat stdout mat hnd <- openFile (combine output_dir ("distance_matrix.txt")) WriteMode printDistMat hnd mat hClose hnd-+ --------------------+ writeFile (combine output_dir ("dendrogram.txt"))+ (show$ fmap treename dendro)+ putStrLn " [finished] Wrote full dendrogram to file dendrogram.txt"+ let plotIt mnameMap = + if True -- do_graph+ then async (do + t0 <- getCurrentTime+ let dot = dotDendrogram cfg "dendrogram" 1.0 dendro mnameMap highlightTrs+ _ <- dotToPDF dot (combine output_dir "dendrogram.pdf") + t1 <- getCurrentTime + putStrLn$ " [finished] Wrote dendrogram diagram to file dendrogram.pdf ("+ ++show(diffUTCTime t1 t0)++")")+ else async (return ()) case dist_thresh of- Nothing -> error "Fully hierarchical cluster output is not finished! Use --editdist."- Just dstThresh -> do+ Nothing -> do a <- plotIt Nothing+ return (M.empty,[],[a])+ -- 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)+ classes = clustsToMap clusts+ -- Cache the lengths:+ wlens = [ (length ls, OneCluster ls) | OneCluster ls <- clusts ]+ sorted0 = reverse$ sortBy (compare `on` fst) wlens -- TODO: Parallel sorting?+ nameMap = clustsToNameMap (map snd sorted0)+ gvizAsync <- plotIt (Just nameMap)+ putStrLn$ "After flattening, cluster sizes are: "++show (map fst sorted0) -- Flatten out the dendogram:- return (clustsToMap clusts, sorted0, [gvizAsync])+ return (classes, sorted0, [gvizAsync]) - reportClusts clust_mode binlist+ unless (null binlist)$ reportClusts clust_mode binlist ---------------------------------------- -- TEST, TEMPTOGGLE: print out edge weights :@@ -212,10 +230,7 @@ -- 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)+-- putStrLn$ "Final number of tree bins: "++ show (M.size classes) unless (null warnings1 && null warnings2) $ writeFile (combine output_dir "WARNINGS.txt")@@ -230,13 +245,13 @@ "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+ async2 <- case clust_mode of+ BinThem -> outputBins binlist output_dir do_graph+ ClusterThem{} -> outputClusters num_taxa binlist output_dir do_graph -- Wait on parallel tasks: putStrLn$ "Waiting for asynchronous tasks to finish..."- mapM_ wait asyncs + mapM_ wait (async2:asyncs) putStrLn$ "Finished." -------------------------------------------------------------------------------- -- End driver@@ -257,12 +272,15 @@ 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+doCluster :: Bool -> Int -> C.Linkage -> [FullTree a] -> IO (DistanceMatrix, C.Dendrogram (FullTree a))+doCluster use_hashrf num_taxa linkage validtrees = do+ t0 <- getCurrentTime+ when use_hashrf$ putStrLn " Using HashRF-style algorithm..." let nwtrees = map nwtree validtrees- numtrees = length validtrees - mat = distanceMatrix nwtrees+ numtrees = length validtrees+ mat = if use_hashrf + then hashRF num_taxa nwtrees+ else fst (distanceMatrix nwtrees) ixtrees = zip [0..] validtrees dist (i,t1) (j,t2) | j == i = 0 -- | i == numtrees-1 = 0 @@ -270,29 +288,29 @@ | 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)+ -- Force the distance matrix:+ V.mapM_ evaluate mat+ t1 <- getCurrentTime+ putStrLn$ "Time to compute distance matrix: "++show(diffUTCTime t1 t0)+ putStrLn$ "Clustering using method "++show linkage return (mat,dendro)- --- reportClusts :: ClustMode -> BinResults StandardDecor -> IO [(Int, StrippedTree, OneCluster StandardDecor)]+++reportClusts :: ClustMode -> [(Int, OneCluster StandardDecor)] -> IO () reportClusts mode binlist = do - let - taxa :: S.Set Int- taxa = S.unions$ map (S.fromList . all_labels . snd3) binlist- binsizes = map fst3 binlist+ let binsizes = map fst 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+ putStrLn$" Up to first 30 bin sizes, excluding singletons:"+ forM_ (zip [1..30] binlist) $ \ (ind, (len, 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+ ClusterThem{} -> hcat [] ] -- | Convert a flat list of clusters into a map from individual trees to clusters. clustsToMap :: [OneCluster StandardDecor] -> BinResults StandardDecor@@ -303,6 +321,14 @@ fn2 theclust acc (FullTree{nwtree}) = M.insert (anonymize_annotated nwtree) theclust acc +-- | Map each tree NAME onto the one-based index (in sorted order) of the cluster it+-- comes from.+clustsToNameMap :: [OneCluster StandardDecor] -> M.Map TreeName Int+clustsToNameMap clusts = foldl' fn M.empty (zip [1..] clusts)+ where+ fn acc (ix,(OneCluster ftrs)) =+ foldl' (\acc' t -> M.insert (treename t) ix acc') acc ftrs+ flattenDendro :: (C.Dendrogram (FullTree DefDecor)) -> OneCluster StandardDecor flattenDendro dendro = case dendro of@@ -327,65 +353,72 @@ -------------------------------------------------------------------------------- -outputClusters binlist output_dir do_graph = do+-- outputClusters :: (Num a1, Ord a1, Show a1) => Int -> [(a1, t, OneCluster a)] -> String -> Bool -> IO ()+outputClusters :: Int -> [(Int, OneCluster a)] -> String -> Bool -> IO (Async ())+outputClusters num_taxa 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+ let consTrs = [ FullTree "consensus" (labelTable $ head ftrees) nwtr + | (_, OneCluster ftrees) <- binlist+ , let nwtr :: NewickTree StandardDecor+ nwtr = annotateWLabLists $ fmap (const (Nothing,0)) $+ consensusTree num_taxa $ map nwtree ftrees ]+ + forM_ (zip3 [1::Int ..] binlist consTrs) $ \ (i, (size, OneCluster ftrees), fullConsTr) -> 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+ writeFile (base i size ++"_consensus.tr") (show (displayStrippedTree fullConsTr) ++ "\n")+ writeFile (base i size ++"_alltrees.tr")+ (unlines [ show (displayStrippedTree ft) | ft <- ftrees ]) - 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+ putStrLn$ " [finished] Wrote contents of each cluster to cluster<N>_<size>.txt"+ putStrLn$ " [finished] Wrote representative (consensus) trees to cluster<N>_<size>_consensus.tr"+ if do_graph then do+ putStrLn$ "Next start the time consuming operation of writing out graphviz visualizations:"+ asyncs <- forM (zip3 [1::Int ..] binlist consTrs) $+ \ (i, (size, OneCluster membs), fullConsTr) -> async$ 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")+ let dot = dotNewickTree ("cluster #"++ show i) 1.0 fullConsTr+ _ <- dotToPDF dot (base i size ++ ".pdf") return ()- putStrLn$ "[finished] Wrote visual representations of trees to "++filePrefix++"<N>_<size>.pdf"-- return ()-+ async $ do+ mapM_ wait asyncs+ putStrLn$ " [finished] Wrote visual representations of consensus trees to "++filePrefix++"<N>_<size>.pdf"+ else async (return ()) +outputBins :: [(Int, OneCluster StandardDecor)] -> String -> Bool -> IO (Async ()) 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 avgs = map (avg_trees . map nwtree . clustMembers . snd) binlist+ forM_ (zip3 [1::Int ..] binlist avgs) $ \ (i, (size, 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"-+ -- Printing the average tree instead of the stripped cannonical one: + writeFile (base i size ++"_avg.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>_avg.tr" + if do_graph then do+ putStrLn$ "Next start the time consuming operation of writing out graphviz visualizations:"+ asyncs <- forM (zip3 [1::Int ..] binlist avgs) $ \ (i, (size, OneCluster membs), avgTree) -> async $ 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 ()+ async $ do+ mapM_ wait asyncs+ putStrLn$ " [finished] Wrote visual representations of trees to "++filePrefix++"<N>_<size>.pdf"+ else async (return ()) -------------------------------------------------------------------------------- -- Monadic mapAccum@@ -436,3 +469,30 @@ type TempDecor = (Double, (Int, Int), Int, [Label]) avg ls = sum ls / fromIntegral (length ls)++-- | Parse trees in addition to the main inputs (for --highlight).+retrieveHighlights :: (String->String) -> LabelTable -> [FilePath] -> IO [[NewickTree ()]]+retrieveHighlights name_hack labelTab ls =+ mapM parseHighlight ls+ where+ parseHighlight file = do + bs <- B.readFile file+ let (lt2,htr) = parseNewick labelTab name_hack file bs+ unless (lt2 == labelTab) $+ error$"Tree given as --highlight includes taxa not present in main tree set: "+++ show(M.keys$ M.difference lt2 labelTab) + return (map (fmap (const())) htr)+++-- | Create a predicate that tests trees for consistency with the set of --highlight+-- (consensus) trees.+matchAnyHighlight :: [[NewickTree ()]] -> NewickTree () -> Bool +-- matchAnyHighlight :: [[NewickTree ()]] -> NewickTree () -> Maybe Int+-- If there is a match, return the index of the highlight that matched.+matchAnyHighlight highlightTrs =+ let matchers = map mkMatcher highlightTrs+ mkMatcher ls = let fns = map compatibleWith ls -- Multiple predicate functions+ in \ tr -> -- Does it match any predicate?+ or$ map (\f -> f tr) fns + in \ newtr -> + any (\f -> f newtr) matchers
Bio/Phylogeny/PhyBin/Binning.hs view
@@ -30,10 +30,8 @@ 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)@@ -179,13 +177,13 @@ -- 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);"+tt = normalize $ annotateWLabLists $ head$ 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));"+norm5 = normalize$ annotateWLabLists$ head$ snd$ parseNewick M.empty id "" "(D,E,C,(B,A));" ----------------------------------------------------------------------------------------------------@@ -332,6 +330,8 @@ annotateWLabLists :: NewickTree DefDecor -> AnnotatedTree annotateWLabLists tr = case tr of NTLeaf (bs,bl) n -> NTLeaf (StandardDecor bl bs 1 [n]) n+ NTInterior (bs,bl) [] -> NTInterior (StandardDecor bl bs 0 []) []+ -- NTInterior (bs,bl) [] -> error "annotateWLabLists: internal invariant broken. Shouldn't have NTInterior with null children." NTInterior (bs,bl) ls -> let children = map annotateWLabLists ls in NTInterior (StandardDecor bl bs@@ -372,7 +372,9 @@ ---------------------------------------------------------------------------------------------------- tre1 :: (LabelTable, NewickTree DefDecor)-tre1 = parseNewick M.empty id "" "(A:0.1,B:0.2,(C:0.3,D:0.4):0.5);"+tre1 = (x,head y)+ where+ (x,y) = 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)))@@ -384,7 +386,7 @@ norm = norm2 . B.pack norm2 :: B.ByteString -> FullTree StandardDecor-norm2 bstr = FullTree "" tbl (normalize $ annotateWLabLists tr)+norm2 bstr = FullTree "" tbl (normalize $ annotateWLabLists$ head tr) where (tbl,tr) = parseNewick M.empty id "test" bstr
Bio/Phylogeny/PhyBin/CoreTypes.hs view
@@ -12,7 +12,7 @@ ClustMode(..), TreeName, -- * Tree operations- displayDefaultTree,+ displayDefaultTree, displayStrippedTree, treeSize, numLeaves, liftFT, get_dec, set_dec, get_children, map_labels, all_labels, foldIsomorphicTrees,@@ -21,7 +21,7 @@ avg_branchlen, get_bootstraps, -- * Command line config options- PhyBinConfig(..), default_phybin_config,+ PhyBinConfig(..), default_phybin_config, -- * General helpers Label, LabelTable,@@ -85,6 +85,9 @@ instance Foldable FullTree where foldMap f (FullTree _ _ tr) = foldMap f tr +instance Functor FullTree where+ fmap f (FullTree n l tr) = FullTree n l $ fmap 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:@@ -114,6 +117,13 @@ Just val -> base <> text ":[" <> text (show val) <> text "]" where base = parens$ sep$ map_but_last (<>text",") $ map loop ls +displayStrippedTree :: FullTree a -> Doc+displayStrippedTree orig = loop tr <> ";"+ where+ (FullTree _ mp tr) = orig -- normalize orig+ loop (NTLeaf _ name) = text (mp M.! name)+ loop (NTInterior _ ls) = parens$ sep$ map_but_last (<>text",") $ map loop ls+ ---------------------------------------------------------------------------------------------------- -- Labels ----------------------------------------------------------------------------------------------------@@ -177,7 +187,7 @@ , labelTable :: LabelTable , nwtree :: NewickTree a }- deriving (Show)+ deriving (Show, Ord, Eq) liftFT :: (NewickTree t -> NewickTree a) -> FullTree t -> FullTree a liftFT fn (FullTree nm labs x) = FullTree nm labs (fn x)@@ -205,6 +215,10 @@ , do_graph :: Bool , do_draw :: Bool , clust_mode :: ClustMode+ , highlights :: [FilePath] -- [FullTree ()]+ , show_trees_in_dendro :: Bool+ , show_interior_consensus :: Bool+ , use_hashrf :: Bool , print_rfmatrix :: Bool , dist_thresh :: Maybe Int , branch_collapse_thresh :: Maybe Double -- ^ Branches less than this length are collapsed.@@ -214,19 +228,28 @@ 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.)"+ , num_taxa = error "must be able to determine the number of taxa expected in the dataset. (Supply it manually with -n.)" , name_hack = id -- Default, no transformation of leaf-labels , output_dir = "./phybin_out/" , inputs = [] , do_graph = False , do_draw = False- , clust_mode = BinThem+-- , clust_mode = BinThem+ , clust_mode = ClusterThem C.UPGMA+#ifdef USE_HASHRF+ , use_hashrf = True+#else+ , use_hashrf = False+#endif+ , highlights = []+ , show_trees_in_dendro = False+ , show_interior_consensus = False , print_rfmatrix = False , dist_thresh = Nothing , branch_collapse_thresh = Nothing } -data ClustMode = BinThem | ClusterThem C.Linkage +data ClustMode = BinThem | ClusterThem { linkage :: C.Linkage } ---------------------------------------------------------------------------------------------------- -- * Simple utility functions for the core types:
Bio/Phylogeny/PhyBin/Parser.hs view
@@ -24,12 +24,24 @@ -- | 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 :: LabelTable -> NameHack -> String -> B.ByteString -> (LabelTable, [NewickTree DefDecor]) parseNewick tbl0 name_hack file input = extractLabelTable tbl0 $ - runB file (newick_parser name_hack) $+ runB file (many1$ newick_parser name_hack) $ B.filter (not . isSpace) input +-- 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))+ -- | 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@@ -41,34 +53,43 @@ parseNewicks name_hack pairs = (labtbl, fullTrs) where fullTrs = [ FullTree (takeBaseName file) labtbl tr- | (file,_) <- pairs- | tr <- trs ]+ | (file,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)+ let (acc',trs) = parseNewick acc name_hack file bstr+ names = if length trs > 1+ then zipWith (\x y -> x++"_"++show y) (repeat file) [0..]+ else [file]+ in (acc', (zip names trs) ++ 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)+extractLabelTable :: LabelTable -> [TempTree] -> (LabelTable, [NewickTree DefDecor])+extractLabelTable tbl0 trs = (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) _)+ -- (_,finMap,finTree) = loop2 (S.fromList (M.elems tbl0)) tbl0 tr+ (_,finMap,finTree) = loop1 flipped tbl0 trs++ -- Mere plumbing:+ loop1 seen acc [] = (seen, acc, [])+ loop1 seen acc (hd:tl) =+ let (seen2,acc2,hd') = loop2 seen acc hd + (seen3,acc3,tl') = loop1 seen2 acc2 tl in+ (seen3,acc3, hd':tl')++ loop2 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) =+ loop2 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+ let (seen3,acc3,x') = loop2 seen2 acc2 x in (seen3, acc3, x':ls)) (seen1,acc1,[]) chlds@@ -173,7 +194,7 @@ -------------------------------------------------------------------------------- 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);"+tre1 = head $ 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
Bio/Phylogeny/PhyBin/RFDistance.hs view
@@ -1,20 +1,36 @@-{-# LANGUAGE ScopedTypeVariables, CPP #-}+{-# LANGUAGE ScopedTypeVariables, CPP, BangPatterns #-} module Bio.Phylogeny.PhyBin.RFDistance- (DenseLabelSet, DistanceMatrix, - allBips, foldBips,- distanceMatrix, printDistMat)+ (+ -- * Types+ DenseLabelSet, DistanceMatrix,++ -- * Bipartition (Bip) utilities+ allBips, foldBips, dispBip,+ consensusTree, bipsToTree, filterCompatible, compatibleWith,++ -- * ADT for dense sets+ mkSingleDense, mkEmptyDense, bipSize,+ denseUnions, denseDiff, invertDense, markLabel,+ + -- * Methods for computing distance matrices+ distanceMatrix, hashRF, ++ -- * Output+ printDistMat) where import Control.Monad+import Control.Monad.ST+import Data.Function (on) import Data.Word import qualified Data.Vector as V-import qualified Data.Vector.Unboxed.Mutable as MV+import qualified Data.Vector.Mutable as MV+import qualified Data.Vector.Unboxed.Mutable as MU 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 System.IO.Unsafe -- import Control.LVish -- import qualified Data.LVar.Set as IS@@ -26,17 +42,33 @@ import Bio.Phylogeny.PhyBin.CoreTypes -- import Data.BitList import qualified Data.Set as S+import qualified Data.List as L import qualified Data.IntSet as SI import qualified Data.Map.Strict as M import qualified Data.Foldable as F+import qualified Data.Traversable as T import Data.Monoid import Prelude as P import Debug.Trace +#ifdef BITVEC_BIPS+import qualified Data.Vector.Unboxed.Bit as UB+import qualified Data.Bit as B+#endif++-- I don't understand WHY, but I seem to get the same answers WITHOUT this.+-- Normalization and symmetric difference do make things somewhat slower (e.g. 1.8+-- seconds vs. 2.2 seconds for 150 taxa / 100 trees)+#define NORMALIZATION+-- define BITVEC_BIPS+ -------------------------------------------------------------------------------- -- A data structure choice -------------------------------------------------------------------------------- +-- type DenseLabelSet s = BitList++ -- | 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). -- @@ -49,27 +81,72 @@ -- -- 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+#ifdef BITVEC_BIPS --- M.write vec lab (B.fromBool True)--- mkEmptyDense size = U.replicate size (B.fromBool False) +# if 1+type DenseLabelSet = UB.Vector B.Bit+markLabel lab = UB.modify (\vec -> MU.write vec lab (B.fromBool True)) +mkEmptyDense size = UB.replicate size (B.fromBool False)+mkSingleDense size ind = markLabel ind (mkEmptyDense size)+denseUnions = UB.unions+bipSize = UB.countBits+denseDiff = UB.difference+invertDense size bip = UB.invert bip+dispBip labs bip = show$ map (\(ix,_) -> (labs M.! ix)) $+ filter (\(_,bit) -> B.toBool bit) $+ zip [0..] (UB.toList bip)+denseIsSubset a b = UB.or (UB.difference b a)+traverseDense_ fn bip =+ U.ifoldr' (\ ix bit acc ->+ (if B.toBool bit+ then fn ix+ else return ()) >> acc)+ (return ()) bip --- markLabel lab set = IS.putInSet lab set --- mkEmptyDense _size = IS.newEmptySet+# else+-- TODO: try tracking the size:+data DenseLabelSet = DLS {-# UNPACK #-} !Int (UB.Vector B.Bit)+markLabel lab (DLS _ vec)= DLS (UB.modify (\vec -> return (MU.write vec lab (B.fromBool True))) ) vec+-- ....+# endif +#else+type DenseLabelSet = SI.IntSet+markLabel lab set = SI.insert lab set +mkEmptyDense _size = SI.empty+mkSingleDense _size = SI.singleton+denseUnions _size = SI.unions +bipSize = SI.size+denseDiff = SI.difference+denseIsSubset = SI.isSubsetOf++dispBip labs bip = "[" ++ unwords strs ++ "]"+ where strs = map (labs M.!) $ SI.toList bip+invertDense size bip = loop SI.empty (size-1)+ where -- There's nothing for it but to iterate and test for membership:+ loop !acc ix | ix < 0 = acc+ | SI.member ix bip = loop acc (ix-1)+ | otherwise = loop (SI.insert ix acc) (ix-1)+traverseDense_ fn bip =+ -- FIXME: need guaranteed non-allocating way to do this.+ SI.foldr' (\ix acc -> fn ix >> acc) (return ()) bip+#endif+ markLabel :: Label -> DenseLabelSet -> DenseLabelSet mkEmptyDense :: Int -> DenseLabelSet+mkSingleDense :: Int -> Label -> 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+-- | Print a BiPartition in a pretty form+dispBip :: LabelTable -> DenseLabelSet -> String +-- | Assume that total taxa are 0..N-1, invert membership:+invertDense :: Int -> DenseLabelSet -> DenseLabelSet +traverseDense_ :: Monad m => (Int -> m ()) -> DenseLabelSet -> m ()++ -------------------------------------------------------------------------------- -- Dirt-simple reference implementation --------------------------------------------------------------------------------@@ -77,20 +154,26 @@ type DistanceMatrix = V.Vector (U.Vector Int) -- | Returns a triangular distance matrix encoded as a vector.-distanceMatrix :: [NewickTree a] -> DistanceMatrix+-- Also return the set-of-BIPs representation for each tree.+distanceMatrix :: [NewickTree a] -> (DistanceMatrix, V.Vector (S.Set DenseLabelSet)) 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))- + mat = V.generate sz $ \ i -> + U.generate i $ \ j ->+ let diff1 = S.size (S.difference (eachbips V.! i) (eachbips V.! j))+ diff2 = S.size (S.difference (eachbips V.! j) (eachbips V.! i))+ in diff1 + diff2+ in (mat, eachbips)+ -- | 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) $+-- trace ("labelbips "++show allLeaves++" "++show size) $+#ifdef NORMALIZATION + fmap (\(a,ls) -> (a,map (normBip size) ls)) $+#endif loop tr where size = numLeaves tr@@ -98,22 +181,28 @@ 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+ sets = map (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 (NTLeaf _ lab) = mkSingleDense size lab leafSet (NTInterior _ ls) = denseUnions size $ map leafSet ls +-- normBip :: DenseLabelSet -> DenseLabelSet -> DenseLabelSet+-- normBip allLeaves bip =+normBip :: Int -> DenseLabelSet -> DenseLabelSet +normBip totsize bip =+ let -- size = bipSize allLeaves+ halfSize = totsize `quot` 2+-- flipped = denseDiff allLeaves bip+ flipped = invertDense totsize bip + in + case compare (bipSize bip) halfSize of+ LT -> bip + GT -> flipped -- Flip it+ EQ -> min bip flipped -- This is a painful case, we need a tie-breaker+ + foldBips :: Monoid m => (DenseLabelSet -> m) -> NewickTree a -> m foldBips f tr = F.foldMap f' (labelBips tr) where f' (_,bips) = F.foldMap f bips@@ -127,6 +216,7 @@ -------------------------------------------------------------------------------- -- First, necessary types: +-- UNFINISHED: #if 0 -- | A collection of all observed bipartitons (bips) with a mapping of which trees -- contain which Bips.@@ -143,12 +233,78 @@ -- | Tree's are identified simply by their order within the list of input trees. -- type TreeID = Int #endif+ --------------------------------------------------------------------------------+-- Alternate way of slicing the problem: HashRF+-------------------------------------------------------------------------------- -- 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. +-- | This version slices the problem a different way. A single pass over the trees+-- populates the table of bipartitions. Then the table can be processed (locally) to+-- produce (non-localized) increments to a distance matrix.+hashRF :: Int -> [NewickTree a] -> DistanceMatrix+hashRF num_taxa trees = build M.empty (zip [0..] trees)+ where+ num_trees = length trees+ -- First build the table:+ build acc [] = ingest acc+ build acc ((ix,hd):tl) =+ let bips = allBips hd+ acc' = S.foldl' fn acc bips+ fn acc bip = M.alter fn2 bip acc+ fn2 (Just membs) = Just (markLabel ix membs)+ fn2 Nothing = Just (mkSingleDense num_taxa ix)+ in + build acc' tl++ -- Second, ingest the table to construct the distance matrix:+ ingest :: M.Map DenseLabelSet DenseLabelSet -> DistanceMatrix+ ingest bipTable = runST theST+ where+ theST :: forall s0 . ST s0 DistanceMatrix+ theST = do + -- Triangular matrix, starting narrow and widening:+ matr <- MV.new num_trees+ -- Too bad MV.replicateM is insufficient. It should pass index. + -- Instead we write this C-style:+ for_ (0,num_trees) $ \ ix -> do + row <- MU.replicate ix (0::Int)+ MV.write matr ix row+ return ()++ unsafeIOToST$ putStrLn$" Built matrix for dim "++show num_trees++ let bumpMatr i j | j < i = incr i j+ | otherwise = incr j i+ incr :: Int -> Int -> ST s0 ()+ incr i j = do -- Not concurrency safe yet:+-- unsafeIOToST$ putStrLn$" Reading at position "++show(i,j)+ row <- MV.read matr i+ elm <- MU.read row j+ MU.write row j (elm+1)+ return ()+ fn bipMembs =+ -- Here we quadratically consider all pairs of trees and ask whether+ -- their edit distance is increased based on this particular BiP.+ -- Actually, as an optimization, it is sufficient to consider only the+ -- cartesian product of those that have and those that don't.+ let haveIt = bipMembs+ -- Depending on how invertDense is written, it could be useful to+ -- fuse this in and deforest "dontHave".+ dontHave = invertDense num_trees bipMembs+ fn1 trId = traverseDense_ (fn2 trId) dontHave+ fn2 trId1 trId2 = bumpMatr trId1 trId2+ in+-- trace ("Computed donthave "++ show dontHave) $ + traverseDense_ fn1 haveIt+ F.traverse_ fn bipTable+ v1 <- V.unsafeFreeze matr+ T.traverse (U.unsafeFreeze) v1++ #if 0 -- | Returns a (square) distance matrix encoded as a vector. distanceMatrix :: [AnnotatedTree] -> IO (U.Vector Word)@@ -157,7 +313,7 @@ -- 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)+ v <- MU.replicate (sz*sz) (0::Word) let fn set () = F.foldrM @@ -181,6 +337,8 @@ #endif --------------------------------------------------------------------------------+-- Miscellaneous Helpers+-------------------------------------------------------------------------------- instance Pretty a => Pretty (S.Set a) where pPrint s = pPrint (S.toList s)@@ -196,3 +354,82 @@ hPutStr h " " hPutStr h "0\n" hPutStrLn h "-----------------------------------------"++-- My own forM for numeric ranges (not requiring deforestation optimizations).+-- Inclusive start, exclusive end.+{-# INLINE for_ #-}+for_ :: Monad m => (Int, Int) -> (Int -> m ()) -> m ()+for_ (start, end) _fn | start > end = error "for_: start is greater than end"+for_ (start, end) fn = loop start+ where+ loop !i | i == end = return ()+ | otherwise = do fn i; loop (i+1)++-- | Which of a set of trees are compatible with a consensus?+filterCompatible :: NewickTree a -> [NewickTree b] -> [NewickTree b]+filterCompatible consensus trees =+ let cbips = allBips consensus in+ [ tr | tr <- trees+ , cbips `S.isSubsetOf` allBips tr ]++-- | Is a tree compatible with a consensus?+-- This is more efficient if partially applied then used repeatedly.+compatibleWith :: NewickTree a -> NewickTree b -> Bool+compatibleWith consensus newTr =+ S.isSubsetOf (allBips consensus) (allBips newTr)++-- Consensus between two trees, which may even have different label maps.+consensusTreeFull (FullTree n1 l1 t1) (FullTree n2 l2 t2) =+ error "FINISHME - consensusTreeFull"++-- | Take only the bipartitions that are agreed on by all trees.+consensusTree :: Int -> [NewickTree a] -> NewickTree ()+consensusTree _ [] = error "Cannot take the consensusTree of the empty list"+consensusTree num_taxa (hd:tl) = bipsToTree num_taxa intersection+ where+ intersection = L.foldl' S.intersection (allBips hd) (map allBips tl)+-- intersection = loop (allBips hd) tl+-- loop :: S.Set DenseLabelSet -> [NewickTree a] -> S.Set DenseLabelSet+-- loop !remain [] = remain+-- -- Was attempting to use foldBips here as an optimization:+-- -- loop !remain (hd:tl) = loop (foldBips S.delete hd remain) tl+-- loop !remain (hd:tl) = loop (S.difference remain (allBips hd)) tl + +-- | Convert from bipartitions BACK to a single tree.+bipsToTree :: Int -> S.Set DenseLabelSet -> NewickTree ()+bipsToTree num_taxa origbip =+-- trace ("Doing bips in order: "++show sorted++"\n") $ + loop lvl0 sorted+ where+ -- We consider each subset in increasing size order.+ -- FIXME: If we tweak the order on BIPs, then we can just use S.toAscList here:+ sorted = L.sortBy (compare `on` bipSize) (S.toList origbip)++ lvl0 = [ (mkSingleDense num_taxa ix, NTLeaf () ix)+ | ix <- [0..num_taxa-1] ]++ -- VERY expensive! However, due to normalization issues this is necessary for now:+ -- TODO: in the future make it possible to definitively denormalize.+ -- isMatch bip x = denseIsSubset x bip || denseIsSubset x (invertDense num_taxa bip)+ isMatch bip x = denseIsSubset x bip ++ -- We recursively glom together subtrees until we have a complete tree.+ -- We only process larger subtrees after we have processed all the smaller ones.+ loop !subtrees [] =+ case subtrees of+ [] -> error "bipsToTree: internal error"+ [(_,one)] -> one+ lst -> NTInterior () (map snd lst)+ loop !subtrees (bip:tl) =+-- trace (" -> looping, subtrees "++show subtrees) $ + let (in_,out) = L.partition (isMatch bip. fst) subtrees in+ case in_ of+ [] -> error $"bipsToTree: Internal error! No match for bip: "++show bip+ ++" out is\n "++show out++"\n and remaining bips "++show (length tl)+ ++"\n when processing orig bip set:\n "++show origbip+ -- loop out tl+ _ -> + -- Here all subtrees that match the current bip get merged:+ loop ((denseUnions num_taxa (map fst in_),+ NTInterior () (map snd in_)) : out) tl+
Bio/Phylogeny/PhyBin/Util.hs view
@@ -27,11 +27,10 @@ import System.FilePath (combine) import System.Directory (doesFileExist, doesDirectoryExist, getDirectoryContents, getCurrentDirectory)-import System.IO (openFile, hClose, IOMode(ReadMode))-import System.Process (system)+import System.IO (openFile, hClose, IOMode(ReadMode), stderr,+ hPutStr, hPutStrLn) import System.Exit (ExitCode(..)) import Test.HUnit ((~:),(~=?),Test,test)-import qualified HSH -- For vizualization: import Text.PrettyPrint.HughesPJClass hiding (char, Style)@@ -89,19 +88,20 @@ --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+ if not exists then do+ error$ "No file or directory found at this path!: "++path+ -- hPutStr stderr$ "Input not a file/directory, assuming wildcard, using 'find' for expansion"+ -- entries <- HSH.run$ "find " ++ path + -- hPutStrLn stderr$ "("++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 "+ hPutStr stderr$ "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+ hPutStrLn stderr$ "("++show (length filtered)++" regular files found): "++ show path return$ filtered else if reg then do return [path]
Bio/Phylogeny/PhyBin/Visualize.hs view
@@ -10,8 +10,11 @@ where import Text.Printf (printf) import Data.List (elemIndex, isPrefixOf)-import Data.Maybe (fromJust)-import Data.Map ((!))+import Data.List.Split (chunksOf)+import Data.Maybe (fromJust, isJust)+import qualified Data.Map as M+import qualified Data.Set as S+import qualified Data.Vector as V import Data.Text.Lazy (pack) import Control.Monad (void) import Control.Concurrent (Chan, newChan, writeChan, forkIO)@@ -19,11 +22,13 @@ 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 Data.GraphViz.Attributes.Colors (Color(RGB)) -- import Test.HUnit ((~:),(~=?),Test,test) import qualified Data.Clustering.Hierarchical as C import Bio.Phylogeny.PhyBin.CoreTypes+import Bio.Phylogeny.PhyBin.RFDistance (filterCompatible, compatibleWith, consensusTree) ---------------------------------------------------------------------------------------------------- -- Visualization with GraphViz and FGL:@@ -33,7 +38,7 @@ toGraph :: FullTree StandardDecor -> G.Gr String Double toGraph (FullTree _ tbl tree) = G.run_ G.empty $ loop tree where- fromLabel ix = tbl ! ix+ fromLabel ix = tbl M.! ix loop (NTLeaf _ name) = do let str = fromLabel name _ <- G.insMapNodeM str@@ -59,24 +64,171 @@ 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+----------------------------------------------------------------------------------------------------+-- Dendrograms+ ----------------------------------------------------------------------------------------------------++-- | Some duplicated code with dotNewickTree.+dotDendrogram :: PhyBinConfig -> String -> Double -> C.Dendrogram (FullTree a) ->+ Maybe (M.Map TreeName Int) -> [[NewickTree ()]] -> Gv.DotGraph G.Node+dotDendrogram PBC{show_trees_in_dendro, show_interior_consensus}+ title edge_scale origDendro mNameMap highlightTrs =+ Gv.graphToDot myparams (G.nmap uid graph)+ where+ (charsDropped, dendro) = truncateNames origDendro+ -- This is ugly, but we modify the name map to match:+ nameMap' = fmap (M.mapKeys (drop charsDropped)) mNameMap + graph :: DendroGraph+ graph = dendrogramToGraph num_taxa dendro++ num_taxa = numLeaves $ nwtree fstLeaf+ FullTree{labelTable} = fstLeaf+ fstLeaf = getFstLeaf dendro+ getFstLeaf (C.Leaf ft) = ft+ getFstLeaf (C.Branch _ l _) = getFstLeaf l ++ uidsToNames = M.fromList $+ map (\nd@NdLabel{uid} -> (uid,nd)) $+-- filter (isJust . tre) $ + map (fromJust . G.lab graph) $ + G.nodes graph++ matchers = map mkMatcher highlightTrs+ mkMatcher ls = let fns = map compatibleWith ls -- Multiple predicate functions+ in \ tr -> -- Does it match any predicate?+ or$ map (\f -> f tr) fns + + wcolors = zip matchers defaultPalette+ findColor tr = loop wcolors+ where loop [] = Nothing+ loop ((f,col):rst) | f tr = Just col+ | otherwise = loop rst+ + -- Map uids to highlight color+ highlightMap = M.map fromJust $+ M.filter isJust $+ M.map (\ (NdLabel _ _ _ ct) -> findColor ct)+ uidsToNames++ 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, UniqueNodeName) -> [GA.Attribute]+ nodeAttrs (_num, uid) =+ let uid' = if show_trees_in_dendro+ then printed_tree++"\n"++uid+ else uid+ printed_tree =+ case M.lookup uid uidsToNames of+ Nothing -> ""+ Just NdLabel{clumpSize,consensus} ->+ (if clumpSize>1 then "size "++show clumpSize++"\n" else "")+ ++ show (displayStrippedTree (FullTree "" labelTable consensus))+ (tag,shp,styl) = -- case eith of+ if isPrefixOf "DUMMY_" uid+ then (if show_trees_in_dendro && show_interior_consensus+ then printed_tree else "",+ if show_interior_consensus+ then GA.BoxShape -- GA.PlainText+ else GA.PointShape,+ [ GA.Color [weighted$ GA.X11Color Gv.Transparent] ]+ )+ else (uid', GA.Ellipse, [ GA.Style [GA.SItem GA.Filled []]])+ highlightColor = + case M.lookup uid highlightMap of+ Nothing -> []+ Just col -> [ GA.Color [weighted col ] ]+ clustColor | not (null highlightTrs) = []+ | otherwise = + case (nameMap', M.lookup uid uidsToNames) of+ (Just nm, Just NdLabel{tre=Just FullTree{treename}}) ->+ case M.lookup treename nm of+ Nothing -> []+ -- Here we color the top TEN clusters:+ Just ind | ind <= 10 -> [ GA.Color [weighted$ defaultPaletteV V.! (ind-1) ] ]+ | otherwise -> []+ -- TODO: shall we color intermediate nodes?+ _ -> []+ in + [ GA.Label$ GA.StrLabel$ pack tag+ , GA.Shape shp+ ] ++ styl ++ highlightColor ++ clustColor++ edgeAttrs = getEdgeAttrs edge_scale+++type UniqueNodeName = String+-- type DendroGraph = G.Gr (Maybe TreeName,UniqueNodeName) Double+type DendroGraph = G.Gr NdLabel Double++-- | When we first convert to a graph representation, there is a bunch of information+-- hanging off of each node.+data NdLabel =+ NdLabel+ { uid :: UniqueNodeName+ , tre :: Maybe (FullTree ())+ , clumpSize :: !Int+ , consensus :: NewickTree ()+ }+ deriving (Show, Ord, Eq)++-- | Create a graph using TreeNames for node labels and edit-distance for edge labels.+dendrogramToGraph :: Int -> C.Dendrogram (FullTree a) -> DendroGraph+dendrogramToGraph num_taxa orig =+ G.run_ G.empty $ void$+ loop (fmap (fmap (const())) orig) where- -- deEither (Left s) = s- -- deEither (Right s) = s- loop node@(C.Leaf FullTree{treename}) = G.insMapNodeM (treename)+-- loop node@(C.Leaf ft) = G.insMapNodeM (NdLabel (treename ft) (Just$ fmap (const()) ft) 1 ft)+ loop node@(C.Leaf ft) =+ let stripped = fmap (const()) ft in+ G.insMapNodeM (NdLabel (treename ft) (Just stripped) 1 (nwtree stripped))+ loop node@(C.Branch 0 left right) = do+ -- As a preprocessing step we collapse clusters that are separated by zero edit distance.+ ----------------------------------------+ let lvs = collapseZeroes left ++ collapseZeroes right+ nms = map (treename) lvs+ lens = map length nms+ total = sum lens+ avg = total `quot` length nms+ -- The goal here is to make an approximately square arrangement:+ -- goal: avg * perline == total / (avg * perline)+ perline = ceiling$ sqrt (fromIntegral total / ((fromIntegral avg)^2))+ chunked = chunksOf perline nms+ fatname = unlines (map unwords chunked)+ G.insMapNodeM (NdLabel fatname (Just (head lvs))+ (length lvs) (consensusTree num_taxa (map nwtree lvs)))+ ---------------------------------------- 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)+ do (_,ll@(NdLabel lid _ s1 c1)) <- loop left+ (_,rr@(NdLabel rid _ s2 c2)) <- loop right+ -- Interior nodes do NOT have their names drawn:+ let ndname = "DUMMY_"++(lid++"_"++rid) -- HACK!+ (midN,mid) <- G.insMapNodeM (NdLabel ndname Nothing (s1+s2) (consensusTree num_taxa [c1,c2]))+ G.insMapEdgeM (ll, mid, dist)+ G.insMapEdgeM (rr, mid, dist) return (midN,mid) + collapseZeroes (C.Leaf tr) = [tr]+ collapseZeroes (C.Branch 0 l r) = collapseZeroes l ++ collapseZeroes r+ collapseZeroes oth = error "dendrogramToGraph: internal error. Not expecting non-zero branch length here." ++-- | The plot looks nicer when the names aren't bloated with repeated stuff. This+-- replaces all tree names with potentially shorter names by removing the common prefix.+-- Returns how many characters were dropped.+truncateNames :: C.Dendrogram (FullTree a) -> (Int, C.Dendrogram (FullTree a))+truncateNames dendro = (prefChars, fmap chopName dendro)+ where + chopName ft = ft{ treename= drop prefChars (treename ft) }+ prefChars = length$ commonPrefix$ S.toList$ allNames dendro+ allNames (C.Leaf tr) = S.singleton (treename tr)+ allNames (C.Branch _ l r) = S.union (allNames l) (allNames r)++----------------------------------------------------------------------------------------------------++ -- | Open a GUI window to displaya tree. -- -- Fork a thread that then runs graphviz.@@ -121,7 +273,7 @@ Gv.graphToDot myparams graph where graph = toGraph2 atree- fromLabel ix = tbl ! ix + fromLabel ix = tbl M.! 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 }@@ -137,6 +289,36 @@ edgeAttrs = getEdgeAttrs edge_scale +-- | Some arbitrarily chosen colors from the X11 set:+defaultPaletteV :: V.Vector GA.Color+defaultPaletteV = V.fromList defaultPalette++defaultPalette :: [GA.Color]+defaultPalette = concat$ replicate 4 $ map GA.X11Color + [ Gv.Aquamarine+ , Gv.PaleVioletRed+ , Gv.MediumPurple+ , Gv.PaleGreen+ , Gv.PapayaWhip+ , Gv.SkyBlue+ , Gv.Yellow+ , Gv.Crimson+ , Gv.Gray+ , Gv.PaleGoldenrod+ ]++-- Grabbed the first couple palettes off a website:+altPalette :: V.Vector GA.Color+altPalette = V.fromList $ concat $ replicate 3 $ + -- http://www.colourlovers.com/palette/2962435/Autumn_Roses+ [ RGB 159 74 81, RGB 217 183 173, RGB 149 91 116, RGB 185 138 148+ -- , RGB 101 69 82 -- too dark+ ] +++ -- http://www.colourlovers.com/palette/2962434/Earthy_warm+ [ RGB 108 74 39, RGB 207 179 83, RGB 180 149 60, RGB 244 242 185+ -- , RGB 61 63 39+ ]+ getEdgeAttrs :: Double -> (t, t1, Double) -> [GA.Attribute] getEdgeAttrs edge_scale = edgeAttrs where @@ -152,8 +334,6 @@ 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 = @@ -161,42 +341,14 @@ -- 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-+weighted c = GC.WC {GC.wColor=c, GC.weighting=Nothing} -- | 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 + fromLabel ix = tbl M.! 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 }@@ -262,3 +414,13 @@ -- TEMP / HACK: prettyPrint' :: Show a => a -> String prettyPrint' = show++-- | Common prefix of a list of lists.+commonPrefix :: Eq a => [[a]] -> [a]+commonPrefix [] = []+commonPrefix ls@(hd:tl)+ | any null ls = []+ | otherwise =+ if all ((== (head hd)) . head) tl+ then head hd : commonPrefix (map tail ls)+ else commonPrefix (map tail ls)
Main.hs view
@@ -1,8 +1,8 @@-{-# LANGUAGE RecordWildCards, TupleSections #-}+{-# LANGUAGE RecordWildCards, TupleSections, NamedFieldPuns #-} {-# OPTIONS_GHC -fwarn-unused-imports -fwarn-incomplete-patterns #-} module Main where-import Data.List (sort)+import Data.List (sort, intersperse, foldl') import qualified Data.ByteString.Lazy.Char8 as B import qualified Data.Map as M import qualified Data.Set as S@@ -16,15 +16,19 @@ import Test.HUnit (runTestTT, Test, test, (~:)) import Control.Applicative ((<$>))+import qualified Data.Vector as V+import Test.HUnit as HU+ import Data.GraphViz (runGraphvizCanvas,GraphvizCommand(Dot),GraphvizCanvas(Xlib)) import Text.PrettyPrint.HughesPJClass hiding (char, Style) import Bio.Phylogeny.PhyBin.CoreTypes import Bio.Phylogeny.PhyBin (driver, binthem, normalize, annotateWLabLists,- unitTests, acquireTreeFiles, deAnnotate)+ unitTests, acquireTreeFiles, deAnnotate,+ retrieveHighlights, matchAnyHighlight) import Bio.Phylogeny.PhyBin.Parser (parseNewick, parseNewicks, parseNewickFiles, unitTests) import Bio.Phylogeny.PhyBin.Visualize (viewNewickTree, dotNewickTree_debug)-import Bio.Phylogeny.PhyBin.RFDistance (distanceMatrix, printDistMat, allBips)+import Bio.Phylogeny.PhyBin.RFDistance import Bio.Phylogeny.PhyBin.PreProcessor import qualified Data.Clustering.Hierarchical as C@@ -49,8 +53,13 @@ | View | TabDelimited Int Int + | Highlight FilePath+ | ShowTreesInDendro | ShowInterior++ | HashRF Bool | SelfTest- | RFMatrix | LineSetDiffMode | PrintNorms | PrintReg+ | RFMatrix | LineSetDiffMode+ | PrintNorms | PrintReg | PrintConsensus | PrintMatching | Cluster C.Linkage | BinningMode | EditDistThresh Int@@ -77,29 +86,43 @@ {- -- 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"+++ , Option ['t'] ["tabbed"] (ReqArg parseTabDelim "NUM1:NUM2")$ "assume the input is a ab-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 [] ["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 [] ["UPGMA"] (NoArg$ Cluster C.UPGMA) $ "Use Unweighted Pair Group Method (average linkage) - DEFAULT mode" , 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)"+ "Combine all clusters separated by DIST or less. Report a flat list of clusters.\n"+++ "Irrespective of whether this is activated, a hierarchical clustering (dendogram.pdf) is produced." +-- , Option [] ["dendogram"] (NoArg DendogramOnly)$ "Report a hierarchical clustering (default)"++ , Option [] [] (NoArg$ error "internal problem") " Select Robinson-Foulds (symmetric difference) distance algorithm:"+ , Option [] ["simple"] (NoArg$ HashRF False)+ ((if hashRF then "" else "(default) ")++ "use direct all-to-all comparison")+ , Option [] ["hashrf"] (NoArg$ HashRF True)+ ((if hashRF then "(default) " else "")++"use a variant of the HashRF algorithm for the distance matrix") , Option [] [] (NoArg NullOpt) "" , Option [] [] (NoArg$ error "internal problem") "----------------------------- Visualization --------------------------------"- , 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 ['g'] ["graphbins"] (NoArg Graph) "use graphviz to produce .dot and .pdf output files"+ -- TODO: Produce the consensus tree as well as the individual trees.+ , Option ['d'] ["drawbins"] (NoArg Draw) "like -g, but open GUI windows to show each bin's tree" , Option ['w'] ["view"] (NoArg View)$ "for convenience, \"view mode\" simply displays input Newick files without binning" + , Option [] ["showtrees"] (NoArg ShowTreesInDendro) "Print (textual) tree topology inside the nodes of the dendrogram"+ , Option [] ["highlight"] (ReqArg Highlight "FILE") $ + "Highlight nodes in the tree-of-trees (dendrogram) consistent with the.\n"+++ "given tree file. Multiple highlights are permitted and use different colors."+ , Option [] ["interior"] (NoArg ShowInterior)+ "Show the consensus trees for interior nodes in the dendogram, rather than just points."+ , Option [] [] (NoArg NullOpt) "" , Option [] [] (NoArg$ error "internal problem") "---------------------------- Tree pre-processing -----------------------------" @@ -141,7 +164,11 @@ , 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"+ , Option [] ["consensus"] (NoArg PrintConsensus) "print a strict consensus tree for the inputs, then exit"+ , Option [] ["matching"] (NoArg PrintMatching) "print a list of tree names that match any --highlight argument" ]+ where+ hashRF = use_hashrf default_phybin_config usage :: String usage = "\nUsage: phybin [OPTION...] files or directories...\n\n"++@@ -181,47 +208,7 @@ 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 @@ -241,6 +228,10 @@ RFMatrix -> return cfg { print_rfmatrix= True } + ShowTreesInDendro -> return cfg { show_trees_in_dendro = True }+ ShowInterior -> return cfg { show_interior_consensus = True }+ Highlight path -> return cfg { highlights = path : highlights cfg }+ LineSetDiffMode -> do bss <- mapM B.readFile files case map (S.fromList . B.lines) bss of@@ -251,10 +242,13 @@ oth -> error $"Line set difference mode expects two files as input, got "++show(length oth) exitSuccess - PrintNorms -> return cfg- PrintReg -> return cfg - + PrintNorms -> return cfg+ PrintReg -> return cfg+ PrintConsensus -> return cfg+ PrintMatching -> return cfg+ Cluster lnk -> return cfg { clust_mode = ClusterThem lnk }+ HashRF bl -> return cfg { use_hashrf = bl } BinningMode -> return cfg { clust_mode = BinThem } EditDistThresh n -> return cfg { dist_thresh = Just n } DendogramOnly -> return cfg { dist_thresh = Nothing }@@ -302,6 +296,25 @@ liftFT (normalize . annotateWLabLists) ft exitSuccess ------------------------------------------------------------+ when (elem PrintConsensus opts) $ do + (_,fts) <- parseNewickFiles (name_hack config) all_inputs+ putStrLn $ "Strict Consensus Tree of "++show (length fts)++" trees:"+ when (null fts) $ error "No trees provided!"+ let ctree = consensusTree (num_taxa config) (map nwtree fts)+ FullTree{labelTable} = head fts+ print$ displayStrippedTree$ FullTree "" labelTable ctree+ exitSuccess + ------------------------------------------------------------+ when (elem PrintMatching opts) $ do + (labs,fts) <- parseNewickFiles (name_hack config) all_inputs+ when (null$ highlights config) $ error "No --highlight given, so --matching makes no sense. Matching what?"+ htrs <- retrieveHighlights (name_hack config) labs (highlights config)+ let isMatch = matchAnyHighlight htrs+ forM_ fts $ \ FullTree{treename,nwtree} ->+ when (isMatch$ fmap (const()) nwtree) $+ putStrLn treename+ exitSuccess + ------------------------------------------------------------ when (View `elem` opts) $ do view_graphs config exitSuccess@@ -316,7 +329,7 @@ str <- B.readFile file putStrLn$ "Parsed: " ++ (B.unpack str) let (tbl,tr) = parseNewick M.empty name_hack file str- (chan, _tr) <- viewNewickTree file (FullTree "" tbl (annotateWLabLists tr))+ (chan, _tr) <- viewNewickTree file (FullTree "" tbl (annotateWLabLists (head tr))) return chan forM_ chans readChan return ()@@ -354,7 +367,50 @@ temp = driver default_phybin_config{ num_taxa=7, inputs=["../datasets/test.tr"] } +--------------------------------------------------------------------------------+-- Aggregated Unit Tests+--------------------------------------------------------------------------------++allUnitTests :: Test+-- allUnitTests = unitTests +++allUnitTests = test + [ Bio.Phylogeny.PhyBin.unitTests+ , Bio.Phylogeny.PhyBin.Parser.unitTests+ , "norm/Bip1" ~: (testNorm prob1)+ , "bipTreeConversion" ~: testBipConversion+ , "t3" ~: t3_consensusTest, "t4" ~: t4_consensusTest, "t5" ~: t5_consensusTest+ ]++-- [2013.07.23] +-- This was INCORRECTLY normalizing to:+-- ((1_, 2_), (7_, (18, 6_)), ((14, 3_), (19, (13, 5_))))+prob1 :: String+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,parseds) = parseNewick M.empty id "test" (B.pack str)+ parsed = head parseds+ 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 "+ ++concat (intersperse " " (map (dispBip labs) (S.toList added))) ++"\n "+ ++concat (intersperse " " (map (dispBip labs) (S.toList removed)))++ -- 112 and 13+rftest :: IO () rftest = do (mp,[t1,t2]) <- parseNewickFiles (take 2) ["tests/13.tr", "tests/112.tr"] putStrLn$ "Tree 13 : " ++ show (displayDefaultTree t1)@@ -369,7 +425,7 @@ putStrLn$ "13 collapsed 0.03: " ++show (disp$ liftFT (collapseBranchLenThresh 0.03) t1) putStrLn$ "112 collapsed 0.03: " ++show (disp$ liftFT (collapseBranchLenThresh 0.03) t2) - let mat = distanceMatrix [nwtree t1, nwtree t2]+ let (mat,_) = distanceMatrix [nwtree t1, nwtree t2] printDistMat stdout mat return () where@@ -377,6 +433,83 @@ let collapsed :: AnnotatedTree collapsed = normalize$ annotateWLabLists tr in displayDefaultTree$ deAnnotate $ FullTree nm labs collapsed++-- | This test was done with --editdist 4 --complete+t3_consensusTest :: IO ()+t3_consensusTest = consensusTest "./tests/t3_consensus/cluster1_284_alltrees.tr"+ "./tests/t3_consensus/cluster1_284_consensus.tr"++-- | This test was done with --editdist 0 --complete+t4_consensusTest :: IO ()+t4_consensusTest = consensusTest "./tests/t4_consensus/cluster1_16_alltrees.tr"+ "./tests/t4_consensus/cluster1_16_consensus.tr"++-- | This test was done with --editdist 1 --complete+t5_consensusTest :: IO ()+t5_consensusTest = consensusTest "./tests/t5_consensus/cluster1_35_alltrees.tr"+ "./tests/t5_consensus/cluster1_35_consensus.tr"++consensusTest :: String -> String -> IO ()+consensusTest alltrees consensus = do + (_,ctree:ftrees) <- parseNewickFiles id [consensus,alltrees]+ let num_taxa = numLeaves (nwtree ctree)+ plainTrs = map nwtree ftrees + eachbips = map allBips plainTrs+ totalBips = foldl' S.union S.empty eachbips+ intersectBips = foldl' S.intersection S.empty eachbips+ FullTree _ labs _ = ctree+ linesPrnt x = unlines (map ((" "++) . dispBip labs) $ S.toList x)+ putStrLn$ "Bips in each: "++ show (map S.size eachbips)+ putStrLn$ "Total bips in all: "++show (S.size totalBips) + putStrLn$ "Bips in common: "++ show (S.size intersectBips)+ putStrLn$ "Bips of first member:\n" ++ linesPrnt (head eachbips)+ putStrLn$ "Some bips in the union that are NOT in the intersection:\n" +++ linesPrnt (S.fromList$ take 20$ S.toList$ S.difference totalBips intersectBips)+ let cbips = allBips $ nwtree ctree+ putStrLn$ "ConsensusBips ("++show (S.size cbips)++"):\n"++linesPrnt cbips+ putStrLn$"Things in the consensus that should NOT be:\n"++linesPrnt (S.difference cbips intersectBips)+ putStrLn$"Things not in the consensus that SHOULD be:\n"++linesPrnt (S.difference intersectBips cbips)++ putStrLn$ "Now recomputing consensus tree for "++show num_taxa++" taxa"+ let ctree2 = consensusTree num_taxa plainTrs+ cbips2 = allBips ctree2+ putStrLn$ "Freshly recomputed consensusBips ("++show (S.size cbips2)++"):\n"++linesPrnt cbips2+ HU.assertEqual "Consensus tree on disk should match computed one:"+ cbips cbips2 -- (allBips$ fmap (const ()) $ nwtree ctree) + + putStrLn " Partial distance matrix WITHIN this cluster:"+ let (mat,_) = distanceMatrix (map nwtree ftrees)+ printDistMat stdout (V.take 30 mat)+ HU.assertBool "Consensus should only include bicbips2ps in the members" (S.isSubsetOf cbips totalBips)+ HU.assertEqual "Consensus tree matches intersected bips" cbips intersectBips + return ()++testBipConversion :: IO ()+testBipConversion = + do (_,trs) <- allTestTrees+ mapM_ testOne trs+ putStrLn "All trees passed bip conversion."+ where+ testOne (FullTree{nwtree}) = do+ let sz = numLeaves nwtree + bips = allBips nwtree+ tr2 = bipsToTree sz bips+ bips2 = allBips tr2+ assertEqual "round trip bips->tree->bips" bips bips2++-- | Read in all test trees which we happen to have put in the repository for testing+-- purposes.+allTestTrees :: IO (LabelTable, [FullTree DefDecor])+allTestTrees =+ parseNewickFiles id $+ [ "./tests/t3_consensus/cluster1_284_alltrees.tr"+ , "./tests/t3_consensus/cluster1_284_consensus.tr"+ , "./tests/t4_consensus/cluster1_16_alltrees.tr"+ , "./tests/t4_consensus/cluster1_16_consensus.tr"+ , "./tests/t5_consensus/cluster1_35_alltrees.tr"+ , "./tests/t5_consensus/cluster1_35_consensus.tr"+ ]+ ---------------------------------------------------------------------------------------------------- -- TODO: expose a command line argument for testing.
phybin.cabal view
@@ -1,5 +1,5 @@ Name: phybin-Version: 0.2.2+Version: 0.2.11 License: BSD3 License-file: LICENSE Stability: Beta@@ -14,9 +14,17 @@ -- 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+-- 0.2.4 -- add consensus trees+-- 0.2.6 -- colorization, hashrf, misc improvements+-- 0.2.7 -- Add command line opt --showtrees+-- 0.2.8 -- Add command line opt --highlight+-- 0.2.9 -- Add command line opt --interior+-- 0.2.10 -- Add command line opt --matching+-- 0.2.11 -- Cleanup and windows compatibility. -- homepage: http://code.haskell.org/phybin-homepage: http://people.csail.mit.edu/newton/phybin/+-- homepage: http://people.csail.mit.edu/newton/phybin/+homepage: http://www.cs.indiana.edu/~rrnewton/projects/phybin/ Copyright: Copyright (c) 2010 Ryan Newton Synopsis: Utility for clustering phylogenetic trees in Newick format based on Robinson-Foulds distance.@@ -26,7 +34,7 @@ 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+ phybin produces 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@@ -40,6 +48,18 @@ type: git location: git://github.com/rrnewton/PhyBin.git +Flag hashrf+ description: Use the HashRF algorithm by default instead of the naive one.+ default: True++Flag bitvec+ description: Use bitvectors rather than IntSets for bipartitions.+ default: False++Flag sequential+ description: Don't use any parallelism at all.+ default: False+ Library Exposed-modules: Bio.Phylogeny.PhyBin Bio.Phylogeny.PhyBin.RFDistance@@ -49,34 +69,50 @@ Bio.Phylogeny.PhyBin.Util Bio.Phylogeny.PhyBin.Visualize -- Other-modules: Data.BitList- Build-Depends: base >= 3 && < 5, directory, process, containers, unix, + Build-Depends: base >= 3 && < 5, directory, process, containers, async, time, filepath, graphviz >= 2999.16, text >= 0.11 && < 0.12, prettyclass, fgl,- HSH, HUnit, bytestring, + HUnit, bytestring, -- For bytestring.lazy support: parsec >= 3.1.0, - bitvec >= 0.1, vector >= 0.10,- hierarchical-clustering >= 0.4+ vector >= 0.10,+ hierarchical-clustering >= 0.4, split >= 0.2 -- lattice-par, GHC-Options: -O2 -funbox-strict-fields -rtsopts + if flag(hashrf)+ CPP-Options: -DUSE_HASHRF+ if flag(bitvec)+ CPP-Options: -DBITVEC_BIPS+ Build-Depends: bitvec >= 0.1+ if flag(sequential)+ CPP-Options: -DSEQUENTIALIZE+ Executable phybin Main-is: Main.hs Build-Depends: phybin -- DUPLICATE:- Build-Depends: base >= 3 && < 5, directory, process, containers, unix, + Build-Depends: base >= 3 && < 5, directory, process, containers, async, time, filepath, graphviz >= 2999.16, text >= 0.11 && < 0.12, prettyclass, fgl,- HSH, HUnit, bytestring, + HUnit, bytestring, -- For bytestring.lazy support: parsec >= 3.1.0, - bitvec >= 0.1, vector >= 0.10,- hierarchical-clustering >= 0.4+ vector >= 0.10,+ hierarchical-clustering >= 0.4, split >= 0.2 -- lattice-par, GHC-Options: -O2 -funbox-strict-fields -rtsopts++ if flag(hashrf)+ CPP-Options: -DUSE_HASHRF+ if flag(bitvec)+ CPP-Options: -DBITVEC_BIPS+ Build-Depends: bitvec >= 0.1+ if flag(sequential)+ CPP-Options: -DSEQUENTIALIZE