生物信息学(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算法罚分规则如下:
(2-1)

以罚分规则为基础,通过以下公式构造一个打分矩阵:
在这里插入图片描述
最后通过回溯得分矩阵,获取全局最优解,回溯的规则是单元从哪里来,就回溯到哪里去。这样说很抽象,接下来举一个例子。

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;		//是否向左回溯
	int score;		//得分矩阵第(i, j)这个单元的分值
};


typedef struct BacktrackingUnit *unitLine;

int max3(int a, int b, int c);
ResUnit traceback(unitLine** item, const int i, const int j, string str1, string str2, string res1, string res2, int n, ResUnit resUnit);
ResUnit NeedlemanWunch(string str1, string str2);

struct SequenceUnit {
	string *str1;	//匹配序列1
	string *str2;	//匹配序列2
	int score;
};

int main() {
	ResUnit s = NeedlemanWunch("GGATCGA","GAATTCAGTTA");
	cout << s.res1 << endl;
	cout << s.res2 << endl;
	system("pause");
}


/**
比较三种路径之间谁最大

f(i-1,j-1),f(i-1,j)+indel,f(i,j-1)+indel
*/
int max3(int a, int b, int c) {
	int temp = a > b ? a : b;
	return temp > c ? temp : c;
}

/**
比较两个字符类型属于什么,match,dismatch,indel
*/
int myCompare(char a, char b) {
	if (a == b)
		return MATCH;
	else if (a == ' ' || b == ' ')
		return INDEL;
	else
		return DIS_MATCH;
}


ResUnit traceback(unitLine** item, const int i, const int j, string str1, string str2, string res1, string res2, int n, ResUnit resUnit) {
	unitLine temp = item[i][j];
	if (resUnit.tag != 1)
	{
		if (!(i || j)) {   // 到矩阵单元(0, 0)才算结束,这代表初始的两个字符串的每个字符都被比对到了

			resUnit.str1 = str1;
			resUnit.str2 = str2;
			resUnit.res1 = res1;
			resUnit.res2 = res2;
			resUnit.tag = 1;
			return resUnit;
		}
		if (temp->goUp) {    // 向上回溯一格
			res1 = str1[i - 1] + res1;
			res2 = INDEL_CHAR + res2;
			resUnit = traceback(item, i - 1, j, str1, str2, res1, res2, n + 1, resUnit);
		}
		if (temp->goLeftUp) {    // 向左上回溯一格 
			res1 = str1[i - 1] + res1;
			res2 = str2[j - 1] + res2;
			resUnit = traceback(item, i - 1, j - 1, str1, str2, res1, res2, n + 1, resUnit);
		}
		if (temp->goLeft) {    // 向左回溯一格
			res1 = INDEL_CHAR + res1;
			res2 = str2[j - 1] + res2;
			resUnit = traceback(item, i, j - 1, str1, str2, res1, res2, n + 1, resUnit);
		}
		return resUnit;
	}
	else
	{
		return resUnit;
	}

}


ResUnit NeedlemanWunch(string str1, string str2) {
	//字符串str1,str2长度
	const int m = str1.length();
	const int n = str2.length();

	int m1, m2, m3, mm;

	unitLine **unit;

	// 初始化
	if ((unit = (unitLine **)malloc(sizeof(unitLine*) * (m + 1))) == NULL) {
		fputs("Error: Out of space!\n", stderr);
		exit(1);
	}
	for (int i = 0; i <= m; i++) {
		if ((unit[i] = (unitLine *)malloc(sizeof(unitLine) * (n + 1))) == NULL) {
			fputs("Error: Out of space!\n", stderr);
			exit(1);
		}
		for (int j = 0; j <= n; j++) {
			if ((unit[i][j] = (unitLine)malloc(sizeof(struct BacktrackingUnit))) == NULL) {
				fputs("Error: Out of space!\n", stderr);
				exit(1);
			}
			unit[i][j]->goUp = 0;
			unit[i][j]->goLeftUp = 0;
			unit[i][j]->goLeft = 0;
		}
	}
	unit[0][0]->score = 0;
	for (int i = 1; i <= m; i++) {
		unit[i][0]->score = INDEL * i;
		unit[i][0]->goUp = 1;
	}
	for (int j = 1; j <= n; j++) {
		unit[0][j]->score = INDEL * j;
		unit[0][j]->goLeft = 1;
	}


	// 动态规划算法计算得分矩阵每个单元的分值
	for (int i = 1; i <= m; i++) {
		for (int j = 1; j <= n; j++) {
			m1 = unit[i - 1][j]->score + INDEL;
			m2 = unit[i - 1][j - 1]->score + myCompare(str1[i - 1], str2[j - 1]);
			m3 = unit[i][j - 1]->score + INDEL;
			mm = max3(m1, m2, m3);
			unit[i][j]->score = mm;
			//判断路径来源
			if (m1 == mm) unit[i][j]->goUp = 1;
			if (m2 == mm) unit[i][j]->goLeftUp = 1;
			if (m3 == mm) unit[i][j]->goLeft = 1;
		}
	}


	//开始回溯
	ResUnit res;
	res.tag = 0;
	res = traceback(unit, m, n, str1, str2, "", "", 0, res);
	res.score = unit[m][n]->score;


	//释放内存
	for (int i = 0; i <= m; i++) {
		for (int j = 0; j <= n; j++) {
			free(unit[i][j]);
		}
		free(unit[i]);
	}
	free(unit);

	//返回值
	return res;
}
4.2 运行结果

在这里插入图片描述

总结

NW算法是双序列比对最基础的算法,且是一个全局比对算法,可以看到,他的回溯路径贯穿了整个矩阵。

同时肉眼可见的一点就是回溯路径上的数字几乎都是正数,而且在打分矩阵除了回溯路径之外的地方有些很大的负数,且完全没有用到。

当然还有,就是回溯路径虽然弯弯曲曲,但是大体上是在斜对称轴上的。

对于回溯路径之外,负数太多的问题,Smith-Waterman给出了解决方案。下一篇即介绍其升级版本Smith-Waterman(SW)算法。而回溯路径在对称轴上这一个规律,在BLAST算法中得到利用。

我的毕设题目是:基于CUDA的多序列比对并行算法的设计与实现,感谢导师赏脸让我成功毕业,所以接下来我会用几篇文章逐步分享我的毕设成果,成果~在这里插入图片描述

包括其他序列比对算法,如SW、BLAST、CLUSTAL,以及如何将CLUSTAL用GPU加速,感兴趣的同学可以移步至同专栏的其他博客。
如果有什么问题,请私信联系我或者在评论区留言
码字不易,若有帮助,给个关注和赞呗

在这里插入图片描述

  • 40
    点赞
  • 136
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 18
    评论
全局比对是一种常用的序列比对方法,它可以比对两个序列的整个长度,并给出它们的相似程度。其中 Needleman-Wunsch 是一种经典的全局比对算法Needleman-Wunsch 算法的本质是动态规划,它将全局比对的问题拆分成若干个子问题,并通过计算每个子问题的得分来求解最终的比对结果。具体实现步骤如下: 1. 初始化得分矩阵(score matrix)并计算第一行和第一列的得分; 2. 逐行逐列计算得分矩阵中的每个元素的得分,根据得分矩阵中每个元素的左上、上、左三个元素的得分值,计算当前元素的得分; 3. 回溯得分矩阵,得到最优比对结果。 以下是 Needleman-Wunsch 算法的 Python 实现: ```python # Needleman-Wunsch 算法实现 import numpy as np def needleman_wunsch(seq1, seq2, match_score=1, mismatch_score=-1, gap_score=-1): # 初始化得分矩阵 m, n = len(seq1), len(seq2) score_matrix = np.zeros((m+1, n+1)) for i in range(1, m+1): score_matrix[i,0] = score_matrix[i-1,0] + gap_score for j in range(1, n+1): score_matrix[0,j] = score_matrix[0,j-1] + gap_score # 填充得分矩阵 for i in range(1, m+1): for j in range(1, n+1): score1 = score_matrix[i-1,j-1] + (match_score if seq1[i-1] == seq2[j-1] else mismatch_score) score2 = score_matrix[i-1,j] + gap_score score3 = score_matrix[i,j-1] + gap_score score_matrix[i,j] = max(score1, score2, score3) # 回溯得分矩阵,得到最优比对结果 align1, align2 = '', '' i, j = m, n while i > 0 and j > 0: score, score1, score2, score3 = score_matrix[i,j], score_matrix[i-1,j-1], score_matrix[i-1,j], score_matrix[i,j-1] if score == score1 + (match_score if seq1[i-1] == seq2[j-1] else mismatch_score): align1 = seq1[i-1] + align1 align2 = seq2[j-1] + align2 i, j = i-1, j-1 elif score == score2 + gap_score: align1 = seq1[i-1] + align1 align2 = '-' + align2 i -= 1 else: align1 = '-' + align1 align2 = seq2[j-1] + align2 j -= 1 while i > 0: align1 = seq1[i-1] + align1 align2 = '-' + align2 i -= 1 while j > 0: align1 = '-' + align1 align2 = seq2[j-1] + align2 j -= 1 return align1, align2 ``` 在上面的代码中,我们使用了 `numpy` 库来创建得分矩阵,这使得代码更加简洁和高效。在实现中,我们还可以通过调整 `match_score`、`mismatch_score` 和 `gap_score` 的值来控制比对的结果。
评论 18
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

大青儿

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值