Biopythonを用いたORFの検索・翻訳と配列の抽出
Biopython（公式サイト）を用いたORF（Open Reading Frame）の検索と翻訳、および翻訳後の配列の抽出を行います。
index
コード全体
動作環境：
- python 3.6.8
- biopython 1.76
| input | output |
|---|---|
| sequence.fasta | orf.csv |
from Bio import SeqIO
import csv
f = open('orf.csv', 'w', newline="")
w = csv.writer(f)
record = SeqIO.read("sequence.fasta","fasta")
table = 11
min_aa_len = 10
def orfs_trans(seq, trans_table, min_aa_length):
out = []
seq_len = len(seq)
for nuc in [seq]:
for frame in range(3):
trans = str(nuc[frame:].translate(trans_table))
trans_len = len(trans)
aa_start = 0
aa_end = 0
while aa_start < trans_len:
aa_end = trans.find("*", aa_start)
if aa_end == -1:
aa_end = trans_len
if aa_end - aa_start >= min_aa_length:
start = frame + aa_start * 3 + 1
end = min(seq_len, frame + aa_end * 3 + 3)
length = int((end - start) / 3 - 1)
out.append((length, start, end, trans[aa_start:aa_end]))
aa_start = aa_end + 1
out.sort()
return out
orf_list = orfs_trans(record.seq, table, min_aa_len)
w.writerows(orf_list)
f.close()
なお、このコードを実行する際には、前提として以下の条件が伴います。
- ここでの「ORF」とは、開始コドンとしてのメチオニン（M）の有無にかかわらず、終止コドンにいたるまで可能な限り長さを延長させたアミノ酸の配列を指す。したがって、開始コドンに関係なく翻訳は開始されるものとする。
- 順方向のみの翻訳とする（公式のクックブックには、逆方向の翻訳の方法も記載されている）。
- エキソン・イントロンは考慮しない（原核生物を対象とする）。
参考として、Biopythonチュートリアル・クックブック 20.1.13 Identifying open reading framesの内容を用いています。
参考先から変更した主な点は以下の通りです。
- 順方向からのみを扱うため、strandに関する記述を削除
- 複数の変数名の変更
- 結果についてアミノ酸の長さで昇順にソートするように変更
重要な箇所を以下に示します。
アミノ酸の翻訳
まずは用いる配列と条件を指定します。
record = SeqIO.read("sequence.fasta","fasta")
table = 11
min_aa_len = 10
BiopythonのSeqIO.readでFASTAファイルを読み込み、配列を取得します。
tableには翻訳テーブルを番号で指定します（参考：ncbi - The Genetic Codes）。今回の場合はバクテリアを指定しました。
min_aa_lenには、取得するアミノ酸の長さの最小値を指定します。今回は10 a.a.を指定しています。
フレームの指定
for frame in range(3):
trans = str(nuc[frame:].translate(trans_table))
ここでは、翻訳時のフレームの位置ごとにループ処理を宣言しています。
DNAからアミノ酸が翻訳されるときには、3つの塩基ごとに1つのコドンが割り当てられています。
例として、以下の配列を翻訳するときには、
AGCTTTTCATTCTGACTGCAACGG
余る塩基を無視する場合には、以下のような読まれ方が考えられます。
- AGC TTT TCA TTC TGA CTG CAA CGG
- (A) GCT TTT CAT TCT GAC TGC AAC (GG)
- (AG) CTT TTC ATT CTG ACT GCA ACG (G)
一番下の読まれ方を、さらに一塩基分ずらすと、一番上の読まれ方から最初のコドンを抜いたものに等しくなります。
そのため、コドンの読まれ方には大きく3通りが存在します。
このアルゴリズムでは、はじめにすべての塩基を翻訳してしまい、その中から終止コドンと終止コドンに挟まれた箇所を探す仕様になっています。
長さ・開始・終止の取得
while aa_start < trans_len:
aa_end = trans.find("*", aa_start)
if aa_end == -1: # アミノ酸配列の最後に達したら
aa_end = trans_len
if aa_end - aa_start >= min_aa_length: # 指定したアミノ酸配列の長さ以上の時
start = frame + aa_start * 3 + 1
end = min(seq_len, frame + aa_end * 3 + 3)
length = int((end - start) / 3)
out.append((length, start, end, trans[aa_start:aa_end]))
aa_start = aa_end + 1
ここでは、まず、aa_start（翻訳が開始されるアミノ酸の位置）について、trans_len（翻訳したアミノ酸（ポリペプチド）の全長）よりも小さい値である間に、while以下の処理を繰り返します。これは、言い換えると「翻訳したポリペプチド全長の最後まで以下の動作を繰り返す」ということになります。
そして、先ほどの「はじめにすべての塩基を翻訳してしまい、その中から終止コドンと終止コドンに挟まれた箇所を探す」作業を行います。
例えば、
AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTCTGA
を翻訳すると、
- SFSF * LQRAICLCVD * KKSV *
となり、ここから
- SFSF
- LQRAICLCVD
- KKSV
を取り出す、ということを繰り返すことになります。
そのために、aa_endを指定します。
aa_end = trans.find("*", aa_start)
ここでは、pythonのfind関数を用いて、aa_start番目の文字列以降から終止コドン「*」のある位置を探し、それをaa_endに代入します。aa_start = 10ならば、翻訳した全アミノ酸配列のうち、10番目のアミノ酸以降における終止コドンを探すことになります。
その後、開始位置・終止位置・長さを得ます。 なお、開始・終止位置はDNA配列上の位置、長さはアミノ酸配列の長さです。
start = frame + aa_start * 3 + 1 # 開始位置
end = min(seq_len, frame + aa_end * 3 + 3) # 終止位置
length = int((end - start) / 3) # アミノ酸配列の長さ
開始位置は、frameによってずらした分を補正するためにframeを加えています。 終止位置は、配列の長さと終止位置の長さを比較して、終止位置の長さのほうが短い場合はそちらを採用します。 長さは、終止位置から開始位置を引いたものを3（塩基 / 1コドン）で割っています。
出力
これにより得られる出力は以下の通りです。例として、大腸菌（Escherichia coli）K-12株 MG1655のゲノム（GenBank: U00096.3）を用いた結果を示します。
10 16 48 LQRAICLCVD
10 608 640 KLSSIRNLPK
10 2154 2186 HQRWGWITGY
10 2553 2585 RKSFALCWQY
10 2841 2873 ASGLMCSGRR
10 2874 2906 HLLMVHCSEM
10 2884 2916 WCIARRCSHG
10 3090 3122 APVPVRWSRR
10 3310 3342 VAVGAGVSGD
10 4108 4140 ASDHSDRDLR
10 8179 8211 YHQGRPVTSP
10 8861 8893 NLPVLQRARL
10 10577 10609 SASGAVSTIF
10 11424 11456 VHRMTRKGDF
.
.
.
このように、アミノ酸の長さによってソートされたアミノ酸配列が得られました。