読者です 読者をやめる 読者になる 読者になる

系列アライメントのアルゴリズムをRubyで実装した

アルゴリズム アルゴリズムデザイン Ruby
  • 動的計画法分かんない
  • しかも、なんか一次元じゃなくて二次元っぽいから怖い><

と思ってgkbrしていた次回アルゴリズムデザインのゼミで担当の系列アライメント(p246 6.6)のアルゴリズムですが、ミスドでホットカフェオレを3杯くらいおかわりしていたら結構簡単にできました。説明用の資料とかはまだ作ってないけど、たぶんどうにかなるでしょう。ということでとりあえずソース貼っておく。

# -*- coding: utf-8 -*-

# アルゴリズムデザイン 6.6 系列アライメント
# 動的計画法による最小コストの算出

class SequenceAlignment
  class Point
    # バックトラック用のクラス
    # ポインタがないので。。。
    attr_reader :x
    attr_reader :y
    def initialize(x = 0, y = 0)
      @x = x
      @y = y
    end
  end

  attr_reader :x
  attr_reader :y
  attr_reader :delta
  def initialize(x, y, delta)
    @x = x.split(//)
    @y = y.split(//)
    @delta = delta
    @a = [] #コストを入れておくための2次元配列
    (0..@x.length).each{|i|
      @a.push []
    }
  end

  def alignment
    # initialization
    (0..@x.length).each{|i|
      @a[i][0] = i * @delta
    }
    (0..@y.length).each{|j|
      @a[0][j] = j * @delta
    }

    # calculate minimum cost for sequence alignment
    # とりあえず表を埋めていけばいい
    (1..x.length).each{|i|
      (1..y.length).each{|j|
        @a[i][j] = min(alpha(i, j) + @a[i-1][j-1], delta + @a[i-1][j], delta + @a[i][j-1])
      }
    }
    return @a
  end

  def min(a, b, c)
    if (a < b) && (a < c)
      return a
    elsif (b < c)
      return b
    else
      return c
    end
  end

  def delta
    return @delta
  end

  def alpha(i, j)
    # iとjのindexは1から始まっているので-1しておく
    a = @x[i-1]
    b = @y[j-1]
    vowel = ["a", "b", "c", "d", "e"]
    if a == b
      return 0
    elsif ( (!vowel.include?(a) && !vowel.include?(b)) || # 母音と母音のマッチング
            (vowel.include?(a) && vowel.include?(b))) # 子音と子音のマッチング
      return 1
    else # 母音と子音のマッチングコスト
      return 3
    end
  end

  def argmin(i, j)
    a = @a[i-1][j]
    b = @a[i-1][j-1]
    c = @a[i][j-1]
    if (a < b) && (a < c)
      return Point.new(i-1, j)
    elsif (b < c)
      return Point.new(i-1, j-1)
    else
      return Point.new(i, j-1)
    end
  end

  def back_track(i = @x.length, j = @y.length, result = [])
    point = argmin(i, j)
    if i==0 && j==0 # 最後まで行きついた
      return result
    elsif result.length == 0 # 再帰の最初
      x = @x[point.x]
      y = @y[point.y]
      result.push [x, y]
      return back_track(point.x, point.y, result)
    else
      prev = result[result.length-1]
      if prev[0] == @x[point.x]
        x = nil 
      else 
        x = @x[point.x]
      end  
      
      if prev[1] == @y[point.y]
        y = nil 
      else
        y = @y[point.y]
      end
      result.push [x, y]
      return back_track(point.x, point.y, result)
    end
  end
  
  def minimum_cost
    return @a[@x.length][@y.length]
  end
end

a = "name"
b = "mean"

puts "Let us consider sequence alignment betweeen \"#{a}\" and \"#{b}\"."

sa = SequenceAlignment.new(a, b, 2)

puts "=" * 30
puts "Array for sequence alignment is as follows."
sa.alignment.reverse.each{|i|
  puts "#{i.join(" ")}"
}
puts "=" * 30

puts "Result is as follows."

sa.back_track.transpose.map{|a|
  a.reverse.map{|b|
    b.nil? ? "-" : b
  }.join(" ")
}.each{|i|
  puts i
}
puts "Minimised cost is #{sa.minimum_cost}."

実行結果。

/Users/syou6162/ruby% ruby sequence_alignment.rb
Let us consider sequence alignment betweeen "name" and "mean".
==============================
Array for sequence alignment is as follows.
8 6 4 5 6
6 4 4 5 4
4 3 2 3 5
2 1 3 5 6
0 2 4 6 8
==============================
Result is as follows.
n a - m e
m e a n -
Minimised cost is 6.

追記

上のソース、バックトラックの付近がおかしくなってました><。argminのindex、iとjがこんがらかってました。iは縦のほうのindex、jは横のindexになっているので、-1されるものが直感と違ったりしまします。

# -*- coding: utf-8 -*-

# アルゴリズムデザイン 6.6 系列アライメント
# 動的計画法による最小コストの算出

require "pp"

class SequenceAlignment
  class Point
    # バックトラック用のクラス
    # ポインタがないので。。。
    attr_reader :x
    attr_reader :y
    def initialize(x = 0, y = 0)
      @x = x
      @y = y
    end
    def to_s
      return "(#{x}, #{y})"
    end
  end

  attr_reader :x
  attr_reader :y
  attr_reader :a
  attr_reader :b
  attr_reader :delta
  def initialize(x, y, delta)
    @x = x.split(//)
    @y = y.split(//)
    @delta = delta
    @a = [] #コストを入れておくための2次元配列
    (0..@x.length).each{|i|
      @a.push []
    }
    @b = [] #バックトラックするための2次元配列
    (0..@x.length).each{|i|
      @b.push []
    }
  end

  def alignment
    # initialization
    (0..@x.length).each{|i|
      @a[i][0] = i * @delta
    }
    (0..@y.length).each{|j|
      @a[0][j] = j * @delta
    }

    (1..@x.length).each{|i|
      @b[i][0] = Point.new(i-1, 0)
    }
    (1..@y.length).each{|j|
      @b[0][j] = Point.new(0, j-1)
    }

    # calculate minimum cost for sequence alignment
    (1..@x.length).each{|i|
      (1..@y.length).each{|j|
        @a[i][j] = [alpha(i, j) + @a[i-1][j-1], delta + @a[i-1][j], delta + @a[i][j-1]].min
        @b[i][j] = argmin(i, j)
      }
    }
    return @a
  end

  def argmin(i, j)
    # transposeになっていることに注意!!
    a = delta + @a[i][j-1] # 左
    b = alpha(i, j) + @a[i-1][j-1]  # 左上
    c = delta + @a[i-1][j] # 上
    if (a < b) && (a < c)
      return Point.new(i, j-1)
    elsif (b < c) 
      return Point.new(i-1, j-1)
    else
      return Point.new(i-1, j)
    end
  end

  def delta
    return @delta
  end

  def alpha(i, j)
    # iとjのindexは1から始まっているので-1しておく
    a = @x[i-1]
    b = @y[j-1]
    vowel = ["a", "b", "c", "d", "e"]
    if a == b
      return 0
    elsif ( (!vowel.include?(a) && !vowel.include?(b)) || # 母音と母音のマッチング
            (vowel.include?(a) && vowel.include?(b))) # 子音と子音のマッチング
      return 1
    else # 母音と子音のマッチングコスト
      return 3
    end
  end
  
  def back_track
    def back_track_prepare(i = @x.length, j = @y.length, result = [])
      if i==0 && j==0 # 最後まで行きついた
        return result
      end
      point = @b[i][j]
      result.push point
      return back_track_prepare(point.x, point.y, result)
    end
    
    prev_x = 1.0/0
    prev_y = 1.0/0
    str_x = []
    str_y = []
    back_track_prepare.reverse.each{|i|
      if i.y == prev_y # 上に移動するとき
        str_x.push @x[i.x] 
        tmp = str_y.pop # 末尾にあったやつを取り出して
        str_y.push nil # nilを入れてから
        str_y.push tmp # 取り出したやつを入れる
      elsif i.x == prev_x # 右に移動するとき
        str_y.push @y[i.y]
        tmp = str_x.pop
        str_x.push nil
        str_x.push tmp
      else # 右上に移動するとき
        str_x.push @x[i.x]
        str_y.push @y[i.y]
      end
      prev_x = i.x
      prev_y = i.y
    }
    return [str_x, str_y]
  end

  def minimum_cost
    return @a[@x.length][@y.length]
  end
end

a = "mean"
b = "name"

puts "Let us consider sequence alignment betweeen \"#{a}\" and \"#{b}\"."

sa = SequenceAlignment.new(a, b, 2)

puts "=" * 30
puts "Array for sequence alignment is as follows."
sa.alignment.reverse.each{|i|
  puts "#{i.join(" ")}"
}
puts "=" * 30

puts "Result is as follows."

puts "Minimised cost is #{sa.minimum_cost}."

sa.back_track.map{|a|
  a.map{|b|
    b.nil? ? "-" : b
  }
}.map{|i|
  i.join("  ")
}.each{|i|
  puts i
}

実行結果。

/Users/syou6162/ruby% ruby sequence_alignment.rb
Let us consider sequence alignment betweeen "mean" and "name".
==============================
Array for sequence alignment is as follows.
8 6 5 4 6
6 5 3 5 5
4 3 2 4 4
2 1 3 4 6
0 2 4 6 8
==============================
Result is as follows.
Minimised cost is 6.
m  e  a  n  -
n  -  a  m  e