Needleman Wunsch 和 Smith Waterman 算法 • CoderXL's Blog

今天学习的两个新算法，Needleman Wunsch 和 Smith Waterman，都使用动态规划思想解决序列匹配问题。

具体来说，给定两个序列

a AACCTATAGCT

b GCGATATA

想要找出两个序列在一定规则下最优的匹配。

其中，“匹配”的定义是分别在 a 和 b 的某些相邻元素之间（或首尾）加入若干空位 -，使之变成长度相等的序列，例如：

a' AACCTATAGCT-

b' --GCGATA--TA

一个匹配的好坏可以使用空位开放罚分、空位扩展罚分、记分矩阵配合使用来评价。例如，如果使用空位开放罚分 $p_o=2$ ，空位扩展罚分 $p_e=2$ ，记分矩阵

\left[\begin{matrix} 3 &-1 &-1 &-1\\ -1 &3 &-1 &-1\\ -1 &-1 &3 &-1\\ -1 &-1 &-1 &3 \end{matrix}\right]

或者简记为 $\displaystyle{f(a_i,b_j)=\left\{\begin{aligned}3&,~~Matched,\\-1&,~~Mismatched.\end{aligned}\right.}$

那么，一个匹配的总分就可以通过对两序列的每个对应位置的分数进行累加得到，进而用于判断其好坏。

而这两个算法要解决的问题就是，找到一个在给定评判标准下，最优的全局匹配和局部匹配（局部匹配找到一对 a 和 b 的子串，这对子串的全局最优匹配分数在所有子串组合中最高）。

先看代码

Needleman-Wunsch.cpp

#include <iostream>
#include <string>
#define N 20
using namespace std;

const int d=2;
struct Dp
{
	int v,ifrm,jfrm;
}dp[N][N];
string sa,sb;
inline int f_cost(char x, char y)
{
	return (x==y?3:-1);
}
int n,m;
int main()
{
    cin>>sa>>sb;
	sa=" "+sa;
	sb=" "+sb;
	m=sa.length();
	n=sb.length();
	for(int i=1;i<m;i++)dp[i][0].v=-i;
	for(int i=1;i<n;i++)dp[0][i].v=-i;
	for(int i=1;i<m;i++)
	{
		for(int j=1;j<n;j++)
		{
			dp[i][j]={dp[i-1][j-1].v+f_cost(sa[i],sb[j]), i-1, j-1};
			if(dp[i-1][j].v-d > dp[i][j].v)
			{
				dp[i][j].v=dp[i-1][j].v-d;
				dp[i][j].ifrm=i-1;
				dp[i][j].jfrm=j;
			}
			if(dp[i][j-1].v-d > dp[i][j].v)
			{
				dp[i][j].v=dp[i-1][j].v-d;
				dp[i][j].ifrm=i;
				dp[i][j].jfrm=j-1;
			}
		}
	}
	for(int i=0;i<m;i++)
	{
		for(int j=0;j<n;j++)
		{
			cout<<dp[i][j].v<<"\t";
		}
		cout<<"\n";
	}
	
	cout<<"\n";
	int curri=m-1,currj=n-1;
	while(curri&&currj)
	{
		cout<<curri<<"\t"<<currj<<"\n";
		int _i=dp[curri][currj].ifrm;
		currj=dp[curri][currj].jfrm;
		curri=_i;
	}
	return 0;
}

cpp

Smith-Waterman.cpp

#include <iostream>
#include <string>
#define N 20
using namespace std;

const int d=2;
struct Dp
{
	int v,ifrm,jfrm;
}dp[N][N];
string sa,sb;
inline int f_cost(char x, char y)
{
	return (x==y?3:-1);
}
int n,m;
int mxdp,mxi,mxj;
int main()
{
    cin>>sa>>sb;
	sa=" "+sa;
	sb=" "+sb;
	m=sa.length();
	n=sb.length();
	for(int i=1;i<m;i++)
	{
		for(int j=1;j<n;j++)
		{
			dp[i][j]={dp[i-1][j-1].v+f_cost(sa[i],sb[j]), i-1, j-1};
			if(dp[i-1][j].v-d > dp[i][j].v)
			{
				dp[i][j].v=dp[i-1][j].v-d;
				dp[i][j].ifrm=i-1;
				dp[i][j].jfrm=j;
			}
			if(dp[i][j-1].v-d > dp[i][j].v)
			{
				dp[i][j].v=dp[i-1][j].v-d;
				dp[i][j].ifrm=i;
				dp[i][j].jfrm=j-1;
			}
			if(dp[i][j].v<0)
			{
				dp[i][j].v=0;
				dp[i][j].ifrm=i;
				dp[i][j].jfrm=j;
			}
			if(dp[i][j].v>mxdp)
			{
				mxdp=dp[i][j].v;
				mxi=i;
				mxj=j;
			}
		}
	}
	for(int i=0;i<m;i++)
	{
		for(int j=0;j<n;j++)
		{
			cout<<dp[i][j].v<<"\t";
		}
		cout<<"\n";
	}
	cout<<"\n";
	cout<<mxi<<"\t"<<mxj<<"\n";
	int curri=mxi,currj=mxj;
	while(dp[curri][currj].v&&curri&&currj)
	{
		cout<<dp[curri][currj].ifrm<<"\t"<<dp[curri][currj].jfrm<<"\n";
		int _i=dp[curri][currj].ifrm;
		currj=dp[curri][currj].jfrm;
		curri=_i;
	}
	return 0;
}

cpp