DNA比对 - 得分是否是附加的？

Question

所以我有一个递归代码,可以为2条DNA链提供最佳的对齐,但问题是它的执行速度非常慢(我需要它递归).然后我在麻省理工学院网站上看到结果是附加的,这对我很好,但后来我想了一下,我发现有一个问题:

网站:http://ocw.mit.edu/courses/electrical-engineering-and-computer-science/6-096-algorithms-for-computational-biology-spring-2005/lecture-notes/lecture5_newest.pdf

麻省理工学院网站说,对于给定的溢出(i,j):first_strand(0,i)和second_strand(0,j)alignment
+ first_strand(i,len)和second_strand(j,len)对齐

等于

first_strand和second strand alignment

但:

GTC GTAA

具有GTA比对的G是G--并且具有A比对的GTA TC是TC并且A-结果= G-TC和GTAA-

真正的最佳结果= GTC-GTAA

谁能解释他们在麻省理工学院网站上的意思？我可能错了!

Answer 1

我假设你在谈论这个链接.

如果是这样,请仔细阅读数百次;-)它是"附加的",因为您只考虑将分割固定在特定(i, j)对上的对齐.

在您假设的反例中,您首先打破了初始G关闭GTC和初始GTA关闭GTAA.然后G--是改变最短的方式GTC进入G.精细.继续使用相同的分割,您仍然需要对齐剩余的右手部分: TCwith A.还好.

这并不是说这是最好的分裂.只有声称它是最好的对齐,因为你只考虑特定的分裂.

这是动态编程方法中的一小步,这是您缺少的部分.仍需要计算所有可能拆分的最佳对齐方式.

动态编程起初很棘手.你不应该期望从盯着电报幻灯片中学到它.阅读一本真正的教科书,或在网上搜索教程.

加快递归版本

注释表明这个"必须"的代码是递归的.那好吧 ;-)

注意:我只是把它放在一起来说明加速合适的递归函数的一般过程.它几乎没有经过测试.

首先是一个完全天真的递归版本:

def lev(a, b):
    if not a:
        return len(b)
    if not b:
        return len(a)
    return min(lev(a[:-1], b[:-1]) + (a[-1] != b[-1]),
               lev(a[:-1], b) + 1,
               lev(a, b[:-1]) + 1)

我将在这里讨论的所有运行中使用"absd31-km"和"ldk3-1fjm"作为参数.

在我的盒子上,使用Python 3,这个简单的函数在大约1.6秒后返回7.它非常缓慢.

最明显的问题是无休止地重复的字符串切片.:索引中的每个都需要与被切片的字符串的当前长度成比例的时间.所以第一个改进是传递字符串索引.由于代码总是切掉字符串的前缀,我们只需要传递"字符串结尾"索引:

def lev2(a, b):
    def inner(j1, j2):
        if j1 < 0:
            return j2 + 1
        if j2 < 0:
            return j1 + 1
        return min(inner(j1-1, j2-1) + (a[j1] != b[j2]),
                   inner(j1-1, j2) + 1,
                   inner(j1, j2-1) + 1)
    return inner(len(a)-1, len(b)-1)

好多了!此版本在"仅"中返回7,大约1.44秒.仍然非常缓慢,但比原来更好.它的优点是会在较长的琴弦上增加,但谁在乎;-)

我们差不多完成了!现在要注意的重要一点是该函数在运行过程中多次传递相同的参数.我们在"备忘录"中捕获那些以避免所有冗余计算:

def lev3(a, b):
    memo = {}
    def inner(j1, j2):
        if j1 < 0:
            return j2 + 1
        if j2 < 0:
            return j1 + 1
        args = j1, j2
        if args in memo:
            return memo[args]
        result = min(inner(j1-1, j2-1) + (a[j1] != b[j2]),
                     inner(j1-1, j2) + 1,
                     inner(j1, j2-1) + 1)
        memo[args] = result
        return result
    return inner(len(a)-1, len(b)-1)

该版本在大约0.00026秒内返回7,比lev2它快5000倍.

现在,如果您已经研究了基于矩阵的算法,并且稍微眯一眼,您将看到lev3() 有效地构建二维矩阵映射索引对以在其memo字典中产生结果.除了递归版本以更复杂的方式构建矩阵之外,它们实际上是相同的.另一方面,递归版本可能更容易理解和推理.请注意,您发现的幻灯片称为memoization aporoach"自上而下",嵌套循环矩阵称为"自下而上".这些都很好描述.

你还没有说过你的递归函数是如何工作的,但如果它遭受任何类似的递归过量,你应该能够使用类似的技术获得类似的加速:-)