双序列比对之重复匹配算法以及C++实现

最新推荐文章于 2023-07-04 09:03:01 发布

隔壁王同学啊

最新推荐文章于 2023-07-04 09:03:01 发布

阅读量2k

点赞数

分类专栏：计算生物学

本文链接：https://blog.csdn.net/weixin_43770577/article/details/105396904

版权

计算生物学专栏收录该内容

23 篇文章 6 订阅

订阅专栏

双序列比对之重复匹配算法

在看此篇文章之前，建议先看一看我之前写的NW算法和SW算法（近期准备将这两篇文章在重新写一下，主要添加一些内容。）
重复匹配算法是迥异于之前介绍的全局比对和局部比对算法。如果我们将要比对的两个序列命名为序列1，它为AWGHEEBAWGRBZ,序列2为PEAWG。

全局比对之后的结果可能为
HE—AWGEBAWGRBZAWG
—PEAWG-------------------

局部比对之后的可能结果为：
H–EAWGEBAWGRBZAWG
PE–AWG-------------------
但是当待比对序列中存在多个重复的片段时候，如我们想找出蛋白质中模体或域的多个拷贝时，上述的算法，显然并不能满足需求。如序列1中AWG重复的出现，序列2中也有AWG的出现，我们需求是只要序列一中出现了包含了AWG的连续片段，我们就让序列2进行联配并显示。
通俗来讲序列1就是母串，序列2就是子串。联配时候字串可以重复的使用，然后将所有匹配上的片段进行打印出来。如：
HEAWGEBAWGRBZAWG
.EAWG.E…AWG…AWG
“.”是间断点，是字串与子串之间的分界点。
核心思想就是一个匹配片段完成后，字串还可以对该片段之后的母串继续进行匹配。这也是它与全局序列和局部序列比对的最大不同之处。
那么如何实现该算法呢。显然我们可以通过将子串与母串先进行对齐，进行一次比对，然后在将子串往后一个位置，再次进行比对，直到推=退到母串的末尾。这是最容易想的算法思路，但同时它的时间复杂度也是很大。想到至此，我们又想到了用KMP算法，对子串进行预处理，得到一个存储前缀长度的next数组，然后再进行匹配，可以降低时间复杂度。但是啊，这个比对太“精准”了，我们要知道整个序列对比的基础就是概率论基础，通过对先验生物学数据进行整合，得到了一个打分矩阵，我们认为当一个比对结果的得分大于某一个阈值时，就代表这结果是可接受。比如字串中的ACG与母串中ACG进行匹配的得分为30，而ACG与-CG得分为25，如果我们将算法中的阈值设置为20，那么表明得分25的匹配结果也是可行的，因为我们要考虑生物进化中的基因突变，即碱基对的插入和缺失等等，-CG表明在进化时发生A缺失了。而KMP算法就是如果字串中是ACG,那么只有母串中出现ACG才算匹配成功，所以说它返回的匹配结果太“精准了”，只有当你特定的想搜索一个特别保守的序列，才考虑使用。或者具有特定的功能的序列，如起始密码子，TATA盒等。
ok,上面的算法思想，你都可以写着试试。其实今天的重复匹配算法还是基于动态规划算法,它和Neddle-Wunch和Smith-Waterman算法是一脉相承的，只不过当形成动态规划表时，重复匹配算法是一列列遍历形成的，而其他两者是一行行遍历形成。因为你想，我如果把子串放列的方向上，母串放在行的方向上，那么我们一列列的遍历是不是就相当于使用子串进行一次又一次匹配，当这次的匹配得分大于阈值时（我们事先自己设置的）我们设置一个间断点，并且这个间断点在第0行，如果有几个间断点，就说明存在几个得分大于阈值的匹配。（如果动态规划不清楚的可以看看我之前关于Neddle-Wunch和Smith-Waterman算法的博客）。
最后我用c++将这个算法复现了，基本思路还是和之前复现Neddle-Wunch和Smith-Waterman算法的思路一致。如果不想看具体的代码的话，我下面简要的叙述一下伪代码思路。结构体数组dp的v值记录了当前匹配的累积最大分值。w值记录了这个最大分值是如何产生的，以dp[2][3]为例（s1字符串是存储了母串，s2存储了子串），如果它的w值是0就代表是s2[2]匹配了一个空格（gap），当动态规划表建立完成时，开始打印匹配情况进行回溯打印，当回溯至该元素时候，就打印出s2[2]和空格匹配，并且下一个回溯至dp[1][3]。如果w值是1，就表示s2[2]与s1[3]进行匹配，回溯打印时，就打印出s2[2]与s1[3]匹配，并且下一个回溯至dp[1][2]。w值是2是s1[3]和空格匹配，并且下一个回溯至dp[2][2]。而动态规划表形成的逻辑见下面的式子。
$\max \left\{ \begin{array} { l } F ( 0,j-1 ) \\ F ( i , j-1 ) - T , i = 1 , \cdots , m \ \end{array} \right.$
(该逻辑是由下面代码中find_max_1函数实现)

$\max \left\{ \begin{array} { l } F ( 0 , j ) \\ F ( i - 1 , j - 1 ) + s \left( x _ { i } , y _ { j } \right) \\ F ( i - 1 , j ) - d \\ F ( i , j - 1 ) - d \end{array} \right.$
(该逻辑是由下面代码中find_max_2函数实现)

其中T代表设置的阈值，d代表空位罚分，s(，)是打分矩阵。并且事先将动态规划表的第零列的元素初始化为0。
下面的代码是重复匹配算法对两条核苷酸序列进行匹配。

#include<iostream>
#include<string>
using namespace std;
string s1,s2;
int  s[5][5];
int  gap=-5;
const int T=5;//可以根据序列的特点进行阈值的设置，这我们暂且将其设置为5。
struct Pair{
int w;//w的允许有三个值，0代表向上的箭头，1代表斜向下的箭头，2代表向右的箭头。
int v;//v是该位置的累计积分。
bool bre;//如果bre为真，那么说明此点是一个间断点，那么此时col值存储间断点指向上一个点的行号，vol存储列号。
int row;
int col; 
}dp[101][101]; 
void get_score(){
  //生成核苷酸打分矩阵
  s[0][0]=10;s[1][0]=-1;s[2][0]=-3;s[2][1]=-5;
	s[1][1]=7;s[3][0]=-4;s[3][1]=-3;s[3][2]=0;
	s[2][2]=9;
	s[3][3]=8;
	for(int i=0;i<4;i++){
		for(int j=0;j<4;j++){
			if(s[i][j]!=0){
				s[j][i]=s[i][j];
			}
		}
	}
}
int tran(char a){
  //将输入的字符串AGCT转化为0，1，2，3 
 	int c;
 	if (a=='A'){
 		c=0;
	}
 	if (a=='G'){
 		c=1;
	}if (a=='C'){
 		c=2;
	}if (a=='T'){
 		c=3;
	}
 	return c;
 }
void find_max_1(Pair dp[101][101],int s,int s2_length,int i){
	int tem_row;
	int tem_col;
	int tem_max;
	for(int k=0;k<=s2_length;k++){
	if(tem_max<=dp[k][s].v-T){
		tem_max=dp[k][s].v-T;
		 tem_row=k;
		 tem_col=s;
		}
	}
	if(tem_max>=dp[i][s].v){
		dp[i][s+1].bre=true;
		dp[i][s+1].row=tem_row;
		dp[i][s+1].col=tem_col;
		dp[i][s+1].v=tem_max;
	}
	else{
		dp[i][s+1].bre=false;
		dp[i][s+1].v=dp[i][s].v;
		dp[i][s+1].w=2;
	}

}
void find_max_2(Pair dp[101][101],int i,int j){
	int cnt=-10000;
	int left=dp[i][j-1].v+gap;
	if(cnt<=left){
		cnt=left;
		dp[i][j].w=2;
		dp[i][j].v=cnt;

	}
	int up=dp[i-1][j].v+gap;
	if(cnt<=up){
		cnt=up;
		dp[i][j].w=0;
		dp[i][j].v=cnt;
	}
	int i_i=tran(s1[i-1]);
	int j_j=tran(s2[j-1]);
	int inc=dp[i-1][j-1].v+s[i_i][j_j];
	if(cnt<=inc){
		cnt=inc;
		dp[i][j].w=1;
		dp[i][j].v=cnt;
	}
}
int main(){
	get_score();
  	cout << "this is seuence_a:";
	cin >> s1 ;//因为开的数组的长度为101，所以输入的字符串应该小于100。
  	cout<< "this is sequence_b:"; 
  	cin>> s2;
	int s1_length=s1.size();
	int s2_length=s2.size();
	
	dp[0][0].v=0;
	dp[0][0].w=0;
	for(int j=0;j<=s2_length;j++){
  		dp[j][0].v=0;
  		dp[j][0].w=0;
  		dp[j][0].bre=false;
	}//对第一列进行初始化。
	for(int j=1;j<=s1_length;j++){
	//按一列一列遍历生成dp表。
  		for(int i=0;i<=s2_length;i++){
	  		if(i==0){
		  	//处理第0行的元素，判断它是否有可能是间断点，如果是间断点的话，就将bre赋值为true。
				find_max_1(dp,j-1,s2_length,i);
	    	}
	  	else{
		  //处理不是第0行的元素。
		  	find_max_2(dp,i,j);
	  		}
  		}
	}
	//以下是将匹配结果打印出来的程序。
	//不同于nw算法从dp表的右下角，开始回溯打印比对的情况，重复匹配算法是从dp表的第零行最后一个格子开始回溯。
	//有几个break点，就有几个比对得分大于我们设置的阈值T。
	string out[101];
	int j=s1_length;
	while(j>0){
		if(dp[0][j].bre==false){
			out[j]="*";//星号表示间断点。
			j-=1;
		}
		else{
			if(out[j]!="*")
			out[j]="*";//开始打印匹配的片段。从下一个点开始跳跃。
			int i=dp[0][j].row;
			j=dp[0][j].col;
			while(i>0){
				out[j]=s2[i-1];
				if(dp[i][j].w==0){
					i-=1;
					out[j]="^";
				}
				if(dp[i][j].w==1){
					out[j]=s2[i-1];
					i-=1;
					j-=1;
				}
				if(dp[i][j].w==2){
					out[j]="-";
					j-=1;
				}
			}
		}
	}
	cout << "this is the result:";
	cout<< endl;
	cout << s1;
	cout << endl;
	for(int i=1;i<=s1_length;i++){
		cout << out[i];
	}
	cout << endl;
	cout<< "this is score matrix:";
	for(int i=0;i<=s2_length;i++){
		cout << endl;
		for(int j=0;j<=s1_length;j++){
			cout << dp[i][j].v<< " ";
		}
	}
	return 0;
}

下面是运行该算法得到的结果截图。

QQ图片.png
下面的代码是重复匹配算法对两条氨基酸序列进行比对，里面是内嵌了BLOSUM62矩阵。其他和上面那个代码想差不多。

#include<iostream>
#include<string>
#include<fstream>
#include<map>
using namespace std;
string s1,s2;
string ss="ARNDCQEGHLKMFPSTWYVBZX";
map <char,int> amino;
int  gap=-4;
int s[25][25];
const int T=2;//可以根据序列的特点进行阈值的设置，这我们暂且将其设置为2。
struct Pair{
int w;//w的允许有三个值，0代表向上的箭头，1代表斜向下的箭头，2代表向右的箭头。
int v;//v是该位置的累计积分。
bool bre;//如果bre为真，那么说明此点是一个间断点，那么此时col值存储间断点指向上一个点的行号，vol存储列号。
int row;
int col; 
}dp[101][101]; 
void get_score(){
  //读入blosum_62矩阵
  ifstream file_1("blosum62.txt");
    for(int i=0;i<23;i++)
    {
        for(int j=0;j<23;j++)
        {
            file_1>>s[i][j];
        }
    }
  file_1.close();
}
void map_char_int(){
    for(int i=0;i<23;i++){
        amino[ss[i]]=i;//将A映射为0，将R映射1.......。
    }
 }
void find_max_1(Pair dp[101][101],int s,int s2_length,int i){
	int tem_row;
	int tem_col;
	int tem_max;
	for(int k=0;k<=s2_length;k++){
	if(tem_max<=dp[k][s].v-T){
		tem_max=dp[k][s].v-T;
		 tem_row=k;
		 tem_col=s;
		}
	}
	if(tem_max>=dp[i][s].v){
		dp[i][s+1].bre=true;
		dp[i][s+1].row=tem_row;
		dp[i][s+1].col=tem_col;
		dp[i][s+1].v=tem_max;
	}
	else{
		dp[i][s+1].bre=false;
		dp[i][s+1].v=dp[i][s].v;
		dp[i][s+1].w=2;
	}

}
void find_max_2(Pair dp[101][101],int i,int j){
	int cnt=-10000;
	int left=dp[i][j-1].v+gap;
	if(cnt<=left){
		cnt=left;
		dp[i][j].w=2;
		dp[i][j].v=cnt;

	}
	int up=dp[i-1][j].v+gap;
	if(cnt<=up){
		cnt=up;
		dp[i][j].w=0;
		dp[i][j].v=cnt;
	}
	int inc=dp[i-1][j-1].v+s[amino[s2[j-1]]][amino[s1[i-1]]];
	if(cnt<=inc){
		cnt=inc;
		dp[i][j].w=1;
		dp[i][j].v=cnt;
	}
}
int main(){
	get_score();
    map_char_int();
  	cout << "this is sequence_a:";
	cin >> s1 ;//因为开的数组的长度为101，所以输入的字符串应该小于100。
  	cout<< "this is sequence_b:"; 
  	cin>> s2;
	int s1_length=s1.size();
	int s2_length=s2.size();
	
	dp[0][0].v=0;
	dp[0][0].w=0;
	for(int j=0;j<=s2_length;j++){
  		dp[j][0].v=0;
  		dp[j][0].w=0;
  		dp[j][0].bre=false;
	}//对第一列进行初始化。
	for(int j=1;j<=s1_length;j++){
	//按一列一列遍历生成dp表。
  		for(int i=0;i<=s2_length;i++){
	  		if(i==0){
		  	//处理第0行的元素，判断它是否有可能是间断点，如果是间断点的话，就将bre赋值为true。
				find_max_1(dp,j-1,s2_length,i);
	    	}
	  	else{
		  //处理不是第0行的元素。
		  	find_max_2(dp,i,j);
	  		}
  		}
	}
	//以下是将匹配结果打印出来的程序。
	//不同于nw算法从dp表的右下角，开始回溯打印比对的情况，重复匹配算法是从dp表的第零行最后一个格子开始回溯。
	//有几个break点，就有几个比对得分大于我们设置的阈值T。
	string out[101];
	int j=s1_length;
	while(j>0){
		if(dp[0][j].bre==false){
			out[j]="*";//星号表示间断点。
			j-=1;
		}
		else{
			if(out[j]!="*")
			out[j]="*";//开始打印匹配的片段。从下一个点开始跳跃。
			int i=dp[0][j].row;
			j=dp[0][j].col;
			while(i>0){
				out[j]=s2[i-1];
				if(dp[i][j].w==0){
					i-=1;
					out[j]="^";
				}
				if(dp[i][j].w==1){
					out[j]=s2[i-1];
					i-=1;
					j-=1;
				}
				if(dp[i][j].w==2){
					out[j]="-";
					j-=1;
				}
			}
		}
	}
	cout << "this is the result:";
	cout<< endl;
	cout << s1;
	cout << endl;
	for(int i=1;i<=s1_length;i++){
		cout << out[i];
	}
	cout << endl;
	cout<< "this is score matrix:";
	for(int i=0;i<=s2_length;i++){
		cout << endl;
		for(int j=0;j<=s1_length;j++){
			cout << dp[i][j].v<< " ";
		}
	}
	return 0;
}

上面就是本片博文全部内容了，如果有什么不懂的话可以详细看看《生物序列分析》。

隔壁王同学啊

关注

0
点赞
踩
9

收藏

觉得还不错? 一键收藏
2
评论
双序列比对之重复匹配算法以及C++实现

双序列比对之重复匹配算法在看此篇文章之前，建议先看一看我之前写的NW算法和SW算法（近期准备将这两篇文章在重新写一下，主要添加一些内容。）重复匹配算法是迥异于之前介绍的全局比对和局部比对算法。如果我们将要比对的两个序列命名为序列1，它为AWGHEEBAWGRBZ,序列2为PEAWG。全局比对之后的结果可能为 HE—AWGEBAWGRBZAWG —PEAWG-----------...
复制链接

扫一扫

专栏目录