DNA双序列全局比对Needleman-Wunsch算法教学设计

张小丹 肖春杨 卫泽刚 杨严硕 刘飞

摘要：文章针对大多数学生学习DNA双序列全局比对算法所面临的困惑，采用循序渐进，由浅入深，理论难点分析，Python代码实现，代码解析，拓展练习的思路对双序列全局比对Needleman-Wunsch算法过程进行详细讲解。首先对比对与回溯过程进行详细建模，将比对过程分解为得分矩阵构建和路径回溯两部分，然后提供了相应的Python代码实现及运行实例，进一步理解比对建模过程，并详细分析了代码难点。最后通过两个提升练习题，加强训练学生解决相关新问题的能力。

关键词：序列对比；全局比对；Needleman-Wunsch；教学设计

中图分类号：TP301      文献标识码：A

文章编号：1009-3044（2024）08-0178-03

开放科学（资源服务）标识码（OSID）

0 引言

序列比对是为了确定两个或多个序列之间的相似性或同源性，将不同序列按照一定规律进行的排列[1-3]。Needleman-Wunsch（NW） 和Smith-Waterman（SW） 是双序列全局比对和局部比对的经典算法[4-5]，均采用动态规划算法策略进行打分，然后通过路径回溯找到两条序列之间的最优比对结果。

在生物信息学相关课程教学过程中，序列比对是重点，同时也是难点。难点主要体现在以下两个方面：1） 序列比对过程中采用的动态规划算法是一种抽象计算思维建模方法，将待求解问题分解成若干子问题，通过对子问题的迭代求解，获得待求解问题的最优解。需要学生有较强的算法理解能力[6]。2） 序列比对需要用程序代码实现出来并成功运行，尽管大多数学生在大一阶段学习了C/C++程序设计基础相关课程，但都是一些基本语法知识，并没有针对实际问题进行求解，缺少一定的算法编程能力。因而在课堂教学中难免会出现学生缺乏积极性、不能跟随老师推导进行思考的问题，这时就需要老师根据序列比对算法的具体特点，将知识讲解得通俗易懂，启发学生们的主动学习兴趣[7]。

对于实现序列比对的编程语言来说，本文结合目前流行的Python语言进行比对过程讲解与代码实现。相对于其他编程语言，Python语法更加简单，简化了很多语法操作，可读性强，便于理解与实现，可使学生更加专注于如何求解比对问题，避免花费很多精力学习编程语法及规则[8]。

1 NW算法比对过程

对于给定的两条序列q和t，其详细过程如下：首先确定好插入（insertion） 、删除（deletion） 、匹配（match） 、替换（substitution或mismatch） 的惩罚得分，然后创建[len（q）+1]×[len（t）+1]大小的打分（得分）矩阵score，其中len（q）和len（t）分别表示序列q和t的长度，然后根据插入、删除、匹配、替换情况，通过选取最优得分计算每个矩阵元素的得分，同时记录得分路径。打分矩阵每个元素的得分计算公式为：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列相应元素位置的得分值，r是插入与删除错误的惩罚得分（默认值为-2） ，又称空位（gap） 罚分，δ函数是判断匹配位置的得分函数，定义为：δ（x， y）=α（若x=y） ，δ（x， y）=β（若x≠y） ，即：如果两个碱基相同，即匹配，则得分为α（默认值2） ，否则为替换错误，得分为β（默认值-2） 。可以观察到，每次计算的得分值是从匹配、替换、插入、删除四种情况下选取的最大值，即每一步计算得分均为当前子序列（子问题）的最优比对结果，当所有最优比对结果计算完毕后，即可得到两条序列的完整比对结果（全局优化）。同时记录每个元素得分值的计算路径（回溯方向，即当前最大得分值的取值方向），最后，从得分矩阵右下角元素开始，通过路径回溯找到最终的序列比对结果，路径回溯的计算公式为：

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

其中qalign与talign分别表示序列q与t的比对结果（字符串），‘-表示空位，符号“+”表示字符串拼接，p（i， j）表示得分矩阵中元素score（i， j）的回溯方向。式（1） 表示，从矩阵中最右下角得分开始，通过路径回溯得到两条序列的详细比对结果，若得分路径来自斜方向（用0表示），则根据序列q和t，按照从后往前的方向（逆序方向）各取出一个碱基，加入qalign与talign，若得分路径来自上方，提取q序列相应位置的碱基加入qalign，将空位‘-加入talign；若路径来自左边，提取q序列相应位置的碱基加入qalign，将空位‘-加入talign，直到回溯到score矩阵的第二行第二列，即可得到序列q與t的最终完整比对结果qalign与talign。

2 NW算法实例演示

接下来通过一个简单例子加深学生对双序列全局对比算法的理解。给定两条序列q=TCCA，t=TCGCA，q长度为4，t长度为5。令匹配、空位、替换相应的惩罚得分为2、-2、-2。如图1所示，首先构建大小为5×6的矩阵，并根据下面公式初始化第一行和第一列：score（1， j） = -gap×j，score（i， 1） = -gap×i，其中j =1，2，…，len（t）+1，i =1，2，…，len（s）+1。然后从第二行第二列开始，根据上述得分计算公式计算得分，如图1（b） 所示，在计算score（2， 2）得分时，最终从三个方向得分中选取最大值作为该单元的分值，如图1（b） 所示，该单元最大值为斜方向的得分：2，即匹配得分，则score（2， 2）=2，同时记录得分路径p（2， 2）=1（1表示对角线方向）。以此类推，最终完成整个矩阵的得分计算，如图1（c） 所示。最后从矩阵中最右下角开始，即p（5， 6）开始，根据公式（3） ，通过路径回溯得到两条序列的相似比对结果：qalign=TC-CA，talign=TCGCA，如图2所示。

3 NW算法Python代码实现

通过上述对双序列全局比对算法的过程分析和实例详解，可将NW双序列全局比对的过程分为两步：（1） 创建打分矩阵和路径矩阵并进行打分计算；（2） 通过回溯输出最终的比对结果。该问题可以用Python编程实现，限于篇幅，代码下载地址为：https：//pan.baidu.com/s/1UuKoZ0O6bulv8Qzemq-IVg？pwd=3fh5，可在Windows或Linux系统下直接运行。

4 代码难点分析

上述代码以函数定义的形式对DNA双序列全局比对算法进行编程实现，直观体现了全局对比算法的过程。其中有以下几处代码是学生难以理解，需要进行详细解释。

1） 创建矩阵。通过第61行代码加载numpy包（简写为np） ，可以快速创建不同大小的矩阵（第5-6行），以上函数通过numpy的zeros函数实现矩阵的创建及初始化（矩阵元素值全部初始化为0） 。

2） 函数定义。Python代码定义函数非常简单，只需通过关键字def即可定义不同的函数及参数传递，函数返回值通过return关键字即可实现。在以上程序代码中，定义了两个函数：打分函数cal_score（第2行）和回溯函数backtrack（第31行）。其中打分函数cal_score需要两个参数：即待比对的两条序列，并返回2个值：打分矩阵scores和路径矩阵path。回溯函数backtrack需要一个参数：打分路径path，并返回3个值：两条序列的比对结果aligned_seq1、aligned_seq1和匹配结果middle。

3） 如何用程序实现每个得分的计算，并选取最大的值。在计算每一个单元的得分值时，需要计算三个方向的得分值（代码第20-28行），然后再选取最大。在计算匹配或替换得分时，首先需要判断对应位置的碱基是否相等（第41行代码），如果相等，需要加上匹配得分，否则加上替换得分。

4） 回溯矩阵如何记录路径。在创建打分矩阵的时候同时创建路径矩阵（第6行），然后根据最大值得分方向记录得分路径（代码第28行）。其中0表示最大值源自左侧方向，1表示最大值源自对角线方向，2表示最大值源自上方方向。

5） Python中字符串操作。由于比对的是DNA序列，在Python直接用字符串表示即可，既方便又简单明了。Python字符串用两个单引号或者双引号表示，然后直接赋值给相应变量即可，如第62-63行代码。在路径回溯过程中，需要根据公式（3） 提取碱基并进行合并，在Python中可直接用“+”符号对字符串进行操作，表示两个字符串的拼接，如上述第39-40行代码所示。

5 加强提升练习

在讲解完NW算法过程后，通过增加练习题让学生课后练习，巩固所学的知识。如果根据新的两条序列可以计算出比对结果，说明学生理解掌握了NW算法。在课后练习阶段，可以通过以下两个问题进行提升：

1） 根据序列比对结果计算两条序列间的相似性，定义为：匹配的碱基数占比对长度的百分比。

2） 相对于NW算法，SW算法是局部序列比对。如何得到局部比对结果，可以在公式（1） 的基础上改动为：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方法一样，然后让学生理解并用Python编程实现。

6 总结

本文首先对Needleman-Wunsch算法的比对与回溯过程进行详细建模，将序列比对过程分解为得分矩阵构建和路径回溯两部分，然后提供了相应的Python编程代码实现及运行实例，进一步理解比对建模过程，并针对代码难点进行详细分析。最后通过两个提升练习题，加强训练学生解决相关新问题的能力。这将为以能力目标导向的生物信息相关课程教学改革提供新思路，为培养学生未来职业发展所需问题解决能力提供教学设计案例。

参考文献：

[1] 卫泽刚.三代测序序列比对与微生物操作单元划分算法[D].西北工业大学，2019.

[2] 罗静初.双序列比对基础和应用实例[J].生物信息学，2023，21（1）：1-19.

[3] 卫泽刚，侯一凡，张小丹，等.微生物操作分類单元划分算法研究[J].宝鸡文理学院学报（自然科学版），2022，42（1）：80-88.

[4] 张小丹，李喆，卫泽刚，等.基于k-mer词频向量的九种DNA序列相似性计算方法比较分析[J].科学技术创新，2023（21）：106-111.

[5] 甘秋云.基于动态规划的Needleman-Wunsch双序列比对算法的分析与研究[J].计算机工程与科学，2021，43（2）：340-346.

[6] 张立臣，王小明.面向问题解决能力培养的算法课程教学设计——以动态规划算法教学为例[J].高教学刊，2021（11）：105-109.

[7] 张小丹，范兴国，卫泽刚，等.多线程排序算法教学设计——以插入排序为例[J].科技视界，2022（25）：121-123.

[8] 卫泽刚，张小丹，赵军娣，等.基于Matlab的选择排序算法教学设计[J].科技视界，2022（22）：86-88.

【通联编辑：王 力】