字符串與模式匹配算法（六）：Needleman–Wunsch算法

一、Needleman-Wunsch 算法

　　尼德曼-翁施算法（英語：Needleman-Wunsch Algorithm）是基於生物信息學的知識來匹配蛋白序列或者DNA序列的算法。這是將動態算法應用於生物序列的比較的最早期的幾個實例之一。該算法是由 Saul B. Needlman和 Christian D. Wunsch 兩位科學家於1970年發明的。本算法高效地解決了如何將一個龐大的數學問題分解為一系列小問題，並且從一系列小問題的解決方法重建大問題的解決方法的過程。該算法也被稱為優化匹配算法和整體序列比較法。時至今日尼德曼-翁施算法仍然被廣泛應用於優化整體序列比較中。

二、初始化得分矩陣

　　首先建立如下的得分矩陣。從第一列第一行的位置起始。計算過程從d _0,0開始，可以是按行計算，每行從左到右，也可以是按列計算，每列從上到下。當然，任何計算過程，只要滿足在計算 d_{(i , j)} 時 d _{(i-1, j)}（上邊）、d _{(i-1 , j-1)} （左上）和 d_{(i, j-1)} （左邊）都已經被計算這個條件即可。在計算 d_{(i , j)}后，需要保存d_{(i , j)}是從 d _{(i-1 , j)} 、d_{(i-1 , j-1)} 或 d_{(i, j-1)} 中的一個推進的，或保存計算的路徑，以便於后續處理。上述計算過程到 d_{(m , n)} 結束，其中m，n各位兩個序列的長度。

		G	C	A	T	G	C	U
	0
G
A
T
T
A
C
A

三、選擇得分體系

　　我們可以看出每個位置配對都有三種可能情況：匹配， 不匹配與錯位（或插入）：

• 匹配： 兩個字母相同

• 不匹配：兩個字母不相同

• 錯位：一個字母與另一個序列中的間隔（空白）相匹配

　　這三種得分情況有很多打分標准，這些情況都總結在得分體系的部分。從現在開始，我們將使用Needleman 和Wunsch創造的簡單體系來進行打分，即匹配得1分，不匹配得-1分，錯位得-1分.

四、填充評分矩陣

　　開始於第二行中的第二列，初始值為0。通過矩陣一行一行移動，計算每個位置的分數。得分被計算為從現有的分數可能的最佳得分（即最高）的左側，頂部或左上方（對角線）。當一個得分從頂部計算，或從左邊這代表在我們的對准的插入缺失。當我們從對角線計算分數這表示兩個字母所得位置匹配的對准。定不存在“向上”或“左上”的位置對第二行，我們只能從現有單元向左添加。因此，我們添加-1的權利，因為這代表了從以前的比分插入缺失。這導致在第一行是0，-1，-2，-3，-4，-5，-6，-7。這同樣適用於第二列，因為我們只具有以上現有分數。因此，我們有：

		G	C	A	T	G	C	U
	0	-1	-2	-3	-4	-5	-6	-7
G	-1
A	-2
T	-3
T	-4
A	-5
C	-6
A	-7

　　第一種情況是存在向三個方向構築矩陣，周圍位置得分如下：

		G
	0	-1
G	-1	?

　　該單元格具有三個可能的候選總和：

• 對角線的左上鄰居的分數為0。G和G的配對是匹配項，因此添加匹配項的分數：0 + 1 = 1
• 正上鄰居的得分為-1，從那里移動代表一個indel，因此將indel的得分相加：（-1）+（-1）=（-2）
• 左鄰居的得分也為-1，代表插入缺失，也產生（-2）。

　　得分最高的單元格也必須記錄下來。在的文章最開始說的那矩陣圖中，它表示為從行和列3中的單元格到行和列2中的單元格的箭頭，?號處根據評分規則最后為1。

　　對接下來的位置，我們有不同的選擇。

		G	C
	0	-1	-2
G	-1	1	X
A	-2	Y

　　X:

• 上: (-2)+(-1) = (-3)

• 左: (+1)+(-1) = (0)

• 左上: (-1)+(-1) = (-2)

　　Y:

• 上: (1)+(-1) = (0)

• 左: (-2)+(-1) = (-3)

• 左上: (-1)+(-1) = (-2)

　　對於X和Y，最高得分均為0。
　　兩個或所有相鄰小格可能會達到最高的候選分數：

	T	G
T	1	1
A	0	X

• Top: (1)+(-1) = (0)

• Top-Left: (1)+(-1) = (0)

• Left: (0)+(-1) = (-1)

　　在這種情況下，必須將達到最高候選得分的所有方向記錄為中並完成圖中可能的起點，例如，在第7行和第7列的單元格中。

　　以這種方式填寫表格會給出分數或所有可能的比對候選，右下角單元格中的分數代表最佳比對的比對分數。

五、追溯本源

　　按照箭頭的方向標記從右下角的單元格到左上角的單元格的路徑。從此路徑開始，將按照以下規則構建序列：

• 對角箭頭表示匹配或不匹配，因此原始單元格的列字母和行字母將對齊。

• 水平或垂直箭頭表示插入/缺失。 水平箭頭將使間隙（“-”）與行的字母（“側邊”序列）對齊，垂直箭頭將使間隙與列的字母（“頂部”序列）對齊。

• 如果有多個箭頭可供選擇，則它們表示路線的分支。 如果兩個或多個分支都屬於從右下角到左上角單元格的路徑，則它們是同等可行的路線。 在這種情況下，請注意將路徑作為單獨的對齊候選。

　　遵循這些規則，上圖中的一個可能的對齊候選者的步驟為：

U → CU → GCU → -GCU → T-GCU → AT-GCU → CAT-GCU → GCAT-GCU
A → CA → ACA → TACA → TTACA → ATTACA → -ATTACA → G-ATTACA
             ↓
    (branch) → TGCU → ...
             → TACA → ...

六、評分系統

　　最簡單的計分方案只是為每個匹配，不匹配和插入缺失給出一個值。 上面的分步指南使用match = 1，mismatch = -1，indel = -1。 因此，對准分數越低，編輯距離越大，因為該計分系統需要較高的分數。 另一個計分系統可能是：

• Match = 0

• Indel = 1

• Mismatch = 1

　　對於此系統，對齊分數將代表兩個字符串之間的編輯距離。 可以針對不同的情況設計不同的評分系統，例如，如果認為差距對您的對齊方式非常不利，則可以使用對差距進行嚴重懲罰的評分系統，例如：

• Match = 0

• Mismatch = 1

• Indel = 10

七、Needleman–Wunsch算法代碼

　　計算評分矩陣，按行推進，根據左、上、左上計算的結果取最大的：

    /**
     * 一種用於計算全局最大相似度矩陣的實用方法。
     *
     * @param sequence1 代表第一個氨基酸序列的字符串。
     * @param sequence1 代表第二個氨基酸序列的字符串。
     * @param matchScore 代表分配給比賽的分數的浮點數。
     * @param mismatchScore 一個浮點數，代表分配給不匹配項的分數。
     * @param indelScore 浮點數表示分配給插入和刪除的分數。
     */
    public static float[][] computeMatrix(String sequence1, String sequence2, float matchScore, float mismatchScore, float indelScore) {
        sequence1 = "-" + sequence1;
        sequence2 = "-" + sequence2;

        float[][] resultMatrix = new float[sequence1.length()][sequence2.length()];

        for (int i = 0; i < sequence1.length(); i++) {
            resultMatrix[i][0] = i * indelScore;
        }
        for (int j = 0; j < sequence2.length(); j++) {
            resultMatrix[0][j] = resultMatrix[0][j] = j * indelScore;
        }

        for (int i = 1; i < sequence1.length(); i++) {
            for (int j = 1; j < sequence2.length(); j++) {
                resultMatrix[i][j] = Math.max(resultMatrix[i - 1][j] + indelScore,
                        Math.max(resultMatrix[i][j - 1] + indelScore, resultMatrix[i - 1][j - 1] +
                                (sequence1.charAt(i) == sequence2.charAt(j) ? matchScore : mismatchScore)));
            }
        }

        return resultMatrix;
    }

　　根據評分矩陣來按最低順序獲得最佳對齊：

    /**
     * 一種計算兩個的最佳全局最低對齊序列的實用方法
     *
     * @param similarityMatrix 2維數組表示以前計算的全局最大相似度矩陣。
     * @param sequence1 代表第一個氨基酸序列的字符串。
     * @param sequence1 代表第二個氨基酸序列的字符串。
     * @param matchScore 代表分配給匹配的分數的浮點數。
     * @param mismatchScore 一個浮點數，代表分配給不匹配項的分數。
     * @param indelScore 浮點數表示分配給插入和刪除的分數。
     */
    public static OptimalAlignment obtainOptimalAlignmentByDownmostOrder(float similarityMatrix[][],
                                                                         String sequence1, String sequence2, double matchScore, double mismatchScore, double indelScore) {

        int i = similarityMatrix.length - 1;
        int j = similarityMatrix[0].length - 1;
        StringBuilder alignment1Builder = new StringBuilder();
        StringBuilder alignment2Builder = new StringBuilder();

        sequence1 = "-" + sequence1;
        sequence2 = "-" + sequence2;

        while (i != 0 && j != 0) {
            if (similarityMatrix[i][j] == (similarityMatrix[i][j - 1] + indelScore)) {
                alignment1Builder.insert(0, "-");
                alignment2Builder.insert(0, sequence2.charAt(j));
                j--;
            } else if (similarityMatrix[i][j] == (similarityMatrix[i - 1][j - 1] + (sequence1.charAt(i) == sequence2.charAt(j) ? matchScore : mismatchScore))) {
                alignment1Builder.insert(0, sequence1.charAt(i));
                alignment2Builder.insert(0, sequence2.charAt(j));
                i--;
                j--;
            } else if (similarityMatrix[i][j] == (similarityMatrix[i - 1][j] + indelScore)){
                alignment1Builder.insert(0, sequence1.charAt(i));
                alignment2Builder.insert(0, "-");
                i--;
            }
        }

        return new OptimalAlignment(alignment1Builder.toString(), alignment2Builder.toString());
    }

　　測試主函數：

　　public static void main(String[] args) {
        String sequence1 = "GCATGCU";
        String sequence2 = "GATTACA";

        float matchScore = 1;
        float mismatchScore = -1;
        float indelScore = -1;

        float[][] computedMatrix = NeedlemanWunsch.computeMatrix(sequence1, sequence2, matchScore, mismatchScore, indelScore);

        System.out.println("Best global downmost alignment: ");
        System.out.println(NeedlemanWunsch.obtainOptimalAlignmentByDownmostOrder(computedMatrix, sequence1, sequence2, matchScore, mismatchScore, indelScore));
    }

　　全局相似度矩陣：

 　　結果：