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を取得する。
記事の最後に、コードの動作確認に利用したテストデータを載せておきます。けっこう崩してしまっても読み込んでくれるんですね。初めて知りました。
参考リスト
- Package SeqIO | Biopython doc
- Biopython Tutorial and Cookbook | Biopython doc
- Genbank形式 | GenomeMatcher project homepage
- DDBJ/GenBank | biopapyrus
テストデータ
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 //