Biopython入門 – 前編 【インストール ~ 簡単な使い方】

Biopythonは生命科学研究に特化したPythonのパッケージであり、さまざまな情報解析をサポートすることができます。

しかしながら日本語で書かれている情報はかなり少ないので、Biopythonをこれから勉強してみようという方向けに3回構成でその基本的な使い方を解説します。

Biopythonとは

データサイエンス界隈で勢い盛んなプログラミング言語Pythonですが、生命科学の世界にもPythonはかなり使われるようになりました。

生命科学における特有のデータ処理 (DNA・RNA・タンパク質配列や各種データベースとの連携など) を容易にする目的で開発されたのがBiopythonパッケージです。

PerlやRuby等の言語におけるBioPerlBioRubyと同じく無料で公開されており簡単にインストールすることができます。

他のPythonのパッケージと同じくpipで入れることができ、ターミナルから

pip install biopython

と打ち込むだけです。

あるいはGoogle Colabにもインストールすることができて、その場合は

! pip install biopython

というように先頭に「!」をつけます。!がついた場所はターミナルのコマンドと同じようになるのです。

Google Colabについての簡単な入門記事はGoogle Colaboratoryでデータサイエンスを始めよう【使い方入門】にまとめています。

Biopythonでできることは

できること (の一部)
ファイルの操作: Blast、Clustalw、FASTA、GenBank、PubMed、UniGeneなど

オンラインデータベースへのアクセス・操作:Blast、Entrez、PubMed、Swiss-Prot、Prositeなど

各種プログラムの実行:Blast、Clustalw、EMBOSS command line tools

このうち、この記事では一番上のファイルの操作、特にBiopythonにおける最も基本操作である文字の入出力を紹介します。

Biopythonによる配列操作入門

DNA配列は文字列で表すことができるので、BiopythonのSeq() を使って配列を作ることができます。

例えばAGTACACTGGTという配列のDNAを作るには

from Bio.Seq import Seq
my_seq = Seq("AGTACACTGGT")
print (my_seq)

200512 1

基本的にPythonの文字列と同じように操作でき、たとえば部分的な配列を取り出すには

# 3塩基目から最後まで
print (my_seq[2:])
# 先頭から4塩基目まで
print (my_seq[:4])
# 逆順に返す
print (my_seq[::-1])

200512 2

この最後については、配列を逆にできるのでしばしば使われるテクニックです。

Pythonにおける文字列の扱いと同じく、複数のDNA配列を結合する事もできます。

sample1 = Seq ('AAAAGG')
sample2 = Seq ('TTCC')
# 2つの配列を結合
concatenate_seq = sample1 + sample2
print (concatenate_seq)

200512 3

Pythonの文字列と同じく、len ()で長さを調べることができます。

len (concatenate_seq)

配列内の文字を数えるcountも使うことができます。例えば、配列の中のGC割合を調べるには

(concatenate_seq.count ('G') + concatenate_seq.count ('C')) / len (concatenate_seq) * 100

実際には、GC割合はBiopythonが用意しているGC関数を使って計算することもできます。

from Bio.SeqUtils import GC
print (GC(concatenate_seq))

このように、配列情報をPythonの文字列と同じ用に扱えることが分かりました。次はもう少し生命科学ででてくる操作を見ていきます。

相補・逆相補、転写と翻訳

与えられた配列を相補鎖や逆相補鎖にしたり、転写やアミノ酸への翻訳などにすることもできます。

相補鎖や逆相補鎖はcomplement()、reverse_complement()でできます。

sequence = Seq ('AGGTCGA')
#相補鎖
print (sequence.complement())
#逆相補鎖
print (sequence.reverse_complement())

200512 4

転写はtranscribe(), 翻訳はtranslate()です。

#転写
print (sequence.transcribe())
#翻訳
print (sequence.translate())

FASTAファイルなどの読み込み

これまでは単純にDNA文字列を扱う方法を見てきましたが、GenbankからダウンロードしたファイルやFASTA形式のものも、Biopythonでは簡単に読み込みできます。

例えば、Genbankでこちらの配列とアノテーションを.gbという拡張子で保存します (右上の「Send to」> Choose Destinationで「File」を選択 > 「Create File」)。

保存したsequence.gbファイルを読み込むには、Biopythonで入出力を担当するSeqIOクラスのSeqIO.read ()を使います。また、配列のIDは.seq, 配列そのものは.seqで取得できます。

# ファイルの読み込み
from Bio import SeqIO
record = SeqIO.read ('sequence.gb', 'genbank')
# 配列のIDは.idで取得できる
print (record.id)
# 配列情報は.seqで取得できる
print (record.seq)

さまざまなアノテーション情報は.annotations内に辞書型で格納されているので、例えば次のようにしていろいろ表示できます。

# アノテーションを表示する
for k in record.annotations.keys ():
  print (k, ':', record.annotations[k])

200512 5

FASTAファイルの書き出し

BiopythonでFASTAファイルに書き出すこともできます。

配列データを作りFASTAに書き込む

FASTAファイルには複数の配列を格納することができますが、それらはSeqRecordクラスを使って作成します。それぞれには配列の他、そのidとdescriptionが必要です。

いくつか作った配列をリストにまとめ、SeqIO.writeでFASTA形式に書き出すことが可能です。

これまでDNAを例に見てきたので、ここではタンパクを扱ってみます。BioAlphabetからgeneric_proteinを使えばOKです。

from Bio.SeqRecord import SeqRecord
from Bio.Alphabet import generic_protein

# 1つ目の配列
rec1 = SeqRecord (
    Seq ('WVDSTPPPGTRVRAVAIYKQSQHMTEVVRRCPHHERCSDSD', generic_protein),
    id= 'AYF55702.1',
    description ='TP53 partial'
)

# 2つ目の配列
rec2 = SeqRecord (
    Seq ('MDFFRVVENQQPPATMPLNVSFTNRNYDLDYDSVQ', generic_protein),
    id= 'NP_002458.2',
    description ='Myc partial'
)

# リスト形式でまとめる
records = [rec1, rec2]
# FASTAファイルへ書き出し
from Bio import SeqIO
SeqIO.write (records, 'example.fa', 'fasta')

作成された’example.fa’をテキストエディタなどで開くか、あるいはUNIX系の簡単なコマンドが使える方はターミナルで

cat example.fa

とすることでも閲覧できます。

そうするとこのように2つの配列がFASTA形式の1つのファイルになっていることが確認できます。
200512 7

Genbank形式のデータをFASTAファイルに書き出す

新たな配列だけでなく、一度Genbank形式で用意したファイルをFASTA形式に変換することもできます。

先程用意したsequence.gbをFASTA形式に変換してみます。やり方はいくつかありますが、SeqIO.parseで読み込んだ後、SeqIO.writeで書き出すという方法をとってみます。

# Genbank形式のファイルを読み込む
records_gen = SeqIO.parse ('sequence.gb', 'genbank')
# FASTA形式にして書き出す
records_fa = SeqIO.write (records_gen, 'sequence.fa', 'fasta')

このように、どのファイル形式なのかをきちんと指定すれば適切に読み書きすることができます。

最後に、きちんとFASTA形式になっているかlessコマンドで確認してみます。lessはcatと同じく画面に表示するUNIXのコマンドですが、一部しか表示されないためデータが大きくなりがちなバイオインフォデータを確認するときによく使われています。

less sequence.fa

200512 8

正しくFASTA形式に変更できているようです。

まとめに代えて

この記事ではBiopythonを紹介し、実際にDNAやタンパクの入出力を中心に使い方を見てきました。

少し長くなったので今回はこれで終わりにしますが、冒頭にも書いたようにいろいろできることがあり、次回は残り2つ、特にBiopythonを使ったBlast等のアラインメントツールや各種データベースへのアクセス法についてまとめます。

Biopython入門 – 中編 【NCBIデータベースやBlastへのアクセス】

また、このチュートリアルは本家のチュートリアル(英語)を参考にして作成していますが、スペースの関係からかなり割愛しています。本家のチュートリアルはたくさんのことがしっかりと書かれているので、こちらも合わせてご覧ください。

関連図書

この記事に関連した内容を紹介している本やサイトはこちらです。

Biopython入門 – 中編 【NCBIデータベースやBlastへのアクセス】

Biopython入門 – 後編 【系統樹・ゲノムデータの可視化】

Biopython Tutorial and Cookbook

Google Colaboratoryでデータサイエンスを始めよう【使い方入門】

今日も【生命医学をハックする】 (@biomedicalhacks) をお読みいただきありがとうございました。当サイトの記事をもとに加筆した月2回のニュースレターも好評配信中ですので、よろしければこちらも合わせてどうぞ