- Nussinov Algorithmは、最もシンプルな二次構造予測アルゴリズムである。
- 「最大いくつの塩基対を作ることができるか」を調べる。
- 動的計画法を用いて効率的に調べるために、以下のようなPseudoknotは無視している。様々な重要なRNAが、Pseudoknotの形をとることが確認されているが、相対的にこのような構造をとる確率は低いので、計算の効率を優先し、これらを無視する。
Algorithm¶
- ※「$i$ から $j$ までの部分塩基配列で最大幾つの塩基対が作れるか」を $\gamma(i,j)$ が保持している。
- 今回の実装では $j-i-1<3$ の時は、$i,j$ が塩基対を作れないものとしている。
- ポイントはRecursionの最後の式で、「$i$ から $k$ までの部分塩基配列」と「$k+1$ から $j$ までの部分塩基配列」を足し合わせる操作を行う。（Bifurcation）
Initialization:¶
$$\begin{aligned} \gamma(i,i)&=0&\text{ for }i=1\text{ to }L\\ \gamma(i,i-1)&=0&\text{ for }i=2\text{ to }L \end{aligned}$$Recursion:¶
$$ \gamma(i,j) = \max \begin{cases} \gamma(i+1,j)\\ \gamma(i,j-1)\\ \gamma(i+1,j-1)+\delta(i,j)\\ \max_{i\leq k\verb|<|j}\left[\gamma(i,k) + \gamma(k+1,j)\right] \end{cases} $$TraceBack¶
基本的なアルゴリズムは同じですが、Bifurcationが起きた際にトレースバックポインターを複数用意する必要があることには注意が必要です。
bp = []
stack = [(0,N-1)]
while(stack):
i,j = stack.pop(0)
delta = 1 if self._is_bp(sequence[i],sequence[j]) and (j-i>3) else 0
if (i>=j): continue
elif gamma[i+1][j] == gamma[i][j]: stack.append((i+1,j))
elif gamma[i][j-1] == gamma[i][j]: stack.append((i,j-1))
elif gamma[i+1][j-1]+delta == gamma[i][j]: bp.append((i,j)); stack.append((i+1,j-1))
else:
# Bifurcation
for k in range(i,j):
if gamma[i][k] + gamma[k+1][j] == gamma[i][j]:
stack.append((k+1,j))
stack.append((i,k))
break
自作のモジュール（kerasy.Bio.structure.Nussinov）を使って実装します。
In [1]:
from kerasy.Bio.structure import Nussinov
In [2]:
model = Nussinov()
In [3]:
# 以下の形のパラメタファイルを受け取ります。
! cat params.json
In [4]:
model.load_params("params.json")
In [5]:
model.params()
In [6]:
# サンプル配列
sequence = "GCGGAAACGGGGTCA"
In [7]:
# 最大ペアの配列
model.predict(sequence, memorize=True, traceback=True)
In [8]:
gamma = model.gamma
In [9]:
# DP行列は以下のようになっています。
model._printAsTerai(gamma, sequence)
$i,j$ が塩基対を組むときの塩基対最大化¶
塩基 $i,j$ が塩基対を組むという制約がある場合を考えます。これは、「$i$ から $j$ までの部分塩基配列を除いた残りで最大幾つの塩基対が作れるか」</b>を保持した $\omega(i,j)$ と $\gamma(i,j)$ を組み合わせることで効率的に求めることができます。
Algorithm¶
- 「$i$ から $j$ までの部分塩基配列で最大幾つの塩基対が作れるか」を $\gamma(i,j)$ が保持している。
- 「$i$ から $j$ までの部分塩基配列を除いた残りで最大幾つの塩基対が作れるか」を $\omega(i,j)$ が保持している。
Initialization:¶
$$\omega(i,L) = 0$$Recursion:¶
$$ \omega(i,j) = \max \begin{cases} \omega(i-1,j)\\ \omega(i,j+1)\\ \omega(i-1,j+1)+\delta(i-1,j+1)\\ \max_{1\leq k\verb|<|i}\left[\omega(k,j) + \gamma(k,i-1)\right]\\ \max_{j\verb|<| k\leq L}\left[\omega(j+1,k) + \gamma(i,k)\right] \end{cases} $$In [10]:
model.outside(sequence, memorize=True)
In [11]:
omega = model.omega
In [12]:
# DP行列は以下のようになっています。
model._printAsTerai(omega, sequence)
$Z(i,j)$ の計算¶
| variable | definition |
|---|---|
| $\gamma\left(i,j\right)$ | the maximum number of base-pairs for subsequence from $i$ to $j$. |
| $\omega\left(i,j\right)$ | the maximum number of base-pairs excluding subsequence from $i$ to $j$. |
上記の定義より、$\gamma,\omega$ を用いて $Z(i,j)$ が以下のように求まることがわかります。
$$Z(i,j) = \begin{cases} \gamma\left(i+1,j-1\right) + 1 + \omega\left(i,j\right) & (\text{if }i\text{ and }j\text{-th nucleotides can form a base-pair})\\ 0 & (\text{otherwise}) \end{cases}$$In [13]:
model.ConstrainedMaximize(sequence, memorize=False)
備忘録¶
Nussinov Algorithmでは、これまでの動的計画法とは異なり、 $i,j$ を斜め方向に舐める必要がありました。
そこで、その際に使用したコードを備忘録として残しておきます。
In [14]:
import numpy as np
In [15]:
N=4
In [16]:
# Nussinov Algorithm (γ)
idx=1
arr = np.zeros(shape=(N,N), dtype=int)
for ini_i in reversed(range(N)):
diff=N-ini_i
for i in reversed(range(ini_i)):
arr[i][i+diff] = idx
idx+=1
print(arr)
In [17]:
# Nussinov outsie (ω)
idx=1
arr = np.zeros(shape=(N,N), dtype=int)
for j in reversed(range(N)):
for i in range(N-j):
arr[i][j+i] = idx
idx+=1
print(arr)
In [ ]: