利用Python实现Smithwaterman算法

def smith_waterman(seq1, seq2, match_score=3, mismatch_penalty=-1, gap_penalty=-2):
# 初始化矩阵
matrix = [[0] * (len(seq2) + 1) for _ in range(len(seq1) + 1)]
max_score = 0
max_i, max_j = 0, 0

# 填充矩阵
for i in range(1, len(seq1) + 1):
    for j in range(1, len(seq2) + 1):
        if seq1[i - 1] == seq2[j - 1]:
            score = matrix[i - 1][j - 1] + match_score
        else:
            score = max(
                matrix[i - 1][j] + gap_penalty,
                matrix[i][j - 1] + gap_penalty,
                matrix[i - 1][j - 1] + mismatch_penalty
            )
        if score > max_score:
            max_score = score
            max_i, max_j = i, j
        matrix[i][j] = max(score, 0)

# 回溯，获取最大得分路径
align1, align2 = '', ''
i, j = max_i, max_j
while i > 0 and j > 0 and matrix[i][j] != 0:
    if matrix[i][j] == matrix[i - 1][j - 1] + (match_score if seq1[i - 1] == seq2[j - 1] else mismatch_penalty):
        align1 = seq1[i - 1] + align1
        align2 = seq2[j - 1] + align2
        i -= 1
        j -= 1
    elif matrix[i][j] == matrix[i - 1][j] + gap_penalty:
        align1 = seq1[i - 1] + align1
        align2 = '-' + align2
        i -= 1
    else:
        align1 = '-' + align1
        align2 = seq2[j - 1] + align2
        j -= 1

return max_score, align1, align2

示例用法

seq1 = “AGCACACA”
seq2 = “ACACACTA”
score, align1, align2 = smith_waterman(seq1, seq2)
print(“Alignment score:”, score)
print(“Sequence 1:”, align1)
print(“Sequence 2:”, align2)