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

【问题标题】：Biopython parsing a GBK file without genome sequenceBiopython解析没有基因组序列的GBK文件
【发布时间】：2014-10-22 11:22:36
【问题描述】：

我编写了一个脚本，该脚本使用 GenBank 文件和 Biopython 从 GBK 文件的序列部分中获取给定基因的序列，我的同事将其用于他们的工作。

我们现在在使用新数据集时遇到了一些问题，结果表明下载的 GBK 文件不包含序列（当您从 NCBI 的 GenBank 网站下载时很容易发生这种情况）。当使用record.seq[start:end] 时，Biopython 不会抛出错误，而是返回一长串 Ns。从一开始就捕获该问题并停止脚本并显示错误消息的最简单方法是什么？

【问题讨论】：

请edit您的问题与MCVE显示您现在正在做什么。请包括 2 个输入 - 一个普通的 GenBank 文件（包含序列）和一个不包含序列的文件。在不确切知道自己在做什么的情况下，很难就如何改变它提出建议。
抱歉，我认为以前使用过 GenBank 和 BioPython 的人不需要输入文件。当然，我应该包含代码。到目前为止，我已经找到了自己的解决方案，所以如果我不尝试在此处找到提供文件的方法，希望您不要介意。无论如何，谢谢，在尝试准备 MCVE 时更容易找到解决方案。

标签： python biopython genbank

【解决方案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!")

【讨论】：

当您说“无序列”或“序列缺失”时，您是指很多“NNNNNNNN”还是只是一个空序列？你能提供一些样本加入吗？如果您只想检查 Seq 是否为 Bio.Seq.UnknownSeq，则可以使用 if isinstance(sequence, Bio.Seq.UnknownSeq) 之类的内容。
如果你去 GenBank 下载一些东西（例如ncbi.nlm.nih.gov/nuccore/U22660.1）并在这样做之前取消选中右侧自定义视图下的“显示序列”框，该文件将不包含任何低于起源特征。不过，如果您使用SeqIO.parse() 读取该文件，然后尝试提取序列，Biopython 将为您提供尽可能多的 N 序列应该是长的。
到目前为止，isinstance() 对我不起作用。这似乎是因为我一直只做from Bio import SeqIO，而从未真正导入Bio。现在可以了，谢谢。 :-)
最好用from Bio.Seq import UnknownSeq导入UnknownSeq，然后只用isinstance(sequence, UnknownSeq)。