基于德布莱英图（De Bruijn graph）的宏基因组de novo拼接

最近在做一批宏基因组数据的拼接工作，这里的拼接主要指从Read到Contig这个水平。然后因为是第一次做，有很多没懂的地方，在学习和实践的过程中，总结了一下。写成几篇博客记录一下，以后又不清楚了可以回来看一下。这篇博客主要记录一下宏基因组拼接的算法之一——基于德布莱英图（De Bruijn graph）的方法。

一、为什么使用德布莱英图？

对于短序列read的拼接有两种基本算法，一种是基于Overlap Graph的算法，另一种是基于德布莱英图的算法。大多数为了桑格测序（Sanger sequencing）得到的测序数据拼接的算法都是基于Overlap Graph的，最典型的就是Overlap-layout-consensus方法，简称为OLC方法，这个方法通过两两之间有overlap区域的reads来构图，每个Read作为图中的点，而拥有足够的overlap则是建立图中的连接的基础。

这种方式在处理序列数很有限，而且序列间overlap区域很明显的时候能取得最好的效果。一些对于二代测序数据设计的assembler使用了OLC方法，但是这种传统的方法计算量实在过大，就算是非常简单的物种的从头组装也需要上百万甚至上千万条序列，所以仍然使用OLC方法显然不太实际。而且二代测序数据的读长片段，现在来说250bp和300bp的读长都是属于新兴事物，想要在茫茫多的短序列中找到可供使用的overlap区域也是OLC方法面临的瓶颈之一。

由于OLC方法跟不上更短的Read长度和更大的Read条数，更多的为二代测序数据拼接的assembler都选择了基于德布莱英图方法的，本文中将其简称为DBG，指de Bruijn graph。DBG方法通过将序列打断成更小的长度为k的小片段（称为k-mer）来降低计算的规模。对于一条长度为N的序列从头到尾（比如从5’-端到3‘-端）依次取k-mer一共可以取到*N-k+1个k-mer，每两个相邻的k-mer之间都有长度为k-1的overlap区域。

通过将序列打断为k-mers，德布莱英图方法降低了大量短序列数据集中的数据冗余性，就好比上图中的这条序列在被切为3-mers后，GAT这个k-mer是重复出现的，并在德布莱英图中被反复利用，也就达到了降低了整体数据量和计算量的目的。对于某一特定的组装来说，最大的有效k-mer size取决于读长和error rate，同时如何选择k的取值也是决定了组装质量的重要条件（我现在也正在调查研究这一点）。实际上，在拼接前是有方法可以先看看合适的k值的，比如kmergenie这款软件就可以给出一个预测的k-mer值。

德布莱英图的另一个很优秀的特性就是在基因组中的一些重复片段会在德布莱英图中被打断，避免了一些错误的组装结果（当然不能很好地解决）。

第一节的内容和图片整理来自：https://banana-slug.soe.ucsc.edu/_media/bioinformatic_tools:abyss_technote_illumina.pdf

二、从德布莱英图到拼接结果

对于某一组测序数据来说，通过相应的算法获得了德布莱英图，怎么将图转化为序列呢？首先我们必须知道，德布莱英图是具有欧拉图（Eulerian Graph）的性质的，也就是说我们需要寻找的通路（path）是一条遍历每个边至少一次的。当然我们也希望它是遍历每个点至少一次，但是哈密顿图（Hamiltonian Graph）作为世界难解的NP-Complete问题，还很难被归化成一个具有概括性的模型。

根据上面的欧拉通路性质，对于下面的这组德布莱英图，我们认为它的拼接结果应该是GATTACATTACAA，或者说图中的ACA-CAT-ATT这个分支环路至少要被遍历一次。

事实上，当序列被打断为短k-mer后，序列的模糊性就变高了，就比如上图中的分支环路到底遍历几次，或者像下图这样的：

好像有很多条通路都可以满足条件。所以在实际情况中，assembler返回的是一组子通路（substring）的集合，通过这些子通路就可以组装出多条通路。比如针对上面这张图，返回的子通路集合是这样的{ACTGA, GACC, GAGTG, GAATG}，其中第一个元素代表这个图的“入口”，第二个元素好代表这个图的“出口”，第三个元素代表的是上面蓝色的分支环路，第四个元素代表的是下面绿色的分支环路，通过按顺序选择分支就可以得到一条完整的路径。

通过拼接得到的contig是不能有分支节点的，也就是说我们想要得到的就是简单的单一路径，或者说单一序列。所以如果我们最终得到的contig是下图这样的，我们就会把它们按分支节点拆开，作为最终的contig结果。

这一部分其实还有大佬讲的更详细一点，参考：https://zhuanlan.zhihu.com/p/54466660。

三、几个值得关注的问题（之后边学习边写）

以上就是我现在比较关心的从read到contig的这样一个拼接过程的简单流程，但是在实用阶段中还有如下几个问题：

1. 以上介绍的都是单k-mer实现序列拼接，但是现在的Assemblers普遍输入k-mer list，也就是说它会有多个k值进行迭代式的拼接。为什么要用multi-k-mer的方法？k值以及k-max、k-min该如何选择？
2. Assembler的参数该如何选择？或者说该如何针对不同的数据情况选择？
3. 德布莱英图只是一个基础，不同的Assembler，比如Megahit、MetaSPAdes、IDBA-UD等都在DBG的方法上进行了自己的改良优化，它们的详细原理是什么（关系到它们的适用场景和参数选择）？
4. 还有很多使用中遇到的问题……

在学习的过程中写了这些笔记，整理出来发成一个文档集。之后应该还会写以上几个问题的学习笔记和几款Assembler的使用说明（弄清楚了咋整的之后……），stay tuned :)

参考文献和学习资料：

1. https://zhuanlan.zhihu.com/p/54466660

2. http://www.cs.bilkent.edu.tr/~calkan/teaching/cs481/pdfslides/20-cs481-intro_to_assembly.pdf

3. http://compbio.charite.de/tl_files/groupmembers/robinson/genome-assembly-1.pdf

4. http://compbio.charite.de/tl_files/groupmembers/robinson/genome-assembly-2.pdf

5. Pasolli E, Asnicar F, Manara S, et al. Extensive Unexplored Human Microbiome Diversity Revealed by Over 150,000 Genomes from Metagenomes Spanning Age, Geography, and Lifestyle. Cell. 2019;176(3):649-662.e20. doi:10.1016/j.cell.2019.01.001

6. Forouzan E, Shariati P, Mousavi Maleki MS, Karkhane AA, Yakhchali B. Practical evaluation of 11 de novo assemblers in metagenome assembly. J Microbiol Methods. 2018 Aug;151:99-105. doi: 10.1016/j.mimet.2018.06.007. Epub 2018 Jun 25. PMID: 29953874.

7. https://banana-slug.soe.ucsc.edu/_media/bioinformatic_tools:abyss_technote_illumina.pdf