如果第一次閱讀本系列文檔請(qǐng)先移步閱讀【ROSALIND】【練Python，學(xué)生信】00 寫在前面 ?謝謝配合~

題目：

編輯距離（Edit Distance）

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

所給：兩個(gè)蛋白質(zhì)序列s和t，長(zhǎng)度都不超過1000個(gè)氨基酸。

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

需得：兩條序列的編輯距離。

?

測(cè)試數(shù)據(jù)

>Rosalind_39

PLEASANTLY

>Rosalind_11

MEANLY

測(cè)試輸出

5

?

生物學(xué)背景

? ? ? ?在06 DNA序列Hamming距離的計(jì)算中，我們已經(jīng)能計(jì)算相同長(zhǎng)度的兩個(gè)同源序列，一個(gè)序列變成另一個(gè)需要的最少點(diǎn)突變次數(shù)。但在現(xiàn)實(shí)生活中，序列間除了有點(diǎn)突變，還會(huì)發(fā)生插入和刪除，使兩條序列長(zhǎng)度不再相同，面對(duì)這種情況，我們需要把之前的方法進(jìn)行推廣，使其可以計(jì)算插入和刪除的情況。

?

數(shù)學(xué)背景

? ? ? ?編輯距離是由俄國(guó)科學(xué)家Vladimir Levenshtein提出的，所以也叫Levenshtein距離。這是度量?jī)蓚€(gè)序列間相似度的指標(biāo)，簡(jiǎn)單來(lái)說就是一個(gè)序列變成另一個(gè)序列所需的最少單字符編輯次數(shù)。

? ? ? ?有關(guān)Levenshtein距離的介紹可以詳見這篇博文（https://blog.csdn.net/asty9000/article/details/81384650）。簡(jiǎn)單來(lái)說，對(duì)于a、b兩個(gè)字符串，a的第i個(gè)字符ai和b中的第j個(gè)字符bj有如下兩個(gè)可能：1）ai和bj相同，不增加a、b兩串的距離2）ai和bj不同，則a、b兩串間的距離要加一，此時(shí)我們又面臨兩個(gè)選擇：是把a(bǔ)1a2…a(i-1)變成b1b2…bj呢？還是把a(bǔ)1a2…ai變成b1b2…b(j-1)呢？這就要看哪種變法距離更小了，選擇距離小的那種方法加一即可。對(duì)a、b中每個(gè)字符都做這種考慮，就形成了一個(gè)矩陣，把右下角的數(shù)字取出來(lái)就是兩個(gè)字符串的Levenshtein距離了。

?

Python背景

? ? ? ? Numpy是Python一個(gè)庫(kù)，常用于科學(xué)計(jì)算。用Numpy中的zeros方法可以快速建立一個(gè)指定形狀和類型的數(shù)組，用“0”填充?；A(chǔ)用法為numpy.zeros(shape，dtype=float，order = 'C')。

?

思路

? ? ? ?本題的求解方法有些類似38 求解最長(zhǎng)公共子序列中用到的動(dòng)態(tài)規(guī)劃算法。我在剛開始時(shí)也希望能借用最長(zhǎng)公共子序列快速得到結(jié)果，不過未能如愿，反而越想越復(fù)雜。實(shí)際上直接用Levenshtein距離算法即可，代碼也非常簡(jiǎn)單。

?

?

代碼

import numpy as np

?

def readfasta(lines):
??? """讀入fasta格式文件的函數(shù)"""
??? 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):
??? """計(jì)算Levenshtein距離的函數(shù)"""
??? ld = np.zeros((len(s1)+1, len(s2)+1)) # 用zeros方法建立一個(gè)矩陣
??? 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]: # 如果兩個(gè)字符相同
??????????????? tep = 0 # 距離不增加
??????????? else:
??????????????? tep = 1

# 字符不同則在點(diǎn)突變、插入、刪除中選使距離最短的
??????????? 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))