驚了，驚了! 原來這個算法還能用在核酸檢測
核酸檢測經(jīng)此一疫已經(jīng)成為了人盡皆知的檢測手段，但是你可曾想過經(jīng)過儀器檢測出來的 DNA序列是如何與正常情況下的 DNA 序列做對比的呢？面對整村的核酸檢測結果，計算工作肯定需要電腦來完成，如果我們想知道一個受試者的 DNA 序列與正常情況下的 DNA 序列有至少幾個個元素差異時，我們就需要一個算法來解決，它叫 萊文斯坦 (levenshtein) 算法，而這個算法背后有一個主要的技術：動態(tài)規(guī)劃！
import?numpy?as?np
本文將只通過 numpy 從頭構建萊文斯坦算法，如果還沒安裝這個庫的讀者們可以通過下面指令在終端安裝：
pip?install?numpy
編輯距離
DNA 序列主要可以用四種英文字母表示，ATCG，至于什么樣的檢測結果要對應到什么樣的英文字母不在這篇文章的討論范圍，我們只討論如何處理這些序列數(shù)據(jù)。假如檢測站的樣本里我們的到了兩組 DNA 序列：
判斷這兩組序列與目標序列之間的差異前，我們需要先定義差異的含義，也就是編輯距離，面對一個序列，我們需要經(jīng)過幾步操作才能夠使其和目標序列變得一模一樣，如果是三步，那么編輯距離就是3。但是廣義地說，修改一個序列使其和另一個序列一樣的方法有無限多種，因此我們需要給另一個約束條件，就是如何用最少的步驟來完成修改的操作，而具體操作動作只能包含下面三種方式：
插入 刪除 替換
為了有效解決這個問題，動態(tài)規(guī)劃就是本篇文章所介紹算法的核心邏輯！
動態(tài)規(guī)劃
動態(tài)規(guī)劃的核心概念就是把一個大的問題拆分成很多個小的問題分別擊破，使得整體問題難度降低，進而更容易找出問題的最優(yōu)答案。面對 DNA 序列問題，大問題就是一串完整的序列到底有哪些不同，而小問題則是對應到每個位置上的字符，到底是應該采取什么樣的動作，先給出計算兩個不同序列編輯距離的計算過程。
為了把 "ATTGTCT" 改成 "GTAGCTT"，根據(jù)矩陣里最右下角的計算結果，最少的操作步驟就是 4 步，至于這四個步驟具體是 插入，刪除，替換 怎么樣的排列組合，那得依據(jù)計算矩陣右下角數(shù)字的過程中小編所標注的框的顏色和對應到一個 2x2 的行為矩陣來判定：
如果我們把行為矩陣里的右下角位置對應到大矩陣里每一個紅色框，紅框的數(shù)字全是從大矩陣中橙色框 + 1 而得到的，而橙色框?qū)诵袨榫仃噭t是 替換 的操作！通過同樣的規(guī)律，我們也可以從藍色的路徑了解到從 "ATTGTCT" 變到 "GTA" 最少需要 5 步，從 "ATT" 變到 "GTAGCTT" 則需要 4 步，至于具體填上數(shù)字的細節(jié)，小編馬上接著說！
1. 原理細節(jié)
在最一開始什么數(shù)字都不存在在矩陣中時，只有 變化前 的序列被放在最上面，和 變化后 的序列被放在最左邊，而序列的最前面需要填上一個代表本來什么都沒有的站位符，我們可以理解為一個字符串在最易開始是什么元素都沒有的，因此不論是變化前后兩個序列最前面都得加上空的占位符。接著，我們就能開始逐一解大問題中的每一個小問題，從空的占位符變成任意長度的序列所需要的步數(shù)肯定等于序列的長度，因此第一個行列就可以理所當然的填上 0 1 2 ... 7。
剩下來還沒填上數(shù)字的格子，我們就需要以 2x2 所排部的四個格子為單位去計算最小的步驟，而這四個格子也就剛好對應了行為矩陣的位置，為了確定小問題的最優(yōu)操作，我們需要知道分別操作了 插入，刪除，替換 之后，最小的成本是多少。
1-1. 不同元素 - 綠色方形塊
如果對應位置的字母需要被改變，例如圖中 A 要換成 G，至少 1 步 (替換) 需要操作，但這是因為一開始問題一目了然，能夠通過直覺看出問題的解，實際上在動態(tài)規(guī)劃的思維中，我們需要發(fā)現(xiàn)一個能夠一直被重復使用的規(guī)律，因此 A 換成 G 是 1 步的答案背后應該要根據(jù) 2x2 范圍里的數(shù)字去比較：
0 + 1 表示替換 A 成 G 所需要的步數(shù) 1 + 1 表示插入 A 成 G 所需要的步數(shù) 1 + 1 表示刪除 A 成 G 所需要的步數(shù)
而最后根據(jù)這三個操作的結果來看，證明了最少步數(shù)是替換操作，因此這步將遵循此操作，并把 0 + 1 的結果寫進 2x2 矩陣范圍的右下角。
1-2. 相同元素 - 黃色方形塊
相反地，如果我們遇到要把 G 換成 G 的情況，就直接把 2x2 矩陣范圍的左上角數(shù)字抄到右下角，因為同樣的字符是不需要做任何變動的，也就沒有任何一個操作合適，所以最少步數(shù)的結論就還是 3 步。了解了處理 不同元素 與 相同元素 的計算規(guī)則之后，就可以逐一填上矩陣里每一個位置對應的數(shù)字，有時候在某幾個特殊的環(huán)節(jié)我們可能會發(fā)現(xiàn)對序列的最少步驟操作可能不是唯一的，這其實是可以接受的，只要兩個不唯一的系列操作最終都是同樣的步驟數(shù)量就行，很明顯的例子就是上圖中大矩陣的藍色路徑，最后的部分，其實先刪除在取代，或者先取代再刪除其實都能達到同樣的效果，雖然序列操作最后有點不同，但都還是得花 5 步來完成序列的修改。
2. 動態(tài)規(guī)劃代碼
根據(jù)上面所解釋的邏輯，我們需要在迭代過程中不斷尋找 2x2 區(qū)域里面的最少步長，如果發(fā)現(xiàn)當前的轉(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 信息的顏色框拿去裝一個單詞，這么一來我們就能分析兩句話之間有幾個單詞不一樣。總的來說，萊文斯坦算法是一個非常輕巧但同時在自然語言處理領域非常好用的一個算法！值得一試！
推薦閱讀
CondInst：性能和速度均超越Mask RCNN的實例分割模型
mmdetection最小復刻版(十一)：概率Anchor分配機制PAA深入分析
MMDetection新版本V2.7發(fā)布，支持DETR，還有YOLOV4在路上!
無需tricks，知識蒸餾提升ResNet50在ImageNet上準確度至80%+
不妨試試MoCo，來替換ImageNet上pretrain模型！
mmdetection最小復刻版(七)：anchor-base和anchor-free差異分析
mmdetection最小復刻版(四)：獨家yolo轉(zhuǎn)化內(nèi)幕
機器學習算法工程師
? ??? ? ? ? ? ? ? ? ? ? ? ??????? ??一個用心的公眾號
?
