基因组组装算法

本文主要内容：

• De Bruijn graph算法
• 基因组组装过程
• Kmer为什么是奇数
• Kmer的作用

在上篇文章cutadapt使用指南中，我们提到可以使用minion来寻找测序数据的接头，但是原理是什么呢？这个问题激发了我的兴趣，于是搜索了原文献：

Unfortunately sometimes adapter information is not properly tracked and attached as metadata to the raw sequencing data and may not be known for a given sample. To alleviate this problem, we have developed Minion, a C-program that uses De Bruijn graphs to detect sequence fragments with characteristics typical of adapter sequence such as high abundance towards the 3′ end of reads.

可是De Bruijn graphs是什么东西呢？参考维基百科：

In graph theory, ann-dimensional De Bruijn graph of m symbols is a directed graphrepresenting overlaps between sequences of symbols.

也就是说，De Bruijn graph是一个展示符号序列之间重叠关系的有方向的示意图。举个栗子，比如说拼接基因组的序列：

De Bruijn graph

但是这样仍然不知道这在干嘛，可以理解为：

De Bruijn graph思路

看到这里，是不是冒出来很多问题？比如，基因组是怎么拼接的？De Bruijn graph在基因组拼接中做了什么？另一个是为什么是开头是AGT？我们来一个个解答。

基因组是怎么组装的？

基因组的组装流程主要可以分为：原始数据→数据过滤→纠错→kmer分析→denovo组装。下图只包括了后三个步骤，并标出了De Bruijn graph的作用。

纠错→kmer分析→denovo组装

无论是一代Sanger法，二代还是三代测序，目前，我们得到的测序数据相对于整个基因组而言都是极小的。就好比你要完成几十个足球场那么大的拼图，但是你已有的都是小拇指大小的碎片，并且存在很多重复的碎片。现在让你拼，你要怎么办呢？是不是脑阔痛又要掉头发了？

科学家们解决拼成上百万、千万、亿万个数据“碎片”（reads）的方法，就是通过碎片中重复的部分（即overlap）构建graph，利用算法从graph中得到最优路径，从而得到最初的contig。

那么多聪明的脑袋，办法多，他们也想了很多算法。目前主流的算法除了De Bruijn graph，还有Overlap-Layout-Consensus，String Graph等，这篇文章我们主要介绍De Bruijn graph。

现在你应该可以回答De Bruijn graph做了什么吧？就是在把“碎片”reads打断成长度为k的核酸片段（kmer）后，利用De Bruijn graph根据kmer间重复的部分构建graph，得到最优化路径从而拼接contig。

为什么是开头是AGT？

在看了更多资料后，发现其实这个问题可以简化为利用De Bruijn graph组装基因组的时候，Kmer为什么必须是奇数？有什么作用？

首先了解下到底什么是Kmer？

所谓Kmer，即为一段短的DNA片段。K为一个奇数，K等于几，就称为几mer。例如一段17bp的序列，则可以称为17mer。
我们可以将一条reads切割成很多小的Kmer片段，从第一个碱基开始，每隔固定距离的碱基开始提取碱基。例如一条100bp长的reads，每隔一个位置取一个17mer的片段。也就是1-17取一个kmer出来，2-18取一个Kmer，3-19取一个Kmer。以此类推，最终84-100为1个kmer。那么最终将会生成将会生成100-17+1，也就是84个Kmer片段。原来一条长片段，变成了很多短的片段，碱基的数目也增加了很多倍。
而且，每次取Kmer是同一条reads正反取两次，也就是对这条reads的反向互补序列再取一次Kmer。

那Kmer为什么是奇数？

根本原因是为了避免正反义链混淆。
目前的NGS测序技术也做不到通测基因组。一般来说都是测出上百万千万亿万个小小的片段（read，长度一般是100bp-300bp）。而且，为了确保准确性，基因组都会被反复测很多层。组装时构建的kmer单位，实际上是对这些read进行的。具体的操作就是按照kmer的长度把这些read切割成更小的、存在重叠关系的片段。那么，此刻当我们构建de Bruijn graph时，如何能够保证正确地把同属于一条read上的Kmer连接起来，就显得极为重要了！我们不能一会儿把A kmer正确地连到它自己所在的read，一会儿又连到它互补链的read上去！

这是什么意思呢？举个例子，4mer的ATAT反向互补后仍是ATAT，而3mer的ATA反向互补是什么？TAT。这样就可以知道这个kmer是在正义链还是互补链了对不~

Kmer分析有什么作用？使用reads拼接不行吗？

优点主要有两个：

我们将测序的reads切割成更短的kmer，形成了一个巨大的kmer的集合。所有的kmer片段长度相同，可以计算kmer出现的频数。也就是把相同的kmer加起来计数。只出现一次的kmer有多少个，出现两次的kmer有多少个，出现三次的kmer有多少。那么最终就会得到一个kmer频数的统计表，绘制一张kmer频数的分布图。其中有一个峰值，这个峰值就是kmer的深度。可以利用这个值估计基因组大小。

在对二代测序进行序列拼接时，为了去除测序错误同时又充分利用数据，所以就采取了将reads切割成kmer的方法。利用kmer可以用于去除测序错误碱基，得到更好的拼接效果。

 

参考：

1. Lecture 7: Assembly - De Bruijn Graph

2. Combining de Bruijn graph, overlap graph and microassembly for de novo genome assembly

3. 什么是kmer，kmer分析有什么用？

4. 利用de Bruijn graph组装基因组的时候，Kmer为什么必须是奇数？