DNA测序技术及其拼接算法综述

李艳慧，张少强

（天津师范大学计算机与信息工程学院，天津300387）

DNA测序技术用以分析特定DNA片段的碱基序列（腺嘌呤（A）、胸腺嘧啶（T）、胞嘧啶（C）与鸟嘌呤（G））的排列方式.DNA序列的测定是分子生物学研究的一项重要且关键的内容，如基因的定位、基因结构和功能的研究、基因表达与调控、基因与疾病的关系等，都需要对DNA结构有详细的了解.因此，DNA测序技术极大推动了生物学和医学的研究和发展.1990年代启动的人类基因组计划，其主要目的就是测定人类DNA的30亿个碱基对序列及识别人类DNA中所有的基因，这直接推动了生命科学和现代医学的发展，也加快了DNA测序技术的日趋成熟.

自DNA测序技术发明以来，人类对各物种基因和基因组结构的研究和探索从未中断.从第一代测序技术到第四代测序技术，测序的读序（reads，即测序仪器产生DNA片段的序列）长度由长到短，然后再由短变长，测序技术取得了长足的发展.目前，第一代测序技术产生的基于荧光标记和Sanger的双脱氧链终止法原理的测序仪慢慢退出了历史舞台，第二代测序技术在市场上占有绝对的优势，但第三代和第四代测序技术也在迅猛发展.第二代测序技术目前是测序的主流，其显著特征是高通量，一次能对几十万条DNA片段进行序列测序.第三代是以单分子测序为特点的测序技术，这种技术无需PCR扩增，更能反应样本的真实情况，测序通量也更高，操作过程更简单，成本更低.第四代测序技术又称纳米孔测序技术，它具有高通量、低成本，并能快速进行基因检测的特点.测序技术的发展历程如图1所示.

测序完成后的第一步也是最重要的一步就是根据读序拼接回贴成完整的序列，其测序与拼接过程如图2所示.拼接算法用于将测出的读序拼接成完整的染色体序列（chromosomesequence）、转录本序列（transcript sequence）或因为测序不完整而形成的支架序列（scaffolds），最终形成基因组、转录组序列或其他功能序列.所有的拼接算法基本上是先将读序根据重叠关系连接成“叠连群”（contigs），然后再将叠连群根据配对（pair-end）读序和其他信息形成有序集合来构建更长的支架序列.

图1 测序技术的发展历程Fig.1 Development of sequencing technology

图2 DNA测序及拼接过程示意图Fig.2 Diagram of DNA sequencing and assembly

1 第一代测序

1.1 第一代测序技术

第一代测序技术起源于1970年代，测序采用均一的单链DNA分子.第一代测序技术的主要方法有化学降解法、双脱氧链终止法、荧光自动测序技术和杂交测序技术.

化学降解法的优点是准确性较好，易掌握，缺点是过程较繁琐.双脱氧链终止法又称Sanger法，该方法简便、快捷，经过后续的不断改良，成为了第一代DNA测序的主流.Sanger法测序的原理是：每个反应含有所有4种脱氧核苷酸三磷酸（dNYTP），使之扩增，并混入限量的一种不同的双脱氧核苷酸三磷酸（ddNTP）使之终止.由于ddNTP缺乏延伸所需要的3′-OH基团，使延长的寡聚核苷酸选择性地在G、A、T或C处终止，终止点由反应中相应的双脱氧核苷酸而定.荧光自动测序技术基于Sanger原理，使用荧光标记代替同位素标记，利用成像系统自动监测，进而提高了测序的准确性和速度.杂交测序技术利用DNA杂交原理，将一系列已知序列的单链寡核苷酸片段固定在基片上，将待测DNA样本片段变性后与其杂交，根据杂交情况排列出样品的序列信息.

1.2 第一代测序算法

第一代测序技术最早应用的拼接算法是OLC（overlap-layout-consensus）算法，最初由 Staden（1980年）[1]开发，随后许多科学家对其进行了扩展和阐述.OLC算法是一个直观的拼接算法，它适用于读序长度较长的测序数据.该算法一般分为3个步骤：（1）找到所有读序之间的重叠（O）；（2）对所有的读序根据信息在图上进行布局（L）；（3）推断出共有序列（C）.随着Sanger测序技术的广泛应用，基于OLC的拼接算法有Arachne[2]、Celera Assembler[3]、CAP3[4]、PCAP[5]、Phrap[6]、Phusion[7].Arachne拼接算法的关键特征包括：寻找重叠读序的高效和灵敏、纠正拼接前的错误实现高度准确性的重叠评分步骤、基于正向-反向拼接的读序合并、通过正向-反向拼接不一致检测重复叠连群.CAP3是第三代CAP序列拼接算法，它具有修剪5’端和3’端低质量区域读序的能力.CAP3和Phrap比较而言，Phrap通常比CAP3产生更长的叠连群，而通常CAP3在共有序列中产生的误差比Phrap少，并且使用正向-反向约束的低通量数据，CAP3比Phrap更容易构建支架序列.Phusion使用读序数据和配对读序数据提高了大型基因组数据拼接结果的高度精确性.

2 第二代测序

2.1 第二代测序技术

近几年，第二代测序技术在基因组学的研究得到了广泛的应用.第二代测序技术也被称为新一代测序技术（next generation sequencing，NGS）.第二代测序技术相比于第一代测序技术，成本大幅度下降.第二代测序技术主要包括罗氏454公司的454级焦磷酸测序技术、Illumina公司的Solexa测序技术和ABI公司的SOLiD测序技术.目前全球使用量最大的第二代测序技术是Illumina公司的Solexa和Hiseq测序技术，这2种测序技术的核心原理是相同的，都采用边合成边测序的方法，它的测序过程主要分为4个步骤：（1）DNA待测文库构建；（2）Flowcell（即将基因组DNA的随机片段附着到光学透明的玻璃表面）；（3）桥式PCR扩增与变性；（4）测序.

第二代测序技术引入表观遗传学后，形成了各种表观遗传学测序的方法.其中，第二代测序技术与染色质免疫共沉淀技术相结合的ChIP-seq技术，已成为功能基因组学，尤其是基因表达调控领域研究的关键技术.当前，ChIP-seq的应用主要包括：（1）DNA序列上转录因子结合位点的识别；（2）表观遗传学领域中基因组DNA甲基化、组蛋白修饰、核小体定位等.目前，第二代测序平台中的Illumina/Solexa测序平台上的RNA-seq技术应用最广.该技术利用高通量测序对组织或细胞中所有RNA反转录而成的cDNA文库进行测序，统计相关读序的数目并计算不同RNA的表达量，进而发现新的转录本，而不再需要构建克隆文库.以亚硫酸氢盐转换法与第二代测序技术相结合的全基因组亚硫酸氢盐测序法是高通量检测基因组中DNA甲基化的标准方法.DNase-seq和FAIRE-seq是2种在全基因组范围内鉴定染色质可接近性、检测与调节活性相关的DNA序列的方法.FAIRE-seq是在DNase-seq的基础上产生的.将第二代测序技术与染色质构象捕获技术（chromosome conformation capture，即3C技术）相结合，用高通量的方式检测染色体的高级结构及基因组位点间的相互作用，产生了如4C-seq、5C-seq、Hi-C 和 ChiA-PET 等技术.

2.2 第二代测序算法

2.2.1 RNA-seq数据转录组拼接算法

第二代测序RNA-seq数据转录组的拼接算法有2种形式（如图3所示）：基于参考基因组的拼接和从头（de novo）拼接（不基于参考基因组的拼接）.基于参考基因组的拼接首先将大量的RNA-seq测序短读序准确地回贴到基因组的相应位置再进行拼接，而从头拼接是根据RNA-seq测序片段之间的重叠关系进行拼接.

图3 RNA-seq数据拼接的2种形式Fig.3 Two categories of RNA-seq data assembly

第二代测序RNA-seq转录组拼接算法有3种类型：剪接图（splicing graph）、重叠图（overlap graph）和De Bruijn图.基于参考基因组拼接大致包括回贴、构图和根据图结构设计算法重构表达转录本.其中，在回贴部分常使用的联配工具有Hisat[8]、Tophat[9]、Tophat2[10]、SpliceMap[11]、MapSplice[12]等，这些工具可以处理可变剪切.回贴结束后，所有测序片段会由于回贴到基因组相应位置而聚类成簇（代表特定基因相应的表达情况），进而形成有向无圈图结构（剪接图或重叠图）.

基于参考基因组的拼接算法包栝StringTie[13]、Bayesember[14]、Cufflinks[15]、 Scripture[16]、IsoLasso[17]、iReckon[18]、CER[19]、Traph[20]、CIDANE[21]和 Scallop[22]等 .比较典型的算法是Cufflinks和StringTie.Cufflinks的主要思想是根据所有回贴后兼容的测序片段的兼容关系建立重叠图（有向无圈），图中的点代表双端测序片段，点之间的有向边代表点之间的测序片段是兼容的.通过最小路覆盖模型[23]，从图中找到数目最少并且能覆盖所有点的路，将找到的路作为最终预测的表达转录本的集合.StringTie的主要思想是根据所有回贴后兼容的测序片段的兼容关系建立剪接图，图中的点代表外显子，点与点之间的有向边表示2个外显子之间发生剪接事件.与Cufflinks的不同之处在于，StringTie设计了一个贪婪算法逐条预测表达的转录本，然后采用最大流模型逐条估计转录本表达量，最后对剪接图进行更新，迭代进行.Scallop算法能精确拼接识别多外显子和低表达水平的转录本，与StringTie相比可以得到更多正确的多外显子转录本.

转录组的从头拼接算法包括ABySS[24]、SOAP de novo-Trans[25]、Oases[26]、IDBA-Tran[27]、BinPacker[28]、Bridger[29]、Trinity[30]等.其中 ABySS、SOAP de novo-Trans、Oases和IDBA-Tran本质上是基因组拼接的方法.Trinity是应用最广、公认度最高的从头转录组拼接软件，也是第一个专门针对转录组拼接开发的软件.Trinity的主要思想是将测序短序打碎成k-mer（连续k长度的片段），并用哈希表存储每个k-mer及其出现的频率，然后选择频率最高的k-mer构建De Bruijn图，最后通过动态规划算法遍历De Bruijn图，根据打分机制筛选出符合条件的路.Bridger将Cufflinks和Trinity两者优势相结合，它的主要思想是构建k-mer哈希表后，选取频率最高的k-mer作为种子直接延伸出整个连通分支，作为一个剪接图，之后在重复选取k-mer延伸连通分支时，被使用过的k-mer根据一定规则被继续使用，最后对每个剪接图应用最小路覆盖算法找出最小路覆盖作为输出.

2.2.2 ChIP-seq数据相关算法

ChIP-seq数据分析流程见图4，其主要阶段包括读序回贴（read mapping）[31]，峰值检测（peak calling）[32-35]和模体识别[36].现在比较流行的回贴软件包括Bowtie[37]、BWA[38]、GSNAP[39]、MAQ[40]等.其中，Bowtie 和 BWA 是基于Burrows-Wheeler的转换法.峰值检测包括3个步骤：寻找sharp peak和broad peak；生成富集倍数（fold enrichment）轨道；生成logLR轨道.峰值检测软件主要有 MACS[41]、CisGenome[42]、SICER[43]、CCAT[44]、ZINBA[45].其中MACS和CisGenome适用于分析尖锐信号，SICER和CCAT适合分析连绵信号，ZINBA和MACS适合分析混合型峰.模体识别是指在给定的一组峰值区域中识别模体并找出转录因子结合位点.其相关算法有DREME[46]、MEME-ChIP[47]、WEEDER[48]等.DREME 专门设计用于发现短核DNA结合模体的真核生物的转录因子，它可以在几分钟内优化分析ChIP-seq数据集.DREME可作为基于模体（motif）的序列分析工具的MEME[49]组件的一部分.MEME-ChIP用于DNA大数据集模体分析，其网络服务设计用于分析ChIP-seq数据峰值区域.WeederWeb[48]是一个Web接口，用于在一组相关的调控DNA序列中自动发现保守模体.

图4 ChIP-seq数据分析流程Fig.4 Analysis process of ChIP-seq data

2.2.3 第二代测序基因组拼接算法

目前基于第二代测序技术产生的短读序序列进行基因拼接的算法，大致分为OLC图、Greedy图和De Bruijin图3种类型.因为第二代测序技术产生的基因片段较短且高通量测序速度较快，使得De Bruijn图算法在第二代测序技术中占据主要地位.

OLC图算法最早提出于第一代测序，之后科学家不断对其进行阐述和扩展.该算法遵循overlap/layout/consensus框架，其原理是利用重叠图，用短读序表示重叠图的结点，读序之间的重叠关系表示重叠图的边，通过计算短读序重叠关系进行短读序基因拼接.其优点是将序列拼接问题转换为图论问题，适用于长读序基因拼接，缺点是算法无法解决重复序列问题，而且对于第二代基因测序数据的特点，算法运算量较大.典型 算 法有 Arachne[2]、CAP3[4]、EDENA[50]、Newbler[51].其中CAP3和Arachne是在第一代测序时提出的，CAP3是经典的OLC算法，该算法能快速找出读序之间的重叠关系，能有效计算和评价重叠关系.EDENA也是经典的OLC算法，其优点是拼接正确的contig概率更大.

Greedy图算法是最早提出的面向第二代测序技术的拼接算法之一，它遵循贪婪算法的原理，给定一个读序或seed，只要它满足条件就对这个seed不断地延长直到所有的数据用完或超过设定阈值.其优点是原理简单，操作方便，缺点是算法容易得到局部最优解，难以得到全局最优解.基于Greedy图的算法有SSAKE[52]、SHARCGS[53]、VCAKE[54]、QSRA[55].SSAKE 和VCAKE是典型的Greedy算法，选择最大重叠进行组装.SHARCGS和QSRA的相似之处在于算法主要有滤除、组装和融合3个步骤.

De Bruijin图算法遵循Euler超路径算法，其原理是将读序打断成长度为k的核苷酸片段（k-mer），利用核苷酸片段间的overlap关系构建DeBruijn图，然后得到叠连群，通过遍历De Bruijn图得到支架序列，进一步操作后得到基因组序列.其优点是能有效解决重复序列问题，使用哈希表存储方便快速查找，缺点是构造De Bruijin图对内存要求较高，消耗计算资源较大.基于 De Bruijn 图的算法有 ABySS[24]、SOAP de novo[25]、EULER-SR[56]、ALLPATHS[57]、Velvet[58].Velvet和 ABySS是比较经典的基于De Bruijn图的拼接算法，其中ABySS最大的特点是采用De Bruijn图的分布式表示，在多台计算机网络中并行计算进行拼接.SOAP de novo是一个短读序的基因拼接算法，该算法能精确拼接大型基因组数据.

3 第三代测序

3.1 第三代测序技术

第三代测序技术是在第二代测序技术基础上增加读序长度，加快测序运行速度，其特点是单分子测序（不经PCR直接进行边合成边测序）.第三代测序技术包括Heliscope测序技术、SMRT（single molecule real time，单分子实时测序）、离子半导体测序技术（ion torrent）等.其中SMRT测序技术是较为成熟的.SMRT技术以SMRT芯片（带有很多ZMW（zero-modewaveguides）孔的厚度为100 nm的金属片）为测序载体进行测序反应，并采用边合成边测序的思想.第三代测序技术原理是：首先将长读序互相比对并识别长读序之间的重叠区，如一个较短的读序完全包含在一个较长的读序内或者2个长读序之间有相同的序列区域；然后通过不同的扩展策略对长读序进行拼接.其优点是更容易匹配读序序列和降低测序成本，缺点是测序方法不够健壮，错误率高，比对时间较长，不能被广泛使用.

3.2 第三代测序算法

3.2.1 第三代测序拼接算法

第三代测序拼接算法利用长读序跨过基因组中大部分重复区，该类方法一般首先将长读序相互比对并识别长读序之间的重叠区，然后通过不同的扩展策略对长读序进行拼接.

目前已有的长读序拼接工具主要包括HGAP[59]、Dazzler[60]、Sprail[61]、MHAP[62]、Sparc[63]、PBcR[64]等 .HGAP首先建立一个加权的有向无环图来表示比对在一起的多个长读序的关联关系，然后通过找到图中的最大权重来确定一致序列，最后利用Celera[65]拼接纠错后的长读序.自拼接的PBcR方法通过BLASR[66]软件识别长读序之间带有噪声的重叠区，利用重叠区对长读序进行纠错，最后同样利用Celera对纠错后的长读序进行拼接.Dazzler方法基于排序和合并策略，设计了一个高敏感度、高效率的过滤器，用来处理长读序之间局部比对所包含的种子点，并提出了一个线性时间的启发式算法用来找到经过种子点的长读序之间的局部比对.Sprail方法首先选择一定比例的任意长读序作为种子读序，利用比对程序将其他长读序比对到种子读序，从而对种子读序进行纠错，最后使用Celera拼接纠错后的长读序.MHAP方法提出了一个概率算法，可以有效地识别带有噪声的长读序之间的重叠区，具体而言，将长读序重叠区识别简化为整数序列比较问题，从而大幅减少长读序重叠识别时间.

3.2.2 第三代测序纠错算法

第三代测序技术解决高错误率策略主要有2个：（1）利用来自另一个测序平台的低错误率的短读序，在保证足够的覆盖倍数下，纠正长读序中的错误，该类方法称为混合方法；（2）利用长读序的错误分布比较均匀的事实，在保证长读序足够的覆盖率下，利用长读序自身的信息进行纠错，该类方法称为长读序自纠错方法.

目前，针对PacBio长读序的纠错算法有：（1）自纠错方法：DAGCon[67]、PBcR[64]、LoRMA[68]；（2）基于短读序的纠错方法：PBcR[64]、LSC[69]、PacBioToCA[70]、proovread[71]、ECTools[72]、Cerulean[73]、LoRDEC[74]、Jabba[75].DAGCon 的主要思想是通过长读序间的多对多比对进行自身纠错.PBcR的主要思想是通过长读序间的局部比对构建多重序列比对，然后生成共有序列来进行纠错.PBcR的最新版本既支持长读序的自纠错，也支持基于短读序的纠错.LSC的主要思想是先将长读序中相同的碱基压缩，然后再将短读序比对到长读序，利用比对结果对长读序进行纠错.PacBioToCA直接将短读序比对到长读序上，利用比对结果对长读序进行纠错.

4 第四代测序技术

第四代测序技术又称为纳米孔测序技术，它不需要对DNA进行生物或化学处理，而是基于物理的原理进行测序.其简要原理是：将经过提取的DNA的一端与一种特殊的蛋白连接，该蛋白主要起到与测序纳米孔的连接及控制序列进入速度的作用，同时将DNA的另一端连接，这样DNA的正义链和反义链可以依次通过纳米孔.通过纳米孔的DNA单链可以产生微弱的电流变化，感应器感受电信号并导出，通过对随时间变化的电信号进行解析，可以识别基因中碱基对的排列顺序.其优点是超长读序、高通量、低成本、更少的测序时间和更为简单的数据分析，缺点是测序准确率较低，纳米技术在制造、测序和集成等方面也存在着诸多挑战.

目前，针对Nanopore长读序的纠错方法包括Nanocorr[76]、NaS[77]、MiRCA[78]混合纠错算法.Nanocorr的主要思想是将短读序比对到长读序上，通过最大递增子序列筛选出最优的短读序比对片段，过滤掉噪音信息.NaS的主要思想是将短读序比对到长读序上，将完全匹配的子序列作为种子，通过种子招募相似的短读序，最后对招募到的短读序进行拼接.MiRCA的主要思想是利用拼接工具将短读序拼接成长的叠连群，然后将叠连群比对到长读序上，利用叠连群和长读序的比对结果，对长读序进行纠错.

5 总结与讨论

随着基因测序技术的不断发展，测序拼接的正确率和速度也在不断提高.第一代测序技术的测序读长可达1 000 bp～2 000 bp，准确率高达99.999%，但存在测序成本高，通量低等方面的缺点.随着基因组研究的不断深入，第一代测序技术已经无法满足人们大规模测序工作的需求，这直接推动了第二代测序技术的产生.由于第二代测序技术产生的读序较短，而且基因组中存在大量复杂的重复区，这些重复区使基于第二代测序技术的短读序序列组装结果比较零散，无法获得完整的基因组序列.第三代测序技术所产生的长读序长度可以达到数万碱基，这些长读序可以跨过基因组中大部分的重复区，进而帮助研究人员获得完整的基因组序列.但是第三代测序技术的测序错误率高达15%左右，这对长读序组装结果的质量和效率造成了较大的影响，特别是错误富集的区域.第四代测序的纳米孔技术，其显著优点为无标签、超长读序（104～106个碱基）、高通量和低材料需求.这些优点使实验过程大大简化，可以很容易地用于DNA测序应用.四代测序技术的比较见表1.

目前第四代技术尚未大面积推广，最新的拼接算法研究思路是将第二三代测序的短读序和长读序相结合，这类算法主要有LINK[79]、SSPACE-LongRead[80]等，另外很多拼接算法存在占用内存资源较多的缺陷.

目前拼接算法面临的主要挑战有：如何有效地分析大型数据集、处理序列错误和正确捕获基因组的重复结构.除此之外，还需要继续开发新的方法来表示和分析拼接图，从而正确地拼接重复序列并正确区分共组装单倍体（co-assembled haplotypes）的基因组变异体.同时，为有效地处理越来越大的数据集，需要对基于De Bruijn图的拼接算法进行并行化.此外，内存效率也是基因组拼接算法的重点研究内容.为减少内存，相关研究在拼接算法中引入了稀疏图、数据结构压缩、布隆过滤器和重叠有效计算的FM索引等.其中，基于布隆过滤器的算法利用了布隆过滤器的空间效率很高的随机数据结构，布隆过滤器利用位数组简洁地表示集合，并判断元素是否属于集合，这使得基于布隆过滤器算法的空间效率和查询时间都远远优于一般的算法，其缺点是有一定的误识别率和删除困难.这些新的高效内存拼接算法允许分析更大的数据集，并且不需要高性能的计算设备.总之，目前拼接算法的应用大部分都处于串行的实现阶段，拼接速度不够高，在序列拼接问题中如何实现并行处理，提高拼接准确率以及减少内存资源占用等方面的研究具有重要意义.

现代测序实验通常会产生配对数据，这些数据可以限制读序对之间的相对方向和距离的信息.分析这些信息可以解决序列重复问题，并可以将单个的重叠群连接成长度较长的支架，这是提升拼接算法效率的一个重要领域.同时，在低覆盖率和重复区域进行局部拼接的算法设计是最困难和亟待解决的问题.

目前，在对拼接问题进行广泛的数学分析的情况下，依然需要依靠启发式算法和其他技术，这是因为难以提出实际的拼接数学模型.应用于拼接的大多数数学公式表明，虽然找到最佳的解决方案是可能的，但仍需要大量工作.

表1 测序技术的比较Tab.1 Comparison of sequencing techniques