食譜

使用 GFF 模組將基因序列轉換為預測的蛋白質。

由 Brad Chapman 貢獻

注意：這需要 Biopython 的提議的 GFF 模組，該模組尚未整合到 Biopython 中。您需要先安裝它才能執行範例

問題

從 GlimmerHMM 以 GFF3 格式的輸出檔案開始，產生預測蛋白質序列的 FASTA 檔案。

解決方案

進行設定，我們導入所需的模組，並使用 SeqIO 將輸入的 FASTA 檔案解析為標準的 python 字典。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 輸出中的每個記錄，然後遍歷該記錄中的每個基因特徵。這些 SeqFeature 中的每一個都包含一個基因預測，外顯子作為頂層特徵的子特徵存在。使用這些外顯子的位置，我們提取當前基因的核苷酸序列，並將其附加到列表中。然後我們將列表加總以獲得完整序列。如果基因位於反向鏈上，則序列將進行反向互補。最後，我們將序列翻譯成蛋白質。產生一個具有序列和預測識別碼的 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 版本更通用，因為它允許多個 contig 同時處理，並且需要較少的自訂程式碼來解決問題。這讓您可以專注於生物學問題，並避免解析程式碼中的錯誤。

感謝 Michael Kuhn 提供範例檔案來執行腳本

另請參閱