python 比对匹配_用Python从头实现Needleman-Wunsch序列比对算法

本文将介绍生物信息学中一个基本的序列比对算法，Needleman-Wunsch比对算法(以下简称NW算法)，并展示如何编写一个简单的Python程序来实现该算法。全文分为上下两个部分：上篇介绍并阐释NW算法的原理，下篇展示如何用Python实现。如果你对该算法不陌生并且仅对实现方法感兴趣，那么可以直接跳到下篇。

——————————————————— 上篇 ————————————————————

序列比对是生物信息学领域的经典问题，出现了很多算法，也衍生出了一大批生物信息学工具，有些已经逐渐淡出历史舞台，有些依然被频繁使用，比如MUSCLE，BLAST，BWA，Bowtie2等。

两序列间的比对(“一对一”)是序列比对话题中最基本的问题，也是本文将要探讨的主题，更复杂的比对问题都可以由此衍生而来，比如，“多对多”：多重序列比对；“一对多”：目标序列对数据库的比对；“多对一”：大量二代测序短序列比对到人类基因组上。

两序列比对，从解决问题的方向上来说，可下分为全局比对(Global Alignment)和局部比对(Local Alignment)两种。前者试图将两条序列从头到尾整体对齐，而后者则只是针对两条序列中高度相似的局部区域进行比对。这两种方式都各有特点，局部比对的意义在于，在实际研究中，我们往往更关注序列间高度保守的区域，因为这暗示着重要的生物学功能。同样的两条序列，采用全局比对和局部比对都可以，但得到的结果可能差别很大：同样的序列，不同的比对方式。上：全局比对；下：局部比对 ##图片来自网络##

NW算法就是一种经典的全局比对算法，由Saul B. Needleman和Christian D. Wunsch于1970年首次提出，当时是为了解决蛋白序列的比对问题，至今依然在使用。值得一提的是，同时期涌现了许多重要的比对算法，比如Smith-Waterman算法。但NW算法是最早的、也是最经典的一个，所以你可能会在很多生物信息学专著中见到它，比如：Richard Durbin 教授的这本 Biological Sequence AnalysisDavid Mount教授的 Bioinformatics: Sequence and Genome Analysis

下面开始介绍NW算法。

第一步，定一个打分策略。这是基础， 定义了两个碱基在各种比对情况下的分值(包括gap)，我们这里将会用到一个简单的：

若是两个碱基一样，即完美匹配，即 +1 分

若是两个碱基不一样，即错配，则 -1 分

若是在任一条链上开一个gap，则 -1 分

我们用 θ(a, b)函数来表示碱基 a 和 b 之间的打分。

第二步，计算分值矩阵。分值矩阵的概念是NW算法的重点，矩阵由格子组成，一个比对方案就是一条从左上角格子移动到右下角格子的路径，一个矩阵的路径包含了两条序列间所有可能的比对情况，我们要做的就是找出这其中的最优的路径，这其中的方法就是根据分值矩阵来搜索。这里我们假设想要比对的序列分别为seq1和seq2：

>seq1

TCATC

>seq2

TCATGGC

在计算矩阵的时候，需要在两条序列之前加上一个 “-”，表示gap。左上角的格子为起点，记为0。接下来看矩阵里剩下的格子，每个格子都由别的格子推算而来，用 C(i, j) 表示第 i 行，第 j 列的格子的得分，计算方法：

其实对应的就是从“上”，”对角线“，”左“ 三个方向分别计算分值，取最大值。第一行和第一列的意义有点特殊，表示序列的起始阶段比对到gap上，因此这里的打分就是开 1 个gap，扣1 分，gap越多，扣的越多。

接下来就是剩下的格子，这就需要比较三个方向计算出的最大值了。以下图为例，要计算T对T这个格子的分值，首先T对T自身是一个完美匹配，+1，从对角线而来就是0+1，得到1；从上面来就是-1-1，得-2；从左边来就是-1-1，得-2。显然取最大值1，同时我们知道了这个T-T位置的1是来自对角线的。取最大值1

剩下的所有格子都可以此类推，每到一个格子干两件事：计算分值，记录路径：红色标示可走的路径