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

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

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

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

文章編號：1009-3044（2024）08-0178-03

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

0 引言

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

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

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

1 NW算法比對過程

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

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

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

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

接下來通過一個簡單例子加深學(xué)生對雙序列全局對比算法的理解。給定兩條序列q=TCCA，t=TCGCA，q長度為4，t長度為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。然后從第二行第二列開始，根據(jù)上述得分計算公式計算得分，如圖1（b） 所示，在計算score（2， 2）得分時，最終從三個方向得分中選取最大值作為該單元的分值，如圖1（b） 所示，該單元最大值為斜方向的得分：2，即匹配得分，則score（2， 2）=2，同時記錄得分路徑p（2， 2）=1（1表示對角線方向）。以此類推，最終完成整個矩陣的得分計算，如圖1（c） 所示。最后從矩陣中最右下角開始，即p（5， 6）開始，根據(jù)公式（3） ，通過路徑回溯得到兩條序列的相似比對結(jié)果：qalign=TC-CA，talign=TCGCA，如圖2所示。

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

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

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

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

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

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

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

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

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

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

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

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

2） 相對于NW算法，SW算法是局部序列比對。如何得到局部比對結(jié)果，可以在公式（1） 的基礎(chǔ)上改動為：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的得分最大值開始回溯到得分為0的位置，回溯方法與NW方法一樣，然后讓學(xué)生理解并用Python編程實(shí)現(xiàn)。

6 總結(jié)

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

參考文獻(xiàn)：

[1] 衛(wèi)澤剛.三代測序序列比對與微生物操作單元劃分算法[D].西北工業(yè)大學(xué)，2019.

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

[3] 衛(wèi)澤剛，侯一凡，張小丹，等.微生物操作分類單元劃分算法研究[J].寶雞文理學(xué)院學(xué)報（自然科學(xué)版），2022，42（1）：80-88.

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

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

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

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

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

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