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

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

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

#方針
最終的に作る関数は、
<img width="931" alt="image.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/595771/1b1ad38f-4be9-8812-55da-79c645035753.png">
とすると、以下のようなアラインメントを返す関数です。
<img width="927" alt="image.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/595771/974e72c3-2527-0ae1-4826-f104d5a592e2.png">


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

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

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

##スコア行列
初期状態を以下のようにします。
<img width="958" alt="image.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/595771/a17b3a52-aa64-21b1-e283-cef6194d49b9.png">

##トレースバック行列
初期状態を以下のようにします。
<img width="959" alt="image.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/595771/83d3dc21-c0da-1aae-031e-ff70d7dcc37b.png">

#行列の埋め方
一定の規則に従って、行列を左上から埋めていきます。規則は以下の通りです。
これから埋めようとするセルの番地を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は、引数の中から最大値を返す関数です。後でコピーペーストする用に、絶対参照を使っています。
<img width="436" alt="image.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/595771/898a5a49-7e22-0f5d-9662-1136872bb657.png">
こっちはちょっとややこしいですね。
<img width="532" alt="image.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/595771/a5011181-e717-a27f-79f2-4e37cd8a87bf.png">

あとは作った式をコピーペーストして行列を埋めます。
<img width="752" alt="image.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/595771/e7a6ce2d-a592-f4b0-f8c4-ab6666c294df.png">

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

<img width="753" alt="image.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/595771/380c2a56-fc57-ef5b-1ce9-e07ba171cca9.png">

これで、2つの配列をどう並べたら良いかがわかったことになりますが、この行列をもとにVBAを使わないでアラインメントを作るのは大変なので、以後はVBAで関数を作ります。
#VBAコード
アルゴリズムさえわかってしまえば、特に難しいところはないですね。ここでは二次元の行列の代わりに、一次元の行列を使っています。幅nの行列の座標(x、y)は、一次元行列では、indexが$x + n*y$になります。

```vb
'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を返すように変更すれば、そういう使い方もできますね。

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