How Do You Extract The Query Title From A Blast Xml Output Using Biopython

I would like to find the alignment title along with the query title and the expect value in a BLAST XML output file with many query sequences. I can get the alignment title and expect value but the query title is eluding me. How do I extract that?

blast_records = NCBIXML.parse(result_handle)
blast_record = blast_records.next()

for blast_record in blast_records:
    for alignment in blast_record.alignments:
        for hsp in alignment.hsps:
            print alignment.title, hsp.expect, blast_record.header.query

blast biopython

• 10,332 views • link

updated 14.9 years ago by David W

• written 14.9 years ago by Farhat

Can you give the cloud link of your output file?

• 0 views • link

written 14.9 years ago by Thaman

Thaman: that won't be of much help since any result file should have its query name in the same place.

• 0 views • link

written 14.9 years ago by Michael Schubert

2 answers

Is blast_record.query what you're looking for?

• 0 views • link

written 14.9 years ago by David W

Yes, that does the job.

• 0 views • link

written 14.9 years ago by Farhat

I'm not a python or a BioPython guy, but if your document is large, it could be a bad idea to upload it in memory. Here is a event-bases (SAX) solution:

import sys
import xml
import xml.sax
from  cStringIO import StringIO
from xml.sax.handler import ContentHandler,DTDHandler,EntityResolver
from xml.sax.xmlreader import InputSource

class BlastHandler(ContentHandler,EntityResolver):
    def __init__(self):
        self.content=None
        self.hitdef=None
        self.evalue=None
        self.query=None
    def startElement(self,name,attrs):
        if(name=="Hit_def" or name=="Hsp_evalue" or name=="BlastOutput_query-def"):
            self.content=""
    def endElement(self,name):
        if(name=="Hsp_evalue"):
            self.evalue=self.content
            print self.hitdef,self.evalue,self.query
            self.evalue=None
        elif(name=="Hit_def"):
            self.hitdef=self.content
        elif(name=="BlastOutput_query-def"):
            self.query=self.content
        self.content=None
    def characters(self,chars):
        if(self.content!=None):
            self.content+=chars
    def notationDecl(self, name, publicId, systemId):
        return None
    def unparsedEntityDecl(self, name, publicId, systemId, ndata):
        return None
    def resolveEntity(self, publicId, systemId):
        input = InputSource()
        input.setByteStream(StringIO(""))
        return input

if __name__=='__main__':
    handler=BlastHandler()
    parser=xml.sax.make_parser()
    parser.setContentHandler(handler)
    parser.setEntityResolver(handler)
    parser.parse(open(sys.argv[1]))

• 0 views • link

written 14.9 years ago by Pierre Lindenbaum

The Biopython parser should only deal with the hits for one query at a time - so big multi-query BLAST XML files are not such a problem.

• 0 views • link

written 14.9 years ago by Peter

I would suggest cElementTree Api instead of SAX

• 0 views • link

written 14.9 years ago by Thaman

Out of interest why? Other bits (newer) of Biopython do use ElementTree, but the BLAST XML parser (which is older) went for SAX.

• 0 views • link

written 14.2 years ago by Peter

Biopython parser does deal with only 1 query at a time so memory is not such an issue.

• 0 views • link

written 14.9 years ago by Farhat