【问题标题】:global alignment sequence function全局比对序列函数
【发布时间】:2018-12-14 17:42:30
【问题描述】:

我正在尝试实现the Needleman-Wunsch algorithm 以获得全局对齐函数中的最低分数,但是当两个序列相等时,我得到的不是最低分数 0,而是 8。

这段代码有什么问题?

alphabet = ["A", "C", "G", "T"] 
score = [[0, 4, 2, 4, 8], \
     [4, 0, 4, 2, 8], \
     [2, 4, 0, 4, 8], \
     [4, 2, 4, 0, 8], \
     [8, 8, 8, 8, 8]]

def globalAlignment(x, y):
#Dynamic version very fast
    D = []
    for i in range(len(x)+1):
        D.append([0]* (len(y)+1))

    for i in range(1, len(x)+1):
        D[i][0] = D[i-1][0] + score[alphabet.index(x[i-1])][-1]
    for i in range(len(y)+1):
        D[0][i] = D[0][i-1]+ score[-1][alphabet.index(y[i-1])]

    for i in range(1, len(x)+1):
        for j in range(1, len(y)+1):
            distHor = D[i][j-1]+ score[-1][alphabet.index(y[j-1])]
            distVer = D[i-1][j]+ score[-1][alphabet.index(x[i-1])]
            if x[i-1] == y[j-1]:
                distDiag = D[i-1][j-1]
            else:
                distDiag = D[i-1][j-1] + score[alphabet.index(x[i-1])][alphabet.index(y[j-1])]

            D[i][j] = min(distHor, distVer, distDiag)

    return D[-1][-1]

x = "ACGTGATGCTAGCAT"
y = "ACGTGATGCTAGCAT"
print(globalAlignment(x, y))

【问题讨论】:

  • 欢迎来到 StackOverflow。请按照您创建此帐户时的建议阅读并遵循帮助文档中的发布指南。 On topic、how to ask 和 ... the perfect question 在此处申请。
  • 具体来说,你没能“让别人轻松帮助你”。你有某种“神奇”的距离计算,你仍然对我们隐藏。一个字母的变量名称,令人眼花缭乱的下标序列——没有文档、解释或调试跟踪。 为什么您希望这些计算为 0? 它是如何 转到8 的?失败的中间步骤是什么?
  • 查看这个可爱的debug 博客寻求帮助。
  • 我从 Ben Langmead 在 coursera 中为我们带来的 DNA 测序算法课程中获得了这段代码。代码在他的机器上运行正常,但在我自己的机器上却不能得到同样的结果。

标签: python dynamic-programming bioinformatics sequence-alignment


【解决方案1】:

只是改变 -> for i in range(len(y)+1): 到 for i in range(1, len(y) + 1): 和 -> distVer = D[i-1][j]+ score[-1][alphabet.index(x[i-1])] 到

distVer = D[i - 1][j] + score[alphabet.index(x[i - 1])][-1]

【讨论】:

    【解决方案2】:

    至少

    distHor = D[i][j-1]+ score[-1][alphabet.index(y[j-1])]
    distVer = D[i-1][j]+ score[-1][alphabet.index(x[i-1])]
    

    是可疑的,因为您没有在初始化中为 [-1] 使用相同的位置, 并且两个距离不太可能在权重中使用相同的方向...
    我想应该是

    score[alphabet.index(x[i-1])][-1]
    

    但它可能不是唯一的错误......

    【讨论】:

    • 嗨,我通过在 socre score = [[0, 4, 2, 4, 8], \ [4, 0, 4, 2 , 8], \ [2, 4, 0, 4, 8], \ [4, 2, 4, 0, 8], \ [0, 0, 0, 0, 0]]
    【解决方案3】:

    我通过在最后一个分数列表中输入 0 而不是 8 来解决问题;

    alphabet = ["A", "C", "G", "T"] 
    score = [[0, 4, 2, 4, 8], \
         [4, 0, 4, 2, 8], \
         [2, 4, 0, 4, 8], \
         [4, 2, 4, 0, 8], \
         [0, 0, 0, 0, 0]]
    
    def globalAlignment(x, y):
    #Dynamic version very fast
    D = []
    for i in range(len(x)+1):
        D.append([0]* (len(y)+1))
    
    for i in range(1, len(x)+1):
        D[i][0] = D[i-1][0] + score[alphabet.index(x[i-1])][-1]
    for i in range(len(y)+1):
        D[0][i] = D[0][i-1]+ score[-1][alphabet.index(y[i-1])]
    
    for i in range(1, len(x)+1):
        for j in range(1, len(y)+1):
            distHor = D[i][j-1]+ score[-1][alphabet.index(y[j-1])]
            distVer = D[i-1][j]+ score[alphabet.index(x[i-1])][-1]
            if x[i-1] == y[j-1]:
                distDiag = D[i-1][j-1]
            else:
                distDiag = D[i-1][j-1] + score[alphabet.index(x[i-1])][alphabet.index(y[j-1])]
    
            D[i][j] = min(distHor, distVer, distDiag)
    
    return D[-1][-1]    
    

    【讨论】:

      猜你喜欢
      • 1970-01-01
      • 1970-01-01
      • 2022-01-23
      • 2012-10-15
      • 2017-03-05
      • 2018-04-30
      • 2012-12-16
      • 1970-01-01
      • 1970-01-01
      相关资源
      最近更新 更多