在Biopython中捕获Genbank文件解析错误

Question

我已经使用genbank Entrez模块下载了一个与this post类似的BioPython文件列表。在随后解析这些文件时，我遇到了一个错误，因为我从Entrez下载的genbank文件是给予基因组不完整的有机体的临时RefSeq的一部分(NZ_CM000913.1)。当我尝试读取这个文件时，我得到一个记录错误，并且我的脚本停止。我正在尝试编写一个函数来避免这些记录。最简单的方法是按大小过滤记录，但我想知道如何更“生物地”地做这件事--测试一个文件是否包含记录，如果不包含，则排除它。当前的ValueError消息被引发，但将停止脚本。

#the error message is something like this    
from Bio import SeqIO
gbkfile = 'NZ_CM000913.1'
SeqIO.read(open(gbkfile), 'gb')

      File "/usr/lib64/python2.6/site-packages/Bio/SeqIO/__init__.py", line 605, in read
        raise ValueError("No records found in handle")

对于我的循环，我可以使用这样的东西：

#filter by length
for gbk in gbklist:
    if len(open(gbk).readlines()) < 50:
        print 'short file: exclude'
    else:
        process_gbk(gbk)

但我想知道是否可以从BioPython中捕获错误消息：

#generate GBK-file exception    
for gbk in gbklist:
    try: 
        SeqIO.read(open(gbk),'gb')
        process_gbk(gbk)
    except BiopythonIO: 
        'this file is not a genbank file'
        pass

Answer 1

你的建议就快到了！这里是为了完成它+一些方便的额外好处(阅读评论)

errorList = []                               # to store your erroneous files for later handling ;)

#generate GBK-file exception 
for gbk in gbklist:
    try: 
        SeqIO.read(open(gbk),'gb')
        process_gbk(gbk)
    except ValueError:                       # handles your ValueError
        print(gbk+' is not a genbank file')  # lets you know the file causing the error "live"
        errorList.append(gbk)                # logs the name of erroneous files in "errorList"
        continue                             # skips straight to the next loop     

Answer 2

实际上，我可以使用Biopython打开NZ_CM000913.1：

>>> from Bio import SeqIO
>>> fname = 'NZ_CM000913.1.gbk'
>>> recs = SeqIO.read(fname, 'genbank')
>>> recs
SeqRecord(seq=UnknownSeq(6760392, alphabet = IUPACAmbiguousDNA(), character = 'N'), id='NZ_CM000913.1', name='NZ_CM000913', description='Streptomyces clavuligerus ATCC 27064 chromosome, whole genome shotgun sequence.', dbxrefs=['Project:47867 BioProject:PRJNA47867'])

是否确定已正确下载文件并且不是空文件？

另外，我注意到您正在向SeqIO.read提供open('NZ_CM000913.1.gbk')。更好(也更容易阅读)的方法是简单地输入文件名(SeqIO.read('NZ_CM000913.1.gbk', 'genbank'))，以防止文件句柄未关闭。