Smith-Waterman算法及其Java实现

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);
    }
}