For example, the DNA of one organism may be
S1= { ACCGGTCGAGTGCGCGGAAGCCGGCCGAA}
S2= {GTCGTTCGGAATGCCGTTGCTCTGTAAA}
one goal of comparing two strands of DNA is to determine how "similar" the two strands are, as some measure of how closely related the two organisms are."
------------------------------------------摘自Introduction to algorithms[算法导论]
现在,我们遇到求公共子序列的时候,公共子序列意味着什么?匹配?相似度?怎么求最长公共子序列(longest common subsequence)?本文将告诉你其算法并给出c/c++ 与Java实现源代码.
[1] 定义:何为最长公共子序列?
eg: X = {A,B,C,B,D,A,B} Y = {B,D,C,A,B,A}此两序列的最长公共子序列是LCS={B,C,B,A}
定义我就不多说了,自己感受咯.
[2] 一个相关的定理:
若X = {x1,x2,...,xm} , Y = {y1,y2,...,yn}的LCS是Z={z1,z2,...,zk}那么有:
1) 如果xm = yn , 则zk = xm = yn 并且 Zk-1 是Xm-1和Yn-1的LCS.
2)如果xm ≠ yn, 那么zk ≠ xm意味着Z 是Xm-1和Y的LCS.
3)如果xm ≠ yn, 那么zk ≠ yn 意味着Z是X和Yn-1的LCS.
这个定理其实不难理解,用反证法可以证明之.但是却暗含递归思想.
let us define c[i,j] to be the length of an LCS of the sequence Xi and Yj.
整理一下.有如下的递归关系:
[3]算法源程序
#include <string.h>
int* LCS_length(char* X, char* Y)
{
int m = strlen(X) + 1;
int n = strlen(Y) + 1;
int (*c)[n] = new int[m][n];
int (*b)[n] = new int[m][n];
for(int i=0;i<m;++i)
for(int j=0;j<n;++j)
b[i][j] = 0;
for(int i=0;i<m;++i)
c[i][0] = 0;
for(int i=0;i<n;++i)
c[0][i] = 0;
for(int i=1;i<m;++i)
for(int j=1;j<n;++j)
{
if (X[i-1] == Y[j-1])
{
c[i][j] = c[i-1][j-1] +1 ;
b[i][j] = 2; // 比较相等之标记
}
else
{
if( c[i-1][j] >= c[i][j-1])
{
c[i][j] = c[i-1][j];
b[i][j] = 1; //
}
else
c[i][j] = c[i][j-1];
}
}
delete[] c;
return (int*)b;
} void LCS(char* X,char* Y,int m, int n,int* b)
{
if(m ==0 || n == 0)
return;
if(b[m* (strlen(Y)+1)+ n] == 2)
{
LCS(X,Y,m-1,n-1,b);
printf("%c",X[m-1]);
}
else if(b[m* (strlen(Y)+1)+ n] == 1)
{
LCS(X,Y,m-1,n,b);
}
else
{
LCS(X,Y,m,n-1,b);
}
} int main()
{
// char* X = "ABCBDAB";
// char* Y = "BDCABA";
char* X = "ACCGGTCGAGTGCGCGGAAGCCGGCCGAA";
char* Y = "GTCGTTCGGAATGCCGTTGCTCTGTAAA";
int m = strlen(X) + 1;
int n = strlen(Y) + 1;
int* c = NULL;
c = LCS_length(X,Y);
LCS(X,Y,strlen(X),strlen(Y),c);
delete[] c;
getchar();
return 0;
}
以上程序在VC6.0下是不能正常编译的.推荐用Dev-C++来调试.我用的版本是4.9.9.0.以上程序或许有点费解,主要是C++对动态的多维数组不支持! 您看,我到求LCS时变二维为一维来求了,但愿你能根据前面的定理看懂我这糟糕的程序咯.
那么Java怎么样?JAVA支持多维数组啦,因此java写起来"好看"多了.^_^
{
final static int NorthWest = 2;
final static int UP = 1;
final static int LEFT = 0;
private static int[][] LCS_length(String X, String Y)
{
int m = X.length()+1;
int n = Y.length()+1;
int[][] c = new int[m][n];
int[][] flag = new int[m][n];
for(int i=0;i<m;++i)
for(int j=0;j<n;++j)
flag[i][j] = LEFT;
for(int i=0;i<m;++i) //递归表达式2
c[i][0] = 0;
for (int i=0;i<n;++i)
c[0][i] = 0;
for(int i=1;i<m;++i)
for(int j=1;j<n;++j)
{
if( X.charAt(i-1) == Y.charAt(j-1))
{
c[i][j] = c[i-1][j-1] +1;//递归表达式3
flag[i][j] = NorthWest; //
}
else
{
if(c[i-1][j] > c[i][j-1])
{
c[i][j] = c[i-1][j];
flag[i][j] = UP; //
}
else
{
c[i][j] = c[i][j-1];
}
}
}
return flag;
}
private static void LCS(String X, int m,int n,int[][] flag)
{
if(m ==0 || n ==0) return;
if(flag[m][n] ==NorthWest)
{
LCS(X,m-1,n-1,flag);
System.out.print(X.charAt(m-1));
}
else if(flag[m][n] == UP)
{
LCS(X,m-1,n,flag);
}
else
{
LCS(X,m,n-1,flag);
}
}
public static void main(String[] args)
{
//String X = "ABCBDAB";
//String Y = "BDCABA";
String X = "ACCGGTCGAGTGCGCGGAAGCCGGCCGAA";
String Y = "GTCGTTCGGAATGCCGTTGCTCTGTAAA";
int m = X.length();
int n = Y.length();
int[][] flag = new int[m+1][n+1];
flag = LCS_length(X,Y);
LCS(X,m,n,flag);
//程序输出:GTCGTCGGAAGCCGGCCGAA
}
}
本文完.如有问题欢迎留言讨论.
参考资料:<Introduction to algorithms>