《基于贪心算法与最短路径的基因组组装最优拼接问题---1411》由会员分享，可在线阅读，更多相关《基于贪心算法与最短路径的基因组组装最优拼接问题---1411（12页珍藏版）》请在金锄头文库上搜索。
1、 基于贪心算法与最小路径的基因组组装优化问题基于贪心算法与最小路径的基因组组装优化问题摘要摘要 随着人类基因组计划的实施和飞速发展，基因组测序拼接作为生物信息学 的核有着极其重要的应用价值。新的测序技术大量涌现，产生的 reads 长度更 短，数量更多，覆盖率更大，能直接读取的碱基对序列长度远小于基因组长度。 本文通过如何在保证组装序列的连续性、完整性和准确性的同时设计耗时短、 内存小的组装算法，建立数学模型来解决基因组组装问题。 针对问题一，首先，利用相应的软件对原基因组 G 进行切割，利用全基因 鸟枪法测序对切割后的短基因进行测序，得到较小的基因组,通过对比多条j iG任意切割后相似的基因
2、组从而找出个别碱基对存在的识别错误。而对于基因j iG组中存在的重复片段可以通过两个 read 之间的 DNA 片段的长度满足一定的分布 规律即 pared end read 来解决。接下来对比任意两个和是否相等，通过 MATLAB 软件建立111mnread322mnreadn m 阶的关联矩阵，最后利用图论中的最短路径方法使更多的基因组能拼接在 一起，尽可能使拼接出来的基因组在原基因组的覆盖率达到最大。针对问题二，先把附件给出的数据提取出来导入 MATLAB 中，再结合问题一 给出的模型对基因组进行重组，从而得到新的基因。最后，基于对基因组组装的研究，为使重组基因能更接近原基因序列，对 问
3、题一提出模型进行合理性的评价。关键词：关键词：基因组组装 全基因鸟枪法测序 贪心算法 最短路径 一、一、 问题的重述问题的重述1.1 问题背景快速和准确地获取生物体的遗传信息对于生命科学研究具有重要的意义。对每个生物体来说，基因组包含了整个生物体的遗传信息，这些信息通常由组成基因组的 DNA 或 RNA 分子中碱基对的排列顺序所决定。获得目标生物基因组的序列信息，进而比较全面地揭示基因组的复杂性和多样性，成为生命科学领域的重要研究内容。1.2 问题提出确定基因组碱基对序列的过程称为测序（sequencing）。测序技术始于 20 世纪 70 年代，伴随着人类基因组计划的实施而突飞猛进。从第一代
4、到现在普遍 应用的第二代，以及近年来正在兴起的第三代，测序技术正向着高通量、低成 本的方向发展。尽管如此，目前能直接读取的碱基对序列长度远小于基因组序 列长度，因此需要利用一定的方法将测序得到的短片段序列组装成更长的序列。 通常的做法是，将基因组复制若干份，无规律地分断成短片段后进行测序，然 后寻找测得的不同短片段序列之间的重合部分，并利用这些信息进行组装。例 如，若有两个短片段序列分别为 ATACCTTGCTAGCGT GCTAGCGTAGGTCTGA 则有可能基因组序列中包含有 ATACCTTGCTAGCGTAGGTCTGA 这一段。当 然，由于技术的限制和实际情况的复杂性，最终组装得到的
5、序列与真实基因组 序列之间仍可能存在差异，甚至只能得到若干条无法进一步连接起来的序列。 对组装效果的评价主要依据组装序列的连续性、完整性和准确性。连续性要求 组装得到的（多条）序列长度尽可能长；完整性要求组装序列的总长度占基因 组序列长度的比例尽可能大；准确性要求组装序列与真实序列尽可能符合。 利用现有的测序技术，可按一定的测序策略获得长度约为 50100 个碱基对 的序列，称为读长（reads）。基因组复制份数约为 50100。基因组组装软件可 根据得到的所有读长组装成基因组，这些软件的核心是某个组装算法。常用的 组装算法主要基于 OLC（Overlap/Layout/Consensus）方
6、法、贪婪图方法、de Bruijn 图方法等。一个好的算法应具备组装效果好、时间短、内存小等特点。 新一代测序技术在高通量、低成本的同时也带来了错误率略有增加、读长较短 等缺点，现有算法的性能还有较大的改善空间。 具体解决问题如下： 问题一：试建立数学模型，设计算法并编制程序，将读长序列组装成基因 组。你的算法和程序应能较好地解决测序中可能出现的个别碱基对识别错误、 基因组中存在重复片段等复杂情况。 问题二：现有一个全长约为 120,000 个碱基对的细菌人工染色体（BAC），采用 Hiseq2000 测序仪进行测序，测序策略以及数据格式的简要说明见附录一 和附录二，测得的读长数据见附录三，测
7、序深度（sequencing depth）约为70，即基因组每个位置平均被测到约 70 次。试利用你的算法和程序进行组装，并使之具有良好的组装效果。二、二、 问题分析问题分析2.1 问题一分析本题要求我们的算法和程序应能较好地解决测序中可能出现的个别碱基对 识别错误、基因组中存在重复片段等复杂情况。故在下列分别对个别碱基识别 错误和基因组中存在重复片段进行分析。2.1.1 个别碱基对识别错误分析read 中每一个碱基都有一个质量值，来表示该碱基被正确测出的概率。一 般来说，5端的碱基正确的概率较大，而 3端 1 到 3 个碱基可能是错误 的。这就要求拼接软件在拼接时能够纠错，但是，可纠错的软件
8、也可能把正确 的碱基当作错 误来纠正。所以不仅要求拼接软件在拼接时能够纠错，尽可能多 的发现真正的错误，而且要求拼接软件尽可能少的将正确的碱基识别成错误的。2.1.2 基因重复片段分析基因组中存在大量重复片段，重复片段可能导致拼接错误，或者导致不连 续的较短 contig 出现。重叠片段类型主要有以下几种，如下图所示。图 1 基因组重叠片段类型图2.2 问题二分析本题题目提供全长约为 120,000 个碱基对的细菌人工染色体，采用新一代 的 Hiseq2000 测序仪进行测序。附件提供了筛选好的定长 reads 数据文件。先 将附件的数据提取出来储存到空文件 A 中，再将之导入到 MATLAB
9、 中。然后使用第一题提出的基于贪心算法与最短路径算法的组装算法的模型中，得出新的基 因组 G，并对结果进行误差分析。三、三、 问题假设问题假设(1)假设测序过程中没有其他因素的干扰； (2)假设题目所给定的序列相对位置的碱基全部遵循 GU-AC 法则； (3)假设题目中所有的序列都是正常可判别的序列，没有出现序列的基因突 变等情况； (4)假设一个完整基因组，打断成 500bp 的片段是随机的； (5)假设基因组每个位置被测到的几率是等可能的； (6)所有片段上的碱基都已经被识别出来，不存在未知碱基。四、四、 模型符号说明模型符号说明j iG原基因进行第 j-1 次复制并对其进行任意切割后的第
10、 i 个基因j iL基因的长度，即有个碱基对j iGj iGj iLK碱基对数量，即有 K 个碱基对1 ijread第一个碱基对到第 K 个碱基对组成的基因j iG3 ijread第-K+1 个碱基对到第个碱基对组成的基因j iGj iLj iL12 ijread第一个碱基对到第-K 个碱基对组成的基因j iGj iL23 ijread第 K+1 个碱基对到第个碱基对组成的基因j iGj iLConting（C）由经过贪心算法和最短路径算法后拼接产生的基因j iG 11mngl从顶点到的一条路的权.也就是的值00g 11mng11m nL 11mngz的父亲点，用以确认最短路的路线. 11mn
11、gS具有永久标号的顶点集.五、五、 模型的建立及求解模型的建立及求解5.1 基因组序列的获取及拼接针对新一代测序数据reads 长度较短、数据海量的特点，全基因组测序方 面的数据分析软件的研发，已成为生物信息学领域最迫切、最重要的研究课题。 基于新一代测序数据的基因组序列拼接，通常分为如下三个阶段：（1）数据的预处理阶段。该阶段通过特定的方法，移除测序数据中的错误 碱基；（2）基因组连续片段（contigs）生成阶段。该阶段将reads 拼接成 contigs；（3）超长序列片段（scaffoldings）组装阶段。该阶段使用配对数据，确 定contigs 之间的方向和位置关系，生成scaff
12、oldings。5.1.1 全基因组鸟枪测序法随着人类基因组计划的完成，人类对自身遗传信息的了解和掌握优乐前所 未有的进步。与此同时，分子水平的基因检测技术平台不断发展和完善，使得 基因检测技术得到了迅猛发展，基因检测效率不断提高。从最初第一代以 Sanger 测序为代表的直接检测技术和以连锁分析为代表的间接测序技术，其最 主要的测序方法是全基因组鸟枪法（WGS）测序。基因组研究的核心目标是获得 生物体的整套遗传密码.其实现的技术途径是大规模 DNA 测序。图 2 WGS 测序的步骤该测序方法的优点在于: 其一,大大缩短侧序周期并降低了侧序成本.由于在构建基因组序列的整个 过程中无需任何物理图
13、谱,节省了大量的时间,以及人工实验所播要花费的大量 成本； 其二,双端侧信息,可以有效的排除r印eat区域的对整个拼接过程的影响:这 样就缩短6nishnig阶段大量人力财力的投入.而缺点在于:拼接过程中计算量明 显加大,对软硬件性能要求较高;双端测信息的引入,需要引入新的算法,并改进 现有的软件。5.1.2 贪心算法贪心算法（又称贪婪算法）是从问题的某一个初始解出发逐步逼近给定的 目标，以尽可能快的地求得更好的解，当达到某算法中的某一步不能再继续前 进时，算法停止。该算法存在问题：1. 不能保证求得的最后解是最佳的；2. 不能用来求最大或最小解问题；3. 只能求满足某些约束条件的可行解的范围
14、。结合该题，贪心策略类型的序列拼接算法主要采用种子迭代扩展的方法， 按一定条件选择初始reads 作为待生成contigs 的种子，通过启发式搜索方式 使得每一步都合并与其具有最多交叠的reads，直至reads 或contigs 两端都不 能再做进一步的扩展。一般而言，reads 的选择是按照拼接质量递减的顺序考 虑的，拼接质量通常用碱基质量和覆盖度来衡量。为避免错拼，有些扩展操作 在发现冲突的信息时就立即停止。SSAKE16、SHARCGS11、VCAKE9即采用 了该类拼接策略。SSAKE 和VCAKE能够处理非完全匹配的reads，SHARCGS适用于 均匀分布、非配对的reads。5
15、.1.3测序数据分析对测序仪测出的数据进行分析，我们发现如下特征： (1) 基因组中有些位置被较多 reads 所覆盖，有些位置被较少 reads 所覆 盖，这些位置是随机的，不可预知。 (2) 每个 reads 的每个碱基都有一个质量值，该质量值能反映该碱基的正 确率。质量值越高，则碱基的正确率越高。 (3) 有些 reads 上某个碱基含量过高（超过 90%），甚至所有碱基都是某一 个碱基，这样的 reads 错误率较高，在拼接过程中应尽力避免使用该类 reads。5.2 问题一的模型建立及求解令=其中 1in,1jm.j iGmnm nmmmnnGGGGGGGGGGGGLMLMMMLL3
16、2122 32 22 111 31 21 1 设，i。,j 且 1，n,1,m。1n2n1m2m1n2n1m2m设 A 为阶空矩阵。mn若=，则 A(n+,n+)=1311mnread122mnread1n1m2n2m若，则 A(n+,n+)=+311mnread122mnread1n1m2n2m则可得关于迪克斯特拉（Dijkstra）算法的邻接矩阵 A。接下来利用最短路径算法和贪心算法求出基因重组的最佳方案。定义为节点图，中的任意一点皆为一个结点，则共有 n m 个点，j iGj iGj iG其中第 n+1 个点表示，即以 n 个点为一周期，以此类推。并从其中随意定义2 1G一点为初始点.00g对每个顶点，定义两个标记（,),其中：= 11mngl 11mngz 11mngl11m nL算法的过程就是在每一步改进这两个标记，使最终为从顶点到 11mngl00g的最短路的权.输入为
