1

我已经使用 BioPython Entrez 模块按照这篇文章的内容下载了 genbank 文件列表。在随后解析这些文件时，我遇到了一个错误，因为我从 Entrez 下载的 genbank 文件是临时 RefSeq 的一部分，该临时 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
4

2 回答 2

2

你的建议差不多了！这是完成它+一些方便的津贴（阅读评论）

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
于 2013-11-23T12:41:22.613 回答
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'])

您确定您已正确下载文件而不是空文件吗？

另外，我注意到你正在喂食open('NZ_CM000913.1.gbk')。SeqIO.read最好（并且更容易阅读）简单地提供文件名（SeqIO.read('NZ_CM000913.1.gbk', 'genbank')）以防止文件句柄未关闭。

于 2012-12-08T01:35:58.940 回答