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