动态编程和基因序列比对

基因组数据库保存了海量的原始数据。人类基因本身就有接近 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。

回页首

序列比对

基因学的一个主要主题就是比较 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]; 
 }

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

• 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

回页首

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

现在请注意 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 算法的这个实现提供的全局比对与前面得到的结果不同，但是得分相同。但是，它们都是最优全局比对。回想一下，在填充表格的时候，单元格的最大值有些时候可能来自多个单元格。如果绘制的箭头指向不同的单元格，最终就会形成不同的比对（但是得分都相同）。

回页首

结束语

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

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

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

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

回页首

致谢

感谢 Sonna Bristle 让我对计算生物学产生了兴趣，感谢 IBM 的 Carlos P. Sosa 审阅了本文并给出了宝贵的建议。

回页首

下载

描述	名字	大小
本文的示例代码	j-seqalign-code.zip	12KB