Smith-Waterman演算法及其Java實現

標籤：stream   賦值   sem   最大   print   new   blog   for   put   

Smith-Waterman演算法是1981年Smith和Waterman提出的一種用來尋找並比較具有局部相似性地區的動態規划算法，很多後來的演算法都是在該演算法的基礎上發展的。這是一種兩序列局部比對演算法，把兩條未知的序列進行排列，通過字母的匹配，刪除和插入操作，使得兩條序列達到同樣長度，在操作的過程中，儘可能保持相同的字母對應在同一個位置。當兩條序列進行比對時，找出待比對序列中的某一子片段的最優比對。這種比對方法可能會揭示一些匹配的序列段，而本來這些序列段是被一些完全不相關的殘基所淹沒的。

其演算法過程簡單描述為：
1) 為每一堿基對或殘基對賦值。相同或類似的賦予正值，對於不同的或有空位的賦予負值；
2) 用0對矩陣邊緣單元初始化；
3) 矩陣中得分值相加，任何小於0的得分值均用0代替；
4) 通過動態規劃的方法，從矩陣中的最大分值單元開始回溯尋找；
5) 繼續，一直到分值為0的單元停止，此回溯路徑的單元即為最優比對序列。
由以上可知，Smith-Waterman演算法主要分兩步．計算得分矩陣和尋找最佳相似片段對。得到得分矩陣以後，用動態規劃回溯的方法找到局部最大相似片段對：先找到得分矩陣中最大的元素．然後按照元素原路徑一步一步往前回溯，直到回溯到0時停止。

 

下面舉例子來說明，這個例子也來源於Smith-Waterman的論文原文。
1) 我們假設需要匹配的兩個序列分別為s1=AAUGCCAUUGACGG,S2=ACAGCCUCGCUUAG。
2) 首先，計算匹配度矩陣H。找到矩陣中得分最大(3.3)的元組H(10,8)，開始回溯的過程。
3) 回溯的思路很簡單，就是檢查位於該元組上方，左方，和左上方的元組，看它的得分是等於上-4/3，還是左-4/3，還是左上+1，還是左上-1/3。簡而言之，就是看看這個元組是“從誰那兒走過來的”。
4) 回溯終止的臨界條件是，某個元組的得分為0，這意味著我們尚未找到匹配這兩個串的子串頭。
5) 整個回溯過程結束後，找到的子串如下：

AAUGCCAUUG
ACAGCC-UCG

下面是用Java語言寫的原始碼：

import java.io.BufferedReader;import java.io.IOException;import java.io.InputStreamReader;import java.util.ArrayList;import java.util.Iterator;import java.util.Stack;public class SWSq {    private int[][] H;    private int[][] isEmpty;    private static int SPACE ;                      //空格匹配的得分      private static int MATCH ;                       //兩個字母相同的得分      private static int DISMACH;                    //兩個字母不同的得分      private int maxIndexM, maxIndexN;    private Stack<Character> stk1, stk2;    public String subSq1, subSq2;                        //相似性最高的兩個子串      public SWSq(){        stk1 = new Stack<Character>();        stk2 = new Stack<Character>();        SPACE = -4;        MATCH = 3;        DISMACH = -1;    }    private int max(int a, int b, int c){        int maxN;        if(a >= b)            maxN = a;        else            maxN = b;        if(maxN < c)            maxN = c;        if(maxN < 0)            maxN = 0;        return maxN;    }    private void calculateMatrix(String s1, String s2, int m, int n){//計算得分矩陣          if(m == 0)            H[m][n] = 0;        else if(n == 0)            H[m][n] = 0;        else{            if(isEmpty[m - 1][n - 1] == 1)                calculateMatrix(s1, s2, m-1, n-1);            if(isEmpty[m][n - 1] == 1)                calculateMatrix(s1, s2, m, n-1);            if(isEmpty[m - 1][n] == 1)                calculateMatrix(s1, s2, m-1, n);            if(s1.charAt(m-1) == s2.charAt(n-1))                H[m][n] = max(H[m - 1][n - 1] + MATCH, H[m][n - 1] + SPACE, H[m - 1][n] + SPACE);            else                H[m][n] = max(H[m - 1][n - 1] + DISMACH, H[m][n - 1] + SPACE, H[m - 1][n] + SPACE);        }        isEmpty[m][n] = 0;    }    private void findMaxIndex(int[][] H, int m, int n){//找到得分矩陣H中得分最高的元組的下標          int curM, curN, i, j, max;        curM = 0;        curN = 0;        max = H[0][0];        for(i = 0; i < m; i++)            for(j = 0; j < n; j++)                if(H[i][j] > max){                    max = H[i][j];                    curM = i;                    curN = j;                }        maxIndexM = curM;        maxIndexN = curN;    }    private void traceBack(String s1, String s2, int m, int n){//回溯 尋找最相似子序列          if(H[m][n] == 0)            return;        if(H[m][n] == H[m-1][n] + SPACE) {            stk1.add(s1.charAt(m-1));            stk2.add(‘-‘);            traceBack(s1, s2, m - 1, n);        }        else if(H[m][n] == H[m][n-1] + SPACE) {            stk1.add(‘-‘);            stk2.add(s2.charAt(n-1));            traceBack(s1, s2, m, n - 1);        }        else {            stk1.push(s1.charAt(m - 1));            stk2.push(s2.charAt(n-1));            traceBack(s1, s2, m - 1, n - 1);        }    }    public String ALtoString(ArrayList<Character> A) {        StringBuilder sb = new StringBuilder();        for (Character a : A) {            sb.append(a.toString());        }        return sb.toString();    }    public void find(String s1, String s2){        //initMatrix(s1.length(), s2.length());          int i, j;        H = new int[s1.length() + 1][s2.length() + 1];        isEmpty = new int[s1.length() + 1][s2.length() + 1];        for(i = 0; i<=s1.length(); i++)            for(j = 0; j<=s2.length(); j++)                isEmpty[i][j] = 1;        calculateMatrix(s1, s2, s1.length(), s2.length());        findMaxIndex(H, H.length, H[0].length);        traceBack(s1, s2, maxIndexM, maxIndexN);        ArrayList<Character> arr1 = new ArrayList<>();        ArrayList<Character> arr2 = new ArrayList<>();        while(!stk1.empty())            arr1.add(stk1.pop());        subSq1 = ALtoString(arr1);        while(!stk2.empty())            arr2.add(stk2.pop());        subSq2 = ALtoString(arr2);    }    public static void main(String[] args) throws IOException {        SWSq x = new SWSq();        String s1 = "AAUGCCAUUGACGG";        String s2 = "ACAGCCUCGCUUAG";        x.find(s1, s2);        System.out.println("----------------------------");        System.out.println(s1);        System.out.println(s2);        System.out.println("----------------------------");        System.out.println(x.subSq1);        System.out.println(x.subSq2);    }}  

　　

 

Smith-Waterman演算法及其Java實現