biopython 多序列比对

本文翻译自 Understanding Partial Order Alignment for Multiple Sequence Alignment，原文链接在

http://simpsonlab.github.io/2015/05/01/understanding-poa/

Jared 开发的 Nanopolish 工具使用 poaV2 工具来对测序序列进行错误修正，poaV2 则使用了偏序比对的算法，其算法描述详见 Lee，Grasso 和 Sharlow 的文章^1,2,3，该方法和 PacBio 的 PBDagCon 使用的方法类似。

基础知识

首篇介绍 POA 的文章指出，序列比对算法有时会产生一些没有生物学意义的结果，这些结果对双序列比对（pairwise alignment）以及高度保守序列之间的多序列比对（multiple alignment）是无害的，但是对更为一般性的多序列比对却会引起很多问题。比如说，考虑如下这些序列：

>seq1
CCGCTTTTCCGC
>seq2
CCGCAAAACCGC

上面的两条序列在选择最佳比对结果时有很多种选择，考虑如下情形，

CCGC----TTTTCGCG   CCGCTTTT----CCGC  CCGC-TT-TT--CGCG   CCGC-T-T-T-TCCGC
CCGCAAAA----CGCG   CCGC----AAAACCGC  CCGCA--A--AACCGC   CCGCA-A-A-A-CCGC

在不考虑生物学意义的情况下，中间的 8 个碱基之间可以任意匹配，那么对于每一条序列来说，就有 $C_{8}^{4}$ 种选择（译者注：将 gap - 看成 4 个球，放到 8 个不同的抽屉里，有 $C_{8}^{4}$ 种放法，gap 位置固定以后，T 或者 A 的位置也固定了），那么两条序列之间就共有 $C_{8}^{4} * C_{8}^{4}$

对于双序列比对，上面不同的比对方案可能相差不大，但对于多序列比对，不同的选择则对最后的结果有比较大的影响。那么，除了上面这种线性的表示方法，是否有更好方法来表达真实的情况？很自然的，我们想到可以用图来表示：

偏序比对图（partial order alignment graph）和字符串表示的比对不同，对于同一个碱基，会有确定的上下游关系。偏序比对图中的每个节点表示一个碱基，和字符串表示的比对类似，图中也会有隐含的上下游顺序关系，每个节点有零个或者多个上游节点和零个或多个下游节点。图中不允许有自连接和环，也就是说，偏序比对图是有向无环图（Directed Acyclic Graph，DAG）。

重复和重排序都可能有生物学意义，有多种比对会考虑这些情况，这样会将问题更一般化，并且有助于组装^4,5。对于以修正测序错误的 Naonpolish 来说，这样的一般化是非必须的。

Smith-Waterman

在考虑如何对偏序图进行比对之前，让我们先回忆一下我们如何对字符串序列进行比对。

对于 Needleman-Wunsh 算法和其变体，我们先假设有两个光标在碱基上，每条序列上有一个光标。让我们考虑光标在每个比对步骤中的位置，如何移动光标使得我们能够得到最佳比对结果？由于全局最优的比对路径必定来自局部最优的“移动”，这样我们就只需将问题转化为如何为光标选择以下三种移动方式：

• 两个光标同时向前移动，表明光标位置上的碱基互相匹配（match）
• 光标 1 向前移动，光标 2 保持不动，将序列 1 中的光标所在位置碱基添加到比对中
• 光标 2 向前移动，光标 1 保持不动，将序列 2 中的光标所在位置碱基添加到比对中

参考下图，对于每一步移动，我们会将该步移动的得分加上之前所有步骤的得分，设置成该位置的得分。

我们可以按照任何喜欢的顺序来计算得分——沿着矩阵的行、列或者对角线——只要我们能够保证所需的上一个步骤的得分被计算在内就可以。

从字符串比对到图比对（String to Graph Alignment）

可能超乎你的想象，将序列比对到 DAG 图上，只会对（序列比对所使用的）动态规划算法增加很小的复杂性。最初的 POA 文章中将动态规划矩阵加上了 3D 的“突起”（来表示节点有多个上游或者下游的情况），这可能导致问题看上去比实际更加复杂。

字符串比对和图比对的动态规划算法的主要不同在于，一个碱基在序列中只有一个上游，而在图中可能有两个或者更多个上游碱基。这样，图比对中光标移动的时候就要考虑上一个步骤中光标有可能来自多个上游中的一个，计算分数的时候要同时考虑这些来自不同上游的分数。

再重申一遍，图比对和字符串序列比对的唯一区别就在与图比对的动态规划算法循环中要考虑多个上游分数的情形以选择合适的插入或者碱基匹配。这就需要为当前碱基增加一个循环来遍历其所有上游碱基，算法的其他部分保持不变。

拓扑排序（Topological Sort）

在开始动态规划循环之前，还有一个需要考虑的步骤。

当对两个序列进行比对时，可以顺序的来对序列的索引进行遍历，这样对于任何正在计算的新位置，必要的上游分数就已经确定了。

但是对于图上的节点来说，并没有序列这样的固有顺序关系。例如，如果考虑节点插入图中的先后顺序，那么插入一个带有新序列的最新节点的时候——该节点有可能成为之前插入的某个节点的上游节点——下游节点开始计算的时候该节点的分数可能还没有准备好。

因此，使用拓扑排序来生成节点的顺序关系，这样可以保证每个节点都在其确定上游节点之后。对于一个有向无环图来说，这样的顺序关系总能找得到，而且事实上可能还能找到不只一套这样的顺序关系。拓扑排序是 make 和类似的工具如何决定在工作流中执行任务的顺序，以及有多少个⁶电子表格程序决定单元格是否需要更新（译者注：这半句不是很懂 … …）。

拓扑排序有两类主要的算法，Kahn（1962）算法和重复深度优先搜索算法。两种算法均适合我们的动态规划算法问题。

如此以来，将一个序列比对到图上就变得简单了：

• 如果图更新了，执行拓扑排序
• 照常运行动态规划步骤：

• 以拓扑排序顺序访问图上的节点，
• 考虑对齐/插入移动的所有有效上游节点

序列的插入（Insertion of aligned sequence）

Consider that we have a graph that so far only contains the sequence CGCTTAT, and our dynamic programming calculation aligning the sequence CGATTACG has given us an alignment that looks like this:
假设我们现在有一个只包含序列 CGCTTAT 的图，利用动态规划算法将该序列和序列 CGATTACG

C G A T T A C G
| | . | | | .
C G C T T A T -

也就是说，序列中的每个碱基，要么和图中的一个碱基进行配对（可能是匹配或者错配），要么就将其插入图中。

我们期望将新的序列插入图中会得到如下图这样的结果：

在上图中我们第一次见到有两种类型的边，加粗的边是有向的（方向未在图中标注，默认未从左到右），表示节点的上下游关系；从左数第三个节点 A 和 C 之间的虚线边表示这两个碱基被比对到了一起，但是是个错配（mismatch），快到结尾处的 C 和 T

我们来同时跟踪上下游节点和所有比对到一起的节点。沿着我们正在插入的序列往前并进行比对计算。如果序列中的某个碱基没有比对到任何地方，或者他直接或间接比对上的节点都没有相同的碱基，那我们就将这个碱基插入图中形成一个新的节点，否则我们直接使用已有的节点并在必要时添加一条边。

更具体的说，我们采取的步骤如下：

• 首先在图中为序列创建一个起始节点
• 将前置位置（previous position）设置为这个起始节点
• 对于比对计算的序列中的每个碱基

• 如果当前的碱基没有比对到图中的任何一个节点，或者该碱基比对到的节点以及其他的节点都不和该碱基相同

• 那么就为这个碱基在图中创建一个新节点，并将该节点设置为当前节点
• 将这个新节点与其比对到的节点（如果有的话）对齐（译者注：按照上一个图片的说明，这里应该是在这两个节点之间添加一条虚线边），其他所有比对到的节点也更新为与这个新节点对齐（译者注：应该也是增加新的虚线边）

• 否则，

• 当前碱基比对到的节点被设置为当前节点

• 如果在前置位置和当前节点之间没有边，那么新添加一条（译者注：这个边应该是加粗边，有方向，从前置位置指到当前节点）
• 将当前序列的标签添加到（前置位置和当前节点之间的）边上；该条边上的标签数依赖与包含这条边和两个节点的所有序列的数目

尽可能的对节点进行融合，以保证 motif 的信息，即在几个序列中类似位置的多次出现的 motif 的信息不被图中的多个路径所掩盖（译者注：这句没太看明白是什么意思）

值得一提的是，通过找到插入图中某序列的起始节点，并沿着标有该序列标签的边遍历图，可以重构插入图中的任何单独序列。

一旦被比对的序列完成插入，并且对图中的节点重新进行拓扑排序，就可以开始进行下一个比对。

公共路径（Consensus paths）

现在你所有的序列都已经被插入到图中，那么怎样从图中得到类似像公共序列（consensus sequence）这样的信息呢？这是另一篇论文²中讨论的内容。

在图中找到单个的最佳路径相对比较简单。事实上，这仍然是一个动态规划问题。现将所有的节点分数置为零，并逐节点对图进行遍历。对于每个节点，选择最佳的边进入该节点——也即被最多序列所包含的边——并且将边的权重和其指向的节点分数相加来得到打分；如果有多条边平行，选择指向得分最高节点的边。

通过最高得分和边的选择会得到图中的最大权重路径。论文中指出，如果边的权重依赖与该边被遍历的概率，那么这就是一个最大似然路径。

然而，你也许想从比对中提取多个公共序列，该序列出现了多次但仍然是一个相对少数的序列（译者注：可能是相对于最佳路径而言）。寻找剩余公共序列的方法需要基于启发式的算法，该论文中有大部分篇幅讨论这个问题。

最基本的想法是去除已经提取的公共序列的边或者降低其权重，然后重复上述发现最佳路径的步骤来寻找其他的公共序列。在该论文中建议的步骤为：

• 通过部分的边或者节点，以及一些其他可能的条件，查找到对应与公共序列的那些序列
• 把和这些序列相关的边的权重降低，一种可能是将权重降至零
• 继续前述的公共序列查找算法

在描述该算法的简单实现里，我们简单的选取了所有大部分碱基都在当前公共序列中的序列，并将相应的边的权重全部去除，重复这个步骤直到没有剩余的序列或者无法找到显著的公共序列。

该论文发现了一中极端的情形，这种情形下公共序列会提早终止，我们允许这种情况的发生。

比对字符串（Alignment strings）

最终，为了交流比对结果，生成一个输入和公共序列的“扁平”比对结果仍然是有用的。

一旦对图进行了拓扑排序，这也是相当简单的。图中的每个节点会按照拓扑排序的顺序，被分配到最终生成表的一列中，彼此对齐的节点被分配到同一列，没有比对到其他节点的节点单独一列。然后把碱基填入，每一条序列（包括公共序列）一行。

因为我们是按照拓扑排序的顺序来将节点分配到列中，那么生成拓扑排序（非唯一）的方法就会影响到比对结果最终的字符串表示形式，即使他们是功能上等同的。Kahn 排序会倾向于将结果序列交替排列，而深度优先排序则必然会按照顺序访问长字符串。这样的话，深度优先生成的比对字串更易读，因此我们在后面的实现中采用的是这种方法。

简单实现

可以在这里找到上述算法的一个简单但是功能完整的实现。在比对阶段，我们给出了两种实现；一种非常容易理解但是运行速度很慢；另外一种用到了 numpy 的向量化来提高效率，这显著的提升了速度，但是读起来更花功夫。

即使是更快的那个实现版本运行起来仍然很慢——比用 C 语言实现的 poaV2 要慢 10 倍，如果 poaV2 使用 -O3 选项进行优化编译，则要慢上 20 倍——但不管怎么说它对于小规模的问题还是有用的。

上述的简单实现可以生成一个 HTML 文件，可以交互式的方式图形化的浏览最终的偏序图；可视化工具在有高性能实现的浏览器上工作良好，但是不支持超过一千个节点的图。

结论

偏序图比对是一个非常有力的工具，它能够生成包含序列比对结构的丰富信息的图，但是相比于其他方法，缺少相关的在线文档和便于浏览的实现。我们希望这个工作能有助于更为广泛的受众更深入的理解这个方法。

References

1. Multiple sequence alignment using partial order graphs (2002) by Lee, Grasso, and Sharlow
2. Generating consensus sequences from partial order multiple sequence alignment graphs (2003) by Lee
3. Combining partial order alignment and progressive multiple sequence alignment increases alignment speed and scalability to very large alignment problems (2004), Grasso and Lee
4. Multiple alignment of protein sequences with repeats and rearrangements (2006) Phouong et al.
5. Cactus: Algorithms for genome multiple sequence alignment (2011) Paten et al.
6. Later versions of excel actually allow circular dependencies in cell calculations.