制限酵素の切断領域(Restriction Sites)を求めるアルゴリズム

一般的に制限酵素の切断領域は回文構造（パリンドローム）になっています。

今回はその回文領域を求めるためのアルゴリズムについて考えていきます。

アルゴリズム

切断領域の定義を4-12のヌクレオチドから構成される回文領域とします。

相補鎖および逆相補鎖の生成

初めに、対象配列の相補鎖および逆相補鎖を用意しておきます。

逆相補鎖は回文構造の探索用、相補鎖は回文の検証に使用します。

from Bio.Seq import Seq
seq = Seq('''GCGCTATAGTCCCCTCCAGTCCCAAACTCGAGGTAGGAGGGTGTGCTAGATTTTGTTTTA
GAAGTAAAGGTCTCGTCGGCAGCTTGTTTGGAGAGCTCCGCTAAATAAGAATCCACTTAA
CAAAGACGCCACACACAAGTACCTCGTACGTAGAAGCACTGCGGAAATACTAAGCTCCAA
TCGCGTTTGAGCAGTCCTCTTCACTGACTAGTTTTCCAGACGGGCCTAGCGCCGCGGCAA
GCCGGTACTATGTGATCACCCCGGTAAGGTGGGCCTTTGCAGGCCAGCATTCGCTCAAGA
ACCGCATCAACGCCTGCACTGCCGTTGAGATTGGATTTCGGATGATCGACGCGAGTATAG
TATCGTATCGTAGGCTGAAGTATCGAGCCATTGAGCTCGAAAGGAGTAACTATATCCCTA
ATGCCCACTAGGCTTAGTTAAAGGACGTCCTTGCGAGGACCTGAGGTGGATGGCATTCTT
TCATGCACACGGTCGGTACACGTTGCTTTGGCCCGCTGAGGTAGCATTATATTGTCCTGG
CAAACGAAGGCACCGTTCATTCTGGTTATCGGACATTATTGCTGTGTCTCCATTTAGCTG
GGATGGAGATCTCACGGAATTACGCCGAACTCTGTTTAAGAAACCACTTATCCTCTATGG
GAATCTCTCGTTTGTCTTGCGCTCCAGCCCGTGGCTGTTAATTGTAGCCGTGCTACGGTA
TCAATTCTCACTGATACCTTCAAGGCGTATACTGGCTAGTGGAAATGGTCGGGCGGTGCC
CATGCTCAATACTTTTCACTTTAGCATATCTATTCTCGAGATATCCGCAACAGCAAACAG
CAACGCGTGCTAAGGCATGTTGATAATGTATAGGGCGGTCGTGCCATGCGCG'''.replace('\n',''))

BiopythonライブラリのSeqモジュールを使用して瞬時に変換します。

seq_c = seq.complement()
seq_r = seq.reverse_complement()

Biopythonのインストールと簡単な塩基配列操作

回文構造の探索

次に、逆相補鎖を利用して、条件を満たすパリンドロームの候補となる配列を探索していきます。

この時、候補となる配列の長さが偶数であることに注意します。

s_l= []
for i in range(len(seq)):
    s_i = seq_r[i:]
    if len(s_i) >= 12:
        for j in range(4,13):
            s = s_i[:j]
            if s in seq and s not in s_l and len(s)%2 ==0 :
                s_l.append(s)
    elif len(s_i) < 12:
        for j in range(4,len(s_i)+1):
            s = s_i[:j]
            if s in seq and s not in s_l and len(s)%2 ==0 :
                s_l.append(s)

探索結果をs_lに代入しました。

回文配列の検証

通常ならs_l中に入っているもののインデックスを片っぱしからfindメソッドで取り出していけばいいのですが、findメソッド及びindexメソッドには最初にマッチしたものを取り出したら仕事を辞めてしまうという欠点があります。

そのため、s_l中の全ての断片が必ずしも回文となると言うことはできません。

そこで、以下のように回文の検証を行ったものを[開始地点,長さ]と言う形で取り出してみたいと思います。

t_l = [[t+1,len(s)] for s in s_l for t in range(len(seq)) if seq[t:t+len(s)] == s and seq_c[t:t+len(s)] == s[::-1]]

結果

回文配列の結果をターミナルに出力します。

for i in sorted(t_l):
    print(str(i[0])+' '+str(i[1]))

実行結果

意外と多いですね。

計算量

見ての通り、二重ループを使用しているため計算量はO(N^2)となることは確実ですが、上記の処理を関数化して文字列数による実行時間の推移をグラフにしてみます。

from Bio.Seq import Seq
import time 
def get_locations_of_restriction_sites(seq):
    seq = Seq(seq.replace('\n',''))
    seq_c = seq.complement()
    seq_r = seq.reverse_complement()
    s_l= []
    for i in range(len(seq)):
        s_i = seq_r[i:]
        if len(s_i) >= 12:
            for j in range(4,13):
                s = s_i[:j]
                if s in seq and s not in s_l and len(s)%2 ==0 :
                    s_l.append(s)
        elif len(s_i) < 12:
            for j in range(4,len(s_i)+1):
                s = s_i[:j]
                if s in seq and s not in s_l and len(s)%2 ==0 :
                    s_l.append(s)
    return sorted([[t+1,len(s)] for s in s_l for t in range(len(seq)) if seq[t:t+len(s)] == s and seq_c[t:t+len(s)] == s[::-1]])

意外とコンパクトにまとまりました。

ここから、推移を確認するためにランダムに塩基配列を増やしていきながら経過を見ていきます。

import random
from matplotlib import pyplot as plt 
from tqdm import tqdm
if __name__ == '__main__':
    fig = plt.figure()
    plt.grid()
    plt.xlabel('bp')
    plt.ylabel('time(s)')
    nucleotides = ['A','C','G','T']
    seq = ''
    k = 10
    for i in tqdm(range(100)):
        s = time.time()
        seq += ''.join(random.choices(nucleotides,k=k))
        get_locations_of_restriction_sites(seq)
        plt.scatter(k*i,time.time()-s,c='blue')
    fig.savefig('a.jpg')

出来上がった画像がこちらです。