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

きくらげ観察日記

好きなことを、適当に。

Needleman-Wunsch法により配列の大域アラインメントを求める

アラインメントとは、2つの配列s, tに欠損記号「-」を挿入して、出来る限り2つの一致度が高くなるように並べたものです。
例えば、2つの文字列"ABCD"と"ACDE"があったとき、これは1つめの文字列から'B'を削除して末尾に'E'を挿入たものだと考えるのが自然でしょう。この時の2つの文字列のアラインメントは、それぞれ"ABCD-"と"A-CDE"となります。

アルゴリズム自体の説明は「Needleman-Wunsch」とかでググれば出てくるので割愛します。
Gaucheによる実装がこちら。
alignment関数は2つの文字列s, tとパラメータweight, gap-penaltyを受け取り、sとtの類似度の最大値と、類似度を最大値にするようなsとtのアラインメントをそれぞれ返します。なお、アラインメントはリストで表され、各要素は文字か欠損値#fで表されます。

(use util.match)
(use gauche.array)
(use gauche.parameter)

(define weight
  (make-parameter
   (^(x y) (if (eqv? x y) 1 -1))))

(define gap-penalty
  (make-parameter 1))

(define arr~ array-ref)

(define (alignment s t)
  (define w (weight))
  (define d (gap-penalty))
  (define s-len (string-length s))
  (define t-len (string-length t))
  (define s-arr (apply array (shape 0 s-len) (string->list s)))
  (define t-arr (apply array (shape 0 t-len) (string->list t)))
  (letrec ((calc
            (^(i j)
              (lazy
               (match (cons i j)
                 ((_ . 0) (* i (- d)))
                 ((0 . _) (* j (- d)))
                 ((_ . _) (max (- (force (arr~ dp (- i 1) j)) d)
                               (- (force (arr~ dp i (- j 1))) d)
                               (+ (force (arr~ dp (- i 1) (- j 1)))
                                  (w (arr~ s-arr (- i 1))
                                     (arr~ t-arr (- j 1))))))
                 ))))
           (dp
            (tabulate-array
             (shape 0 (+ s-len 1) 0 (+ t-len 1)) calc)))
    (receive (s-align t-align) (split (path s-arr t-arr dp))
      (values
       (force (arr~ dp s-len t-len))
       s-align t-align)
      )))

(define (path s-arr t-arr dp)
  (define w (weight))
  (define d (gap-penalty))
  (let loop ((i (array-length s-arr 0))
             (j (array-length t-arr 0))
             (result ()))
    (if (and (= i 0) (= j 0))
        result
        (let1 dp-val (force (arr~ dp i j))
          (cond
           ((= dp-val (- (force (arr~ dp (- i 1) j)) d))
            (loop (- i 1) j
                  (cons (cons (arr~ s-arr (- i 1)) #f)
                        result)))
           ((= dp-val (- (force (arr~ dp i (- j 1))) d))
            (loop i (- j 1)
                  (cons (cons #f (arr~ t-arr (- j 1)))
                        result)))
           (else
            (loop (- i 1) (- j 1)
                  (cons (cons (arr~ s-arr (- i 1)) (arr~ t-arr (- j 1)))
                        result)))
           )))))

(define (split lst)
  (match lst
    (() (values () ()))
    (((x . y) . rest) (receive (xs ys) (split rest)
                        (values (cons x xs) (cons y ys))))))

実行結果:

gosh> (alignment "ABCD" "ACDE")
1
(#\A #\B #\C #\D #f)
(#\A #f #\C #\D #\E)

動的計画法は遅延評価のある言語だと漸化式を書くだけで綺麗に解ける、ということを示したかったのですが、Schemeだとそこまで綺麗になりませんでした。HaskellだとSTモナドを使ってdpテーブルを破壊的に更新するよりも100倍綺麗に書けます。

なお、普通の書き方をするとalignment関数は以下のようになります。

(define (alignment-destructive s t)
  (define w (weight))
  (define d (gap-penalty))
  (define s-len (string-length s))
  (define t-len (string-length t))
  (define s-arr (apply array (shape 0 s-len) (string->list s)))
  (define t-arr (apply array (shape 0 t-len) (string->list t)))
  (define dp
    (make-array
     (shape 0 (+ s-len 1) 0 (+ t-len 1)) #f))
  (define (calc! i j)
    (define val (arr~ dp i j))
    (or val
        (let1 x (match (cons i j)
                  ((_ . 0) (* i (- d)))
                  ((0 . _) (* j (- d)))
                  ((_ . _) (max (- (calc! (- i 1) j) d)
                                (- (calc! i (- j 1)) d)
                                (+ (calc! (- i 1) (- j 1))
                                   (w (arr~ s-arr (- i 1))
                                      (arr~ t-arr (- j 1))))))
                  )
          (array-set! dp i j x)
          x
          )))
  (calc! s-len t-len)
  (receive (s-align t-align) (split (path s-arr t-arr dp))
    (values
     (arr~ dp s-len t-len)
     s-align t-align)
    ))