基于混合策略的de,novo序列拼接算法构造
时间：2023-02-15 14:35:08 来源：千叶帆 本文已影响人
肖存威,石海鹤,王 岚,程柏良
(江西师范大学计算机信息工程学院,江西 南昌 330022)
以人类基因组计划(human genome project,HGP)的实施为界,可以将生物信息学(Bioinformatics)的发展大致分为2个阶段,分别为前基因组时代和后基因组时代.HGP的主要贡献是提供了一项新的技术——基因测序[1].测序技术是进行生物研究最基本的技术之一,近几年来,随着测序技术的高速发展,这既降低了测序成本,又推动了生物信息学的发展,对精准医疗和基因检测等研究产生了促进作用.
在测序技术的发展过程中,一共产生了3代基因测序技术.第1代基因测序技术在分子生物学研究中发挥了重要作用,比较知名的基因测序方法是双脱氧法[2].第2代基因测序技术是在HGP的研究过程中被逐步开发出来,主要有Roche公司的454[3]、Illumina公司的Solexa[4]和ABI公司的Solid[5]等技术方法.第2代基因测序技术是通过对DNA片段复制扩增,把碱基的信号强度放大,加入dNTP,在每个接头上聚合新的碱基,就可以得到在整个接头上的DNA序列信息[6].目前出现了以单分子测序为特点的第3代基因测序技术,使得在测序成本和测序通量上有了较大的改善.第3代基因测序技术主要以生物科学公司的HeliScope[7]、牛津纳米孔技术公司的Nanopore[8]以及太平洋生物科学公司的PacBio[9]等技术方法为主流[10].与第2代基因测序技术不同的是,第3代基因测序技术既无需打断整条DNA链,又无需对PCR扩增,就可实现对DNA分子的单独测序[11].测序读取的最大长度可达到12～15 kbp,从而减少了拼接次数,降低了成本.不进行PCR扩增,避免了在扩增过程中出现错误,然而相比于测序的准确率来说,第3代基因测序技术的平均错误率约为15%,这远远高于第2代基因测序技术的平均错误率(1%).
因为当前测序读取的最大长度的序列也无法将整个基因组的遗传信息表达准确,所以需要将经测序得到的短基因序列通过相关的技术拼接起来,得到一条完整的基因序列,这就是序列拼接(sequences assembly),其流程如图1所示.
图1 序列拼接流程图
按照是否有参考序列参与的标准,序列拼接可分为2类:第1类是de novo(从头)拼接,不需要参考序列,仅依靠短序列之间的重叠关系拼接成长序列;第2类是参考序列拼接,需要本物种或相似物种的基因组序列作为参考序列,以降低拼接过程的复杂度和成本.在实际应用中对一个新物种进行的测序往往没有合适的参考序列,在这种情况下只能采取de novo方式进行序列拼接.
随着物种的范围越来越大,序列拼接的准确度和长度要求越来越高,这给序列拼接工作带来挑战.在拼接序列检测方面,目前拼接得到的大多数序列是基于de novo拼接的,没有参考基因组,无法从基因组层面对结果进行比对检测,即使有参考基因组参与比对验证,在拼接过程中碱基也可能会发生突变,从而降低拼接的准确度.在测序方面,目前的测序技术无法同时满足低错误率和读长较长的要求,如第2代基因测序准确度较高,但是读长较短、测序时间长、费用高,而第3代基因测序读长较长、时间短、测序成本低,但准确度较低.在序列拼接算法方面,大多算法是基于单一策略开发的,具有一定的缺陷,如贪婪策略算法是利用启发式搜索方式[12],选择与种子reads重叠最多碱基的reads进行拼接,最后得到局部最优解或次优解,却无法保证得到全局最优解.基于OLC(overlap-layout-consensus)策略和DBG(De Bruijngraph,De Bruijn图)策略开发的算法,在单独使用时存在缺陷.对于复杂度较高的基因组,若单独采用OLC策略,则杂合reads之间的同源性降低,容易导致重叠;若单独采用DBG策略,则所形成的有向图会存在较多的分叉,导致contigs长度减小,无法进一步拼接.因此,可以尝试研究在混合策略下的序列拼接问题,期望得到更好的序列拼接结果.
本文在分析同一拼接策略算法实现原理的基础上,进一步研究在相同策略下的不同算法在具体实现过程中的差异性,得到算法在领域内的共同特征和差异特征.利用形式化方法及形式化平台方面的优势,结合领域分析建模和产生式编程的方法,形式化开发算法在领域内的序列拼接算法构件.对于共同特征可以生成单种策略下的必须执行构件,差异特征可以生成在单种策略下的可选执行构件[13].根据用户对可选构件的选择和输入文件的特征自动组装序列拼接算法.根据(OLC+DBG)→OLC混合模式将不同策略算法在领域内生成的拼接算法进行组装,构造得到基于混合策略的ODO拼接算法.ODO算法在综合不同策略优势、提高拼接效率和准确度的同时也保证了结果的可靠性.
1.1 单一策略拼接
基于贪婪策略开发的序列拼接算法是在整个序列拼接研究过程中最早提出的一类算法,也被称为贪心算法,是从问题的某一个初始值出发,逐步从两边逼近目标值,以尽快地求得一个较好的解,当算法运行到某一步不能继续前进时,算法停止.运行结束后会得到一个局部最优解或次优解,在运行过程中这一类算法所做的每个选择都是在当前状态下最好的选择(贪心选择),贪心策略拼接原理示意图如图2(a)所示.常见的基于贪婪策略开发的序列拼接算法有SSAKE[14]、VCAKE[15]、SHARCGS[16]、NPSCARF[17]等.
OLC策略主要基于reads之间的重叠关系进行拼接,是最能直观体现序列拼接原理思想的策略.OLC策略在序列拼接时将所有待拼接的reads构造成一个有向图,图中的每个节点都代表特定的reads;若2个节点之间存在边,则说明这2个reads的重叠部分的碱基数量大于设置的阈值.对图进行简化和分割后,可以在一系列子图上寻找经过每个节点一次且仅一次的路径,最后得到所需的序列.由此把序列拼接问题转化为Hamilton路径问题,OLC策略拼接原理示意图如图2(b)所示.基于OLC策略开发的序列拼接算法有CABOG[18]、Phrap[19]、Celera Assembly[20]、Edena[21]等.
DBG策略的特点是不需要进行序列比对,可大幅度减少系统内存的消耗,De Bruijn图的大小只与基因组大小和算法复杂度有关.在DBG算法中,把reads按照长度k分割成k-mer,并将这些k-mer存入Hash表中,不同k-mer在Hash表中只存储1次,然后将这些k-mer作为点,构建De Bruijn图.对De Bruijn图进行简化和分割,并在图中找到1条欧拉路径构成contig,再通过Pair-end文库可以确定contig之间的位置关系,指导contig形成长度更长的scaffold,DBG策略拼接原理示意图如图2(c)所示.基于DBG策略开发的序列拼接算法有SOAPdenovo[22]、Velvet[23]、ABySS[24]、EULER[25]、ALLPATHS[26]等.
图2 单一策略拼接原理示意图
1.2 多策略混合拼接
贪婪策略在拼接序列时总是做出在当前情况下最好的选择,不从整体加以考虑,且并非对所有问题都适用.OLC策略在使用时因PCR扩增而会存在大量的重复数据,并且对长reads拼接有偏向性.DBG策略在使用时将短reads重新分割成长度更短的k-mer进行序列拼接,增加了计算时间和拼接的复杂度.由此可见这些策略在单独使用时均存在不足.
为了解决单一策略在序列拼接时易出现的问题,综合各个策略的优点,这可以使用混合模式来实现策略的混合.如(贪婪策略+DBG)→OLC、贪婪策略+OLC、贪婪策略+DBG、DBG+OLC、(OLC+DBG)→DBG、(OLC+DBG)→OLC、(OLC+OLC)→DBG以及(DBG+DBG)→OLC等模式.A.K. Kusuma等[27]对(OLC+DBG)→OLC进行了研究,在实验过程中,首先使用Velvet算法拼接reads,然后使用Edena算法对相同的reads再次拼接,最后使用Minimus assembler对拼接产生的contigs再次拼接可得到最终结果.使用上述方法,A.K. Kusuma等[27]在3个真实数据集和4个模拟数据集上进行了实验.该方法既避免了DBG策略拼接导致的错误率过高的现象发生,又解决了OLC策略拼接带来的重复序列和占用内存过大的问题.
本文充分利用在软件开发方法上及其在平台方面具有的优势,对(OLC+DBG)→OLC混合模式的算法进行了研究.使用形式化、产生式编程等方法和技术,分别构造了基于OLC策略和基于DBG策略的算法域构件库,借助在PAR平台上的Apla→Java程序生成系统生成了基于OLC的拼接算法OLC_assembly_1、OLC_assembly_2和基于DBG的拼接算法DBG_assembly,进一步拼装出混合模式的ODO算法,从而将尽可能多的创造性劳动转化为非创造性劳动,提高了算法开发过程的可理解性、算法的可靠性以及算法构件库在相关领域中的可重用性,ODO算法流程图如图3所示.
图3 ODO算法流程图
2 算法构造
实现相同策略功能的算法可以是不同的,这是因为可按照使用需求生成同种单一策略下的不同算法.软件形式化是实现软件自动化的关键,PAR是一种统一的软件形式化方法,包含了自定义泛型算法设计语言Radl、泛型抽象程序设计语言Apla、统一的算法程序设计方法以及PAR平台(Radl到Apla程序生成系统,Apla到Java、C++等系列程序生成系统)和其他关键技术[28].根据序列拼接算法的需求设计出相应的程序规约,并在此基础上进行程序设计,得到最终的算法和程序.
首先在领域层次上建立领域模型,然后利用产生式编程对算法构件进行交互设计,再形式化实现算法的构件库,最后在不同策略的构件库中自动生成不同单一策略下的拼接算法.利用生成的序列拼接算法按照(OLC+DBG)→OLC模式组装后可得ODO算法,同时利用形式化方法保证了算法构造的高效性和可靠性.下面叙述了基于DBG策略拼接算法的构造过程,阐述了基于OLC策略拼接算法的构造过程.
2.1 DBG策略拼接算法生成
基于DBG策略开发的算法有SOAPdenovo、Velvet、ABySS、EULER等,将这些序列拼接算法联合起来,构成了基于DBG策略的算法拼接领域.通过对在领域内算法的分析,找出算法的共同特征和可变特征并进行分类,选择和定义需要分析和解决的领域,收集相关的领域信息,并形成领域模型[29].根据领域分析的结果,抽取算法共性,用参数表示算法的差异性,设计出Apla语言描述的拼接算法构件,为实现完整的算法构件库,进一步分析不同构件之间的交互模型和构件之间的约束关系.
2.2 OLC策略拼接算法生成
2.2.1 领域特征模型分析 将基于OLC策略的序列拼接算法联合起来,构成OLC算法领域,在领域内的序列拼接算法在进行序列拼接时的过程大致可以分为3步:首先,通过序列比对找到reads之间重叠长度超过设置阈值的重叠关系;其次,通过这些重叠信息将reads建立一种新的组合关系,可得到重叠群contigs,再将contigs进一步排列,生成较长的scaffold重叠群;最后,在经过简化和分割后的有向图中找到一条最优的路径(Hamilton路径)所对应的序列,这即是需要的拼接序列.接下来对在OLC算法领域中的算法CABOG、Celera Assembly和Edena进行分析,找出这3个算法的共同特性和可变特性.如在实现序列核查、序列比对、构建有向图等功能中,这3种算法使用的方法和步骤是类似的.但在纠正测序错误这个步骤中,Edena使用快速检测虚假reads来实现对测序错误的纠正,CABOG使用rocks和stones技术来实现对测序错误的纠正,Celera Assembly使用PacBioToCA自纠算法来实现对测序错误的纠正.在对这3个算法进行领域分析后,找到了3个算法的共同特性和不同特性,归纳得到了OLC领域算法流程图和领域特征模型(见图4和图5).
图4 OLC领域算法流程图
图5 OLC领域算法特征模型
2.2.2 领域构件设计 根据对OLC领域的分析和建模,可抽取出其中的特征,用参数表示其中的差异性,设计出可使用Apla语言描述的算法构件[29].在构建得到ODO算法后,测序所得的序列作为输入,算法将会对输入的序列做核查,避免由硬件问题而出现字符错误的现象发生,在纠正测序错误后,再进行序列比对,构建有向图,去除错误分支,生成contigs和scaffolds.因此,这种OLC算法的主要构件是序列核查构件(Seq_check)、纠错构件(Error_correction)、序列比对构件(Sequence_alignment)、有向图生成构件(Graph_structure)、图化简去分支构件(Graph_simplification)、Contigs拼接构件(Contigs_assembly)、Scaffolds拼接构件(Scaffolds_assembly)和最长路径构件(Longest_path).算法构件关系图如图6所示,这里使用Radl语言为Seq_check构件的功能做形式化描述.
图6 算法构件交互模型
Seq_check构件:此构件对输入的基因序列进行核查,如在DNA基因序列中仅含有4种碱基(A,G,C,T),若序列中含有其他碱基,则会停止运行并提示错误.
|[in A_seq:list of char;B_seq:list of char;n:integer;out_flag:boolean]|;
AQ:n=2∧|A_seq|>0∧|B_seq|=4∧perm(
B_seq,{A,G,C,T});
(perm表示标准核查序列B_seq和序列{ A,G,C,T}互为置换);
AR:flag=(∀i:0≤i≤|A_seq|:(∃j:0≤j≤3:A_seq [i]=B_seq [j])).
在上述的形式化规约中,输入(in)和输出(out)是在PAR平台中定义的2个关键字标识符,序列类型(list)、整型(integer)和布尔型(boolean)是在PAR平台中预定义的数据类型,AQ和AR分别是构件的前置条件和后置条件.
2.2.3 Apla形式化实现 抽象程序设计语言Apla可以直接使用抽象数据类型和抽象过程编写程序,可抽象描述算法,并进行正确性验证,这就保证了生成算法的可靠性,并且易于生成各种可执行程序设计语言程序Java、C++等.在Apla中的标识符、关键字、符号表达方式和在Radl语言中涉及的一致.
在使用形式化方法开发图6中算法构件的Apla实现后,利用PAR平台Apla→Java程序生成系统得到对应的Java代码.
下面仅列出Seq_check构件的Apla实现.
function Seq_check(A_seq:list(char);B_seq:list(char,4)):boolean;
vari,j,sum:integer;
begin
i,j,sum:=0,0,0;
do (0 do (j<3)→if(A_seq[i]==B_seq[j])→ sum:=sum+1; fi; j:=j+1; od; i:=i+1; od; if(sum==|A_seq|) →Seq_check:=true []→Seq_check:=false; fi; end. 再得到相应的Java程序如下: public boolean Seq_check(String A_seq,String B_seq){ inti=0; intj=0; int sum=0; for (i=0;i for(j=0;j<3;j++){ if(A_seq[i]==B_seq[j]){ sum++; } } } if(sum == A_seq.length ()){ return true; }else{ return false; } }. 2.2.4 OLC策略拼接算法组装 在构件设计完成后,根据使用需求和输入文件的特征可选择相对应的算法构件组装ODO算法.本文选择Seq_check、Error _correction(rocks and stones)、Sequence _alignment、Contigs_assembly等构件构造出了OLC_assembly_1算法.按照相同的方法,从OLC和DBG领域算法构件库中选择符合使用需求的算法构件,生成OLC_assembly_2和DBG_ assembly算法,进一步组装得到ODO算法,由此将尽可能多的创造性劳动转化为非创造性劳动. 本次实验选择了3种单一策略(贪婪策略、OLC策略和DBG策略)的算法(SSAKE、Edena和Velvet)与ODO构造算法进行对比实验. 本文从GenBank中下载了3种基因组较小的样本来进行实验,实验样本信息如表1所示,基因名称分别是Clostridiumperfringens(产气荚膜梭菌)、Acetobacteraceti(醋酸杆菌)、Escherichiacoli(大肠杆菌),该实验在如表2所示的实验平台上运行. 表1 实验样本信息 表2 实验运行平台 coverage depth增加表示基因组复制的次数增加,经过鸟枪法(shot gun)剪切,可得到更多不同的reads,建立更完善的Hash表和更完整的有向图.但若无限的增加coverage depth的值,则过多的reads会导致文库过大,消耗大量的内存和计算时间,同时还会增加在有向图中的分支和错误节点,降低拼接准确度.本次实验将coverage depth的值设为50. Velvet算法和ODO算法在构建De Bruijn图时会将文库中的reads切割成长度为k的k-mer.k值越大越容易得到唯一的序列,但在拼接过程中越容易产生gaps;随着k值减小,可以增加De Bruijn图的连通性,但会提高算法的复杂度,增加处理时间.同时,为了避免正反链混淆产生回文序列的现象发生,k值应选奇数.本次实验所使用的k值取为31. 构建了OLC和DBG领域构件库,生成了基于OLC的拼接算法OLC_assembly_1、OLC_assembly_2和基于DBG的拼接算法DBG_ assembly,进一步组装出(OLC+DBG)→OLC混合模式算法,简称ODO算法.为评估实验效果,在表2的实验平台上用SSAKE(贪婪策略)、Edena(OLC策略)、Velvet(DBG策略)和构造生成的ODO算法分别对在表1中的3种基因样本(Clostridiumperfringens、Acetobacteraceti、Escherichiacoli)进行序列拼接实验. 由实验结果(见表3)可知:构造的ODO算法无论是在N50大小上还是在基因组覆盖度(coverage depth)上都优于其他的单一策略.如对于Clostridiumperfringens这一样本基因,ODO算法拼接出来的结果N50为1 625 bp,基因覆盖度为93.9%,优于单一策略进行序列拼接的结果. 表3 4种策略算法在不同样本上的拼接结果 在评估指标中N50表示在拼接得到的Contigs长度从大到小排列后最中间Contigs长度的值,N50的值越大说明拼接得到的Contigs长度越长,即拼接效果越好.在对Escherichiacoli实验样本进行拼接时,SSAKE、Edena、Velvet和ODO得到的N50值分别为328 bp、545 bp、673 bp和894 bp.由实验结果可知:相比于其他的单一策略下的算法,ODO算法取得了较优的结果.产生这一结果的主要原因有2个方面:(i)先对输入的reads文库进行2次拼接,然后将2次拼接的结果混合再次作为第3阶段拼接算法的输入,得到最终的Contigs输出;(ii)在ODO算法的构造过程中使用形式化方法开发算法领域构件,可根据用户的需求和输入文件的特征选择不用的构件组装得到策略拼接算法,在提高算法准确度和效率的同时,保证了结果的可靠性. 在评估指标中Contigs number表示拼接结果的Contigs数量,在基因组碱基数量大致相同的情况下,Contigs数量越少表示每条Contigs的长度越长,即拼接结果越好.在对Acetobacteraceti实验样本进行拼接时,SSAKE、Edena、Velvet和ODO得到的Contigs number值分别为106、87、54和67.由评估指标可知:ODO算法的结果比SSAKE和Edena算法的结果更优,但是离Velvet算法的结果还有一定的差距.产生这一结果的可能原因是:在序列中还有较多的gaps,当碱基比对遇到gap时算法会默认当前匹配不成功,并停止匹配,输出结果,这将导致Contigs number反向增大. 本文利用形式化方法及形式化平台方面的优势,结合领域分析建模和产生式编程的方法,形式化开发了序列拼接算法构件,并构造了(OLC+DBG)→OLC混合模式算法,简称ODO算法.实验结果显示:在对3种基因样本进行序列拼接对比实验时,构造生成的ODO算法比单一策略算法所取得的结果在N50和Coverage等参数上有一定的优势.ODO算法在Escherichiacoli样本的序列拼接实验中基因组覆盖度达到97.1%,N50的值达到894 bp.研究结果表明:利用产生式编程和形式化方法开发的ODO算法在综合不同策略优势、提高拼接效率和准确度的同时也保证了结果的可靠性,为深入研究算法构造在序列拼接上的影响提供了参考.3.1 实验材料
3.2 coverage depth和k值对序列拼接的影响
3.3 实验结果与分析
