今回はNGSによって得られた無数の重複DNA断片から一つの塩基配列を決定するアルゴリズムについて考えます。
全域探索をした際、最悪の場合O(N^2)になる可能性があるため、あらかじめDNA断片をソートして組み立てを開始することにします。
目次
アルゴリズム
「DNA断片を組み立てて一つの配列にする」ということをアルゴリズム風に言い換えると、「全断片の中で共有部分が最長である２配列を組み合わせていき最終的に一つの長い配列を作る」ということになります。
1.DNA断片のソート
組み立て前のDNA断片を辞書順にソートしたものを、fragmentsに代入します。
fragments = sorted(['ATTAGACCTG','CCTGCCGGAA','AGACCTGCCG','GCCGGAATAC'])2.最長共有配列の探索
最長共有配列(Longest common substring)を探す関数lcsを定義します。
def lcs(seq1,seq2):
t_1 = ''
t_2 = ''
for i in range(len(seq1)):
if seq1[:i] == seq2[-i:]:
t_1 = seq1[:i]
for i in range(len(seq2)):
if seq2[:i] == seq1[-i:]:
t_2 = seq2[:i]
if len(s_1) > len(s_2):
return t_1
elif len(s_1) < len(s_2):
return t_2
elif len(s_1) == len(s_2):
return t_1t_1はseq1の前方とseq2の後方がオーバーラップした場合、t_2はseq2の前方とseq1の後方がオーバーラップした場合の共通配列で、lcsはこのうち長い方の配列を返します。
3.配列加工
s_nに２配列の共有部分の重複を取り除いた配列を代入します。
この操作は後で多用するため、関数seq_modifyとして定義します。
def seq_modify(s_1,s_2,seq_overlapped):
if seq_overlapped == s_1[:len(seq_overlapped)]:
s_n = s_2 + s_1[len(seq_overlapped):]
elif seq_overlapped == s_2[:len(seq_overlapped)]:
s_n = s_1 + s_2[len(seq_overlapped):]
return s_nここからがアルゴリズムの本番ですが、リストの中身が１つになるまで以下の操作を繰り返します。
import time
start = time.time()
while len(fragments) !=1:
lcs_len = [0]
for i in range(len(fragments)):
t_i = fragments[i]
for j in range(len(fragments)):
t_j = fragments[j]
if j!=i and len(lcs(t_i,t_j))>lcs_len[-1]:
lcs_len.append(len(lcs(t_i,t_j)))
s_1 = t_i
s_2 = t_j
else:
pass
fragments.remove(s_1)
fragments.remove(s_2)
fragments.append(seq_modify(s_1,s_2,lcs(s_1,s_2)))
print(fragments)
print(fragments[0])
process_time = time.time() - start
print(f'{process_time}s')lcs_lenに最長共有部分の長さを代入していき、この長さを超える新たな共有部分が出てきた場合のみにリストの更新を行えば一回分の全域探索を終えると最長共有部分を取得していることになります。
このループが終了すると、fragmentsには一つの配列のみが格納されています。これが今回の目的である塩基配列です。
4.検証
STEP.3まででほとんど終わりのようなものですが、念の為出来上がった配列に全てのDNA断片が含まれているか検証します。
sequence にfragments[0]を代入します。
sequence = fragments[0]DNA断片を格納したfragmentsには破壊的処理を施したため、fragments_に新たにDNA断片を格納したリストを代入します。
#検証用に変数mismatchesをカウンターとします。
mismatches = 0
for seq in fragments_:
if seq not in fragments_:
mismatches +=1
else:
pass
if mismatches != 0:
print('Assembly failed.')
elif mismatches == 0:
print('Assembly succeeded.')アルゴリズムを実行しターミナルに”Assembly succeeded.”と表示されれば、組み立ての成功です。
