xml2x 0.2 → 0.4
raw patch · 5 files changed
+159/−135 lines, 5 filesdep ~bionew-uploader
Dependency ranges changed: bio
Files
- README +57/−15
- src/Html.hs +11/−7
- src/Options.hs +7/−2
- src/Xml2X.hs +82/−109
- xml2x.cabal +2/−2
README view
@@ -3,7 +3,8 @@ xml2x - convert blast output in XML format, either to a (csv) table suitable for e.g. importing into Excel or OOCalc, or- to HTML. Optionally annotating the output with GO terms.+ to HTML. Optionally annotating the output with GO terms+ and KEGG KOs. INSTALLATION @@ -14,19 +15,14 @@ xml2x [options] xmlfile1 xmlfile2... - Options include --annotations to specify the mapping between- UniProt accessions and GO terms. This file is usually called- "gene_association.goa_uniprot", and is available from the GO- consortium. The file is several GB, you may want to consider- trimming it down a bit by filtering out the automatic (IEA)- annotations.+ Use -v if you are on an interactive terminal to keep track of+ progress. - Additionally, you can use --go-defs to specify the description of- the GO terms, and the output will then be somewhat more- meaningful. The file is usually called "GO.terms_and_ids",- similarly available.+ Output format is specified with -C (CSV) or -H (HTML), with -C+ being the default. Note that only one output format can be used+ at a time. - Output format is specified with -C or -H, with -C being the default.+CSV OUTPUT For CSV output, the following modes are supported @@ -34,14 +30,60 @@ --top - output only the top hit for each input sequence --region - output top hit for regions that overlap <50% - Use -v on an interactive terminal to keep track of progress.+ Use -o to specify an output file, the default is to output to+ standard out. -BUGS+HTML OUTPUT - HTML isn't implemented yet.+ For HTML output, a directory called "blast.d" is created (or+ re-used if already present), and an index is constructed in a file+ named "index.html" in the current directory. The index lists some+ information about the highest scoring blast hit, and links to the+ file displaying the alignment. + The directory contains one HTML file per input+ sequence, and uses a HTML table to rendering the alignments.+ Color codes indicate level of identity (not total match score or+ E-value!), so short, brightly red matches may have lower score than long gray+ ones. Frame (for BLASTX) or strand (for BLASTN) is indicated as+ text for each match.++ The files are named consistently, so if you run BLAST in both+ directions (i.e. swapping -i and -d), you should be able to go+ back and forth by clicking on the sequence names.++ANNOTATIONS++ Options include --annotations to specify the mapping between+ UniProt accessions and GO terms. This file is usually called+ "gene_association.goa_uniprot", and is available from the GO+ consortium [1]. The file is several GB, you may want to consider+ trimming it down a bit by filtering out the automatic (IEA)+ annotations - however, xml2x will first scan the blast output to+ extract only relevant GO annotations, so keeping it all in memory+ is not necessary.++ Additionally, you can use --ontology to specify the description of+ the GO terms, and the output will then be somewhat more+ meaningful. The file is usually called "gene_ontology.obo",+ similarly available [2].++ You can also add KEGG annotations with the -k (or --kegg-organism)+ option. This option takes a file prefix as a parameter, and+ for a prefix $P, expects to find files $P_uniprot.list and+ $P_ko.list. These files are read, and used to mapp KEGG KOs to+ each UniProt hit. Available from [3].++BUGS+ XML parsing is slow, but ndm said he'd look into it. Must be compiled with -smp to avoid huge memory requirements, but the plus side is that with -smp, we use a lot less RAM than AutoFact.++REFERENCES++ [1] http://www.geneontology.org/ontology/gene_ontology.obo+ [2] ftp://ftp.ebi.ac.uk/pub/databases/GO/goa/UNIPROT/+ [3] ftp://ftp.genome.jp/pub/kegg/genes/organisms/
src/Html.hs view
@@ -52,7 +52,9 @@ type Link = String --- generate a result file, and return a link to it+-- | Generate two results: +-- 1. a new file displaying the aligned matches (via genBrfile)+-- 2. a row in the index file, with a link to the former (using writeRow) mkHtml :: Handle -> (BlastRecord -> [[String]]) -> BlastRecord -> IO () mkHtml h writer br1 = do let ls@((q:_):_) = writer br1@@ -68,11 +70,13 @@ record (t,w,rs) = table -- ! [ border 1 ] << (tr << (th << t +++ th << ruler w +++ map hit rs))- where ruler i = table ! [ border 1, cellspacing 0 ] << tr << go i- go i | i > 100 = td ! [ width "97" ] - << font ! [ size "-2" ] << "|" +++ go (i-100)- | otherwise = td ! [ width (show i)] << font ! [ size "-2" ] << ""+ where ruler i = table ! [ border 1, cellspacing 0, width (show i) ] << tr << go i+ go i | i > 100 = td ! [ width "100" , small_font ] + << "|" +++ go (i-100)+ | otherwise = td ! [ width (show i), small_font ] << "" +small_font :: HtmlAttr+small_font = thestyle "font-size:50%" default_bg = "#d0e0ff" -- | Format a BlastHit as a <tr>@@ -90,8 +94,8 @@ -- $ map (\b -> let (f,t) = (q_from b,q_to b) -- in (f,t,bits b/fromIntegral (t-f),aux b)) $ bs - mycell (Cell w c fr) = td ! [width (show w), bgcolor (makeColor c)] - << font ! [size "-2"] << (maybe "" showFrame fr)+ mycell (Cell w c fr) = td ! [width (show w), bgcolor (makeColor c), small_font] + << (maybe "" showFrame fr) showFrame (Frame Plus n) = "+"++show n showFrame (Frame Minus n) = "-"++show n
src/Options.hs view
@@ -29,6 +29,7 @@ ++ "optionally adding GO annotations.\n" data Options = Opts { godef, goann :: Maybe FilePath+ , ko :: [FilePath] -- multiple species allowed , format :: Format , output :: String -> IO () , inputs :: [FilePath]@@ -45,11 +46,14 @@ data Select = All | Top | Reg deriving Eq options :: [OptDescr (Options -> IO Options)]-options = [Option [] ["go-defs"] (ReqArg (\arg opt -> return opt { godef = Just arg }) "File")- "Go ids and defs"+options = [Option [] ["ontology"] (ReqArg (\arg opt -> return opt { godef = Just arg }) "File")+ "GO ontology definition (.obo file)" ,Option [] ["annotations"] (ReqArg (\arg opt -> return opt { goann = Just arg }) "File") "GOA associations" + ,Option ['k'] ["kegg-organism"] (ReqArg (\arg opt -> return opt { ko = arg : ko opt }) "File")+ "Organism prefix for KEGG"+ -- Verbosity ,Option ['v'] ["verbose"] (NoArg (\opt -> return opt {verb = True, count = \n msg -> countIO msg " ..done!\n" n }))@@ -74,6 +78,7 @@ defaultopts :: Options defaultopts = Opts { godef = Nothing, goann = Nothing, format = Csv+ , ko = [] , output = putStr, verb = False, inputs = [] , count = \_ _ -> return . id, select = All }
src/Xml2X.hs view
@@ -2,112 +2,51 @@ module Main where -import Data.ByteString.Lazy.Char8 (ByteString,unpack,pack)+import Data.ByteString.Lazy.Char8 (unpack) import qualified Data.ByteString.Lazy.Char8 as B-import Data.Map (Map) import qualified Data.Map as M-import Data.Set (Set)-import qualified Data.Set as S-import Data.List (intersperse,sortBy,partition,group,nub)+import Data.Map (Map) import Control.Monad import Data.Maybe import System.IO import System.Directory (doesDirectoryExist, createDirectory) import Data.Array.Unboxed-import Data.Word+import Text.Printf -import Bio.Sequence.GOA import qualified Bio.Alignment.BlastData as BD import Bio.Alignment.BlastData (BlastResult)-import Bio.Alignment.BlastFlat-import Bio.Alignment.BlastXML +import Bio.Alignment.BlastXML as X -- new BlastFlat exports this too+import Data.List (intersperse,sortBy,partition,nub)+import Bio.Sequence.GeneOntology+import Bio.Alignment.BlastFlat ++import Annotate import Options import Html main :: IO ()-main = getOptions >>= annotate -- >>= tabulate---- a sequence name and a list of annotation and bit? scores-type AnnotatedSequence = (ByteString, [(Annotation,Int)])--{- | Read files and return list of annotated sequences- 1. build a map of GO terms to description (25K lines - tiny)- 2. collect all proteins names from the XML file (1 pass over .xml)- 3. build a map of proteins to GO terms (1 pass over gene_assoc..), and add GO term desc- 4. revisit .xml, output: sequence, hits?, go terms/summary,...whatever--}-annotate :: Options -> IO () -- [AnnotatedSequence]-annotate opts = do -- (Opts (Just go) (Just goa),[xml],[])- let readXMLs fs = return . concat =<< mapM readXML fs+main = do+ opts <- getOptions when (null (inputs opts)) $ fail "No input files specified!"-- -- goterms are not necessary when we don't also use goanns- (goterms,goanns) <- case goann opts of- Nothing -> return (M.empty,M.empty)- Just goa -> do- ts <- case godef opts of- Nothing -> return M.empty- Just godb -> do- t <- return . goTerms =<< (count opts) 1000 "reading GO terms: "- =<< readGO godb- when (isNothing $ M.lookup (GO 0) t) $ return () -- strictness hack- return t-- proteins <- getProts opts- when (S.member (B.empty) proteins == False) $ return () -- strictness- mt <- readGOA goa >>= (count opts) 1000 "reading GO annotations: "- >>= return . protTerms proteins- when (isNothing $ M.lookup B.empty mt) $ return () -- strictness, again- return (ts,mt)+ (goterms,goanns) <- read_go opts+ keggs <- read_kegg opts - let tabulator = tabulate opts goterms goanns+ let tabulator = tabulate opts goterms goanns keggs+ readXMLs fs = return . concat =<< mapM X.readXML fs+ case format opts of- Csv -> readXMLs (inputs opts) >>= csvize opts tabulator . concatMap results+ Csv -> readXMLs (inputs opts) >>= csvize opts tabulator Html -> readXMLs (inputs opts) >>= htmlize opts tabulator --- ------------------------------------------------------------------------- Build tables of annotation information--- -------------------------------------------------------------------------- | Parse GO.terms_and_ids into a map-goTerms :: [GoDef] -> Map GoTerm GoDef-goTerms = M.fromList . map go2keyval- where go2keyval gd@(GoDef gt _ _) = (gt,gd)---- | Read blast xml output and collect UniProt hits-getProts :: Options -> IO (Set UniProtAcc)-getProts opts = do- hs <- mapM (\f -> (count opts) 100 ("reading accessions from '"++f++"': ")- . concatMap BD.hits . concatMap BD.results =<< readXML f) (inputs opts)- return . S.fromList . map (chop . BD.subject) $ concat hs---- | Convert "UniProt_ABCDEF Blah blah blah" to "ABCDEF"--- Warning: only tested on this format, may break!-chop :: ByteString -> ByteString-chop = B.drop 1 . B.dropWhile (/='_') . head . B.words---- | Using GO definitions and the set of uniprot terms, scan the associations file and--- snarf all relevant data. Requires the annotations to be grouped by protein.-protTerms :: (Set UniProtAcc) -> [Annotation] -> Map UniProtAcc (UArray Int Word16)-protTerms ps as = M.fromAscList . map toArray . partitions . filter isMember $ as- where toArray (Ann up (GO gt1) _:xs) = let gs = nub (gt1:[ gt | Ann _ (GO gt) _ <- xs ])- a = listArray (0,length gs-1) $ map fromIntegral gs- in a `seq` (pack . unpack $ up,a)- toArray [] = error "cannot read an empty list of annotations!"- isMember (Ann up _ _) = up `S.member` ps- partitions xs@(Ann x _ _:_) = let (p1,pps) = span (\(Ann y _ _) -> y==x) xs- in p1:partitions pps- partitions [] = []- -- ------------------------------------------------------------------------ | Implement CSV output+-- Implement CSV output -- --------------------------------------------------------------------- -csvize :: Options -> ([String], BlastRecord -> [[String]]) -> [BlastRecord] -> IO ()-htmlize :: Options -> ([String], BlastRecord -> [[String]]) -> [BlastResult] -> IO ()+-- | build CSV and HTML output from a header and a function to apply to blast records+csvize, htmlize :: Options -> ([String], BlastRecord -> [[String]]) -> [BlastResult] -> IO () csvize opts (header,writer) brs =- do hs <- count opts 10 "Generating output: " brs+ do hs <- count opts 10 "Generating output: " $ concatMap results brs let csvFormat = unlines . map (concat . intersperse "," . map maybeQuote) maybeQuote = filter (/=',') output opts $ csvFormat (header : concatMap writer hs)@@ -123,42 +62,47 @@ hPutStr h htmlfooter -- | read file and return a header, and a function for generating output-tabulate :: Options -> Map GoTerm GoDef -> Map UniProtAcc (UArray Int Word16)- -> ([String], BlastRecord -> [[String]])-tabulate opts gds pts =+tabulate :: Options -> GoDefinitions -> GoAnnotations -> KoAnnotations -> ([String], BlastRecord -> [[String]])+tabulate opts gds pts kos = let header = ["Query","from","to","Target","Description","from","to","ident","bitscore","E-value","direction"]- ++ if isJust (goann opts) then ["GO terms"]- ++ case select opts of All -> []- _ -> ["indirect GO"]- else []+ ++ (if isJust (goann opts) then ["GO terms"]+ ++ case select opts of + All -> case format opts of + Html -> ["Indirect","Hierarchy"]+ _ -> ["Hierarchy"]+ _ -> ["Indirect","Hierarchy"]+ else [])+ ++ (if (not . null . ko) opts then ["KEGG"] else []) wr = case format opts of Html -> showTop Csv -> case select opts of All -> showAll Top -> showTop Reg -> showReg- in (header, wr (isJust (goann opts)) gds pts)+ in (header, wr gds pts kos) -- a writer generates one or more output records from a BlastRecord-type Writer = Bool -> Map GoTerm GoDef -> Map UniProtAcc (UArray Int Word16) -> (BlastRecord -> [[String]])+type Writer = GoDefinitions -> GoAnnotations -> KoAnnotations -> (BlastRecord -> [[String]]) showAll, showTop, showReg :: Writer --- nb: the returns are due to the type of flatten :: [BlastRecord] -> ..-showAll go gds pts = map show1 . flatten . return- where show1 bf = showFlat bf ++ if go then [showGo gds pts [subject bf]] else []+-- nb: the 'return' is due to flatten needing a *list* of BlastRecord+showAll gds pts ks = map show1 . flatten . return+ where show1 bf = showFlat bf + ++ (if not (M.null pts) then [showGo gds pts [bf],showGoHier gds pts [bf]] else [])+ ++ (if not (M.null ks) then [showKegg ks [bf]] else []) -showTop go g p = showTop' go g p . flatten . return+showTop g p k = showTop' g p k . flatten . return -showTop' :: Bool -> Map GoTerm GoDef -> Map UniProtAcc (UArray Int Word16) -> [BlastFlat]- -> [[String]]-showTop' go gds pts = map show1 . select_first+showTop' :: GoDefinitions -> GoAnnotations -> KoAnnotations -> [BlastFlat] -> [[String]]+showTop' gds pts ks = map show1 . select_first where- show1 (bf,go1,go2) = showFlat bf ++ if go then [go1,go2] else []+ show1 (bf,go1,go2,go3,k1) = showFlat bf ++ (if not (M.null pts) then [go1,go2,go3] else [])+ ++ (if not (M.null ks) then [k1] else []) select_first [] = [] -- take all hits for this query select_first (x:xs) = let (as,zs) = span (\y -> query x == query y) xs (bs,ys) = span (\y -> subject x == subject y) as- ysubs = map head $ group $ map subject ys- in (merge (x:bs), showGo gds pts [subject x],showGo gds pts ysubs)+ -- ysubs = map head $ group $ map subject ys -- NB!+ in (merge (x:bs), showGo gds pts [x], showGo gds pts ys, showGoHier gds pts (x:ys), showKegg ks (x:ys)) : select_first zs -- merge all hits against the same subject@@ -171,7 +115,7 @@ in x { q_from = q1, q_to = q2, h_from = s1, h_to = s2 } merge [] = error "needs at least one blast hit to generate output" -showReg go gds pts = concatMap (showTop' go gds pts) . select_region . flatten . return+showReg gds pts ks = concatMap (showTop' gds pts ks) . select_region . flatten . return where select_region [] = [] select_region (x:xs) = let (this' ,next) = span ((query x==).query) xs@@ -203,11 +147,40 @@ dir bf = case aux bf of Strands p p' -> if p==p' then "Fwd" else "Rev" Frame d _ -> if d==Minus then "Rev" else "Fwd" -showGo :: Map GoTerm GoDef -> Map UniProtAcc (UArray Int Word16) -> [UniProtAcc] -> String-showGo gds pts gs = showGD . unions . map fromJust . filter isJust- . map (flip M.lookup pts) . map chop $ gs+-- | Render KEGG inforation from a set of UniProtAcc's.+showKegg :: KoAnnotations -> [BlastFlat] -> String+showKegg ks fs = unwords $ map show $ concatMap (flip M.lookup ks) $ map chop $ map subject fs++-- | Render GO information from a set of UniProtAcc's.+showGo :: GoDefinitions -> GoAnnotations -> [BlastFlat] -> String+showGo gds pts = showGD . nub . concatMap (extractGOs pts) . map subject where- showGD = concatMap (\gd -> show (GO $ fromIntegral gd) ++ showGT (M.lookup (GO $ fromIntegral gd) gds))- showGT (Just (GoDef _ str cls)) = " ("++show cls++": "++unpack str++") "- showGT Nothing = " "- unions = nub . concatMap elems+ showGD = concatMap (\gt -> show gt ++ showGT (liftM fst $ M.lookup gt gds))+ showGT (Just (GoDef _ str cls)) = " ("++show cls++": "++unpack str++") "+ showGT Nothing = " "++-- | Simlar to "ShowGo", but retain bit scores and extract hierarchical information by+-- recursive lookups.+showGoHier :: GoDefinitions -> GoAnnotations -> [BlastFlat] -> String+showGoHier gds pts = showGD . sortBy (compare `on` (negate . snd)) . M.toList . M.unionsWith max . map getAll+ where+ getAll :: BlastFlat -> Map GoTerm Double+ getAll b = let (s,x) = (subject b, bits b)+ gs = nub $ extractGOs pts s+ in M.fromList $ zip (concatMap (collect gds) $ gs) (repeat x)+ showGD = concatMap (\(gt,x) -> show gt ++ printf " [%.1f]" x ++ showGT (liftM fst $ M.lookup gt gds))+ showGT (Just (GoDef _ str cls)) = " ("++show cls++": "++unpack str++") "+ showGT Nothing = " "++ on f g x y = f (g x) (g y)++-- recursively extract all GoTerms from a starting point+collect :: GoDefinitions -> GoTerm -> [GoTerm]+collect gds t = t : case M.lookup t gds of+ Just (_,[]) -> []+ Nothing -> []+ Just (_,ts) -> concatMap (collect gds) ts++-- convert a UniProtAcc into a list of GoTerms+extractGOs :: GoAnnotations -> UniProtAcc -> [GoTerm] +extractGOs pts = map (GO . fromIntegral) . maybe [] elems . flip M.lookup pts . chop
xml2x.cabal view
@@ -1,5 +1,5 @@ Name: xml2x-Version: 0.2+Version: 0.4 License: GPL License-File: LICENSE @@ -16,7 +16,7 @@ Build-Type: Simple Tested-With: GHC==6.8.2 -Build-Depends: base>3, bio >= 0.3.2, containers, bytestring, array, xhtml, directory+Build-Depends: base>3, bio >= 0.3.4, containers, bytestring, array, xhtml, directory -- Build-Depends: base, bio - GHC 6.6 -- Build-Depends: base, bio, fps - GHC 6.4