基于reads引导的基因组序列拼接算法

曾培龙

摘 要：针对新一代测序技术数据读取片段reads长度短、准确度低、数据海量等特点，本文提出了基于reads引导的基因序列拼接算法（SRGA），以整条reads为拼接单位，并首次提出了基于数据特征和拼接信息累计的评分机制。选取常用测试集，将本文中的算法与序列拼接领域中的经典算法进行对比和分析，取得了较好的效果。

关键词：生物信息学；新一代测序技术；基因组序列拼接

中图分类号：TP391 文献标识码：A 文章编号：2095-2163（2015）03-

GENOME ASSEMBLY GUIDED BY READS

ZENG Peilong

（China Ship Development and Design Center， WuHan 430064，China）

Abstract：Due to next generation sequencing data of mass， short length and relatively low precision， this paper proposes a new genome assembly guided by reads， regarding one entire reads sequences as assembly unit. This algorithm firstly invents a scoring mechanism based on accumulated assembly information and data charactistics. Then the paper gives the metrics results of several algorithms on the test set， the proposed （SRGA） and several classical algorithm of genome assembly. Experimental results show SRGA can obtain satisfactory stereo matching results.

Key words：Bioinformatics， Next-generation Sequencing， Genome Assembly

0 引 言

新一代测序技术促进了生命科学的快速发展，但其产生的基因读取片段reads具有长度短、准确度低、数据海量等特点[1-2]，这就对序列拼接算法提出了相当严峻的挑战，传统的序列拼接软件已不再适用[3]。为此，即需针对新一代测序的数据特点，从实际应用需求出发，研发新的优质高效的序列拼接软件。

本文针对新一代测序数据的数据特点，提出了基于reads引导的基因组序列拼接算法（SRGA），并以整条reads为拼接单位，首次提出了基于数据特征和拼接信息累计的评分机制，从而减少不必要的重复计算，同时也提高了基因组序列拼接的质量和速度。

1 reads数据预处理

在序列拼接过程中，reads数据的预处理具有重要意义。由于新一代测序数据精确度较低，以及数据海量，造成reads中含有大量的碱基错误。理论上讲，De Bruijn图构建时，图的大小规模只与基因组的大小相关，与reads数据量无关。但在碱基错误的影响下，De Bruijn图[4]实际大小会随着reads数据的增加呈现几何型增长。拼接时，reads中的碱基错误也极易导致错误拼接。因此，拼接前需要对数据进行预处理，去除初始数据中的错误碱基。

（1）新一代测序数据的准确率较低，错误碱基主要分布在reads 3端，并且越靠近3端错误率越高，reads 3端错误率更高，接近20%，而在5端则非常准确[5]，如图1所示。为降低错误碱基对拼接的影响，拼接前需过滤掉出错率较高的碱基数据。处理方法为：以靠近3端二分之一reads长度的碱基序列为基准，计算该区域碱基序列的质量平均值，若该值小于15，则过滤掉该条reads。该平均值对应碱基的错误率，计算公式为：

（1）

其中，Q为碱基质量值， 为碱基出错率。

图1 Solexa数据错误分布

Fig.1 Error display of Solexa data

（2）测序过程中往往会产生许多人工数据[6]，这些reads数据会有许多标识为A的碱基序列，需要去除。处理方法为：若某条reads中A含量>=0.9，则该reads被过滤掉。

（3）测序过程中有时会产生一些没有被测出来的通常表示为“N”或“.”未知碱基[7]，需要去除。处理方法为：如果某条reads中含有未知碱基，则该条reads被过滤掉。

2 De Bruijn图构建

SRGA以reads为输入，经过数据预处理、De Bruijn图构建、长序列片段contigs拼接等，生成基因组序列。以reads为拼接单位，考虑所有正在参与拼接的reads信息，采用基于数据特征和拼接累计信息的评分机制，选择最佳的路径进行contigs扩展，直到contigs扩展结束。

2.1 De Bruijn图的建立

建立De Bruijn图时，为了节省存储空间及便于信息的快速查找，分两次解析数据文件，包括reads的碱基序列和碱基质量数据。De Bruijn图只存储碱基序列，碱基质量数据用来数据预处理。Reads上的连续子碱基串称之为kmer，通常用K值表示kmer的长度。

第一次解析reads初始数据文件：统计kmer数量，构建无边的De Bruijn图。依次读取reads数据文件中的每条reads上的kmer，通过哈希函数将kmer映射到De Bruijn图中，并更新图中kmer的数量。

第二次解析reads初始数据文件：存储kmer碱基序列，并更新kmer位置、标识等信息。依次解析数据文件中的每条reads， 并取得相应的kmer碱基序列。对于每个kmer，首先计算其在De Bruijn图中对应的地址，然后将其碱基序列存储在对应的数组中。

2.2 De Bruijn图的基本操作

序列拼接时，De Bruijn图的主要操作包括kmer的查找、kmer的删除、kmer的更新、kmer所在reads数据信息的更新等。kmer的碱基序列在De Bruijn图中采用位存储的方式存储，其查找、删除、更新等与位操作一一对应。拼接时，需要将kmer的碱基序列经过二进制编码转换成整数，碱基序列所对应的整数经过位操作也可以转换成kemr的碱基序列。位操作的方式一方面节省了kmer碱基序列存储所需要的内存空间，另一方面也提高了kmer的查找、删除等操作的效率。

3 contigs拼接

拼接过程中，综合考虑kmer的碱基质量、拼接信息累计等各项信息，以整条reads为拼接单位，整体研究拼接问题，设计合理的评分机制，选择得分最高的kmer进行contigs扩展，直至拼接结束，是基于reads引导的基因组序列拼接算法SRGA的主要思想。

3.1 kmer评分机制

Kmer的评分机制基于kmer的碱基质量及所在reads的拼接信息累计。针对某特定的kmer，其得分计算方法为：首先计算拼接时当前kmers所在的每条reads的得分，并将得分累加，即得该k-er的正向得分；然后取该kmer的反向互补序列，并计算其得分，两个得分相加，即为该特定kmer的综合得分。

将reads按kmer所出现的位置及reads的方向划分成5个分值区域，其中的k-mers按照k-mer的位置和该read的方向划分成5个区域，如图2所示，“+”代表reads的正向得分区域，“-”代表reads的反向得分区域。

图2 reads的得分区域示意

Fig.2 Score display of reads

kmer所在reads的得分计算公式为：

（2）

kmer综合得分的计算公式为：

（3）

序列拼接过程中，一般会优先选择拼接即将完成的reads，但考虑到reads的数据区域特征，正向reads越靠近末端，其碱基的错误率越高。为提高序列拼接的拼接质量，将kmer拼接次数和碱基的错误率结合起来，以评分的方法进行量化，为kmer的导航选择提供依据。

Kmer评分机制的数学模型如图3所示，图中阴影部分的面积即为kmer的得分。该得分综合了reads的累计拼接信息和kmer所在reads区域的碱基质量，reads拼接次数越多，碱基序列所在的区域分值越高，对kmer的得分贡献越大。

图3 kmer得分的数学模型

Fig.3 Mathematical model of kmer score

3.2 contigs的拼接

基于reads引导的contigs拼接过程如下：

step 1：选择De Bruijn图中出现次数为所有kmer出现次数平均值的kmer作为初始kmer，并将该kmer作为初始的contigs，开始第一轮拼接；

step 2：取得contigs扩展待选的所有kmer，根据reads拼接信息及kmer所对应的区域得分，对这些kmer进行评分，选择分值最高的kmer进行拼接；

step 3：如果待选kmer为空，且拼接状态处于第二轮拼接，则停止拼接，转向step 6；如果仍处于第一轮拼接，则拼接状态标记为第二轮拼接，取contigs的反向互补，在其3端选择kmer，进行第二轮拼接；

Step 4：更新reads拼接信息，删除拼接失败或成功的reads信息；

Step 5：删除De Bruijn图中所有拼接成功reads中的kmer信息，并读取reads的碱基序列，添加到contigs；

Step 6：拼接结束，保存contig及拼接成功的reads。

contigs分两轮进行，采用双向扩展。第一轮拼接从contig3端开始，直至扩展终止；第二轮取contig的反向互补，其3端重新进行扩展，直至扩展终止。contigs拼接过程如图4所示。

图4 contig双向扩展

Fig.4 Two directions extension of contigs

4 实验结果及分析

选取大肠杆菌和酵母菌的基因组作为测试数据，对SRGA进行测试，将拼接结果与velvet[8]、Edena[9]和SOAPdenovo[10]等常用软件进行对比，并用Mauve Assembly Metrics软件对各自的拼接结果进行指标分析、客观评测。

选定的大肠杆菌E.coil（ERR001665）基因组数据，其数据量为4.2G ，reads长度为36bp，目标基因组序列为4.6M。

选定的酵母菌 S.pombe（ERR018885）基因组数据，其数据量为1.3G ，reads长度为54bp，目标基因组序列为12.6M。

两组数据的测评结果分别如表1和表2所示。

表1 大肠杆菌的测评结果

Tab.1 Metrics of E.coil

如表1所示，针对大肠杆菌基因组，SRGA的Miscalled bases、N50、Missing bases、DCJ distance等指标上优于SOAPdenovo与velve，N90、N50、DCJ distance等优于Edena。

表2 酵母菌的测评结果

Tab.2 Metrics of S.pombe

如表2所示，针对酵母菌基因组，SRGA的N90、N50、DCJ distance、Missing bases，等指标Edena及velvet，Missing bases、DCJdistance、 Miscalled bases等优于SOAPdenovo。

通过全面分析两组基因组数据的测评结果可知，采用SRGA进行的序列拼接，无论是拼接精度还是拼接的contigs长度均表现不俗，且SRGA适用于大型基因组的序列拼接。

5 结束语

本文首次在基于reads引导的序列拼接算法中，提出了基于基因组序列数据特征和拼接信息累计的评分导航机制，避免了传统的基于reads拼接的计算模型中数据质量信息丢失的问题，为contigs的快速拼接提供了精确导航。然而，基因组序列拼接是一个NP难问题，需要不断地创新和探索。本文的研究工作只是解决该难题的一小步，期待在此基础有更大的进展。

参考文献：

[1] SANGER F， NICKLEN S， COULSON A R. DNA sequencing with chain-terminating inhibitors. 1977. [J]. Biotechnology， 1992， 24（12）：5463-5467.

[2] MAXAM A M， GILBERT W. A new method for sequencing DNA[J]. Proceedings of the National Academy of Sciences， 1977， 74（2）：560-564.

[3] POP M， SALZBERG S L. Bioinformatics challenges of new sequencing technology[J]. Trends Genet， 2008， 24（3）：142-149.

[4] PEVZNER P A， TANG H， WATERMAN M S. An Eulerian path approach to DNA fragment assembly[J]. Proceedings of the National Academy of Sciences of the United States of America， 2001， 98（17）：9748-9753.

[5] ERLICH Y， MITRA P P， Delabastide M， et al. Alta-Cyclic： a self-optimizing base caller for next-generation sequencing[J]. Nat Methods，2008， 5：679-682.

[6] SCHEIBYE-ALSING K， HOFFMANN S， FRANKEL A， et al. Sequence assembly[J]. Comput Biol Chem， 2009， 33（2）：121-136.

[7] DOHM J C， LOTTAZ C， BORODINA T， et al. Sharcgs， a fast and highly accurate short-read assembly algorithm for De Novo genomic sequencing[J]. Genome Research， 2007， 17（11）：1697-1706.

[8] ZERBINO D R， BIRNEY E. Velvet： Algorithms for De Novo short read assembly using De Bruijn graphs[J]. Genome Research， 2008， 18（5）：821-829.

[9] HERNANDEZ D， FRANCOIS P， FARINELLI L， et al. De Novo bacterial genome sequencing： millions of very short reads assembled on a desktop computer[J]. Genome Research 2008， 18（5）：802-809.

[10] LI R， ZHU H， RUAN J， et al. De Novo assembly of human genomes with massively parallel short read sequencing[J]. Genome Research， 2009， 20（2）：265-272.