本文介紹基於最長公共子序列的文本比較算法——Needleman/Wunsch算法。還是以實例說明:字符串A=kitten,字符串B=sitting那他們的最長公共子序列為ittn(注:最長公共子序列不需要連續出現,但一定是出現的順序一致),最長公共子序列長度為4。
和LD算法類似,Needleman/Wunsch算法用的都是動態規划的思想,兩者十分相似。
舉例說明:A=GGATCGA,B=GAATTCAGTTA,計算LCS(A,B)。
第一步:初始化動態轉移矩陣
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
狀態轉移方程是:若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算法矩陣,通過回溯,能找到匹配字串
第一步:定位在矩陣的右下角
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算法的回溯路徑是一樣的。這樣找到的匹配字串也是一樣的。