phybin 0.1.2.1 → 0.1.2.4
raw patch · 3 files changed
+315/−1207 lines, 3 filesdep +textdep ~graphviz
Dependencies added: text
Dependency ranges changed: graphviz
Files
- Bio/Phylogeny/PhyBin/Main.hs +0/−1203
- Main.hs +300/−0
- phybin.cabal +15/−4
− Bio/Phylogeny/PhyBin/Main.hs
@@ -1,1203 +0,0 @@-{-# LANGUAGE ScopedTypeVariables, RecordWildCards, TypeSynonymInstances, CPP #-}-{-# OPTIONS_GHC -fwarn-incomplete-patterns #-}----module Bio.Phylogeny.PhyBin.Main where-module Main where--import Text.Printf-import Text.Parsec-import Text.Parsec.Char-import Text.Parsec.Combinator-import Text.Parsec.ByteString.Lazy --import Numeric-import Data.Function-import Data.List-import Data.Maybe-import Data.Char-import qualified Data.ByteString.Lazy.Char8 as B-import qualified Data.Map as M-import qualified Data.Set as S--import Control.Monad-import Control.Exception hiding (try)---import Control.DeepSeq-import Control.Concurrent---- Wow, I actually couldn't figure out how to open a file (and get a--- HANDLE) so that I could then use getFileAttributes under--- System.Win32. Giving up because I think I can just use the--- OS-independent System.Directory--import System.Directory-import System.FilePath--import System.Environment-import System.Directory-import System.Console.GetOpt-import System.Exit-import System.IO--import Test.HUnit-import HSH hiding (run)-import qualified HSH ---- For vizualization:-import Data.Graph.Inductive as G hiding (run) -import Data.GraphViz as Gv hiding (Label, toLabel) -import qualified Data.Graph.Inductive as G -import qualified Data.GraphViz as Gv ---import Data.Graph.Inductive.Query.DFS--import Text.PrettyPrint.HughesPJClass hiding (char, Style)--import Debug.Trace--phybin_version = "0.1.2.1" -- NOTE: Remember to keep me in sync with the .cabal file--------------------------------------------------------------------------------------------------------- Type definitions-------------------------------------------------------------------------------------------------------type BranchLen = Double---- Even though the Newick format allows it, ignoring interior node--- labels. (They are not commonly used.)-data NewickTree a = - NTLeaf a Label- | NTInterior a [NewickTree a]- deriving (Show, Eq, Ord)--{---- [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 Pretty (NewickTree dec) where - pPrint (NTLeaf dec name) = text (fromLabel name)- pPrint (NTInterior dec ls) = - --parens$ commasep ls- (parens$ sep$ map_but_last (<>text",") $ map pPrint ls)----- Experimental: toggle this to change the representation of labels:--------------------------------------------type Label = Atom; (toLabel, fromLabel) = (toAtom, fromAtom)------------------------------------------type Label = String; (toLabel, fromLabel) = (id, id)------------------------------------------fromLabel :: Label -> String--------------------------------------------------------------------------------------------------------- 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 = doesFileExist-is_directory = doesDirectoryExist -file_exists path = - do f <- doesFileExist path- d <- doesDirectoryExist path- return (f || d)--------------------------------------------------------------------------------------------------------- General helper/utility functions:---------------------------------------------------------------------------------------------------------commacat ls = hcat (intersperse (text ", ") $ map pPrint ls)-commasep ls = sep (intersperse (text ", ") $ map pPrint ls)--map_but_last fn [] = []-map_but_last fn [h] = [h]-map_but_last fn (h:t) = fn h : map_but_last fn t--fst3 (a,_,_) = a-snd3 (_,b,_) = b-thd3 (_,_,c) = c--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 ls [] = ls-demerge [] ls = error$ "demerge: first list did not contain all of second, remaining: " ++ show ls-demerge l@(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 Nothing ls = ls-maybeCons (Just x) ls = x : ls--maybeInsert fn Nothing ls = ls-maybeInsert fn (Just x) ls = insertBy fn x ls---maybeInsert Nothing ls = ls---maybeInsert (Just x) ls = insert (x) ls -- Relies on ORD--------------------------------------------------------------------------------------------------------- Newick file format parser definitions:-------------------------------------------------------------------------------------------------------tag l s =- case s of - NTLeaf _ n -> NTLeaf l n- NTInterior _ ls -> NTInterior l ls---- This parser ASSUMES that whitespace has been prefiltered from the input.-newick_parser :: Parser (NewickTree BranchLen)-newick_parser = - do x <- subtree- l<-len- char ';'- return$ tag l x--subtree :: Parser (NewickTree BranchLen)-subtree = internal <|> leaf--leaf :: Parser (NewickTree BranchLen)-leaf = do n<-name; return$ NTLeaf 0.0 (toLabel n)--internal :: Parser (NewickTree BranchLen)-internal = do char '(' - bs <- branchset- char ')' - nm <- name -- IGNORED- return$ NTInterior 0.0 bs--branchset :: Parser [NewickTree BranchLen]-branchset =- do b <- branch <?> "at least one branch"- rest <- option [] $ try$ do char ','; branchset- return (b:rest)--branch :: Parser (NewickTree BranchLen)-branch = do s<-subtree; l<-len; - return$ tag l s--len :: Parser Double-len = option 0.0 $ do char ':'; number--number :: Parser Double-number = - do sign <- option "" $ string "-"- fst <- many1 digit- snd <- option "0" $ try$ do char '.'; many1 digit- return (read (sign ++ fst++"."++snd) :: Double)--name :: Parser String-name = option "" $ many1 (letter <|> digit <|> oneOf "_.-")---------------------------------------------------------------------------------------------------------- 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.-------------------------------------------------------------------------------------- A common type of tree is "AnnotatedTree", which contains three things:--- (1) branch length from parent to "this" node--- (2) subtree weights for future use--- (defined as number of LEAVES, not counting intermediate nodes)--- (3) sorted lists of labels for symmetry breaking-type AnnotatedTree = NewickTree (BranchLen, Int, [Label])---annotateWLabLists :: NewickTree BranchLen -> AnnotatedTree-annotateWLabLists tr = case tr of - NTLeaf bl n -> NTLeaf (bl,1,[n]) n- NTInterior bl ls -> - let children = map annotateWLabLists ls in - NTInterior (bl, sum $ map get_weight children,- foldl1 merge $ map get_label_list children)- children--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)--map_labels fn (NTLeaf dec lab) = NTLeaf dec $ fn lab-map_labels fn (NTInterior dec ls) = NTInterior dec$ map (map_labels fn) ls--all_labels (NTLeaf _ lab) = [lab]-all_labels (NTInterior _ ls) = concat$ map all_labels ls--get_dec (NTLeaf dec _) = dec-get_dec (NTInterior dec _) = dec---- Set all the decorations to a constant:-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 (NTLeaf _ _) = []-get_children (NTInterior _ ls) = ls---- Number of LEAVES contained in subtree:-get_weight = snd3 . get_dec---- Sorted list of leaf labels contained in subtree-get_label_list = thd3 . get_dec--add_weight (l1,w1,sorted1) node = - let (_,w2,sorted2) = get_dec node in - (l1, w1+w2, merge sorted1 sorted2)---- Remove the influence of one subtree from the metadata of another.-subtract_weight (l1,w1,sorted1) node = - let (_,w2,sorted2) = get_dec node in - (l1, w1-w2, demerge sorted1 sorted2)---- Turn on for extra invariant checking:-debug = False- ------------------------------------------------------------------------------------- I ran into a nasty bug as a result of "deriving Ord". But I didn't end up getting rid of it.---instance Ord AnnotatedTree where --- compare (NTLeaf _ _) (NTInterior _ _) = LT--- compare (NTLeaf _ _) (NTLeaf _ _) = EQ---compare_nodes :: AnnotatedTree -> AnnotatedTree -> Ordering--- Our sorting criteria for the children of interior nodes:--compare_childtrees node1 node2 = - case get_weight node1 `compare` get_weight 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 get_label_list node1 `compare` get_label_list node2 of- --EQ -> error "FIXME EQ"- 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----- 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.- -- 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 context = case node of- NTLeaf dec@(l,w,sorted) name -> - (node, - -- If the leaf becomes the root... we could introduce another node:- NTInterior (add_context (0,w,sorted) context) $- (verify_sorted "1" id$ maybeInsert compare_childtrees context [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@(l,w,_) 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 context- inverted_children = verify_sorted "2" id$ maybeInsert compare_childtrees context 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 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 = normalize $ annotateWLabLists $ run newick_parser "(A,(C,D,E),B);"--tt0 = drawNewickTree "tt0" $ annotateWLabLists $ run newick_parser "(A,(C,D,E),B);"--tt2 = toGraph tt-tt3 = drawNewickTree "tt3" tt--norm4 = norm "((C,D,E),B,A);"-tt4 = drawNewickTree "tt4"$ trace ("FINAL: "++ show (pPrint norm4)) $ norm4--norm5 = normalize$ annotateWLabLists$ run newick_parser "(D,E,C,(B,A));"-tt5 = drawNewickTree "tt5"$ norm5--tt5' = prettyPrint' $ dotNewickTree "norm5" 1.0 norm5--ttall = do tt3; tt4; tt5--------------------------------------------------------------------------------------------------------- Equivalence classes:-------------------------------------------------------------------------------------------------------data BinEntry = BE {- members :: [String], - trees :: [AnnotatedTree]-}- deriving Show ---- We index the results of binning by topology-only trees that have their decorations removed.--- (But we leave the weights on and leave the type as AnnotatedTree so as to acces Ord.)-type BinResults = M.Map AnnotatedTree BinEntry---- Takes labeled trees, classifies labels into equivalence classes.---binthem :: [(String, NewickTree BranchLen)] -> M.Map (NewickTree ()) BinEntry-binthem :: [(String, NewickTree BranchLen)] -> BinResults-binthem ls = binthem_normed normalized- where- normalized = map (\ (lab,tree) -> (lab, normalize $ annotateWLabLists tree)) ls---- This version accepts trees that are already normalized:-binthem_normed :: [(String, AnnotatedTree)] -> BinResults-binthem_normed normalized = --- foldl (\ acc (lab,tree) -> M.insertWith update tree (BE{ members=[lab] }) acc)- foldl (\ acc (lab,tree) -> M.insertWith update (anonymize_annotated tree) (BE [lab] [tree]) acc)- M.empty normalized- --(map (mapSnd$ fmap (const ())) normalized) -- still need to STRIP them- where - --(++)--- update new old = BE{ members= (members new ++ members old) }- update new old = BE (members new ++ members old) (trees new ++ trees old)- --strip = fmap (const ())---- Remove branch lengths and labels but leave weights-anonymize_annotated :: AnnotatedTree -> AnnotatedTree-anonymize_annotated = fmap (\ (bl, w, labs) -> (0, 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);--}---------------------------------------------------------------------------------------------------------- Visualization with GraphViz and FGL:--------------------------------------------------------------------------------------------------------- First we need to be able to convert to FGL graphs:-toGraph :: AnnotatedTree -> Gr String Double-toGraph tree = run_ G.empty $ loop tree- where- loop (NTLeaf (l,_,_) name) = - do let str = fromLabel name- G.insMapNodeM str- return str- loop (NTInterior (l,_,sorted) ls) =- do let bigname = concat$ map fromLabel sorted- names <- mapM loop ls- G.insMapNodeM bigname- mapM_ (\x -> insMapEdgeM (bigname, x, 0.0)) names- return bigname---- This version uses the tree nodes themselves as labels.-toGraph2 :: AnnotatedTree -> Gr AnnotatedTree Double-toGraph2 tree = run_ G.empty $ loop tree- where- loop node@(NTLeaf _ _) = - do G.insMapNodeM node - return ()- loop node@(NTInterior (l,_,sorted) ls) =- do mapM_ loop ls- G.insMapNodeM node - mapM_ (\x -> insMapEdgeM (node, x, fst3$ get_dec x)) ls- return ()----- The channel retuned will carry a single message to signal--- completion of the subprocess.-drawNewickTree :: String -> AnnotatedTree -> IO (Chan (), AnnotatedTree)-drawNewickTree title tree =- do chan <- newChan- let dot = dotNewickTree title (1.0 / avg_branchlen [tree])- tree- runit = do runGraphvizCanvas default_cmd dot 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 = Neato---- Show a float without scientific notation:-myShowFloat :: Double -> String--- showFloat weight = showEFloat (Just 2) weight ""-myShowFloat fl = printf "%.4f" fl---dotNewickTree :: String -> Double -> AnnotatedTree -> DotGraph G.Node-dotNewickTree title edge_scale tree = - --trace ("EDGE SCALE: " ++ show edge_scale) $- graphToDot myparams graph- where - graph = toGraph2 tree- myparams :: GraphvizParams AnnotatedTree Double () AnnotatedTree- myparams = defaultParams { globalAttributes= [GraphAttrs [Gv.Label (StrLabel title)]],- fmtNode= nodeAttrs, fmtEdge= edgeAttrs }- nodeAttrs :: (Int,AnnotatedTree) -> [Attribute]- nodeAttrs (num,node) =- let children = get_children node in - [ Gv.Label (StrLabel$ concat$ map fromLabel$ thd3$ get_dec node)- , Shape (if null children then {-PlainText-} Ellipse else PointShape)- , Style [SItem Filled []]- ]-- -- TOGGLE:- -- edgeAttrs (_,_,weight) = [ArrowHead noArrow, Len (weight * edge_scale + bump), Gv.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) $- [ArrowHead noArrow, Gv.Label (StrLabel$ myShowFloat weight)] ++ -- TEMPTOGGLE- --[ArrowHead noArrow, Gv.Label (StrLabel$ show draw_weight)] ++ -- TEMPTOGGLE- if weight == 0.0- then [Color [X11Color Red], Len minlen]- else [Len draw_weight]- 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)----- This version shows the ordered/rooted structure of the normalized tree.-dotNewickTree_debug :: String -> AnnotatedTree -> DotGraph G.Node-dotNewickTree_debug title tree = graphToDot myparams graph- where - graph = toGraph2 tree- myparams :: GraphvizParams AnnotatedTree Double () AnnotatedTree- myparams = defaultParams { globalAttributes= [GraphAttrs [Gv.Label (StrLabel title)]],- fmtNode= nodeAttrs, fmtEdge= edgeAttrs }- nodeAttrs :: (Int,AnnotatedTree) -> [Attribute]- nodeAttrs (num,node) =- let children = get_children node in - [ Gv.Label (if null children - then StrLabel$ concat$ map fromLabel$ thd3$ get_dec node- else RecordLabel$ take (length children) $ - -- This will leave interior nodes unlabeled:- map (PortName . PN) $ 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)- )- , Shape Record- , Style [SItem Filled []]- ]-- edgeAttrs (num1,num2,weight) = - let node1 = fromJust$ lab graph num1 - node2 = fromJust$ lab graph num2 - ind = fromJust$ elemIndex node2 (get_children node1)- in [TailPort$ LabelledPort (PN$ show$ 1+ind) (Just South)]----------------------------------------------------------------------------------------------------------- Utilities and UNIT TESTING-------------------------------------------------------------------------------------------------------parseNewick :: String -> B.ByteString -> NewickTree BranchLen-parseNewick file input = - runB file newick_parser $- B.filter (not . isSpace) input--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--runPr prs str = print (run prs str)-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"--cnt :: NewickTree a -> Int-cnt (NTLeaf _ _) = 1-cnt (NTInterior _ ls) = 1 + sum (map cnt ls)--tr1 = run newick_parser "(A:0.1,B:0.2,(C:0.3,D:0.4):0.5);"-tr1draw = drawNewickTree "tr1"$ annotateWLabLists tr1-tr1dot = putStrLn$ prettyPrint' $ dotNewickTree "" 1.0 $ annotateWLabLists tr1---norm = normalize . annotateWLabLists . run newick_parser-norm2 = normalize . annotateWLabLists . parseNewick "test"-tests = - let - ntl s = NTLeaf 0.0 (toLabel s)- in - 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 "A"- , "subtree" ~: ntl "A" ~=? run subtree "A"-- -- These are not allowed:- , "null branchset" ~: errortest$ run branchset ""-- , "internal" ~: NTInterior 0.0 [ntl "A"] ~=? run internal "(A);"-- , "example: no nodes are named" ~: NTInterior 0 [ntl "", ntl "",NTInterior 0 [ntl "", ntl ""]]- ~=? run newick_parser "(,,(,));"- , "example: leaf nodes are named" ~: 6 ~=? cnt (run newick_parser "(A,B,(C,D));")- , "example: all nodes are named" ~: 6 ~=? cnt (run newick_parser "(A,B,(C,D)E)F;")-- , "example: all but root node have a distance to parent" ~: 6 ~=? cnt (run newick_parser "(:0.1,:0.2,(:0.3,:0.4):0.5);")- , "example: all have a distance to parent" ~: 6 ~=? cnt (run newick_parser "(:0.1,:0.2,(:0.3,:0.4):0.5):0.6;")- , "example: distances and leaf names (popular)" ~: 6 ~=? cnt tr1- , "example: distances and all names" ~: 6 ~=? cnt (run newick_parser "(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 ~=? cnt (run newick_parser "((B:0.2,(C:0.3,D:0.4)E:0.5)F:0.1)A;")-- , "merge" ~: [1,2,3,4,5,6] ~=? merge [1,3,5] [2,4,6]-- , "demerge" ~: [2,4,6] ~=? demerge [1,2,3,4,5,6] [1,3,5]- , "demerge" ~: [1,3,5] ~=? demerge [1,2,3,4,5,6] [2,4,6]-- , "annotateWLabLists" ~: - --NTInterior (0.0,[A,B,C,D]) [NTLeaf (0.1,[A]) A,NTLeaf (0.2,[B]) B,NTInterior (0.5,[C,D]) [NTLeaf (0.3,[C]) C,NTLeaf (0.4,[D]) D]]- map toLabel ["A","B","C","D"] -- ORD on atoms is expensive... it must use the whole string.- ~=? get_label_list (annotateWLabLists tr1)-- -- Make sure that all of these normalize to the same thing.- , "normalize1" ~: "(C, D, E, (A, B))" ~=? show (pPrint$ 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);"--- , "dotConversion" ~: True ~=? 100 < length (prettyPrint' $ dotNewickTree "" 1.0$ norm "(D,E,C,(B,A));") -- 444-- - , "phbin: these 3 trees should fall in the same category" ~: - 1 ~=? (length $ M.toList $- binthem [("one", run newick_parser "(A,(C,D,E),B);"),- ("two", run newick_parser "((C,D,E),B,A);"),- ("three", run newick_parser "(D,E,C,(B,A));")])-- ]--run_tests = runTestTT tests-t = run_tests-- -------------------------------------------------------------------------------------------------------- Driver to put the pieces together (parse, normalize, bin)--------------------------------------------------------------------------------------------------------- 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 :: Label -> Label- , output_dir :: String- , inputs :: [String]- , do_graph :: Bool- , do_draw :: Bool- }--default_phybin_config = - PBC { verbose = False- , num_taxa = error "must be able to determine the number of taxa expected in the dataset"- , name_hack = id -- Default, no transformation of leaf-labels- , output_dir = "./"- , inputs = []- , do_graph = False- , do_draw = False- }---driver :: PhyBinConfig -> IO ()-driver PBC{..} =- do - --------------------------------------------------------------------------------- -- First, find out where we are and open the files:- --------------------------------------------------------------------------------- cd <- getCurrentDirectory - --putStrLn$ "PHYBIN RUNNING IN DIRECTORY: "++ cd-- 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-- let files = concat all -- take 10 $ concat all- num_files = length files-- putStrLn$ "Parsing "++show num_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.- --------------------------------------------------------------------------------- --- -- results contains: num-nodes, parsed, warning-files- results <- forM files $ \ file -> - do --stat <- getFileStatus file - reg <- is_regular_file file- if not reg then return (0,[],[(-1, file)]) else do-- h <- openFile file ReadMode - bstr <- B.hGetContents h-- -- Clip off the first three characters:- let - parsed = map_labels name_hack $ parseNewick file bstr- annot = annotateWLabLists parsed- normal = normalize annot- weight = get_weight annot-- -- TEMPTOGGLE- when False $ do putStrLn$ "DRAWING TREE"; drawNewickTree "Annotated" annot; drawNewickTree "Normalized" normal- putStrLn$ "WEIGHTS OF NORMALIZED' CHILDREN: "++ show (map get_weight$ get_children normal)-- if not$ weight == num_taxa- then do --putStrLn$ "\n WARNING: file contained an empty or single-node tree: "++ show file- when verbose$ putStrLn$ "\n WARNING: file contained unexpected number of leaves ("- ++ show weight ++"): "++ show file- return (0,[], [(weight, file)])- else do - when verbose$ putStr "."-- --evaluate$ deepseq$ runB newick_parser bstr- --evaluate$ cnt$ runB newick_parser bstr- num <- evaluate$ cnt parsed- --num <- evaluate$ cnt normal-- hClose h- --return$ (num, [normal])- return$ (num, [parsed], [])-- putStrLn$ "\nNumber of input trees: " ++ show num_files- putStrLn$ "Number of VALID trees (correct # of leaves/taxa): " ++ show (length$ concat$ map snd3 results)- putStrLn$ "Total tree nodes contained in valid trees: "++ show (sum$ map fst3 results)-- --------------------------------------------------------------------------------- -- Do the actual binning:- ---------------------------------------------------------------------------------- putStrLn$ "Creating equivalence classes (bins), bin sizes:"-- let classes = --binthem_normed$ zip files $ concat$ map snd3 results- binthem$ zip files $ concat$ map snd3 results- binlist = reverse $ sortBy (compare `on` fst3) $- map (\ (tr,ls) -> (length (members ls), tr,ls)) $ M.toList classes- numbins = length binlist- taxa = S.unions$ map (S.fromList . all_labels . snd3) binlist- warnings = concat $ map thd3 results- base i size = combine output_dir ("bin" ++ show i ++"_"++ show size)-- ----------------------------------------- -- 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.- ---------------------------------------------------------------------------------- forM_ binlist $ \ (len, tr, _) -> do- when (len > 1) $ -- Omit that long tail of single element classes...- putStrLn$ " "++ show (pPrint tr) ++" members: "++ show len-- putStrLn$ "\nTotal unique taxa ("++ show (S.size taxa) ++"): "++ - show (sep $ map (text . fromLabel) $ S.toList taxa)-- putStrLn$ "Final number of tree bins: "++ show (M.size classes)- forM_ (zip [1..] binlist) $ \ (i, (size, tr, bentry)) -> do- --putStrLn$ (" WRITING " ++ combine output_dir ("bin" ++ show i ++"_"++ show size ++".txt"))- writeFile (base i size ++".txt") (concat$ map (++"\n") (members bentry))- writeFile (base i size ++".tr") (show (pPrint tr) ++ ";\n")--- writeFile (base i size ++".rawtree") (show tr ++ ";\n") -- TempToggle-- when (not$ null warnings) $- writeFile (combine output_dir "bin_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 (\ (n,file) -> - (if n == -1 - then "Not a regular/readable file: "++ file - else "Wrong number of taxa ("++ show n ++"): "++ file)- ++"\n") - warnings))- putStrLn$ "Wrote contents of each bin to bin<N>_<binsize>.txt"- putStrLn$ "Wrote representative trees to bin<N>_<binsize>.tr"-- when (do_graph) $ do- forM_ (zip [1..] binlist) $ \ (i, (size, tr, bentry)) -> do- let - dot = dotNewickTree ("bin #"++ show i) (1.0 / avg_branchlen (trees bentry))- --(annotateWLabLists$ fmap (const 0) tr)- -- TEMP FIXME -- using just ONE representative tree:- (--trace ("WEIGHTED: "++ show (head$ trees bentry)) $ - --(head$ trees bentry) )- (avg_trees$ trees bentry))- when (size > 1 || numbins < 100) $ do - runGraphvizCommand default_cmd dot Pdf (base i size ++ ".pdf")- return ()- putStrLn$ "Wrote visual representations of trees to bin<N>_<binsize>.pdf"-- --putStrLn$ "Wrote representative tree to bin<N>_<binsize>.tr"- putStrLn$ "Finished."- --------------------------------------------------------------------------------- -- End driver- ------------------------------------------------------------------------------------- Average branch length across several trees.-avg_branchlen :: [AnnotatedTree] -> Double-avg_branchlen ls = fst total / snd total- where- total = sum_ls $ map sum_tree ls- sum_ls ls = (sum$ map fst ls, sum$ map snd ls)-{-- sum_tree (NTLeaf (l,_,_) _) | l < 0 = - trace ("!!! GOT NEGATIVE BRANCH LENGTH: "++ show l) $- (0,0)--}- sum_tree (NTLeaf (0,_,_) _) = (0,0)- sum_tree (NTLeaf (l,_,_) _) = (abs l,1)- sum_tree (NTInterior (l,_,_) ls) = - let (x,y) = sum_ls$ map sum_tree ls in- if l == 0 then (x, y) else ((abs l) + x, 1+y)--{--nonzero_blens :: AnnotatedTree -> Int-nonzero_blens node = - let children = sum $ map nonzero_blens $ get_children node in- if (fst3 $ get_dec node) == 0 - then children- else children + 1--}---- 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 ls = --summed -- TEMPTOGGLE- fmap (\ (blen,w,ls) -> (blen / count, w,ls)) summed- where- summed = foldl1 sum_2trees ls- count = fromIntegral$ length ls-- sum_2trees a b = case (a,b) of- (NTLeaf (l1,w,ls) nm, NTLeaf (l2,_,_) _) -> NTLeaf (l1+l2,w,ls) nm- (NTInterior (l1,w,ls) ls1, - NTInterior (l2,_,_) ls2) -> NTInterior (l1+l2,w,ls) $ - map (uncurry sum_2trees) $ zip ls1 ls2- _ -> error "avg_trees: applied to non-isomorphic trees"---bump = 0.00001 -- for DIRTY HACKS--{- - ----------------------------------------- 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.-- -}---------------------------------------------------------------------------------------------------------- MAIN script: Read command line options and call the program.--------------------------------------------------------------------------------------------------------- Note: ORDER is important here, we process options in this order:-data Flag - = Verbose - | Version - | Output String- | NumTaxa Int- | NullOpt- | Graph | Draw- | Force - | View- | TabDelimited Int Int-- | SelfTest-- | NameCutoff String- | NamePrefix Int- | NameTable String -- Must come after Cutoff/Prefix-- deriving (Show, Eq, Ord) -- ORD is really used.--parseTabDelim str = - TabDelimited 8 9- -options :: [OptDescr Flag]-options =- [ 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 [] ["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") "----------------------------- Visualization --------------------------------"- , Option ['g'] ["graphbins"] (NoArg Graph) "use graphviz to produce .dot and .pdf output files named bin1.*, bin2.*, etc"- , 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!!- , Option ['n'] ["numtaxa"] (ReqArg (NumTaxa . read) "NUM") "expect NUM taxa for this dataset"--{- -- 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."-- , Option [] [] (NoArg NullOpt) ""- , Option ['s'] ["namesep"] (ReqArg NameCutoff "STR") $ --"\n"++- "An alternative to --nameprefix, STR provides a set of delimeter characters,\n"++- "for example '-' or '0123456789'. The taxa name is then a variable-length\n"++- "prefix of each gene name up to but not including any character in STR."-- , Option [] [] (NoArg NullOpt) ""- , Option ['m'] ["namemap"] (ReqArg NameTable "FILE") $ - "Even once prefixes are extracted it may be necessary to use a lookup table\n"++- "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."- ]--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"++-- "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"++-- "\nOptions include:\n"-defaultErr errs = error $ "ERROR!\n" ++ (concat errs ++ usageInfo usage options)---main = - do argv <- getArgs -- (opts,files) <- - case getOpt Permute options argv of- (o,n,[] ) -> return (o,n)- (_,_,errs) -> defaultErr errs-- 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 run_tests; exitSuccess-- Output s -> return cfg { output_dir= s }-- NumTaxa n -> return cfg { num_taxa= n }- Graph -> return cfg { do_graph= True } - Draw -> return cfg { do_draw = True } - View -> return cfg -- Handled below-- TabDelimited _ _ -> error "tabbed option not handled yet"- 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 }-- -- 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 } - (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--view_graphs :: PhyBinConfig -> IO ()-view_graphs PBC{..} = - do chans <- forM inputs $ \ file -> do - putStrLn$ "Drawing "++ file ++"...\n"- str <- B.readFile file- putStrLn$ "Parsed: " ++ (B.unpack str)- (chan, tr) <- drawNewickTree file $ - annotateWLabLists$ - map_labels name_hack $ - parseNewick file str- return chan- forM_ chans readChan - return ()-------------------------------------------------------------------------------------- Every dataset it seems needs a new hack on the names.-name_table_reader file = - do contents <- readFile file- let mp = M.fromList $ - map (\ls -> case ls of - [a,b] -> (a,b)- _ -> error$ "Each line of "++file++"must contain two whitespace free strings: "++ unwords ls) $ - filter (not . null) $- map words $ - lines contents- - return$ - \ name_to_hack -> - case M.lookup name_to_hack mp of -- Could use a trie- Just x -> x- Nothing -> name_to_hack--temp = driver default_phybin_config{ num_taxa=7, inputs=["../datasets/test.tr"] }--------------------------------------------------------------------------------------------------------- TODO: expose a command line argument for testing.--- The below test exposed my normalization bug relating to "deriving Ord".--- I need to transform it into one or more proper unit tests.--main_test = - withArgs ["-w","~/newton_and_newton_local/datasets/yersinia/yersinia_trees/111.dnd","-m","../datasets/yersinia/name_table_hack_yersinia.txt"]- main ----a :: AnnotatedTree--- annotateWLabLists$ -a :: NewickTree Double-a = set_dec 1 $ - NTInterior () [NTInterior () [NTLeaf () "RE",NTInterior () [NTLeaf () "SD",NTLeaf () "SM"]],NTInterior () [NTLeaf () "BB",NTLeaf () "BJ"],NTInterior () [NTLeaf () "MB",NTLeaf () "ML"]]--b :: NewickTree Double-b = set_dec 1 $ - NTInterior () [NTInterior () [NTLeaf () "BB",NTLeaf () "BJ"],NTInterior () [NTLeaf () "MB",NTLeaf () "ML"],NTInterior () [NTLeaf () "RE",NTInterior () [NTLeaf () "SD",NTLeaf () "SM"]]]--ls = [("a",a),("b",b)]---- This is one:-num_binned = M.size $ binthem ls--a_ = ("980.dnd",NTInterior 0.0 [NTInterior 5.697e-2 [NTLeaf 3.95e-2 "SM",NTLeaf 5.977e-2 "SD"],NTLeaf 0.13143 "RE",NTInterior 0.13899 [NTInterior 9.019e-2 [NTLeaf 0.11856 "BB",NTLeaf 0.13592 "BJ"],NTInterior 0.13194 [NTLeaf 0.19456 "MB",NTLeaf 0.16603 "ML"]]])--b_ = ("999.dnd",NTInterior 0.0 [NTInterior 6.527e-2 [NTInterior 0.13734 [NTLeaf 2.975e-2 "SM",NTLeaf 3.002e-2 "SD"],NTLeaf 0.18443 "RE"],NTInterior 6.621e-2 [NTLeaf 0.16184 "MB",NTLeaf 0.15233 "ML"],NTInterior 0.23143 [NTLeaf 9.192e-2 "BB",NTLeaf 0.10125 "BJ"]])---- But THIS is two: ack!-num2 = M.size $ binthem [a_,b_]---- Here's the test that breaks things:-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 = 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"]]]--d1 = drawNewickTree "" a_norm-d2 = drawNewickTree "" b_norm--d1_ = forkIO $ do runGraphvizCanvas Dot (dotNewickTree_debug "" a_norm) Xlib; return ()-d2_ = forkIO $ do runGraphvizCanvas Dot (dotNewickTree_debug "" b_norm) Xlib; return ()
+ Main.hs view
@@ -0,0 +1,300 @@+{-# LANGUAGE RecordWildCards #-}+{-# OPTIONS_GHC -fwarn-unused-imports #-}++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 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 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 Bio.Phylogeny.PhyBin.Visualize (viewNewickTree, dotNewickTree_debug)++import Version++----------------------------------------------------------------------------------------------------+-- MAIN script: Read command line options and call the program.+----------------------------------------------------------------------------------------------------++-- Note: ORDER is important here, we process options in this order:+data Flag + = Verbose + | Version + | Output String+ | NumTaxa Int+ | BranchThresh Double + | NullOpt+ | Graph | Draw+ | Force + | View+ | TabDelimited Int Int++ | SelfTest++ | NameCutoff String+ | NamePrefix Int+ | NameTable String -- Must come after Cutoff/Prefix++ deriving (Show, Eq, Ord) -- ORD is really used.++parseTabDelim :: String -> Flag+parseTabDelim _str = + TabDelimited 8 9+ +options :: [OptDescr Flag]+options =+ [ 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 [] ["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") "----------------------------- Visualization --------------------------------"+ , Option ['g'] ["graphbins"] (NoArg Graph) "use graphviz to produce .dot and .pdf output files named bin1.*, bin2.*, etc"+ , 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 ['n'] ["numtaxa"] (ReqArg (NumTaxa . read) "NUM") "expect NUM taxa for this dataset"++ , Option ['b'] ["branchcut"] (ReqArg (BranchThresh . read) "LEN") "collapse branches less than LEN"+ ++{- -- 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."++ , Option [] [] (NoArg NullOpt) ""+ , Option ['s'] ["namesep"] (ReqArg NameCutoff "STR") $ --"\n"+++ "An alternative to --nameprefix, STR provides a set of delimeter characters,\n"+++ "for example '-' or '0123456789'. The taxa name is then a variable-length\n"+++ "prefix of each gene name up to but not including any character in STR."++ , Option [] [] (NoArg NullOpt) ""+ , Option ['m'] ["namemap"] (ReqArg NameTable "FILE") $ + "Even once prefixes are extracted it may be necessary to use a lookup table\n"+++ "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."+ ]++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"++++ "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"++++ "\nOptions include:\n"++defaultErr :: [String] -> t+defaultErr errs = error $ "ERROR!\n" ++ (concat errs ++ usageInfo usage options)+++allUnitTests :: Test+-- allUnitTests = unitTests +++allUnitTests = test + [ Bio.Phylogeny.PhyBin.unitTests+ , Bio.Phylogeny.PhyBin.Parser.unitTests+ ]+-- Bio.Phylogeny.PhyBin.Parser.unitTests++main :: IO ()+main = + do argv <- getArgs ++ (opts,files) <- + case getOpt Permute options argv of+ (o,n,[] ) -> return (o,n)+ (_,_,errs) -> defaultErr errs++ 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+ + Output s -> return cfg { output_dir= s }++ NumTaxa n -> return cfg { num_taxa= n }+ BranchThresh n -> return cfg { branch_collapse_thresh= Just n }+ Graph -> return cfg { do_graph= True } + Draw -> return cfg { do_draw = True } + View -> return cfg -- Handled below++ TabDelimited _ _ -> error "tabbed option not handled yet"+ 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 }++ -- 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 } + (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++view_graphs :: PhyBinConfig -> IO ()+view_graphs PBC{..} = + do chans <- forM inputs $ \ file -> do + 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+ return chan+ forM_ chans readChan + return ()+++--------------------------------------------------------------------------------+-- Every dataset it seems needs a new hack on the names!++name_table_reader :: String -> IO (String -> String)+name_table_reader file = + do contents <- readFile file+ let mp = M.fromList $ + map (\ls -> case ls of + [a,b] -> (a,b)+ _ -> error$ "Each line of "++file++" must contain two whitespace free strings: "++ unwords ls) $ + filter (not . null) $+ map tokenize $ + lines contents+ return$ + \ name_to_hack -> + case M.lookup name_to_hack mp of -- Could use a trie+ Just x -> x+ Nothing -> name_to_hack+ where+ tokenize :: String -> [String]+ tokenize line =+ case words line of+ [] -> []+ [one] -> error$"Name table contained bad line: "++ show one+ [one,two] -> [one,two]+ (one:rest) -> [one, unwords rest]++temp :: IO ()+temp = driver default_phybin_config{ num_taxa=7, inputs=["../datasets/test.tr"] }++----------------------------------------------------------------------------------------------------+-- TODO: expose a command line argument for testing.+-- The below test exposed my normalization bug relating to "deriving Ord".+-- I need to transform it into one or more proper unit tests.++main_test :: IO ()+main_test = + withArgs ["-w","~/newton_and_newton_local/datasets/yersinia/yersinia_trees/111.dnd","-m","../datasets/yersinia/name_table_hack_yersinia.txt"]+ main ++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"]]++pb :: NewickTree DefDecor+pb = set_dec (Nothing,1) $ + NTInterior () [NTInterior () [NTLeaf () "BB",NTLeaf () "BJ"],NTInterior () [NTLeaf () "MB",NTLeaf () "ML"],NTInterior () [NTLeaf () "RE",NTInterior () [NTLeaf () "SD",NTLeaf () "SM"]]]++ls1 :: [(String, NewickTree DefDecor)]+ls1 = [("a",pa),("b",pb)]++-- This is one:+num_binned :: Int+num_binned = M.size $ binthem ls1++a_ :: (String, NewickTree DefDecor)+a_ = ("980.dnd",+ fmap (\x -> (Nothing,x)) $ + NTInterior 0.0 [NTInterior 5.697e-2 [NTLeaf 3.95e-2 "SM",NTLeaf 5.977e-2 "SD"],NTLeaf 0.13143 "RE",NTInterior 0.13899 [NTInterior 9.019e-2 [NTLeaf 0.11856 "BB",NTLeaf 0.13592 "BJ"],NTInterior 0.13194 [NTLeaf 0.19456 "MB",NTLeaf 0.16603 "ML"]]])++b_ :: (String, NewickTree DefDecor)+b_ = ("999.dnd",+ fmap (\x -> (Nothing,x)) $ + NTInterior 0.0 [NTInterior 6.527e-2 [NTInterior 0.13734 [NTLeaf 2.975e-2 "SM",NTLeaf 3.002e-2 "SD"],NTLeaf 0.18443 "RE"],NTInterior 6.621e-2 [NTLeaf 0.16184 "MB",NTLeaf 0.15233 "ML"],NTInterior 0.23143 [NTLeaf 9.192e-2 "BB",NTLeaf 0.10125 "BJ"]])++-- But THIS is two: ack!+num2 :: Int+num2 = M.size $ binthem [a_,b_]++-- Here's the test that breaks things:+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"]]]++b_norm :: NewickTree StandardDecor+b_norm = fmap (\ (bl,w,ls) -> StandardDecor bl Nothing w ls) b_norm_++d1 :: IO (Chan (), NewickTree StandardDecor)+d1 = viewNewickTree "" a_norm++d2 :: IO (Chan (), NewickTree StandardDecor)+d2 = viewNewickTree "" b_norm++d1_ :: IO ThreadId+d1_ = forkIO $ do runGraphvizCanvas Dot (dotNewickTree_debug "" a_norm) Xlib; return ()+ +d2_ :: IO ThreadId+d2_ = forkIO $ do runGraphvizCanvas Dot (dotNewickTree_debug "" b_norm) Xlib; return ()+++-- | A of a tree with _____ weights attached to it:+withBootstrap :: String+withBootstrap = "((((A8F330_:0.01131438136322714984,(G0GWK2_:0.00568050636963043226,(Q92FV4_:0.00284163304504484121,((B0BVQ5_:0.00319487112504297311,A8GU65_:0.00000122123005994819)74:0.00279881991324161267,(C3PM27_:0.00560787769333294297,C4K2Z0_:0.00559642713265556899)15:0.00000122123005994819)4:0.00000122123005994819)56:0.00276851661606284868)60:0.00283144414216590342)76:0.00886304965525876697,(A8GQC0_:0.05449879836105625541,(A8F0B2_:0.04736199885985507840,Q4UJN9_:0.02648399728559588939)64:0.00905997055810744446)28:0.00323255855543533657)29:0.02237505187863457132,(Q1RGK5_:0.00000122123005994819,A8GYD7_:0.00000122123005994819)100:0.28299884298270094884)100:0.05776841634437222123,(Q9ZC84_:0.00000122123005994819,D5AYH5_:0.00000122123005994819)99:0.00951976341375833368,Q68VM9_:0.04408933524904214141);"++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);"
phybin.cabal view
@@ -1,11 +1,17 @@ Name: phybin-Version: 0.1.2.1+Version: 0.1.2.4 License: BSD3 License-file: LICENSE Stability: Beta Author: Ryan Newton <rrnewton@gmail.com> Maintainer: Ryan Newton <rrnewton@gmail.com> ++-- Version history:+-- 0.1.2 -- first significant release+-- 0.1.2.1 -- +-- 0.1.2.4 -- new release, several new features+ -- homepage: http://code.haskell.org/phybin homepage: http://people.csail.mit.edu/newton/phybin/ @@ -30,15 +36,20 @@ -- location: http://code.haskell.org/phybin/repo/ Executable phybin- Main-is: Bio/Phylogeny/PhyBin/Main.hs+ Main-is: Main.hs+-- Main-is: Bio/Phylogeny/PhyBin/Main.hs Build-Depends: base >= 3 && < 5, directory, process, containers, unix, - stringtable-atom, filepath, graphviz, prettyclass, fgl,+ stringtable-atom, filepath, +-- graphviz == 2999.11.0.0,+ graphviz >= 2999.12,+ text >= 0.11 && < 0.12,+ prettyclass, fgl, HSH, HUnit, bytestring, -- For bytestring.lazy support: parsec >= 3.1.0 -- deepseq, -- extensions: CPP- extensions: ScopedTypeVariables, RecordWildCards, TypeSynonymInstances+-- extensions: ScopedTypeVariables, RecordWildCards, TypeSynonymInstances GHC-Options: -O2 -- GHC-Options: -O2 -threaded