我编写了一个脚本,它使用 GenBank 文件和 Biopython 从 GBK 文件的序列部分中获取给定基因的序列,我的同事将其用于他们的工作。
我们现在在使用新数据集时遇到了一些问题,结果表明下载的 GBK 文件不包含序列(当您从 NCBI 的 GenBank 网站下载时很容易发生这种情况)。Biopython 在使用record.seq[start:end]
. 从一开始就捕获该问题以停止带有错误消息的脚本的最简单方法是什么?
对了,我找到了办法。如果我计算序列中的 N 并检查序列是否与序列一样多,我知道序列丢失了:
import sys
from Bio import SeqIO
for seq_record in SeqIO.parse("sequence.gb", "genbank"):
sequence = seq_record.seq
if len(sequence) == sequence.count("N"):
sys.exit("There seems to be no sequence in your GenBank file!")
我更喜欢检查序列类型的解决方案,因为空序列是Bio.Seq.UnknownSeq
,而不是Bio.Seq.Seq
真正的序列,如果有人能在那个方向提出建议,我将不胜感激。
更新
@xbello 让我再次尝试检查序列类型,现在这也有效:
import sys, Bio
from Bio import SeqIO
for seq_record in SeqIO.parse("sequence.gb", "genbank"):
sequence = seq_record.seq
if isinstance(sequence, Bio.Seq.UnknownSeq):
sys.exit("There seems to be no sequence in your GenBank file!")