Aritalab:Lecture/Bioinformatics/Alignment

• このページで説明するアルゴリズムの Java コード

Needleman-Wunsch (-Sellers) アルゴリズム

1970年、分子生物学者の Saul B. NeedlemanとChristian D. Wunsch は、大域アライメントのアルゴリズムを Journal of Molecular Biology誌 (1970) 48: 443-453 に発表しました。いまでは誤って Needleman-Wunsch アルゴリズムと呼ばれるようになりましたが、論文は効率が悪い手法、O(n³) 時間、を紹介しています。ここで紹介されるアルゴリズムの原型を作ったのは Peter Sellers (1974) です。したがって正確には Needleman-Wunsch-Sellers アルゴリズムと呼ぶべきでしょう 。

基本は最長共通部分列(LCS)アルゴリズムと同じで、動的計画法によるLCSのページで紹介した以下の漸化式に基づいてスコア行列を計算します。ギャップとミスマッチに対してそれぞれペナルティスコア (gap penalty σ, mismatch penalty μ)が与えられています。

$D_{i,j} = max \begin{cases} D_{i-1,j} - \sigma \\ D_{i,j-1} - \sigma \\ D_{i-1,j-1} - \mu & \mbox{if } s_i \ne t_j \\ D_{i-1,j-1} +1 & \mbox{if } s_i = t_j \end{cases}$

ここではギャップペナルティを 2、ミスマッチペナルティを 1 として、２つの配列 gctagg と aattgaagg の間のアライメントを考えてみましょう。スコア行列は右のようになります。 1行目と1列目は 0, -2, -4, -6, ... -18 とギャップペナルティに基づく等差数列ができています。これは左隣（または上）のセルからギャップペナルティ 2 を順次マイナスしてできる部分です。大域アライメントでは配列の最初と最後を揃えるために、初期状態としてギャップの数に比例するペナルティを与えます。 残りのセルは動的計画法に基づいてスコアが埋められ、最適アライメントは以下のようになります。 1:aattgaagg:9 | | || 1:gct--a-gg:6 右図のトレースバックと見比べてみてください。大域アライメントのスコアは マッチ 4 × 1 + ミスマッチ 2 × (-1) + ギャップ 3 × (-2) = -4 となっています。大域アライメントにおけるギャップの数は配列長の差になります。具体的なアルゴリズムとJavaコードはこのページを参照してください。

g c t a g g 0 -2 -4 -6 -8 -10 -12 a -2 -1 -3 -5 -5 -7 -9 a -4 -3 -2 -4 -4 -6 -8 t -6 -5 -4 -1 -3 -5 -7 t -8 -7 -6 -3 -2 -4 -6 g -10 -7 -8 -5 -4 -1 -3 a -12 -9 -8 -7 -4 -3 -2 a -14 -11 -10 -9 -6 -5 -4 g -16 -13 -12 -11 -8 -5 -4 g -18 -15 -14 -13 -10 -7 -4

Smith-Waterman アルゴリズム

バイオインフォマティクスの草分けである Temple F. Smith と Michael S. Waterman が 1981 年に Journal of Molecular Biology誌 (147: 195–197) に発表した局所アライメントを求めるアルゴリズムです。MS Waterman は動的計画法を RNA 二次構造の解析などに応用して活躍した数学者ですが、TF Smithは配列アライメントに関する実績もなく、このアルゴリズムをWalther Goad から聞いて論文にしたといわれる曰くつきのアルゴリズムです。FASTAパッケージ（後述）の中では SSEARCH の名前で実装されていました。局所アライメントは以下のように定義されます。

局所アライメント = 与えられた２つの配列それぞれの全ての部分配列のペア中で、大域アライメントのスコアを最大にするもの

この定義だとすべての部分配列を列挙して大域アライメントを施す必要が思えますが、そうではありません。動的計画法において常にスコア 0 を候補に置くだけで、同じ計算時間で局所アライメントを求められるところが大変エレガントなアルゴリズムです。オリジナル論文はアフィンギャップ関数（後述）という、生物学的に意味のあるギャップペナルティを採用しています。

$D_{i,j} = max \begin{cases}0 \\ D_{i-1,j} - \sigma \\ D_{i,j-1} - \sigma \\ D_{i-1,j-1} - \mu & \mbox{if } s_i \ne t_j \\ D_{i-1,j-1} +1 & \mbox{if } s_i = t_j \end{cases}$

ここではギャップペナルティを 2、ミスマッチペナルティを 1 として、２つの配列 gctagg と aattgaagg の間のアライメントを考えてみましょう。スコア行列は右のようになります。 どこからアライメントを初めてもよいため 1 行目と 1 列目はすべて 0 にします。 残りのセルは動的計画法に基づいてスコアが埋められ、最適アライメントは以下のようになります。この例は説明の都合上、大域アライメントと同じ配列のため、たまたま右下部分が最大値になっています。アルゴリズム自体は最適値がスコア行列の任意の箇所に出てもよいことに注意してください。 7:agg:9 ||| 4:agg:6 具体的なアルゴリズムとJavaコードはこのページを参照してください。

g c t a g g 0 0 0 0 0 0 0 a 0 0 0 0 1 0 0 a 0 0 0 0 1 0 0 t 0 0 0 1 0 0 0 t 0 0 0 1 0 0 0 g 0 1 0 0 0 1 1 a 0 0 0 0 1 0 0 a 0 0 0 0 1 0 0 g 0 1 0 0 0 2 1 g 0 1 0 0 0 1 3

アフィンギャップと Gotoh アルゴリズム

上記のアルゴリズムでは、ギャップとして欠失する塩基数に比例したペナルティが課されますが、配列進化の観点では不適切なスコアです。 ギャップは配列が切断されることを意味しますが、その間に挿入される配列の長さは配列の切断というイベントに比較して重要ではありません。そのため、実用のアルゴリズムでは

ギャップ開始 (open) ペナルティ σ_o

ギャップ伸長 (extension) ペナルティ σ_x

を導入し、長さ n のギャップに対し σ_o + (n-1) σ_x のコストとします。ギャップ伸長ペナルティの値は開始ペナルティよりもずっと小さくします。ギャップ長に比例する値でコストが増加するペナルティ方式を（平行移動のアフィン変換にちなんで）アフィンギャップと呼びます。計算量を O(n²) に抑えたままアフィンギャップコストを計算するアルゴリズムを最初に示したのは日本の後藤修（ごとうおさむ）先生です。

この関数を使うと、アライメントスコアは

(マッチ数) − μ(ミスマッチ数) − σ_o(ギャップ開始数) − σ_x(ギャップスペース数)

と書けます。

まずスコア行列と同じサイズのギャップ管理用行列 P (横方向用)、Q (縦方向用) を用意します。それらの初期値は位置 ( i , 0 ), ( 0 , j ) に対してそれぞれ −σo − i σx と −σo − j σx とします。ここでは σo = 2, σx = 0.1 で計算してあります。そして以下の漸化式を用います。 $\begin{align} D_{i,j} &= max \begin{cases}0 \\ P_{i-1,j} \\ Q_{i,j-1} \\ D_{i-1,j-1} - \mu & \mbox{if } s_i \ne t_j \\ D_{i-1,j-1} +1 & \mbox{if } s_i = t_j \end{cases} \\ P_{i,j} &= max \begin{cases} D_{i-1, j} - \sigma_o \\ P_{i-1,j} - \sigma_x \end{cases} \\ Q_{i,j} &= max \begin{cases} D_{i, j-1} - \sigma_o \\ Q_{i,j-1} - \sigma_x \end{cases} \\ \end{align}$

初期行列の値（スコア用、ギャップ管理用とも共通） -2.0 -2.1 -2.2 -2.3 ... j -2.0 -2.1 -2.2  : i si, j 位置 (i , j) にいるとき、スコアの求め方を考えましょう。普通に考えるとセル (i , j) の左側全部と上側全部のセル（薄い影がかかっている領域）全ての位置から、ギャップが入る可能性を検討しなくてはなりません。それを素直に実行すると各セルの値を計算するのに O(n) 時間かかるため、全体で O(n3) 時間の計算量になります。その手間を省くために行列 P, Q を導入します。 行列 P, Q はそれぞれ位置 (i , j) が横方向および縦方向のギャップに含まれる際のスコアを記録しておく部分です。P と Q の漸化式をみると、それまで続いていたギャップを延長するか (P なら位置 i をスキップし、Q なら位置 j をスキップ）、新たに位置 i , j からギャップを開くか判断しています。これらの行列に前の方からギャップを伸長した場合のスコアを記録してあり、通常のスコアとの比較をおこないます。

参考文献