读取多个blast文件(生物工程)

Question

我试图读取通过向NCBI blast网站提交多个序列而生成的XML文件列表。从每个文件中，我想打印某些信息行。我想要读的文件都给出了后缀"_recombination.xml"。

for file in glob.glob("*_recombination.xml"):
    result_handle= open(file)
    blast_record=NCBIXML.read(result_handle)
    for alignment in blast_record.alignments:
        for hsp in alignment.hsps:
            print "*****Alignment****"
            print "sequence:", alignment.title
            print "length:", alignment.length
            print "e-value:", hsp.expect
            print hsp.query
            print hsp.match
            print hsp.sbjct

脚本首先查找所有带有"_recombination.xml"后缀的文件，然后我希望它读取每个文件，并打印某些行(这几乎是BioPython烹饪书中的一个直接副本)，它似乎就是这样做的。但我得到了以下错误：

Traceback (most recent call last):
File "Scripts/blast_test.py", line 202, in <module>
blast_record=NCBIXML.read(result_handle)
File "/Library/Python/2.7/site-packages/Bio/Blast/NCBIXML.py", line 576, in read
first = iterator.next()
File "/Library/Python/2.7/site-packages/Bio/Blast/NCBIXML.py", line 643, in parse
expat_parser.Parse("", True) # End of XML record
xml.parsers.expat.ExpatError: no element found: line 3106, column 7594 

我不太确定问题出在哪里。我不确定它是否试图遍历它已经读取的文件--例如，关闭这些文件似乎有帮助：

for file in glob.glob("*_recombination.xml"):
    result_handle= open(file)
    blast_record=NCBIXML.read(result_handle)
    for alignment in blast_record.alignments:
        for hsp in alignment.hsps:
            print "*****Alignment****"
            print "sequence:", alignment.title
            print "length:", alignment.length
            print "e-value:", hsp.expect
            print hsp.query
            print hsp.match
            print hsp.sbjct
    result_handle.close()
    blast_record.close()

但它也给了我另一个错误：

Traceback (most recent call last): 
File "Scripts/blast_test.py", line 213, in <module> blast_record.close() 
AttributeError: 'Blast' object has no attribute 'close'

Answer 1

我通常使用解析方法而不是read，也许它可以帮助您：

for blast_record in NCBIXML.parse(open(input_xml)):
    for alignment in blast_record.alignments:
        for hsp in alignment.hsps:
            print "*****Alignment****"
            print "sequence:", alignment.title
            print "length:", alignment.length
            print "e-value:", hsp.expect
            print hsp.query
            print hsp.match
            print hsp.sbjct

并确保在查询blast中使用-outfmt 5生成xml。

Answer 2

我会在Biogeek的回答中添加一个评论，但我不能(还没有足够的声誉)。实际上，他是对的，你应该用

NCBIXML.parse(open(input_xml))

而不是NCBIXML.read(打开(Input_xml))，因为您“试图读取一个XML文件列表”，而对于XML文件列表，您需要解析而不是读取。它解决了你的问题吗？