# 基因序列比對演算法
\begin{gather*}Needleman-Wunsch Algorithm\\基因序列比對演算法\end{gather*}
一般而言，我們會把DNA和蛋白質分別看成是由4和20個英文字母所組成的序列或字串，因為他們分別是由4種核苷酸和20種胺基酸所組成的。對DNA而言，突變是非常平常的事情，也是自然的演化過程。藉由基因的突變，生物可以適應自然環境的改變。

## DNA突變的類型
常見的DNA突變有3種，分別是取代、插入及刪除。
|DNA突變行為|意義|
| ---- | -------------------------- |
|取代(配錯)|把一個字母用另一個字母取代|
|插入|在DNA序列的某一個位置插入或刪除一個新的或舊的字母|
|刪除|在DNA序列的某一個位置插入或刪除一個新的或舊的字母|

## 編輯距離
通常生物學家會利用所謂的**編輯距離**，來衡量兩條DNA序列之間的相異程度。生命總是朝著最短路徑進行演化，所以兩條序列之間的編輯距離被定義為：把其中一條序列編輯轉成另外一條序列，所需最少的編輯運算個數。兩條DNA序列之間的編輯距離越小，代表它們之間的相似程度越高。從演化的觀點來說，這意味著它們演化自同一個祖先（即所謂的同源），所以彼此間應該會有相似的結構及功能。

生物學家喜歡利用比對來衡量兩條序列之間的相似程度。
拿GACGGATAG和GATCGGAATAG這兩條DNA序列來說，乍看之下這兩條長度不同的DNA序列似乎不太相似。但是，當我們如下圖一般，把它們重疊在一起，並在第1條序列的第2個和第3個字母之間與第6個和第7個字母之間分別插入一個空白字，就可發現其實這兩條DNA序列還蠻相像的。

| 重疊前 | 重疊後 |
| ----------- | ------ |
| GACGGATAG | GA_CGGA_TAG |
| GATCGGAATAG | GATCGGAATAG |
這種序列重疊的方式，就稱為序列的比對。

我們可以在兩條序列的任意位置上插入一個或多個空白字，目的是讓相同或相似的字母能夠儘量對齊。但要特別注意的是，不能讓兩個插入的空白字對齊在一起，因為這樣對衡量序列之間的相似程度並無幫助。

因此字母之間對齊的方式就只有2種:
字母與字母的對齊、字母與空白字的對齊( 即所謂的開gap）

---
## 序列比對的演算法
我們可以利用動態規劃（dynamic pro-gramming）技巧，設計出一個簡單又有效率的演算法來解決兩條序列比對的問題

首先我們定義一些符號，以便能夠簡單及清楚地介紹兩條序列比對問題的動態規劃演算法

1. 令A＝a1a2...am和B＝b1b2...bn分別表示長度各為m和n的兩條序列
2. 令Align（i,j）表示兩條序列片段Ai＝a1a2...ai和Bj=b1b2...bj之間分數（或相似度）最高的對齊方式，其中i和j的範圍分別是1≦i≦m與1≦j≦n
3. 令S（i,j）表示Align（i,j）的分數

我們的目的是:
**求得兩條序列之間最佳的分數S（m,n）及對齊方式Align（m,n）**。

---

## 欄位的對齊
在這之前，先把焦點放在S（i,j）和Align（i,j）的身上，因為若有辦法對任意的i和j都能計算出S（i,j）和Align（i,j），自然就可以求得S（m,n）和Align（m,n）。首先，我們把Align（i,j）切成左右兩半：右半部只含Align（i,j）的最後一個欄位，左半部則包含剩下的欄位。在這種情況之下，S（i,j）等於左半部的分數加上右半部的分數。

不管Align（i,j）的對齊方式為何，最後一個欄位的對齊方式就只有3種：
1. 字母對齊字母（即ai對齊bj）
2. 字母對齊空白字（即ai對齊－）
3. 空白字對齊字母（即－對齊bj）

若是第1種情形，那麼右半部的分數等於ai對齊bj的分數（用δ（ai,bj）表示），而左半部的分數等於Ai－1＝a1a2...ai－1與bj－1＝b1b2...bj－1之間對齊的最高分數（即等於S（i-1,j-1））。換句話說，S（i,j）= S（i-1,j-1）+δ（ai,bj）。
若是第2種情形，同理可證，則S（i,j）= S（i-1,j）+δ（ai,－）
若是第3種情形，則S（i,j）= S（i,j-1）+δ（－,bj），其中δ（ai,－）和δ（－,bj）分別表示「ai對齊空白字」與「空白字對齊bj」的分數。

綜合這3種情形＠我們可以把S（i,j）表示成一個遞迴函式如下：


例如我們可以把計算S（m,n）的問題分解成計算S（m-1,n-1），S（m-1,n）和S（m,n-1）的3個小問題，然後再把這3個小問題的解代入遞迴函式，就可求出S（m,n）。同樣的道理，可以把計算S（m-1,n-1）的小問題，分解成計算S（m-2,n-2）、S（m-2,n-1）及S（m-1,n-2）這3個更小的問題。如此遞迴下去，我們可以把整個計算S（m,n）過程用一顆樹狀圖來描述。


## 為何使用動態配置
然而，這種由上而下計算S（m,n）的方式會十分耗時，因為整顆樹會有指數個節點，且每一個節點都需要花時間去計算其遞迴函式值。但是，若仔細觀察S（m,n）的樹狀圖，我們會發現有許多的節點是在做相同的事情。越接近樹的底層，重複的節點就越多。因此若有辦法對相同的節點只做一次計算，那麼整體的執行速度不就可以加快了嗎？

於是我們想到利用由下而上的方式計算S（m,n）：就是對所有相同的節點，只需對第1次遇到的節點做計算，然後把計算出來的值儲存起來，待下次碰到同樣的節點時，只要呼叫或參考這個節點的儲存值即可，如此便可省去重複計算的時間。

這個主意讓我們可以把上述的樹狀圖簡化成一個（m+1）×（n+1）的二維矩陣圖，其中座標（i, j）的節點等於S（i,j）。根據遞迴關係，S（i,j）只跟它相鄰的上、左上及左節點有關係，我們必須先知道這3個相鄰節點的分數之後才能夠求出S（i,j）

## 節點與節點的關係
其中節點之關係分為以下三種:
| 方向 | 關係式 | 意義 |
| ---- | -------- | -------- |
|直線|S（i,j）= S（i-1,j）+ δ（ai,－）|ai與空白字對齊|
|斜對線|S（i,j）= S（i-1,j-1）+ δ（ai,bj）|且ai和bj對齊|
|橫線|S（i,j）= S（i,j-1）+ δ（－,bj）|空白字和b對齊|

為了求得S（m,n），我們可以採用以列為導向的方式計算矩陣中每一個節點分數：就是以由左到右的方式先計算出第1列節點的分數（開始時令S（0,0）= 0），然後第2列、第3列......一直到求出最後一列的S（m,n）為止。當然也可以採取以欄或斜對角線為導向的方式，計算矩陣中節點的分數。換句話說，整個計算S（m,n）的過程就好像在表格（即二維矩陣）上填寫數字（即分數）一般。所以，動態規劃的演算法又被稱為填表格的方法。

---

## 動態規劃表
|
| |_|A|C|G|T|A|T|A|
|-|-|-|-|-|-|-|-|-|
|_|0|-1|-2|-3|-4|-5|-6|-7|
|A|-1|1|0|-1|-2|-3|-4|-5|
|G|-2|0|0|1|0|-1|-2|-3|
|G|-3|-1|-1|1|0|-1|-2|-3|
|A|-4|-2|-2|0|0|1|0|-1|
|C|-5|-3|-1|-1|-1|0|0|-1|
|T|-6|-4|-2|-2|0|-1|1|0|
|A|-7|-5|-3|-3|-1|1|0|2|


δ分數配置
| 配對正確 | 配對錯誤 | 填補gap(空白) |
| -------- | -------- | -------- |
| 1分|-1分|-1分|

以比對AGGACTA與ACGTATA這兩條DNA序列的二維矩陣為例,採用的配分函數是:配對給1分,配錯或開gap給-1分。除此之外,也利用箭頭表示每一個節點的分數是透過哪一個相鄰節點所得到的。從最右下角的節點分數,可以知道AGGACTA與ACGTATA 之間最佳的比對分數是2分。但是這兩條DNA序列之間的字母該怎麼對齊才能得到這個最佳分數呢?答案是利用在二矩陣上所留下的箭頭,從最右下角的節點(7,7)開始,沿著捐進來的箭頭逆推回至最左上角的飾點(0,0).如此,可以得到一條(甚至數條)由(0.0)到(7,7)的最短路徑,然後再根據這條路徑上每一箭頭的方向,寫出相對應字母之間對齊的方式,最後再把這些對齊的字母合併起來即可

---

不難發現，我們總共有（m+1）×（n+1）個節點要計算，而且每一個節點計算最多只花3個加法運算及2次的比較運算。因此，整個計算過程所花的運算時間等於5×（m+1）×（n+1），其中m和n是影響整個時間複雜度最關鍵的因素，所以整個時間複雜度又被簡寫成O（m×n），意味著時間複雜度與兩條序列長度乘積成正比。

## 資料來源
專題報導:生物資訊:基因序列比對的演算法動態規劃動態規劃
(交通大學生物資訊研究所-盧錦隆)
http://ezphysics.nchu.edu.tw/15_course/Bioinformation/RF05.pdf
Sequence Alignment Algorithms
(Notes from Dr. Bino John and Dr. Takis Benos )
https://www.cs.cmu.edu/~02710/Lectures/SeqAlign2015.pdf
計算生物學-Sequence Alignment and Dynamic Programming投影片
https://ocw.mit.edu/courses/electrical-engineering-and-computer-science/6-096-algorithms-for-computational-biology-spring-2005/lecture-notes/lecture5_newest.pdf
維基百科:Needleman–Wunsch_algorithm
https://en.wikipedia.org/wiki/Needleman%E2%80%93Wunsch_algorithm

---
## 演算法實驗課作業
Questions:
1. 在一般的實作當中，f（m，n）可以找到最佳對齊方式之一。 如何找出所有最佳比對？
In normally implementation, f(m,n) can find one of optimal alignment. How to find out all of optimal alignments?

(找任意最佳解、找所有最佳解02)建構table算法

d ← Gap penalty score
for i = 0 to length(A)
    F(i,0) ← d * i
for j = 0 to length(B)
    F(0,j) ← d * j
for i = 1 to length(A)
    for j = 1 to length(B)
    {
        Match ← F(i−1, j−1) + S(Ai, Bj)
        Delete ← F(i−1, j) + d
        Insert ← F(i, j−1) + d
        F(i,j) ← max(Match, Insert, Delete)
    }

(找任意最佳解)backtrack算法

AlignmentA ← ""
AlignmentB ← ""
i ← length(A)
j ← length(B)
while (i > 0 or j > 0)
{
    if (i > 0 and j > 0 and F(i, j) == F(i−1, j−1) + S(Ai, Bj))
    {
        AlignmentA ← Ai + AlignmentA
        AlignmentB ← Bj + AlignmentB
        i ← i − 1
        j ← j − 1
    }
    else if (i > 0 and F(i, j) == F(i−1, j) + d)
    {
        AlignmentA ← Ai + AlignmentA
        AlignmentB ← "−" + AlignmentB
        i ← i − 1
    }
    else
    {
        AlignmentA ← "−" + AlignmentA
        AlignmentB ← Bj + AlignmentB
        j ← j − 1
    }
}

(找所有最佳解02)backtrack算法
(以遞歸形式重寫backtrack算法)

i ← length(A)
j ← length(B)
backtrack(i,j){
    if (i == 0 and j == 0)
    {
            AlignmentA ← ""
            AlignmentB ← ""
    }
    //改成if
    if (i > 0 and j > 0 and F(i, j) == F(i−1, j−1) + S(Ai, Bj))
    {
            AlignmentA ← Ai + AlignmentA
            AlignmentB ← Bj + AlignmentB
            i ← i − 1
            j ← j − 1
            backtrack(i,j);
            //補回i和j值
            i++;
            j++;
    }
    //改成if
    if (i > 0 and F(i, j) == F(i−1, j) + d)
    {
        AlignmentA ← Ai + AlignmentA
        AlignmentB ← "−" + AlignmentB
        i ← i − 1
        backtrack(i,j);
        //補回i值
        i++;
    }
    //改成if
    if
    {
        AlignmentA ← "−" + AlignmentA
        AlignmentB ← Bj + AlignmentB
        j ← j − 1
        backtrack(i,j);
    }
    //補回j值
    j++;
}

(找所有最佳解01)建構table算法

d ← Gap penalty score
for i = 0 to length(A)
    F(i,0) ← d * i
for j = 0 to length(B)
    F(0,j) ← d * j
for i = 1 to length(A)
    for j = 1 to length(B)
    {
        Match ← F(i−1, j−1) + S(Ai, Bj)
        Delete ← F(i−1, j) + d
        Insert ← F(i, j−1) + d
        passValue(stackMatrix[i],[j],i,j, Match, Delete，Insert)
    }
passValue(stack<int>& Stack,int i,int j,int M,int I,int D){
    sameChar ← (Ai == Bj)
    maxValue ← max(M,I,D)
    F(i,j) ← maxValue
    if(maxValue == M)//往左上走
        Stack.push(1)
    if(maxValue == I)//往上走
        Stack.push(2)
    if(maxValue == D)//往左走
        Stack.push(3)
}
stack<int> stackMatrix[length(A)][length(B)]

(找所有最佳解01)backtrack算法

stack<char> AlignmentA ← ""
stack<char> AlignmentB ← ""
void backtrack(i,j){
    if(i == 0 and j == 0){
        output AlignmentA
        output AlignmentB
        return
    }
    else{
        while(!stackMatrix(i,j).empty()){
            switch(stackMatrix(i,j).pop()){
                case(3):
                    AlignmentA.push(Ai)
                    AlignmentB.push('-')
                    backtrack(i,j-1)//向左走
                    AlignmentA.pop()
                    AlignmentB.pop()
                case(2):
                    AlignmentA.push('-')
                    AlignmentB.push(Bj)
                    backtrack(i-1,j)//向上走
                    AlignmentA.pop()
                    AlignmentB.pop()
                case(1):
                    AlignmentA.push(Ai)
                    AlignmentB.push(Bj)
                    backtrack(i-1,j-1)//向左上走
                    AlignmentA.pop()
                    AlignmentB.pop()
            }
        }
        return
    }
}

2. 可能存在多少個最佳比對？ 請構建一組輸入以解釋您的答案。
How many optimal alignments may exist? Please construct a set of input to explain your answer.
最多$3^{mn}$個

| |-|T|T|T|
|-|-|-|-|-|
|-|0|0|0|0|
|T|0|1|1|1|
|T|0|1|2|2|
|T|0|1|2|3|

3. 假設A和B都很長，我們無法在內存中儲存所有m×n個值。 請找出一個只需緩存n個值的方法。
Suppose both A and B are very long, that we can’t maintain all m×n scores in memory. Please find the way which only cache n values.
4. 分析Q1和Q2在最佳情況和最壞情況下的空間複雜度，時間複雜度。
Analyze space complexity, time complexity in best case and worst case in Q1 and Q2.

1.時間複雜度:因為每一個方向都至多有3種backtrack方向，而這樣的運算在沒有DP的情況下最壞會持續m*n次，因此，複雜度為O($3^{mn}$)
2.空間複雜度:因為只需要儲存m列n行，複雜度為O(mn)
5. 完成 Solve http://oj.csie.ndhu.edu.tw/problem/ALG04C.

#include <bits/stdc++.h>

using namespace std;

int main()
{
    string a, b;
    bool first = true;
    while(getline(cin,a)){
        getline(cin,b);
        if(!first) cout << endl;
        first = false;
        vector<vector<int>> dp(a.length()+1,vector<int>(b.length()+1));
        set<string> res;
        int best = 0;
        // longest common substring dp
        for(int i=0;i<a.length();i++){
            for(int j=0;j<b.length();j++){
                if(a[i] == b[j]){
                    dp[i+1][j+1] = dp[i][j]+1;
                    if(dp[i+1][j+1] > best){
                        best = dp[i+1][j+1];
                        res.clear();
                    }
                    if(dp[i+1][j+1] == best){
                        res.insert(a.substr(i-best+1,best));
                    }
                }
            }
        }
        if(res.empty()) cout << "No common sequence." << endl;
        else {
            for(auto& str : res)
                cout << str << endl;
        }
        getline(cin, a); //junk
    }
}