生物信息学（1）——双序列比对之Needleman-Wunsch(NW)算法详解及C++实现

生物信息学系列博客索引

生物信息学（1）——双序列比对之Needleman-Wunsch(NW)算法详解及C++实现
生物信息学（2）——双序列比对之Smith-Waterman(SW)算法详解
生物信息学（3）——双序列比对之BLAST算法简介
生物信息学（4）——多序列比对之CLUSTAL算法详解及C++实现
生物信息学（5）——基于CUDA的多序列比对并行算法的设计与代码实现
项目gitee地址，内附源码、论文等文档信息

1. 什么是序列比对

所谓的序列比对，就是两个或者多个序列按照碱基排列进行比较，从而反映片段之间的相似性和阐明序列的同源性。这里主要是将未知功能的序列与已知序列进行比对，从而确定序列分析。序列比对的基本思想是，基于生物学中序列决定结构，结构决定功能的普遍规律，将核酸序列和蛋白质一级结构上的序列都看成由基本字符组成的字符串，检测序列之间的相似性，发现生物序列中的功能、结构和进化的信息

2. 引入序列比对的原因

对于两种相似的序列，DNA复制一共有三种情况可能导致两个序列不同:
(1)SNP(单核苷酸多态性)，即为碱基的替换，这是出现频率最高的，如将ACCT复制成了ACGT;
(2)INSERT，即为多复制一个碱基，如将AGCT复制成了AAGCT;
(3)DELETION，即为少复制了碱基，如将AGCT复制成了AGC。
由于DNA复制有三种可能会出现的出错情况，属于相同的种类的基因并不完全相同，即序列比对算法并不是简单的检测序列完全相等的情况，而是综合评判两个序列的相似程度，则就需要序列比对来对两个序列进行同源检测或遗传规律的检测。基于有这三种出错情况，那么两个序列对应的存在可能性一共有三种:
(1)MATCH:上下匹配;
(2)MISMATCH:出现SNP，上下不匹配;
(3)INDEL:出现INSERT或DELETION，导致两个序列有一个出现了空缺。
通常情况下，MATCH表示一个正确的匹配，而MISMATCH表示错误的匹配，INDEL相对于MISMATCH，虽然也是错误匹配，但由于INDEL一位，后面的碱基对则全部会错，所以INDEL错误比MISMATCH错误具有更大的破坏性，在序列比对中，INDEL错误会比MISMATCH有更高的罚分。

3. Needleman-Wunsch(NW)算法详解

3.1 公式

Needleman-wunsch算法是双序列比对算法里最经典的算法，其创新的将动态规划的思想引入到生物信息学中，开辟了一个新的纪元。Needleman-wunsch算法定义三种情况，MATCH、DISMATCH与INDEL，每一种对应不同的罚分，NW算法罚分规则如下:

以罚分规则为基础，通过以下公式构造一个打分矩阵:

最后通过回溯得分矩阵，获取全局最优解，回溯的规则是单元从哪里来，就回溯到哪里去。这样说很抽象，接下来举一个例子。

3.2 举个栗子

给定序列

Seq1 = GGATCGA

Seq2 = GAATTCAGTTA

(1) 初始化打分矩阵

首先将矩阵的第0行与第0列分别用Seq1与Seq2填充，填充时，注意预留出两个字符，令(1，1)= 0，然后每个纵向与横向的单元均比前一个单元加上一个INDEL的罚分

(2) 通过公式构建整个打分矩阵

(3) 回溯

从矩阵最右下角开始，向上回溯

(4) 根据回溯路径得出结果

其中路径朝向左上，即 MATCH/DISMATCH,路 径朝左为 seq1 出现 INDEL 情况,路径朝上为 seq2 出现 INDEL 情况，使用 ’-’ 代 替。最后得出如图结果，为全局最优解

4. Needleman-Wunsch(NW)算法C++实现

4.1 代码

#include <stdio.h>
#include<iostream>
#include <string>
#include<algorithm>
#include <vector>                          
#include <iterator> 
#include<ctime>
#include <windows.h>
//声明命名空间std
using namespace std;
#define MATCH 1
#define DIS_MATCH -1
#define INDEL -3
#define INDEL_CHAR '-'

class ResUnit {
   		//一次双序列比对后的结果
public:
	string str1;	//原始序列1
	string str2;	//原始序列2
	string res1;	//结果序列1
	string res2;	//结果序列2
	int score;		//序列总得分，反映两个序列的相似程度
	int tag;		//禁止迭代多次
};

class SingleSeq {
   	//一个序列被整合后的样子
public:
	string str;		//一个序列的原始序列
	string res;		//一个序列被整合后的样子

};

struct BacktrackingUnit {
   
	int goUp;		//是否向上回溯
	int goLeftUp;	//是否向左上回溯
	int goLeft