全局比對與動態規划

前言

學過生信的肯定知道Needleman–Wunsch算法和Smith–Waterman 算法，一個用來進行全局比對，一個用來進行局部比對。單純看算法抽象后的算法公式的話也不復雜，用兩個短序列來套公式計算的話，畫過箭頭的同學都知道不難。本文簡單講述Needleman–Wunsch算法是如何利用動態規划來尋找兩序列最大相似度。
注：有些數學公式沒渲染出來可點擊鏈接查看：https://www.jianshu.com/p/002bbebcaaef
注：代碼實現見 https://www.cnblogs.com/huanping/p/11273391.html

動態規划

動態規划是數學上的優化方法，也是計算機編程的一種方法。對於一名生物狗，從數學層面去看動態規划，就太過硬核了。我們就從計算機編程角度理解動態規划就好了。
動態規划是一種解決問題的方法，將一個復雜大問題拆分為幾個小的問題，從幾個小問題的最優解推出大問題的最優解。當然這不是說將一個大問題拆分了就能用動態規划了。同時問題也需要滿足兩個性質，這兩個性質后面再講。先回到我的主要探討問題全局比對上：對於兩條蛋白質或者核酸序列來說，尋找它們的最大相似度，如何拆分成幾個小問題呢？（提示：矩陣、算法公式）

拆分問題

就這樣憑空想，如何拆分，我覺得我是想破腦袋也想不了T_T. 所以我們按照，學全局比對的套路手動來一遍序列比對操作，從中看看如何拆分問題。

為方便計算，得分矩陣就弄個簡單點的：

N	A	T	G	C	-
A	10	5	5	5	-5
T	5	10	5	5	-5
G	5	5	10	5	-5
C	5	5	5	10	-5
-	-5	-5	-5	-5	0

也就說完全匹配得10分，不全得5分，匹配到空位-5分，空位對空位一開始用來初始化時，才有意義。在之后比對中空位對空位沒有意義。
先給出公式（d為gap罰分，即鹼基比對空格的得分）：

現假設我們有兩條序列Seq1: CATCCCCCAAA 和 Seq2: CAAAA。長度分別為11，5。根據它們的長度畫一個合適的二維矩陣，填充好初始數據：

加粗的分別是 Seq1,Seq2，其中 -代表空格。一開始空格對空格，我們初始化，記$F(0,0)=0$，同時記空格分別對應的行和列是第0行，第0列。然后第1行，第1列分別是兩條首字母對應的行列。用$(i，j)$代表（行，列），$i,j$能有多大，取決於序列對的長度。$s(x_i,y_i)$的分數是一條序列中第$i$個鹼基和另一條序列中第$j$個鹼基比對的分數。放在矩陣中就是了第$i$行的鹼基和第$j$列的鹼基比對結果。

按照$F(i,j)$函數公式，$F(i,j)$的值由左上的值($F(i-1,j-1)$),上方值($F(i-1,j)$)和左方值($F(i,j-1)$)確定。
（注:這里的左上，上，左是一個相對位置，在當前設定矩陣坐標系中才有意義，這里是將$（0,0）$位置放在了左上方，向右和向下延伸。但其實，你也可以將出放在其它位置，比如右下方）

極端比對

我們看下第0行，它們只能由左方的值得到。我們算幾個值：
$F(0,1)=F(0,1-1)-5=0-5$
$F(0,2)=F(0,2-1)-5=-5-5=-10$
$F(0,3)=F(0,3-1)-5=-10-5=-15$
第0列，只能由上方的值得到，計算同理，不再演示了。
因為第0行和0列，很容易就得出來了，所以一般畫矩陣時，就順便一起畫出來了。

你應該知道，往右或往下延申，則比對到了空格。那這是為啥？
就拿第0行的比對情況舉個例子，第0行中，$i=0$已經固定了，它代表一條序列（后面稱為$S1$序列）的第0元素。$j=0，1,2,3...$,分別代表另一條序列（后面稱為$S2$序列）第0，1，2，3...元素。序列第一個元素是序列的第一個鹼基。第0個是啥玩意，第0個啥也不是，那我們就把它當做空格吧。

1. 序列比對從兩條序列的第$（0,0）$元素開始，兩個空格比對，得分為0。
2. $i=0$固定，一直表示$S1$序列第0個元素，但是$j$是一直增加的，所以$j=1$時，就是用$S2$序列的的第$j$個元素和$S1$序列的第$i=0$個元素比對，但是$S1$序列第$i=0$元素已經和$S2$序列的第$j=0$元素比對過了，你不能重復比對啊，你也不能私自和$S1$序列的$i=1$元素比對，所以沒辦法，就只能和一個空格比對，產生空位罰分，在之前得分基礎上上加上空位罰分。
3. $S2$序列的第$2,3,4...$元素都是如此，只能比對空格，產生空位罰分。

比對往右延申的過程就是行數$i$固定了，沒有增加，但是列數$j$增加，又兩條序列的中的元素不能重復比對，也不能比對序列其它的元素，所以就只能比對空格了。（向下延申同理。）
我們看兩個極端的例子：

##黑色路線
CATCCCCCAAA-----
-----------CAAAA
##紅色路線
-----CATCCCCCAAA
CAAAA-----------

你們說這第0行，0列的有啥用處，放這么多負分而產生極端比對。要排除這些不合適比對結果也要花費許多資源。

比對開始

上面進行了一次小探討。比對直接向右或向下會產生gap。因為行列有一方被固定，沒法提供可比對元素。如果要正確比對的話，要使得$i$和$j$同時增加，即向右下方延申。除了第0行和0列的元素，只有單一來源，矩陣其它的元素都有三個來源，即上方，左方，左上。如何確定$F(i,j)$的值，按照公式，選擇其中最大的值。$F(m,n)$($m,n$為兩條序列的長度)為最終序列比對的相似得分。
既然我們想使$F(m,n)$最大，肯定不能過多的產生空位罰分吧。那么我們從一開始就一直向右下方延申，直至不能繼續，能到得到最大的相似得分呢？我們試一下吧：

我們用箭頭表示從哪個方向得到的最大值，兩個箭頭則可以從兩個方向得到。從圖中紅色箭頭的路線可以看出來，一直選擇往右下延申的比對，似乎可以達到最大值。而且再比對三對元素就可以到達終點了。現在看紅色箭頭，似乎是全場全場最靚的仔。然而最終結果如何？我們繼續把剩余幾個各自補全看看。

現在再看全場最靚的仔似乎被綠了？？而且是一個綠一個的。簡直就是螳螂捕蟬，黃雀在后，然后漁翁得利。。。

貪心算法：紅色路線選擇的策略，從一開始，三個方向的選擇，每一步，它都選擇三個方向中得分最大的那個，但最后結果卻不是最好的。這種在每一步都做出當前看來是最好選擇的方法，而不考慮整體最優的方法，叫做貪心算法。這種方法在這里得不到最優結果。
窮舉搜索：這是一種很直觀的方法，就是把所有可能的情況都列舉出來，從中找最優的。但是用於序列比對是是非常差的一種方法。從矩陣左上方到達矩陣右下方，有多少種路線？矩陣每一個小格有三種可能路線，矩陣大小$m*n$,即有$3^{mn}$種路線。其中$m,n$稍微取大一點的值，數值就爆炸了。完全不可取。

emm，寫到這里，還沒講到全局比對是如何拆分問題的，現在我們就開始討論這個問題吧。。
首先，在矩陣中當比對到第$（m,n）$格時，代表着比對終止。因為兩條序列的每個元素都參加了比對過程，無論序列中的元素是比對到了空格，還是另一條序列的某個元素。因為序列中的元素都比對完了，在比對下去就只有空格對空格了，沒什么意義。所以$（m,n）$是終點。
我們就從終點開始探討如何分割大問題為小問題吧。
$$F(m,n)=max\begin{cases}F(m,n)+s(x_i,y_j)\F(m-1,n)+d\F(m,n-1)+d\end{cases}$$
其中$s(x_i,y_j)$和$d$是知道的,那問題$F(m,n)$則可以通過求出$F（m-1,n-1）$,$F（m,n-1）$和$F（m-1,n）$的分數得到。即一個大問題，分成了3個子問題了。三個子問題也不知道分數，就繼續分子問題唄。如圖：

分治算法：分治的思想也是將一個問題分為幾個問題來解決。這與動態規划的分離子問題有啥區別？從上面的圖我們可以看出來由$F(m,n)$分割的子問題們是存在重復問題的。分治分離的子問題一般獨立的不會相互重疊，而動態規划的子問題則一般發生重疊現象。

兩個性質

之前空着沒講的兩個性質，現在來講一下。一個是最優子結構，一個是無后效性
最優子結構：如果一個問題的最優解，可以將其分為幾個子問題，然后能從子問題的最優解推出大問題的最優解。就稱其有最優子結構性質。這一性質很容易從函數公式中看出來。
無后效性：給定某一階段的狀態，這一階段以后的發展過程不受這階段以前各段狀態的影響（來自阮行止的回答）這話是什么意思呢，以我們序列比對為例。比如：當你確定了前面某個階段的比對，比對情況如下面的比對方案的前面一部分，這時比對分數已經確定了。無后效性就是前面的比對情況對后面的比對結果的分數不會有任何影響。前面部分的比對方案的任意變動也不會影響后面部分的結果。

##1
CAT  CCCCC-----AAA
CA-  ----------AAA
##2
CAT  CCCCCAAA-----
CA-  AA----------A
##3
CAT  CCCCCAAA-----
CA-  ----------AAA

無后效性，可以讓我們放心大膽的在矩陣每個格子存放最優的結果，舍棄許多不是最優的分數。因為不存在前面階段的最差比對方案，還能與后面階段的最差方案組合成為全局最優比對方案。

復雜度

我們來看一下Needleman–Wunsch算法的復雜度。從二維矩陣中，可以看到我們需要計算每一格的所占最優分值，同時也為了之后的計算方便也得儲存每一個的分值。故時間復雜度和空間復雜度都是$O(mn)$。

其它

為了理解動態規划，我在知乎上這個問題什么是動態規划（Dynamic Programming）？動態規划的意義是什么？的答案下看了許多回答。許多答主都給了自己的理解。他們的回答有相似之處，也有不同。不能說誰對誰錯（小聲bb:也不敢說誰對誰錯）
就全局比對我也談談我的看法：
動態規划和其它方法有許多相似的地方。要是將其它方法的定義定廣一點，也能把動態規划說成其它方法。比如：
窮舉搜索： 之前上文曾提到過，不過那時我們說要窮舉的對象是每一個比對方案。現在我們換一個觀點想一想：窮舉每一個矩陣中的格子，即窮舉每一個最優的$F(i,j)$。這里我們的對象就從比對路徑而換成了$F(i,j)$（序列比對到第$(i,j)$的最大相似度）
貪心算法: 之前提到的貪心算法，我們狹隘的將最優，認為是向下、向右或向右下而能得到的分數。現在我們將最優看成相鄰$F(i,j)$中最大的那個。即按照這個貪心選取，它會可能不斷的對比相鄰的$F(i,j)$而選擇最優。直至，到最后剩下三條路徑到達$F(m,n)$,選中最優的那個。
總之，不管黑貓白貓，能幫助解決問題就是好貓！選擇你認為最好的一種理解。在以后多次遇到動態規划后，在慢慢完善你的認知。