📂アルゴリズム

ニードルマン・ワンシュアルゴリズム：グローバルシークエンスアラインメント

概要

二つの塩基配列の共通部分が最も多くなるように整列を見つけることをグローバルアライメントと呼び、その方法として最も広く使われているのがニードルマン-ウンシュアルゴリズムである。

配列アライメントにはあまりにも多くのケースがあるため、ダイナミックプログラミングを通じて効率的かつ迅速に計算する必要がある。

アルゴリズム¹

入力

二つの文字列$\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)$のように表され、特にギャップを表すキャラクターは$-$として表すことにする。

---

ステップ1. 初期化

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

$$ s_{i,0} \gets i \cdot \delta \left( v_{i},- \right) \\ s_{0,j} \gets j \cdot \delta \left( -, w_{j} \right) $$

---

ステップ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) \end{cases}$$

---

出力

$s_{mn}$から行と列が右下から左上に向かって一致する要素を選ぶことで、二つの文字列と最も類似したグローバルアライメントを得ることができる。そのスコアは$s_{mn}$である。

例

入力

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$を受ける。

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

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

---

ステップ2. 各セルでスコアを計算する。

C3（3行3列）に次の関数を入力して、すべてのセルにコピーすればよい。

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

---

出力. 右下からトレースする。

右下の隅から次のルールに従って逆方向にたどる。

• 赤 ↖: 行と列がマッチしていないので、どちらかを選んでも構わない。
• オレンジ ↖: 行と列がマッチしているので、その塩基を選択する。
• グレー ← または ↑: ギャップ、つまりインデルを意味する。上と左のどちらか得点が高い方を選ぶ。

トレース過程でわかるように、最適な配列アライメントは必ずしもユニークではないかもしれない。「A-T-TAGCCAC」、「A-T-TAGCGAC」、さらには「ACT-TAGCGAC」も可能だが、見つかった可能な配列アライメントのスコアは全て22点で同じである。

コード

ダイナミックプログラミングであるため、上記のようにエクセルでも簡単に実装できる。例で見たように、勉強目的であれば問題ないが、これをずっと再利用することは明らかに無駄である。次はジュリアで同じ作業を迅速かつ簡単に行うコードだ。

using BioAlignments

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

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

実行結果

using BioAlignments

julia> s1 = "ATTAGCTCAC";

julia> s2 = "ACTATAGCGCAG";

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

julia> global_result = pairalign(GlobalAlignment(), s1, s2, scoremodel)
PairwiseAlignmentResult{Int64,String,String}:
  score: 22
  seq:  1 A-T-TAGCTCAC 10
          | | |||| || 
  ref:  1 ACTATAGCGCAG 12