【ROSALIND】【练Python,学生信】57 使编辑距离最小的序列比对

如果第一次阅读本系列文档请先移步阅读【ROSALIND】【练Python,学生信】00 写在前面 谢谢配合~

题目:
编辑距离(Edit Distance Alignment)
Given: Two protein strings s and t in FASTA format (with each string having length at most 1000 aa).
所给:以fasta格式给出两条蛋白质序列s和t,长度均不超过1000个氨基酸。
Return: The edit distance dE(s,t) followed by two augmented strings s’ and t’ representing an optimal alignment of s and t.
需得:编辑距离dE(s,t),以及s和t比对后得到的s’和t’。
测试数据
>Rosalind_43
PRETTY
>Rosalind_97
PRTTEIN
测试输出
4
PRETTY--
PR-TTEIN
生物学背景
在06 DNA序列Hamming距离的计算中,我们知道了如何在序列只发生点突变时计算两序列间的Hamming距离。而在45 两条不等长序列间的Levenshtein距离中,我们知道了如何在考虑点突变、插入、删除的情况下计算两条不等长序列的Levenshtein距离。
在考虑插入和删除时,我们需要引入空位(‘-’)作为占位符,帮助我们实现不等长序列的比对。假设有两个序列s和t,在添加空位进行比对后,这两个序列会分别变成序列s’和t’,s’和t’长度一定相同,并且相同位上肯定不都为空位,即假如s’[j]是空位,t’[j]肯定不为空位。
因为s’和t’已具有相同的长度,我们可以把这两个序列写在一起,直接一一比对每一位的变化。比如s’[j]和t’[j]不一致,相当于由s到t此位上出现一个点突变;s’[j]为空位,说明由s到t出现了一个插入;t’[j]为空位说明,由s到t出现了一个删除。
我们可以把由s’到t’经过的变化数量定义为dH(s’, t’),则dE(s,t)=min dH(s’, t’),我们称这个使编辑数量最小的比对方法为最佳联配(optimal alignment)。
思路
我是在24 求解最长递增子序列代码的基础上写的这道题。在此再做简单解释。
动态规划的整个过程可以划为四个部分:
1)存储子问题的最优化的动态规划矩阵;
2)最优化的递归计算方法;
3)给出子问题最优解的矩阵填充过程;
4)寻找最优化比对路径的回溯方法。
以测试数据为例:
1) 存储子问题的最优化的动态规划矩阵;
[['', 0], ['-', 0], ['P', 0], ['R', 0], ['T', 0], ['T', 0], ['E', 0], ['I', 0], ['N', 0]]
[['-', 0], ['', 0], ['←', 1], ['←', 2], ['←', 3], ['←', 4], ['←', 5], ['←', 6], ['←', 7]]
[['P', 0], ['↑', 1], ['', 0], ['', 0], ['', 0], ['', 0], ['', 0], ['', 0], ['', 0]]
[['R', 0], ['↑', 2], ['', 0], ['', 0], ['', 0], ['', 0], ['', 0], ['', 0], ['', 0]]
[['E', 0], ['↑', 3], ['', 0], ['', 0], ['', 0], ['', 0], ['', 0], ['', 0], ['', 0]]
[['T', 0], ['↑', 4], ['', 0], ['', 0], ['', 0], ['', 0], ['', 0], ['', 0], ['', 0]]
[['T', 0], ['↑', 5], ['', 0], ['', 0], ['', 0], ['', 0], ['', 0], ['', 0], ['', 0]]
[['Y', 0], ['↑', 6], ['', 0], ['', 0], ['', 0], ['', 0], ['', 0], ['', 0], ['', 0]]
2) 最优化的递归计算方法;
这里的核心思路是在序列中找一个匹配方法使一个序列变成另一个序列需要增加的编辑次数最少,每一步都如此考虑,由此形成递归。具体方法详见代码注释。
3) 给出子问题最优解的矩阵填充过程;
每一步都记下使编辑次数增加最小的填充方法以及编辑次数,如下表。
[['', 0], ['-', 0], ['P', 0], ['R', 0], ['T', 0], ['T', 0], ['E', 0], ['I', 0], ['N', 0]]
[['-', 0], ['', 0], ['←', 1], ['←', 2], ['←', 3], ['←', 4], ['←', 5], ['←', 6], ['←', 7]]
[['P', 0], ['↑', 1], ['↖', 0], ['←', 1], ['←', 2], ['←', 3], ['←', 4], ['←', 5], ['←', 6]]
[['R', 0], ['↑', 2], ['↑', 1], ['↖', 0], ['←', 1], ['←', 2], ['←', 3], ['←', 4], ['←', 5]]
[['E', 0], ['↑', 3], ['↑', 2], ['↑', 1], ['↖', 1], ['↖', 2], ['↖', 2], ['←', 3], ['←', 4]]
[['T', 0], ['↑', 4], ['↑', 3], ['↑', 2], ['↖', 1], ['↖', 1], ['←', 2], ['↖', 3], ['↖', 4]]
[['T', 0], ['↑', 5], ['↑', 4], ['↑', 3], ['↖', 2], ['↖', 1], ['↖', 2], ['↖', 3], ['↖', 4]]
[['Y', 0], ['↑', 6], ['↑', 5], ['↑', 4], ['↑', 3], ['↑', 2], ['↖', 2], ['↖', 3], ['↖', 4]]
4) 寻找最优化比对路径的回溯方法。
[['', 0], ['-', 0], ['P', 0], ['R', 0], ['T', 0], ['T', 0], ['E', 0], ['I', 0], ['N', 0]]
[['-', 0], ['', 0], ['←', 1], ['←', 2], ['←', 3], ['←', 4], ['←', 5], ['←', 6], ['←', 7]]
[['P', 0], ['↑', 1], ['↖', 0], ['←', 1], ['←', 2], ['←', 3], ['←', 4], ['←', 5], ['←', 6]]
[['R', 0], ['↑', 2], ['↑', 1], ['↖', 0], ['←', 1], ['←', 2], ['←', 3], ['←', 4], ['←', 5]]
[['E', 0], ['↑', 3], ['↑', 2], ['↑', 1], ['↖', 1], ['↖', 2], ['↖', 2], ['←', 3], ['←', 4]]
[['T', 0], ['↑', 4], ['↑', 3], ['↑', 2], ['↖', 1], ['↖', 1], ['←', 2], ['↖', 3], ['↖', 4]]
[['T', 0], ['↑', 5], ['↑', 4], ['↑', 3], ['↖', 2], ['↖', 1], ['↖', 2], ['↖', 3], ['↖', 4]]
[['Y', 0], ['↑', 6], ['↑', 5], ['↑', 4], ['↑', 3], ['↑', 2], ['↖', 2], ['↖', 3], ['↖', 4]]
可以看到,两个序列分别为
PRET-TY
PRTTEIN
编辑次数为4
注意:此处答案与范例给的序列不一致,但因为编辑次数仍为4,所以能够通过。可以看到,本算法追求的并不是比对上的字符数量最多,而仅能使编辑次数最小。如果希望得到最多的比对字符,稍微更改下罚分方式即可。
代码