Biopython 解析没有基因组序列的 GBK 文件

Question

我写了一个脚本，使用GenBank文件和Biopython来获取指定基因的序列，这些序列是我的同事们在工作中需要用到的。

最近我们在处理一个新的数据集时遇到了一些问题，结果发现下载的GBK文件里没有包含序列（这在从NCBI的GenBank网站下载时很常见）。当我们用record.seq[start:end]来获取序列时，Biopython并不会报错，而是返回一长串的N。这种情况下，有什么简单的方法可以在一开始就捕捉到这个问题，并让脚本停止并显示错误信息呢？

Answer 1

好的，我找到了一种方法。如果我数一下序列中的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!")