サムネがコーヒーの記事は書きかけです。

ORF(オープンリーディングフレーム）の抽出

SHARE

先日ORF(Open Reading Frame）について考える機会があり、配列からの抽出のアルゴリズムを知りたくなったので試しに一からコードを書いてみました。（普通はORFが欲しい時はORFinder等のwebサービスを利用するのが無難です。）

ORFinder NCBI https://www.ncbi.nlm.nih.gov/orffinder/

以下の記事は、本記事より高速化したORF取得アルゴリズムについて解説しています。

ORFをもっと効率よく求めてみる

抽出にあたって、ORFを以下のように定義します。

[box01 title=”ORFの定義”]1.ORFは開始コドン(ATG)から始まること。

2.ORFは終止コドン（TAA,TAG,TGA)で終了すること。

3.メチオニン(M)単体でもORFとみなす。

4.ComplimentarySeqについても読み込みを行う[/box01]

本来生体内では全てのORFが発現するわけではありませんが、今回の抽出操作では全メチオニンを翻訳開始の可能性がある場所とみなし、全てのORF候補配列を網羅して翻訳まで行います。

また、真核生物はスプライシングを行うため、ゲノムDNAからORFを抽出することは不可能です。

そのため、今回の抽出はゲノムDNAではなくスプライスされたmRNAを想定して行います。（なぜか塩基はチミンを使用していますが、、、）

ORF抽出のアウトライン

抽出操作全体を以下の３つに分割します。

転写、翻訳などの配列操作はBiopythonを使えば一発で行えますが、パッケージのインポートがうまくいかなかったので、標準ライブラリのみを使用して配列操作を行います。

(1)下準備

ここでの目標は、FASTA形式で取得した塩基配列を加工して、ATGCの４文字のみを使用した文字列を取得することです。

from functions import transcribe
#functionsから転写用の関数を取得します。

seq_raw = """>FASTA_0210
GCATGAGTGCCCACGAGGAAGTCGTACCGACAAAGCCATCATGATCCTCGATCGATGCCA
GTGGTCGGTGGTCCTCCTCTAGGTACGTGTAGCTGTGAATGTTGTAGAGAGGATGACCAC
CAGTTCTAGGCGGGCCTTTGCGGCAACCAATGATTCATTCAGGACGTGAACCCGTTGAAA"""
#seq_rawにFASTA形式の塩基配列を代入します。

seq = seq_raw[8:].replace("\n","")
#seq にFASTA形式からATGC以外の文字を取り除いた文字列を代入します。

seq_c = transcribe(seq)[::-1]
#seq_cにseqを転写及びリバースした配列を代入します。

以上の操作で、通常塩基配列seqとcomplimentaryseqの２配列を取得することができました。

def transcribe(seq):
    for nucleotide in range(len(seq)):
        target = seq[nucleotide]
        if target == "A":
            target.replace("A","T")
        elif target == "T":
            target.replace("T","A")
        elif target == "G":
            target.replace("C","G")
        elif target == "C":
            target.replace("G","C")
    return seq

(2)関数の作成

ORF抽出のメイン部分だと言っても過言ではないです。

2.1 codon_2_Aminoacid を作成します

本来なら、biopythonを使い一瞬で翻訳することができますが、ORF抽出時に不都合があったので自力でやっています。確かに、biopythonのtranslate()を使用すれば一瞬で翻訳できて良さそうなのですが、ORFの定義にしたがって翻訳をしたときに終止コドンを見落としてしまう可能性を危惧したため、通常翻訳+終止コドンをXとして出力する仕様にしています。

こうすれば、後からXをマーカーとして配列加工が容易にできます。

関数の内容はシンプルで、protein にコドンに対応するアミノ酸を代入していき、返値として出力します。おそらく変数名はproteinではなくaminoacidの方が良かったです、、、

import math
def codon_2_Aminoacid(codon):
        protein = ""
        if codon == "TTT" or codon == "TTC":
            protein += "F"
        elif codon == "TTA" or codon == "TTG" or codon == "CTT" or codon == "CTC" or codon == "CTA" or codon == "CTG":
            protein += "L"
        elif codon == "ATT" or codon == "ATC" or codon == "ATA":
            protein += "I"
        elif codon == "GTT" or codon == "GTC" or codon == "GTA" or codon == "GTG":
            protein += "V"
        elif codon == "TCT" or codon == "TCC" or codon == "TCA" or codon == "TCG":
            protein += "S"
        elif codon == "CCT" or codon == "CCC" or codon == "CCA" or codon == "CCG":
            protein += "P" 
        elif codon == "ACT" or codon == "ACC" or codon == "ACA" or codon == "ACG":
            protein += "T"
        elif codon == "GCT" or codon == "GCC" or codon == "GCA" or codon == "GCG":
            protein += "A"
        elif codon == "TAT" or codon == "TAC":
            protein += "Y"
        elif codon == "CAT" or codon == "CAC":
            protein += "H"
        elif codon == "CAA" or codon == "CAG":
            protein += "Q"
        elif codon == "AAT" or codon == "AAC":
            protein += "N"
        elif codon == "AAA" or codon == "AAG":
            protein += "K"
        elif codon == "GAT" or codon == "GAC":
            protein += "D"
        elif codon == "GAA" or codon == "GAG":
            protein += "E"
        elif codon == "TGT" or codon == "TGC":
            protein += "C"
        elif codon == "TGG":
            protein += "W"
        elif codon == "CGT" or codon == "CGC" or codon == "CGA" or codon == "CGG":
            protein += "R"
        elif codon == "AGT" or codon == "AGC":
            protein += "S"
        elif codon == "AGA" or codon == "AGG":
            protein += "R"
        elif codon == "GGT" or codon == "GGC" or codon == "GGA" or codon == "GGG":
            protein += "G"
        elif codon == "ATG":
            protein += "M"
        elif codon == "TAA" or codon == "TGA" or codon == "TAG":
            protein += "X"
        return protein

 

2.2関数ORFを作成します

def ORF(seq):
#引数としてseqを使用します。
    index_list = []
    for i in range(len(seq)-2):
        codon_i = seq[i] + seq[i+1] + seq[i+2]
        if codon_i == "ATG":
            index_list.append(i)
            print("-------------------------------------------------")
            print(index_list)
            print("-------------------------------------------------")
#全域をスキャンし、メチオニンに該当するコドンのAのインデックスをindex_listに格納します。
    protein_list1 = []
    for j in range(len(index_list)):
        protein = ""
        start_location = index_list[j]
        for k in range(math.floor((len(seq)-start_location)/3)):
            codon_k = seq[start_location+3*k] + seq[start_location+3*k+1] + seq[start_location+3*k+2]
＃作成した関数codon_2_Aminoacidを使用し、変数proteinにアミノ酸を代入していき鎖を伸ばしていきます。
            protein += codon_2_Aminoacid(codon_k)
            print(protein)
        protein_list1.append(protein)
    print(protein_list1)
＃ここまでで、終止マーカーXを含む全ORFを網羅的に取得することができました。
#次に、終止マーカーを削除していく処理を行います。
    protein_list = []
    for target in range(len(protein_list1)):
        target_protein = protein_list1[target]
        counterX = target_protein.count("X")
        if counterX != 0:
            protein_list.append(target_protein)
        else:
            pass
    protein_modified_list = []
    for target_to_cut in range(len(protein_list)):
        target_to_cut = protein_list[target_to_cut]
        protein_modified = ""
        for target_amino_acid in range(len(target_to_cut)):
            target_amino_acid = target_to_cut[target_amino_acid]
            if target_amino_acid == "X":
                break 
            else:
                protein_modified += target_amino_acid
        protein_modified_list.append(protein_modified)
    return protein_modified_list

これで無事、全ORF候補が網羅的に格納されたprotein_modified_listを取得する関数ができました。

ここまでで終了したいところですが、今回はターミナルに出力し確認するところまで行きたいのでもう少し頑張ります。

2.3ターミナル展開用の関数

＃これらの関数は単にprotein_modified_listをターミナルに展開するために作成しています。
def get_answer(protein_list):
    print(protein_list)
    print("----------------------------------")
    for i in range(len(protein_list)):
        print(protein_list[i])
    print("----------------------------------")
def unify_lists(list1,list2):
    list3 = list1 + list2
    list_ = list(set(list3))
    for result in range(len(list_)):
        print(list_[result])

ここまでが、関数の作成です。

ORFの格納を正規配列とリバースの２つに分けて行ったので、unify_listsで統合して出力できるようにします。

(3)処理の実行

(1)および(2)で処理の実行に必要な要素を全て揃えたので、いよいよORF抽出を開始します。

import time
#処理時間を計測します。
from functions import codon_2_Aminoacid , ORF ,get_answer,transcript,unify_lists
＃全関数をインポートします。
from seq import seq , seq_c
#作成した２配列をインポートします。
start = time.time()

modified_protein_list1 = ORF(seq)
modified_protein_list2 = ORF(seq_c)

＃この部分がターミナルにORFのリストを一要素ずつ展開します。
unify_lists(modified_protein_list1,modified_protein_list2)


process_time = time.time() - start
print(process_time)

実行結果

以下がORF抽出の実行結果です。

[box05 title=”実行結果”]MPVVGGPPLGTCSCECCREDDHQF
MILDRCQWSVVLL
MSAHEEVVPTKPS
ML
MLKEHP
MDLLLVAGDRS
MCMDLLLVAGDRS
MTTSSRRAFAATNDSFRT
0.001405954360961914[/box05]

上記のコードは全てコピペ改変等大丈夫ですので、自分のパソコンで処理を実行してみたいという人はやってみてください。

FASTA形式の塩基配列を拾ってきて、自分でORFを取得してみるのも面白いです。

アミノ酸配列から全コドンパターンの逆算
Genome assembly(DNA配列決定）アルゴリズム

 

1 COMMENT

コメントを残す

メールアドレスが公開されることはありません。 ※ が付いている欄は必須項目です

前の記事

通常処理とBiopythonの処理速度の比較

次の記事

アミノ酸配列から全コドンパターンの逆算