第7章の練習，投票

再帰版のアラインメント

def max(x,y)
  if y < x	
    x
  else
    y
  end		
end

def g()
  -2
end

def q(x,y)
  if x==y
    2
  else
    -1
  end
end

def align_sub(s,t,i,j)
  if i==0 || j==0
    i*g() + j*g()
  else
    v0 = align_sub(s,t,i,  j-1) + g()
    v1 = align_sub(s,t,i-1,j-1) + q(s[i-1],t[j-1])
    v2 = align_sub(s,t,i-1,j)   + g()
    max(v0,max(v1,v2))
  end
end

def align_rec(s,t)
  align_sub(s,t,s.length(),t.length())
end

align_rec("ATAG","AAC")
align_rec("GACG","GCAG")

再帰版の実行時間の計測

require "benchmark"
Benchmark.measure{ align_rec("012345","543210") }.total()
Benchmark.measure{ align_rec("0123456","6543210") }.total()
Benchmark.measure{ align_rec("01234567","76543210") }.total()
Benchmark.measure{ align_rec("012345678","876543210") }.total()
Benchmark.measure{ align_rec("0123456789","9876543210") }.total()

ここは?

1. 1/16
2. 3/16
3. 5/16
4. 7/16
5. 9/16

動的計画法を用いたアラインメント

def max(x,y)
  if y < x	
    x
  else
    y
  end		
end

def make1d(w)
  a=Array.new(w)
  for i in 0..(w-1)
    a[i]=0
  end
  a
end
def make2d(h,w)
  a=Array.new(h)
  for i in 0..(h-1)
    a[i]=make1d(w)
  end
  a
end

def g()
  -2
end

def q(c,d)
   if c == d
       return 2
   else
       return -1
   end
end

def align(s,t)
  m = s.length()
  n = t.length()
  a = make2d(m+1,n+1)
  a[0][0] = 0
  for j in 1..n
    a[0][j] = a[0][j-1] + g()
  end
  for i in 1..m
    a[i][0] = a[i-1][0] + g()
  end
  for i in 1..m 
    for j in 1..n
      v2 = a[i][j-1] + (ここを埋める) # 上から
      v1 = a[i-1][j-1] + (ここを埋める) # 左上から
      v0 = a[i-1][j] + g() # 左からのコスト
      a[i][j]=max(v0,max(v1,v2)) # これらの最大のものを返す
    end
  end
  a
end

def traceback(a,s,t)
  u = ""
  v = ""
  i = s.length()
  j = t.length()
  while i>0 || j>0
    if j>0 && a[i][j] == a[i][j-1] + g()
      # 上から求められた場合
      u = "-"           + u
      v = t[j-1 .. j-1] + v
      j = j - 1 # go left
    else
      if i>0 && j>0 && a[i][j] == a[i-1][j-1] + q(s[i-1], t[j-1])
          # 左上から求められた場合
        u = (ここを埋める) + u
        v = (ここを埋める) + v
        j = j - 1
        i = i - 1
      else
        if i>0 && a[i][j] == a[i-1][j] + g()
          # 左から求められた場合
          u = (ここを埋める) + u
          v = (ここを埋める) + v
          i = i - 1
        end
      end
    end
  end
  [u,v]
end

def align_dp(s,t)
  traceback(align(s,t),s,t)
end

ここは?

1. 1
2. 2
3. 3
4. 4
5. 5
6. -1
7. -2
8. -3
9. -4

練習

• align.rb を完成させて、ATAG と AAC のアラインメントを求めよう。
• RNase_P.rb をロードすると、文字列を返す（引数なしの）関数 seq0() と seq1() が定義される。これらのアラインメントを求めよう。

align_dp(seq0(),seq1())

進捗状況の確認

1. align の結果に対して traceback が正しく動いた時点で投票してください。
2. align は動いたが、traceback の結果が正しくない。
3. align はできたが、うまく動かない。
4. align がまだできていない。