哇。生物专业的复习月可不是盖的,十几门的专业课,并且是纯背的。。。。提前一个月开始复习,复习吐了,写个博客调节调节。。
本次介绍的局部双序列对比算法,也就是SW算法。基本上的话,只要学了NW算法也就差不多理解了SW算法了。这个算法顾名思义,就是找到最佳的局部片段。注意这个最佳,怎末实现。
这个最佳可以翻译为,从所有的比对情况中,截取得分最高的片段。当一个比对情况形成时,也就是这个情况的计分值不断累加时,当进行至本位点比对时,累计的分值变为了负值,那么将这个负值直接赋值为0,然后进行下一个位点的比对。
最后dp表形成后,遍历找到最大值的格子,从此处进行回溯,直到回溯到0值。
#include<iostream>
#include<string>
using namespace std;
struct Pair{
int w;//w的允许有三个值,0代表向下的箭头,1代表斜向下的箭头,2代表向右的箭头。
int v;//v是该位置的累计积分。
};
int tran(char a){
int c;
if (a=='A'){
c=0;
}
if (a=='G'){
c=1;
}
if (a=='C'){
c=2;
}
if (a=='T'){
c=3;
}
return c;
}
int main(){
string s1;
string s2;
cin >> s1 >> s2;//输入要比对的核苷酸的序列,双序列
int m,n;
m=s1.size();
n=s2.size();
//将输入的字符串AGCT转化为0,1,2,3
int s[5][5];
int gap=-5;
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];
}
}
}//输入计分矩阵
Pair dp[1000][1000];
dp[0][0].v=0;
dp[0][0].w=0;
for(int i=0;i<=m;i++){
for(int j=0;j<=n;j++){
if(i==0&&j!=0){
dp[i][j].v=dp[i][j-1].v+gap;
dp[i][j].w=2;
}
if(j==0&&i!=0){
dp[i][j].v=dp[i-1][j].v+gap;
dp[i][j].w=0;
}//以上将第一行第一列分别进行初始化。
}
}
for(int i=1;i<=m;i++){
for(int j=1;j<=n;j++){
int x=tran(s1[i-1]);
int y=tran(s2[j-1]);
if(dp[i-1][j].v+gap>dp[i-1][j-1].v+s[x][y]){
dp[i][j].v=dp[i-1][j].v+gap;
dp[i][j].w=0;
}
else{
dp[i][j].v=dp[i-1][j-1].v+s[x][y];
dp[i][j].w=1;
}
if(dp[i][j].v<dp[i][j-1].v+gap){
dp[i][j].w=2;
dp[i][j].v=dp[i][j-1].v+gap;
}
if(dp[i][j].v<0){
dp[i][j].v=0;
}//当此格子的累计分值小于0时,将其赋值为0。达到局部匹配的效果
}
}
cout << "最终得分矩阵为" << endl;
for(int i=0;i<=m;i++){
for(int j=0;j<=n;j++){
cout << dp[i][j].v << " ";
if(j==n){
cout << endl;
}
}
}
cout << "输出标记的路径:" << endl;
for(int i=0;i<=m;i++){
for(int j=0;j<=n;j++){
cout << dp[i][j].w<< " ";
if(j==n){
cout << endl;
}
}
}
string s3,s4;//存放排队的方案
int temp=0;
int i_max;
int j_max;
for(int i=0;i<=m;i++){
for(int j=0;j<=n;j++){
if(dp[i][j].v>temp){
i_max=i;
j_max=j;
temp=dp[i][j].v;
}
}
}//遍历找出最大累计分值
int xx=0;
int i=i_max;
int j=j_max;
while(dp[i][j].v!=0){
if(dp[i][j].w==1){
s3[xx]=s1[i-1];
s4[xx]=s2[j-1];
i-=1;
j-=1;
}
if(dp[i][j].w==0){
s3[xx]=s1[i-1];
s4[xx]=' ';
i-=1;
}
if(dp[i][j].w==2){
s3[xx]=' ';
s4[xx]=s2[j-1];
j-=1;
}
xx+=1;
}//以下输出局部序列比对的结果
for(int i=xx-1;i>=0;i--){
cout << s3[i];
}
cout << endl;
for(int i=xx-1;i>=0;i--){
cout << s4[i];
}
return 0;
}