Metagenome Assembly - Part2：[文献阅读] Why are de Bruijn graphs useful for genome assembly?

发表日期：2011年11月8日

发表期刊：Nature Biotechnology 29, 987-991(2011)

作者：

Phillip E. C. Compeau,

Department of Mathematics, University of California San Diego, La Jolla, CA, USA

Pavel A. Pevzner, and

Department of Computer Science and Engineering, University of California San Diego, La Jolla, CA, USA

Glenn Tesler

Department of Mathematics, University of California San Diego, La Jolla, CA, USA

将大量短序列拼接为连续的基因组是一道难题

针对二代测序数据的算法思路的发展可以追溯到300年前的普鲁士柯尼斯堡（Prussian city of Königsberg，现在的俄罗斯的加里宁格勒，Kaliningrad），那时的柯尼斯堡城被城中河分为4个部分，4个城区之间共有7座桥相邻接（如图1a所示）。柯尼斯堡居民非常享受在城中漫步，于是他们想到了如下的问题：有没有可能从一个城区开始逛遍整个城市，过程中经过每座桥梁恰好1次，最后返回到出发的城区。伟大的数学家莱昂哈德·欧拉（Leonhard Euler）在1735年对这个问题（我们把它称作经典的七桥问题，Bridges of Königsberg Problem）做出了概念上的总结和突破，同时这项突破也让今天的生物学上的短序列拼接成为了可能。

欧拉的思想首先将每个城区看做点（point，called a node），每座桥则看做点之间的连线（called an edge），这样就构成了图（graph）——由点和点与点之间的连接组成的网络（如图1b所示）。通过描述任意一个图中是否含有欧拉环路（Eulerian cycle，图中经过每条连接并回到起点的路径）的步骤，欧拉不仅解决了七桥问题，同时也标志了数学中相当重要的分支——*图论（graph theory）*的开端。

基于对齐拼接的计算问题

为了说明图对于基因组拼接的有利之处，我们会用从一个小的环形基因组（ATGGCGTGCA, 图2a）中取出的5个简单序列（CGTGCAA, ATGGCGT, CAATGGC, GGCGTGC, TGCAATG)来举例。现在的NGS方法生成的读长种类很多，最流行的读长约在100bp（本文2011年发表，现在的250bp、300bp也很普遍了）。将短序列拼接成长的连续序序列的一个简单直接的方法就是把每条短序列看做点，然后将有重叠区域（overlap）的点用箭头连接起来（箭头具有方向性，也说明了这样的图是有向图，directed edge）组成图，2001的人类基因组拼接与很多其他用Sanger测序方法得到的测序数据的拼接都是基于这样的方法。如图2b里展示的那样，至少有5个碱基重叠的两个序列就会被箭头连起来。

假设现在有一只蚂蚁在沿着图上的边爬行会更方便我们理解图算法。在基因组拼接的过程中，蚂蚁爬过了的路径就是可能的基因组拼接结果。特别地，如果小蚂蚁爬过的路径是：ATGGCGT → GGCGTGC → CGTGCAA → TGCAATG → CAATGGC → ATGGCGT，那么它的路径在图中就构成了一个哈密顿环路，它指的是一类经过每个点恰好1次并最终返回起点的路径，也就是说在基因组拼接中，每条序列只会被使用一次。

现在的拼接软件是利用一系列k-mers来工作的，k-mer通常来说比read序列短。100bp的序列可以依次生成46个55-mers。接下来我们可以用这些k-mers来寻找哈密顿环路：首先，对于一套reads来说，首先将所有序列生成的k-mer作为节点，然后对于每个k-mer，规定它的*前缀（prefix）是除了最后一个碱基之外的前面k-1个碱基序列，它的后缀（suffix）*是除了第一个碱基之外的后面k-1个碱基序列。如果前一个k-mer的后缀和后一个k-mer的前缀相同，那么就将这两个k-mer的节点连接起来。之后，寻找一个可以代表可能基因组的哈密顿环路，这个环路使用了每个检测到的k-mer，并且具有最小长度，因为它遍历每个k-mer仅1次。

然而，这个过程并不如上面描述的这么简单。且不谈Illumina产生的数据量之巨大，现在对million级（更别说billion级）数量的节点的哈密顿环路整理并没有有效的算法。哈密顿方法的计算压力实在是过高，以至于基于NGS的数据项目都不得不放弃哈密段环路方法。事实上，这确实是计算科学给基因组序列处理的设下的一个限制：寻找哈密顿环路在实际上属于计算科学中的一类叫做NP-Complete的问题。时至今日，世界上顶尖的计算机科学家也没有找到一个很好的能解决这类问题的方法。

利用德布莱英图的拼接

前面我们已经提到了，想找到遍历每个节点恰好1次的哈密顿环路比较难，但是我们马上就讲到，想要找到遍历每个连接（边）恰好1次的环路是相对来说简单的。这为我们的拼接提供了新的思路，我们不再把k-mer考虑成图中的节点，而我们把在read中定位的k-mer考虑为图中的连接，这就构建了德布莱英图（de Bruijn graph），简写E。首先，对于每个k-mer的前缀和后缀，把它们考虑为节点，也就是说每个给定的k-1长度的序列只能在图中作为节点出现仅1次；然后，如果某些k-mer具有前缀x和后缀y，那么就把节点x和节点y用有向箭头连接起来，然后将这个连接标作该k-mer（如图2d所示）。

现在我们再想象小蚂蚁的路径，它不再是遍历所有点，而是尝试去遍历图E的每条边。这个问题听起来熟悉吗？没错，这实际上就把问题转化为了七桥问题，而我们要寻找的环路就是欧拉环路。由于遍历了每条边，也就是遍历了所有k-mer，那么生成的环路也就是可能的基因组序列。

欧拉证明了：对于任何无向连通图（图中任意两点间都至少有一条路径），仅当图中每个节点的度（节点具有的连接数）为偶数时，该图中就一定有欧拉环路。至此，我们可以知道，七桥问题的答案是否定的，因为该问题中的4个节点的度全是奇数。

上述结论在有向图中也是类似的。对于有向图中的任意一点v，定义指向v的连接的数目为v的入度（indegree），由v指向其他节点的连接的数目为v的出度（outdegree）。如果某一有向图的所有节点的出度和入度相等，则称该图为平衡（balanced）。**欧拉的理论说明有向连接图当且仅当它是平衡时存在欧拉环路。**特别地，该理论还可以进一步说明。对于德布莱英图E，只要我们找到了所有可能的k-mer，它就一定是存在欧拉环路的，因为每个k-mer都是由独特的前缀节点和后缀节点组建的，对于每个节点来说，他的出入度就表示了(k-1)-mer作为节点被使用的次数。

其实很容易知道一个含有欧拉环路的有向图是平衡的，因为每当小蚂蚁遍历欧拉环路的某一特定的顶点，它就走上了环路的一条边并从下一条边离开，这就将接到每个顶点的两条边配对起来了，说明指向某一顶点的边和从该点离开的边数目各占一半。想看有向图中有没有欧拉环路就相对难一些。为了证明，欧拉也派了小蚂蚁去随机探索一张图，小蚂蚁的旅行只有一条限制，它不能走已经走过的边。那么，小蚂蚁迟早会被困在某一个节点，此时此刻它面前的所有边都是之前被走过的，欧拉注意到以为图是平衡状态的，这个“禁止通行”的节点正是小蚂蚁旅行的起点。这就说明小蚂蚁已经完成了一个环路的遍历，如果这个环路正好遍历了所有边，那么相当于小蚂蚁找到了一条欧拉环路。否则，欧拉就会再派一只小蚂蚁取随机遍历前一只小蚂蚁没有走过的边，这样来找到第二条环路。之后，欧拉又说明了，两只小蚂蚁走过的环路实际上可以合并为一个环路，如果这个环路包括了图中所有的边，那么两只小蚂蚁就合力找到了一个欧拉环路。如果还不行，欧拉就还会派出更多蚂蚁……最终，欧拉环路会被寻得。

在现在的计算机上，这样的算法可以很有效地在含有亿万节点的图中找到欧拉环路，这样也避开了由NP-Completeness带来的困难。像这样简单地将一个问题的框架稍加改动就能将序列拼接问题转化为一个更容易处理的问题，这也是在计算科学中常常用到的策略。

计算机完成的欧拉算法过程运行时间和德布莱英图中的节点数目大致成比例。在哈密顿方法中，这个运行时间很可能会长很多。

不过，德布莱英图也不是万能的。通过我们的说明，做出了几个假设（还需要很多工作来证明和解决）。不过对于序列拼接的很多应用，利用德布莱英图的一些类似形式将涉及哈密顿环路的问题转化为欧拉环路问题解决的案例已经非常丰富且颇有成效。此外，德布莱英图的一些衍生方法也已经用于其他生物信息学问题的解决，如抗体序列问题，synteny block重建、RNA拼接等。

随着新的测序技术的出现和进步，对于基因组组装的最好方案也可能会发生变化。影响算法选择的因素包括数据量（读长、覆盖度）、数据质量（错误率）、基因组结构（重复区段的数目与大小、GC含量）。短序列测序技术产生的大量reads，正好有利于德布莱英图的使用。由于基于overlap方法的拼接算法分辨不了长于读长的重复区段，德布莱英图也可以很好地解决这类重复区段问题。然而，如果未来的测序技术会生成读长长而序列数少的高质量read，组装算法选择的钟摆可能又会摇回到基于overlap方法了（是不是预言了三代测序的长序列！？算是神预言吗？）。

译者：阅读和翻译感想

这篇文章最初是国科大的生物信息学课上，王秀杰老师课件上推荐的一篇文献，当时没仔细看，就是按老师说的理解了一下基因组组装这件事。自己开始真的做拼接了才知道这个背后的算法细节和过程非常复杂，其实单看这篇文章也只能了解大概，不过可以为脑子里构思图的构建画一张蓝图。其实读过了几个Assembler的文章之后，发现它们的算法虽说是基于德布莱英图，但也各有各的特点和改进优化。真正困难的是，对于不同类型的数据，同一Assembler设置不同参数，怎么才能让assembly有最好的结果？从理论上来讲，这些结果需要对算法的充分理解才能预测。

本文仅为对原文的翻译，翻译水平有限，如果感兴趣欢迎评论或者参阅原文。

Reference

Compeau, P., Pevzner, P. & Tesler, G. How to apply de Bruijn graphs to genome assembly. Nat Biotechnol 29, 987–991 (2011). https://doi.org/10.1038/nbt.2023

本文同步发表在本人个人博客上，链接：Metagenome Assembly – Part2：[文献阅读] Why are de Bruijn graphs useful for genome assembly?