※実装時のパラメータ¶
PairHMMによるアラインメントを実装していく。なお、（kerasy.Bio.alignment.PairHMM）に定義してあります。
In [1]:
from kerasy.Bio.alignment import PairHMM
In [2]:
# サンプル配列
seqX = "CCAGAGCTGTGGCAGACAGTGGCT"
seqY = "CCAGCTGTGCAGACACTGGCTT"
In [3]:
model = PairHMM()
In [4]:
! cat params.json
In [5]:
model.set_params(path="params.json")
Pair HMM とその利点
- HMMのために開発されたアルゴリズムを使える。（らしい←正解ラベルの無いアラインメントでBaum-Welchアルゴリズムをどのように使うのかは現在考察中です。）
- アラインメントの確率的な解釈が可能になる。
グローバルアラインメントを出力する Pair-HMM¶
- Mは x と y がマッチした状態。（xi と yj がアライン）
- Xは x に挿入が起きている状態。（xi と − がアライン）
- Yは x に欠失が起きている状態。（− と yj がアライン）
Pair HMMのためのViterbiアルゴリズム¶
2本数の配列 x,y を与えると確率最大のアラインメントを求めるアルゴリズム。
Initialization:¶
{VM(0,0)=1,VM(i,0)=VM(0,j)=0 for i≠0,j≠0VX(0,0)=0,VX(1,0)=δqx1VX(i,0)=εqxiVX(i−1,0) for i=2,…,nVY(0,0)=0,VY(0,1)=δqy1VY(0,j)=εqyjVY(0,j−1) for j=2,…,mRecursion: i=1,…,n,j=1,…,m¶
VM(i,j)=pxiyjmax{(1−2δ−τ)VM(i−1,j−1)(1−ε−τ)VX(i−1,j−1)(1−ε−τ)VY(i−1,j−1)VX(i,j)=qximax{δVM(i−1,j)εVX(i−1,j)VY(i,j)=qyjmax{δVM(i,j−1)εVX(i,j−1)Termination:¶
ve=τmax(VM(n,m),VX(n,m),VY(n,m))In [6]:
P_optimal = model.align(seqX,seqY,return_score=True)
In [7]:
print(f"π*={P_optimal}")
アラインメントの信頼性¶
配列 x,y が与えられた時、アラインメント π が得られる条件付き確率を考えることができ、それがアラインメントの信頼度を表していると考ることができる。
P(π|x,y)=P(x,y,π)P(x,y)P(x,y)=∑alignments πP(x,y,π)※ なお、ここで P(x,y) は x,y に関する全ての可能なアラインメントについて確率を足し合わせたものであり、Forwardアルゴリズムで計算することができる。
Pair HMMのためのForwardアルゴリズム¶
※ Viterbiの max を sum に変更しただけ。
Initialization:¶
{FM(0,0)=1,FM(i,0)=FM(0,j)=0 for i≠0,j≠0FX(0,0)=0,FX(1,0)=δqx1FX(i,0)=εqxiFX(i−1,0) for i=2,…,nFY(0,0)=0,FY(0,1)=δqy1FY(0,j)=εqyjFY(0,j−1) for j=2,…,mRecursion: i=1,…,n,j=1,…,m¶
FM(i,j)=pxiyjsum{(1−2δ−τ)FM(i−1,j−1)(1−ε−τ)FX(i−1,j−1)(1−ε−τ)FY(i−1,j−1)FX(i,j)=qxisum{δFM(i−1,j)εFX(i−1,j)FY(i,j)=qyjsum{δFM(i,j−1)εFX(i,j−1)Termination:¶
Fe=τ sum(FM(n,m),FX(n,m),FY(n,m))In [8]:
F,P_all = model.forward(seqX,seqY)
In [9]:
print(f"𝐹^𝐸(𝑛,𝑚)=𝑃(𝑥,𝑦)=Σalignments𝜋𝑃(𝑥,𝑦,𝜋)={P_all}")
In [10]:
print(f"P(π*|x, y)={P_optimal/P_all}")
xi と yj がアラインされる確率¶
- 次の式で表される。 P(xi♢yi|x,y)=∑π∈ΩijP(x,y,π)P(x,y)=FM(i,j)⋅BM(i,j)P(x,y)
※ Ωij：xi と yj がアラインされるアラインメントの集合
Pair HMMのためのBackwardアルゴリズム¶
Needleman-Wunsh-Gotoh のBackwardアルゴリズムとの類似性に注目！
Initialization:¶
{BM(n,m)=BX(n,m)=BY(n,m)=τBX(i,m)=εqxi+1BX(i+1,m) for i=n−1,…,1BY(i,m)=εqyj+1BX(n,j+1) for j=m−1,…,1BM(i,m)=δqxi+1BX(i+1,m) for i=n−1,…,1BM(n,j)=δqyj+1BY(n,j+1) for j=m−1,…,1Recursion: i=n−1,…,1,j=m−1,…,1¶
BM(i,j)=sum{(1−2δ−τ)pxiyjBM(i+1,j+1)δqxi+1B(i+1,j)δqyj+1BY(i,j+1)BX(i,j)=(1−ε−τ)pxi+1,yj+1BM(i+1,j+1)+εqxi+1BX(i+1,j)BY(i,j)=(1−ε−τ)pxi+1,yj+1BM(i+1,j+1)+εqyj+1BY(i,j+1)In [11]:
B = model.backward(seqX,seqY)
In [12]:
Pij = (F*B)/P_all
In [13]:
import seaborn as sns
import matplotlib.pyplot as plt
In [14]:
sns.heatmap(Pij)
plt.show()
In [15]:
# 可視化のために pandas を使用。
import pandas as pd
r,c = Pij.shape
pd.set_option('display.max_rows', r)
pd.set_option('display.max_columns', c)
In [16]:
df = pd.DataFrame(Pij)
In [17]:
df
Out[17]:
In [ ]: