菜谱

使用 GFF 模块从基因序列预测蛋白质。

由 Brad Chapman 贡献

注意：这需要提议的 Biopython GFF 模块，该模块尚未集成到 Biopython 中。您需要先安装它才能运行示例

问题

从 GlimmerHMM 输出文件（GFF3 格式）开始，生成预测蛋白质序列的 FASTA 文件。

解决方案

设置此操作，我们导入所需的模块并将输入 FASTA 文件解析为标准 python 字典，使用 SeqIO。SeqIO 也用于写入输出文件。我们向输出函数传递一个 python 生成器，SeqIO 将一次将生成器转换为 FASTA 文件。这允许脚本泛化到任意大的数据集，而不会遇到内存问题。

from __future__ import with_statement
import sys
import os
import operator

from Bio import SeqIO
from Bio.SeqRecord import SeqRecord

from BCBio import GFF


def main(glimmer_file, ref_file):
    with open(ref_file) as in_handle:
        ref_recs = SeqIO.to_dict(SeqIO.parse(in_handle, "fasta"))

    base, ext = os.path.splitext(glimmer_file)
    out_file = "%s-proteins.fa" % base
    with open(out_file, "w") as out_handle:
        SeqIO.write(protein_recs(glimmer_file, ref_recs), out_handle, "fasta")

接下来，我们提供一个函数来解析 GlimmerHMM 输出。由于 GFF 库，这非常容易。我们一次迭代一组行以处理非常大的输出文件，并提供我们之前使用 SeqIO 解析的参考记录作为输入

def glimmer_predictions(in_handle, ref_recs):
    """Parse Glimmer output, generating SeqRecord and SeqFeatures for predictions"""
    for rec in GFF.parse(in_handle, target_lines=1000, base_dict=ref_recs):
        yield rec

这是程序的核心。我们遍历 GlimmerHMM 输出中的每个记录，然后遍历该记录中的每个基因特征。每个 SeqFeatures 都包含一个基因预测，其中外显子作为顶级特征的 sub_features 存在。使用这些外显子的位置，我们提取当前基因的核苷酸序列并将其追加到列表中。然后，我们将列表加在一起以获得完整序列。如果基因位于反向链上，则序列会被反向互补。最后，我们将序列翻译成蛋白质。生成一个包含序列和预测标识符的 SeqRecord，提供我们之前用来写入输出文件的生成器

def protein_recs(glimmer_file, ref_recs):
    """Generate protein records from GlimmerHMM gene predictions."""
    with open(glimmer_file) as in_handle:
        for rec in glimmer_predictions(in_handle, ref_recs):
            for feature in rec.features:
                seq_exons = []
                for cds in feature.sub_features:
                    seq_exons.append(
                        rec.seq[cds.location.nofuzzy_start : cds.location.nofuzzy_end]
                    )
                gene_seq = reduce(operator.add, seq_exons)
                if feature.strand == -1:
                    gene_seq = gene_seq.reverse_complement()
                protein_seq = gene_seq.translate()
                yield SeqRecord(protein_seq, feature.qualifiers["ID"][0], "", "")

讨论

比较完整脚本

这突出了使用标准库和标准化格式如何节省编码，而不是自己编写。GFF3 版本更加通用，因为它允许同时处理多个重叠群，并且需要更少的自定义代码来解决问题。这使您能够专注于生物学问题，并避免解析代码中的错误。

示例文件可用于运行脚本，感谢 Michael Kuhn 提供

另请参阅