分类：
2010-07-15 15:58:32
[基因序列相似度,字符串相似度]
LCS问题就是求两个字符串最长公共子串的问题。解法就是用一个矩阵来记录两个字符串中所有位置的两个字符之间的匹配情况，若是匹配则为1，否则为0。然后求出对角线最长的1序列，其对应的位置就是最长匹配子串的位置。
下面是字符串21232523311324和字符串312123223445的匹配矩阵，前者为X方向的，后者为Y方向的。不难找到，红色部分是最长的匹配子串。通过查找位置我们得到最长的匹配子串为：21232
0 0 0 1 0 0 0 1 1 0 0 1 0 0 0
0 1 0 0 0 0 0 0 0 1 1 0 0 0 0
1 0 1 0 1 0 1 0 0 0 0 0 1 0 0
0 1 0 0 0 0 0 0 0 1 1 0 0 0 0
1 0 1 0 1 0 1 0 0 0 0 0 1 0 0
0 0 0 1 0 0 0 1 1 0 0 1 0 0 0
1 0 1 0 1 0 1 0 0 0 0 0 1 0 0
1 0 1 0 1 0 1 0 0 0 0 0 1 0 0
0 0 0 1 0 0 0 1 1 0 0 1 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 1 0
0 0 0 0 0 0 0 0 0 0 0 0 0 1 0
0 0 0 0 0 1 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
但是在0和1的矩阵中找最长的1对角线序列又要花去一定的时间。通过改进矩阵的生成方式和设置标记变量，可以省去这部分时间。下面是新的矩阵生成方式：
0 0 0 1 0 0 0 1 1 0 0 1 0 0 0
0 1 0 0 0 0 0 0 0 2 1 0 0 0 0
1 0 2 0 1 0 1 0 0 0 0 0 1 0 0
0 2 0 0 0 0 0 0 0 1 1 0 0 0 0
1 0 3 0 1 0 1 0 0 0 0 0 1 0 0
0 0 0 4 0 0 0 2 1 0 0 1 0 0 0
1 0 1 0 5 0 1 0 0 0 0 0 2 0 0
1 0 1 0 1 0 1 0 0 0 0 0 1 0 0
0 0 0 2 0 0 0 2 1 0 0 1 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 1 0
0 0 0 0 0 0 0 0 0 0 0 0 0 1 0
0 0 0 0 0 1 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
不用多说，你大概已经看出来了。当字符匹配的时候，我们并不是简单的给相应元素赋上1，而是赋上其左上角元素的值加一。我们用两个标记变量来标记矩阵中值最大的元素的位置，在矩阵生成的过程中来判断当前生成的元素的值是不是最大的，据此来改变标记变量的值，那么到矩阵完成的时候，最长匹配子串的位置和长度就已经出来了。
这样做速度比较快，但是花的空间太多。我们注意到在改进的矩阵生成方式当中，每生成一行，前面的那一行就已经没有用了。因此我们只需使用一维数组即可。最终的代码如下：(其中是s1是横向，s2纵向）
Private Function LCS(ByVal str_1 As String, ByVal str_2 As String) As String
If str_1 = "" Or str_2 = "" Then Return ""
Dim c(str_1.Length) As Integer
Dim max, maxj, i, j As Integer
maxj = 0 : max = 0 '这两个是标志变量
For i = 0 To str_2.Length - 1
For j = str_1.Length - 1 To 0 Step -1
If str_2.Chars(i) = str_1.Chars(j) Then
If i = 0 Or j = 0 Then
c(j) = 1
Else
c(j) = c(j - 1) + 1
End If
Else
c(j) = 0
End If
If c(j) > max Then '把>改成>=则返回最后一个最长匹配子串
max = c(j) : maxj = j '更新标志变量
End If
Next
Next
If max = 0 Then Return ""
Return str_1.Substring(maxj - max + 1, max) '直接从标志变量得出结果
End Function
//////////////////////////////////////////////////////////
附一个c++版的：
string LCS(string s1,string s2)
{
if(s1==""||s2=="")
return "";
int c[s1.size()];
int max,maxj,i,j;
maxj=0;max=0;
for(i=0;i<=s2.size()-1;i++)
for(j=s1.size()-1;j>=0;j--)
{
if(s2[i]==s1[j])
{
if(i==0||j==0)
c[j]=1;
else
c[j]=c[j-1]+1;
}
else
c[j]=0;
if(c[j]>max)
{
max=c[j];
maxj=j;
}
}
if(max==0)
return "";
string ss="";
for(j=maxj-max+1;j<=maxj;j++)
ss+=s1[j];
return ss;
}
=================================================
（非连续的）
首先将要看到如何运用动态编程查找两个 DNA 序列的最长公共子序列（longest common subsequence，LCS）。发现了新的基因序列的生物学家通常想知道该基因序列与其他哪个序列最相似。查找 LCS 是计算两个序列相似程度的一种方法：LCS 越长，两个序列越相似。
子序列中的字符与子字符串中的字符不同，它们不需要是连续的。例如，ACE 是 ABCDE 的子序列，但不是它的子字符串。请看下面两个 DNA 序列：
这两个序列的 LCS 是 GCCAG。（请注意，这仅是一个 LCS，而不是唯一的 LCS，因为可能存在其他长度相同的公共子序列。这种最优化问题和其他最优化问题的解可能不止一个。）
首先，考虑如何递归地计算 LCS。令：
有三个递归子问题：
结果表明（而且很容易使人相信）原始问题的解就是下面三个子序列中最长的一个：
（基线条件（base case）是 S1 或 S2 为长度为零的字符串的情形。在这种情况下，S1 和 S2 的 LCS 显然是长度为零的字符串。）
但是，就像计算斐波纳契数的递归过程一样，这个递归解需要多次计算相同的子问题。可以证明，这种递归解法需要耗费指数级的时间。相比之下，这一问题的动态编程解法的运行时间是 Θ(mn)，其中 m 和 n 分别是两个序列的长度。
为了用动态编程有效地计算 LCS，首先需要构建一个表格，用它保存部分结果。沿着顶部列出一个序列，再沿着左侧从上到下列出另一个序列，如图 2 所示：
这种方法的思路是：将从上向下、从左到右填充表格，每个单元格包含一个数字，代表该行和该列之前的两个字符串的 LCS 的长度。也就是说，每个单元格包含原始问题的一个字问题的解。例如，请看第 6 行第 7 列的单元格：它在 GCGCAATG 序列第二个 C 的右侧，在 GCCCTAGCG 的 T 的下面。这个单元格最终包含的数字就是 GCGC 和 GCCCT 的 LCS 的长度。
首先看一下表格的第二行中应该是什么条目。这一行的单元格保存的 LCS 长度对应的是序列 GCGCAATA 的零长前端和序列 GCCCTAGCG 的 LCS。显然，这些 LCS 的值都是 0。类似的，沿着第二列向下的值也都是 0，这与递归解的基线条件对应。现在表格如图 3 所示：
接下来，要实现与递归算法中递归子问题对应的情景，但这时使用的是表格中已经填充的值。在图 4 中，我已经填充了一半左右的单元格：
在填充单元格时，需要考虑以下条件：
下面三个值分别对应着我在前面列出的三个递归子问题返回的值。
在空单元格中填充下面 3 个数字中的最大值：
请注意，我在图中还添加了箭头，指向当前单元格值的来源。后面的 “回溯” 一节将用这些箭头建立实际的 LCS（与仅仅发现 LCS 长度相反）。
现在填充图 4 中接下来的空单元格 — 在 GCCCTAGCG 中第三个 C 下面和 GCGCAATG 第二个 C 的右侧的单元格。它上面的值是 2，左侧的值是 3，左上侧的值是 2。这个单元格上面的字符和左侧的字符相等（都是 C），所以必须选择 2、3 和 3（左上侧单元格中的 2 + 1）的最大值。所以，这个单元格的值为 3。绘制一个箭头，从该单元格指向其中的值的源单元格。在这个示例中，新的数值可能来自不止一个单元格，所以可以任选一个：例如左上侧单元格。
作为练习，您可以尝试填充表格的余下部分。如果在关联过程中，一直按照左上侧-上侧-左侧的顺序选择单元格，那么会得到如图 5 所示的表格。（当然，如果在关联过程中做了不同的选择，那么箭头就会不同，但是数字是相同的。）
回想一下，任何单元格中的数字都是该单元格所在行之上和列之前的字符串的 LCS 长度。所以，表格右下角的数字就是字符串 S1 和 S2 （在本例中是 GCCCTAGCG 和 GCGCAATG）的 LCS 的长度。所以，这两个序列的 LCS 长度是 5。
这是在所有动态编程算法中都要牢记的关键点。表格中的每个单元格都包含单元格所在行上面和所在列左侧序列前端问题的解。
接下来要做的就是寻找实际的 LCS。使用单元格箭头进行回溯可以完成。在构建表格的时候，请记住，如果箭头指向左上侧的单元格，那么当前单元格中的值要比左上侧单元格的值大 1，这意味着左侧单元格和上面单元格中的字符相等。构建 LCS 时，这会将相应的字符添加到 LCS 中。所以，构建 LCS 的途径就是从右下角的单元格开始，沿着箭头一路返回。每当沿着对角箭头回到左上角的单元格而且 该单元格的值比当前单元格的值小 1 时，就要将对应的公共字符添加到 正在构建的 LCS 的前端。请注意，之所以将字符放在 LCS 前端，是因为我们是从 LCS 末端开始的。（在 图 5 的示例中，右下角的 5 与要添加的第 5 个字符对应。）
依此类推，继续构建 LCS。从右下侧的单元格开始，看到单元格指针指向左上侧的单元格，而且当前单元格的值（5）比其左上侧单元格的值（4）大 1。所以将字符 G 添加到最初的零长度的字符串之前。下一个箭头，从包含 4 的单元格开始，也指向左上侧，但是值没有变。接着这个箭头也是如此。下一个单元格的箭头还是指向左上侧，但是这次值从 3 变为 4。这意味着需要将这一行和这一列中的公共字符 A 添加到 LCS 中。所以现在的 LCS 是 AG。接下来，沿着指针向左（对应着跳过上面的 T）到达另一个 3。然后有一个对角指针指向 2。因此，又添加了在当前行和当前列中的公共字符 C，生成的 LCS 为 CAG。继续使用这种方式，直到最后到达 0。图 6 显示了整个回溯过程：
通过这种回溯方法，得到的 LCS 为 GCCAG
python的代码实现：
def lcs_similar(x, y):
x = unicode(x, 'utf-8', 'ignore')
y = unicode(y, 'utf-8', 'ignore')
m = len(x) + 1
n = len(y) + 1
if m == 1 and n == 1:
return 100
elif m==1 or n == 1:
return 0
lcs_len = [[0 for a in range(n)] for b in range(m)]
for i in range(1, m):
for j in range(1, n):
if(x[i-1] == y[j-1]):
lcs_len[i][j] = lcs_len[i-1][j-1] + 1
else:
max_len = lcs_len[i-1][j]
if lcs_len[i][j-1] > lcs_len[i-1][j]:
max_len = lcs_len[i][j-1]
lcs_len[i][j] = max_len
i = m - 2
j = n - 2
common_substr = u''
while True:
if i == -1 or j == -1:
break
if x[i] == y[j]:
common_substr = u"%s%s" % (x[i], common_substr)
i = i - 1
j = j -1
else:
if lcs_len[i][j+1] > lcs_len[i+1][j]:
i = i-1
else:
j = j-1
print common_substr
//////////////////////////////////////////////////////////////////
添加一个C++版本的：
string LCS_2(string s1,string s2)
{
if(s1==""||s2=="")
return "";
int m=s1.size()+1;
int n=s2.size()+1;
int lcs[m][n],i,j;
for(i=0;i
for(i=1;i
if(s1[i-1]==s2[j-1])
lcs[i][j]=lsc[i-1][j-1]+1;
else
lcs[i][j]=lcs[i-1][j]>=lcs[i][j-1]?lcs[i-1][j]:lcs[i][j-1];//取上面或左侧最大值
}
i=m-2;
j=n-2;
string ss=""
while(i!=-1&&j!=-1)
{
if(x[i]==y[j])
{
ss+=s[i];
i--;
j--;
}
else
{
if(lcs[i+1][j+1]==lcs[i][j])
{
i--;
j--;
}
else
{
if(lcs[i][j+1]>=lcs[i+1][j])
i--;
else
j--;
}
}
}
reverse(ss);//逆转ss
return ss;
}
