在生物应用中,经常需要比较两个(或多个)不同生物体的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(" s3: " + 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); } }
最后运行结果: