Needleman/Wunsch算法介绍以及JAVA代码编写
首先声明:网上关于Needleman/Wunsch算法有很很多,也有更详细的介绍,但是关于用JAVA如何去实现两段文本字符串的比较就比较少,本文是结合大量同行的的代码的同时加以总结,希望能对后来阅读的人有所帮助,若对某哥的代码有所抄袭,请见谅:
part 1:Needleman/Wunsch算法介绍
本文介绍基于最长公共子串的文本比较算法——Needleman/Wunsch算法。
还是以实例说明:字符串A=kitten,字符串B=sitting
那他们的最长公共子串为ittn(注:最长公共子串不需要连续出现,但一定是出现的顺序一致),最长公共子串长度为4。
定义:
LCS(A,B)表示字符串A和字符串B的最长公共子串的长度。很显然,LSC(A,B)=0表示两个字符串没有公共部分。
Rev(A)表示反转字符串A
Len(A)表示字符串A的长度
A+B表示连接字符串A和字符串B
性质:
LCS(A,A)=Len(A)
LCS(A,"")=0
LCS(A,B)=LCS(B,A)
0≤LCS(A,B)≤Min(Len(A),Len(B))
LCS(A,B)=LCS(Rev(A),Rev(B))
LCS(A+C,B+C)=LCS(A,B)+Len(C)
LCS(A+B,A+C)=Len(A)+LCS(B,C)
LCS(A,B)≥LCS(A,C)+LCS(B,C)
LCS(A+C,B)≥LCS(A,B)+LCS(B,C)
为了讲解计算LCS(A,B),特给予以下几个定义
A=a1a2……aN,表示A是由a1a2……aN这N个字符组成,Len(A)=N
B=b1b2……bM,表示B是由b1b2……bM这M个字符组成,Len(B)=M
定义LCS(i,j)=LCS(a1a2……ai,b1b2……bj),其中0≤i≤N,0≤j≤M
故: LCS(N,M)=LCS(A,B)
LCS(0,0)=0
LCS(0,j)=0
LCS(i,0)=0
对于1≤i≤N,1≤j≤M,有公式一
若ai=bj,则LCS(i,j)=LCS(i-1,j-1)+1
若ai≠bj,则LCS(i,j)=Max(LCS(i-1,j-1),LCS(i-1,j),LCS(i,j-1))
计算LCS(A,B)的算法有很多,下面介绍的Needleman/Wunsch算法是其中的一种。和LD算法类似,Needleman/Wunsch算法用的都是动态规划的思想。在Needleman/Wunsch算法中还设定了一个权值,用以区分三种操作(插入、删除、更改)的优先级。在下面的算法中,认为三种操作的优先级都一样。故权值默认为1。
举例说明:A=GGATCGA,B=GAATTCAGTTA,计算LCS(A,B)
第一步:初始化LCS矩阵
G | A | A | T | T | C | A | G | T | T | A | ||
---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | |
G | 0 | |||||||||||
G | 0 | |||||||||||
A | 0 | |||||||||||
T | 0 | |||||||||||
C | 0 | |||||||||||
G | 0 | |||||||||||
A | 0 |
第二步:利用公式一,计算矩阵的第一行
G | A | A | T | T | C | A | G | T | T | A | ||
---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | |
G | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 |
G | 0 | |||||||||||
A | 0 | |||||||||||
T | 0 | |||||||||||
C | 0 | |||||||||||
G | 0 | |||||||||||
A | 0 |
第三步:利用公式一,计算矩阵的其余各行
G | A | A | T | T | C | A | G | T | T | A | ||
---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | |
G | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 |
G | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 2 | 2 | 2 | 2 |
A | 0 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
T | 0 | 1 | 2 | 2 | 3 | 3 | 3 | 3 | 3 | 3 | 3 | 3 |
C | 0 | 1 | 2 | 2 | 3 | 3 | 4 | 4 | 4 | 4 | 4 | 4 |
G | 0 | 1 | 2 | 2 | 3 | 3 | 3 | 4 | 5 | 5 | 5 | 5 |
A | 0 | 1 | 2 | 3 | 3 | 3 | 3 | 4 | 5 | 5 | 5 | 6 |
则,LCS(A,B)=LCS(7,11)=6
可以看出,Needleman/Wunsch算法实际上和LD算法是非常接近的。故他们的时间复杂度和空间复杂度也一样。时间复杂度为O(MN),空间复杂度为O(MN)。空间复杂度经过优化,可以优化到O(M),但是一旦优化就丧失了计算匹配字串的机会了。由于代码和LD算法几乎一样。这里就不再贴代码了。
还是以上面为例A=GGATCGA,B=GAATTCAGTTA,LCS(A,B)=6
他们的匹配为:
A:GGA_TC_G__A
B:GAATTCAGTTA
如上面所示,蓝色表示完全匹配,黑色表示编辑操作,_表示插入字符或者是删除字符操作。如上面所示,蓝色字符有6个,表示最长公共子串长度为6。
利用上面的Needleman/Wunsch算法矩阵,通过回溯,能找到匹配字串
第一步:定位在矩阵的右下角
G | A | A | T | T | C | A | G | T | T | A | ||
---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | |
G | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 |
G | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 2 | 2 | 2 | 2 |
A | 0 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
T | 0 | 1 | 2 | 2 | 3 | 3 | 3 | 3 | 3 | 3 | 3 | 3 |
C | 0 | 1 | 2 | 2 | 3 | 3 | 4 | 4 | 4 | 4 | 4 | 4 |
G | 0 | 1 | 2 | 2 | 3 | 3 | 3 | 4 | 5 | 5 | 5 | 5 |
A | 0 | 1 | 2 | 3 | 3 | 3 | 3 | 4 | 5 | 5 | 5 | 6 |
第二步:回溯单元格,至矩阵的左上角
若ai=bj,则回溯到左上角单元格
G | A | A | T | T | C | A | G | T | T | A | ||
---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | |
G | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 |
G | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 2 | 2 | 2 | 2 |
A | 0 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
T | 0 | 1 | 2 | 2 | 3 | 3 | 3 | 3 | 3 | 3 | 3 | 3 |
C | 0 | 1 | 2 | 2 | 3 | 3 | 4 | 4 | 4 | 4 | 4 | 4 |
G | 0 | 1 | 2 | 2 | 3 | 3 | 3 | 4 | 5 | 5 | 5 | 5 |
A | 0 | 1 | 2 | 3 | 3 | 3 | 3 | 4 | 5 | 5 | 5 | 6 |
若ai≠bj,回溯到左上角、上边、左边中值最大的单元格,若有相同最大值的单元格,优先级按照左上角、上边、左边的顺序
G | A | A | T | T | C | A | G | T | T | A | ||
---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | |
G | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 |
G | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 2 | 2 | 2 | 2 |
A | 0 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
T | 0 | 1 | 2 | 2 | 3 | 3 | 3 | 3 | 3 | 3 | 3 | 3 |
C | 0 | 1 | 2 | 2 | 3 | 3 | 4 | 4 | 4 | 4 | 4 | 4 |
G | 0 | 1 | 2 | 2 | 3 | 3 | 3 | 4 | 5 | 5 | 5 | 5 |
A | 0 | 1 | 2 | 3 | 3 | 3 | 3 | 4 | 5 | 5 | 5 | 6 |
若当前单元格是在矩阵的第一行,则回溯至左边的单元格
若当前单元格是在矩阵的第一列,则回溯至上边的单元格
G | A | A | T | T | C | A | G | T | T | A | ||
---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | |
G | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 |
G | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 2 | 2 | 2 | 2 |
A | 0 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
T | 0 | 1 | 2 | 2 | 3 | 3 | 3 | 3 | 3 | 3 | 3 | 3 |
C | 0 | 1 | 2 | 2 | 3 | 3 | 4 | 4 | 4 | 4 | 4 | 4 |
G | 0 | 1 | 2 | 2 | 3 | 3 | 3 | 4 | 5 | 5 | 5 | 5 |
A | 0 | 1 | 2 | 3 | 3 | 3 | 3 | 4 | 5 | 5 | 5 | 6 |
依照上面的回溯法则,回溯到矩阵的左上角
第三步:根据回溯路径,写出匹配字串
若回溯到左上角单元格,将ai添加到匹配字串A,将bj添加到匹配字串B
若回溯到上边单元格,将ai添加到匹配字串A,将_添加到匹配字串B
若回溯到左边单元格,将_添加到匹配字串A,将bj添加到匹配字串B
搜索晚整个匹配路径,匹配字串也就完成了
可以看出,LD算法和Needleman/Wunsch算法的回溯路径是一样的。这样找到的匹配字串也是一样的。
不过,Needleman/Wunsch算法和LD算法一样,若要找出匹配字串,空间的复杂度就一定是O(MN),在文本比较长的时候,是极为耗用存储空间的。故若要计算出匹配字串,还得用其他的算法,留待后文介绍。
part 2:Java实现
-
首先了解一下这个算法,这是动态规划中的一个经典例子,对于递归的概念的深刻理解,将一个大问题变成小问题,并逐步求解。由第一个字符,我们假设为缺失开始,每一加一个字符,都有三种情况,mismatch,match,deletion/insert,对应的进行打分。高分值的为最优解,逐步延伸至最后一个字符。如有不懂,可以对这个模型在进行深刻理解。
-
该代码对回溯步骤进行简化,仅返回一个最优解,但是该方法简单,适合初级编程者,自学所用。
源代码:
/*Needleman Wunsch是一个全局比对算法,可用于DNA和蛋白质序列的全局比对*/
public class Needleman_Wunsch {
/*全局变量用于回溯是的指针*/
static int l=0;
public static void main(String[] args) {
/*比对的两列字符串*/
String t="GCGCAATG";
String p ="GCCCTAGCG";
/*创建H矩阵用于打分,成为打分矩阵,创建D矩阵用于回溯,成为指针矩阵或者方向矩阵*/
int tlen=t.length();
int plen=p.length();
int [][] h=new int[tlen+1][plen+1];
int [][] d=new int[tlen+1][plen+1];
/*初始化矩阵,第一列或行为deletion后者insert,扣分2*/
for(int i=0;i<1;i++){
for(int j=0;j<plen+1;j++){
h[i][j]=-2*j;
d[i][j]=3;
}
}
for(int j=0;j<1;j++){
for(int i=0;i<tlen+1;i++){
h[i][j]=-2*i;
d[i][j]=1;
}
}
/*动态规划用于打分*/
for(int i=1;i<tlen+1;i++){
for(int j=1;j<plen+1;j++){
/*分值:mismatch(失配)-1,deletion(缺失)/inserting(插入)-2,match(匹配)1,*/
int s1=-2,s2=0,s3=-2;
if(t.charAt(i-1)==p.charAt(j-1)){
s2=1;
}else{
s2=-1;
}
h[i][j]=maximum(h[i-1][j]+s1,h[i-1][j-1]+s2,h[i][j-1]+s3);
d[i][j]=l;
}
}
/*输出打分矩阵*/
System.out.println("score matrix:");
for(int i=0;i<tlen+1;i++){
for(int j=0;j<plen+1;j++){
System.out.printf("%4d",h[i][j]);
if(j!=0&&j%plen==0){
System.out.println();
}
}
}
/*输出索引矩阵*/
System.out.println("index matrix:");
for(int i=0;i<tlen+1;i++){
for(int j=0;j<plen+1;j++){
System.out.print(d[i][j]+" ");
if(j!=0&&j%plen==0){
System.out.println();
}
}
}
/*输出结果*/
System.out.print("Target sequence:");
String result =get_back(t,p,d);
for (int i=0;i<result.length();i++){
System.out.print(result.charAt(i)+" ");
}
System.out.println();
System.out.print("Source sequence:");
for (int i=0;i<p.length();i++){
System.out.print(p.charAt(i)+" ");
}
}
/*求最大值的方法*/
public static int maximum(int a,int b,int c){
int max =a;
l=1;
if(a<b){
max=b;
l=2;
if(b<c){
max=c;
l=3;
}
}else if(a<c){
max=c;
l=3;
}
if(max==a&&max==b){
l=4;
}else if(max==a&&max==c){
l=5;
}else if(max==b&&max==c){
l=6;
}
if(max==a&&max==b&&max==c){
l=7;
}
return max;
}
/*回溯方法*/
public static String get_back(String t,String p,int[][] d){
int i=t.length();
int j=p.length();
StringBuffer sb = new StringBuffer();
while(i>=0&&j>0){
int start = d[i][j];
switch(start){
case 1:sb.insert(0, t.charAt(i-1));i=i-1;break;
case 2:sb.insert(0, t.charAt(i-1));i=i-1;j=j-1;break;
case 3:sb.insert(0, '-');j=j-1;break;
case 4:sb.insert(0, t.charAt(i-1));i=i-1;j=j-1;break;
case 5:sb.insert(0, t.charAt(i-1));i=i-1;break;
case 6:sb.insert(0, t.charAt(i-1));i=i-1;j=j-1;break;
case 7:sb.insert(0, t.charAt(i-1));i=i-1;j=j-1;break;
}
}
String result =sb.toString();
return result;
}
}
-
结果分析:
Target sequence:G C G C - A A T G
Source sequence:G C C C T A G C G
此结果为最优解之一。