一、问题描述
尝试找到两个完整的序列 S1 和 S2 之间的最佳比对。如S1=GCCCTAGCG S2=GCGCAATG 如果设定每个匹配字符为1分,每个空格为-2分,每个不匹配为-1分,则下面的比对就是全局最优比对:S1’=GCCCTAGCG S2’=GCGC_AATG,连字符“_”代表空格。在 S2’ 中有五个匹配字符,一个空格(或者反过来说,在 S1’ 中有一个插入项),有三个不匹配字符。这样得到的分数是 (5×1) + (1×-2) + (3×-1) = 0,这是能够实现的最佳结果。
二、设计算法
1. 算法策略
Needleman-Wunsch全局序列比对算法:
它的思路与 LCS 算法相似。这个算法也使用二维表格,一个序列沿顶部展开,一个序列沿左侧展开。而且也能通过以下三个途径到达每个单元格:1.来自上面的单元格,代表将左侧的字符与空格比对。2.来自左侧的单元格,代表将上面的字符与空格比对。3.来自左上侧的单元格,代表与左侧和上面的字符比对(可能匹配也可能不匹配)。该单元格的值来自于以下3个中的最大值:(1)上方的值-2 (2)左边的值-2 (3)如果该单元格所在的行于所在的列对应的字符相等,则为左上值加1,否则为左上值-1。
2. 数据结构
① 得分矩阵
二维数组dp表示得分矩阵,s1位于dp数组的上方,s2位于dp数组的左方其中dp[i][j] 表示s1[0:j]与s2[0:i]的最高得分。
② 父节点矩阵
二维数组status记录当前节点的父节点。status[i][j] 记录dp[i][j] 的父节点。1代表左上角,2代表左边,4代表上边,若有多解也可表示,3=1+2,5=1+4,6=2+4,7=1+2+4。
3. 求解步骤
① 初始化得分矩阵
初始化第0行:所有节点的得分依次-2,且父节点均来自左边,status[0][j] = 2;
初始化第0列:所有节点的得分依次-2,且父节点均来自上边,status[i][0] = 4。
② 计算得分矩阵的每一项
状态转移方程(得分体系):dp[i][j] = max{ 1. dp[i-1][j-1] +/- 1 左上角
2. dp[i][j-1] - 2 左边s2为空
4. dp[i-1][j] - 2 上边s1为空
}
并填充status[i][j]。
③ 回溯,得到最优解
通过status数组进行回溯构造所有解,进行字符串拼接,最后反转得出所有正确解。
④ 输出打印
最优值(最高得分) = dp[n][m],打印由③所得的最优解s1和s2。
三、实现算法
1. Main.java
package com.全局序列比对;
import java.util.Scanner;
/**
* @author fanb
* @email fanb@nwafu.edu.cn
* 问题描述
* 尝试找到两个完整的序列 S1 和 S2 之间的最佳比对。
* 如果设定每个匹配字符为1分,每个空格为-2分,每个不匹配为-1分,
*
* 求解策略
* Needleman-Wunsch(尼德曼-翁施)算法:
* 它的思路与 LCS 算法相似。这个算法也使用二维表格,一个序列沿顶部展开,一个序列沿左侧展开
* s1位于dp数组的上方,s2位于dp数组的左方
* dp[i][j]得分矩阵 表示s1[0:j] 与 s2[0:i]的最高得分
* 状态转移方程(得分体系) dp[i][j] = max{
* 1. dp[i-1][j-1] +/- 1 左上角
* 2. dp[i][j-1] - 2 左边s2为空
* 4. dp[i-1][j] - 2 上边s1为空
* }
* 回溯找出最优解 status[i][j] 表示dp[i][j] 的父节点 1.左上角 2.左边 4.上边
* 若有多解也可表示:3=1+2,5=1+4,6=2+4,7=1+2+4
*
* 时间复杂度 O(m*n)
* 空间复杂度 O(m*n)
*/
public class Main {
public static void main(String[] args) {
Scanner scanner = new Scanner(System.in);
System.out.print("输入s1:");
String s1 = scanner.next();
System.out.print("输入s2:");
String s2 = scanner.next();
scanner.close();
Solution solution = new Solution(s1, s2);
int ans = solution.needlemanWunsch(); // 求解最优值
System.out.println("最优值:"+ans);
solution.print(); // 打印最优解
}
}
class Solution{
String s1,s2; // s1上侧 s2左侧
int m, n; // m列 n行
int[][] dp; // 得分矩阵
int[][] status; // 记录父节点,1(2^0)代表左上,2(2^1)代表来自左边(左空格),4(2^2)代表来自上边(上空格)
// 构造器
public Solution(String s1, String s2) {
this.s1 = s1;
this.s2 = s2;
this.m = s1.length();
this.n = s2.length();
this.dp = new int[n+1][m+1];
this.status = new int[n+1][m+1];
status[0][0] = 0;
dp[0][0] = 0;
}
// Needleman-Wunsch 全局序列比对算法
public int needlemanWunsch(){
// 初始化第0列
for (int i = 1; i <= n; i++) {
dp[i][0] = dp[i-1][0] - 2;
status[i][0] = 4;
}
// 初始化第0行
for (int i = 1; i <= m; i++) {
dp[0][i] = dp[0][i-1] - 2;
status[0][i] = 2;
}
// 动态规划,填满dp数组
for (int i = 1; i <= n; i++) { // 行
for (int j = 1; j <= m; j++) { // 列
int leftTop = dp[i-1][j-1] + (s1.charAt(j-1) == s2.charAt(i-1) ? 1 : -1); // 左上角
int left = dp[i][j-1] - 2; // 来自左边 左空格
int top = dp[i-1][j] - 2; // 来自上边 上空格
int state = 0;
int max = Math.max(leftTop, left); // 当前位置最高得分
max = Math.max(top,max);
// 左上角
if(leftTop == max) state += 1;
// 来自左边 左空格
if(left == max) state += 2;
// 来自上边 上空格
if(top == max) state += 4;
status[i][j] = state;
dp[i][j] = max;
// System.out.print(max + " ");
}
// System.out.println();
}
return dp[n][m];
}
// print 回溯dp二维数组输出最优解
public void print(){
StringBuilder s11 = new StringBuilder();
StringBuilder s22 = new StringBuilder();
boolean kind = false; // 标志位,判断是否有多种情况
int i, j;
while(true){
s11.append(' ');
s22.append(' ');
kind = false;
i = n; // 行
j = m; //列
while(status[i][j] != 0){
switch (status[i][j]){
case 1: // 只来自左上角
s11.append(s1.charAt(j-1));
s22.append(s2.charAt(i-1));
i--;
j--;
break;
case 2: // 只来自左空格
s11.append(s1.charAt(j-1));
s22.append('_');
j--;
break;
case 3: // 来自 1 + 2
// 若kind已经为true,说明已经走了一个分支,此时不需要-2
if(!kind){
status[i][j] -= 2;
kind = true;
}
s11.append(s1.charAt(j-1));
s22.append('_');
j--;
break;
case 4: // 只来自上空格
s11.append('_');
s22.append(s2.charAt(i-1));
i--;
break;
case 5: // 来自 1 + 4
case 6: // 来自 2 + 4
case 7: // 来自 1 + 2 + 4
// 若kind已经为true,说明已经走了一个分支,此时不需要-4
if(!kind){
status[i][j] -= 4;
kind = true;
}
s11.append('_');
s22.append(s2.charAt(i-1));
i--;
break;
}
}
if(!kind) break;
}
System.out.println("s1:" + s11.reverse());
System.out.println("s2:" + s22.reverse());
}
}
2. 测试
测试用例
S1=GCCCTAGCG
S2=GCGCAATG
测试结果(存在多解)
最优值:0
s1: GCCCTAGCG GCCCTAGCG GCCCTAGCG
s2: GCGC_AATG GCGCAA_TG GCGCAAT_G
四、复杂度分析
(1)时间复杂度
初始化第0行T(m),初始化第0列T(n),填充得分矩阵T(mn)
T = T(mn) + T(m) + T(n)
所以时间复杂度:O(m*n)
(2)空间复杂度
得分矩阵dp[n+1][m+1]占(n+1)(m+1),父节点矩阵status[n+1][m+1]占(n+1)(m+1)
S = 2mn + 2m + 2n + 2
所以空间复杂度:O(m*n)