塩基配列のアライメントの仕組みを知るためのコード(python)

こんにちは。早雲です。

この記事はすこし専門的な内容になっており、『生物学に特化した短編小説とエッセイ』とは別のマガジンに組み込まれています。

遺伝学とプログラミングに興味がある方はぜひ試してみてください。すこし不親切な記事ですが……

塩基配列のアライメントの仕組みを知るためのコード

勉強がてら、遺伝子のアライメントの簡単なコードを書いてみましたので共有します。初歩的ですが、どのような挙動でアライメントが行われるか知るのには役に立つと思います。
コードはpython3.7で実行し確認しました。

sequence_s="AGCT"
sequence_t="ACGCT"

sequence_s_list = [x for x in sequence_s]
sequence_t_list = [y for y in sequence_t]


m=len(sequence_s)
n=len(sequence_t)
d=-2

def comparison_site(si,tj):
   if si == tj:
       return 1
   else:
       return -3
       
score_mat = {}
alin_mat = {}

for i in range(m+1):
   point_i=str(i)+",0"
   score_mat[point_i] = i*d
   alin_mat[point_i] = "lateral"
for j in range(n+1):
   point_j="0,"+str(j)
   score_mat[point_j] = j*d
   alin_mat[point_j] = "down"

for i in range(m):
   li = i+1
   for j in range(n):
       lj = j+1
       pre_score_a = score_mat[str(li-1)+","+str(lj-1)]
       pre_score_b = score_mat[str(li-1)+","+str(lj)]
       pre_score_c = score_mat[str(li)+","+str(lj-1)]
       score_candidate_list = [pre_score_a + comparison_site(sequence_s_list[i],sequence_t_list[j]), pre_score_b + d, pre_score_c + d]         
       score = max(score_candidate_list)
       score_chosed_index = score_candidate_list.index(max(score_candidate_list))
       if score_chosed_index == 0:
           alin_direction = "Diagonally" 
       elif score_chosed_index == 1:
           alin_direction = "lateral"
       elif score_chosed_index == 2:
           alin_direction = "down"
       score_mat[str(li)+","+str(lj)] = score
       alin_mat[str(li)+","+str(lj)] = alin_direction
alin_mat['0,0'] = "non"

M = m
N = n
sequence_s_alin = ""
sequence_t_alin = ""
flag = True

while flag:
   trace_direction = alin_mat[str(M)+","+str(N)]
   if trace_direction == "Diagonally" :
       sequence_s_alin = sequence_s_list[M-1] + sequence_s_alin
       sequence_t_alin = sequence_t_list[N-1] + sequence_t_alin
       M = M - 1
       N = N - 1
   elif trace_direction == "lateral":
       sequence_t_alin = "-" + sequence_t_alin 
       sequence_s_alin = sequence_s_list[M-1] + sequence_s_alin
       M = M - 1
           
   elif trace_direction == "down":            
       sequence_s_alin = "-" + sequence_s_alin              
       N = N - 1
       sequence_t_alin = sequence_t_list[N-1] + sequence_t_alin 
   elif trace_direction == "non":
       flag = False
       
print("sequence_s : ",sequence_s_alin)
print("sequence_t : ",sequence_t_alin)
print("score : ",score_mat[str(m)+","+str(n)])

うーん。なんかもっとスマートな書き方がありそうなんですが、とりあえずこんな感じで。

また、『生物学に特化した短編小説とエッセイ』もよろしくお願いいたします。

ではまた。

2022/10/26内容修正

掲載したコードのインデントがなぜか消えているのに気が付きましたので、修正します。大変失礼しました。誰か参考にしている方がいるかもしれないので。。

それにしても、このコードを書いて二年くらい経って見直しているのですが、さすがにひでえな。。関数全然使わないじゃん。。モジュール化のモの字もないじゃん。。と思いました。

上で「なんかもっとスマートな書き方がありそうなんですが、とりあえずこんな感じで」って言ってますが、「とりあえずこんな感じで」じゃねえよ。。