【ROSALIND】【练Python，学生信】45 两条不等长序列间的Levenshtein距离

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

题目：

编辑距离（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距离算法即可，代码也非常简单。

代码

import numpy as np

def readfasta(lines):
    """读入fasta格式文件的函数"""
    seq = []
    index = []
    seqplast = ""
    numlines = 0
    for i in lines:
        if '>' in i:
            index.append(i.replace("\n", "").replace(">", ""))
            seq.append(seqplast.replace("\n", ""))
            seqplast = ""
            numlines += 1
        else:
            seqplast = seqplast + i.replace("\n", "")
            numlines += 1
        if numlines == len(lines):
            seq.append(seqplast.replace("\n", ""))
    seq = seq[1:]
    return index, seq


def Levenshtein(s1, s2):
    """计算Levenshtein距离的函数"""
    ld = np.zeros((len(s1)+1, len(s2)+1)) # 用zeros方法建立一个矩阵
    for i in range(1,len(s1)+1): # 填充第一列
        ld[i][0] = i
    for i in range(1,len(s2)+1): # 填充第一行
        ld[0][i] = i

    for i in range(1,len(s1)+1):
        for j in range(1,len(s2)+1):
            if s1[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)

    return ld[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))