エクセルVBAで動的計画法でDNA配列のアラインメントを作成 Needleman-Wunschアルゴリズムの実装

はじめに

この記事は、DNA配列のアラインメントを返すエクセルの関数についての記事です。

Needleman-Wunschによるアルゴリズムによるアラインメントですので関数名はnwa()としました。

方針

最終的に作る関数は、

とすると、以下のようなアラインメントを返す関数です。

エクセルシート上で、アルゴリズムを示した後で、関数の作成に取り掛かります。

エクセルシート上でアルゴリズムをチェック

2つの行列を用意します。アラインしたい2つのDNA配列の長さをmとnとすれば、縦m x 横nの大きさの行列です。行列は、スコア行列と、トレースバック行列です。初期状態は以下の通りです。

2つのDNA配列は、一文字ずつ区切ってセルに入れ、行と列のラベルとして使います。

スコア行列

初期状態を以下のようにします。

トレースバック行列

初期状態を以下のようにします。

行列の埋め方

一定の規則に従って、行列を左上から埋めていきます。規則は以下の通りです。
これから埋めようとするセルの番地をx,yとした時に、以下の3つのスコアのうち、最も高いスコアをセル(x,y)の値とする(2つ以上が同じ値で、最高値になることもある)。

1. そのセルに対応する塩基が同じ場合(セルの列ラベルと行ラベルを比べて同じ場合)、セル(x-1,y-1)の値に1を加えた値、そうでなければ、セル(x-1,y-1)の値から1を減じた値
2. セル(x-1,y)の値から1を減じた値
3. セル(x,y-1)の値から1を減じた値

要するにセル(x-1,y-1)は左上のセル、セル(x-1,y)は左のセル、セル(x,y-1)は上のセルです。
またこの時、トレースバック行列の対応するセルにも値を入れます。規則は、

1. 最高スコアを示したのが「1」だったら、LEFTUP
2. 最高スコアを示したのが「2」だったら、LEFT
3. 最高スコアを示したのが「3」だったら、UP、
   とします。

エクセル上では、こんな風になります。MAXは、引数の中から最大値を返す関数です。後でコピーペーストする用に、絶対参照を使っています。

こっちはちょっとややこしいですね。

あとは作った式をコピーペーストして行列を埋めます。

トレースバック行列の右下のセルから、記しておいた通りに、さかのぼっていきます。LEFTUPと書いてあったら左上に飛ぶ、というようにします。辿ったセルの背景を赤くしておきました。

これで、2つの配列をどう並べたら良いかがわかったことになりますが、この行列をもとにVBAを使わないでアラインメントを作るのは大変なので、以後はVBAで関数を作ります。

VBAコード

アルゴリズムさえわかってしまえば、特に難しいところはないですね。ここでは二次元の行列の代わりに、一次元の行列を使っています。幅nの行列の座標(x、y)は、一次元行列では、indexが$x + n*y$になります。

'DNA配列を2つ受け取ります。DNAXとDNAYとします。

Public Function nwa(DNAX As String, DNAY As String) As String

Dim n As Integer
Dim m As Integer
Dim u As Integer
Dim l As Integer
Dim ul As Integer
Dim score As Integer
Dim alignmentX As String
Dim alignmentY As String
Dim alignmentLine As String

alignmentX = ""
alignmentY = ""
alignmentLine = ""
n = Len(DNAX)
m = Len(DNAY)

Dim scm() As Integer 'スコアマトリックスscm
Dim trm() As String 'トレースバックマトリックス trm
ReDim scm(0 To n * m)
ReDim trm(0 To n * m)


'マトリックスの初期化
For K = 1 To n - 1
    scm(K) = -2 * K
    trm(K) = "LEFT"
Next K

For J = 1 To m - 1
    scm(J * n) = -2 * J
    trm(J * n) = "UP"
Next J

'規則従って2つのマトリックスを埋める
For Y = 1 To m - 1
    For X = 1 To n - 1
        u = scm(X + (Y - 1) * n) - 1
        l = scm(X - 1 + Y * n) - 1
        If Mid(DNAX, X, 1) = Mid(DNAY, Y, 1) Then
            ul = scm(X - 1 + (Y - 1) * n) + 1
            Else
            ul = scm(X - 1 + (Y - 1) * n) - 1
        End If

        If ul >= u And ul >= l Then
                scm(X + Y * n) = ul
                trm(X + Y * n) = "LEFTUP"
        ElseIf u >= ul And u >= l Then
                scm(X + Y * n) = u
                trm(X + Y * n) = "UP"
        ElseIf l >= ul And l >= u Then
                scm(X + Y * n) = l
                trm(X + Y * n) = "LEFT"
        End If
    Next X
Next Y


'スコア算出と、アラインメントの作成
'アラインメントの各行を別々に作って最後に足し合わせる。
'ギャップ部分は - を挿入
score = 0
X = n - 1
Y = m - 1
While X <> 0 Or Y <> 0
        q = trm(X + Y * n)
        If q = "LEFT" Then
            alignmentX = Mid(DNAX, X, 1) & alignmentX
            alignmentY = "-" & alignmentY
            alignmentLine = " " & alignmentLine
            X = X - 1
        ElseIf q = "UP" Then
           alignmentX = "-" & alignmentX
            alignmentY = Mid(DNAY, Y, 1) & alignmentY
            alignmentLine = " " & alignmentLine
            Y = Y - 1
        ElseIf q = "LEFTUP" Then
            If Mid(DNAX, X, 1) = Mid(DNAY, Y, 1) Then
                score = score + 1
                alignmentLine = "|" & alignmentLine
                Else
                alignmentLine = " " & alignmentLine
            End If
            alignmentX = Mid(DNAX, X, 1) & alignmentX
            alignmentY = Mid(DNAY, Y, 1) & alignmentY
            X = X - 1
            Y = Y - 1
        End If
Wend

'VBA関数の作り方のルールに従って、最後に関数名=返り値とする。
nwa = alignmentX & vbNewLine & alignmentLine & vbNewLine & alignmentY
End Function

終わりに

最後にscoreを返すように変更すれば、そういう使い方もできますね。

以前作った関数ですが、ここで公開することにしました。かつて読んだ本にあった解説がなかなか理解できず、煩悶したことがあります。ここでの記事が、悩める人の理解の助けになれば幸いです。