基因組數據庫保存了海量的原始數據。人類基因本身就有接近 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個斐波納契數的遞歸函數
|
1
2
3
4
5
6
7
8
9
|
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. 從下往上計算斐波納契數
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
|
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類（部分清單）
|
1
2
3
4
5
6
|
public class Cell {
private Cell prevCell;
private int score;
private int row;
private int col;
}
|
這些算法的第一步都是初始化表格中的值，有時還需要初始化表格中的指針。設置的初始值和指針因算法不同而不同。正因如此，清單 4 所示的 DynamicProgramming類定義了兩個抽象方法：
清單 4. DynamicProgramming初始化代碼
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
|
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代碼
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
|
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()方法
|
1
|
abstract protected Object getTraceback();
|
LCS 的 Java 實現
現在可以編寫 LCS 算法的 Java 實現了。
初始化單元格中的值的方法很簡單：只需將它們的初值全部設為 0（后面還會重設其中的一些單元格），如清單 7 所示：
清單 7. LCS 的初始化方法
|
1
2
3
|
protected int getInitialScore(int row, int col) {
return 0;
}
|
清單 8 顯示了向表格中的初始單元格填充值和指針的代碼：
清單 8. 填充單元格值和指針的 LCS 方法
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
|
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 回溯代碼
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
|
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 = GCCCTAGCGS1''= GCGS2''= GCGS2 = 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 的初始化代碼
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
|
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 代碼
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
|
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 共同的回溯代碼
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
|
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 算法的回溯代碼
|
1
2
3
4
5
6
7
|
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 算法的初始化
|
1
2
3
4
5
6
7
|
protected int getInitialScore(int row, int col) {
return 0;
}
protected Cell getInitialPointer(int row, int col) {
return null;
}
|
第二，在填充表格時，如果某個得分為負，那么就用 0 代替，只對得分為正的單元格添加返回指針。請注意在清單 15 中仍然在跟蹤哪個單元格的得分高，在回溯時要使用這個記錄：
清單 15. 填充單元格的 Smith-Waterman 代碼
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
|
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 算法的回溯代碼
|
1
2
3
4
5
6
7
|
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 = GCCCTAGCGS1''= GCGS2''= GCGS2 = GCGCAATG
BioJava
BioJava 是一個開源項目，目的是開發一個處理生物學數據的 Java 框架。它的功能包括：操縱生物學序列的對象，進行序列分析的 GUI 工具，以及包含動態編程工具包的分析和統計例程（請參閱 參考資料）。
清單 17 演示了如何根據本章前面的示例使用的序列和得分方案，運行 Needleman-Wunsch 和 Smith-Waterman 算法的 BioJava 實現：
清單 17. BioJava 代碼中的序列比對代碼（基於 Andreas Dräger 的 BioJava 示例代碼）
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
|
// 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. 輸出
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
|
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)
這算法是用來比較兩個字符串差異的。LD值越大，兩個字符串差異越大，LD值越小，兩個字符串差異越小。