diff --git a/README b/README --- a/README +++ b/README @@ -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/ diff --git a/src/Html.hs b/src/Html.hs --- a/src/Html.hs +++ b/src/Html.hs @@ -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 @@ -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 diff --git a/src/Options.hs b/src/Options.hs --- a/src/Options.hs +++ b/src/Options.hs @@ -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 } diff --git a/src/Xml2X.hs b/src/Xml2X.hs --- a/src/Xml2X.hs +++ b/src/Xml2X.hs @@ -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 diff --git a/xml2x.cabal b/xml2x.cabal --- a/xml2x.cabal +++ b/xml2x.cabal @@ -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