Needleman-Wunsch 算法和Smith-Waterman算法

全局和局部序列比对:

      全局序列比对  尝试找到两个完整的序列 S1 和 S2 之间的最佳比对。如S1=GCCCTAGCG S2=GCGCAATG 如果设定每个匹配字符为1分，每个空格为-2分，每个不匹配为-1分，则下面的比对就是全局最优比对：S1'=GCCCTAGCG S2'=GCGC_AATG,连字符“_”代表空格。在 S2' 中有五个匹配字符，一个空格（或者反过来说，在 S1' 中有一个插入项），有三个不匹配字符。这样得到的分数是 (5×1) + (1×-2) + (3×-1) = 0，这是能够实现的最佳结果。

 

      局部序列比对 不必对两个完整的序列进行比对；可以在每个序列中使用某些部分来获得最大得分。使用同样的序列 S1 和 S2，以及同样的得分方案，可以得到以下局部最优比对 S1'' 和 S2''：S1=GCCCTAGCG S2=GCGCAATG   S1''=GCG S2''=GCG，虽然这个局部比对恰好没有不匹配字符或空格，但是一般情况下，局部比对可能存在不匹配字符或空格。这个局部比对的得分是 (3×1) + (0×-2) + (0×-1) = 3。（最佳局部比对的得分要大于或等于最佳全局比对的得分，这是因为全局比对也属于局部比对)

 

 Needleman-Wunsch 算法：

            该算法用来计算全局比对。它的思路与 LCS 算法相似。这个算法也使用二维表格，一个序列沿顶部展开，一个序列沿左侧展开。而且也能通过以下三个途径到达每个单元格：1.来自上面的单元格，代表将左侧的字符与空格比对。2.来自左侧的单元格，代表将上面的字符与空格比对。3.来自左上侧的单元格，代表与左侧和上面的字符比对（可能匹配也可能不匹配）。该单元格的值来自于一下3个中的最大值：（1）上方的值-2 （2）左边的值-2 （3）如果该单元格所在的行于所在的列对应的字符相等，则为左上值加1，否则为左上值-1。

        首先，必须初始化表格。这意味着填充第二行和第二列的分数和指针。填充第二行的操作意味着使用位于顶部的第一个序列中的字符，并使用空格，而不是使用左侧从上到下的序列中的第一个字符。空格的扣分是 -2，所以每次使用空格的时候，就给以前的单元格加了 -2 分。以前的单元格是左侧的单元格。这就说明了在第二行中为什么得到了 0, -2, -4, -6, ... 这样的序列。用相似的方式得到第二列的得分和指针.

       接下来，需要填充余下的单元格。同 LCS 算法一样，对于每个单元格，都有三个选择，要从中选择最大的。可以从上面、左侧、左上侧到达每个单元格。假设 S1 和 S2 是要比对的字符串，S1' 和 S2' 是生成的比对中的字符串。从上面到达单元格相当于将左面的字符从 S2 加入 S2'，跳过上面的 S1 中的当前字符，并在 S1'中加入一个空格。因为一个空格的分数是 -2，所以当前单元格的得分要从上面的单元格得分减 2 得到。类似的，将左边的单元格得分减 2，可以从左侧到达空单元格。

          最后，可以将上面的字符加入到 S1' 中，将左边的字符加入到 S2' 中。这就相当于从左上侧进入空白单元格。这两个字符将会匹配，在这种情况下，新的得分就是左上侧单元格的得分减 1。在这三种可能性当中，选择得分最大的一个（如果得分相等，可以从得分高的单元格中从任选一个）。接下来，需要得到实际的比对字符串S1' 和S2'以及比对的得分。右下角单元格中的得分包含 S1 和 S2 的最大比对得分，就像在 LCS 算法中包含 LCS 的长度一样。而且，与 LCS 算法类似，要获得 S1' 和 S2'，要从右下角单元格开始沿着指针回溯，反向构建 S1' 和 S2'。从表格的构建过程可知，从上向下对应着将左侧字符从 S2 加入到 S2' 中，将空格加入 S1' 中；从左向右对应着将上面的字符从 S1 加入到 S1' 中，将空格加入 S2' 中；而向下和向右移动意味着分别将来自S1 和 S2 的字符加入 S1' 和 S2'。

 

 

Smith-Waterman 算法：在 Smith-Waterman 算法中，不必比对整个序列。两个零长字符串即为得分为 0 的局部比对，这一事实表明在构建局部比对时，不需要使用负分。这样会造成进一步比对所得到的分数低于通过 “重设” 两个零长字符串所能得到的分数。而且，局部比对不需要到达任何一个序列的末端，所以也不需要从右下角开始回溯：可以从得分最高的单元格开始回溯。

          这导致 Smith-Waterman 算法与 Needleman-Wunsch 算法存在着三个区别。

           首先，在初始化阶段，第一行和第一列全填充为 0（而且第一行和第一列的指针均为空）。

           第二，在填充表格时，如果某个得分为负，那么就用 0 代替，只对得分为正的单元格添加返回指针。

           最后，在回溯的时候，从得分最高的单元格开始，回溯到得分为 0 的单元格为止。除此之外，回溯的方式与 Needleman-Wunsch 算法完全相同。

每个单元格的值来自于一下的最大值：：（1）上方的值-2 （2）左边的值-2 （3）如果该单元格所在的行于所在的列对应的字符相等，则为左上值加1，否则为左上值-1。但是切记，如果最大值为负数，则该单元格的值用0代替，不指向任何单元格。

 

两个算法的Java代码实现：

[java]  view plain copy

1. public abstract class SequenceAlignment extends DynamicProgramming {  
2.     protected int match; //当匹配时的分值  
3.     protected int mismatch; // 当不匹配时的分值  
4.     protected int space; // 当为空格时的分值  
5.     protected String[] alignments;  
6.     public SequenceAlignment(String sequence1, String sequence2) {  
7.         this(sequence1, sequence2, 1, -1, -1);  
8.     }  
9.     public SequenceAlignment(String sequence1, String sequence2, int match,  
10.             int mismatch, int gap) {  
11.         super(sequence1, sequence2);  
12.         this.match = match;  
13.         this.mismatch = mismatch;  
14.         this.space = gap;  
15.     }  
16.     /* 
17.      * (non-Javadoc) 两个算法通用的回溯 
18.      * 
19.      * @see com.ibm.combio.DynamicProgramming#getTracebace() 
20.      */  
21.     protected Object getTracebace() {  
22.         StringBuffer align1Buf = new StringBuffer();  
23.         StringBuffer align2Buf = new StringBuffer();  
24.         Cell currentCell = getTracebackStartingCell();  
25.         while (traceBackIsNotDone(currentCell)) {  
26.             if (currentCell.getRow() - currentCell.getPrevCell().getRow() == 1) {  
27.                 align2Buf.insert(0, sequence2.charAt(currentCell.getRow() - 1));  
28.             } else {  
29.                 align2Buf.insert(0, '-');  
30.             }  
31.             if (currentCell.getCol() - currentCell.getPrevCell().getCol() == 1) {  
32.                 align1Buf.insert(0, sequence1.charAt(currentCell.getCol() - 1));  
33.             } else {  
34.                 align1Buf.insert(0, '-');  
35.             }  
36.             currentCell = currentCell.getPrevCell();  
37.         }  
38.         String[] alignments = new String[] { align1Buf.toString(),  
39.                 align2Buf.toString() };  
40.         return alignments;  
41.     }  
42.     public int getAlignmentScore() {  
43.         if (alignments == null) {  
44.             getAlignment();  
45.         }  
46.         int score = 0;  
47.         for (int i = 0; i < alignments[0].length(); i++) {  
48.             char c1 = alignments[0].charAt(i);  
49.             char c2 = alignments[1].charAt(i);  
50.             if (c1 == '-' || c2 == '-') {  
51.                 score += space;  
52.             } else if (c1 == c2) {  
53.                 score += match;  
54.             } else {  
55.                 score += mismatch;  
56.             }  
57.         }  
58.         return score;  
59.     }  
60.     public String[] getAlignment() {  
61.         ensureTableIsFillIn();  
62.         alignments = (String[]) getTracebace();  
63.         return alignments;  
64.     }  
65.     protected abstract boolean traceBackIsNotDone(Cell currentCell); // 回溯停止判断  
66.     protected abstract Cell getTracebackStartingCell(); // 取得回溯起始单元  
67. }  
68. 2.NeedlemanWunsch 算法类  
69. public class NeedlemanWunsch extends SequenceAlignment {  
70.     public NeedlemanWunsch(String sequence1, String sequence2) {  
71.         super(sequence1, sequence2);  
72.     }  
73.     public NeedlemanWunsch(String sequence1, String sequence2, int match,  
74.             int mismatch, int gap) {  
75.         super(sequence1, sequence2, match, mismatch, gap);  
76.     }  
77.     protected void fillInCell(Cell currentCell, Cell cellAbove,  
78.             Cell cellToLeft, Cell cellAboveLeft) {  
79.         int rowSpaceScore = cellAbove.getScore() + space;  
80.         int colSpaceScore = cellToLeft.getScore() + space;  
81.         int matchOrMismatchScore = cellAboveLeft.getScore();  
82.         if (sequence2.charAt(currentCell.getRow() - 1) == sequence1  
83.                 .charAt(currentCell.getCol() - 1)) {  
84.             matchOrMismatchScore += match;  
85.         } else {  
86.             matchOrMismatchScore += mismatch;  
87.         }  
88.         if (rowSpaceScore >= colSpaceScore) {  
89.             if (matchOrMismatchScore >= rowSpaceScore) {  
90.                 currentCell.setScore(matchOrMismatchScore);  
91.                 currentCell.setPrevCell(cellAboveLeft);  
92.             } else {  
93.                 currentCell.setScore(rowSpaceScore);  
94.                 currentCell.setPrevCell(cellAbove);  
95.             }  
96.         } else {  
97.             if (matchOrMismatchScore >= colSpaceScore) {  
98.                 currentCell.setScore(matchOrMismatchScore);  
99.                 currentCell.setPrevCell(cellAboveLeft);  
100.             } else {  
101.                 currentCell.setScore(colSpaceScore);  
102.                 currentCell.setPrevCell(cellToLeft);  
103.             }  
104.         }  
105.     }  
106.     @Override  
107.     protected boolean traceBackIsNotDone(Cell currentCell) {  
108.         return currentCell.getPrevCell() != null;  
109.     }  
110.     @Override  
111.     protected Cell getTracebackStartingCell() {  
112.         return scoreTable[scoreTable.length - 1][scoreTable[0].length - 1];  
113.     }  
114.     /* 
115.      * (non-Javadoc) 
116.      * 
117.      * @see java.lang.Object#toString() 
118.      */  
119.     @Override  
120.     public String toString() {  
121.         return "[NeedlemanWunsch: sequence1=" + sequence1 + ", sequence2="  
122.                 + sequence2 + "]";  
123.     }  
124.     /* (non-Javadoc) 
125.      * 初始化单元格的前驱指针 
126.      * @see com.ibm.combio.DynamicProgramming#getInitialPointer(int, int) 
127.      */  
128.     protected Cell getInitialPointer(int row, int col) {  
129.         if (row == 0 && col != 0) {  
130.             return scoreTable[row][col - 1];  
131.         } else if (col == 0 && row != 0) {  
132.             return scoreTable[row - 1][col];  
133.         } else {  
134.             return null;  
135.         }  
136.     }  
137.     /* (non-Javadoc) 
138.      * 初始化单元格的值 
139.      * @see com.ibm.combio.DynamicProgramming#getInitialScore(int, int) 
140.      */  
141.     protected int getInitialScore(int row, int col) {  
142.         if (row == 0 && col != 0) {  
143.             return col * space;  
144.         } else if (col == 0 && row != 0) {  
145.             return row * space;  
146.         } else {  
147.             return 0;  
148.         }  
149.     }  
150.     @Override  
151.     public Object getTracebace() {  
152.         return super.getTracebace();  
153.     }  
154. }  
155. 3.SmithWaterman 类  
156. public class SmithWaterman extends SequenceAlignment {  
157.     private Cell highScoreCell; // 最高分值得单元  
158.     public SmithWaterman(String sequence1, String sequence2) {  
159.         super(sequence1, sequence2);  
160.     }  
161.     public SmithWaterman(String sequence1, String sequence2, int match,  
162.             int mismatch, int gap) {  
163.         super(sequence1, sequence2, match, mismatch, gap);  
164.     }  
165.     @Override  
166.     protected void initialize() {  
167.         // TODO Auto-generated method stub  
168.         super.initialize();  
169.         highScoreCell = scoreTable[0][0];  
170.     }  
171.     @Override  
172.     protected Cell getTracebackStartingCell() {  
173.         return highScoreCell;  
174.     }  
175.     @Override  
176.     protected boolean traceBackIsNotDone(Cell currentCell) {  
177.         return currentCell.getScore() != 0;  
178.     }  
179.     @Override  
180.     public Object getTracebace() {  
181.         return super.getTracebace();  
182.     }  
183.     @Override  
184.     protected void fillInCell(Cell currentCell, Cell cellAbove,  
185.             Cell cellToLeft, Cell cellLeftAbove) {  
186.         int rowSpaceScore = cellAbove.getScore() + space;  
187.         int colSpaceScore = cellToLeft.getScore() + space;  
188.         int matchOrMismatchScore = cellLeftAbove.getScore();  
189.         if (sequence1.charAt(currentCell.getCol() - 1) == sequence2  
190.                 .charAt(currentCell.getRow() - 1)) {  
191.             matchOrMismatchScore += match;  
192.         } else {  
193.             matchOrMismatchScore += mismatch;  
194.         }  
195.         if (rowSpaceScore >= colSpaceScore) {  
196.             if (matchOrMismatchScore >= rowSpaceScore) {  
197.                 if (matchOrMismatchScore > 0) {  
198.                     currentCell.setScore(matchOrMismatchScore);  
199.                     currentCell.setPrevCell(cellLeftAbove);  
200.                 }  
201.             } else {  
202.                 if (rowSpaceScore > 0) {  
203.                     currentCell.setScore(rowSpaceScore);  
204.                     currentCell.setPrevCell(cellAbove);  
205.                 }  
206.             }  
207.         } else {  
208.             if (matchOrMismatchScore >= colSpaceScore) {  
209.                 if (matchOrMismatchScore > 0) {  
210.                     currentCell.setScore(matchOrMismatchScore);  
211.                     currentCell.setPrevCell(cellLeftAbove);  
212.                 }  
213.             } else {  
214.                 if (colSpaceScore > 0) {  
215.                     currentCell.setScore(colSpaceScore);  
216.                     currentCell.setPrevCell(cellToLeft);  
217.                 }  
218.             }  
219.         }  
220.         if (currentCell.getScore() > highScoreCell.getScore()) {  
221.             highScoreCell = currentCell;  
222.         }  
223.     }  
224.     @Override  
225.     protected Cell getInitialPointer(int row, int col) {  
226.         return null;  
227.     }  
228.     @Override  
229.     protected int getInitialScore(int row, int col) {  
230.         return 0;  
231.     }  
232. }  
233. 4.编写测试类  
234. public class test {  
235.     public static void main(String[] args) {  
236.         String s1 = "GCCCTAGCG";  
237.         String s2 = "GCGCAATG";  
238.         SmithWaterman sw = new SmithWaterman(s1, s2);  
239.         sw.getScoreTable();  
240.         String[] s = (String[]) sw.getTracebace();  
241.         for (String str : s)  
242.             System.out.println(str);  
243.         NeedlemanWunsch nw = new NeedlemanWunsch(s1, s2);  
244.         nw.getScoreTable();  
245.         s = (String[]) sw.getTracebace();  
246.         for (String str : s)  
247.             System.out.println(str);  
248.     }  
249. }  

1.两个算法通用的部分定义为一个抽象类，该类继承在LCS算法中定义的动态编程框架抽象类未进行优化，查找全局和局部比对的时间为 O(mn)，LCS及这两个算法都是动态编程的应用。动态编程还可用于矩阵链相乘、装配线规划和计算机象棋程序。使用动态编程能够解决的LCS及以上两个问题，这些问题都有共同的特征：每个问题的解都能用递归关系表示;用递归方法对这些递归关系的直接实现会造成解决方案效率低下，因为其中包含了对子问题的多次计算;一个问题的最优解能够用原始问题的子问题的最优解构造得到.

 

 

 

http://en.wikipedia.org/wiki/Needleman%E2%80%93Wunsch_algorithm

 

http://en.wikipedia.org/wiki/Smith%E2%80%93Waterman_algorithm

 

http://www.w3school.com.cn/php/php_ref_mysql.asp

 

http://hi.baidu.com/%CF%FE%D4%C2%B7%C9%B7%C9/blog/item/693e1481c53417dd9123d93d.html