全局序列比对动态规划找到矩阵中的最小值

Question

我有 2 个序列，AACAGTTACC 和 TAAGGTCA，我正在尝试查找全局序列比对。我设法创建了一个二维数组并创建了矩阵，我什至用半动态方法填充了它。

这是我填充矩阵的代码：

void process() {
    for (int i = 1; i <= sequenceA.length; i++) {
        for (int j = 1; j <= sequenceB.length; j++) {
            int scoreDiag = opt[i-1][j-1] + equal(i, j);
            int scoreLeft = opt[i][j-1] - 1;
            int scoreUp = opt[i-1][j] - 1;
            opt[i][j] = Math.max(Math.max(scoreDiag, scoreLeft), scoreUp);
        }
    }
}

private int equal(int i, int j) {
    if (sequenceA[i - 1] == sequenceB[j - 1]) {
        return 1;
    } else {
        return -1;
    }
}

我的主要问题是这段代码生成了这个输出：

 0   -1   -2   -3   -4   -5   -6   -7   -8     
-1   -1    0   -1   -2   -3   -4   -5   -6      
-2   -2    0    1    0   -1   -2   -3   -4     
 -3   -3   -1    0    0   -1   -2   -1   -2      
 -4   -4   -2    0   -1   -1   -2   -2    0      
 -5   -5   -3   -1    1    0   -1   -2   -1       
-6   -4   -4   -2    0    0    1    0   -1      
-7   -5   -5   -3   -1   -1    1    0   -1       
-8   -6   -4   -4   -2   -2    0    0    1      
 -9   -7   -5   -5   -3   -3   -1    1    0      
-10   -8   -6   -6   -4   -4   -2    0    0 

但我希望它看起来像这样（我只关心图片中的数字）：

我必须应用惩罚：每不匹配 1 和每间隙 2，如果匹配 0。

Answer 1

有几处需要修改：

1. 请注意，在您给我们的图像中，对齐方式是从右下角到左上角。所以在这张图片中，他们并没有真正对齐AACAGTTACC 和TAAGGTCA，而是CCATTGACAA 和ACTGGAAT。
2. 你说你想要一个global alignment，但你实际上计算了一个local alignment。主要区别在于序列开始时的处罚。在全局对齐中，您必须计算第一行和第一列的插入和删除。
3. 第三，你没有正确应用你提到的惩罚。取而代之的是，您总是以 -1 惩罚并以 +1 奖励。
4. 在示例图像中，它们不是在每个位置上取最大值，而是取最小值（这是因为你的惩罚是正数，而奖励是 0，而不是相反，所以你想最小化这些值）。

完整的解决方案是：

// Note that these sequences are reversed!
String sequenceA ="CCATTGACAA";
String sequenceB = "ACTGGAAT";

// The penalties to apply
int gap = 2, substitution = 1, match = 0;

int[][] opt = new int[sequenceA.length() + 1][sequenceB.length() + 1];

// First of all, compute insertions and deletions at 1st row/column
for (int i = 1; i <= sequenceA.length(); i++)
    opt[i][0] = opt[i - 1][0] + gap;
for (int j = 1; j <= sequenceB.length(); j++)
    opt[0][j] = opt[0][j - 1] + gap;

for (int i = 1; i <= sequenceA.length(); i++) {
    for (int j = 1; j <= sequenceB.length(); j++) {
        int scoreDiag = opt[i - 1][j - 1] +
                (sequenceA.charAt(i-1) == sequenceB.charAt(j-1) ?
                    match : // same symbol
                    substitution); // different symbol
        int scoreLeft = opt[i][j - 1] + gap; // insertion
        int scoreUp = opt[i - 1][j] + gap; // deletion
        // we take the minimum
        opt[i][j] = Math.min(Math.min(scoreDiag, scoreLeft), scoreUp);
    }
}

for (int i = 0; i <= sequenceA.length(); i++) {
    for (int j = 0; j <= sequenceB.length(); j++)
        System.out.print(opt[i][j] + "\t");
    System.out.println();
}

结果和你给我们的例子一样（但是相反，记住！）：

0   2   4   6   8   10  12  14  16  
2   1   2   4   6   8   10  12  14  
4   3   1   3   5   7   9   11  13  
6   4   3   2   4   6   7   9   11  
8   6   5   3   3   5   7   8   9   
10  8   7   5   4   4   6   8   8   
12  10  9   7   5   4   5   7   9   
14  12  11  9   7   6   4   5   7   
16  14  12  11  9   8   6   5   6   
18  16  14  13  11  10  8   6   6   
20  18  16  15  13  12  10  8   7

所以最终的对齐分数位于opt[sequenceA.length()][sequenceB.length()] (7)。如果您确实需要显示图像中的反转矩阵，请执行以下操作：

for (int i = sequenceA.length(); i >=0; i--) {
    for (int j = sequenceB.length(); j >= 0 ; j--)
        System.out.print(opt[i][j] + "\t");
    System.out.println();
}

Answer 2

看看http://en.wikipedia.org/wiki/Longest_common_substring，代码几乎是几种语言的复制粘贴，并且很容易适应并告诉您对齐索引。我不得不做类似的事情并最终得到https://github.com/Pomax/DOM-diff/blob/rewrite/rewrite/rewrite.html#L103

（它返回的 SubsetMapping 基本上是一个简单的结构，它为两个上下文提供索引，https://github.com/Pomax/DOM-diff/blob/rewrite/rewrite/rewrite.html#L52）