基因组序列拼接 - 图文

基因组组装模型论文

2014年成都理工大学 校内数学建模竞赛论文

题 目 编 号 队编号 参 赛 队 员 姓 名 张萌立 何理 张玲玲 C题-基因组组装 9 学号 201313030206 201305090108 201308050710 专业 计科 空间 财务管理

二0一四年五月二十五日

1

基因组组装模型论文 基因组组装

队编号：9 摘要：本文所要研究的就是全基因组的从头测序的组装问题。

首先，本文简要介绍了测序技术及测序策略，认真分析了基因系列拼装所面临的主要挑战，比如reads数据海量、可能出现的个别碱基对识别错误、基因组中存在重复片段等复杂情况，探讨了当前基因组序列拼接所采用的主要策略，即OLC（Overlap/Layout/Consensus）方法、de Bruijn图方法，且深入探讨了de Bruijn图方法。

其次，针对题中问题，以一条reads为基本单位，分为reads拼接和contig组装两个阶段，其中contig是由reads拼接生成的长序列片段。Reads的拼接阶段主要包括数据预处理、de-Bruijn 图、contig构建等，而contig的组装阶段主要包括序列的相对位置的确定以及重叠部分overlap的检测，用序列比对的方法来提高拼接的精度。

最后，进行了算法的验证与性能的评价，并且针对问题2，进行了组装分析与验证，结果表明，得到的拼接基因组序列在小范围内与原基因组序列大致吻合。

关键词：基因组系列拼接； reads；de Bruijn图；contig组装；k-mer片段；

2

基因组组装模型论文

一.问题重述

基因组组装

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

确定基因组碱基对序列的过程称为测序（sequencing）。测序技术始于20世纪70年代，伴随着人类基因组计划的实施而突飞猛进。从第一代到现在普遍应用的第二代，以及近年来正在兴起的第三代，测序技术正向着高通量、低成本的方向发展。尽管如此，目前能直接读取的碱基对序列长度远小于基因组序列长度，因此需要利用一定的方法将测序得到的短片段序列组装成更长的序列。通常的做法是，将基因组复制若干份，无规律地分断成短片段后进行测序，然后寻找测得的不同短片段序列之间的重合部分，并利用这些信息进行组装。例如，若有两个短片段序列分别为

ATACCTTGCTAGCGT GCTAGCGTAGGTCTGA

则有可能基因组序列中包含有ATACCTTGCTAGCGTAGGTCTGA这一段。当然，由于技术的限制和实际情况的复杂性，最终组装得到的序列与真实基因组序列之间仍可能存在差异，甚至只能得到若干条无法进一步连接起来的序列。对组装效果的评价主要依据组装序列的连续性、完整性和准确性。连续性要求组装得到的（多条）序列长度尽可能长；完整性要求组装序列的总长度占基因组序列长度的比例尽可能大；准确性要求组装序列与真实序列尽可能符合。

利用现有的测序技术，可按一定的测序策略获得长度约为50–100个碱基对的序列，称为读长（reads）。基因组复制份数约为50–100。基因组组装软件可根据得到的所有读长组装成基因组，这些软件的核心是某个组装算法。常用的组装算法主要基于OLC（Overlap/Layout/Consensus）方法、贪婪图方法、de Bruijn图方法等。一个好的算法应具备组装效果好、时间短、内存小等特点。新一代测序技术在高通量、低成本的同时也带来了错误率略有增加、读长较短等缺点，现有算法的性能还有较大的改善空间。

问题一：试建立数学模型，设计算法并编制程序，将读长序列组装成基因组。你的算法和程序应能较好地解决测序中可能出现的个别碱基对识别错误、基因组中存在重复片段等复杂情况。

问题二：现有一个全长约为120,000个碱基对的细菌人工染色体（BAC）， 采用Hiseq2000测序仪进行测序，测序策略以及数据格式的简要说明见附录一和附录二，测得的读长数据见附录三，测序深度（sequencing depth）约为70×，即基因组每个位置平均被测到约70次。试利用你的算法和程序进行组装，并使之具有良好的组装效果。

附录一：测序策略

测序策略如下图所示。DNA分子由两条单链组成，在图中表现为两条平行直

3

基因组组装模型论文

线，两条直线上相对位置的两个碱基相互结合形成碱基对（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，用加权的边来判断到底要和哪个片段进行拼接，我们在查阅资料后，发现可以通过de bruijn图并对其进行相应的改进后来建立数学模型对问题进行求解。

4

基因组组装模型论文

设想一本杂志被复制成多份,将每份杂志均以不同的方式剪切,将多份剪切的杂志放在一起。在剪切的过程中,一些碎片丢失,一些碎片被污渍浸染,一些碎片存在着重叠现象。根据上述情况来寻找恢复原始杂志的方法。这是DNA序列拼接问题的现实模型描述。

基于de Bmijn图的序列拼接原理主要是通过构造并简化de Bmijn图结构来实现整个序列拼接的过程。

三. 基于De Bruijn图的序列拼接技术分析与比较

二十世纪八十年代末，Pevzner等人提出基于de bruijn图的算法，并首次将该算法用于DNA序列拼接。基于de bruijn图的算法的核心思是将序列拼接问题转换为人们所熟悉的欧拉路径问题。Pevzner等人认为传统的

overlap-layout-consensus算法导致了将DNA序列拼接问题转换为Hamilton路径问题，他们受到杂交测序方法SBH(Sequencing by Hybridization)的启发，创造性地提出了在de Bruijn图中寻找欧拉路径的构想，尽管杂交测序方法SBH从未在测序工程中实际应用过，但它直接引发了基因芯片工业的诞生。

构造de Bruijn图的方法如下所述：

（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作为de Bruijn图的顶点。

（2）对于给定的两条k-mer x和y，如果在某read ri中存在一条长度为k+1的子串，且该子串的前k个碱基与k-mer x(或y)精确匹配，同时该子串的后k个碱基与k-mer y(或x)精确匹配，那么该算法认为两条k-mer x和y之间存在一条公共边。

将采用上述方法构造的de Bruijn图记作G。对于read集合R={r1，r2，?，rn}中的任意一条read ri，若在de Bruijn图G中存在一条路径P，且该路径P访问ri中的每一条k-mer仅一次，则欧拉路径问题便可理解为：给定某一de Bruijn图G以及G中的路径集合P，在de Bruijn图G中确定某一条欧拉路径Q，使得路径集合P中的每一个元素都是欧拉路径Q的子路径。利用欧拉路径算法进行DNA序列拼接的主要步骤如下所述：首先利用纠错软件修正read中测序错误的碱基；然后按照上述方法构建de Bruijn图；构建deBruijn图之后，应将read集合中的所有read排列在de Bruijn图中，在deBruijn图中，每一条read均被视作一条路径；最后在de Bruijn图中寻找一条欧拉路径，使得该路径包含de Bruijn图中所有read所对应的路径。

在OLC中,在Overlap步骤中,采用了序列比对算法来寻找read之间的重叠信息,该算法的时间复杂度为0(?2),其中,《SDNA序列中read的数量。当前DNA测序数据序列越来越短,对同一个物种进行测序,其产生的read数量大大增加,这使得OLC的计算量增加;而基于deBruijn图原理的序列拼接中,抛弃了 OLC中序列比对算法,而是采用以k-mer为图中顶点构建图,从而减少了序列比对算法所消耗的时间,提高了算法的效率与overlap-layout-consensus算法相比，基于de bruijn图的算法有更低的时间复杂度，这是因为欧拉路径问题实际上是一个线性时间的问题。利用欧拉路径思想的拼接算法有EULER-SR、ALLPATHS、Velvet

5

基因组组装模型论文

首先，随机选取若干条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.2 contigs相对位置的确定

在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.3 contigs连接

contigs相对位置的确定指的是contigs排放位置的确定，并没真正的连接在一起。contigs连接时需要先计算contigs之间的交叠overlap，由于contigs

11

基因组组装模型论文

之间的相对位置已经确定，我们可以根据相邻两条contigs上配对的reads来计算contigs之间的距离distance，从而确定overlap的大小。此时的distance大小可以大于0，也可以小于0。若计算出的distance值小于0，则表明这两条contigs之间可能有overlap存在；若计算出的distance值大于0，则表明这两条contigs之间可能有真正的gaps存在。

4.3.4 gaps的填充

contigs之间的相对位置确定后，在contigs连接阶段，有交叠的contigs被连接到一起，合并成长度更长的contigs，没有被连接的相邻contigs之间有空隙gaps的存在，需要进行gaps填充。gaps填充就是通过局部序列拼接的方式将contigs之间缺少的那部分碱基序列构造出来。配对数据往拼接生成的contigs上映射时，有这样的一部分配对的reads：它们只有其中的一条找到了匹配位置，而与之配对的另一条没有找到，我们称配对reads中在contigs上找到匹配位置的那些reads为悬浮hanging reads。针对每两条相邻的contigs，我们将悬浮的那些reads保存在一个表中，该表记为悬浮reads表HRL。与此同时，收集那些与hang reads对应的在contigs上未找到匹配位置的reads，并用它们构建出一个De Bruijn图。有了悬浮reads表HRL和所构建出来的De Bruijn图，我们就可以很快速的对contigs之间的gaps进行填充。如图4-5-4.1[8]所示。

图4-5-4.1gaps填充示意图

如图4-5-4.2[9]所示，选取contigs A末端READ_LEN长的序列处开始拼接，每次扩展一个碱基。待拼接kmer有效的条件：悬浮reads表中存在与kmer所在read的配对read，并且该kmer所在的read的方向与contigs的扩展方向相反以及该配对reads的距离在配对文库的误差范围之内。最优kmer的确定方法与reads拼接阶段的kmer选取方法一样，就是通过打分机制，对待选的kmer进行打分，取得分最大者。在reads拼接阶段，如果下一个kmer为空时，则停止contigs的首轮扩展，并进行contigs的反向扩展，而在填充gaps时，当下一个kmer为空时，我们把N作为扩展字符，只进行单向扩展。

12

基因组组装模型论文

图4-5-4.2局部拼接填充gaps

为了加快拼接速度，我们在实际操作时选择从gaps两端同时相向拼接，即在相邻contigs对应的两末端同时开始在悬浮reads表和De Bruijn图的引导下对gaps进行填充。在填充gaps的同时，需要时刻检测这两条contigs之间的gaps信息，一旦gaps的大小接近0时就停止填充，以防止过填充情形的出现。

4.5 算法的小结

本节详细介绍了基于reads引导的基因组序列拼接算法，该算法的实施具体包括reads拼接和contigs组装两个阶段。reads拼接阶段主要有数据的预处理、De Bruijn图的建立和contigs构建等几部分，其中重点介绍了拼接时下一个kmer的选择策略，即以配对reads数据为拼接导航、综合考虑整条reads参与拼接的累计信息，并结合reads数据的区域特征提出打分机制，对待选的所有kmer进行打分，取得分最大者。contigs组装阶段主要介绍了contigs相对位置的确定、contigs连接和gaps填充等几方面的内容，就contigs之间序列交叠overlap和间隙gaps的界定，以及contigs末端出错碱基的处理及纠正等做了深入的探究。高效的数据结构，巧妙的拼接组装策略等，为该算法的成功实施提供了重要保障。

五.模型中算法的验证

在基于reads引导的基因组序列拼接算法理论的指导下，我们结合实际应用的需要，开发出一套基因组序列拼接系统，该序列拼接系统被命名为SRGA，经数据测试取得了不错的效果，从而验证了算法的正确性。

图5-1是用相关软件做出的基因组比对相似度的一部分。

13

基因组组装模型论文

图5-1.基因组的相似度比对

5.1算法的输入与输出

SRGA序列拼接系统主要包括两个大的功能模块，reads拼接模块和contigs组装模块。这两个模块相对独立，作用于输入数据的不同处理阶段，reads拼接模块生成的contigs是contigs组装模块的输入。在算法实现时，我们选择了C语言作为系统开发语言，并将软件在linux环境下发布。

14

基因组组装模型论文

新一代测序技术的快速发展，产生了大量的配对reads数据，而这些配对reads的数据文件就是系统的初始输入。如图4-1所示，reads配对数据通常分布在两个数据文件，在每个单独的文件里，每4行作为一个基本单位，记录着一条reads的序列信息和碱基质量信息。其中，第1行和第3行是测序信息，第2行是reads的碱基序列，第4行是reads每个碱基的质量值，reads的碱基序列和碱基质量值是基因组序列拼接所需要的数据信息。配对reads所分布的两个文件里对应的reads相互配对，即第一个文件里的第1条reads与第二个文件的第1条reads配对，第一个文件里的第2条reads与第二个文件的第2条reads配对，依次类推。

就系统的初始输入而言，除了配对reads数据文件，还有一个重要的参数需要人为设定，即K值。K值表征kmer的大小，在不同K值条件下，基因组序列的拼接效果往往相差很大。K值的选择还没有通用的方法，一般是通过综合考虑reads长度、数据覆盖度、目标基因组长度等因素而设定。

图5-2reads拼接阶段生成的contigs示意

初始输入的配对reads经过SRGA的reads拼接功能模块处理后，会生成大量的长度更长的序列片段contigs。如图4-2所示，生产的contigs包括contigs的编号，contigs的长度及contigs的碱基序列等。contigs只是序列拼接的中间结果，将与初始输入的配对reads一起作为系统contigs组装模块的输入。

15

本文来源：https://www.bwwdw.com/article/6iqt.html