对基因组拼接算法的深入分析和探讨

张磊

摘要：研究基因组拼接算法，进而更高效地实现全基因组拼接、获取生物体遗传信息，对于生命科学研究具有重要的意义。 拼接之前对数据进行一定的预处理。拼接采用经典的DBG算法，分析了DBG算法的实现步骤以及需要注意的问题。根据read切割K-mer，然后采用位运算压缩存入Hash表。以有效的K-mer为顶点，相邻的K-mer之间连有向边，建立DBG图，然后在该图中寻找Euler路径。每一条Euler路径对应一条contig。最后进行参数组合实验，并将结果与专业软件CLC的结果进行对比。两者效率相近，说明了算法的高效性。

关键词：拼接；DBG图；K-mer；contig；scaffold

中图分类号：Q812；TP311 文献标识码：A 文章编号：1009-3044（2015）33-0123-02

快速和准确地获取生物体的遗传信息对于生命科学研究具有重要的意义。对每个生物体来说，基因组包含了整个生物体的遗传信息，这些信息通常由组成基因组的DNA或RNA分子中碱基对的排列顺序所决定。获得目标生物基因组的序列信息，进而比较全面地揭示基因组的复杂性和多样性，成为生命科学领域的重要研究内容。

本文深入探讨了基因组拼接算法的实现细节，为实现基因组组装算法小平台化提供了参考。

测试数据：一个全长约为120000个碱基对的细菌人工染色体（BAC）， 采用Hiseq2000测序仪进行测序，测序深度（sequencing depth）约为70×。

1 背景介绍

基因组组装是指以基因组测序所得的序列片段read数据为输入，利用序列片段之间的交叠关系，再加上配对数据、碱基质量等辅助信息，通过一定的方法手段重建目标基因组序列的过程。基于新一代测序数据的拼接算法，通常包括如下三个阶段：

1） 数据的预处理：主要是对测序数据中的碱基进行纠错；

2） 拼接contig：短read片段被拼接成长度较长的contig；

3） 组装scaffold：利用配对数据，确定contig之间的相对位置关系和距离，对contig进行组装。

2 拼接算法

本文讨论的是第二个阶段——拼接contig。

2.1 DBG算法简介

目前，基于DBG图策略的拼接算法被最广泛应用到新一代测序数据的处理中，比如华大基因的SOAPdenovo [5]。基于DBG图的拼接算法，非常巧妙地将具有交叠关系的read映射到一起，降低了计算交叠时的复杂度，减少了内存消耗。本文也采用DBG算法进行拼接。

在本文中，主要用到的经典算法有Euler路径套圈算法，以及Hash表的插入算法。读者可以参考相关的算法书籍。

DBG图以K-mer为顶点，相邻的K-mer之间连有向边，配对文件中的任一条read对应一个K-mer路径，将寻找contig的问题转化为在该图中寻找一条Euler路径。

2.2 拼接算法大致流程

笔者实现的DBG算法，具体步骤如下：

S1：用户输入参数：K-mer的大小、K-mer有效的最低频次、contig有效的最短长度。

S2：读入并保存预处理后的文件中的read，并生成对应的反向互补串。

S3：Hash表初始化。根据指定的K-mer，将所有的read及其反向互补串连续地切割成指定长度的串，相邻的两个K-mer交叠K-1个碱基。利用位运算压缩存储，插入Hash表。数组p[c][i][j]记录该K-mer对应的去重后编号。

S4：根据数组p[c][i][j]构建数组opp[]。

S5：将Hash表去除冗余，存入数组kmer[]。

S6：根据数组p[c][i][j]构建DBG图。如果去重后编号为p[c][i][j]和p[c][i][j+1]的K-mer都是有效的，则在两个K-mer对应的顶点之间连一条有向边，由p[c][i][j]指向p[c][i][j+1]，记录在数组map[][]里面。同时记录每个顶点i的入度数组ind[]和出度数组outd[]。

S7：从每个有效而且ind[i]=0，outd[i]>0的顶点出发，寻找Euler路径，每一条Euler路径对应一个contig。寻找的时候，总是出现频次最大的后继K-mer对应的路径优先走。每经过一条边，便在DBG图中删除该边以及相应的反向互补串对应的边。每找完一个contig，记录其长度，并根据contig对应的K-mer序列，输出该contig。

S8：根据各个contig的长度，并计算N50以及拼接序列总长度。

2.3 拼接算法的细节处理

关于DBG算法，需要说明以下几点：

1）S2中生成反向互补串，这是由DNA实际结构决定的：DNA的两条链之间的关系为反向互补。

2）S3中位运算采用二进制编码，碱基与碱基编码对应关系为：A-00，C-01，T-10，G-11。采用二进制编码，将存储字符转化为存储二进制数，大大减少了内存的消耗。一般来说，对整数的位操作远远快于对序列碱基的字符操作。其中，K-mer的第k个碱基对应的int的下标p1计算如下：

[p1=ints-[（K-mer-k+1）×2-1]/bits] 其中：ints为用来存储K-mer的int数据的个数，bits为一个int用来存储碱基信息的位数（int数据的最高位为符号位，所以bits最多为30）。式中除号为整除。

K-mer的第k个碱基对应int的数据位p2，p2+1，p2计算如下：

[p2=bits-[（K-mer-k+1）×2-1]%bits]其中：%表示求余运算。

3）S3中采用的Hash函数为：

[p=（i=1intsbit[i]/ints）%hashsize]

其中：bit[i]为该K-mer对应的第i个int数据，hashsize为Hash表的大小。

4）经过实验，采用上述Hash函数时，当Hash表的大小为最大可能值的3-4倍时，冲突便可以控制在理想的范围以内。

5）S3中p[c][i][j]，c=1表示原串，c=2表示反向互补串，i表示read序号，j表示K-mer在read中第一碱基的位置。该数组是为了计算数组opp[]和构建DBG图。

6）S4中opp[i]表示去重后编号为i的K-mer的反向互补串对应的去重后编号。

7）S6中map[i][j]，i表示某个K-mer的去重后编号，j为某个碱基对应的二进制编码。若map[i][j]=0，表示去重后编号为i的K-mer后面不可能紧邻j对应的碱基。若map[i][j]>0，则表示去重后编号为i的K-mer去掉第一个碱基，后面紧邻j对应的碱基形成的K-mer的去重后编号。

8）S7中条件ind[i]=0和outd[i]>0保证了i一定是某个read（包括反向互补串）的初始K-mer。

9）S7中选择出现频次最大的后继K-mer对应的路径优先走，是因为出现频次越大，该K-mer正确的可能性越大。这样可以有效地避免遇到错误碱基而导致的拼接中断。

10）S7中如果不删除相应的反向互补串对应的边，可能会出现两个contig或者其子串之间存在反向互补关系的现象，而这样显然是冗余的。

11）本文中未考虑重复片段的影响。位于重复区域的read会被集中拼接到一个区域，而在另一个区域却无法拼接。

12）基于DBG图的拼接算法的一个关键参数是K值，一般认为对于原核生物的基因组拼接，K值在21-35之间是合适的[6-7]。提供的测试数据为某细菌的，而细菌属于原核生物。

4 结论

目前的几乎所有组装软件都需要较大内存的计算平台，但是本文的程序在普通的PC机上即可实现，具有很大的可推广价值。但是，对测序错误导致的tip、bubble结构以及repeat的处理，没有进行深入考虑。另一方面，需要用户进行参数组合实验，不能自动计算最优参数组合。如果进一步进行研究，可以从采用多种配对文库、长序列与短序列混合、并行化拼接等方面进行考虑。

参考文献：

[1] Lucian I， Farideh F， Silvana I. HiTEC： accurate error correction in high-throughput sequencing data[J/OL]. （2010-07-28）[2012-05-28].

[2] Jan S， Heiko S， Simon J P， et al. SHREC： a short-read error correction method[EB/OL]. Oxford：OxfordJournal，2009[2012-06-03].http：//bioinformatics.oxfordjournals.org/content.

[3] Yang X， Dorman K S， Aluru S. Reptile： representative tiling for short read error correction[EB /OL]. Bioinformatics， 2010， 285： 313-314[2012-06-10].

[4] Hernandez D， Francois P， Farinelli L. De Novo Bacterial Genome Sequencing： Millions of Very Short Reads Assembled on a Desktop Computer[J]. Genome Research 2008， 18（5）：802-809.

[5] Zerbino DR， Birney E. Velvet： Algorithms for De Novo Short Read Assembly Using De Bruijn Graphs[J]. Genome Research 2008， 18（5）：821-829.

[6] 赵玮. 基于配对信息的 DNA 从头测序组装算法研究[D]. 哈尔滨： 哈尔滨工业大学， 2011.

[7] 韩东涛. 基于概率模型的基因组从头测序算法研究[D]. 哈尔滨： 哈尔滨工业大学， 2012.

[8] 曾培龙. 基于 reads 引导的基因组序列拼接[D]. 哈尔滨： 哈尔滨工业大学， 2012.

[9] 黄勇. 基于高通量测序的微生物基因组学研究[D]. 北京： 中国人民解放军军事医学科学院， 2013.

[10] 王东阳， 任世军， 王亚东. DNA 序列拼接中 de Bruijn 图结构的研究[J]. 智能计算机与应用， 2011， 1（2）： 20-25.