最长公共子序列问题(LCS)(生物信息学中常用算法)
子序列的概念: 设X = <x1,x2,┅, xm>,若有1≤i1<i2<┅ <ik≤m,使得Z=< z1, z2,┅, zk> =<xi1, xi2,┅, xik>,则称Z是X的子序列,记为Z<X。e.g.X=<A,B,C,B,D,A,B>, Z=<B,C,B,A>, 则有Z<X。
公共子序列的概念:
设X,Y是两个序列,且有Z<X和Z<Y,则称Z是X和Y 的公共子序列。
最长公共子序列的概念:若Z<X,Z<Y,且不存在比Z更长的X和Y 的公共子序列,则称Z是X和Y 的最长公共子序列,记为ZÎLCS(X , Y)。
最长公共子序列往往不止一个。e.g. X=<A,B,C,B,D,A,B>, Y=<B,D,C,A,B,A>, 则Z=<B,C,B,A>, Z’=<B,C,A,B>, Z’’=<B,D,A,B>均属于LCS(X , Y),即X,Y有3个LCS。
一.分析
记Xi=﹤x1,…,xi﹥即X序列的前i个字符 (1≤i≤m)(前缀),Yj=﹤y1,…,yj﹥即Y序列的前j个字符 (1≤j≤n)(前缀)
假定Z=﹤z1,…,zk﹥∈LCS(X , Y)。若xm=yn(最后一个字符相同),则不难用反证法证明:
该字符必是X与Y的任一最长公共子序列Z(设长度为k)的最后一个字符,即有zk = xm = yn。且显然有Zk-1∈LCS(Xm-1 ,Yn-1)即Z的前缀Zk-1是Xm-1与Yn-1的最长公共子序列。
若xm≠yn,则亦不难用反证法证明:要么Z∈LCS(Xm-1, Y),要么Z∈LCS(X , Yn-1)。
由于zk≠xm与zk≠yn其中至少有一个必成立,因此:若zk≠xm则有Z∈LCS(Xm-1 ,Y),
若zk≠yn 则有Z∈LCS(X , Yn-1)。
∴若xm=yn,则问题化归成求Xm-1与Yn-1的LCS(LCS(X , Y)的长度等于LCS(Xm-1, Yn-1)的长度加1)
若xm≠yn,则问题化归成求Xm-1与Y的LCS及X与Yn-1的LCS。
LCS(X , Y)的长度为:Max {LCS(Xm-1, Y)的长度, LCS(X , Yn-1)的长度}求LCS(Xm-1, Y)的长度与LCS(X , Yn-1)的长度这两个问题不是相互独立的:
∵两者都需要求LCS(Xm-1,Yn-1)的长度,因而具有重叠性。此外,两个序列的LCS中包含了两个序列的前缀的LCS,故问题具有最优子结构性质考虑用动态规划法。
根据以上分析可知:最长公共子序列问题具有性质:(1)经过分解后的子问题具有高度重复性;(2)具有最优子结构性质。因此可以采用动态规划法求解问题。
二.算法设计
为了构造出LCS,还需要使用一个二维数组b[m][n],b[i][j]记录C[i][j]是通过哪个子问题的值求得的,以决定搜索的方向,欲求出所有的LCS,定义数组b如下:
设1表示“↖”对角线方向,2表示“↑”向上,3表示“←”向左,4表示“←↑”向上或向左,
若X[i]=Y[j],b[i][j] = 1,
若C[i-1][j]>C[i][j-1], 则b[i][j] = 2,
若C[i-1][j]<C[i][j-1], 则b[i][j] = 3,
若C[i-1][j]=C[i][j-1], 则b[i][j] = 4,
根据以上辅助数组C和b的定义,算法首先需要求出这两个数组,C[m][n]中记录的最长公共子序列的长度,b中记录了查找子序列元素的搜索方向。
利用C和b 数组的内容,outputLCS用递归的方式输入全部LCS。笔者实现的方式要用到一个字符数组lcs[]用来存放结果,长度为C[m][n],即LCS的长度lcs_length。还用到一个变量current_length表示当前的位置,初始值为lcs_length。基本思路为:从数组b的右下角往前搜索即outputLCS(b,X,i,j,current_length,lcs_length),其中b是之前求得的b数组,X是X序列,i、j分别是X、Y序列的长度。如果b[i][j]=1,则xi为LCS中的元素,xi存入lcs[current_length],同时将current_length减1,然后沿对角线方向继续搜索,即递归进入outputLCS(b,X,i-1,j-1,current_length,lcs_length);如果b[i][j]=2,则沿向上方向搜索,递归进入outputLCS(b,X,i-1,j,current_length,lcs_length);如果b[i][j]=3,则沿向左方向搜索,递归进入outputLCS(b,X,i,j-1,current_length,lcs_length);如果b[i][j]=4,则沿向上和向左两个方向搜索,先后递归执行outputLCS(b,X,i-1,j,current_length,lcs_length)和outputLCS(b,X,i,j-1,current_length,lcs_length)。递归结束的标志是搜索到了边界即i=0或j=0时,此时输出lcs[]数组的内容即为LCS中的一个,然后结束本层递归,返回上一层。上面的每一次递归进入时都要先判断是否搜索到了边界。因为算法在沿不同路径搜索时可能会出现相同的LCS序列,因此要进行去重处理,本算法的处理方法是,将每一次的lcs存入到一个ArrayList中,如果ArrayList中不存在lcs,则将其加入,否则,不加入。最后输出ArrayList中的所有元素即可得到所有的LCS,ArrayList的长度是LCS的个数。
源代码如下:
/**
* 功能:最长公共子序列(LCS)问题的代码实现
*/
import java.util.ArrayList;
public class LCS {
//定义存放两个序列的字符数组
private char seqX[]=null;
private char seqY[]=null;
//定义求解过程用到的两个二维数组C数组和b数组
private int cArray[][]=null;
private int bArray[][]=null;
private char[] lcs=null;
//定义用于存放最终结果的ArrayList对象
private ArrayList<String> lcsArray=new ArrayList<String>();
//构造函数,传入两个字符串对象
public LCS(String str1,String str2){
this.seqX=str1.toCharArray();
this.seqY=str2.toCharArray();
}
//初始化两个数组
private void initArray(){
int row=seqX.length;
int col=seqY.length;
int[][] cArray=new int[row+1][col+1];
int[][] bArray=new int[row+1][col+1];
for(int i=0;i<row+1;i++){
for(int j=0;j<col+1;j++){
cArray[i][j]=0;
bArray[i][j]=0;
}
}
this.cArray=cArray;
this.bArray=bArray;
}
//初始化存放结果的数组
private char[] initRusult(int length){
char[] lcs=new char[length];
for(int i=0;i<length;i++){
lcs[i]='%';
}
return lcs;
}
/**功能:动态规划法求两个数组
*参数:
* seqX:序列X
* seqY:序列Y
**/
private void getLcs(char[] seqX,char[] seqY){
int row=seqX.length;
int col=seqY.length;
for(int i=1;i<row+1;i++){
for(int j=1;j<col+1;j++){
if(seqX[i-1]==(seqY[j-1])){
cArray[i][j]=cArray[i-1][j-1]+1;
bArray[i][j]=1;//1表示“↖”,2表示“↑”,3表示“←”,4表示“←↑”,
}else{
if(cArray[i-1][j]>cArray[i][j-1]){
cArray[i][j]=cArray[i-1][j];
bArray[i][j]=2;//1表示“↖”,2表示“↑”,3表示“←”,4表示“←↑”,
}else if(cArray[i-1][j]<cArray[i][j-1]){
cArray[i][j]=cArray[i][j-1];
bArray[i][j]=3;//1表示“↖”,2表示“↑”,3表示“←”,4表示“←↑”,
}else{
cArray[i][j]=cArray[i-1][j];
bArray[i][j]=4;//1表示“↖”,2表示“↑”,3表示“←”,4表示“←↑”,
}
}
}
}
}
/**功能:递归得到所有LCS
*参数:
* bArray:标志数组
* seqX:序列x
* i:序列X的长度
* j:序列Y的长度
* current_len:当前长度
* lcs_max_len:最长公共子序列长度
**/
private void outputLCS(int[][] bArray,char[] seqX,int i ,int j,int current_len,int lcs_max_len){
if(i==0 || j==0){
String result="";
for(int k=0;k<lcs_max_len;k++)
{
result+=lcs[k];
}
if(result!=null){
if(!lcsArray.contains(result))//如果ArrayList中不存在此lcs则加入
lcsArray.add(result);
}
return;
}
if(bArray[i][j]==1){
current_len--;
lcs[current_len]=seqX[i-1];
outputLCS(bArray, seqX, i-1, j-1,current_len,lcs_max_len);
}else if(bArray[i][j]==2){
outputLCS(bArray, seqX, i-1, j,current_len,lcs_max_len);
}else if(bArray[i][j]==3){
outputLCS(bArray, seqX, i, j-1,current_len,lcs_max_len);
}else if(bArray[i][j]==4){
outputLCS(bArray, seqX, i-1, j,current_len,lcs_max_len);
outputLCS(bArray, seqX, i, j-1,current_len,lcs_max_len);
}
}
private void outputLCS(int[][] bArray,char[] seqX,int i ,int j){
if(i==0 || j==0){
System.out.print("|");
return;
}
if(bArray[i][j]==1){
outputLCS(bArray, seqX, i-1, j-1);
System.out.print(seqX[i-1]);
}else if(bArray[i][j]==2){
outputLCS(bArray, seqX, i-1, j);
}else if(bArray[i][j]==3){
outputLCS(bArray, seqX, i, j-1);
}else if(bArray[i][j]==4){
outputLCS(bArray, seqX, i-1, j);
outputLCS(bArray, seqX, i, j-1);
}
}
/**
* 功能:调用上面的方法输出最后的结果,作为对外接口供用户使用
*/
public void outResult(){
//初始化数组
initArray();
//求数组C和b
getLcs(seqX, seqY);
int lcs_max_length=cArray[seqX.length][seqY.length];
lcs=initRusult(lcs_max_length);
outputLCS(bArray,seqX,seqX.length,seqY.length,lcs_max_length,lcs_max_length);
System.out.println("公共子序列共有:"+lcsArray.size());
for(int i=0;i<lcsArray.size();i++){
System.out.println(lcsArray.get(i));
}
}
}
/**
* 功能:LCS类测试
*/
public class LCSTest {
public static void main(String[] args) {
String str1=" ABCBDAB ";
String str2=" BDCABA ";
System.out.println("序列X:"+str1);
System.out.println("序列Y:"+str2);
//实例化LCS类对象
LCS lcsTest=new LCS(str1,str2);
//调用类提供的接口输出结果
lcsTest.outResult();
}
}
三.运行测试
1.测试用例1
Stringstr1="ABCBDAB";
Stringstr2="BDCABA";
(1) 运行结果为:
(2) 结果分析:
上图为数组C和数组b中的内容。
以上两张数组中的内容可以转换为如下的信息:
递归调用outputLCS(),第一次b[7][6]=4,则要分别向上和向左两个方向搜索,先向上,b[6][6]=1,则C[6][6] 为公共子序列中的元素,存入数组中,此时current_length=4,则lcs[3]=’A’,然后逐步递归,下图中红色圆圈圈出的位置的元素会依次存入lcs[]中,最后到达C[0][0]时,i=0结束本次递归,此时lcs中的内容为:’B’、’C’、’B’、’A’,则”BCBA”被存入到ArrayList中。
回溯到C[7][6]位置,向左搜索,b[7][5]=1,则C[7][5] 为公共子序列中的元素,存入数组中,此时current_length=4,则lcs[3]=’B’,然后逐步递归。当搜索到C[5][3]时,即下图中蓝色圆圈圈中的位置,b[5][3]=1,又要分为两个方向搜索,此时lcs[3]=’B’,lcs[2]=’A’, 此current_length的值为2。先往上搜索,到达‘边界’时,lcs[1]=’C’,lcs[0]=’B’。然后”BCAB”被存入到ArrayLIst中。回溯到C[5][3]位置,向左搜索,此时current_length的值为2,那么后面的C[5][2]、C[4][1]位置的元素会分别被存入lcs[1]、lcs[0]中,在C[3][0]处到达边界,此时lcs中为”BDAB”,存入ArrayList中。整个递归过程结束。
从上述的分析中可以看出,本算法很好的处理了搜索过程中的分支问题。在没有用此方法之前,我用的是如下的方法输出结果:
private void outputLCS(int[][] bArray,char[] seqX,int i ,int j){
if(i==0 || j==0){
System.out.print("|");//分割输出结果
return;
}
if(bArray[i][j]==1){
outputLCS(bArray, seqX, i-1, j-1);
System.out.print(seqX[i-1]);
}else if(bArray[i][j]==2){
outputLCS(bArray, seqX, i-1, j);
}else if(bArray[i][j]==3){
outputLCS(bArray, seqX, i, j-1);
}else if(bArray[i][j]==4){
outputLCS(bArray, seqX, i-1, j);
outputLCS(bArray, seqX, i, j-1);
}
}
输出的结果为:
(其中的‘|’是用来分割结果的,如果不分开的话结果是一串字符,无法辨认)。其中BCBA和BDAB两个结果是正确的,而中间只有BC两个字符,通过对问题和算法的分析,发现,由于C[5][3]位置产生了分支,回溯后沿向左方向搜索,C[6][4]和C[7][5]位置的元素只被输出了一次,但是这两个元素是共有的,却只输出一次,所以导致结果不够清晰。当然,对这个长串进行分割、合并、去重处理,也能得到正确的结果。
2.测试用例2
String str1="asjfdljzifjqernljcv";
String str2="ajfkljewiurzjvc";
(1) 运行结果为:
(2) 结果分析:
下图给出了数组C和数组b中的内容。由于测试用例较长,本实例不对全部进行详细分析,只取出部分进行分析。
这个测试用例是我在键盘上随机敲击形成的,旨在测试程序对较长序列的适用性。可以观测测试结果是正确的。
依然用未改进的算法输出:得到的是一个相当长的结果:|ajf|ajflj|ajf|ajflj|ajf|ajflji|ajf|ajflj|ajf|ajflj|ajf|ajflji|ajf|ajflj|ajf|ajflj|ajf|ajflji|ajf|ajflj|ajf|ajflj|ajf|ajflji|ajf|ajflj|ajf|ajflj|ajf|ajflji|ajf|ajfljer|ajf|ajflj|ajf|ajflj|ajf|ajflji|ajf|ajflj|ajf|ajflj|ajf|ajflji|ajf|ajflj|ajf|ajflj|ajf|ajflji|ajf|ajflj|ajf|ajflj|ajf|ajflji|ajf|ajflj|ajf|ajflj|ajf|ajflji|ajf|ajfljer|ajf|ajflj|ajf|ajflj|ajf|ajflji|ajf|ajflj|ajf|ajflj|ajf|ajflji|ajf|ajflj|ajf|ajflj|ajf|ajflji|ajf|ajflj|ajf|ajflj|ajf|ajflji|ajf|ajflj|ajf|ajflj|ajf|ajflji|ajf|ajfljerjc|ajf|ajflj|ajf|ajflj|ajf|ajflji|ajf|ajflj|ajf|ajflj|ajf|ajflji|ajf|ajflj|ajf|ajflj|ajf|ajflji|ajf|ajflj|ajf|ajflj|ajf|ajflji|ajf|ajflj|ajf|ajflj|ajf|ajflji|ajf|ajfljer|ajf|ajflj|ajf|ajflj|ajf|ajflji|ajf|ajflj|ajf|ajflj|ajf|ajflji|ajf|ajflj|ajf|ajflj|ajf|ajflji|ajf|ajflj|ajf|ajflj|ajf|ajflji|ajf|ajflj|ajf|ajflj|ajf|ajflji|ajf|ajfljer|ajf|ajflj|ajf|ajflj|ajf|ajflji|ajf|ajflj|ajf|ajflj|ajf|ajflji|ajf|ajflj|ajf|ajflj|ajf|ajflji|ajf|ajflj|ajf|ajflj|ajf|ajflji|ajf|ajflj|ajf|ajflj|ajf|ajflji|ajf|ajfljerjv。其中很多的字串是重复的,由实例1的分析可知,在进行搜索时肯定产生了很多的分支。
取出数组中的一个片段进行分析,如下图所示。
可以看出,图中蓝色方框圈中的部分C数组的元素都是7,从右下角到左上角每一次都有分支,也就是在第五部分时间复杂度分析中提到的最坏情况,当然只是局部的,此处是一个3×4的矩阵,路径共有C(2+3,2) =10条,仅仅此出的一个片段就用10条路径,也就是输出时此片段后面的LCS的前半部分子串会被输出10次。而且从C数组和b数组的结果可以观测,矩阵中类似的片段有很多,因此输出结果那么长也就不足为奇了。
3.测试用例3
String str1="128aql934&*7&*c(&a%1%&&89%&%27e5%$^&@2859";
String str2="28&^%93503a4$*&c*(682c%$@*(@908d5690e238";
(1) 运行结果为:
(2) 结果分析:
由于本测试实例较长,C数组和b数组都较大,故此处不再给出。
此测试用例也是随机产生的,而且含有不同类型的字符,包含字母、数字、符号,旨在测试程序对不同类型字符的适用性,从结果可以看出,本程序依然能完美的解决此问题。
四.时间复杂度分析
有上述的对outputLCS()算法的分析可知,求出所有的LCS即是根据搜索的方向信息遍历所有的路径找出满足条件的元素集合。求出辅助数组C和b所用的时间为O(mn+m+n)。求LCS时,在求某一个LCS的过程中,遇到b[i][j]=1时,要将序列中的字符放入数组中,时间为O(1),这样的操作要做lcs_length次,到达边界时,要把lcs[]数组中的结果放入ArrayList中,在常数时间内可以完成。这些操作可以分为两类,一是在搜索过程中,碰到每个点的操作,二是到达边界时的操作。利用会计方法将每个点中的操作的代价放入到达边界时操作的代价里面,这样达到一次边界要花费的代价是lcs_length+1,因此只要求出遍历的路径数就可以得到整体的花费时间。而路径数是由搜索方向决定的。显然算法在最好的情况下,即m=n时且b中每次搜索都沿着对角线方向,时间复杂度为O(n)。
于此相反,数组b中所有的值都指示要沿着两个方向(向上或向左)搜索,这种情况下每一个点都要沿着两个方向分别调用outputLCS(),遇到i=0或j=0时返回,直到搜索完所有的路径才能结束。而且这种情况下,每次点上没有操作,这样每搜索一次只有在边界上的操作,时间为O(1),因此路径数即为最坏情况的时间。最坏情况下的搜索矩阵如下图所示(C矩阵中元素全为0)。由此可知,确定最坏情况的时间复杂度,可以转化为计算从点(m,n)到i=0或j=0(除(i,j)=(0,0)外)的所有路径。如图建立坐标系,对于坐标系中的点可以沿向上或向左两个方向前进设点S(m,n),所求问题即为点 S(m,n)到x轴上点或y轴上点(除原点O(0,0)外)所有路径之和。
现假设x 轴上的点 均可以向上前进, y轴上的系列坐标点 均可以向左前进。则这些点到原点的路径是唯一的。因此 经过x轴上点 或y轴上点到达原点 的路径点S(m,n)与 到达x轴上点 或y轴上点的路径数相同。由此,所求问题有转化为了,从点S(m,n)开始只能向上或向左方向前进(边界上只有一个方向)到达原点O(0,0)的路径数。求解如下:
把m次向上的移动记为:Um、Um-1、….、U2、U1,把n次向左的移动记为:Ln、Ln-1、……、L2、L1。将上述(m+n)个不同的元素,按从大到小的顺序混合排列起来,即为上述问题的答案。有排列组合知识可知,从(m+n)个位置中选出n个给L,即有 C(m+n,n)种方法,所以所求的排列数为C(m+n,n) 。所以算法最坏情况的时间复杂度为 C(m+n,n)。

被折叠的 条评论
为什么被折叠?



