尼德曼-翁施演算法

尼德曼-翁施演算法（Needleman-Wunsch Algorithm）是基於生物資訊學的知識來匹配蛋白序列或者DNA序列的演算法。這是將動態演算法應用於生物序列的比較的最早期的幾個實例之一。該演算法是由 Saul B. Needlman和 Christian D. Wunsch 兩位科學家於1970年發明的。本演算法高效地解決了如何將一個龐大的數學問題分解為一系列小問題，並且從一系列小問題的解決方法重建大問題的解決方法的過程。該演算法也被稱為優化匹配演算法和整體序列比較法。時至今日尼德曼-翁施演算法仍然被廣泛應用於優化整體序列比較中。

## 新手指南

本演算法可以用於比較任意兩個序列。本指南將會使用兩個短小的DNA序列作為例子圖1來演示該演算法

```
GCATGCU
GATTACA
```

### 初始化得分矩陣

首先建立如圖一的得分矩陣。從第一列第一行的位置起始。計算過程從d 0 , 0開始，可以是按行計算，每行從左到右，也可以是按列計算，每列從上到下。當然，任何計算過程，只要滿足在計算d i , j時d i-1 , j、d i-1 , j-1和d i, j-1都已經被計算這個條件即可。在計算d i , j後，需要保存d i , j是從d i-1 , j、d i-1 , j-1或di, j-1中的哪一個推進的，或保存計算的路徑，以便於後續處理。上述計算過程到d m , n結束。

### 選擇得分體系

接下來我們需要決定如何給每個獨立配對打分。從我們的序列中，你可能就能發現最好的序列配對之一：

```
GCATG-CU
G-ATTACA
```

我們可以看出每個位置配對都有三種可能情況：匹配，不匹配與錯位（或插入）：

- 匹配：兩個字母相同
- 不匹配：兩個字母不相同
- 錯位：一個字母與另一個序列中的間隔（空白）相匹配

這三種得分情況有很多打分標準，這些情況都總結在得分體系的部分。從現在開始，我們將使用Needleman 和Wunsch創造的簡單體系來進行打分，即匹配得一分，不匹配得-1分，錯位得-1分.

### 填充得分矩陣

開始於第二行中的第二列，初始值為0。通過矩陣一行一行移動，計算每個位置的分數。得分被計算為從現有的分數可能的最佳得分（即最高）的左側，頂部或左上方（對角線）。當一個得分從頂部計算，或從左邊這代表在我們的對準的插入缺失。當我們從對角線計算分數這表示兩個字母所得位置匹配的對準。定不存在「向上」或「左上」的位置對第二行，我們只能從現有單元向左添加。因此，我們添加-1的權利，因為這代表了從以前的比分插入缺失。這導致在第一行是0，-1，-2，-3，-4，-5，-6，-7。這同樣適用於第二列，因為我們只具有以上現有分數。因此，我們有：

G｜C｜A｜T｜G｜C｜U

0｜-1｜-2｜-3｜-4｜-5｜-6｜-7

G｜-1

A｜-2

T｜-3

T｜-4

A｜-5

C｜-6

A｜-7

第一種情況是存在向三個方向構築矩陣，周圍位置得分如圖

G

0｜-1

G｜-1｜?

對接下來的位置，我們有不同的選擇。

G｜C

-1｜-2

G｜1｜0

C｜A

A｜0｜1

T｜-1｜0

在某些位置，可能有兩個或者三個方向得分相同，如圖所示：

T｜G

T｜1｜1

A｜0｜0

完整的得分矩陣如下圖所示

所以根據該矩陣就可以方便的出我們想要的匹配路徑了。

```
结果：

序列         最佳匹配
GCATGCU      GCATG-CU      GCA-TGCU      GCAT-GCU
GATTACA      G-ATTACA      G-ATTACA      G-ATTACA

解读：

你可以从最左位置开始解读匹配

匹配:  +GCATGCU
      +GATTACA
分数： 0 没有任何分数<div><br>
</div><div>匹配：  +GCATGCU
       +GATTACA
分数： 0 1 间隔 分数 -1</div><div>匹配：   +GCATGCU
       +GATTACA</div><div>分数： 0 2 间隔 分数 -2</div><div>匹配:     +GCATGCU
       +GATTACA
分数:  0  3 间隔 分数 -3

以此类推
</div>
```

### 追溯本源

注意：可能有多個路徑都能得到最高的分數，在這裡我們只展示其中一種

根據指向上一個位置的箭頭來回溯取得最佳匹配的路徑。然後根據從起始位置指向最終位置的路徑來完成匹配的建立。方法如下：

- 箭頭表示匹配或者不匹配，不管哪一種，都表明這兩個位置上的字母相關聯。如以下例子所示：

```
空--> G
空 --> G
```

- 如果有一個水平的箭頭表示會有一個縱向序列的間隔。如下圖例子所示：

```
G --> GC
G --> G-
```

- 如果有一個垂直的箭頭表示會有一個橫向序列的間隔。如下圖例子所示：

```
GCA --> GCA-
G-A --> G-AT
```

- 注意可能會有多個箭頭指向同一個位置，這代表了匹配的路徑有很多不同的分支，如果這些路徑的得分相同，代表這些匹配方法實質上是等效的

根據以上規則我們可以得出其中一種可能的匹配路徑如下:

```
G →  GC → GCA → GCA- → GCA-T → GCA-TG → GCA-TGC →  '''GCA-TGCU'''
G →  G- → G-A → G-AT → G-ATT → G-ATTA → G-ATTAC →  '''G-ATTACA'''
```

## 得分體系

### 基礎得分標準

最簡單的打分體系只考慮匹配，不匹配與間隔三種情況。在新手指南中使用的打分體系是匹配=1，不匹配=-1，間隔=-1。在這種情況下序列越長，得分越低，在本體系中我們期望好的匹配取得高分。另一種打分體系如下：

- 匹配=0
- 間隔=1
- 不匹配=1

在本體系中得分將會代表兩序列中匹配的距離。不同的得分體系適用於不同的情況。例如，如果在你的匹配中間隔被認為有很大的負面影響，你可以使用間隔扣分很多的打分體系：

- 匹配=0
- 間隔=10
- 不匹配=1

### 相似度矩陣

更複雜的打分體系不僅僅考慮選擇的不同類型，也考慮每個不同的字母的影響。例如，字母A與A匹配將得到4分。為了表示所有字母匹配得分的情況，我們將會使用相似度矩陣來幫助計算

A｜G｜C｜T

A｜1｜-1｜-1｜-1

G｜-1｜1｜-1｜-1

C｜-1｜-1｜1｜-1

T｜-1｜-1｜-1｜1

每個分數代表從一個字母轉換到下一個位置的所有可能的情況的得分。

A｜G｜C｜T

A｜1｜-1｜-1｜-1

G｜-1｜1｜-1｜-1

C｜-1｜-1｜1｜-1

T｜-1｜-1｜-1｜4

不同的得分矩陣適用於不同的情況。該得分體系對於蛋白序列的匹配比對尤為重要。以下為兩種廣為引用的蛋白序列相似度矩陣

- PAM
- BLOSUM

### 間隔補償

當比對序列往往存在差距（即插入缺失），有時路數。生物學上，有很大的差距更容易，而不是多個單一缺失的發生為一體的大型刪除。因此，我們應該進球兩個小的插入缺失是有過之而無不及一個大的。簡單和常見的方式做，這是通過一個大的差距，比分開始一個新的插入缺失和一個較小的間隙延伸得分為每個擴展了插入缺失信。例如，新插入缺失可能會花費-5和擴展，插入缺失可能會花費-1。以這種方式取向，例如：

```
GAAAAAAT
G--A-A-T
```

該序列有多種匹配方式，類似的匹配方式有：

```
GAAAAAAT
GAA----T
```

或者任意有四個間隔的匹配方式

## 進階矩陣

例如，如果相似度矩陣如下圖所示：

A｜G｜C｜T

A｜10｜-1｜-3｜-4

G｜-1｜7｜-5｜-3

C｜-3｜-5｜9｜0

T｜-4｜-3｜0｜8

那麼匹配結果如下：

```
AGACTAGTTAC
CGA---GACGT
```

如果間隔懲罰因子是-5，那麼該匹配的得分為：

= -3 + 7 + 10 - (3 × 5) + 7 + (-4) + 0 + (-1) + 0 = 1

計算序列相似程度時還應該考慮序列長度的影響。令S(s, t) 表示兩條長度分別為m和n的序列的相似性得分，假設S(s, t) = 99，兩條序列有99個字符一致。如果 m=n=100，則可以說這兩條序列非常相似，幾乎一樣。但是，如果m=n=1000，則僅有10% 的字符相同。所以，在實際序列比較時，使用相對長度的得分就更有意義，定義：

（4）

以sim（s，t）作為衡量序列相似性的指標。

從所使用的數據結構d i , j本身及其計算過程來看，序列兩兩比對基本演算法的空間複雜度和時間複雜度都是O (mn)。如果兩條序列的長度相近，空間和時間複雜度就為O (n2)。

首先將要看到如何運用動態編程查找兩個 DNA 序列的 最長公共子序列（longest common subsequence，LCS）。發現了新的基因序列的生物學家通常想知道該基因序列與其他哪個序列最相似。查找 LCS 是計算兩個序列相似程度的一種方法：LCS 越長，兩個序列越相似。

子序列中的字符與子字串中的字符不同，它們不需要是連續的。例如，ACE是 ABCDE的子序列，但不是它的子字串。請看下面兩個 DNA 序列：

•S1= GCCCTAGCG

•S2= GCGCAATG

$$F_{ij} = \max(F_{i-1,j-1} + S(A_{i}, B_{j}), \; F_{i,j-1} + d, \; F_{i-1,j} + d)$$

計算程序如下

```
'''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(A<sub>i</sub>, B<sub>j</sub>)
    Delete ← F(i-1, j) + d
    Insert ← F(i, j-1) + d
    F(i,j) ← '''max'''(Match, Insert, Delete)
  }
```

這兩個序列的 LCS 是 GCCAG。（請注意，這僅是 一個LCS，而不是 唯一的LCS，因為可能存在其他長度相同的公共子序列。這種最優化問題和其他最優化問題的解可能不止一個。）

```
AlignmentA &#x2190; ""
AlignmentB &#x2190; ""
i &#x2190; '''length'''(A)
j &#x2190; '''length'''(B)
'''while''' (i > 0 '''or''' j > 0)
{
  '''if''' (i > 0 '''and''' j > 0 '''and''' F(i,j) == F(i-1,j-1) + S(A<sub>i</sub>, B<sub>j</sub>))
  {
    AlignmentA &#x2190; A<sub>i</sub> + AlignmentA
    AlignmentB &#x2190; B<sub>j</sub> + AlignmentB
    i &#x2190; i - 1
    j &#x2190; j - 1
  }
  '''else''' '''if''' (i > 0 '''and''' F(i,j) == F(i-1,j) + d)
  {
    AlignmentA &#x2190; A<sub>i</sub> + AlignmentA
    AlignmentB &#x2190; "-" + AlignmentB
    i &#x2190; i - 1
  }
  '''else'''
  {
    AlignmentA &#x2190; "-" + AlignmentA
    AlignmentB &#x2190; B<sub>j</sub> + AlignmentB
    j &#x2190; j - 1
  }
}
```

## 歷史與演算法發展

通過Needleman和Wunsch的描述的演算法的最初目的是找到在兩種蛋白質的胺基酸序列的相似性。

Needleman和Wunsch的描述他們的演算法為明確的情況下，當對準被匹配和不匹配僅僅懲罰，並差距沒有懲罰（D =0）。從1970年最初發布提示遞歸。$F_{ij} = \max_{h<i,k<j} \{ F_{h,j-1}+S(A_{i},B_{j}), F_{i-1,k}+S(A_i,B_j) \}$.

對應的動態規劃演算法需要立方時間。該報還指出，遞歸可容納任意差距處罰公式：

懲罰因子，減去對於由每個間隙，可被評估為阻擋在允許的間隙。懲罰因子可以是間隙的大小和/或方向的函數。頁444

一個更好的動態規劃演算法二次運行時間為同樣的問題（無間隙罰款）最早是由大衛Sankoff1972年類似的二次時間演算法推出了由TK Vintsyuk於1968年獨立發現了一種語音處理（「時間扭曲」），由羅伯特·華格納和麥可·J·菲舍爾在1974年的字串匹配。

Needleman和Wunsch的配製最大化相似性的問題。另一種可能性是最小化序列之間的編輯距離，弗拉基米爾的Levenshtein引入。彼得·塞勒斯在1974年發現的兩個問題是等價的。

Needleman-Wunsch演算法仍是廣泛使用的最佳全局比對，特別是當全局比對的質量即結果是最重要的。然而，該演算法是相對於時間和空間的昂貴，正比於兩個序列，因而不適合長度特別長的序列比對。

全局性配對和局部性配對的比較全局性配對和局部性配對的基本區別是在局部性配對中，你將嘗試用你訴求的序列的片段或局部來進行配對。而在全局性配對中，你將嘗試使用整個序列（從起始端到結束端）來進行配對，也即是說，在全局性配對中，當你的比對序列長度不一致時，在配對過程中可能會產生間隔或者不匹配的情況。注意，在局部性配對中也有可能產生間隔。

局部性配對

<image>

全局性配對

<image>

以下是些基於兩種演算法不同點的比較。

在Needleman-Wunsch（全局性）演算法中，分數跟蹤系統是基於得分矩陣的右下角開始的。然而在Smith-Waterman（局部性）演算法中，得分矩陣的起始點是基於分數最高的坐標位置。

全局性配對主要應用於比較同源性的基因的相似性而局部性配對主要應用於在非同源性的基因序列中尋找具有相似性的同源的基因域。

全局性配對是當你將全部的序列納入考慮範圍之內後進行的序列配對，而局部性配對則是將序列的一小部分進行優先比對。以下是為了幫助理解該概念而舉得一個例子：

假設你有一個較長的序列作為參考，假設該序列長為2000鹼基。你所要配對的序列長度較短，約為100鹼基左右。假設在你的參考序列中有與配對序列極為接近的鹼基序列片段出現，當你進行局部性配對時，你將得到一個非常好的匹配結果和相當高的匹配分數。然而，當你進行全局性配對時，匹配程度將會很低。因為該配對方法要求將兩個序列的所有鹼基都進行配對，所以當兩序列長度相差很大時，將會造成配對結果含有很長的間隔。儘管在某處你的配對序列和參考序列的匹配程度相當高，在全局性配對時，將不會把這兩個片段直接進行配對，而是儘量嘗試讓參考序列內的所有鹼基都參與匹配，造成結果與局部性配對方法有極大差別。

因此，在選擇配對演算法的時候，必須參考兩序列的具體長度和配對的實際意義和需要，選擇適當的方法，才能得到有生物學意義的結果。

## 程序代碼

本演算法可以在多個平台上實現運算，常見的程式語言均有公開的代碼可以參考使用。

以下提供在C語言中演算法的實現usage: ./bin/needleman_wunsch OPTIONS seq1 seq2

```
     Needleman-Wunsch optimal global alignment (maximises score).  
     Takes a pair of sequences on the command line, or can read from a
     file and from sequence piped in.  Can read gzip files, FASTA and FASTQ.

     OPTIONS:
       --file <file>        同时读取比对的两序列文件
       --files <f1> <f2>    
       --stdin              读取标准文件头

       --case_sensitive     使用大小写字母区分
       --match <score>      [默认: 1]
       --mismatch <score>   [默认: -2]
       --gapopen <score>    [默认: -4]
       --gapextend <score>  [默认: -1]

       --scoring <PAM30|PAM70|BLOSUM80|BLOSUM62>
       --substitution_matrix <file>  
       --substitution_pairs <file>   

       freestartgap       间隔不扣分
       freeendgap         序列结束不扣分

       printscores        打印优化后的比对分数
       zam               
       printmatrices      打印动态矩阵
       printfasta         打印fasta格式的文件头
       pretty             打印描述行
       colour             打印颜色

     可选优化项：
       nogapsin1          被比对序列不允许有间隔
       nogapsin2          比对序列不允许有间隔
       nogaps             所有序列均不允许有间隔
       nomismatches       不允许有误差配对
```

## 使用Needleman-Wunsch演算法的工具

- EMBOSS Needle and EMBOSS Stretcher Global Alignment Tools
- Needleman-Wunsch alignment for two nucleotide sequences
- MathWorks - Globally align two sequences using Needleman-Wunsch algorithm

應用——如何使用該演算法生物資訊學國家中心網站提供了包括各種生物資訊學計算工具的優秀資源，其中就有進行全局性配對計算的網頁連結。

NCBI 全局性配對

1 進入NCBI 主頁 (ncbi.nlm.nih.gov)

2 尋找Popular Resources, 點擊Blast

- 選擇Needleman-Wunsch演算法你可以選擇自行輸入任意序列對，下滑頁面至底部並點擊Align，比對的結果將提供鹼基的配對結果，相似的百分比和最終得分等資訊。

或者，你可以選擇搜索基因名稱來選擇要比對的序列。

1 進入NCBI 主頁2 尋找到Popular Resources，點擊Gene 3 輸入你要搜索的基因名稱和種族名稱4 點擊結果中顯示的基因名稱5 選擇Needleman-Wunsch演算法並運行

## 生物資訊學以外的應用

### 三維立體畫

三維立體畫是利用人眼立體視覺現象製作的繪畫作品。普通繪畫和攝影作品，包括電腦製作的三維動畫，只是運用了人眼對光影、明暗、虛實的感覺得到立體的感覺，而沒有利用雙眼的立體視覺，一隻眼看和兩隻眼看都是一樣的。充分利用雙眼立體視覺的立體畫，將使你看到一個精彩的世界。人有兩隻眼，兩隻眼有一定距離，這就造成物體的影象在兩眼中有一些差異，見右圖，由圖可見，由於物體與眼的距離不同，兩眼的視角會有所不同，由於視角的不同所看到是影象也會有一些差異，大腦會根據這種差異感覺到立體的景象。三維立體畫就是利用這個原理，在水平方向生成一系列重複的圖案，當這些圖案在兩隻眼中重合時，就看到了立體的影象。參見下圖，這是一幅不能再簡單的立體畫了。圖中最上一行圓最遠，最下一行圓最近，請注意：最上一行圓之間距離最大，最下一行圓之間距離最小。重複圖案的距離決定了立體影象的遠近，生成三維立體畫的程序就是根據這個原理，依據三維影象的遠近，生成不同距離的重複圖案。

立體畫有兩種形式：第一種是由相同的圖案在水平方向以不同間隔排列而成，看起來是遠近不同的物體。另一種立體畫較複雜，在這種立體畫上你不能直接看到物體的形象，畫面上只有雜亂的圖案，製作這樣的立體畫只有使用程序了。兩種作品看法是一樣的，原理都是使左眼看到左眼的影象，讓右眼看到右眼的影象。

### 電腦立體視覺效果

立體匹配是從一對立體圖像的三維重建進程的一個重要步驟。當圖像已被糾正，打個比方可以對準核苷酸/蛋白序列和匹配像素屬於掃描線，因為這兩個任務，著眼於人物的兩個字串之間建立最佳的對應關係繪製。實際上，一對立體聲的「正確」的圖像可以被看作是「左」圖像的突變版本：噪聲和單個攝像機的靈敏度改變的像素值（即字符取代）;和不同的視角揭示以前閉塞的數據，並引入新的閉塞（即插入和字符刪除）。作為結果，中的Needleman-Wunsch演算法的小的修改，使其適合於立體匹配。儘管該演算法在準確性方面性能都不算是最先進的，該但演算法的相對簡單性允許其在嵌入式系統上執行，使其有廣泛的應用.

.

## 其他連接

- Smith-Waterman演算法
- 序列比較

## 外部連結

- NW-align: A protein sequence-to-sequence alignment program by Needleman-Wunsch algorithm (online server & source code)
- Needleman-Wunsch Algorithm as Haskell Code
- A live Javascript-based demo of Needleman-Wunsch
- An interactive Javascript-based visual explanation of Needleman-Wunsch Algorithm
- B.A.B.A. — an applet (with source) which visually explains the algorithm.
- A clear explanation of NW and its applications to sequence alignment
- Sequence Alignment Techniques at Technology Blog
- Biostrings R package implementing Needleman–Wunsch algorithm among others
- Parallel Needleman-Wunsch Algorithm for Grid Implementation by Tahir Naveed, Imitaz Saeed Siddiqui and Shaftab Ahmed - Bahria University
- Needleman-Wunsch Algorithm as Haxe Code

## 基本資料

- 車輛類別：序列比對
- 時間：$O(mn)$
- 容積、佔地面積：$O(mn)$