Biopythonを用いた開始コドンとSD配列を有するORFの抽出
Biopython（公式サイト）を用いたSD配列を上流に有するORF（Open Reading Frame）のうち、開始コドンを最初のアミノ酸として持つ配列の検索と抽出を行います。
- SD配列 …（Shine-Dalgarno：シャイン・ダルガノ）原核生物のmRNAにおいて、翻訳開始のためにリボソームが結合する領域
- 開始コドン …ポリペプチドが翻訳され始めるコドン（ATG:メチオニン（M）が多いがGTGやTTGなども開始点となる場合がある）
なお、大部分のコードはBiopythonを用いたORFの検索・翻訳と配列の抽出と同様です。そのため、説明が重複する箇所は省略します。
ここで記載するORFについての注意事項を再掲します。
- 順方向のみの翻訳とする（公式のクックブックには、逆方向の翻訳の方法も記載されている）。
- エキソン・イントロンは考慮しない（原核生物を対象とする）。
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)
seq_list = [seq]
base_list = list(seq_list[0])
sd_list = [index for index, base in enumerate(base_list) if ''.join(base_list[index:index+5]) == 'GGAGG']
for nuc in seq_list:
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
start_codon = trans[aa_start]
sd2start = nuc[start-11:start-1]
end = min(seq_len, frame + aa_end * 3 + 3)
length = int((end - start) / 3)
if start_codon == 'M' and start-11 in sd_list:
out.append((length, start, end, sd2start, 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()
参考として、Biopythonチュートリアル・クックブック 20.1.13 Identifying open reading framesの内容を用いています。
重要な箇所を以下に示します。
配列リスト取得
base_list = list(seq_list[0])
sd_list = [index for index, base in enumerate(base_list) if ''.join(base_list[index:index+5]) == 'GGAGG']
今回は検索対象とする配列は一つだけであるため、seq_listとしてリスト化する必要はありませんが、seq_list[0]として配列リスト中の最初の配列を指定しています。配列リスト中に複数の配列がある場合には、この指定をループさせる必要があります。
base_listでは、ある配列中の塩基を一つずつリストの要素としています。
たとえば、ATGATGAという配列がある場合には、base_listは
['A','T','G','A','T','G','A']
というように、一塩基ずつ分解されてリストに入れられます。
その後、sd_listとして、内包的表現でSD配列が始まる位置のリストを取得しています。
その際、enumerate()を用いることで、配列のインデックスと要素を同時に取得しています。
出力条件
if start_codon == 'M' and start-11 in sd_list:
アミノ酸配列の最初が開始コドンMであること、およびSD配列の開始点がsd_listと一致したときにcsvファイルにoutputするようにします。
また、のちにSD配列の一塩基目から開始コドンの直前の塩基までの配列を取得するために、
sd2start = nuc[start-11:start-1]
を指定しています。
出力
これにより得られる出力は以下の通りです。例として、大腸菌（Escherichia coli）K-12株 MG1655のゲノム（GenBank: U00096.3）を用いた結果を示します。
50 649216 649368 GGAGGTATAA MRVTTPWRMCSTSGAWQSLRELDARVRFLKPIFAMTSITILTVRSPPRKA
条件に完全に一致する配列は、大腸菌ゲノム中より上記の配列のみ見つかりました。
なお、SD配列はGGAGGに完全に一致せずとも、また開始コドンより上流5塩基を挟まずとも翻訳は開始されますので、目的の配列を幅広く見出すためにはより緩い条件で検索する必要があります。