用Python实现Needleman-Wunsch和Smith-Waterman算法

Python实现

1 简单介绍

Needleman-Wunsch和Smith-Waterman算法都是Pairwise alignment算法，用于双序列比对。

• Needleman-Wunsch是全局序列比对算法，最终会得到在全局上2条序列最佳的匹配结果
• Smith-Waterman是局部序列比对算法，最终得到的是2条序列在局部的最佳匹配片段，即挑选出得分最高的比对片段，而不一定是两个序列全部比对。

:::success 为什么要用局部序列比对算法呢？
功能相关的蛋白之间虽然整体的序列相差甚远，但是常常具有相同的功能域。对于这种情况，仅靠全局比对算法是不行的；另一方面，70年代内含子的发现，使得在进行核酸序列比对时，必须要能正确处理内含子导致的大片段差异。局部序列比对算法因此诞生。 :::

2 算法实现介绍

两个算法都是用了动态规划算法，什么叫动态规划呢。先思考下，若我们想要得到两条序列间最佳比对结果的话，你会怎么解决，穷举法？但是此方法需要穷举的次数随着序列长度的增加而不断增大，而且要考虑两条序列在哪里插入空位，比对结果的长度范围为n(两条序列完全匹配上)-2n(两条序列均对到空位)。

比对序列的长度越长，需要穷举的次数越多，这对我们来说将会耗时很多，这是不太现实的。所以我们需要一种新的方法进行比对——动态规划。动态规划算法是一种求取最优解的方法，该方法通过将问题拆解为若干子问题，求取子问题的最优解，所有子问题的最优解合并即为全局最优解。

下面具体介绍下两种算法的细节。

这两个算法都是先建立一个打分矩阵，考虑每个cell从各个方向的来源得分，记录最大得分的来源，最后按照规则回溯，从而得到比对结果。

2.1 来源方向

来源方向有三个：横着来的horizontal，竖着来的vertical，对角线来的diagonal。

其中horizontal，vertical对应gap，而diagonal对应mismatch和match（diagonal来源的比对可以是确实比对上，也可以是强行比对）。

2.2 打分矩阵

对于全局比对的Needleman-Wunsch算法而言：就考虑来自三个方向的分数情况，记录最大得分的来源

而对于局部比对的Smith-Waterman算法而言：还需要考虑得分>=0，如果得分少于0，则不记录方向。

:::success Note

• 对于diagonal方向的得分，有两种方式
  • 一种就是看两个碱基match不match，match就+x分，mismatch就-y分
  • 一种是替代矩阵M，M[A][B]为得分，经典的替代矩阵PAM250和BLOSUM62，其中PAM250是基于高同源序列构建，BLOSUM62基于远程同源序列构建，因此寻找远程同源序列一般用BLOSUM62。
• 对于vertical和horizonal方向的罚分，也有流行的两种方式
  • 一种是线性罚分，出现一个gap，我就扣k分
    • 举例，设置gap_open=-2。连续出现三个gap的罚分就是
  • 另一种叫affine gap penalty，翻译叫仿射罚分，即在打开第一个gap的时候引入gap open罚分，而在该gap的基础上进行延续则采用gap extension罚分。之所以这样做，相比从生物学原理来说，应该更倾向于仅出现一段长gap而不是多个短gap。
    • 举例，设置gap_open=-2，gap_extension=-1。连续出现三个gap的罚分就是%5Ctimes-1%3D-4#card=math&code=-2%2B%283-1%29%5Ctimes-1%3D-4&id=ZUKux)

:::

2.3 回溯

Needleman-Wunsch就是头对头、尾对尾的算法，所以会从最右下角进行回溯

Smith-Waterman只考虑找出最大得分的比对结果，所以会找出最高分的位置，开始回溯，直到cell的分数为0结束

:::success Note：
两个算法回溯起始的cell的score值就是最终序列的比对得分。 :::

Question： 回溯时，横着和竖着的方向代表存在gap，那怎么知道gap是在哪条序列上呢？ Answer：看方向，垂直箭头的那个序列，就在对应碱基前加一个gap。比如横着的箭头，说明左侧的seq存在gap，竖着的箭头，则是顶部的seq存在gap。

Question： 如果考虑affine gap penalty，正好一个块的上一个块就有两个不同的来源，那这时候，如果这个块是算它连续插入/删除还是不连续插入/删除的分数呢 ？ Answer：根据下面张图应该是2-1

总结