算法笔记学习000——Smith-Waterman算法寻找两个字符串中匹配度最高的子串

<span style="font-family: Arial, Helvetica, sans-serif; background-color: rgb(255, 255, 255);"></span>

<span style="font-family: Arial, Helvetica, sans-serif; background-color: rgb(255, 255, 255);"> </span>

<span style="font-family: Arial, Helvetica, sans-serif; background-color: rgb(255, 255, 255);">//第一次在CSDN上写笔记，是想着身为一个程序员，有些时候有些知识是转瞬即逝的，现在理解了说不定会很快忘记，所以记录一下加深印象，也有助于以后复习~</span>

<span style="font-family: Arial, Helvetica, sans-serif; background-color: rgb(255, 255, 255);">//本文参考了著名经典论文<em>Identification of Common Molecular Subsequences,</em>【其实根本就是照着人家论文学的嘛喂！</span>

<span style="font-family: Arial, Helvetica, sans-serif; background-color: rgb(255, 255, 255);">Smith-Waterman算法与1981年提出，出发点是信息生物学，意欲寻找两个碱基序列中匹配度最高的两个片段，以提取出有价值的生物学信息。</span>

Smith-Waterman算法的基本思路是，对于给定的长度为m的串s1(a1,a2,a3,...,am)和长度为n的串s2(b1,b2,...,bn)，构造一个(m+1)*(n+1)的得分矩阵H。H中的元组H(i,j)的含义是，匹配子串（a1,a2,..,ai）和（b1,b2,...bj）的最高匹配度。

下面我们来讨论匹配度的定义。

我们知道，当匹配进行到(ai,bj)元组时，我们有三种选择:

1. 在bj元组前插入一个空格，用这个空格去匹配ai，匹配后的形式如下，此种匹配方式得分为-4/3

...,ai,...

   ...,--,bj...

2.在ai元组前插入一个空格，用这个空格去匹配bj，匹配后的形式如下，此种匹配方式得分为-4/3

       ...,--,ai,...

   ...,bj,...

3.直接匹配(ai,bj)，如果ai与bj相同，则得分为1，否则得分为-1/3。

那么，匹配度H(i,j)，就应该是在保证匹配度大于0的前提下这三种匹配方式中得分最高者。于是我们有了一个递归的定义：

H(i,j) = max{H(i-1,j-1) + MATCH(或者DISMATCH), H(i-1,j) - 4, H(i, j-1) - 4,0}

其中0<=i<=m, 0<=j<=n.

矩阵的第一行和第一列都为0，即H(0,j)和H(i,0)都等于零，相当于尚未开始匹配的匹配度。

有了匹配度矩阵H，我们的目标就是找到H中得分最大的元组H(i,j)，并用回溯的方法找到一个轨迹，按照上述的评分规则沿轨迹从头到尾走一遍，就可以得到那个最高匹配度H(i,j)，也就是找到s1,s2的匹配过程，自然也就找到了两个匹配度最高的子序列。

下面举例子来说明，这个例子也来源于Smith-Waterman的论文原文。我们假设需要匹配的两个序列分别为s1=AAUGCCAUUGACGG,S2=ACAGCCUCGCUUAG。

首先，计算匹配度矩阵H：

找到矩阵中得分最大(3.3)的元组H(10,8)，开始回溯的过程。

其实回溯的思路很简单，就是检查位于该元组上方，左方，和左上方的元组，看它的得分是等于上-4/3，还是左-4/3，还是左上+1，还是左上-1/3。简而言之，就是看看这个元组是“从谁那儿走过来的”。

回溯终止的临界条件是，某个元组的得分为0，这意味着我们尚未找到匹配这两个串的子串头。

整个回溯过程结束后，找到的子串如下：

下面是用java语言写的源代码，大家可以到自己机器上试着跑一下~

import java.io.BufferedReader;import java.io.IOException;import java.io.InputStreamReader;import java.util.ArrayList;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 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<Character>();        ArrayList<Character> arr2 = new ArrayList<Character>();        while(!stk1.empty())            arr1.add(stk1.pop());        subSq1 = arr1.toString();        while(!stk2.empty())            arr2.add(stk2.pop());        subSq2 = arr2.toString();    }    public static void main(String[] args) throws IOException {        SWSq x = new SWSq();        String s1,s2;        BufferedReader br = new BufferedReader(new InputStreamReader(System.in));        System.out.println("Please enter one String:");        s1 = br.readLine();        System.out.println("Enter another String:");        s2 = br.readLine();        x.find(s1, s2);        System.out.println("The subsequences with greatest similarity are " + x.subSq1 + " and " + x.subSq2);    }}