文本比较算法:Needleman/Wunsch算法

本文介绍基于最长公共子序列的文本比较算法——Needleman/Wunsch算法。还是以实例说明:字符串A=kitten,字符串B=sitting那他们的最长公共子序列为ittn(注:最长公共子序列不需要连续出现,但一定是出现的顺序一致),最长公共子序列长度为4。

和LD算法类似,Needleman/Wunsch算法用的都是动态规划的思想,两者十分相似。

举例说明:A=GGATCGA,B=GAATTCAGTTA,计算LCS(A,B)。

第一步:初始化动态转移矩阵

Needleman/Wunsch算法矩阵
  GAATTCAGTTA
  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算法矩阵
  GAATTCAGTTA
  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算法矩阵
  GAATTCAGTTA
  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

状态转移方程是:若A(i)=B(j),LCS(i,j)=LCS(i-1,j-1)+1;否则LCS(i,j)=max(LCS(i-1,j-1),LCS(i,j-1),LCS(i-1,j))=max(LCS(i,j-1),LCS(i-1,j))。程序实现:

/*
 *侯凯,2014-9-15
 *功能:最长子序列
 */
#include<iostream>
using namespace std;

int CalTheDistance(string A,string B)
{
    int **ptr = new int*[ A.size()+ 1];
    for(int i = 0; i < A.size() + 1 ;i++)
    {
        ptr[i] = new int[B.size() + 1];
    }

    for(int i=0;i<A.size()+1;i++)
    {
        ptr[i][0] = 0;
    }
    for(int i=0;i<B.size()+1;i++)
    {
        ptr[0][i] = 0;
    }
    for(int i=0;i<A.size();i++)
    {
        for(int j=0;j<B.size();j++)
        {
            if(A[i]==B[j])
                ptr[i+1][j+1]=ptr[i][j]+1;
            else
                ptr[i+1][j+1]=max(ptr[i+1][j],ptr[i][j+1]);
        }
    }
    int result = ptr[A.size()][B.size()];
    for(int i = 0; i < A.size() + 1 ;i++)
    {
        delete [] ptr[i];
        ptr[i] = NULL;
    }
    delete[] ptr;
    ptr = NULL;
    return result;
}

int main()
{
    string str1 = "GGATCGA";
    string str2 = "GAATTCAGTTA";
    //最长子序列为6
    int distance = CalTheDistance(str1,str2);
    cout<<distance<<endl;
    system("Pause");
}

以上面为例A=GGATCGA,B=GAATTCAGTTA,LCS(A,B)=6

他们的匹配为:

A:GGA_TC_G__A

B:GAATTCAGTTA

如上面所示,蓝色表示完全匹配,黑色表示编辑操作,_表示插入字符或者是删除字符操作。如上面所示,蓝色字符有6个,表示最长公共子串长度为6。

利用上面的Needleman/Wunsch算法矩阵,通过回溯,能找到匹配字串

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

 

Needleman/Wunsch算法矩阵
  GAATTCAGTTA
  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,则回溯到左上角单元格

 

Needleman/Wunsch算法矩阵
  GAATTCAGTTA
  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,回溯到左上角、上边、左边中值最大的单元格,若有相同最大值的单元格,优先级按照左上角、上边、左边的顺序

 

Needleman/Wunsch算法矩阵
  GAATTCAGTTA
  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算法矩阵
  GAATTCAGTTA
  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算法的回溯路径是一样的。这样找到的匹配字串也是一样的。

posted @ 2014-09-15 17:05  侯凯  阅读(5433)  评论(0编辑  收藏  举报