張小丹 肖春楊 衛(wèi)澤剛 楊嚴(yán)碩 劉飛

摘要：文章針對(duì)大多數(shù)學(xué)生學(xué)習(xí)DNA雙序列全局比對(duì)算法所面臨的困惑，采用循序漸進(jìn)，由淺入深，理論難點(diǎn)分析，Python代碼實(shí)現(xiàn)，代碼解析，拓展練習(xí)的思路對(duì)雙序列全局比對(duì)Needleman-Wunsch算法過(guò)程進(jìn)行詳細(xì)講解。首先對(duì)比對(duì)與回溯過(guò)程進(jìn)行詳細(xì)建模，將比對(duì)過(guò)程分解為得分矩陣構(gòu)建和路徑回溯兩部分，然后提供了相應(yīng)的Python代碼實(shí)現(xiàn)及運(yùn)行實(shí)例，進(jìn)一步理解比對(duì)建模過(guò)程，并詳細(xì)分析了代碼難點(diǎn)。最后通過(guò)兩個(gè)提升練習(xí)題，加強(qiáng)訓(xùn)練學(xué)生解決相關(guān)新問(wèn)題的能力。

關(guān)鍵詞：序列對(duì)比；全局比對(duì)；Needleman-Wunsch；教學(xué)設(shè)計(jì)

中圖分類號(hào)：TP301? ? ? 文獻(xiàn)標(biāo)識(shí)碼：A

文章編號(hào)：1009-3044（2024）08-0178-03

開(kāi)放科學(xué)（資源服務(wù)）標(biāo)識(shí)碼（OSID）

0 引言

序列比對(duì)是為了確定兩個(gè)或多個(gè)序列之間的相似性或同源性，將不同序列按照一定規(guī)律進(jìn)行的排列[1-3]。Needleman-Wunsch（NW） 和Smith-Waterman（SW） 是雙序列全局比對(duì)和局部比對(duì)的經(jīng)典算法[4-5]，均采用動(dòng)態(tài)規(guī)劃算法策略進(jìn)行打分，然后通過(guò)路徑回溯找到兩條序列之間的最優(yōu)比對(duì)結(jié)果。

在生物信息學(xué)相關(guān)課程教學(xué)過(guò)程中，序列比對(duì)是重點(diǎn)，同時(shí)也是難點(diǎn)。難點(diǎn)主要體現(xiàn)在以下兩個(gè)方面：1） 序列比對(duì)過(guò)程中采用的動(dòng)態(tài)規(guī)劃算法是一種抽象計(jì)算思維建模方法，將待求解問(wèn)題分解成若干子問(wèn)題，通過(guò)對(duì)子問(wèn)題的迭代求解，獲得待求解問(wèn)題的最優(yōu)解。需要學(xué)生有較強(qiáng)的算法理解能力[6]。2） 序列比對(duì)需要用程序代碼實(shí)現(xiàn)出來(lái)并成功運(yùn)行，盡管大多數(shù)學(xué)生在大一階段學(xué)習(xí)了C/C++程序設(shè)計(jì)基礎(chǔ)相關(guān)課程，但都是一些基本語(yǔ)法知識(shí)，并沒(méi)有針對(duì)實(shí)際問(wèn)題進(jìn)行求解，缺少一定的算法編程能力。因而在課堂教學(xué)中難免會(huì)出現(xiàn)學(xué)生缺乏積極性、不能跟隨老師推導(dǎo)進(jìn)行思考的問(wèn)題，這時(shí)就需要老師根據(jù)序列比對(duì)算法的具體特點(diǎn)，將知識(shí)講解得通俗易懂，啟發(fā)學(xué)生們的主動(dòng)學(xué)習(xí)興趣[7]。

對(duì)于實(shí)現(xiàn)序列比對(duì)的編程語(yǔ)言來(lái)說(shuō)，本文結(jié)合目前流行的Python語(yǔ)言進(jìn)行比對(duì)過(guò)程講解與代碼實(shí)現(xiàn)。相對(duì)于其他編程語(yǔ)言，Python語(yǔ)法更加簡(jiǎn)單，簡(jiǎn)化了很多語(yǔ)法操作，可讀性強(qiáng)，便于理解與實(shí)現(xiàn)，可使學(xué)生更加專注于如何求解比對(duì)問(wèn)題，避免花費(fèi)很多精力學(xué)習(xí)編程語(yǔ)法及規(guī)則[8]。

1 NW算法比對(duì)過(guò)程

對(duì)于給定的兩條序列q和t，其詳細(xì)過(guò)程如下：首先確定好插入（insertion） 、刪除（deletion） 、匹配（match） 、替換（substitution或mismatch） 的懲罰得分，然后創(chuàng)建[len（q）+1]×[len（t）+1]大小的打分（得分）矩陣score，其中l(wèi)en（q）和len（t）分別表示序列q和t的長(zhǎng)度，然后根據(jù)插入、刪除、匹配、替換情況，通過(guò)選取最優(yōu)得分計(jì)算每個(gè)矩陣元素的得分，同時(shí)記錄得分路徑。打分矩陣每個(gè)元素的得分計(jì)算公式為：score（i， j）=max[score（i-1， j-1）+δ（qi-1， tj-1）， score（i， j-1）+r， score（i-1， j）+r]，其中score（i， j）是得分矩陣score中第i行，第j列相應(yīng)元素位置的得分值，r是插入與刪除錯(cuò)誤的懲罰得分（默認(rèn)值為-2） ，又稱空位（gap） 罰分，δ函數(shù)是判斷匹配位置的得分函數(shù)，定義為：δ（x， y）=α（若x=y） ，δ（x， y）=β（若x≠y） ，即：如果兩個(gè)堿基相同，即匹配，則得分為α（默認(rèn)值2） ，否則為替換錯(cuò)誤，得分為β（默認(rèn)值-2） 。可以觀察到，每次計(jì)算的得分值是從匹配、替換、插入、刪除四種情況下選取的最大值，即每一步計(jì)算得分均為當(dāng)前子序列（子問(wèn)題）的最優(yōu)比對(duì)結(jié)果，當(dāng)所有最優(yōu)比對(duì)結(jié)果計(jì)算完畢后，即可得到兩條序列的完整比對(duì)結(jié)果（全局優(yōu)化）。同時(shí)記錄每個(gè)元素得分值的計(jì)算路徑（回溯方向，即當(dāng)前最大得分值的取值方向），最后，從得分矩陣右下角元素開(kāi)始，通過(guò)路徑回溯找到最終的序列比對(duì)結(jié)果，路徑回溯的計(jì)算公式為：

[qalign='-'，if pi，j==0qi，otherwise+ qaligntalign='-'，if pi，j==2ti，otherwise + talign]? ? （1）

其中qalign與talign分別表示序列q與t的比對(duì)結(jié)果（字符串），‘-表示空位，符號(hào)“+”表示字符串拼接，p（i， j）表示得分矩陣中元素score（i， j）的回溯方向。式（1） 表示，從矩陣中最右下角得分開(kāi)始，通過(guò)路徑回溯得到兩條序列的詳細(xì)比對(duì)結(jié)果，若得分路徑來(lái)自斜方向（用0表示），則根據(jù)序列q和t，按照從后往前的方向（逆序方向）各取出一個(gè)堿基，加入qalign與talign，若得分路徑來(lái)自上方，提取q序列相應(yīng)位置的堿基加入qalign，將空位‘-加入talign；若路徑來(lái)自左邊，提取q序列相應(yīng)位置的堿基加入qalign，將空位‘-加入talign，直到回溯到score矩陣的第二行第二列，即可得到序列q與t的最終完整比對(duì)結(jié)果qalign與talign。

2 NW算法實(shí)例演示

接下來(lái)通過(guò)一個(gè)簡(jiǎn)單例子加深學(xué)生對(duì)雙序列全局對(duì)比算法的理解。給定兩條序列q=TCCA，t=TCGCA，q長(zhǎng)度為4，t長(zhǎng)度為5。令匹配、空位、替換相應(yīng)的懲罰得分為2、-2、-2。如圖1所示，首先構(gòu)建大小為5×6的矩陣，并根據(jù)下面公式初始化第一行和第一列：score（1， j） = -gap×j，score（i， 1） = -gap×i，其中j =1，2，…，len（t）+1，i =1，2，…，len（s）+1。然后從第二行第二列開(kāi)始，根據(jù)上述得分計(jì)算公式計(jì)算得分，如圖1（b） 所示，在計(jì)算score（2， 2）得分時(shí)，最終從三個(gè)方向得分中選取最大值作為該單元的分值，如圖1（b） 所示，該單元最大值為斜方向的得分：2，即匹配得分，則score（2， 2）=2，同時(shí)記錄得分路徑p（2， 2）=1（1表示對(duì)角線方向）。以此類推，最終完成整個(gè)矩陣的得分計(jì)算，如圖1（c） 所示。最后從矩陣中最右下角開(kāi)始，即p（5， 6）開(kāi)始，根據(jù)公式（3） ，通過(guò)路徑回溯得到兩條序列的相似比對(duì)結(jié)果：qalign=TC-CA，talign=TCGCA，如圖2所示。

3 NW算法Python代碼實(shí)現(xiàn)

通過(guò)上述對(duì)雙序列全局比對(duì)算法的過(guò)程分析和實(shí)例詳解，可將NW雙序列全局比對(duì)的過(guò)程分為兩步：（1） 創(chuàng)建打分矩陣和路徑矩陣并進(jìn)行打分計(jì)算；（2） 通過(guò)回溯輸出最終的比對(duì)結(jié)果。該問(wèn)題可以用Python編程實(shí)現(xiàn)，限于篇幅，代碼下載地址為：https：//pan.baidu.com/s/1UuKoZ0O6bulv8Qzemq-IVg？pwd=3fh5，可在Windows或Linux系統(tǒng)下直接運(yùn)行。

4 代碼難點(diǎn)分析

上述代碼以函數(shù)定義的形式對(duì)DNA雙序列全局比對(duì)算法進(jìn)行編程實(shí)現(xiàn)，直觀體現(xiàn)了全局對(duì)比算法的過(guò)程。其中有以下幾處代碼是學(xué)生難以理解，需要進(jìn)行詳細(xì)解釋。

1） 創(chuàng)建矩陣。通過(guò)第61行代碼加載numpy包（簡(jiǎn)寫為np） ，可以快速創(chuàng)建不同大小的矩陣（第5-6行），以上函數(shù)通過(guò)numpy的zeros函數(shù)實(shí)現(xiàn)矩陣的創(chuàng)建及初始化（矩陣元素值全部初始化為0） 。

2） 函數(shù)定義。Python代碼定義函數(shù)非常簡(jiǎn)單，只需通過(guò)關(guān)鍵字def即可定義不同的函數(shù)及參數(shù)傳遞，函數(shù)返回值通過(guò)return關(guān)鍵字即可實(shí)現(xiàn)。在以上程序代碼中，定義了兩個(gè)函數(shù)：打分函數(shù)cal_score（第2行）和回溯函數(shù)backtrack（第31行）。其中打分函數(shù)cal_score需要兩個(gè)參數(shù)：即待比對(duì)的兩條序列，并返回2個(gè)值：打分矩陣scores和路徑矩陣path?；厮莺瘮?shù)backtrack需要一個(gè)參數(shù)：打分路徑path，并返回3個(gè)值：兩條序列的比對(duì)結(jié)果aligned_seq1、aligned_seq1和匹配結(jié)果middle。

3） 如何用程序?qū)崿F(xiàn)每個(gè)得分的計(jì)算，并選取最大的值。在計(jì)算每一個(gè)單元的得分值時(shí)，需要計(jì)算三個(gè)方向的得分值（代碼第20-28行），然后再選取最大。在計(jì)算匹配或替換得分時(shí)，首先需要判斷對(duì)應(yīng)位置的堿基是否相等（第41行代碼），如果相等，需要加上匹配得分，否則加上替換得分。

4） 回溯矩陣如何記錄路徑。在創(chuàng)建打分矩陣的時(shí)候同時(shí)創(chuàng)建路徑矩陣（第6行），然后根據(jù)最大值得分方向記錄得分路徑（代碼第28行）。其中0表示最大值源自左側(cè)方向，1表示最大值源自對(duì)角線方向，2表示最大值源自上方方向。

5） Python中字符串操作。由于比對(duì)的是DNA序列，在Python直接用字符串表示即可，既方便又簡(jiǎn)單明了。Python字符串用兩個(gè)單引號(hào)或者雙引號(hào)表示，然后直接賦值給相應(yīng)變量即可，如第62-63行代碼。在路徑回溯過(guò)程中，需要根據(jù)公式（3） 提取堿基并進(jìn)行合并，在Python中可直接用“+”符號(hào)對(duì)字符串進(jìn)行操作，表示兩個(gè)字符串的拼接，如上述第39-40行代碼所示。

5 加強(qiáng)提升練習(xí)

在講解完NW算法過(guò)程后，通過(guò)增加練習(xí)題讓學(xué)生課后練習(xí)，鞏固所學(xué)的知識(shí)。如果根據(jù)新的兩條序列可以計(jì)算出比對(duì)結(jié)果，說(shuō)明學(xué)生理解掌握了NW算法。在課后練習(xí)階段，可以通過(guò)以下兩個(gè)問(wèn)題進(jìn)行提升：

1） 根據(jù)序列比對(duì)結(jié)果計(jì)算兩條序列間的相似性，定義為：匹配的堿基數(shù)占比對(duì)長(zhǎng)度的百分比。

2） 相對(duì)于NW算法，SW算法是局部序列比對(duì)。如何得到局部比對(duì)結(jié)果，可以在公式（1） 的基礎(chǔ)上改動(dòng)為：score（i， j）=max[score（i-1， j-1）+δ（qi-1， tj-1）， score（i， j-1）+r， score（i-1， j）+r， 0]。與NW方法主要不同之處為：SW矩陣第一行與第一列全部初始化為0，回溯方法是從score的得分最大值開(kāi)始回溯到得分為0的位置，回溯方法與NW方法一樣，然后讓學(xué)生理解并用Python編程實(shí)現(xiàn)。

6 總結(jié)

本文首先對(duì)Needleman-Wunsch算法的比對(duì)與回溯過(guò)程進(jìn)行詳細(xì)建模，將序列比對(duì)過(guò)程分解為得分矩陣構(gòu)建和路徑回溯兩部分，然后提供了相應(yīng)的Python編程代碼實(shí)現(xiàn)及運(yùn)行實(shí)例，進(jìn)一步理解比對(duì)建模過(guò)程，并針對(duì)代碼難點(diǎn)進(jìn)行詳細(xì)分析。最后通過(guò)兩個(gè)提升練習(xí)題，加強(qiáng)訓(xùn)練學(xué)生解決相關(guān)新問(wèn)題的能力。這將為以能力目標(biāo)導(dǎo)向的生物信息相關(guān)課程教學(xué)改革提供新思路，為培養(yǎng)學(xué)生未來(lái)職業(yè)發(fā)展所需問(wèn)題解決能力提供教學(xué)設(shè)計(jì)案例。

參考文獻(xiàn)：

[1] 衛(wèi)澤剛.三代測(cè)序序列比對(duì)與微生物操作單元?jiǎng)澐炙惴╗D].西北工業(yè)大學(xué)，2019.

[2] 羅靜初.雙序列比對(duì)基礎(chǔ)和應(yīng)用實(shí)例[J].生物信息學(xué)，2023，21（1）：1-19.

[3] 衛(wèi)澤剛，侯一凡，張小丹，等.微生物操作分類單元?jiǎng)澐炙惴ㄑ芯縖J].寶雞文理學(xué)院學(xué)報(bào)（自然科學(xué)版），2022，42（1）：80-88.

[4] 張小丹，李喆，衛(wèi)澤剛，等.基于k-mer詞頻向量的九種DNA序列相似性計(jì)算方法比較分析[J].科學(xué)技術(shù)創(chuàng)新，2023（21）：106-111.

[5] 甘秋云.基于動(dòng)態(tài)規(guī)劃的Needleman-Wunsch雙序列比對(duì)算法的分析與研究[J].計(jì)算機(jī)工程與科學(xué)，2021，43（2）：340-346.

[6] 張立臣，王小明.面向問(wèn)題解決能力培養(yǎng)的算法課程教學(xué)設(shè)計(jì)——以動(dòng)態(tài)規(guī)劃算法教學(xué)為例[J].高教學(xué)刊，2021（11）：105-109.

[7] 張小丹，范興國(guó)，衛(wèi)澤剛，等.多線程排序算法教學(xué)設(shè)計(jì)——以插入排序?yàn)槔齕J].科技視界，2022（25）：121-123.

[8] 衛(wèi)澤剛，張小丹，趙軍娣，等.基于Matlab的選擇排序算法教學(xué)設(shè)計(jì)[J].科技視界，2022（22）：86-88.

【通聯(lián)編輯：王 力】