1

私はBioPythonEntrezモジュールを使用して、この投稿の行に沿ってgenbankファイルのリストをダウンロードしました。その後、これらのファイルを解析すると、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
4

2 に答える 2

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.readSeqIO.read('NZ_CM000913.1.gbk', 'genbank')ファイルハンドルが閉じられないままになるのを防ぐために、ファイル名()を単にフィードする方が良い(そして読みやすい) 。

于 2012-12-08T01:35:58.940 に答える
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 に答える