Needleman–Wunsch アルゴリズム

導入

グローバルアライメントを求めるアルゴリズムである！
詳しくはこちらをご覧ください。
https://bi.biopapyrus.jp/seq/alignment/needleman%E2%80%93wunsch.html

プログラム

日本語で書いてある事を、そのまま実装しただけです。
といっても書き方は人それぞれですので、これが正解ということはありません。
実装の一例というだけです。

(define nil '())

(define f (make-hash-table 'equal?))
(define p (make-hash-table 'equal?))

(define score-match 2)
(define score-mismatch -1)
(define score-gap -2)
(define d (* -1 score-gap))

(define (score w1 w2)
  (cond ((or (null? w1) (null? w2)) score-gap)
        ((eq? w1 w2) score-match)
        (else score-mismatch)))

(define (s x y)
  (score (hash-table-get f (make-key x -1) nil)
         (hash-table-get f (make-key -1 y) nil)))

(define (make-key a b)
  (string-append (number->string a) "," (number->string b)))

(define (get-x key)
  (if (null? key)
      nil
      (string->number (car (string-split key ",")))))

(define (get-y key)
  (if (null? key)
      nil
      (string->number (cadr (string-split key ",")))))

(define (init word1 word2 lim-x lim-y)
  (define (iter-x x w)
    (if (> x lim-x)
        'done
        (begin
          (hash-table-put! f (make-key x 0) (* x -1 d))
          (hash-table-put! f (make-key x -1) (car w))
          (hash-table-put! p (make-key x 0) (list (list 3 (make-key (- x 1) 0))))
          (iter-x (+ x 1) (cdr w)))))
  (define (iter-y y w)
    (if (> y lim-y)
        'done
        (begin
          (hash-table-put! f (make-key 0 y) (* y -1 d))
          (hash-table-put! f (make-key -1 y) (car w))
          (hash-table-put! p (make-key 0 y) (list (list 1 (make-key 0 (- y 1)))))
          (iter-y (+ y 1) (cdr w)))))
  (iter-x 1 word1)
  (iter-y 1 word2)
  (hash-table-put! f (make-key 0 0) 0))

(define (evaluate-scores lim-x lim-y)
  (define (iter-x x y)
    (if (> x lim-x)
        'done
        (let* ((key (make-key x y))
               (up-key (make-key x (- y 1)))
               (up-left-key (make-key (- x 1) (- y 1)))
               (left-key (make-key (- x 1) y))
               (up-left (+ (hash-table-get f up-left-key) (s x y)))
               (up (+ (hash-table-get f up-key) (* -1 d)))
               (left (+ (hash-table-get f left-key) (* -1 d)))
               (p0 (sort `((,up 1 ,up-key) (,up-left 2 ,up-left-key) (,left 3 ,left-key)) > car))
               (p1 (car p0))
               (p2 (cadr p0))
               (p3 (caddr p0)))
          (hash-table-put! f key (car p1))
          (cond ((= (car p1) (car p2) (car p3))
                 (hash-table-put! p key (sort (cons (cdr p1) (hash-table-get p key nil)) < car))
                 (hash-table-put! p key (sort (cons (cdr p2) (hash-table-get p key nil)) < car))
                 (hash-table-put! p key (sort (cons (cdr p3) (hash-table-get p key nil)) < car)))
                ((= (car p1) (car p2))
                 (hash-table-put! p key (sort (cons (cdr p1) (hash-table-get p key nil)) < car))
                 (hash-table-put! p key (sort (cons (cdr p2) (hash-table-get p key nil)) < car)))
                (else
                 (hash-table-put! p key (sort (cons (cdr p1) (hash-table-get p key nil)) < car))))
          (iter-x (+ x 1) y))))
  (define (iter-y y)
    (if (> y lim-y)
        'done
        (begin
          (iter-x 1 y)
          (iter-y (+ y 1)))
        ))
  (iter-y 1))

(define (make-global-alignment key)
  (if (equal? key (make-key 0 0))
      nil
      (let* ((val (hash-table-get p key nil))
             (up (assv 1 val))
             (up-left (assv 2 val))
             (left (assv 3 val)))
        (cons key
              (list
               (if up (make-global-alignment (cadr up)) nil)
               (if up-left (make-global-alignment (cadr up-left)) nil)
               (if left (make-global-alignment (cadr left)) nil))))))

(define (make-word paths)
  (define (iter p xx yy w1 w2)
    (if (null? p)
        (begin
          (print (reverse w1))
          (print (reverse w2))
          (print ""))
        (let ((x (get-x (car p)))
              (y (get-y (car p))))
          (cond ((and (> x xx) (> y yy))
                 (iter (cdr p) x y (cons (hash-table-get f (make-key x -1)) w1) (cons (hash-table-get f (make-key -1 y)) w2)))
                ((> x xx)
                 (iter (cdr p) x yy (cons (hash-table-get f (make-key x -1)) w1) (cons "-" w2)))
                ((> y yy)
                 (iter (cdr p) xx y (cons "-" w1) (cons (hash-table-get f (make-key -1 y)) w2)))))))
  (iter paths 0 0 nil nil))

(define (print-score len1 len2)
  (print "SCORE: " (hash-table-get f (make-key len1 len2))))

(define (print-words lst paths)
  (let* ((nodes (cdr lst))
         (up (car nodes))
         (up-left (cadr nodes))
         (left (caddr nodes)))
    (if (and (null? up) (null? up-left) (null? left))
        (make-word (cons (car lst) paths))
        (begin
          (and (not (null? up)) (print-words (car nodes) (cons (car lst) paths)))
          (and (not (null? up-left)) (print-words (cadr nodes) (cons (car lst) paths)))
          (and (not (null? left)) (print-words (caddr nodes) (cons (car lst) paths)))
          'done))))

(define (needleman-wunsch word-lst1 word-lst2)
  (hash-table-clear! f)
  (hash-table-clear! p)

  (let ((len1 (length word-lst1))
        (len2 (length word-lst2)))
    (init word-lst1 word-lst2 len1 len2)
    (evaluate-scores len1 len2)
    (print-score len1 len2)
    (print-words (make-global-alignment (make-key len1 len2)) nil)))

実行

gosh> (needleman-wunsch '(A T T G C) '(A T G C))
SCORE: 6
(A T T G C)
(A - T G C)

(A T T G C)
(A T - G C)

done

で、具体的に何の役に立つの？

それを言っちゃあ、おしまいよ。