在生物应用中,经常需要比较两个(或多个)不同生物体的DNA,
例如:某种生物的DNA可能为S1=ACCGGTCGAGTGCGCGGAAGCCGGCCGAA,
另一种生物的DNA可能为S2=GTCGTTCGGAATGCCGTTGCTCTGTAAA
我们比较两个DNA串的一个原因是希望确定它们的相似度,作为度量两种生物的近似程度指标
寻找第三个串S3,它所有碱基也都出现在S1和S2中,且三个串中的顺序都相同,但在S1和S2中不要求连续出现。
可以找到的S3越长,就可以认为S1和S2的相似度越高。在这个例子中最长的S3为GTCGTCGGAAGCCGGCCGAA
我们定义C[i, j]表示Xi和Yj的LCS长度。如果i = 0或j = 0,即一个序列长度为0,那么LCS的长度为0
根据LCS问题的最优子结构性质,可得如下公式:
C[i, j] = 0,若i = 0 或 j = 0
C[i, j] = C[i - 1, j - 1] + 1,若i,j > 0 且 Xi = Yj
C[i, j] = max(C[i, j - 1], C[i - 1, j]) ,若i, j > 0且Xi != Yj
代码如下:
package 动态规划;
/**
* Lcs即最长公共子序列问题(longest common subsequence problem)
* @author wangdong20
*
*/
public class Lcs {
public static final int empty = 0;
public static final int upLeft = 1;
public static final int up = 2;
public static final int left = 3;
public static int[][][] lcsLength(String x, String y){
int m = x.length();
int n = y.length();
int[][][] result = new int[2][m + 1][n + 1]; // result[0]表示子序列长度 result[1]表示LCS矩阵方向
for(int i = 0; i < m + 1; i++){
result[0][i][0] = 0;
result[1][i][0] = empty;
}
for(int j = 0; j < n + 1; j++){
result[0][0][j] = 0;
result[1][0][j] = empty;
}
for(int i = 1; i <= m; i++){
for(int j = 1; j <= n; j++){
if(x.charAt(i - 1) == y.charAt(j - 1)){
result[0][i][j] = result[0][i - 1][j - 1] + 1;
result[1][i][j] = upLeft;
}
else if(result[0][i - 1][j] >= result[0][i][j - 1]){
result[0][i][j] = result[0][i - 1][j];
result[1][i][j] = up;
}
else{
result[0][i][j] = result[0][i][j - 1];
result[1][i][j] = left;
}
}
}
return result;
}
public static void printLcs(int[][][] b, String x, int i, int j){
if(i == 0 || j == 0)
return;
if(b[1][i][j] == upLeft){
printLcs(b, x, i - 1, j - 1);
System.out.print(x.charAt(i - 1));
}
else if(b[1][i][j] == up){
printLcs(b, x, i - 1, j);
}
else{
printLcs(b, x, i, j - 1);
}
}
/**
* @param args
*/
public static void main(String[] args) {
// TODO 自动生成方法存根
String s1 = "ACCGGTCGAGTGCGCGGAAGCCGGCCGAA";
String s2 = "GTCGTTCGGAATGCCGTTGCTCTGTAAA";
String s3 = "amputation";
String s4 = "spanking";
System.out.println("s1: " + s1);
System.out.println("s2: " + s2);
System.out.println("最长公共子序列: ");
int result[][][] = lcsLength(s1, s2);
printLcs(result, s1, s1.length(), s2.length());
System.out.println("\ns3: " + s3);
System.out.println("s4: " + s4);
System.out.println("最长公共子序列: ");
int result2[][][] = lcsLength(s3, s4);
printLcs(result2, s3, s3.length(), s4.length());
}
}
实质上lcsLength(s3, s4)返回的是两个二维数组组成的三维数组
代码中result[0][i][j]保存的是图中显示的到字符串Xi, Yj目前的LCS长度
result[1][i][j]保存的是图中显示的字符串Xi, Yj的指引方向关系
得到这幅图我们就可以从中得出表b[m, n]
为了得出最后的LCS字符串,只需要从b[m, n]开始,按照箭头方向追踪下去即可。
当b[i, j]遇到upLeft左上时,意味着Xi = Yj是LCS的一个元素.
按照这种方法可以逆序依次构造出LCS的所有元素
public static void printLcs(int[][][] b, String x, int i, int j){
if(i == 0 || j == 0)
return;
if(b[1][i][j] == upLeft){
printLcs(b, x, i - 1, j - 1);
System.out.print(x.charAt(i - 1));
}
else if(b[1][i][j] == up){
printLcs(b, x, i - 1, j);
}
else{
printLcs(b, x, i, j - 1);
}
}
最后运行结果: