Biopython：Genbank file parsing cheat sheet

Genbankファイルのparsingはfastaの場合に比べるとやっかいですが、多くの情報を抜き出すことができるので便利です。

ここでは、セクションごとに、情報にアクセスする記法と、どのような返り値(class)が得られるかをまとめます。あとは使い方次第です。

下準備

以下のようにして、与えられたgenbankファイルにあるrecordを順に読んでいくものとします。

import sys
from Bio import SeqIO

gb_file = sys.argv[1]

# LOOP-(1)
for rec in SeqIO.parse(gb_file, 'genbank'):
    # do something

LOCUS

この一行から、以下の情報を取得することができます。

アクセス方法	得られる値のclass
rec.name	str
rec.annotations['topology']	str
rec.annotations['data_file_division']	str
rec.annotations['date']	str

DEFINITION

ここからはただ一つ。

アクセス方法	得られる値のclass
rec.description	str

ACCESSION

バージョンなしのaccessionです。バージョンも含めたい時にはこれではなく、次のVERSIONからの情報を取りましょう。keyが示しているとおり、複数形（'accessions'）なので、中身がひとつでもリストが返ることに注意しましょう。

アクセス方法	得られる値のclass
rec.annotations['accessions']	list

VERSION

バージョン込みのaccessionと、gi、バージョン部分も個別に取得できます。

アクセス方法	得られる値のclass	memo
rec.id	str	accession with version
rec.annotations['sequence_version']	int
rec.annotations['gi']	int

KEYWORDS

これも複数形。リストで返ります。

アクセス方法	得られる値のclass
rec.annotations['keywords']	list

SOURCE

配列が由来する生物種の情報たち（<--日本語の複数形表現）。

アクセス方法	得られる値のclass
rec.annotations['source']	str
rec.annotations['organism']	str
rec.annotations['taxonomy']	list

REFERENCES

複数のreferenceがあることが多く、それらがリストにはいっています。そのひとつづつにアクセスすることになります。ひとつのreferenceはひとつの'Bio.SeqFeature.Reference' objectとして情報を保持しています。このobjectにアクセスすることになります。

アクセス方法	得られる値のclass	memo
rec.annotations['references']	list	list of 'Bio.SeqFeature.Reference' objects

まず、referenceのobjectひとつづつにアクセスするために、以下のようにforループを回すとします。（このループは、先のforループ、LOOP-(1)のもとで回します）

# LOOP-(1.1)
for ref in rec.annotations['references']:
    # access info under the object

このループのもとで、以下のものにアクセスできます。
位置情報がリンクされています。そうです、annotationは徹底的に配列中のある位置に対して行うものだからです。関連referenceも例外ではないんですね。'FeatureLocation' objectに関しては下で述べます。

アクセス方法	得られる値のclass	memo
ref.location	list	list of 'Bio.SeqFeature.FeatureLocation' objects
ref.authors	str	'.authors'なのですが、リストではなくstrです。
ref.journal	str
ref.medline_id	str
ref.pubmed_id	str

FEATURES

いよいよです。
REFERENCES sectionと同様、通常複数のfeatureが記載されていて、一つのfeatureが'Bio.SeqFeature.SeqFeature' objectとして情報を保持しています。これらをリストとして取得した後、ループを回しながらobjectの持つ情報にアクセスしていきます。

まずは、feature objectsの取得から。

アクセス方法	得られる値のclass	memo
rec.features	list	list of 'Bio.SeqFeature.SeqFeature' objects

LOOP-(1)の下で、次のようにループを回します。

# LOOP-(1.2)
for feat in rec.features:
    # access info under feat obj

どのような情報にアクセスできるかは、以下のとおりです。

アクセス方法	得られる値のclass	memo
feat.type	str	'source', 'gene', 'CDS' etc. featureの名前です
feat.location	'Bio.SeqFeature.FeatureLocation' object	annotation位置に関する情報
feat.qualifiers	dict	qualifier:[value]が辞書に入っています

少し覚えにくく、慣れるまで混乱します。ひとつづつ見ていきましょう。

feat.type

上のLOOP-(1.2)のもとで、feat.typeの返り値(str)として得られます。'source', 'gene', 'CDS', 'mRNA', 'exon'などいろいろあります。すべてのfeatureにひとつづつ必ず存在します。

feat.location

ひとつのfeatureにつき、ひとつ存在します。annotationは、配列のある特定の位置に対してつけていくものだからです。genbank fileを開いて見てみると、各featureの最初の行、feature typeの行に位置情報があるのが見えます。feat.locationで得られるのは、'Bio.SeqFeature.FeatureLocation' objectで、以下のようにしてobject内の情報にアクセスします。基本的にはstartとendなのですが、file中に示されている配列とは逆側（相補鎖側）を意味する場合もあるので、strand情報も用意されています。

アクセス方法	得られる値のclass	memo
feat.location.start	下記参照
feat.location.start.position	int
feat.location.end	下記参照
feat.location.end.position	int
feat.location.strand	int	1 for 'top', -1 for 'bottom(complement)' strand

start, endに関しては、type()と、print時の結果とともに、下に説明します。

# 1..5860
feat.location.start          # 'Bio.SeqFeature.ExactPosition', 0
feat.location.start.position # 'int', 0
feat.location.end            # 'Bio.SeqFeature.ExactPosition', 5860
feat.location.end.position   # 'int', 5860

# <2012..3445
feat.location.start          # 'Bio.SeqFeature.BeforePosition', <2011
feat.location.start.position # 'int', 2011
feat.location.end            # 'Bio.SeqFeature.ExactPosition', 3445
feat.location.end.position   # 'int', 3445

# complement(736..>1488)
feat.location.start          # 'Bio.SeqFeature.ExactPosition', 735
feat.location.start.position # 'int', 735
feat.location.end            # 'Bio.SeqFeature.AfterPosition', >1488
feat.location.end.position   # 'int', 1488

ポイントは以下の通りです。

• feat.location.start（あるいは.end）はそのままprint()に入れると数字を表示する。このとき、"<, >"がある場合はそれもつけたものを表示する。
• feat.location.start.positionとすると、"<, >"は除かれたものが得られる。
• 何よりの注意点は、start, endの返す値（.positionをつけても同様）。1..5860では、0, 5860が得られる。これは部分配列（スライス）をとるときにそのまま使うための配慮。例えば、普通の会話でいうところの、1から5860まで、は、スライスでいえば、seq[0:5860]となる。したがって、

s = feat.location.start
e = feat.location.end
fragment = seq[s:e]

として正しい位置のfragmentを取ることができるのだが、同僚に渡すファイルに

s = feat.location.start
e = feat.location.end
print("start..end", str(s) + '..' + str(e))

としたりすると、あれ？startの位置がちょっとちがうぞ、といわれることになる。

ということで、start、endにどのようなものが返るか、については、少し注意が必要です。

feat.qualifiers

feat.qualifiersは、辞書を返します。ファイル中の、"/"で始まる部分が辞書のkey、"="以降がvalueです。注意点は、valueがリストになっていることです。これまでのところ、複数のelementがリストに入っているものに出くわしたことがないのですが、複数あるとしてコードを書かないと思ったように動きません。

探したいvalueがあるならば、辞書のkeyについて順にループを回すのが簡単です。

# never catches what I am looking for
for k in feat.qualifiers.keys():
    if feat.qualifiers[k] == "what_I_am_looking_for':
        # do what I want to do...

# this works for me
for k in feat.qualifiers.keys():
    if "what_I_am_looking_for" in feat.qualifiers[k]:
        # do what I want to do...

ORIGIN

ここは配列です。

アクセス方法	得られる値のclass
rec.seq	'Bio.Seq.Seq' object

返る値は'Bio.Seq.Seq' objectですが、そのままprint()に入れると配列が得られ、stringのようにふるまいます。スライスもとることができます。

# get from 10 to 20 (= 11 bases)
rec.seq[9:20]

# get 50 bases from position 100 (= 100..149)
rec.seq[99:149]

example codeはまた別のポストで。

May26, 2017追記
BiopythonでGenbankファイル操作：CDSをすべてFASTAに書き出すを書きました。