计算DNA序列python的编辑距离答案

【问题标题】：Computing edit distance of DNA sequence python计算DNA序列python的编辑距离
【发布时间】：2017-02-21 20:03:52
【问题描述】：

因此，我的任务是在 2 个 DNA 序列之间对齐最低成本。失败的输入之一是：

CGCAATTCTGAAGCGCTGGGGAAGACGGGT & TATCCCATCGAACGCCTATTCTAGGAT

正确对齐的成本为 24，但我得到的成本为 23。

我必须阅读转换基地的成本，比如 A -> T、G -> C 等等……我得到的成本文件是：

*,-,A,T,G,C
-,0,1,2,1,3
A,1,0,1,5,1
T,2,1,0,9,1
G,1,5,9,0,1
C,3,1,1,1,0

我制作了一个 python 字典，可以完美地链接所有基础，如下所示：

{'AT': '1', '-C': '3', 'TG': '9', '-G': '1', 'AC': '1', 'C-': '3', 'CA': '1', 'AA': '0', '--': '0', 'TA': '1', 'T-': '2', 'CG': '1', '-T': '2', 'CC': '0', 'GG': '0', 'A-': '1', 'CT': '1', 'AG': '5', 'GC': '1', 'GT': '9', 'GA': '5', 'G-': '1', '-A': '1', 'TC': '1', 'TT': '0'}

目前我的实现适用于某些情况，而在其他情况下我偏离了 +/- 1。

这是我的代码的 sn-p：

def align(one_Line, costBook):
    Split_Line = one_Line.split(",")
    array = [[0] * len(Split_Line[1]) for i in Split_Line[0]]  # Zero fill array,

    xLine = Split_Line[0]
    yLine = Split_Line[1]

    for i in range(1, len(xLine)):
        array[i][0] = array[i - 1][0] + int(costBook['-' + xLine[i - 1]])
    for i in range(1, len(yLine)):
        array[0][i] = array[0][i - 1] + int(costBook[yLine[i - 1] + '-'])

    for i in range(1, len(xLine)):
        for j in range(1, len(yLine)):
            array[i][j] = min(array[i - 1][j] + diff('-', xLine[i], costBook), 
                              array[i][j - 1] + diff(yLine[j], '-', costBook), 
                              array[i - 1][j - 1] + diff(yLine[j], xLine[i], costBook))

diff函数在哪里：

def diff(x, y, cost):
    if x == y:
        return 0
    keyStr = x + y
    return int(costBook[keyStr])

我哪里出错了？是在实际的数组中填充本身，还是我做错了基本情况？

编辑：这是一个半工作版本，至少编辑成本是正确的：

AGTTGTGAAAGAACAAGCGCACAATATTGCCGCGCCGAAAGCT,TTCTTTCATTATTCAA‌ATGTATAGTTTAGAGCGTTA‌A

【问题讨论】：

这应该是 Needleman-Wunsch 算法吗？
另外，你怎么知道正确的对齐成本应该是多少？
类似但不完全是，每个序列的长度可以是任何值。匹配或不匹配或差距的成本可能会有所不同。哦，并给出了适当的费用

标签： python dynamic-programming bioinformatics dna-sequence

【解决方案1】：

我认为您的算法失败的原因是您没有考虑数组中的“间隙”行（评分矩阵）。

考虑两个序列A 和B，每个序列的长度为n。如果我们查看Needleman-Wunsch 和Smith-Waterman 算法的维基百科文章，我们会看到在它们各自的“评分”矩阵中，开头有一个额外的行。这表示A 或B 的第一个字符与空格配对。我建议您快速查看这些页面以了解我的意思

当您将数组定义为：

array = [[0] * len(Split_Line[1]) for i in Split_Line[0]]

你没有包括这个额外的行。

您需要修改align 函数以添加此额外行，并更改计算分数的逻辑。即：

def align(one_Line, costBook):
    Split_Line = one_Line.split(",")

    # Defining an array with an extra row and col to represent a leading gap
    array = [[0] * (len(Split_Line[1])+1) for i in range(len(Split_Line[0])+1)]  # Zero fill array,

    xLine = Split_Line[0]
    yLine = Split_Line[1]

    # Changed so we include our extra line in the loop
    for i in range(1, len(xLine)+1):
        array[i][0] = array[i - 1][0] + int(costBook['-' + xLine[i - 1]])
    for i in range(1, len(yLine)+1):
        array[0][i] = array[0][i - 1] + int(costBook[yLine[i - 1] + '-'])

    # Changed so we include our extra row/col in the loop
    for i in range(1, len(xLine)+1):
        for j in range(1, len(yLine)+1):
            # The references to the original string now need -1 (i.e. i-1)
            array[i][j] = min(array[i - 1][j] + diff('-', xLine[i-1], costBook), 
                              array[i][j - 1] + diff(yLine[j-1], '-', costBook), 
                              array[i - 1][j - 1] + diff(yLine[j-1], xLine[i-1], costBook))
    return array

【讨论】：