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=a₁a₂……a_N，表示A是由a₁a₂……a_N这N个字符组成，Len(A)=N

　　B=b₁b₂……b_M，表示B是由b₁b₂……b_M这M个字符组成，Len(B)=M

　　定义LCS(i,j)=LCS(a₁a₂……a_i,b₁b₂……b_j)，其中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，有公式一

　　若a_i=b_j，则LCS(i,j)=LCS(i-1,j-1)+1

　　若a_i≠b_j，则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矩阵

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
G	0
A	0
T	0
C	0
G	0
A	0

　　第二步：利用公式一，计算矩阵的第一行

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
A	0
T	0
C	0
G	0
A	0

 　　第三步：利用公式一，计算矩阵的其余各行

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

　　则，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算法矩阵，通过回溯，能找到匹配字串

　　第一步：定位在矩阵的右下角

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

　　第二步：回溯单元格，至矩阵的左上角

　　　　若a_i=b_j，则回溯到左上角单元格

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

 　　　　若a_i≠b_j，回溯到左上角、上边、左边中值最大的单元格，若有相同最大值的单元格，优先级按照左上角、上边、左边的顺序

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

　　　　若当前单元格是在矩阵的第一行，则回溯至左边的单元格

　　　　若当前单元格是在矩阵的第一列，则回溯至上边的单元格

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

　　　　依照上面的回溯法则，回溯到矩阵的左上角

　　第三步：根据回溯路径，写出匹配字串

　　　　若回溯到左上角单元格，将a_i添加到匹配字串A，将b_j添加到匹配字串B

　　　　若回溯到上边单元格，将a_i添加到匹配字串A，将_添加到匹配字串B

　　　　若回溯到左边单元格，将_添加到匹配字串A，将b_j添加到匹配字串B

　　　　搜索晚整个匹配路径，匹配字串也就完成了

 

　　可以看出，LD算法和Needleman/Wunsch算法的回溯路径是一样的。这样找到的匹配字串也是一样的。

　　不过，Needleman/Wunsch算法和LD算法一样，若要找出匹配字串，空间的复杂度就一定是O(MN)，在文本比较长的时候，是极为耗用存储空间的。故若要计算出匹配字串，还得用其他的算法，留待后文介绍。

part 2:Java实现

1. 首先了解一下这个算法，这是动态规划中的一个经典例子，对于递归的概念的深刻理解，将一个大问题变成小问题，并逐步求解。由第一个字符，我们假设为缺失开始，每一加一个字符，都有三种情况，mismatch，match，deletion/insert,对应的进行打分。高分值的为最优解，逐步延伸至最后一个字符。如有不懂，可以对这个模型在进行深刻理解。

   
2. 2

   该代码对回溯步骤进行简化，仅返回一个最优解，但是该方法简单，适合初级编程者，自学所用。

   源代码：

   /*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;    

        }    

   }

3. 3

   结果分析：

   Target sequence:G C G C - A A T G 

   Source sequence:G C C C T A G C G

   此结果为最优解之一。