LCS最长子序列匹配算法讲解.doc

首先将要看到如何运用动态编程查找两个 DNA 序列的最长公共子序列（longest common subsequence，LCS）。发现了新的基因序列的生物学家通常想知道该基因序列与其他哪个序列最相似。查找 LCS 是计算两个序列相似程度的一种方法：LCS 越长，两个序列越相似。 子序列中的字符与子字符串中的字符不同，它们不需要是连续的。例如，ACE 是 ABCDE 的子序列，但不是它的子字符串。请看下面两个 DNA 序列： S1 = DEGCCCTAGCGDE S2 = DEGCGCAATGDE 这两个序列的 LCS 是 GCCAG。（请注意，这仅是一个 LCS，而不是唯一的 LCS，因为可能存在其他长度相同的公共子序列。这种最优化问题和其他最优化问题的解可能不止一个。） LCS 算法 首先，考虑如何递归地计算 LCS。令： C1 是 S1 最右侧的字符 C2 是 S2 最右侧的字符 S1 是 S1 中 “切掉” C1 的部分 S2 是 S2 中 “切掉” C2 的部分 有三个递归子问题： L1 = LCS(S1, S2) L2 = LCS(S1, S2) L3 = LCS(S1, S2) 结果表明（而且很容易使人相信）原始问题的解就是下面三个子序列中最长的一个： L1 L2 如果 C1 等于 C2，则为 L3 后端加上 C1 ，如果 C1 不等于 C2，则为 L3。 （基线条件（base case）是 S1 或 S2 为长度为零的字符串的情形。在这种情况下，S1 和 S2 的 LCS 显然是长度为零的字符串。） 但是，就像计算斐波纳契数的递归过程一样，这个递归解需要多次计算相同的子问题。可以证明，这种递归解法需要耗费指数级的时间。相比之下，这一问题的动态编程解法的运行时间是 Θ(mn)，其中 m 和 n 分别是两个序列的长度。 为了用动态编程有效地计算 LCS，首先需要构建一个表格，用它保存部分结果。沿着顶部列出一个序列，再沿着左侧从上到下列出另一个序列，如图 2 所示： 图 2. 初始 LCS 表格 这种方法的思路是：将从上向下、从左到右填充表格，每个单元格包含一个数字，代表该行和该列之前的两个字符串的 LCS 的长度。也就是说，每个单元格包含原始问题的一个字问题的解。例如，请看第 6 行第 7 列的单元格：它在 GCGCAATG 序列第二个 C 的右侧，在 GCCCTAGCG 的 T 的下面。这个单元格最终包含的数字就是 GCGC 和 GCCCT 的 LCS 的长度。 首先看一下表格的第二行中应该是什么条目。这一行的单元格保存的 LCS 长度对应的是序列 GCGCAATA 的零长前端和序列 GCCCTAGCG 的 LCS。显然，这些 LCS 的值都是 0。类似的，沿着第二列向下的值也都是 0，这与递归解的基线条件对应。现在表格如图 3 所示： 图 3.填充了基线条件的 LCS 表格 接下来，要实现与递归 HYPERLINK /base/31 \o 算法与数据结构知识库 \t _blank 算法中递归子问题对应的情景，但这时使用的是表格中已经填充的值。在图 4 中，我已经填充了一半左右的单元格： 图 4.填充了一半的 LCS 表格 在填充单元格时，需要考虑以下条件： 它左侧的单元格 它上面的单元格 它左上侧的单元格 下面三个值分别对应着我在前面列出的三个递归子问题返回的值。 V1 = 左侧单元格的值 V2 = 上面单元格的值 V3 = 左上侧单元格的值 在空单元格中填充下面 3 个数字中的最大值： V1 V2 如果 C1 等于 C2 则为 V3 + 1，如果 C1 不等于 C2，则为 V3 ，其中 C1 是当前单元格上面的字符，C2 是当前单元格左侧的字符 请注意，我在图中还添加了箭头，指向当前单元格值的来源。后面的 “回溯” 一节将用这些箭头建立实际的 LCS（与仅仅发现 LCS 长度相反）。 现在填充图 4 中接下来的空单元格 — 在 GCCCTAGCG 中第三个 C 下面和 GCGCAATG 第二个 C 的右侧的单元格。它上面的值是 2，左侧的值是 3，左上侧的值是 2。这个单元格上面的字符和左侧的字符相等（都是 C），所以必须选择 2、3 和 3（左上侧单元格中的 2 + 1）的最大值。所以，这个单元格的值为 3。绘制一个箭头，从该单元格指向其中的值的源单元格。在这个示例中，新的数值可能来自不止一个单元格，所以可以任选一个：例如左上侧单元格。 作为练习，您可以尝试填充表格的余下部分。如果在关联过程中，一直按照左上侧-上侧-左侧的顺序选择单元格，那么会得到如图 5 所示的表格。（当然，如果在关联过程中做了不同的选择，那么箭头就会不同，但是数字是相同的。） 图 5.填充好的 LCS 表格 回想一