diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,3 @@ +#!/usr/bin/env runhaskell +import Distribution.Simple +main = defaultMain diff --git a/blastxml.cabal b/blastxml.cabal new file mode 100644 --- /dev/null +++ b/blastxml.cabal @@ -0,0 +1,25 @@ +Name: blastxml +Version: 0.2 +Synopsis: Library for reading Blast XML output +Description: This library contains a data structure and functions for + reading output from the BLAST sequence alignment program. +Homepage: http://biohaskell.org/ +License: LGPL +Author: Ketil Malde +Maintainer: ketil@malde.org +Stability: Experimental +Category: Bioinformatics +Build-type: Simple +Cabal-version: >=1.6 + +Source-Repository this + Type: darcs + Location: http://malde.org/~ketil/biohaskell/blastxml + Tag: 0.2 + +Library + Exposed-modules: Bio.BlastXML + Other-modules: Bio.BlastData + Build-depends: base >= 3 && < 5, bytestring, parallel, tagsoup + Hs-Source-Dirs: src + Ghc-Options: -Wall diff --git a/src/Bio/BlastData.hs b/src/Bio/BlastData.hs new file mode 100644 --- /dev/null +++ b/src/Bio/BlastData.hs @@ -0,0 +1,61 @@ +{- | + This module implements a hierarchical data structure for BLAST results. + + BLAST is a tool for searching in (biological) sequences for + similarity. This library is tested against NCBI-blast version + 2.2.14. There exist several independent versions of BLAST, so expect some + incompatbilities if you're using a different BLAST version. + + For parsing BLAST results, the XML format (blastall -m 7) is by far the most + robust choice, and is implemented in the "Bio.Alignment.BlastXML" module. + + The format is straightforward (and non-recursive). + For more information on BLAST, check + +-} + +module Bio.BlastData where + +import Data.ByteString.Lazy.Char8 (ByteString) + +-- ------------------------------------------------------------ + +-- | The sequence id, i.e. the first word of the header field. +type SeqId = ByteString + +-- | The 'Strand' indicates the direction of the match, i.e. the plain sequence or +-- its reverse complement. +data Strand = Plus | Minus deriving (Read,Show,Eq) + +-- | The Aux field in the BLAST output includes match information that depends +-- on the BLAST flavor (blastn, blastx, or blastp). This data structure captures +-- those variations. +data Aux = Strands !Strand !Strand -- ^ blastn + | Frame !Strand !Int -- ^ blastx + deriving (Show,Eq) + +-- | A 'BlastResult' is the root of the hierarchy. +data BlastResult = BlastResult + { blastprogram, blastversion, blastdate :: !ByteString + , blastreferences :: !ByteString + , database :: !ByteString + , dbsequences, dbchars :: !Integer + , results :: [BlastRecord] } + deriving Show + +-- | Each query sequence generates a 'BlastRecord' +data BlastRecord = BlastRecord { query :: !SeqId, qlength :: !Int + , hits :: [BlastHit] } deriving Show + +-- | Each match between a query and a target sequence (or subject) +-- is a 'BlastHit'. +data BlastHit = BlastHit { subject :: !SeqId, slength :: !Int + , matches :: [BlastMatch] } deriving Show +-- | A 'BlastHit' may contain multiple separate matches (typcially when +-- an indel causes a frameshift that blastx is unable to bridge). +data BlastMatch = BlastMatch { bits :: !Double, e_val :: !Double + , identity :: (Int,Int) + , q_from, q_to, h_from, h_to :: !Int + , qseq, hseq :: !ByteString + , aux :: !Aux } deriving Show + diff --git a/src/Bio/BlastXML.hs b/src/Bio/BlastXML.hs new file mode 100644 --- /dev/null +++ b/src/Bio/BlastXML.hs @@ -0,0 +1,128 @@ +{-# Language OverloadedStrings #-} + +{- | Parse blast XML output. + + If you use a recent version of NCBI BLAST and specify XML output (blastall -m 7), + this module should be able to parse the result into a hierarchical 'BlastResult' + structure. + + While the process may consume a bit of memory, the parsing is lazy, + and file sizes of several gigabytes can be parsed (see e.g. the + xml2x tool for an example). To parse XML, we use + 'Text.HTML.TagSoup'. +-} + +module Bio.BlastXML ( readXML + , module Bio.BlastData) where + +import Bio.BlastData + +import qualified Data.ByteString.Lazy.Char8 as B +import Text.HTML.TagSoup +import Control.Monad +import Control.Parallel + +type STag = Tag B.ByteString + +-- | Parse BLAST results in XML format +readXML :: FilePath -> IO BlastResult +readXML fp = do + fc <- B.readFile fp + when (not $ B.isPrefixOf " isTagOpenName "Iteration" t || isTagOpenName "Hit" t) ts + return $ xml2br h iters + +-- | breaks p = groupBy (const (not.p)) +breaks :: (a -> Bool) -> [a] -> [[a]] +breaks p (x:xs) = let first = x : takeWhile (not.p) xs + rest = dropWhile (not.p) xs + in rest `par` first : if null rest then [] else breaks p rest +breaks _ [] = [] + +getFrom :: [STag] -> B.ByteString -> B.ByteString +getFrom list tag = let xs = sections (isTagOpenName tag) list + in if null xs || null (head xs) || (null . drop 1 . head) xs + then error ("Couldn't find tag '"++show tag++"' in\n"++showSome list) + else case xs !! 0 !! 1 of + TagText s -> s + x -> error ("Unexpeced tag: "++ show x) + +-- Use pattern match since 'length' is strict, defeating the purpose. +showSome :: [STag] -> String +showSome a@(_:_:_:_:_:_:_) = (init . show . take 5 $ a)++" ... ]" +showSome a = show a + +xml2br :: [STag] -> [[STag]] -> BlastResult +xml2br h is = BlastResult { blastprogram = get "BlastOutput_program" + , blastversion = bv + , blastdate = bd + , blastreferences = get "BlastOutput_reference" + , database = get "BlastOutput_db" + , dbsequences = 0 + , dbchars = 0 + , results = map iter2rec $ breaks (isTagOpenName "Iteration" . head) is + } + where (bv,bd) = B.break (=='[') $ get "BlastOutput_version" + get = getFrom h + +iter2rec :: [[STag]] -> BlastRecord +iter2rec (i:hs) = BlastRecord + { query = get "Iteration_query-def" + , qlength = readI $ get "Iteration_query-len" + , hits = map hit2hit hs + } + where get = getFrom i + +iter2rec [] = error "iter2rec: got empty list of sections!" + +hit2hit :: [STag] -> BlastHit +hit2hit hs = BlastHit + { subject = get "Hit_def" + , slength = readI $ get "Hit_len" + , matches = map hsp2match $ partitions (isTagOpenName "Hsp") hs + } + where get = getFrom hs + + +readI :: B.ByteString -> Int +readI x = case B.readInt x of + Just (n,_) -> n + _ -> error ("Couldn't read an Int from string: '"++B.unpack x++"'") + +readF :: B.ByteString -> Double +readF = read . B.unpack + +hsp2match :: [STag] -> BlastMatch +hsp2match ms = BlastMatch + { bits = readF $ get "Hsp_bit-score" + , e_val = readF $ get "Hsp_evalue" + , q_from = readI $ get "Hsp_query-from" + , q_to = readI $ get "Hsp_query-to" + , h_from = readI $ get "Hsp_hit-from" + , h_to = readI $ get "Hsp_hit-to" + , identity = (readI $ get "Hsp_identity", readI $ get "Hsp_align-len") + , qseq = get "Hsp_qseq" + , hseq = get "Hsp_hseq" + -- blastx has query-frame ±1..3 + -- tblastn has hit-frame + -- blastn has both hit and query + -- tblastx has query-frame = 1, hit-frame ±1..3 + , aux = case sections (isTagOpenName "Hsp_hit-frame") ms of + [] -> mkFrame $ get "Hsp_query-frame" + [(_o:TagText hf:_c)] -> case sections (isTagOpenName "Hsp_query-frame") ms of + [] -> mkFrame hf + [(__o:TagText qf:__c)] -> mkStrands hf qf + e -> error ("hsp2match: should be tagopen/text/close:\n"++show e) + e -> error ("hsp2match: failed to determine frame:\n"++show e) + } + where get = getFrom ms + mkFrame f = Frame (strand' $ signum $ readI f) (abs $ readI f) + mkStrands h q = Strands (strand' $ readI h) (strand' $ readI q) + -- ignore frame also for tblastx hits (it can be reconstructed from location) + strand' :: Int -> Strand + strand' s = if s > 0 then Plus else Minus + +