核酸やアミノ酸の配列情報を扱うためのFASTAファイルは、フォーマットがシンプルなので、わざわざBiopythonのモジュールを使わなくとも、たいていのことができると思いますが、知っていると便利にできることがあります。今回は、SeqIOを用いて、簡単な使い方の紹介をします。
準備
環境：python 3.4.2, biopython 1.68 （かならずしもこのとおりでなくてよい）
まず、pythonがインストールされてパスが通っている必要があります。コマンドラインからwhich pythonと打ってパスが表示されればOKです。pythonと打ってみましょう。
Python 3.4.2 (v3.4.2:ab2c023a9432, Oct 6 2014, 22:15:05) [MSC v.1600 32 bit (Intel)] on win32
Type "help", "copyright", "credits" or "license" for more information.
>>>
ここで、pythonのバージョンが表示されています。2.7.xが表示される場合、python3と打つと、python3が立ち上がるかもしれません。Macでは標準で2.7.xが入っていると思います。この紹介ではpython3のほうを使います。Python自体はここからダウンロードできます。
Biopythonは、こちらからダウンロードできます。お使いのpython versionにあったものをインストールしてください。（以前はNumpyを別途インストールする必要があり、そこでコケたりすることがありましたが、今はインストーラーにふくまれています。ありがたいですね）
Biopythonのインストールが完了したら、コマンドラインからpythonと打つか、IDLEを立ち上げて（python -m idlelib.idleと打つか、エイリアス・ショートカットアイコンをダブルクリック）、
from Bio import SeqIO
と打ってみます。
>>> from Bio import SeqIO
>>>
のように何事も起こらなければうまくいっています。exit()と打ってpython promptをいったん終了します。
FASTA format
フォーマットについては下のreference pagesにあるページが参考になります。ここではSeqIOを使う上で理解しておくべきことを書きます。
いま少々世間を騒がせているジカ熱ウイルスのゲノム配列を例にとります。下の例がfasta formatです。ジカウイルスはRNA(ssRNA(+))をゲノムにもちますので、あれ？この配列はDNAじゃないの？と思った方は、鋭い！のですが、このデータベース配列ではUがTとして表記されています。
>KX893855.1 Zika virus strain Zika virus/Homo sapiens/VEN/UF-2/2016, complete genome
AGTTGTTACTGTTGCTGACTCAGACTGCGACAGTTCGAGTTTGAAGCGAAAGCTAGCAACAGTATCAACA
GGTTTTATTTTGGATTTGGAAACGAGAGTTTCTGGTCATGAAAAACCCAAAAAAGAAATCCGGAGGATTC
CGGATTGTCAATATGCTAAAACGCGGAGTAGCCCGTGTGAGCCCCTTTGGGGGCTTGAAGAGGCTGCCAG
===途中略===
TGTGGATCTCCAGAAGAGGGACTAGTGGTTAGAGGAGACCCCCCGGAAAACGCAAAACAGCATATTGACG
CTGGGAAAGACCAGAGACTCCATGAGTTTCCACCACGCTGGCCGCCAGGCACAGATCGCCGAATAGCGGC
GGCCGGTGTGGGGAAATCCATGGGTCTT
さて、fastaフォーマットでは、二種類の行しかありません。ひとつ目は>から始まる１行目。残りの行はすべて配列行です。以上。
といいたいところなのですが、biopythonを使う上で、もう少し知っておくことがあります。Description lineあるいはdefinition line (defline)と呼ばれる一行目についてです。>のあとから最初のスペースまでは.id で取り出すことができます。（accession id/numberがここに置かれることが多い。）>を除くdefline全体は、.descriptionで得られます。
配列部分は.seqで取得することができます。
なにはともあれ、実際にコードを書いてどう動くかみてみます。
# parsingFasta1.py
import sys
from Bio import SeqIO
fasta_in = sys.argv[1]
for record in SeqIO.parse(fasta_in, 'fasta'):
id_part = record.id
desc_part = record.description
seq = record.seq
print('id:', id_part)
print('desc:', desc_part)
print('seq:', seq)
出力。
$ python fastaMethods.py KX893855.fasta
id: KX893855.1
desc: KX893855.1 Zika virus strain Zika virus/Homo sapiens/VEN/UF-2/2016, complete genome
seq: AGTTGTTACTGTTGCTGACTCAGA...省略
id, def line, sequenceを取り出すことができました。
Fasta format fileには上で述べたようにdef lineとsequenceの二種類の行しかなく、取得することのできる情報は、
- id
- def line全体
- sequence
の３つです。biopythonを使わずに同様のことをしようとすると、
fileを開いて
一行ずつみていき
">"で始まる行なら
">" より後ろをスペースで一回だけsplitする
splitの結果[0]がid
">"より後ろがdescription
それ以外の行は全て`+= その行`のようにconcatenateしていって最終的にseqとする
のようにすればできそうです。pythonの練習課題としては楽しいかもしれませんが、biopython を使った方がわかりやすく書けます。
Biopythonを使う意味
今回のようなsimple fasta parsingでは、biopythonをわざわざ使う意味が見えにくいかもしれませんが、特にgenbank fileを扱い始めると便利さが身にしみます。今回はSeqIOの"I"の方しか使いませんでしたが、SeqIOの"O"の方も便利になってきます。少しづつポストしていきます。
- reverse complement
- 15 character systemを難なく扱えて便利（例：S->S, M->K, B->V）。
- multiple fasta
- 複数のレコードが入ったfastaの扱い。
- multiple fasta fileの内容をsingle fasta filesに分割する。
- deflineの一斉書き換え(prefixを足したり、配列長を追記したり)
- ある条件に一致するレコードのみ選択してファイルに書き出す。
- genbank file parsing
- これはpythonのみで一般的に使えるものを作るとするとかなり面倒。便利です。
- genbank -> fasta format変換
- genbak中の必要な情報をfastaのdef lineに入れたりするのに便利。
parsingFasta1.pyの補足説明
fasta_in = sys.argv[1]
コマンドラインから入力fastaへのパスを受けます。
for record in SeqIO.parse(fasta_in, 'fasta'):
この部分は毎回登場する決まった書き方です。Fileにはひとつのレコードしかないのになぜfor loop？はごもっともですがこうしておくと複数レコードを含むファイル(multi fasta)にもそのまま対応できます。SeqIO.parse()にはfile pathとformat をこの順で渡します。
さて、record変数ですが、これはSeqRecord objectです。先ほどのfastaMethods.pyを少し変更してこのobjectが持つmethodとpropertyを覗いてみます。
# fastaMethods_peekVars.py
import sys
from Bio import SeqIO
fasta_in = sys.argv[1]
for record in SeqIO.parse(fasta_in, 'fasta'):
id_part = record.id
desc_part = record.description
seq = record.seq
print('type(record):', type(record))
print('dir(record):', dir(record))
print()
print('type(seq):', type(seq))
print('dir(seq):', dir(seq))
出力は以下の通り。
$ python3 fastaMethods.py KX893855.fasta
type(record): <class 'Bio.SeqRecord.SeqRecord'>
dir(record): ['__add__', '__bool__', '__class__', '__contains__', '__delattr__', '__dict__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattribute__', '__getitem__', '__gt__', '__hash__', '__init__', '__iter__', '__le__', '__le___', '__len__', '__lt__', '__module__', '__ne__', '__new__', '__nonzero__', '__radd__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__', '__weakref__', '_per_letter_annotations', '_seq', '_set_per_letter_annotations', '_set_seq', 'annotations', 'dbxrefs', 'description', 'features', 'format', 'id', 'letter_annotations', 'lower', 'name', 'reverse_complement', 'seq', 'upper']
type(seq): <class 'Bio.Seq.Seq'>
dir(seq): ['__add__', '__class__', '__contains__', '__delattr__', '__dict__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattribute__', '__getitem__', '__gt__', '__hash__', '__init__', '__le__', '__len__', '__lt__', '__module__', '__ne__', '__new__', '__radd__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__', '__weakref__', '_data', '_get_seq_str_and_check_alphabet', 'alphabet', 'back_transcribe', 'complement', 'count', 'endswith', 'find', 'lower', 'lstrip', 'reverse_complement', 'rfind', 'rsplit', 'rstrip', 'split', 'startswith', 'strip', 'tomutable', 'tostring', 'transcribe', 'translate', 'ungap', 'upper']
record変数、seq変数はそれぞれ、SeqRecord, Seq objectであることがわかります。出力を右の方へスクロールすると、たしかにSeqRecordには 'description', 'id'が見えます。ほかにも面白そうなものが見えますので、上のscript中、print(record.xxx)などとして何が得られるか試してみるのもいいと思います。
seq変数は、例えばそのままprint()関数に入れられるなど、stringsのようにふるまいますが、stringではありません。辞書のkeyとしてそのまま使ったりしない方が無難です。上に書いたように、script中でprint(seq.xxx)などとして、どんな結果がでるか見てみましょう。
print(seq.reverse_complement)
print(seq.reverse_complement()) # see the difference
print(seq.transcribe())
record.id, record.descriptionとして得ている、id_part, desc_part変数はともにstringです。