📂アルゴリズム

スミス-ウォーターマン アラインメント：ローカル シーケンス アラインメント

概要

二つの塩基配列から、最も似ている部分のアラインメントを見つけることを局所アラインメント^{local Alignment}と言い、その方法として最も広く使われているのがスミス-ウォーターマンアルゴリズム^{smith-Waterman Algorithm}を紹介する。

配列アラインメントにはあまりにも多くのケースがあるため、動的計画法を通じて効率的かつ迅速に計算する必要がある。

アルゴリズム¹

Input

二つの文字列 $\textbf{v}, \textbf{w}$ と代替行列 $\delta$ が与えられているとして、次のように表示されているとしよう。

$$\textbf{v} := \left( v_{1} , \cdots , v_{m} \right) \\ \textbf{w} := \left( w_{1} , \cdots , w_{n} \right)$$

$\textbf{v}$ の長さは $m$、$\textbf{w}$ の長さは $n$ だ。二つのキャラクター $x,y$ のペナルティは $\delta (x,y)$ のように表され、特にギャップを表すキャラクターは $-$ と表すことにしよう。

---

Step 1. 初期化

空の行列 $\left( s_{i,j} \right)_{i=0, \cdots, m \\ j = 0, \cdots , n}$ を作り、次のように最初の行、最初の列を初期化する。

$$s_{i,0} \gets 0 \\ s_{0,j} \gets 0$$

---

Step 2. 動的計画法

$s_{i-1,j}$ と $s_{i,j-1}$ がすでに計算されていれば、$s_{i,j}$ を次のように計算する。（もちろん、$s_{i-1,j-1}$ は計算されているだろう。）

$$s_{i,j} \gets \max \begin{cases} s_{i-1,j} + \delta \left( v_{i}, - \right) \\ s_{i,j-1} + \delta \left( -, w_{j} \right) \\ s_{i-1, j-1} + \delta \left( v_{i} , w_{j} \right) \\ 0 \end{cases}$$

---

Ouput

行列で最もスコアが高い要素を $s_{\text{max}}$ としよう。$s_{\text{max}}$ から始めて左、左上、上の要素から最もスコアが高い要素を選び、0に遭遇するまで止まる。二つの文字列の一部と最も似ている局所アラインメントを得る。そのスコアは $s_{\text{max}}$ だ。

例

Input

julia> s1 = "ATTAGCTCAC";

julia> s2 = "ACTATAGCGCAG";

上のように二つの塩基配列をグローバルアラインメントしよう。

julia> EDNAFULL
SubstitutionMatrix{BioSymbols.DNA,Int64}:
     A  C  M  G  R  S  V  T  W  Y  H  K  D  B  N
  A  5 -4  1 -4  1 -4 -1 -4  1 -4 -1 -4 -1 -4 -2
  C -4  5  1 -4 -4  1 -1 -4 -4  1 -1 -4 -4 -1 -2
  M  1  1 -1 -4 -2 -2 -1 -4 -2 -2 -1 -4 -3 -3 -1
  G -4 -4 -4  5  1  1 -1 -4 -4 -4 -4  1 -1 -1 -2
  R  1 -4 -2  1 -1 -2 -1 -4 -2 -4 -3 -2 -1 -3 -1
  S -4  1 -2  1 -2 -1 -1 -4 -4 -2 -3 -2 -3 -1 -1
  V -1 -1 -1 -1 -1 -1 -1 -4 -3 -3 -2 -3 -2 -2 -1
  T -4 -4 -4 -4 -4 -4 -4  5  1  1 -1  1 -1 -1 -2
  W  1 -4 -2 -4 -2 -4 -3  1 -1 -2 -1 -2 -1 -3 -1
  Y -4  1 -2 -4 -4 -2 -3  1 -2 -1 -1 -2 -3 -1 -1
  H -1 -1 -1 -4 -3 -3 -2 -1 -1 -1 -1 -3 -2 -2 -1
  K -4 -4 -4  1 -2 -2 -3  1 -2 -2 -3 -1 -1 -1 -1
  D -1 -4 -3 -1 -1 -3 -2 -1 -1 -3 -2 -1 -1 -2 -1
  B -4 -1 -3 -1 -3 -1 -2 -1 -3 -1 -2 -1 -2 -1 -1
  N -2 -2 -1 -2 -1 -1 -1 -2 -1 -1 -1 -1 -1 -1 -1
(underlined values are default ones)

代替行列のギャップペナルティは、線形に $-5$ 点としよう。すると、使われるキャラクターは四つの主要な塩基だけで、キャラクターが一致するとき$5$点をもらい、間違ったとき$-4$点、ギャップペナルティで$-5$をもらう。

理解を助けるため、エクセルで簡単に実装してみよう。

Step 1. クエリシーケンスを行に、リファレンスシーケンスを列にして初期化する。

---

Step 2. 各セルでスコアを計算する。

C3（3行3列）から次のように関数を入力して、全てのマスにコピーすればいい。

=MAX(B2+IF(C$1=$A3,5,-4), C2-5,B3-5,0)

---

Output. 右下からトレースバックする。

右下の隅から次のルールに従って逆に登っていく。

• 赤 ↖: 行と列がマッチしていないので、どちらかを選んでも構わない。
• オレンジ ↖: 行と列がマッチしているので、その塩基を選ぶ。
• 灰色 ← もしくは ↑: ギャップ、インデル^indelを意味する。上と左のうち、相対的にスコアが高い方を選ぶ。

コード

動的計画法として、例で見たようにエクセルで簡単に実装することもできる。勉強目的ならいいが、これをずっと再利用するのは明らかに愚かなことだ。次はジュリアで同じ作業を簡単かつ迅速に実行するコードだ。

using BioAlignments

s1 = "ATTAGCTCAC";
s2 = "ACTATAGCGCAG";

scoremodel = AffineGapScoreModel(EDNAFULL, gap_open=-5, gap_extend=0);
local_result = pairalign(LocalAlignment(), s1, s2, scoremodel)

実行結果

julia> using BioAlignments

julia> s1 = "ATTAGCTCAC";

julia> s2 = "ACTATAGCGCAG";

julia> scoremodel = AffineGapScoreModel(EDNAFULL, gap_open=-5, gap_extend=0);

julia> local_result = pairalign(LocalAlignment(), s1, s2, scoremodel)
PairwiseAlignmentResult{Int64,String,String}:
  score: 26
  seq:  1 ATTAGCTCA  9
          || ||| ||
  ref:  4 AT-AGCGCA 11