Biopython：Genbank file parsing cheat sheet

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

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

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

```py3

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)のもとで回します）

```py3
# 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)の下で、次のようにループを回します。

```py3
# 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時の結果とともに、下に説明します。

```py3
# 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]`となる。したがって、

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

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

```py3
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について順にループを回すのが簡単です。

```py3

# 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のようにふるまいます。スライスもとることができます。

```py3
# 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に書き出す](http://qiita.com/sato32ha6/items/b9d85fad2d9de53374bf)を書きました。