动态编程和基因序列比对 计算机科学助力分子生物学

分子生物学越来越多地将计算机科学算法作为研究工具。本文将介绍生物信息学 — 用计算机解决生物学问题。学习动态编程 的基本原理，这是一种高级的计算技术，您将发现它在许多编程项目中都很有用。

<!----><!----><!---->

基因组数据库保存了海量的原始数据。人类基因本身就有接近 30 亿个 DNA 碱基对。为了查遍所有数据并找到其中有意义的关系，分子生物学家们越来越依赖于高效的计算机科学字符串算法。本文将介绍三个这方面的算法，它们都利用动态编程 技术，这是解决最优化问题的一种高级的算法技术，它从下向上寻找子问题的最优解。本文将使用这些算法的 Java™ 实现，还将学习一个用于处理生物学数据的开源的 Java 框架。

基因和字符串算法

基因资料 — DNA 和 RNA — 链是称为核苷酸 的小单元组成的序列。为了回答某些重要的研究问题，研究人员把基因串看作计算机科学的字符串 — 也就是说，可以忽略基因串的物理和化学性质，而将其想像成字符的序列（虽然严格地讲，它们的化学性质通常被编码为字符串算法的参数，您将在本文看到）。

本文的示例使用 DNA，DNA 由腺嘌呤（A ）、胞嘧啶（C ）、胸腺嘧啶（T ）和鸟嘌呤（G ）组成的核苷酸双螺旋组成。DNA 的双螺旋彼此反向互补。A 和 T 是互补的碱基对，C 和 G 也是互补的碱基对。这意味着一个链中的 A 与另一个链中的 T 组成一对（反之亦然），一个链中的 C 与另一个链中的 G 组成一对（反之亦然）。所以，如果知道一个链中的 A 、C 、T 和 G 的顺序，那么就能推导出另一个链中的顺序。所以，可以将一条 DNA 链想像成由字母 A 、C 、G 和 T 组成的字符串。

动态编程

动态编程 是在序列分析中经常使用的一种算法技术。在可以使用递归，但因为递归重复解决相同的子问题造成效率低下的时候，则可以采用动态编程。例如，请看斐波纳契（Fibonacci）序列：0, 1, 1, 2, 3, 5, 8, 13, ... 第一个和第二个斐波纳契数字分别定义为 0 和 1。第 n 个斐波纳契数字是前两个斐波纳契数字的和。所以，可以用清单 1 中的递归函数计算第 n 个斐波纳契数：

清单 1.计算第 n 个斐波纳契数的递归函数

                

public int fibonacci1(int n) {

   if (n == 0) {

      return 0;

   } else if (n == 1) {

      return 1;

   } else {

      return fibonacci1(n - 1) + fibonacci1(n - 2);

   }

}

但是清单 1 的代码效率不高，因为它重复地解决相同的递归子问题。例如，考虑一下 fibonacci1(5) 的计算，如图 1 所示：

图 1.斐波纳契数的递归计算

从图 1 可以看到，fibonacci1(2) 被计算了 3 次。如果从下往上计算斐波纳契数，而不是从上往下计算，效率会更高，如清单 2 所示：

清单 2. 从下往上计算斐波纳契数

                

public int fibonacci2(int n) {

   int[] table = new int[n + 1];

   for (int i = 0; i < table.length; i++) {

      if (i == 0) {

         table[i] = 0;

      } else if (i == 1) {

         table[i] = 1;

      } else {

         table[i] = table[i - 2] + table[i - 1];

      }

   }



   return table[n];

}

清单 2 将中间结果保存在表格中，这样就能重复使用它们，而不是在抛弃之后再重复计算多次。确实，存储表格的内存使用效率较低，因为一次只使用表格中的两个条目，但是现在暂且将这件事放在一边。我们将在本文中使用相似的表格（不过是二维表格）处理比清单 1 复杂得多的示例。清单 2 中的实现所需的时间比清单 1 短许多，清单 2 的运行时间为 O(n) ，而清单 1 中的递归实现的运行时间是 n 的指数。

这就是动态编程的工作原理。遇到一个能用递归算法从上向下解决的问题，然后用从下向上迭代的方式解决。将中间结果保存在表格中供后续步骤使用；否则，需要重复计算它们 — 这显然是一种效率很低的算法。但是动态编程通常被用于最优化问题（比如本文后面的示例），而不是像斐波纳契数这样的问题。下面的示例是一个字符串算法，与计算生物学中经常使用的算法相似。

最长公共子序列问题

首先将要看到如何运用动态编程查找两个 DNA 序列的最长公共子序列（longest common subsequence，LCS） 。发现了新的基因序列的生物学家通常想知道该基因序列与其他哪个序列最相似。查找 LCS 是计算两个序列相似程度的一种方法：LCS 越长，两个序列越相似。

子序列中的字符与子字符串中的字符不同，它们不需要是连续的。例如，ACE 是 ABCDE 的子序列，但不是它的子字符串。请看下面两个 DNA 序列：

• S1 = GCCCTAGCG
• S2 = GCGCAATG

这两个序列的 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 表格

接下来，要实现与递归算法中递归子问题对应的情景，但这时使用的是表格中已经填充的值。在图 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 表格

回想一下，任何单元格中的数字都是该单元格所在行之上和列之前的字符串的 LCS 长度。所以，表格右下角的数字就是字符串 S1 和 S2 （在本例中是 GCCCTAGCG 和 GCGCAATG）的 LCS 的长度。所以，这两个序列的 LCS 长度是 5。

这是在所有动态编程算法中都要牢记的关键点。表格中的每个单元格都包含单元格所在行上面和所在列左侧序列前端问题的解。

使用回溯方法寻找实际的 LCS

接下来要做的就是寻找实际的 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 显示了整个回溯过程：

图 6.在填满的 LCS 表格上进行回溯

通过这种回溯方法，得到的 LCS 为 GCCAG。

动态编程的 Java 语言实现

下面将使用 Java 语言实现动态编程算法 — 首先实现 LCS 算法，然后实现其他两个进行 序列比对（align） 的算法。在每个示例中，都会比较两个序列，而且将使用二维表格存储子问题的解。我们将定义一个抽象类 DynamicProgramming ，它包含所有算法的公共代码。本文的所有示例代码都可以 下载 获得。

开始之前，需要一个类来表示表格中的单元格，如清单 3 的示：

清单 3. Cell 类（部分清单）

                

public class Cell {

   private Cell prevCell;

   private int score;

   private int row;

   private int col;

}

这些算法的第一步都是初始化表格中的值，有时还需要初始化表格中的指针。设置的初始值和指针因算法不同而不同。正因如此，清单 4 所示的 DynamicProgramming 类定义了两个抽象方法：

清单 4. DynamicProgramming 初始化代码

                

protected void initialize() {



   for (int i = 0; i < scoreTable.length; i++) {

      for (int j = 0; j < scoreTable[i].length; j++) {

         scoreTable[i][j] = new Cell(i, j);

      }

   }

   initializeScores();

   initializePointers();



   isInitialized = true;

}



protected void initializeScores() {

   for (int i = 0; i < scoreTable.length; i++) {

      for (int j = 0; j < scoreTable[i].length; j++) {

         scoreTable[i][j].setScore(getInitialScore(i, j));

      }

   }

}



protected void initializePointers() {

   for (int i = 0; i < scoreTable.length; i++) {

      for (int j = 0; j < scoreTable[i].length; j++) {

         scoreTable[i][j].setPrevCell(getInitialPointer(i, j));

      }

   }

}



protected abstract int getInitialScore(int row, int col);



protected abstract Cell getInitialPointer(int row, int col);

接下来，要用值和指针填充表格的每个单元格。同样，填充方法也视算法不同而不同。所以使用抽象方法 fillInCell(Cell, Cell, Cell, Cell) 。清单 5 显示了 DynamicProgramming 类中用于填充表格的方法：

清单 5. 填充表格的 DynamicProgramming 代码

                

protected void fillIn() {

   for (int row = 1; row < scoreTable.length; row++) {

      for (int col = 1; col < scoreTable[row].length; col++) {

         Cell currentCell = scoreTable[row][col];

         Cell cellAbove = scoreTable[row - 1][col];

         Cell cellToLeft = scoreTable[row][col - 1];

         Cell cellAboveLeft = scoreTable[row - 1][col - 1];

         fillInCell(currentCell, cellAbove, cellToLeft, cellAboveLeft);

      }

   }

}



protected abstract void fillInCell(Cell currentCell, Cell cellAbove,

      Cell cellToLeft, Cell cellAboveLeft);

最后，进行回溯。回溯方法因算法不同而不同。清单 6 显示了 DynamicProgramming.getTraceback() 方法：

清单 6. DynamicProgramming.getTraceback() 方法

                

abstract protected Object getTraceback();

LCS 的 Java 实现

现在可以编写 LCS 算法的 Java 实现了。

初始化单元格中的值的方法很简单：只需将它们的初值全部设为 0（后面还会重设其中的一些单元格），如清单 7 所示：

清单 7. LCS 的初始化方法

                

protected int getInitialScore(int row, int col) {

   return 0;

}

清单 8 显示了向表格中的初始单元格填充值和指针的代码：

清单 8.填充单元格值和指针的 LCS 方法

                

protected void fillInCell(Cell currentCell, Cell cellAbove, Cell cellToLeft,

      Cell cellAboveLeft) {

   int aboveScore = cellAbove.getScore();

   int leftScore = cellToLeft.getScore();

   int matchScore;

   if (sequence1.charAt(currentCell.getCol() - 1) == sequence2

         .charAt(currentCell.getRow() - 1)) {

      matchScore = cellAboveLeft.getScore() + 1;

   } else {

      matchScore = cellAboveLeft.getScore();

   }

   int cellScore;

   Cell cellPointer;

   if (matchScore >= aboveScore) {

      if (matchScore >= leftScore) {

         // matchScore >= aboveScore and matchScore >= leftScore

         cellScore = matchScore;

         cellPointer = cellAboveLeft;

      } else {

         // leftScore > matchScore >= aboveScore

         cellScore = leftScore;

         cellPointer = cellToLeft;

      }

   } else {

      if (aboveScore >= leftScore) {

         // aboveScore > matchScore and aboveScore >= leftScore

         cellScore = aboveScore;

         cellPointer = cellAbove;

      } else {

         // leftScore > aboveScore > matchScore

         cellScore = leftScore;

         cellPointer = cellToLeft;

      }

   }

   currentCell.setScore(cellScore);

   currentCell.setPrevCell(cellPointer);

}

最后，用回溯的方式构建实际的 LCS：

清单 9. LCS 回溯代码

                

protected Object getTraceback() {

   StringBuffer lCSBuf = new StringBuffer();

   Cell currentCell = scoreTable[scoreTable.length - 1][scoreTable[0].length - 1];

   while (currentCell.getScore() > 0) {

      Cell prevCell = currentCell.getPrevCell();

      if ((currentCell.getRow() - prevCell.getRow() == 1 && currentCell

            .getCol()

            - prevCell.getCol() == 1)

            && currentCell.getScore() == prevCell.getScore() + 1) {

         lCSBuf.insert(0, sequence1.charAt(currentCell.getCol() - 1));

      }

      currentCell = prevCell;

   }



   return lCSBuf.toString();

}

很容易发现，该算法需要花费 Θ(mn) 的时间（和空间）来进行运算，其中 m 和 n 是两个序列的长度。填充每个单元格需要的时间是恒定的 — 只需有限数量的相加和比较 — 必须填充 mn 个单元格。而且，回溯算法的用时为 O(m + n) 。

下面两个 Java 示例用于实现序列的比对算法： Needleman-Wunsch 和 Smith-Waterman。

序列比对

同源性（Homology） 同源性 是个重要的生物学概念。如果两个物种具有共同的进化祖先，就认为它们是同源的。同源的物种的 DNA 中的许多部分是相同的。反过来，如果两个物种有相似的 DNA 子串，那么通常可以断定它们来自共同的祖先。序列比对算法可用来查找这种类似的 DNA 子串。

基因学的一个主要主题就是比较 DNA 序列并尝试找出两个序列的公共部分。如果两个 DNA 序列有类似的公共子序列，那么这些两个序列很可能是同源的 （请参阅 “同源性 ” 侧栏）。在比对两个序列时，不仅要考虑完全匹配的字符，还要考虑一个序列中的空格或间隙（或者，相反地，要考虑另一个序列中的插入部分）和不匹配，这两个方面都可能意味着突变。在序列比对中，需要找到最优的比对（最优比对大致是指要将匹配的数量最大化，将空格和不匹配的数量最小化）。如果要更正式些，您可以确定一个分数，为匹配的字符添加分数、为空格和不匹配的字符减去分数。

全局和局部序列比对

全局序列比对 尝试找到两个完整的序列 S1 和 S2 之间的最佳比对。以下面两个 DNA 序列为例：

• S1 = GCCCTAGCG
• S2 = GCGCAATG

如果为每个匹配字符一分，一个空格扣两分，一个不匹配字符扣一分，那么下面的比对就是全局最优比对：

• S1' = GCCCTAGCG
• S2' = GCGC-AATG

连字符（-）代表空格。在 S2' 中有五个匹配字符，一个空格（或者反过来说，在 S1' 中有一个插入项），有三个不匹配字符。这样得到的分数是 (5 * 1) + (1 * -2) + (3 * -1) = 0，这是能够实现的最佳结果。

使用局部序列比对 ，不必对两个完整的序列进行比对；可以在每个序列中使用某些部分来获得最大得分。使用同样的序列 S1 和 S2 ，以及同样的得分方案，可以得到以下局部最优比对 S1'' 和 S2'' ：

• S1      = GCCCTAGCG
• S1'' =            GCG
• S2'' =            GCG
• S2      =            GCGCAATG

虽然这个局部比对恰好没有不匹配字符或空格，但是一般情况下，局部比对可能存在不匹配字符或空格。这个局部比对的得分是 (3 * 1) + (0 * -2) + (0 * -1) = 3。（最佳局部比对的得分要大于或等于最佳全局比对的得分，这是因为全局比对 也属于 局部比对。）

Needleman-Wunsch 算法

Needleman-Wunsch 算法 用来计算全局比对。它的思路与 LCS 算法相似。这个算法也使用二维表格，一个序列沿顶部展开，一个序列沿左侧展开。而且也能通过以下三个途径到达每个单元格：

• 来自上面的单元格，代表将左侧的字符与空格比对。
• 来自左侧的单元格，代表将上面的字符与空格比对。
• 来自左上侧的单元格，代表与左侧和上面的字符比对（可能匹配也可能不匹配）

我首先给出完整的表格（参见图 7），在解释如何填充表格的时候可以返回来查看它：

图 7.带有回溯的填充好的 Needleman-Wunsch 表格

首先，必须初始化表格。这意味着填充第二行和第二列的分数和指针。填充第二行的操作意味着使用位于顶部的第一个序列中的字符，并使用空格，而不是使用左侧从上到下的序列中的第一个字符。空格的扣分是 -2，所以每次使用空格的时候，就给以前的单元格加了 -2 分。以前的单元格是左侧的单元格。这就说明了在第二行中为什么得到了 0, -2, -4, -6, ... 这样的序列。用相似的方式得到第二列的得分和指针。清单 10 展示了 Needleman-Wunsch 算法的初始化代码：

清单 10. Needleman-Wunsch 的初始化代码

                

protected Cell getInitialPointer(int row, int col) {

   if (row == 0 && col != 0) {

      return scoreTable[row][col - 1];

   } else if (col == 0 && row != 0) {

      return scoreTable[row - 1][col];

   } else {

      return null;

   }

}



protected int getInitialScore(int row, int col) {

   if (row == 0 && col != 0) {

      return col * space;

   } else if (col == 0 && row != 0) {

      return row * space;

   } else {

      return 0;

   }

}

接下来，需要填充余下的单元格。同 LCS 算法一样，对于每个单元格，都有三个选择，要从中选择最大的。可以从上面、左侧、左上侧到达每个单元格。假设 S1 和 S2 是要比对的字符串，S1' 和 S2' 是生成的比对中的字符串。从上面到达单元格相当于将左面的字符从 S2 加入 S2' ，跳过上面的 S1 中的当前字符，并在 S1' 中加入一个空格。因为一个空格的分数是 -2，所以当前单元格的得分要从上面的单元格得分减 2 得到。类似的，将左边的单元格得分减 2，可以从左侧到达空单元格。最后，可以将上面的字符加入到 S1' 中，将左边的字符加入到 S2' 中。这就相当于从左上侧进入空白单元格。这两个字符将会匹配，在这种情况下，新的得分就是左上侧单元格的得分减 1。在这三种可能性当中，选择得分最大的一个（如果得分相等，可以从得分高的单元格中从任选一个）。观察 图 7 中的指针，能够找到这三种可能性的示例。

清单 11 显示了填充空白单元格的代码：

清单 11.填充表格的 Needleman-Wunsch 代码

                

protected void fillInCell(Cell currentCell, Cell cellAbove, Cell cellToLeft,

      Cell cellAboveLeft) {

   int rowSpaceScore = cellAbove.getScore() + space;

   int colSpaceScore = cellToLeft.getScore() + space;

   int matchOrMismatchScore = cellAboveLeft.getScore();

   if (sequence2.charAt(currentCell.getRow() - 1) == sequence1

         .charAt(currentCell.getCol() - 1)) {

      matchOrMismatchScore += match;

   } else {

      matchOrMismatchScore += mismatch;

   }

   if (rowSpaceScore >= colSpaceScore) {

      if (matchOrMismatchScore >= rowSpaceScore) {

         currentCell.setScore(matchOrMismatchScore);

         currentCell.setPrevCell(cellAboveLeft);

      } else {

         currentCell.setScore(rowSpaceScore);

         currentCell.setPrevCell(cellAbove);

      }

   } else {

      if (matchOrMismatchScore >= colSpaceScore) {

         currentCell.setScore(matchOrMismatchScore);

         currentCell.setPrevCell(cellAboveLeft);

      } else {

         currentCell.setScore(colSpaceScore);

         currentCell.setPrevCell(cellToLeft);

      }

   }

}

接下来，需要得到实际的比对字符串 — S1' 和 S2' — 以及比对的得分。右下角单元格中的得分包含 S1 和 S2 的最大比对得分，就像在 LCS 算法中包含 LCS 的长度一样。而且，与 LCS 算法类似，要获得 S1' 和 S2' ，要从右下角单元格开始沿着指针回溯，反向构建 S1' 和 S2' 。从表格的构建过程可知，从上向下对应着将左侧字符从 S2 加入到 S2' 中，将空格加入 S1' 中；从左向右对应着将上面的字符从 S1 加入到 S1' 中，将空格加入 S2' 中；而向下和向右移动意味着分别将来自 S1 和 S2 的字符加入 S1' 和 S2' 。

Needleman-Wunsch 中使用的回溯代码与 Smith-Waterman 中局部比对的回溯代码基本相同，区别只是开始的单元格以及如何知道何时结束回溯。清单 12 显示了两个算法共享的代码：

清单 12. Needleman-Wunsch 和 Smith-Waterman 共同的回溯代码

                

protected Object getTraceback() {

   StringBuffer align1Buf = new StringBuffer();

   StringBuffer align2Buf = new StringBuffer();

   Cell currentCell = getTracebackStartingCell();

   while (traceBackIsNotDone(currentCell)) {

      if (currentCell.getRow() - currentCell.getPrevCell().getRow() == 1) {

         align2Buf.insert(0, sequence2.charAt(currentCell.getRow() - 1));

      } else {

         align2Buf.insert(0, '-');

      }

      if (currentCell.getCol() - currentCell.getPrevCell().getCol() == 1) {

         align1Buf.insert(0, sequence1.charAt(currentCell.getCol() - 1));

      } else {

         align1Buf.insert(0, '-');

      }

      currentCell = currentCell.getPrevCell();

   }



   String[] alignments = new String[] { align1Buf.toString(),

         align2Buf.toString() };



   return alignments;

}



protected abstract boolean traceBackIsNotDone(Cell currentCell);



protected abstract Cell getTracebackStartingCell();

清单 13 显示了特定于 Needleman-Wunsch 算法的回溯代码:

清单 13. Needleman-Wunsch 算法的回溯代码

                

protected boolean traceBackIsNotDone(Cell currentCell) {

   return currentCell.getPrevCell() != null;

}



protected Cell getTracebackStartingCell() {

   return scoreTable[scoreTable.length - 1][scoreTable[0].length - 1];

}

原始的 Needleman-Wunsch 算法 严格地讲，我还没有演示 Needleman-Wunsch 算法。Needleman-Wunsch 发表的原始算法运行时间达到三个数量级，而且已经不再使用了。但是，这里讨论的二次数量级的算法通常依然称为 Needleman-Wunsch 算法。

通过回溯能够得到本节开始时提到的最优全局比对：

• S1' = GCCCTAGCG
• S2' = GCGC-AATG

显然，这个算法的运行时间为 O(mn) 。

Smith-Waterman 算法

在 Smith-Waterman 算法中，不必比对整个序列。两个零长字符串即为得分为 0 的局部比对，这一事实表明在构建局部比对时，不需要使用负分。这样会造成进一步比对所得到的分数低于通过 “重设” 两个零长字符串所能得到的分数。而且，局部比对不需要到达任何一个序列的末端，所以也不需要从右下角开始回溯：可以从得分最高的单元格开始回溯。

这导致 Smith-Waterman 算法与 Needleman-Wunsch 算法存在着三个区别。首先，在初始化阶段，第一行和第一列全填充为 0（而且第一行和第一列的指针均为空）。清单 14 显示了 Smith-Waterman 算法的初始化代码：

清单 14. Smith-Waterman 算法的初始化

                

protected int getInitialScore(int row, int col) {

   return 0;

}



protected Cell getInitialPointer(int row, int col) {

   return null;

}

第二，在填充表格时，如果某个得分为负，那么就用 0 代替，只对得分为正的单元格添加返回指针。请注意在清单 15 中仍然在跟踪哪个单元格的得分高，在回溯时要使用这个记录：

清单 15.填充单元格的 Smith-Waterman 代码

                

protected void fillInCell(Cell currentCell, Cell cellAbove, Cell cellToLeft,

      Cell cellAboveLeft) {

   int rowSpaceScore = cellAbove.getScore() + space;

   int colSpaceScore = cellToLeft.getScore() + space;

   int matchOrMismatchScore = cellAboveLeft.getScore();

   if (sequence2.charAt(currentCell.getRow() - 1) == sequence1

         .charAt(currentCell.getCol() - 1)) {

      matchOrMismatchScore += match;

   } else {

      matchOrMismatchScore += mismatch;

   }

   if (rowSpaceScore >= colSpaceScore) {

      if (matchOrMismatchScore >= rowSpaceScore) {

         if (matchOrMismatchScore > 0) {

            currentCell.setScore(matchOrMismatchScore);

            currentCell.setPrevCell(cellAboveLeft);

         }

      } else {

         if (rowSpaceScore > 0) {

            currentCell.setScore(rowSpaceScore);

            currentCell.setPrevCell(cellAbove);

         }

      }

   } else {

      if (matchOrMismatchScore >= colSpaceScore) {

         if (matchOrMismatchScore > 0) {

            currentCell.setScore(matchOrMismatchScore);

            currentCell.setPrevCell(cellAboveLeft);

         }

      } else {

         if (colSpaceScore > 0) {

            currentCell.setScore(colSpaceScore);

            currentCell.setPrevCell(cellToLeft);

         }

      }

   }

   if (currentCell.getScore() > highScoreCell.getScore()) {

      highScoreCell = currentCell;

   }

}

最后，在回溯的时候，从得分最高的单元格开始，回溯到得分为 0 的单元格为止。除此之外，回溯的方式与 Needleman-Wunsch 算法完全相同。清单 16 显示了 Smith-Waterman 算法的回溯代码：

清单 16. Smith-Waterman 算法的回溯代码

                

protected boolean traceBackIsNotDone(Cell currentCell) {

   return currentCell.getScore() != 0;

}



protected Cell getTracebackStartingCell() {

   return highScoreCell;

}

图 8 演示了在本文使用的 S1 和 S2 序列上运行 Smith-Waterman 算法的情况：

图 8.带有回溯的填充好的 Smith-Waterman 表格

同 Needleman-Wunsch 算法一样，运行 Smith-Waterman 算法得到（或者查阅图 8 得到的）的最优局部比对是：

• S1      = GCCCTAGCG
• S1'' =            GCG
• S2'' =            GCG
• S2      =            GCGCAATG

ALIGN、FASTA 和 BLAST 本文介绍了 Needleman-Wunsch 和 Smith-Waterman 算法的基本实现，未进行优化，查找全局和局部比对的时间为 O(mn) 。实际的研究人员通常不是比较两个序列，而是试图找出与特定序列相似的所有序列。如果发现某个相似序列具有已知的生物学功能，那么原来的序列也很有可能具有相似的功能，因为相似的序列通常会有相似的功能。ALIGN、FASTA 和 BLAST（Basic Local Alignment Search Tool）是用于查找全局（ALIGN）和局部（FASTA 和 BLAST）比对的专业应用程序。BLAST 在大型序列数据库中搜索与用户输入的序列相似（和可能同源的）的序列，并根据相似程度给出结果的等级（请参阅 参考资料 ）。BLAST 不直接使用 Smith-Waterman 算法，因为即使使用具有二次数量级的运行时间，将输入序列与超大型基因序列数据库的每个序列（每个序列可能还包含 30 亿或者更多个碱基对）进行比对还是太慢。BLAST 首先通过一个称为 播种(seeding） 的流程寻找 种子 ，种子就是可能的匹配的开始位置。然后 BLAST 用动态编程算法将找到的可能匹配扩展为与输入序列匹配的实际的局部比对。最后再确定哪些匹配在统计意义上是重要的，并对其进行评级。这种带有启发性的过程没有 Smith-Waterman 算法那么灵敏 （精确），但要快得多。

BioJava

BioJava 是一个开源项目，目的是开发一个处理生物学数据的 Java 框架。它的功能包括：操纵生物学序列的对象，进行序列分析的 GUI 工具，以及包含动态编程工具包的分析和统计例程（请参阅 参考资料 ）。

清单 17 演示了如何根据本章前面的示例使用的序列和得分方案，运行 Needleman-Wunsch 和 Smith-Waterman 算法的 BioJava 实现：

清单 17. BioJava 代码中的序列比对代码（基于 Andreas Dräger 的 BioJava 示例代码）

                

// The alphabet of the sequences. For this example DNA is chosen.

FiniteAlphabet alphabet = (FiniteAlphabet) AlphabetManager

      .alphabetForName("DNA");

// Use a substitution matrix with equal scores for every match and every

// replace.

int match = 1;

int replace = -1;

SubstitutionMatrix matrix = new SubstitutionMatrix(alphabet, match,

      replace);

// Firstly, define the expenses (penalties) for every single operation.

int insert = 2;

int delete = 2;

int gapExtend = 2;

// Global alignment.

SequenceAlignment aligner = new NeedlemanWunsch(match, replace, insert,

      delete, gapExtend, matrix);

Sequence query = DNATools.createDNASequence("GCCCTAGCG", "query");

Sequence target = DNATools.createDNASequence("GCGCAATG", "target");

// Perform an alignment and save the results.

aligner.pairwiseAlignment(query, // first sequence

      target // second one

      );

// Print the alignment to the screen

System.out.println("Global alignment with Needleman-Wunsch:\n"

      + aligner.getAlignmentString());



// Perform a local alignment from the sequences with Smith-Waterman.

aligner = new SmithWaterman(match, replace, insert, delete, gapExtend,

      matrix);

// Perform the local alignment.

aligner.pairwiseAlignment(query, target);

System.out.println("\nLocal alignment with Smith-Waterman:\n"

      + aligner.getAlignmentString());

BioJava 方法有一个共同之处。首先，请注意 SubstitutionMatrix 的用法。目前为止的示例都假定 DNA 碱基对之间不匹配的扣分应该相等 — 例如，认为 G 变异为 A 与变异为 C 的可能性相当。但在真正的生物学序列中并不是这样，尤其是在蛋白质的氨基酸中。少见的不匹配的扣分要比常见不匹配的扣分多。使用置换矩阵能够根据每对符号单独分配匹配得分。在某种意义上，替换矩阵是对化学属性的编码。例如，BLAST 搜索中经常使用蛋白质的 BLOSUM（BLOcks SUbstitution Matrix）矩阵；BLOSUM 矩阵的值就是根据实际经验确定的。

然后，请注意，insert 和 delete 分数的使用，而不仅仅使用空格得分。就像我说的，可以将空格想象成在没有空格的序列中执行了一次插入，或在有空格的序列中执行了一次删除。一般来说，比较两个序列时还有两种互补的方法。我们刚才一直用 “静态” 的方式考察序列以及序列间的差异。除此之外，还可以查找最少数量的插入、删除和对单独的符号所作的更改来比较序列，从而将一个序列转换成另一个序列。最小数量的更改称为编辑距离 。在计算编辑距离的时候，可以给插入和删除分配不同的得分。例如，可能插入更常见，所以插入扣的分要比删除少。

Perl 在生物信息学中的应用 许多大型服务器生物信息学软件是用 C 或 C++ 编写的。BLAST 最初是用 C 编写的，现在有一个 C++ 版本。但是研究人员 — 在许多情况下，专业生物学家才是他们的第一职业，其次才是程序员 — 编写的许多小应用程序都是用 Perl 编写的。这可能是因为最大的开源生物信息学库 Bioperl 是用 Perl 编写的。如果想进行生物信息学编程，那么可能需要学一些 Perl 和 Bioperl 知识。

现在请注意 gapExtend 变量。从技术上讲，间隙是最大的连续空格序列。但是，某些文献使用 间隙 这个术语时，实际指的是空格。在本文的示例中，所有空格的得分都是相等的，即使这些空格处在更大的间隙内。但是，在实际中，一旦一个间隙开始，那么它包含多个空格的机率就比只包含一个空格的机会要大。所以，为了得到有意义的结果，对于间隙中的后续空格的扣分应该比初始空格的扣分少。这就是 gapExtend 变量的用途。请记住，从算法上讲，所有这些得分方案存在一些主观因素，但是您显然希望计算的字符串编辑距离尽可能地符合自然界的进化距离。（针对不同情况使用不同的得分方案本身就是一个相当有趣和复杂的子领域。）

最后，insert 、delete 和 gapExtend 变量的值都为正，而不像前面那样为负，这是因为它们被定义为开支（成本或扣除）。

在运行 清单 17 中的代码时，得到以下输出：

清单 18. 输出

                

Global alignment with Needleman-Wunsch:



  Length:9

  Score:-0.0

  Query:query,Length:9

  Target:target,Length:8



Query:  1 gccctagcg 9

          || | |  |

Target: 1 gcgc-aatg 8





Local alignment with Smith-Waterman:



  Length:3

  Score:3.0

  Query:query,Length:9

  Target:target,Length:8



Query:  7 gcg 9

          |||

Target: 1 gcg 3

对于局部和全局比对，得到的得分与前面的得分相同。Smith-Waterman 算法的这个实现给出的局部比对与前面得到的结果相同。Needleman-Wunsch 算法的这个实现提供的全局比对与前面得到的结果不同，但是得分相同。但是，它们都是最优全局比对。回想一下，在填充表格的时候，单元格的最大值有些时候可能来自多个单元格。如果绘制的箭头指向不同的单元格，最终就会形成不同的比对（但是得分都相同）。

结束语

分享这篇文章……

 

	提交到 Digg
	发布到 del.icio.us
	Slashdot 一下！

本文介绍了使用动态编程能够解决的三个问题示例。这些问题都有共同的特征：

• 每个问题的解都能用递归关系表示。
• 用递归方法对这些递归关系的直接实现会造成解决方案效率低下，因为其中包含了对子问题的多次计算。
• 一个问题的最优解能够用原始问题的子问题的最优解构造得到。

动态编程还可用于矩阵链相乘、装配线规划和计算机象棋程序。解决编程竞赛的困难问题时经常需要使用动态编程。对动态编程感兴趣的读者可以参考图书 Introduction to Algorithms （请参阅 参考资料 ），了解关于使用动态编程的时机，以及通常如何证明动态编程算法的正确性的更多细节。

动态编程可能是计算机科学在生物学上最重要的应用，但肯定不是唯一的应用。生物信息学和计算生物学是交叉学科领域，正在迅速成为具有专门的学术程序的独立学科。许多分子生物学家现在都掌握了一些编程技术，对于能够了解一些生物学理论的程序员来说，还有许多非常有趣和重要的工作等着他们去做。如果想了解更多内容，请参阅 参考资料 ，获得可能有用的资料。