Hirschberg 法 - シーケンス・アラインメント用

最新版の GNU Diffutils のテキスト差分のロジック部分は Gnulib の lib/diffseq.h に切り出されています。 lib/diffseq.h を読むと、 ロジックの基本は Myers An O(ND) Difference Algorithm and Its Variations (1986) の 4 節に記載されている Hirschberg法 (1975) と組み合わせて線形作業領域サイズで済むように手直ししたもので、 それに加えてワーストケースに陥らないようにヒューリスティックな手法を組み込んだものになっています。 なお、 日本語のウェブサイトで Myers 法として紹介されることが多い手法は 4 節の改良版ではなく、 その前の 3 節で説明してある基本となっている方です。

元になっている Hirschberg 法は、 分子生物学のために発見された手法で、 核酸たんぱく質の 2 本のシーケンスを比較するのに使われます。このシーケンス・アラインメントは、 比較したいシーケンスを上下に並べ挿入・削除されている区間をハイフンで埋めたギャップで表します。

例:
AGTACGCA
--TATGC-

Hirschberg 法はさらに元になっている長さ N と長さ M の記号列のシーケンス・アラインメントを O(NM) の作業領域を使う Needleman-Wunsch 法 (1970)を分割統治法で改良したものです。 Needleman-Wunsch 法は、 編集グラフ格子と同等の N 行 M 列のスコア格子を使い、 O(NM) 時間で計算をおこないます。 この手の編集グラフ格子を扱う手法の系統は原則として Hirschberg 法で作業領域のサイズを O(M) または O(N) に減らして処理する手法に変換することができます。

Needleman-Whunsch 法は編集距離を求める Wagner-Fischer 法 (1974) よりも前に発見されたもので、 記号同士にスコアを行列で指定することで複数の候補の中からギャップの入れ方を選択できるようになっています。 編集グラフを格子で求める点は Wagner-Fischer 法と同じです*1。 違いは、 Needleman-Whunsch 法ではスコアが大きい方が編集距離が短いことで、 編集グラフの格子を求めるときに min ではなく、max を使います。

def needleman_whunsch(a, b)
  s = lambda{|x, y| x == y ? +2 : -1 }  # 変更なしスコア +2 と、 修正スコア -1
  d = -2                                # 挿入スコア -2 と削除スコア -2
  # 格子を求めます
  grid = Array.new(a.size + 1) { Array.new(b.size + 1, 0) }
  (0 .. a.size).each {|i| grid[i][0] = d * i }
  (0 .. b.size).each {|j| grid[0][j] = d * j }
  (1 .. a.size).each do |i|
    (1 .. b.size).each do |j|
      match_score = grid[i - 1][j - 1] + s[a[i - 1], b[j - 1]]
      delete_score = grid[i][j - 1] + d
      insert_score = grid[i - 1][j] + d
      grid[i][j] = [match_score, delete_score, insert_score].max
    end
  end
  # 格子からシーケンス・アラインメントを求めます
  a_alignment = ""
  b_alignment = ""
  i = a.size
  j = b.size
  while i > 0 or j > 0
    if i > 0 and j > 0 and grid[i][j] == grid[i - 1][j - 1] + s[a[i - 1], b[j - 1]]
      a_alignment = a[i - 1] + a_alignment
      b_alignment = b[j - 1] + b_alignment
      i = i - 1
      j = j - 1
    elsif i > 0 and grid[i][j] == grid[i - 1][j] + d
      a_alignment = a[i - 1] + a_alignment
      b_alignment = "-" + b_alignment
      i = i - 1
    else #if j > 0 and grid[i][j] == grid[i][j - 1] + d
      a_alignment = "-" + a_alignment
      b_alignment = b[j - 1] + b_alignment
      j = j - 1
    end
  end
  [a_alignment, b_alignment]
end

needleman_whunsch("AGTACGCA", "TATGC").each {|s| puts s }

これに Hirschberg 法を摘要して作業領域を削ったコードに書き直します。 作業領域が O(MN) 必要だったスコアの格子からは、 中央の一行分だけを求めて使います。 2 つの記号列の一方を 2 分割して先頭側から順方向に求めた格子の行と、 末尾側から逆方向に求めた格子の行を比較して、 スコアの合計が最大の桁を探します。 この桁でもう一方の記号列を 2 分割し、 分割されたそれぞれの記号列のうち先頭側からの一組と末尾側の一組をそれぞれ再帰でシーケンス・アラインメントを求めて結合します。 下記のコードでは、 分割位置を求めるための末尾側から行を求めるところを、 横着して文字列の逆転メソッド reverse を使っていますが、 本来ならスコアの格子を逆順に求める対になるルーチンを書くべきでしょう。

def hirschberg(a, b)
  if a.size == 0
    ['-' * b.size, b.dup]
  elsif b.size == 0
    [a.dup, '-' * a.size]
  elsif a.size == 1
    i = b.rindex(a) || b.size - 1
    ["-" * i + a + "-" * (b.size - i - 1), b.dup]
  elsif b.size == 1
    j = a.rindex(b) || a.size - 1
    [a.dup, "-" * j + b + "-" * (a.size - j - 1)]
  else
    amid = a.size / 2
    a0, a1 = a[0 ... amid], a[amid ... a.size] || ""
    forw_score = needleman_whunsch_score(a0, b)
    back_score = needleman_whunsch_score(a1.reverse, b.reverse)
    bmid = (0 .. b.size).max_by {|j| forw_score[j] + back_score[b.size - j] }
    b0, b1 = b[0 ... bmid], b[bmid ... b.size] || ""
    a0_alignment, b0_alignment = hirschberg(a0, b0)
    a1_alignment, b1_alignment = hirschberg(a1, b1)
    [a0_alignment + a1_alignment, b0_alignment + b1_alignment]
  end
end

def needleman_whunsch_score(a, b)
  s = lambda{|x, y| x == y ? +2 : -1 }
  d = -2
  prev = Array.new(b.size + 1, 0)
  curr = Array.new(b.size + 1) {|j| d * j }
  (1 .. a.size).each do |i|
    prev, curr = curr, prev
    curr.fill(0)
    (1 .. b.size).each do |j|
      match_score = prev[j - 1] + s[a[i - 1], b[j - 1]]
      delete_score = curr[j - 1] + d
      insert_score = prev[j] + d
      curr[j] = [match_score, delete_score, insert_score].max
    end
  end
  curr
end

hirschberg("AGTACGCA", "TATGC").each {|s| puts s }

*1:厳密には、1970 年の Needleman-Wunsch 法では 3 次元格子を使って計算します。 後に 2 次元格子を使うように改良され、 現代の慣習では 2 次元格子を使う方を Needleman-Wunsch 法と呼ぶようです。