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

きくらげ観察日記

好きなことを、適当に。

Smith-Waterman法で配列の局所アラインメントを求める

今度は局所アラインメントを求めます。
前回の大域アラインメントと違って、今度は最も一致度の高い部分文字列を求めるアルゴリズムになります。

type alignment = char option list

let maximums_by (key : 'a -> 'b) (lst : 'a list) : 'a list =
  let rec loop max_val maxs xs =
    match xs with
    | [] -> maxs
    | x :: xs' ->
       let x_val = key x
       in if max_val < x_val then loop x_val [x] xs'
          else if max_val = x_val then loop max_val (x :: maxs) xs'
          else loop max_val maxs xs'
  in match lst with
     | [] -> raise (Failure "maximums_by: empty list")
     | max :: xs -> loop (key max) [max] xs

let maximum (lst : 'a list) : 'a =
  List.hd (maximums_by (fun x -> x) lst)

let fromto (from : int) (to_ : int) : int list =
  let rec loop i acc =
    if i < from then acc
    else loop (i - 1) (i :: acc)
  in loop to_ []

let prod (xs : 'a list) (ys : 'b list) : ('a * 'b) list =
  List.flatten (List.map (fun x -> List.map (fun y -> (x, y)) ys) xs)

module type ALIGNMENT =
  sig
    val align : string -> string -> (alignment * alignment) list
  end

module type WEIGHT =
  sig
    val weight : char -> char -> int
    val gap_penalty : int
  end

module Alignment (W : WEIGHT) : ALIGNMENT =
  struct
    let path s t m n calc : alignment * alignment =
      let rec loop i j acc =
        match i, j with
        | 0, 0 -> List.split acc
        | _, _ ->
           let dp_val = calc i j
           in if dp_val = 0 then
                List.split acc
              else if dp_val = calc (i - 1) j - W.gap_penalty then
                loop (i - 1) j ((Some s.[i - 1], None) :: acc)
              else if dp_val = calc i (j - 1) - W.gap_penalty then
                loop i (j - 1) ((None, Some t.[j - 1]) :: acc)
              else loop (i - 1) (j - 1) ((Some s.[i - 1], Some t.[j - 1]) :: acc)
         in loop m n []

    let align s t =
      let s_len = String.length s in
      let t_len = String.length t in
      let dp = Array.make_matrix (s_len + 1) (t_len + 1) None in
      let rec calc i j =
        let value = dp.(i).(j) in
        match value with
        | Some (v) -> v
        | None ->
           let new_value =
             match i, j with
             | 0, _ -> i * -W.gap_penalty
             | _, 0 -> j * -W.gap_penalty
             | _, _ -> maximum [
                           0
                         ; calc (i - 1) j - W.gap_penalty
                         ; calc i (j - 1) - W.gap_penalty
                         ; calc (i - 1) (j - 1) + W.weight s.[i - 1] t.[j - 1]
                         ]
           in dp.(i).(j) <- Some(new_value);
              new_value in
      let indices = maximums_by (fun (i, j) -> calc i j)
                                (prod (fromto 0 s_len) (fromto 0 t_len)) in
      List.map (fun (m, n) -> path s t m n calc) indices
  end

module Default_weight : WEIGHT =
  struct
    let weight c d = if c = d then 5 else -3
    let gap_penalty = 4
  end

module Default_alignment = Alignment(Default_weight)

Default_weightのweightとgap_penaltyの値は、以下のサイトを参考にしました。
vlab.amrita.edu


実行結果:

# Default_alignment.align "GACTTAC" "CGTGAATTCAT";;
- : (alignment * alignment) list =
[([Some 'G'; Some 'A'; Some 'C'; Some 'T'; Some 'T'; Some 'A'; Some 'C'],
  [Some 'G'; Some 'A'; Some 'A'; Some 'T'; Some 'T'; None; Some 'C']);
 ([Some 'G'; Some 'A'; Some 'C'; Some 'T'; Some 'T'; None; Some 'A'],
  [Some 'G'; Some 'A'; Some 'A'; Some 'T'; Some 'T'; Some 'C'; Some 'A'])]

得られた局所アラインメントは、

GACTTAC
GAATT-C

GACTT-A
GAATTCA

の2つです。先ほどのサイトの結果と一致しているので、アルゴリズムは正しく動いているようです。

余談ですが、このアルゴリズムにおけるweight, gap_penaltyのような、引数にするほどでもないけど自由に設定したい値を渡すのにFunctorを使うと便利です。他のオブジェクト指向言語ならインスタンスのメンバに割り当てたり、Schemeならparameterを使ったりするのでしょうが、Haskellだとなかなかこういうことができなくて少し不便です。