BiopythonでGenbankファイルからアミノ酸配列を抽出する

Genbankファイル中のCDS情報から遺伝子のアミノ酸配列を抽出し、遺伝子ごとのアミノ酸配列が記載されたFASTAファイルを作ったまとめです。

Genbankファイルの扱いはBiopythonを利用すると簡単です。公式のチュートリアルに詳しい説明がありますが、英語だったのでこちらの方の記事を参考にさせていただきました。

大変詳しく説明されているのですが、自分なりにまとめ直したものも備忘録として残しておきます。誰かの理解の助けになれば幸いです。

とりあえず、完成形のコードを下に示します。詳しい解説はこの後で。

#!/usr/bin/env python3

import sys
from Bio import SeqIO

def extract_cds(gbk):
    for record in SeqIO.parse(gbk, "genbank"):
        for feature in record.features:
            if feature.type != "CDS":
                continue
            try:
                product     = feature.qualifiers["product"][0]
                protein_id  = feature.qualifiers["protein_id"][0]
                fas_desc    = "{0}\t{1}".format(product, protein_id)
                translation = feature.qualifiers["translation"][0]
                print(">{0}\n{1}".format(fas_desc, translation))
            except KeyError:
                continue

if __name__ == '__main__':
    genbank = sys.argv[1]
    extract_cds(genbank)

以下、解説です。

from Bio import SeqIO
for record in SeqIO.parse([genbank_file path], "genbank"):
    # 何らかの処理
    ...

BiopythonでGenbankファイルを読み込むときは、上記の形が基本となります。ちなみに、ファイルを読み込みための関数はparseの他にもう1種「SeqIO.read」がありますが、readの方はファイル中に1つの配列情報しか記載されていないときに使います。parseとは異なり、複数の配列情報が記載されているとエラーが発生するので確認が必要な状況では便利です。

上記の基本例に沿って、これ以降record変数からGenbankファイル中の様々な情報を取得します。以下に一部ですがよく使うものをまとめました。

表記	得られる情報
record.id	配列のアクセッションナンバー
record.description	配列のDEFINITION
record.seq	配列そのもの(Seqオブジェクト)
record.annotations	アノテーション情報が格納された辞書
record.annotations["organism"]	生物種名
record.annotations["taxonomy"]	生物分類のリスト
record.features	FEATURESの情報が格納されたリスト

今回、遺伝子のアミノ酸配列を抽出するにあたり重要なのは最終行のrecord.featuresです。以下のように利用します。

# 上からの続き
for feature in record.features:
    # 何らかの処理
    ...

record.featuresにはFEATURESの項目の情報が上から順に格納されています。これをループで処理することで、特定のfeature keyの情報のみを抽出します。今回の場合はCDSですね。上の例に沿って進めると、feature keyの種類はfeature.typeに､そのfeature keyが示す領域の詳しい情報については、Pythonの辞書形式でfeature.qualifiersに格納されています。CDSの場合のfeature.qualifiersの例を一部ですが下にまとめます。

表記	得られる情報
feature.qualifiers["protein_id"]	タンパク質のアクセッションナンバーのリスト
feature.qualifiers["product"]	タンク質のDEFINITIONのリスト
feature.qualifiers["translation"]	アミノ酸配列のリスト

一点注意しなければならないのは、qualifiersではすべての値がリストに格納されていることです。そのため、そこに格納された情報を利用する際にはインデックスをつける必要があります。

ここまでをまとめて、GenbankファイルからCDS情報のみを抜き出す基本形は以下の通りです。

# 上からの続き
for feature in record.features:
    # feature keyの種類がCDSでない場合は除外
    if feature.type != "CDS":
        continue
    qualifiers = feature.qualifiers
    # 何らかの処理
    ...

以上の情報をまとめると、初めのコードの挙動は次のように説明できます。

SeqIO.parseでGenbankファイルを読み込みrecord変数に代入
FEATURESの情報をfeature変数に代入
feature keyがCDSのものだけを対象に、protein_idとproductとtranslationを取得する。

記事の最後に、コードの動作確認に利用したテストデータを載せておきます。けっこう崩してしまっても読み込んでくれるんですね。初めて知りました。

参考リスト

テストデータ

LOCUS       Test
DEFINITION  Dummy
FEATURES             Location/Qualifiers
     CDS             190..255
                     /product="thr operon leader peptide"
                     /protein_id="AAC73112.1"
                     /translation="MKRISTTITTTITITTGNGAG"
     CDS             complement(16751..16960)
                     /product="regulatory protein MokC"
                     /protein_id="AAC73129.1"
                     /translation="MLNTCRVPLTDRKVKEKRAMKQHKAMIVALIVICITAVVAALVTRKDLCEVHIRTGQTEVAVFTAYESE"
     CDS             complement(16751..16903)
                     /product="protein HokC"
                     /protein_id="AAT48122.1"
                     /translation="MKQHKAMIVALIVICITAVVAALVTRKDLCEVHIRTGQTEVAVFTAYESE"
ORIGIN      
        1 agcttttcat tctgactgca acgggcaata tgtctctgtg tggattaaaa aaagagtgtc
//

ほぼ中立ブログ

少しだけ趣味に偏った雑記ブログ

BiopythonでGenbankファイルからアミノ酸配列を抽出する

参考リスト