最长公共子序列问题

最长公共子序列问题（LCS）（生物信息学中常用算法）

子序列的概念：   设X = <x₁,x₂,┅, x_m>，若有1≤i₁<i₂<┅ <i_k≤m，使得Z=< z₁, z₂,┅, z_k> =<x_i1, x_i2,┅, x_ik>，则称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。

 

一.分析

记X_i=﹤x₁，…，x_i﹥即X序列的前i个字符 (1≤i≤m)（前缀）,Y_j=﹤y₁，…，y_j﹥即Y序列的前j个字符 (1≤j≤n)（前缀）

假定Z=﹤z₁，…，z_k﹥∈LCS(X , Y)。若x_m=y_n（最后一个字符相同），则不难用反证法证明：

该字符必是X与Y的任一最长公共子序列Z（设长度为k）的最后一个字符，即有z_k = x_m = y_n。且显然有Z_k-1∈LCS(X_m-1 ,Y_n-1)即Z的前缀Z_k-1是X_m-1与Y_n-1的最长公共子序列。

若x_m≠y_n，则亦不难用反证法证明：要么Z∈LCS(X_m-1, Y)，要么Z∈LCS(X , Y_n-1)。

由于z_k≠x_m与z_k≠y_n其中至少有一个必成立，因此：若z_k≠x_m则有Z∈LCS(X_m-1 ,Y)，

若z_k≠y_n 则有Z∈LCS(X , Y_n-1)。  

 ∴若x_m=y_n，则问题化归成求X_m-1与Y_n-1的LCS（LCS(X , Y)的长度等于LCS(X_m-1, Y_n-1)的长度加1）

若x_m≠y_n，则问题化归成求X_m-1与Y的LCS及X与Y_n-1的LCS。

LCS(X , Y)的长度为：Max {LCS(X_m-1, Y)的长度, LCS(X , Y_n-1)的长度}求LCS(X_m-1, Y)的长度与LCS(X , Y_n-1)的长度这两个问题不是相互独立的：

∵两者都需要求LCS(X_m-1，Y_n-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,则沿向上和向左两个方向搜索，