【ROSALIND】【練Python,學(xué)生信】57 使編輯距離最小的序列比對

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

題目:
編輯距離(Edit Distance Alignment)
?
Given: Two protein strings s and t in FASTA format (with each string having length at most 1000 aa).
所給:以fasta格式給出兩條蛋白質(zhì)序列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’。
?
測試數(shù)據(jù)
>Rosalind_43
PRETTY
>Rosalind_97
PRTTEIN
測試輸出
4
PRETTY--
PR-TTEIN
?
生物學(xué)背景
? ? ? ? 在06 DNA序列Hamming距離的計算中,我們知道了如何在序列只發(fā)生點突變時計算兩序列間的Hamming距離。而在45 兩條不等長序列間的Levenshtein距離中,我們知道了如何在考慮點突變、插入、刪除的情況下計算兩條不等長序列的Levenshtein距離。
? ? ? ??在考慮插入和刪除時,我們需要引入空位(‘-’)作為占位符,幫助我們實現(xiàn)不等長序列的比對。假設(shè)有兩個序列s和t,在添加空位進行比對后,這兩個序列會分別變成序列s’和t’,s’和t’長度一定相同,并且相同位上肯定不都為空位,即假如s’[j]是空位,t’[j]肯定不為空位。
? ? ? ??因為s’和t’已具有相同的長度,我們可以把這兩個序列寫在一起,直接一一比對每一位的變化。比如s’[j]和t’[j]不一致,相當(dāng)于由s到t此位上出現(xiàn)一個點突變;s’[j]為空位,說明由s到t出現(xiàn)了一個插入;t’[j]為空位說明,由s到t出現(xiàn)了一個刪除。
? ? ? ??我們可以把由s’到t’經(jīng)過的變化數(shù)量定義為dH(s’, t’),則dE(s,t)=min dH(s’, t’),我們稱這個使編輯數(shù)量最小的比對方法為最佳聯(lián)配(optimal alignment)。
??
思路
? ? ? ??我是在24 求解最長遞增子序列代碼的基礎(chǔ)上寫的這道題。在此再做簡單解釋。
? ? ? ??動態(tài)規(guī)劃的整個過程可以劃為四個部分:
? ? ? ??1)存儲子問題的最優(yōu)化的動態(tài)規(guī)劃矩陣;
? ? ? ??2)最優(yōu)化的遞歸計算方法;
? ? ? ??3)給出子問題最優(yōu)解的矩陣填充過程;
?? ? ? ?4)尋找最優(yōu)化比對路徑的回溯方法。?
? ? ? ??以測試數(shù)據(jù)為例:
1)? 存儲子問題的最優(yōu)化的動態(tài)規(guī)劃矩陣;
[['', 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)? 最優(yōu)化的遞歸計算方法;
? ? ? ??這里的核心思路是在序列中找一個匹配方法使一個序列變成另一個序列需要增加的編輯次數(shù)最少,每一步都如此考慮,由此形成遞歸。具體方法詳見代碼注釋。
3)? 給出子問題最優(yōu)解的矩陣填充過程;
? ? ? ??每一步都記下使編輯次數(shù)增加最小的填充方法以及編輯次數(shù),如下表。
[['', 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)? 尋找最優(yōu)化比對路徑的回溯方法。
[['', 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
? ? ? ??編輯次數(shù)為4
注意:此處答案與范例給的序列不一致,但因為編輯次數(shù)仍為4,所以能夠通過??梢钥吹剑舅惴ㄗ非蟮牟⒉皇潜葘ι系淖址麛?shù)量最多,而僅能使編輯次數(shù)最小。如果希望得到最多的比對字符,稍微更改下罰分方式即可。
?
代碼