重叠比对（易错读段重叠）：现有代码故障排查与修复

问题定义 #

寻找两个字符串的最高得分重叠比对（Overlap Alignment）。

• 输入：匹配得分m、错配罚分μ、空位罚分σ，以及两个DNA字符串s和t。
• 输出：s与t的重叠比对最大得分，以及达到该得分的重叠比对结果。

背景说明 #

生物学家利用重叠读段组装基因组，但读段中的错误增加了该问题的复杂度。为找到易错读段间的重叠，我们定义字符串s = s₁…sₙ与t = t₁…tₘ的重叠比对为s的一个后缀与t的一个前缀的全局比对。最优重叠比对是在所有i和j中，最大化s的i后缀（sᵢ…sₙ）与t的j前缀（t₁…tⱼ）的全局比对得分。

输入输出格式 #

输入格式 #

• 第一行：m、μ、σ（空格分隔）
• 第二行：DNA字符串s
• 第三行：DNA字符串t

输出格式 #

• 第一行：重叠比对的最大得分
• 第二行：插入适当空位的s的后缀
• 第三行：插入适当空位的t的前缀

代码现状 #

以下Python代码可通过部分测试用例，但在部分测试用例中失败：

def align(m, mu, sigma, s, t):
    # initialize the scoring matrix
    score = [[0 for _ in range(len(t) + 1)] for _ in range(len(s) + 1)]
    backtrack = [[0 for _ in range(len(t) + 1)] for _ in range(len(s) + 1)]

    for i in range(1, len(s) + 1):
        backtrack[i][0] = "↓"
    for j in range(1, len(t) + 1):
        score[0][j] = score[0][j - 1] - sigma
        backtrack[0][j] = "→"

    # fill in the scoring matrix
    for i in range(1, len(s) + 1):
        for j in range(1, len(t) + 1):
            if s[i - 1] == t[j - 1]:
                match = score[i - 1][j - 1] + m
            else:
                match = score[i - 1][j - 1] - mu
            delete = score[i - 1][j] - sigma
            insert = score[i][j - 1] - sigma
            score[i][j] = max(match, delete, insert)
            if score[i][j] == match:
                backtrack[i][j] = "↘"
            elif score[i][j] == delete:
                backtrack[i][j] = "↓"
            else:
                backtrack[i][j] = "→"

    # reconstruct the alignment
    max_score = max(score[-1])
    i = len(s)
    j = score[-1].index(max_score)
    aligned_s, aligned_t = "", ""
    while i > 0 and j > 0:
        if backtrack[i][j] == "↘":
            aligned_s = s[i - 1] + aligned_s
            aligned_t = t[j - 1] + aligned_t
            i -= 1
            j -= 1
        elif backtrack[i][j] == "↓":
            aligned_s = s[i - 1] + aligned_s
            aligned_t = "-" + aligned_t
            i -= 1
        else:
            aligned_s = "-" + aligned_s
            aligned_t = t[j - 1] + aligned_t
            j -= 1
    print(score[len(s)][len(t)])
    print(aligned_s)
    print(aligned_t)


if __name__ == "__main__":
    m, mu, sigma = map(int, input().strip().split())
    s, t = [input().strip() for _ in range(2)]
    align(m, mu, sigma, s, t)

通过的测试用例 #

1. 输入：

1 1 2
GAGA
GAT

输出：

2
GA
GA

2. 输入：

1 1 1
CCAT
AT

输出：

2
AT
AT

3. 输入：

1 5 1
GAT
CAT

输出：

1
-AT
CAT

4. 输入：

1 5 1
ATCACT
AT

输出：

1
ACT
A-T

失败的测试用例 #

1. 输入：

1 1 5
ATCACT
ATG

预期输出：

0
CT
AT

实际输出：

0
"" # empty string
"" # empty string

2. 输入：

3 2 1
CAGAGATGGCCG
ACG

预期输出：

5
-CG
ACG

实际输出：

5
ATGGCCG
A----CG

3. 输入：

2 3 1
CTT
AGCATAAAGCATT

预期输出：

0
--CT-T
AGC-AT

实际输出：

0
"" # empty string
"" # empty string

内容的提问来源于stack exchange，提问作者Ali Alsawad