尼德曼-翁施算法

尼德曼-翁施算法（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)$