BioPythonによるバイオデータ分析入門
BioPythonパッケージを使うと、Pythonで生物データ分析を簡単に取り扱うことができます。
この記事ではその入門編として、BioPythonを用いたDNA配列データの基本操作について整理しました。
BioPythonとは
データ分析を目的としたプログラミングにはPythonが広く活用されるようになってきました。
BioPythonは、生物データの分析をPythonで容易に行うために開発されたパッケージです。
使用するコンピュータにPythonをインストールした後、pipコマンドなどで無料で簡単にインストールすることができます。
適する開発環境は人それぞれですが、jupyter notebookでの実行がおすすめです。
DNA配列の定義
今回は、適当な塩基配列データを定義して、それを色々と編集してみます。
1 2 3 |
import Bio from Bio.Seq import Seq test_seq = Seq("ATGCATGCAAAAAAA") |
まずはBioPythonをimportしたら、Seq関数で括って指定する文字列が遺伝子配列であることを宣言します。
それを今回はtest_seqという変数に格納しました。このtest_seqを色々といじってみましょう。
DNA配列の操作
配列の長さの取得
配列の長さを知りたい場合は、Pythonのリスト操作と同様にlenを用いれば一瞬です。
1 |
len(test_seq) |
1 |
15 |
配列の一部分の取得
配列の一部分を取り出す方法もPythonのリスト操作と同一です。
例えば、3文字目から6文字目を取り出したかったら、
1 |
test_seq[3:6] |
1 |
Seq('CAT') |
初めから6文字を取り出したかったら、
1 |
test_seq[:6] |
1 |
Seq('ATGCAT') |
お尻から9文字を取り出したかったら、
1 |
test_seq[-9:] |
1 |
Seq('GCAAAAAAA') |
というように書けば、配列の任意の箇所を引っ張ってくることができます。
逆方向配列の取得
逆方向の配列もPythonのリスト操作と同様の書き方をすることが可能で、以下のように取得します。
1 |
test_seq[::-1] |
1 |
Seq('AAAAAAACGTACGTA') |
相補差配列の取得
BioPythonには相補鎖配列を取得する独自のメソッドが用意されています。
1 |
test_seq.complement() |
1 |
Seq('TACGTACGTTTTTTT') |
また、上記の結果を反転させても良いのですが、以下のようにすれば逆相補鎖も取得できます。
1 |
test_seq.reverse_complement() |
1 |
Seq('TTTTTTTGCATGCAT') |
転写されたRNA配列の取得
RNAに転写された配列を取得するには以下のようにします。
1 |
test_seq.transcribe() |
1 |
Seq('AUGCAUGCAAAAAAA', RNAAlphabet()) |
T(チミン)がU(ウラシル)に置き換わっていることが分かります。
Seq関数の第2引数にRNAであることが明示されていますが、特に気にする事なくこの変数は扱っていけます。
翻訳されたアミノ酸配列の取得
タンパク質(アミノ酸配列)に翻訳された配列も簡単に取得できます。
1 |
test_seq.translate() |
1 |
Seq('MHAKK', ExtendedIUPACProtein()) |
アミノ酸の慣用名はIUPAC命名法に従ってアルファベット1文字で表現されています。
GC含有量の計算
GC含有量を計算するメソッドもあります。
これには、下記のようにSeqUtilsというクラスをimportしておく必要があります。
1 2 |
from Bio.SeqUtils import * GC(test_seq) |
1 |
26.666666666666668 |
test_seqのGC含有量は約26.7%だと分かりました。
また、コドンの場所ごとのGC含有量を計算するGC123というメソッドもあります。
1 |
GC123(test_seq) |
1 |
(26.666666666666668, 40.0, 20.0, 20.0) |
数値が4つ出てきましたが、これらはそれぞれ、
・配列全体のGC含有量
・塩基の第1ポジションのみのGC含有量
・塩基の第2ポジションのみのGC含有量
・塩基の第3ポジションのみのGC含有量
となっています。
構成比率の計算
さらに細かく配列の構成情報を見たい場合は、Python自体の様々なメソッドを活用すれば色々と調べることができます。
例えば、Gだけの数を知りたいなら、
1 |
test_seq.count("G") |
1 |
2 |
とします。上述したlenと組み合わせれば含有割合の計算もできます。
一気にすべての要素数を知りたければ、collectionsモジュールを使用する方法もあります。
1 2 |
from collections import Counter Counter(test_seq) |
1 |
Counter({'A': 9, 'T': 2, 'G': 2, 'C': 2}) |
4つくらいならcountを4回使っても大した手間では無いですが、扱う配列がアミノ酸配列ならばこちらの方が良さそうです。
以上のように、DNA情報の基本的な操作・計算は全てシンプルに行えるようになっています。
生物データの分析を行うならば、まず何はともあれBioPythonをインストールすると良いかと思います。