系統樹の仕組みを確かめるための初歩的なコード(python)

こんにちは。早雲です。

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

遺伝学とプログラミングに興味がある方はぜひ試してみてください。

・塩基配列から系統樹の作成

系統樹がどのように作られているのかを知りたいなと思い、プログラムを書いてみました。
有名なのは近隣接合法ですが、ひとまず一番簡単なUPGMA法を実装してみました。
ちなみに系統樹を作る方法ですが、各塩基配列の距離行列を算出して、それを階層的クラスタリングする、と考えるとわかりやすかもしれません。
では、コードです。
コードはpython3.7で挙動を確認しました。

import math
import itertools
import pandas as pd
import numpy as np


#Jules-Cantor genetic distance

def jules_cantor(sequence_a,sequence_b):    
   sequence_a_gap_del = ""
   sequence_b_gap_del = ""
   for a,b in zip(sequence_a,sequence_b):
       if "-" in a:
           continue
       elif "-" in b:
           continue
       else:
           sequence_a_gap_del = sequence_a_gap_del + a
           sequence_b_gap_del = sequence_b_gap_del + b
   
   q_counter=0
   for a_gap_del, b_gap_del in zip(sequence_a_gap_del,sequence_b_gap_del):
       if a_gap_del == b_gap_del:
           q_counter+=1
   p = 1 - q_counter/len(sequence_a_gap_del)
   d = -(3/4)*(math.log(1-((4*p)/3)))
   return d

def UPGMA(seq_dic):
   combinations_list = list(itertools.combinations(list(seq_dic.keys()), 2))
   dist_mat_dic = {}
   
   for i in combinations_list:
       dist_mat_dic[str(i[0])+","+str(i[1])]=jules_cantor(seq_dic[i[0]],seq_dic[i[1]])
   
   dataframe_index_colonm_tmp = []
   for key,value in dist_mat_dic.items():
       tmp = key.split(",")
       dataframe_index_colonm_tmp.append(tmp)
   dataframe_index_colonm_list = [e for x in dataframe_index_colonm_tmp for e in x]
   dataframe_index_colonm = list(set(dataframe_index_colonm_list))
   dataframe_index_colonm.sort()
   
   dataframe_list=[]
   for e in dataframe_index_colonm:
       tmp_dic = {}
       for keys,values in dist_mat_dic.items():
           tmp_list = keys.split(",")
           if e in tmp_list:
               tmp_list.remove(e)
               tmp_dic[tmp_list[0]] = values
       row = []
       for c in dataframe_index_colonm:
           try:
               row.append(tmp_dic[c])
           except KeyError:
               row.append(np.nan)
       dataframe_list.append(row)
   
   df = pd.DataFrame(dataframe_list, index = dataframe_index_colonm, columns = dataframe_index_colonm)        
   
   #最小値をとる
   num_dic = {}
   while len(df.columns.values) > 1:
       min_list = ["first",100]
       new_columns_list = list(df.columns.values)
       for e in new_columns_list:
           tmp_list = [df[e].idxmin(),df[e].min()]
           if tmp_list[1] < min_list[1]:
               min_list = tmp_list
               final_list = [e,min_list[0],min_list[1]]
       try:
           pre_dis_0 = num_dic[final_list[0]]
       except KeyError:
           pre_dis_0 = 0
       try:
           pre_dis_1 = num_dic[final_list[1]]
       except KeyError:
           pre_dis_1 = 0
       dis_0 = (final_list[2]/2) - pre_dis_0
       dis_1 = (final_list[2]/2) - pre_dis_1
       df_l = (df[final_list[0]]+df[final_list[1]])/2
       df_1 = df.drop(final_list[0],axis = 1)
       df_2 = df_1.drop(final_list[1],axis = 1)
   
       df_2["("+final_list[0]+":"+str(dis_0)+","+final_list[1]+":"+str(dis_1)+")"]=df_l
               
       df_m = (df_2.loc[final_list[0]]+df_2.loc[final_list[1]])/2
       df_3 = df_2.drop(final_list[0],axis = 0)
       df_4 = df_3.drop(final_list[1],axis = 0)
       df_4.loc["("+final_list[0]+":"+str(dis_0)+","+final_list[1]+":"+str(dis_1)+")"]=df_m
       df = df_4
       num_dic["("+final_list[0]+":"+str(dis_0)+","+final_list[1]+":"+str(dis_1)+")"] = final_list[2]/2
       final_list=[]
       min_list = ["first",100]
   return str(df.columns.values[0])+";"

seq_dic = {"sequence_a":"AA-ACATTT", "sequence_b":"AACCC-TCC", "sequence_c":"AT-AC-TAT", "sequence_d":"AAAAC-TAT"}
newick = UPGMA(seq_dic)
print(newick)

長いですね。。。興味がある方は確認してみてください。

それと、『生物学に特化した短編小説とエッセイ』では生物学を題材にした小説やエッセイを書いています。

興味がある方はぜひ、読んでください。
ではまた。