基因组序列拼接.docx
- 文档编号:24482769
- 上传时间:2023-05-28
- 格式:DOCX
- 页数:21
- 大小:1.23MB
基因组序列拼接.docx
《基因组序列拼接.docx》由会员分享,可在线阅读,更多相关《基因组序列拼接.docx(21页珍藏版)》请在冰豆网上搜索。
基因组序列拼接
2014年成都理工大学
校内数学建模竞赛论文
题目
编号
C题-基因组组装
队编号
9
参
赛
队
员
姓名
学号
专业
张萌立
201313030206
计科
何理
201305090108
空间
张玲玲
201308050710
财务管理
二0一四年五月二十五日
基因组组装
摘要:
本文所要研究的就是全基因组的从头测序的组装问题。
首先,本文简要介绍了测序技术及测序策略,认真分析了基因系列拼装所面临的主要挑战,比如reads数据海量、可能出现的个别碱基对识别错误、基因组中存在重复片段等复杂情况,探讨了当前基因组序列拼接所采用的主要策略,即OLC(Overlap/Layout/Consensus)方法、deBruijn图方法,且深入探讨了deBruijn图方法。
其次,针对题中问题,以一条reads为基本单位,分为reads拼接和contig组装两个阶段,其中contig是由reads拼接生成的长序列片段。
Reads的拼接阶段主要包括数据预处理、de-Bruijn图、contig构建等,而contig的组装阶段主要包括序列的相对位置的确定以及重叠部分overlap的检测,用序列比对的方法来提高拼接的精度。
最后,进行了算法的验证与性能的评价,并且针对问题2,进行了组装分析与验证,结果表明,得到的拼接基因组序列在小范围内与原基因组序列大致吻合。
关键词:
基因组系列拼接;reads;deBruijn图;contig组装;k-mer片段;
一.问题重述
基因组组装
快速和准确地获取生物体的遗传信息对于生命科学研究具有重要的意义。
对每个生物体来说,基因组包含了整个生物体的遗传信息,这些信息通常由组成基因组的DNA或RNA分子中碱基对的排列顺序所决定。
获得目标生物基因组的序列信息,进而比较全面地揭示基因组的复杂性和多样性,成为生命科学领域的重要研究内容。
确定基因组碱基对序列的过程称为测序(sequencing)。
测序技术始于20世纪70年代,伴随着人类基因组计划的实施而突飞猛进。
从第一代到现在普遍应用的第二代,以及近年来正在兴起的第三代,测序技术正向着高通量、低成本的方向发展。
尽管如此,目前能直接读取的碱基对序列长度远小于基因组序列长度,因此需要利用一定的方法将测序得到的短片段序列组装成更长的序列。
通常的做法是,将基因组复制若干份,无规律地分断成短片段后进行测序,然后寻找测得的不同短片段序列之间的重合部分,并利用这些信息进行组装。
例如,若有两个短片段序列分别为
ATACCTTGCTAGCGT
GCTAGCGTAGGTCTGA
则有可能基因组序列中包含有ATACCTTGCTAGCGTAGGTCTGA这一段。
当然,由于技术的限制和实际情况的复杂性,最终组装得到的序列与真实基因组序列之间仍可能存在差异,甚至只能得到若干条无法进一步连接起来的序列。
对组装效果的评价主要依据组装序列的连续性、完整性和准确性。
连续性要求组装得到的(多条)序列长度尽可能长;完整性要求组装序列的总长度占基因组序列长度的比例尽可能大;准确性要求组装序列与真实序列尽可能符合。
利用现有的测序技术,可按一定的测序策略获得长度约为50–100个碱基对的序列,称为读长(reads)。
基因组复制份数约为50–100。
基因组组装软件可根据得到的所有读长组装成基因组,这些软件的核心是某个组装算法。
常用的组装算法主要基于OLC(Overlap/Layout/Consensus)方法、贪婪图方法、deBruijn图方法等。
一个好的算法应具备组装效果好、时间短、内存小等特点。
新一代测序技术在高通量、低成本的同时也带来了错误率略有增加、读长较短等缺点,现有算法的性能还有较大的改善空间。
问题一:
试建立数学模型,设计算法并编制程序,将读长序列组装成基因组。
你的算法和程序应能较好地解决测序中可能出现的个别碱基对识别错误、基因组中存在重复片段等复杂情况。
问题二:
现有一个全长约为120,000个碱基对的细菌人工染色体(BAC),采用Hiseq2000测序仪进行测序,测序策略以及数据格式的简要说明见附录一和附录二,测得的读长数据见附录三,测序深度(sequencingdepth)约为70×,即基因组每个位置平均被测到约70次。
试利用你的算法和程序进行组装,并使之具有良好的组装效果。
附录一:
测序策略
测序策略如下图所示。
DNA分子由两条单链组成,在图中表现为两条平行直线,两条直线上相对位置的两个碱基相互结合形成碱基对(bp),并且与碱基A结合的碱基必为T,与碱基C结合的碱基必为G。
将一个含120,000个bp的完整基因组,随机打断成500bp的片段,然后对500bp的片段进行测序。
测序方法如第3步所示,分别从500bp片段的两端,对两条单链进行测序,测得的读长记为reads1,reads2。
reads1,reads2的长度均为88bp,且该对reads相距500bp。
●图1测序策略示意图
附录二:
数据格式
读长数据格式为fastq格式:
每4行表示一条reads
第一行:
@序列ID,包含index序列及read1或read2标志;
第二行:
碱基序列,大写“ACGTN”;
第三行:
“+”,省略了序列ID;
第四行:
质量值序列:
字符的ASCII码值-64=质量值。
附录三:
读长数据
测序得到的读长数据存放于两个fastq文件中(见附件一),其中McMc_BAC_1.fq.gz.clean.dup.clean和McMc_BAC_2.fq.gz.clean.dup.clean分别存放reads1和reads2的数据。
二.问题分析
正如上面问题所描述的一样,我们要解决的是要将基因小序列read组装成连续的基因大序列乃至最终的完整基因序列,而这就要将两个read1和read2片段进行比较与拼接,比较的时候,因为相似片段的长短问题而不能确定拼接正确性,因此可以用两片段相似的权值来判断拼接的合理性,这样,若用点来代替read,用加权的边来判断到底要和哪个片段进行拼接,我们在查阅资料后,发现可以通过debruijn图并对其进行相应的改进后来建立数学模型对问题进行求解。
设想一本杂志被复制成多份,将每份杂志均以不同的方式剪切,将多份剪切的杂志放在一起。
在剪切的过程中,一些碎片丢失,一些碎片被污渍浸染,一些碎片存在着重叠现象。
根据上述情况来寻找恢复原始杂志的方法。
这是DNA序列拼接问题的现实模型描述。
基于deBmijn图的序列拼接原理主要是通过构造并简化deBmijn图结构来实现整个序列拼接的过程。
三.基于DeBruijn图的序列拼接技术分析与比较
二十世纪八十年代末,Pevzner等人提出基于debruijn图的算法,并首次将该算法用于DNA序列拼接。
基于debruijn图的算法的核心思是将序列拼接问题转换为人们所熟悉的欧拉路径问题。
Pevzner等人认为传统的overlap-layout-consensus算法导致了将DNA序列拼接问题转换为Hamilton路径问题,他们受到杂交测序方法SBH(SequencingbyHybridization)的启发,创造性地提出了在deBruijn图中寻找欧拉路径的构想,尽管杂交测序方法SBH从未在测序工程中实际应用过,但它直接引发了基因芯片工业的诞生。
构造deBruijn图的方法如下所述:
(1)在read集合R={r1,r2,…,rn}中,首先将每一条read分割成若干k-mer(长度更短的DNA片段),分割方法如图1-1所示。
假定集合R中任意一条read的长度均为l,k-mer长度值设为k,那么集合R中的任意一条read均可被分为l–k+1条k-mer,并且这些k-mer作为deBruijn图的顶点。
(2)对于给定的两条k-merx和y,如果在某readri中存在一条长度为k+1的子串,且该子串的前k个碱基与k-merx(或y)精确匹配,同时该子串的后k个碱基与k-mery(或x)精确匹配,那么该算法认为两条k-merx和y之间存在一条公共边。
将采用上述方法构造的deBruijn图记作G。
对于read集合R={r1,r2,…,rn}中的任意一条readri,若在deBruijn图G中存在一条路径P,且该路径P访问ri中的每一条k-mer仅一次,则欧拉路径问题便可理解为:
给定某一deBruijn图G以及G中的路径集合P,在deBruijn图G中确定某一条欧拉路径Q,使得路径集合P中的每一个元素都是欧拉路径Q的子路径。
利用欧拉路径算法进行DNA序列拼接的主要步骤如下所述:
首先利用纠错软件修正read中测序错误的碱基;然后按照上述方法构建deBruijn图;构建deBruijn图之后,应将read集合中的所有read排列在deBruijn图中,在deBruijn图中,每一条read均被视作一条路径;最后在deBruijn图中寻找一条欧拉路径,使得该路径包含deBruijn图中所有read所对应的路径。
在OLC中,在Overlap步骤中,采用了序列比对算法来寻找read之间的重叠信息,该算法的时间复杂度为0(?
2),其中,《SDNA序列中read的数量。
当前DNA测序数据序列越来越短,对同一个物种进行测序,其产生的read数量大大增加,这使得OLC的计算量增加;而基于deBruijn图原理的序列拼接中,抛弃了OLC中序列比对算法,而是采用以k-mer为图中顶点构建图,从而减少了序列比对算法所消耗的时间,提高了算法的效率与overlap-layout-consensus算法相比,基于debruijn图的算法有更低的时间复杂度,这是因为欧拉路径问题实际上是一个线性时间的问题。
利用欧拉路径思想的拼接算法有EULER-SR、ALLPATHS、Velvet和EULER等。
四.模型建立
4.1.1模型的假设
1.假设模型中的read片段都是由一条完整的DNA经过测序而来,它们进过拼接后可以形成一个大片段。
2.模型中出现的各个序列中DNA的双链都准确
3.模型中read在拼接时合理地去掉的公共部分在误差允许的范围内。
4.由于总会在测序中出现read的碱基错误,因此,假设这少量的错误在模型求解时时在误差允许的范围内的。
5.在基因组的剪切过程中未发生基因的丢失,DNA改变,基因的重叠等
4.1.2数据在拼接的预处理
Reads在拼接时,由于新一代序列数据很多,准确度较低,导致reads中含有大量错误碱基。
在这种错误下,deBruijn图的实际大小会随着reads数据量的增加呈现指数型增长,并且容易造成错误拼接。
因此,在此之前需要对reads进行预处理,修正或消除初始reads中的碱基错误。
(1)新一代测序数据错误率高,且主要分布在靠近reads3’端部分,并且越靠近3’出错率越高,而5’端比较正确,如图3-1[2]所示。
为减少错误,我们的方法是:
计算3’端reads长度一般的碱基的平均质量,过滤掉该区域平均质量小于15的reads。
该质量值对应碱基的出错率,计算公式如下:
Q=-10×lg(ε)。
其中,Q为碱基的质量只值,ε为碱基的出错率。
图4-1solexa测序数据的错误特征
(2)测序数据中会有许多全A或者基本上全A的reads,这些数据很可能是solexa测序过程中的人工数据,需要去除。
方法为:
该定A的含量阀值为0.9,过滤掉A含量大于0.9的reads。
(3)测序数据中含有一些未知的碱基,通常同“N”或“.”表示,代表没有被测出来的碱基,对拼接有不利的影响,因此含有未知碱基的reads需要过滤掉。
4.2算法的步骤
4.2.1deBruijn图策略
基于deBruijn图策略的拼接算法的做法是:
(1)构建deBruijn图:
将reads分割成一系列连续的子串k-mers(一般用k只表示kmer碱基数目的大小),作为图中的边,相邻的两个k-mers交叠(k-1)个碱基;
(2)化简deBruijn图:
方法是合并路径出度入度唯一的节点,按照一定的规则去除图中的尖端(tips)和泡状结构(bubbles),如图3-2所示[3];
(3)构建contigs:
在deBruijn图或其子图中寻找一条最优的欧拉路径(经过每条边一次且仅一次的路径),该路径对应的碱基序列即为contigs;
(4)生成scaffolding:
利用配对数据,确定contigs之间的相对方向与位置关系,对contigs进行组装,并填充contigs之间的gaps,最终得到scaffolds序列。
图4-2deBruijn图的简化
4.2.2deBruijn图的构建
与其他基于deBruijn图拼接算法不同的是,基于reads引导的基因序列拼接以deBruijn图架构为基础,这里所建立的deBruijn图只有孤立的顶点,而没有边的存在。
基于deBruijn图的数据存储方法具有快速高效的优点,然而内存开销比较大,内存是计算机的宝贵资源,直接制约着拼接的正常进行,因此,设计新的内存高效的存储方法,对拼接具有重要的影响。
(1)能够灵活操作内存,我们采用C语言作为拼接算法的实现工具。
(2)一条reads,不额外生成其反向互补,而是在拼接过程中,自动检测该reads是正向还是反向。
这样就避免了数据的二倍增长,以及不必要的数据冗余和计算资源的浪费。
(3)设计高效的数据结构,deBruijn图的主要数据结果如图4-3[4]所示。
图(a)表示的是如何在一个read上获取kmers,kmer是reads上的连续子串,其中一般用k值表示kmer的大小。
从reads上的首个碱基开始,依次截取k个碱基。
由于reads数据的高覆盖度及k值一般比较小,这就造成在获取kmer的过程中会出现大量一样的kmer。
因此,就没个kmer而言,需要记录出现的reads编号和reads上几号位置,还需要对其出现的次数累加。
图(b)展现了reads1和reads2在deBruijn图中的存在形式,每个read在deBruijn图中都对应着一个kmer通路。
图(c)表示的是deBruijn图的数据结构。
图4-3deBruijn图
4.2.3contigs构建
综合考虑参与拼接的reads的各项信息,从reads整体角度研究拼接问题,是基于reads引导的全基因组段序列拼接算法的基本思想。
其主要做法是:
将正在参与拼接的reads的拼接信息记录到决策表中,综合考虑决策表中reads累计拼接信息和数据的出错特征,设计合理的评分方法,选择得分最大的kmer进行contig的扩展。
根据选择的kmer更新决策表中reads的信息,contigs每次扩展一个碱基,第一轮扩展结束后,取其反向互补,进行反向扩展。
基于reads引导的contigs构建过程如下:
1.轮拼接,初始化决策表为空,在deBruijn图中选择拼接的初始kmer作为初始拼接的contigs,将参与拼接的reads信息添加进决策表;
2.扩展候选的k-mers,根据决策表中的reads信息,对这些kmer评分,选择得分最高的k-mer;
3.候选k-mers为空,并且已经处于第二轮拼接,则拼接停止,转到第6步;如果仍处于第一轮拼接,则标记为第二轮拼接,初始化决策表为空,将contigs反向互补,在新的3’端重新选择k-mers进行拼接;
4.选择的k-mer更新决策表中的reads的拼接信息,删除拼接成功或者失败的reads;
5.在deBruijn图中标记拼接成功的reads的所有k-mers为删除状态,将成功的reads添加进contigs;
6.拼接结束,则保存contig及成功的reads。
该算法的拼接示意图如图4-4[5]所示。
经过打分导航机制确定了一个带扩展的新kmer。
该kmer在rid为3的当前read中没有出现,因此它的mismatches增1,lastpos变为0.而该kmer在rid为6的kmer中出现了,则它的matches增1,lastpos增1,变为14.如果此时rid为6的lastpos不是13而是0,则last由0变为14.rid为2和8的read是新引进决策表的read。
图4-4contigs拼接示意(read长为36bp,k值取12)
4.2.4contigs组装
本阶段以reads拼接阶段生成的contigs和配对数据文件作为输入,通过配对数据确定contigs之间的相对位置,连接contigs,填充相邻contigs间的gaps区域等过程,最终生成长度更长的scaffolds,如图4-5[6]所示。
图4-5contigs组装示意图
4.3.1模型配对文库参数的校正:
在contigs组装阶段,配对数据常用语确定相邻两个contigs的相对方向和位置关系,校正contigs中的序列重排,以及填充相邻的两个contigs间gaps区域,最终生成质量更高的超长序列scaffolds。
在该阶段,reads被映射到contigs上,以获得配对reads在contigs中的位置。
如果配对的两条reads分别出现在不同的contigs中,则说明这两个contigs在目标基因组中是相邻的。
Contigs组装的关键是利用配对数据确定contigs的相对方向和位置关系,通常采用基于图的方法解决该问题。
在图中,contigs作为顶点,contigs之间的配对连接作为边,配对俩接的数量作为边上的权重,这样就建立了任意两条contigs之间的配对关系,该问题就转化为在图中或子图中寻找遍历每个节点的最优路径。
通常与给定的contigs具有配对关系的contigs数量是很少的,因此,该图并不复杂,处理起来相对简单。
然而,在实际中,由于目标基因组中的重复片段以及可能测序错误的影响,contigs之间的连接有一些是不真实的(或错误的),为了保证结果的准确性,拼接算法通常忽略配对数量低于一定阀值的contigs之间的连接。
通常,配对文库的片段大小服从正态分布,并且标准差大约在均值的10%左右。
在实际应用中,往往只给出配对文库的均值,而标准差是未知的。
然而,配对文库的片段尺寸的分布通常与事先给定的分布并不完全相同,总会有一些差距。
因此,我们对配对数据的分布参数进行了校正。
首先,随机选取若干条contigs作为参数校正的contigs样本并对数据文件的reads进行编号。
一般而言,配对数据会存储在两个文件里,对数据文件1的reads从1开始用奇数编号,对数据文件2的reads从2开始用偶数编号,这样,1、2号是配对的,3、4号是配对的,以此类推。
然后,对选取的contigs按read的长度建立索引。
把两个文件里的reads往contigs上映射。
通常,如果reads数据文件太大的话,可以截取文件1和文件2相对应的各一部分,再往选取的contigs上映射。
最后,分析映射结果,计算配对数据的分布参数,提取reads数据的映射信息,并按reads编号由小到大排列,这样便于分析映射到每一条contigs上的配对数据之间的相关信息。
图3-5-1[7]展现了对两对配对数据映射到contigs上的情形。
图4-5-1配对数据映射到contigs
4.3.2contigs相对位置的确定
在contigs组装之前,需要先无额定contigs之间的相对位置和方向。
将配对数据往reads拼接阶段生成的contigs上映射,通常会有若干配对的reads映射到不同的contigs上,如果被配对的reads映射上的两条contigs之间有交叠,则看配对的两条reads在它们上的位置之间的距离是否在合理范围内。
如果距离在合理范围内且配对reads数目达到设定的阀值,则认为两条contigs是相邻的,应该连接到一起。
如果这两个contigs没有交叠,则认为它们之间存在间隙gaps,需要后续的gaps填充操作。
在contigs的构建阶段,拼接错误通常发生在contigs的末端,这是因为contigs是因找不到符合田间的下一个kmer才结束拼接的,contigs靠近末端的若干碱基正确率较低。
因此,我们摒弃了contigs组装通常采用的把配对reads数据往整条contigs上映射的做法,而只取contigs的两端各Lbp长的序列片段,这样既减少了内存消耗,也有利于提高计算速度。
如果contigs的长度小于L,则取其整个序列。
L值的确定于配对文库的平均插入距离有关,一般取值要大于配对文库的平均插入距离,比如如果所用配对文库的平均插入距离为200bp的话,那么L就取值300。
(1)建立contigs索引;
(2)映射配对reads和contigs之间的匹配信息;
(3)确定contigs的相对位置。
4.3.3contigs连接
contigs相对位置的确定指的是contigs排放位置的确定,并没真正的连接在一起。
contigs连接时需要先计算contigs之间的交叠overlap,由于contigs之间的相对位置已经确定,我们可以根据相邻两条contigs上配对的reads来计算contigs之间的距离distance,从而确定overlap的大小。
此时的distance大小可以大于0,也可以小于0。
若计算出的distance值小于0,则表明这两条contigs之间可能有overlap存在;若计算出的distance值大于0,则表明这两条contigs之间可能有真正的gaps存在。
4.3.4gaps的填充
contigs之间的相对位置确定后,在contigs连接阶段,有交叠的contigs被连接到一起,合并成长度更长的contigs,没有被连接的相邻contigs之间有空隙gaps的存在,需要进行gaps填充。
gaps填充就是通过局部序列拼接的方式将contigs之间缺少的那部分碱基序列构造出来。
配对数据往拼接生成的contigs上映射时,有这样的一部分配对的reads:
它们只有其中的一条找到了匹配位置,而与之配对的另一条没有找到,我们称配对reads中在contigs上找到匹配位置的那些reads为悬浮hangingreads。
针对每两条相邻的contigs,我们将悬浮的那些reads保存在一个表中,该表记为悬浮reads表HRL。
与此同时,收集那些与hangreads对应的在contigs上未找到匹配位置的reads,并用它们构建出一个DeBruijn图。
有了悬浮reads表HRL和所构建出来的DeBruijn图,我们就可以很快速的对contigs之间的gaps进行填充。
如图4-5-4.1[8]所示。
图4-5-4.1gaps填充示意图
如图4-5-4.2[9]所示,选取contigsA末端READ_LEN长的序列处开始拼接,每次扩展一个碱基。
待拼接kmer有效的条件:
悬浮reads表中存在与kmer所在read的配对read,并且该kmer所在的read的方向与contigs的扩展方向相反以及该配对reads的距离在配对文库的误差范围之内。
最优kmer的确定方法与reads拼接阶段的kmer选取方法一样,就是通过打分机制,对待选的kmer进行打分,取得分最大者。
在reads拼接阶段,如果下一个kmer为空时,则停止contigs的首轮扩展,并进行contigs的反向扩展,而在填充gaps时,当下一个kmer为空时,我们把N作为扩展字符,只进行单向扩展。
图4-5-4.2局部拼接填充gaps
为了加快拼接速度,我们在实际操作时选择从gaps两端同时相向拼接,即在相邻contigs对应的两末端同时开始在悬浮reads表和DeBruijn图的引导下对g
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 基因组 序列 拼接