基因组序列拼接

2014年成都理工大学校内数学建模竞赛论文二0一四年五月二十五日摘要：本文所要研究的就是全基因组的从头测序的组装问题。

首先，本文简要介绍了测序技术及测序策略，认真分析了基因系列拼装所面临的主要挑战，比如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片段；一.问题重述基因组组装快速和准确地获取生物体的遗传信息对于生命科学研究具有重要的意义。

对每个生物体来说，基因组包含了整个生物体的遗传信息，这些信息通常由组成基因组的DNA或RNA分子中碱基对的排列顺序所决定。

获得目标生物基因组的序列信息，进而比较全面地揭示基因组的复杂性和多样性，成为生命科学领域的重要研究内容。

确定基因组碱基对序列的过程称为测序（sequencing）。

测序技术始于20世纪70年代，伴随着人类基因组计划的实施而突飞猛进。

从第一代到现在普遍应用的第二代，以及近年来正在兴起的第三代，测序技术正向着高通量、低成本的方向发展。

尽管如此，目前能直接读取的碱基对序列长度远小于基因组序列长度，因此需要利用一定的方法将测序得到的短片段序列组装成更长的序列。

通常的做法是，将基因组复制若干份，无规律地分断成短片段后进行测序，然后寻找测得的不同短片段序列之间的重合部分，并利用这些信息进行组装。

例如，若有两个短片段序列分别为ATACCTT GCTAGCGTGCTAGCGT AGGTCTGA则有可能基因组序列中包含有ATACCTT GCTAGCGT AGGTCTGA这一段。

当然，由于技术的限制和实际情况的复杂性，最终组装得到的序列与真实基因组序列之间仍可能存在差异，甚至只能得到若干条无法进一步连接起来的序列。

对组装效果的评价主要依据组装序列的连续性、完整性和准确性。

连续性要求组装得到的（多条）序列长度尽可能长；完整性要求组装序列的总长度占基因组序列长度的比例尽可能大；准确性要求组装序列与真实序列尽可能符合。

利用现有的测序技术，可按一定的测序策略获得长度约为50–100个碱基对的序列，称为读长（reads）。

基因组复制份数约为50–100。

基因组组装软件可根据得到的所有读长组装成基因组，这些软件的核心是某个组装算法。

常用的组装算法主要基于OLC（Overlap/Layout/Consensus）方法、贪婪图方法、de Bruijn 图方法等。

一个好的算法应具备组装效果好、时间短、内存小等特点。

新一代测序技术在高通量、低成本的同时也带来了错误率略有增加、读长较短等缺点，现有算法的性能还有较大的改善空间。

问题一：试建立数学模型，设计算法并编制程序，将读长序列组装成基因组。

你的算法和程序应能较好地解决测序中可能出现的个别碱基对识别错误、基因组中存在重复片段等复杂情况。

问题二：现有一个全长约为120,000个碱基对的细菌人工染色体（BAC），采用Hiseq2000测序仪进行测序，测序策略以及数据格式的简要说明见附录一和附录二，测得的读长数据见附录三，测序深度（sequencing depth）约为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，用加权的边来判断到底要和哪个片段进行拼接，我们在查阅资料后，发现可以通过de bruijn图并对其进行相应的改进后来建立数学模型对问题进行求解。

设想一本杂志被复制成多份,将每份杂志均以不同的方式剪切,将多份剪切的杂志放在一起。

在剪切的过程中,一些碎片丢失,一些碎片被污渍浸染,一些碎片存在着重叠现象。

根据上述情况来寻找恢复原始杂志的方法。

这是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和EULER等。

四．模型建立4.1.1模型的假设1.假设模型中的read片段都是由一条完整的DNA经过测序而来，它们进过拼接后可以形成一个大片段。

2.模型中出现的各个序列中DNA的双链都准确3.模型中read在拼接时合理地去掉的公共部分在误差允许的范围内。

4.由于总会在测序中出现read的碱基错误，因此，假设这少量的错误在模型求解时时在误差允许的范围内的。

5.在基因组的剪切过程中未发生基因的丢失，DNA改变，基因的重叠等4.1.2数据在拼接的预处理Reads在拼接时，由于新一代序列数据很多，准确度较低，导致reads中含有大量错误碱基。

在这种错误下，de Bruijn图的实际大小会随着reads数据量的增加呈现指数型增长，并且容易造成错误拼接。

因此，在此之前需要对reads 进行预处理，修正或消除初始reads中的碱基错误。

（1）新一代测序数据错误率高，且主要分布在靠近reads3’端部分，并且越靠近3’出错率越高，而5’端比较正确，如图3-1[2]所示。

为减少错误，我们的方法是：计算3’端reads长度一般的碱基的平均质量，过滤掉该区域平均质量小于15的reads。