基于贪心算法与最短路径的基因组组装最优拼接问题---1411.doc

基于贪心算法与最小路径的基因组组装优化问题摘要 随着人类基因组计划的实施和飞速发展，基因组测序拼接作为生物信息学的核有着极其重要的应用价值。新的测序技术大量涌现，产生的reads长度更短，数量更多，覆盖率更大，能直接读取的碱基对序列长度远小于基因组长度。本文通过如何在保证组装序列的连续性、完整性和准确性的同时设计耗时短、内存小的组装算法，建立数学模型来解决基因组组装问题。针对问题一，首先，利用相应的软件对原基因组G进行切割，利用全基因鸟枪法测序对切割后的短基因进行测序，得到较小的基因组,通过对比多条任意切割后相似的基因组从而找出个别碱基对存在的识别错误。而对于基因组中存在的重复片段可以通过两个read之间的DNA片段的长度满足一定的分布规律即pared end read来解决。 接下来对比任意两个和是否相等，通过MATLAB软件建立nm阶的关联矩阵，最后利用图论中的最短路径方法使更多的基因组能拼接在一起，尽可能使拼接出来的基因组在原基因组的覆盖率达到最大。 针对问题二，先把附件给出的数据提取出来导入MATLAB中，再结合问题一给出的模型对基因组进行重组，从而得到新的基因。 最后，基于对基因组组装的研究，为使重组基因能更接近原基因序列，对问题一提出模型进行合理性的评价。关键词：基因组组装 全基因鸟枪法测序 贪心算法 最短路径 一、 问题的重述1.1问题背景快速和准确地获取生物体的遗传信息对于生命科学研究具有重要的意义。对每个生物体来说，基因组包含了整个生物体的遗传信息，这些信息通常由组成基因组的DNA或RNA分子中碱基对的排列顺序所决定。获得目标生物基因组的序列信息，进而比较全面地揭示基因组的复杂性和多样性，成为生命科学领域的重要研究内容。1.2问题提出确定基因组碱基对序列的过程称为测序（sequencing）。测序技术始于20世纪70年代，伴随着人类基因组计划的实施而突飞猛进。从第一代到现在普遍应用的第二代，以及近年来正在兴起的第三代，测序技术正向着高通量、低成本的方向发展。尽管如此，目前能直接读取的碱基对序列长度远小于基因组序列长度，因此需要利用一定的方法将测序得到的短片段序列组装成更长的序列。通常的做法是，将基因组复制若干份，无规律地分断成短片段后进行测序，然后寻找测得的不同短片段序列之间的重合部分，并利用这些信息进行组装。例如，若有两个短片段序列分别为ATACCTTGCTAGCGTGCTAGCGTAGGTCTGA则有可能基因组序列中包含有ATACCTTGCTAGCGTAGGTCTGA这一段。当然，由于技术的限制和实际情况的复杂性，最终组装得到的序列与真实基因组序列之间仍可能存在差异，甚至只能得到若干条无法进一步连接起来的序列。对组装效果的评价主要依据组装序列的连续性、完整性和准确性。连续性要求组装得到的（多条）序列长度尽可能长；完整性要求组装序列的总长度占基因组序列长度的比例尽可能大；准确性要求组装序列与真实序列尽可能符合。利用现有的测序技术，可按一定的测序策略获得长度约为50100个碱基对的序列，称为读长（reads）。基因组复制份数约为50100。基因组组装软件可根据得到的所有读长组装成基因组，这些软件的核心是某个组装算法。常用的组装算法主要基于OLC（Overlap/Layout/Consensus）方法、贪婪图方法、de Bruijn图方法等。一个好的算法应具备组装效果好、时间短、内存小等特点。新一代测序技术在高通量、低成本的同时也带来了错误率略有增加、读长较短等缺点，现有算法的性能还有较大的改善空间。具体解决问题如下：问题一：试建立数学模型，设计算法并编制程序，将读长序列组装成基因组。你的算法和程序应能较好地解决测序中可能出现的个别碱基对识别错误、基因组中存在重复片段等复杂情况。问题二：现有一个全长约为120,000个碱基对的细菌人工染色体（BAC）， 采用Hiseq2000测序仪进行测序，测序策略以及数据格式的简要说明见附录一和附录二，测得的读长数据见附录三，测序深度（sequencing depth）约为70，即基因组每个位置平均被测到约70次。试利用你的算法和程序进行组装，并使之具有良好的组装效果。二、 问题分析2.1 问题一分析本题要求我们的算法和程序应能较好地解决测序中可能出现的个别碱基对识别错误、基因组中存在重复片段等复杂情况。故在下列分别对个别碱基识别错误和基因组中存在重复片段进行分析。2.1.1个别碱基对识别错误分析read 中每一个碱基都有一个质量值，来表示该碱基被正确测出的概率。一般来说，5端的碱基正确的概率较大，而 3端 1 到 3 个碱基可能是错误的。这就要求拼接软件在拼接时能够纠错，但是，可纠错的软件也可能把正确的碱基当作错 误来纠正。所以不仅要求拼接软件在拼接时能够纠错，尽可能多的发现真正的错误，而且要求拼接软件尽可能少的将正确的碱基识别成错误的。2.1.2基因重复片段分析 基因组中存在大量重复片段，重复片段可能导致拼接错误，或者导致不连续的较短contig出现。重叠片段类型主要有以下几种，如下图所示。图1 基因组重叠片段类型图2.2问题二分析本题题目提供全长约为120,000个碱基对的细菌人工染色体，采用新一代的Hiseq2000测序仪进行测序。附件提供了筛选好的定长reads数据文件。先将附件的数据提取出来储存到空文件A中，再将之导入到MATLAB中。然后使用第一题提出的基于贪心算法与最短路径算法的组装算法的模型中，得出新的基因组G，并对结果进行误差分析。三、 问题假设(1)假设测序过程中没有其他因素的干扰；(2)假设题目所给定的序列相对位置的碱基全部遵循GU-AC法则；(3)假设题目中所有的序列都是正常可判别的序列，没有出现序列的基因突变等情况；(4)假设一个完整基因组，打断成500bp的片段是随机的；(5)假设基因组每个位置被测到的几率是等可能的；(6)所有片段上的碱基都已经被识别出来，不存在未知碱基。四、 模型符号说明原基因进行第j-1次复制并对其进行任意切割后的第i个基因 基因的长度，即有个碱基对K碱基对数量，即有K个碱基对第一个碱基对到第K个碱基对组成的基因第-K+1个碱基对到第个碱基对组成的基因第一个碱基对到第-K个碱基对组成的基因第K+1个碱基对到第个碱基对组成的基因Conting（C）由经过贪心算法和最短路径算法后拼接产生的基因从顶点到的一条路的权.也就是的值的父亲点，用以确认最短路的路线.S 具有永久标号的顶点集.五、 模型的建立及求解5.1 基因组序列的获取及拼接 针对新一代测序数据reads 长度较短、数据海量的特点，全基因组测序方面的数据分析软件的研发，已成为生物信息学领域最迫切、最重要的研究课题。基于新一代测序数据的基因组序列拼接，通常分为如下三个阶段： （1）数据的预处理阶段。该阶段通过特定的方法，移除测序数据中的错误碱基； （2）基因组连续片段（contigs）生成阶段。该阶段将reads 拼接成contigs； （3）超长序列片段（scaffoldings）组装阶段。该阶段使用配对数据，确定contigs 之间的方向和位置关系，生成scaffoldings。5.1.1 全基因组鸟枪测序法随着人类基因组计划的完成，人类对自身遗传信息的了解和掌握优乐前所未有的进步。与此同时，分子水平的基因检测技术平台不断发展和完善，使得基因检测技术得到了迅猛发展，基因检测效率不断提高。从最初第一代以Sanger测序为代表的直接检测技术和以连锁分析为代表的间接测序技术，其最主要的测序方法是全基因组鸟枪法（WGS）测序。基因组研究的核心目标是获得生物体的整套遗传密码.其实现的技术途径是大规模DNA测序。图2 WGS测序的步骤该测序方法的优点在于:其一,大大缩短侧序周期并降低了侧序成本.由于在构建基因组序列的整个过程中无需任何物理图谱,节省了大量的时间,以及人工实验所播要花费的大量成本；其二,双端侧信息,可以有效的排除r印eat区域的对整个拼接过程的影响:这样就缩短6nishnig阶段大量人力财力的投入.而缺点在于:拼接过程中计算量明显加大,对软硬件性能要求较高;双端测信息的引入,需要引入新的算法,并改进现有的软件。5.1.2 贪心算法 贪心算法（又称贪婪算法）是从问题的某一个初始解出发逐步逼近给定的目标，以尽可能快的地求得更好的解，当达到某算法中的某一步不能再继续前进时，算法停止。该算法存在问题： 1. 不能保证求得的最后解是最佳的； 2. 不能用来求最大或最小解问题； 3. 只能求满足某些约束条件的可行解的范围。 结合该题，贪心策略类型的序列拼接算法主要采用种子迭代扩展的方法，按一定条件选择初始reads 作为待生成contigs 的种子，通过启发式搜索方式使得每一步都合并与其具有最多交叠的reads，直至reads 或contigs 两端都不能再做进一步的扩展。一般而言，reads 的选择是按照拼接质量递减的顺序考虑的，拼接质量通常用碱基质量和覆盖度来衡量。为避免错拼，有些扩展操作在发现冲突的信息时就立即停止。SSAKE16、SHARCGS11、VCAKE9即采用了该类拼接策略。SSAKE 和VCAKE能够处理非完全匹配的reads，SHARCGS适用于均匀分布、非配对的reads。5.1.3测序数据分析对测序仪测出的数据进行分析，我们发现如下特征：(1) 基因组中有些位置被较多reads所覆盖，有些位置被较少reads所覆盖，这些位置是随机的，不可预知。(2) 每个reads的每个碱基都有一个质量值，该质量值能反映该碱基的正确率。质量值越高，则碱基的正确率越高。(3) 有些reads上某个碱基含量过高（超过90%），甚至所有碱基都是某一个碱基，这样的reads错误率较高，在拼接过程中应尽力避免使用该类reads。5.2 问题一的模型建立及求解令=其中1in,1jm.设，i。,j且1，n,1,m。设A为阶空矩阵。若=，则A(n+,n+)=1若，则A(n+,n+)=+则可得关于迪克斯特拉（Dijkstra）算法的邻接矩阵A。接下来利用最短路径算法和贪心算法求出基因重组的最佳方案。定义为节点图，中的任意一点皆为一个结点，则共有nm个点，其中第n+1个点表示，即以n个点为一周期，以此类推。并从其中随意定义一点为初始点. 对每个顶点，定义两个标记（,),其中：= 算法的过程就是在每一步改进这两个标记，使最终为从顶点到的最短路的权.输入为带权邻接矩阵W.（1） 赋初值：令S=，=0.=VS,令=W（，）,=,.（2） 更新，:=VS,若+,则令=+,=（3） 设是使取得最小值的中的顶点，则令S=S，.（4） 若，转（2）；否则，停止.3.2问题二的模型建立及求解根据新一代全基因组鸟枪法测序过程，本文建立了如下模型：分别从500bp片段的两端，对两条单链进行测序，测得的读长记为reads1，reads2。reads1，reads2的长度均为88bp，且该对reads相距500bp。由reads2通过碱基互补配对原则可得reads3,再把碱基ATCG分别记为1、2、3、4，则可构成188的矩阵。通过reads1与reads3的矩阵相减可判断该条单链是否可以重组，另一条单链可由碱基互补配对获得，再对数据进行拟合，通过拟合的效果可判断组装的正确率。七、 模型的评价与改进7.1 模型的评价7.1.1 模型的优点(1) 本模型的算法容易推广到实际的基因组组装中，具有一定的实际应用价值；(2) 本模型针对新一代测序技术出现的问题逐一进行了解决。新一代分割的读长较短，不易发现reads之间的重叠关系，本文放弃了传统的重叠图算法。把基因重组问题成功的转化为Dijkstra的问题，配合全基因组鸟枪法测序实现了在较短的时间内对大量reads数据进行比对处理。7.1.2 模型的不足(1)忽略了碱基存在的内环境因素及其生化结构的影响；(2)在实际中，基因组组装是一个复杂的数学问题，存在着大量的不确定性；(3)一些新的想法缺乏足够的理论根据，所以有些问题的解决带有一定的主观性；(4)该模型处理重复片段的能力较弱。当有reads拼接失败时，意味着小于reads长度的重复片段被检验出来，但无法处理大于reads长度的重复片段。7.2 模型的改进 （1）对于处理重复片段弱的情况，尝试使用覆盖度分析，建立相应对的概率模型，必要时使用配对reads信息，来发现大于reads长度的重复片段。 （2）基因测序仍是一个探究性问题，虽然测序方法比较多，但对于基因组而言，现有的方法在精度上还不够令人满意，需要更多研究人员来投入时间去完善。 参考文献1孔波. 几种新化学/生物分析方法及分析数据建模算法的研究与应用D.湖南大学,2010.2徐静怡. 基因组序列拼接算法及ncRNA新基因的发现D.中国科学院研究生院（计算技术研究所）,2004.3逯雯雯,卢志远,王亚旭,孙啸. 面向新一代基因组测序技术的序列拼接算法J. 生物信息学,2010,03:248-253.4曾培龙. 基于reads引导的基因组序列拼接D.哈尔滨工业大学,2012.5罗春清,杨焕明. 微生物全基因组鸟枪法测序J.遗传,2002,03:310-314.附录package cn.data;import java.io.*;import java.sql.*;public class DataInput /*A=00*C=01*G=10*T=11*/public static void main(String args) throws SQLException,ClassNotFoundException/ TODO Auto-generated method stubFile file=new File(E:基因数据,data1.txt); /*File file2=new File(D:基因数组,oo.txt);*/String a = new String60000;String b = new String60000;String c = new String60000;DataInput myData = new DataInput();tryFileReader inOne=new FileReader (file);BufferedReader inTwo=new BufferedReader(inOne); /*FileWriter outOne=new FileWriter(file2);*/*BufferedWriter outTwo=new BufferedWriter(outOne);*/String s=null;int i;int j;int k;int p;for(j=0,k=0,p=0,i=1;(s=inTwo.readLine()!=null;i+) if(i%4 = 1)aj = s;j+;else if(i%4 = 2)bk = myData.StringToByte(s);k+;else if(i%4 = 0)cp = s;p+;/* System.out.println(i+:+s);outTwo.write(i+:+s);/write()函数用法outTwo.newLine();*/inOne.close();inTwo.close();/* outTwo.flush();/outOne.close();outTwo.close();*/catch(Exception e)System.out.println(e);for(int i = 0;i100;i+)System.out.println(a+i+ai);for(int i = 0;i100;i+)System.out.println(bi);System.out.println(a40000+a40000);try Class.forName(com.mysql.jdbc.Driver);String url = jdbc:mysql:/localhost/genmdata?user=cyh&password=1011325;Connection conn =DriverManager.getConnection(url); /Statement stmt = conn.createStatement();String sq