核酸檢測經(jīng)此一疫已經(jīng)成為了人盡皆知的檢測手段，但是你可曾想過經(jīng)過儀器檢測出來的 DNA序列是如何與正常情況下的 DNA 序列做對比的呢？面對整村的核酸檢測結(jié)果，計算工作肯定需要電腦來完成，如果我們想知道一個受試者的 DNA 序列與正常情況下的 DNA 序列有至少幾個個元素差異時，我們就需要一個算法來解決，它叫 萊文斯坦 (levenshtein) 算法，而這個算法背后有一個主要的技術(shù)：動態(tài)規(guī)劃！

import?numpy?as?np

本文將只通過 numpy 從頭構(gòu)建萊文斯坦算法，如果還沒安裝這個庫的讀者們可以通過下面指令在終端安裝：

pip?install?numpy

編輯距離

DNA 序列主要可以用四種英文字母表示，ATCG，至于什么樣的檢測結(jié)果要對應(yīng)到什么樣的英文字母不在這篇文章的討論范圍，我們只討論如何處理這些序列數(shù)據(jù)。假如檢測站的樣本里我們的到了兩組 DNA 序列：

判斷這兩組序列與目標(biāo)序列之間的差異前，我們需要先定義差異的含義，也就是編輯距離，面對一個序列，我們需要經(jīng)過幾步操作才能夠使其和目標(biāo)序列變得一模一樣，如果是三步，那么編輯距離就是3。但是廣義地說，修改一個序列使其和另一個序列一樣的方法有無限多種，因此我們需要給另一個約束條件，就是如何用最少的步驟來完成修改的操作，而具體操作動作只能包含下面三種方式：

插入
刪除
替換

為了有效解決這個問題，動態(tài)規(guī)劃就是本篇文章所介紹算法的核心邏輯！

動態(tài)規(guī)劃

動態(tài)規(guī)劃的核心概念就是把一個大的問題拆分成很多個小的問題分別擊破，使得整體問題難度降低，進而更容易找出問題的最優(yōu)答案。面對 DNA 序列問題，大問題就是一串完整的序列到底有哪些不同，而小問題則是對應(yīng)到每個位置上的字符，到底是應(yīng)該采取什么樣的動作，先給出計算兩個不同序列編輯距離的計算過程。

為了把 "ATTGTCT" 改成 "GTAGCTT"，根據(jù)矩陣?yán)镒钣蚁陆堑挠嬎憬Y(jié)果，最少的操作步驟就是 4 步，至于這四個步驟具體是 插入，刪除，替換 怎么樣的排列組合，那得依據(jù)計算矩陣右下角數(shù)字的過程中小編所標(biāo)注的框的顏色和對應(yīng)到一個 2x2 的行為矩陣來判定：

如果我們把行為矩陣?yán)锏挠蚁陆俏恢脤?yīng)到大矩陣?yán)锩恳粋€紅色框，紅框的數(shù)字全是從大矩陣中橙色框 + 1 而得到的，而橙色框?qū)?yīng)了行為矩陣則是 替換 的操作！通過同樣的規(guī)律，我們也可以從藍色的路徑了解到從 "ATTGTCT" 變到 "GTA" 最少需要 5 步，從 "ATT" 變到 "GTAGCTT" 則需要 4 步，至于具體填上數(shù)字的細節(jié)，小編馬上接著說！

1. 原理細節(jié)

在最一開始什么數(shù)字都不存在在矩陣中時，只有 變化前 的序列被放在最上面，和 變化后 的序列被放在最左邊，而序列的最前面需要填上一個代表本來什么都沒有的站位符，我們可以理解為一個字符串在最易開始是什么元素都沒有的，因此不論是變化前后兩個序列最前面都得加上空的占位符。接著，我們就能開始逐一解大問題中的每一個小問題，從空的占位符變成任意長度的序列所需要的步數(shù)肯定等于序列的長度，因此第一個行列就可以理所當(dāng)然的填上 0 1 2 ... 7。

剩下來還沒填上數(shù)字的格子，我們就需要以 2x2 所排部的四個格子為單位去計算最小的步驟，而這四個格子也就剛好對應(yīng)了行為矩陣的位置，為了確定小問題的最優(yōu)操作，我們需要知道分別操作了 插入，刪除，替換 之后，最小的成本是多少。

1-1. 不同元素 - 綠色方形塊

如果對應(yīng)位置的字母需要被改變，例如圖中 A 要換成 G，至少 1 步 (替換) 需要操作，但這是因為一開始問題一目了然，能夠通過直覺看出問題的解，實際上在動態(tài)規(guī)劃的思維中，我們需要發(fā)現(xiàn)一個能夠一直被重復(fù)使用的規(guī)律，因此 A 換成 G 是 1 步的答案背后應(yīng)該要根據(jù) 2x2 范圍里的數(shù)字去比較：

0 + 1 表示替換 A 成 G 所需要的步數(shù)
1 + 1 表示插入 A 成 G 所需要的步數(shù)
1 + 1 表示刪除 A 成 G 所需要的步數(shù)

而最后根據(jù)這三個操作的結(jié)果來看，證明了最少步數(shù)是替換操作，因此這步將遵循此操作，并把 0 + 1 的結(jié)果寫進 2x2 矩陣范圍的右下角。

1-2. 相同元素 - 黃色方形塊

相反地，如果我們遇到要把 G 換成 G 的情況，就直接把 2x2 矩陣范圍的左上角數(shù)字抄到右下角，因為同樣的字符是不需要做任何變動的，也就沒有任何一個操作合適，所以最少步數(shù)的結(jié)論就還是 3 步。了解了處理 不同元素 與 相同元素 的計算規(guī)則之后，就可以逐一填上矩陣?yán)锩恳粋€位置對應(yīng)的數(shù)字，有時候在某幾個特殊的環(huán)節(jié)我們可能會發(fā)現(xiàn)對序列的最少步驟操作可能不是唯一的，這其實是可以接受的，只要兩個不唯一的系列操作最終都是同樣的步驟數(shù)量就行，很明顯的例子就是上圖中大矩陣的藍色路徑，最后的部分，其實先刪除在取代，或者先取代再刪除其實都能達到同樣的效果，雖然序列操作最后有點不同，但都還是得花 5 步來完成序列的修改。

2. 動態(tài)規(guī)劃代碼

根據(jù)上面所解釋的邏輯，我們需要在迭代過程中不斷尋找 2x2 區(qū)域里面的最少步長，如果發(fā)現(xiàn)當(dāng)前的轉(zhuǎn)換前轉(zhuǎn)換后字符是一樣的，那就直接拷貝 2x2 區(qū)域里面的左上角值接著迭代：

def?Levenshtein(astr,?bstr):
????#?讀取輸入的兩個字符串并創(chuàng)建一個空的矩陣
????rows,?cols?=?(len(bstr)?+?1,?len(astr)?+?1)
????matrix?=?np.array([None]?*?(rows*cols)).reshape(rows,?cols)
????
????#?加上第一行與第一列的初始值
????matrix[0,?:]?=?np.arange(cols)
????matrix[:,?0]?=?np.arange(rows)
????
????#?對序列中的所有字符逐一遍歷
????for?i,?b?in?enumerate(bstr):
????????for?j,?a?in?enumerate(astr):
????????????#?如果兩個字符是一樣的,?...
????????????if?a?==?b:
????????????????#?...?直接拷貝左上角的值
????????????????matrix[i+1,?j+1]?=?matrix[i,?j]
????????????else:
????????????????#?否則就找?2x2?矩陣中最少步長
????????????????matrix[i+1,?j+1]?=?np.min([matrix[i+0,?j+1],
???????????????????????????????????????????matrix[i+0,?j+0],
???????????????????????????????????????????matrix[i+1,?j+0]])?+?1
????return?matrix


print(Levenshtein('ATTGTCT',?'GTAGCTT'))

輸出:
[[0,?1,?2,?3,?4,?5,?6,?7],
?[1,?1,?2,?3,?3,?4,?5,?6],
?[2,?2,?1,?2,?3,?3,?4,?5],
?[3,?2,?2,?2,?3,?4,?4,?5],
?[4,?3,?3,?3,?2,?3,?4,?5],
?[5,?4,?4,?4,?3,?3,?3,?4],
?[6,?5,?4,?4,?4,?3,?4,?3],
?[7,?6,?5,?4,?5,?4,?4,?4]]

其他類似的序列替換問題同樣可以通過這個簡單的動態(tài)規(guī)劃算法來計算出兩個序列之間的 編輯距離，例如我們可以把本來用來裝 DNA 信息的顏色框拿去裝一個單詞，這么一來我們就能分析兩句話之間有幾個單詞不一樣。總的來說，萊文斯坦算法是一個非常輕巧但同時在自然語言處理領(lǐng)域非常好用的一個算法！值得一試！

推薦閱讀

CondInst：性能和速度均超越Mask RCNN的實例分割模型

PyTorch 源碼解讀之 torch.autograd

centerX: 用新的視角的方式打開CenterNet

mmdetection最小復(fù)刻版(十一)：概率Anchor分配機制PAA深入分析

MMDetection新版本V2.7發(fā)布，支持DETR，還有YOLOV4在路上!

CNN：我不是你想的那樣

TF Object Detection 終于支持TF2了！

無需tricks，知識蒸餾提升ResNet50在ImageNet上準(zhǔn)確度至80%+

不妨試試MoCo，來替換ImageNet上pretrain模型！

重磅！一文深入深度學(xué)習(xí)模型壓縮和加速

從源碼學(xué)習(xí)Transformer！