python不等长编码有哪些_【ROSALIND】【练Python，学生信】45 两条不等长序列间的Levenshtein距离...

题目：

编辑距离（Edit Distance）

Given: Two protein strings s and t in FASTA format (each of length at most 1000 aa).

所给：两个蛋白质序列s和t，长度都不超过1000个氨基酸。

Return: The edit distance dE(s,t).

需得：两条序列的编辑距离。

测试数据

>Rosalind_39

PLEASANTLY

>Rosalind_11

MEANLY

测试输出

5

生物学背景

在06 DNA序列Hamming距离的计算中，我们已经能计算相同长度的两个同源序列，一个序列变成另一个需要的最少点突变次数。但在现实生活中，序列间除了有点突变，还会发生插入和删除，使两条序列长度不再相同，面对这种情况，我们需要把之前的方法进行推广，使其可以计算插入和删除的情况。

数学背景

编辑距离是由俄国科学家Vladimir Levenshtein提出的，所以也叫Levenshtein距离。这是度量两个序列间相似度的指标，简单来说就是一个序列变成另一个序列所需的最少单字符编辑次数。

有关Levenshtein距离的介绍可以详见这篇博文（https://blog.csdn.net/asty9000/article/details/81384650）。简单来说，对于a、b两个字符串，a的第i个字符ai和b中的第j个字符bj有如下两个可能：1）ai和bj相同，不增加a、b两串的距离2）ai和bj不同，则a、b两串间的距离要加一，此时我们又面临两个选择：是把a1a2…a(i-1)变成b1b2…bj呢？还是把a1a2…ai变成b1b2…b(j-1)呢？这就要看哪种变法距离更小了，选择距离小的那种方法加一即可。对a、b中每个字符都做这种考虑，就形成了一个矩阵，把右下角的数字取出来就是两个字符串的Levenshtein距离了。

Python背景

Numpy是Python一个库，常用于科学计算。用Numpy中的zeros方法可以快速建立一个指定形状和类型的数组，用“0”填充。基础用法为numpy.zeros(shape，dtype=float，order = 'C')。

思路

本题的求解方法有些类似38 求解最长公共子序列中用到的动态规划算法。我在刚开始时也希望能借用最长公共子序列快速得到结果，不过未能如愿，反而越想越复杂。实际上直接用Levenshtein距离算法即可，代码也非常简单。

代码

importnumpyasnp

defreadfasta(lines):"""读入fasta格式文件的函数"""seq = []

index = []

seqplast =""numlines =0foriinlines:if'>'ini:

index.append(i.replace("\n","").replace(">",""))

seq.append(seqplast.replace("\n",""))

seqplast =""numlines +=1else:

seqplast = seqplast + i.replace("\n","")

numlines +=1ifnumlines ==len(lines):

seq.append(seqplast.replace("\n",""))

seq = seq[1:]returnindex, seqdefLevenshtein(s1, s2):"""计算Levenshtein距离的函数"""ld = np.zeros((len(s1)+1,len(s2)+1))#用zeros方法建立一个矩阵foriinrange(1,len(s1)+1):#填充第一列ld[i][0] = iforiinrange(1,len(s2)+1):#填充第一行ld[0][i] = iforiinrange(1,len(s1)+1):forjinrange(1,len(s2)+1):ifs1[i-1] == s2[j -1]:#如果两个字符相同tep =0#距离不增加else:

tep =1

#字符不同则在点突变、插入、删除中选使距离最短的ld[i][j] =min(ld[i-1][j-1] + tep, ld[i][j-1] +1, ld[i-1][j] +1)

returnld[i][j]#返回最短距离f =open('input.txt','r')

lines = f.readlines()

f.close()

[index, seq] = readfasta(lines)

seq1 = seq[0]

seq2 = seq[1]

res = []

res = Levenshtein(seq1, seq2)print(int(res))