基因序列拼接数学建模

1、基于欧拉路径的基因测序重组模型摘要本文基于新一代基因测序原理，优化了传统基因测序拼接模型并根据所给的数据，依据欧拉路径拼接算法建立了正确、高效的优化型基因测序重组模型，拼接了细菌人工染色体，并根据泊松过程，验证了模型的适用性与有效性。问题一，由于在测序时可能出现个别碱基对识别错误、基因组中存在重复片段等复杂情况，导致传统基因测序拼接模型重组的基因组存在不可靠、不准确、不连续和不正确等问题。为解决这些问题，构建准确高效的基因测序重组模型，做到了以下几点：（1）建立质量对比纠错模型，利用基因测序仪检测出的格式数据中的，将所有的碱基质量与标准值进行对比，剔除存在碱基的质量值不达标的错误的；（2）建立

2、重复检查剔除模型，将经过质量对比纠错模型处理的放在集中，一一对比集中的数据，找出数据中的重复片段，将多余个体剔除；（3）建立基于欧拉路径的基因测序重组模型，依据欧拉路径拼装算法，通过编程，建立基于欧拉路径的基因测序重组模型。具体过程如下图所示：既利用拼接成，再由拼接成过程的。问题二，仔细分析附录所给的数据，根据问题一建立的建立基于欧拉路径的基因测序重组模型重组基因组，利用泊松路径近似原理，对覆盖深度的基因测序数据，给出了一个岛中所含片段数的概率计算公式：以及岛长的期望的计算公式：对重组后的基因组进行精确的计算，岛中所含的片段数以及岛的长度与基于欧拉路径的基因测序重组模型完全一致，证明了模型的可

3、靠性与有效性。 关键词：基因测序、欧拉路径、纠错模型、重复剔除模型、泊松过程1、问题一重述快速和准确地获取生物体的遗传信息对于生命科学研究具有重要的意义。对每个生物体来说，基因组包含了整个生物体的遗传信息，这些信息通常由组成基因组的或分子中碱基对的排列顺序所决定。获得目标生物基因组的序列信息，进而比较全面地揭示基因组的复杂性和多样性，成为生命科学领域的重要研究内容。确定基因组碱基对序列的过程称为测序。测序技术始于20世纪70年代，伴随着人类基因组计划的实施而突飞猛进。从第一代到现在普遍应用的第二代，以及近年来正在兴起的第三代，测序技术正向着高通量、低成本的方向发展。尽管如此，目前能直接读取的碱

4、基对序列长度远小于基因组序列长度，因此需要利用一定的方法将测序得到的短片段序列组装成更长的序列。通常的做法是，将基因组复制若干份，无规律地分断成短片段后进行测序，然后寻找测得的不同短片段序列之间的重合部分，并利用这些信息进行组装。例如，若有两个短片段序列分别为 则有可能基因组序列中包含有这一段。当然，由于技术的限制和实际情况的复杂性，最终组装得到的序列与真实基因组序列之间仍可能存在差异，甚至只能得到若干条无法进一步连接起来的序列。对组装效果的评价主要依据组装序列的连续性、完整性和准确性。连续性要求组装得到的（多条）序列长度尽可能长；完整性要求组装序列的总长度占基因组序列长度的比例尽可能大；准确

5、性要求组装序列与真实序列尽可能符合。利用现有的测序技术，可按一定的测序策略获得长度约为50100个碱基对的序列，称为读长。基因组复制份数约为50100。基因组组装软件可根据得到的所有读长组装成基因组，这些软件的核心是某个组装算法。常用的组装算法主要基于方法、贪婪图方法、图方法等。一个好的算法应具备组装效果好、时间短、内存小等特点。新一代测序技术在高通量、低成本的同时也带来了错误率略有增加、读长较短等缺点，现有算法的性能还有较大的改善空间。问题一：试建立数学模型，设计算法并编制程序，将读长序列组装成基因组。你的算法和程序应能较好地解决测序中可能出现的个别碱基对识别错误、基因组中存在重复片段等复杂

6、情况。问题二：现有一个全长约为120,000个碱基对的细菌人工染色体， 采用测序仪进行测序，测序策略以及数据格式的简要说明见附录一和附录二，测得的读长数据见附录三，测序深度约为70×，即基因组每个位置平均被测到约70次。试利用你的算法和程序进行组装，并使之具有良好的组装效果。2、模型的假设1、在进行基因组测序时测序仪都是好的，并且都无故障。2、假设题目所给定的序列相对位置的碱基全部遵循GU-AC法则;3、假设题目中所有的序列都是正常可判别的序列，没有出现序列的基因突变等情况;4、假设一个完整基因组，打断成500bp的片段是随机的;5、假设基因组每个位置被测到的几率是等可能的;6、所有

7、片段上的碱基都已经被识别出来，不存在未知碱基。3、符号说明 第个测序片段的中第个碱基的质量值第个测序片段的中第个碱基的质量值碱基质量的标准值利用测序得到的长度为88个碱基对的序列（读长）以为单位，单倍体基因组的长度以为单位，随机测得的的长度测序所得的的数目比对过程中确定两个片段相邻关系最小重叠长度两个相邻相同的碱基序列占长度的比例 比对时，不需要进行匹配的比例，即 表示两个相邻片段起始位置的间距 表示目标序列上最左端未被片段覆盖的部分的长度 表示表示目标序列上最右端未被片段覆盖的部分的长度 称为每个片段的有效覆盖度 称为有效冗余度 比对过程中确定两个片段相邻关系最小重叠长度4、问题一模型分析、

8、建立与求解4.1 问题一分析根据新一代基因测序原理，要想构建一个能够快速和准确地组合基因测序片段，还原基因组真实碱基序列的模型，必须要解决以下三个问题：由基因测序仪检测出的格式数据存在错误，个别碱基对在识别是较为模糊，导致基因片段值出现误差，最终将影响重组基因组碱基序列的准确性。 ：为保证原基因组碱基序列能够被完美覆盖，基因测序仪检测出的具有数量大的特点，本题测序深度约为70×，即基因组每个位置平均被测到约70次，这使得测序数据中必将存在大量重复片段，这将导致拼接错误，或者导致连续不断的较短出现。：拼接程序的构建要满足：高效、准确和连续这三个要求。对于，可以利用基因测序仪检测出的格式

9、数据中的，计算出中每个碱基的质量值，然后制定一个质量标准值，将与进行比较，剔除存在碱基的质量值不在质量标准值内的错误的。对于，为解决的重复问题，首先将通过碱基互补配对的原则转化成，使和处于同一条单链，然后将和放到同一个集中，进行两两对比，并剔除重复的。对于，以经过前两步错误与重复数据剔除的集为基础，建立基于欧拉路径的基因测序重组模型。4.2 问题一模型建立问题一模型建立思路图如下图所示：图4-1 问题一建模思路示意图4.2.1 基因测序原理测序系统是一种高通量测序技术，其测序原理和测序系统相似，仍然是采用可逆终止法的边合成边测序技术。这种测序技术通过将基因组的随机片断附着到光学透明的表面，这些

10、片断通过延长和桥梁扩增，形成了具有数以亿计的，每个具有约1000拷贝的相同的模板，然后用4种末端被封闭的不同荧光标记的碱基进行边合成边测序。这种新方法确保了高精确度和真实的一个碱基接一个碱基的测序，排除了序列方面的特殊错误，能够测序同聚物和重复序列。测序策略如下图所示。分子由两条单链组成，在图中表现为两条平行直线，两条直线上相对位置的两个碱基相互结合形成碱基对，并且与碱基结合的碱基必为，与碱基结合的碱基必为。将一个含120000个的完整基因组，随机打断成的片段，然后对的片段进行测序。测序方法如第3步所示，分别从片段的两端，对两条单链进行测序，测得的读长记为，。，的长度均为，且该对相距。图4-2

11、 基因测序理论图4.2.2 格式以附录中第一条测序的为例，数据如表一：表4-1 格式数据展示表格FCC2J1CACXX:3:1101:2153:2476#41_112/1CGGTTAGGGTTAGCAATAAAACAAAAATCTAAAAAGAAAAGAAAAATTAGGGTTAATCAAGACCTAGATTAGATCGATTTTTTCTTTT+cddeegfbXWcffcgcfededffghcWaeffM_c_bZ_cdZZGMUZZTGZ_YGTZYYYaaababbabbbb其中第一行以开头，后面是的以及其他信息，例如上例中后的代表设备名称，代表中的第三个，代表第三个中的第个，代表该在

12、该中的坐标信息；，该标志代表样本的编号，用来区分个样本中的； 代表中的前一个。第二行为的序列。紧接着下面两行代表该的质量。第三行以开头，跟随者该的名称（一般于后面的内容相同），但有时可以省略，但一定不能省。第四行代表的质量。这一行可以详细说一下！测序仪是按照荧光信号来判断所测序的碱基是哪一种的，例如红黄蓝绿分别对应，那么一旦出现一个紫色的信号该怎么判断呢，因此对每个结果都有一个概率的问题。起初中心用来衡量该中每个碱基的质量，既，其中代表该碱基被测序错误的概率，如果该碱基测序出错的概率为，则 应该为，那么，那么对应的码为，则在第四行中该碱基对应的质量代表值即为，一般地，碱基质量从0-40，既码为

13、从 到。以上是中心采用记录测序质量的方法，起初没有完全依照中心的方法来定义测序质量，而是把换成了。 其他完全按照的定义来做。但是他这形式在某些情况下是不准确的，可以看出当测序质量很高的情况下两种形式几乎没区别，但低质量的碱基则有区别了。4.2.3 质量对比纠错模型根据附录所给格式数据，给所有测序片段按顺序编号，利用编程，将的数据转化成十进制数据，根据基因组学与应用生物学1中的基因组测序新标准，基因测序得到的碱基出错的概率为以上时既定为识别错误，即碱基质量的标准值为：设定大于的碱基为识别正确的碱基，小于的碱基为识别错误的碱基，将所有碱基的质量归纳整合，以碱基数为纵坐标，碱基质量为横坐标，作碱基质

14、量统计图，如图4-2所示：图4-3 碱基质量统计图由上图可知，基因测序得出的数据中存在识别错误，这些存在错误碱基的数据认为是不正确的，会给基因组重组造成错误。编程剔除附录所给格式数据中所有包含识别错误碱基的测序片段数据。再处理剔除错误测序片段后的数据，求得处理后的碱基的质量，并对所有碱基以质量为标准制图，如图4-3：图4-4 改后碱基质量统计图 由图可知，质量对比纠错模型很好的剔除了存在识别错误的碱基的测序片段，处理后的数据中已经没有了识别错误概率在0.1以上的碱基。（具体程序见附录一）4.2.4 重复检查剔除模型将经过质量对比纠错模型处理的测序数据片段归纳为R集，R集中的都是成功的，并存在大

15、量重复测序，由于基因测序是对一条测序片段两端相隔的两条单链上分别测取的，所以与不能直接比较。根据碱基互补配对原则。利用编程，使得转化成，且一一对应。以第一条为例如下图所示： 图4-5 碱基互补配对原则实例图在同一条单链上，只要组合就能得到原基因组的一条单链，然后再通过碱基互补配对原则就能得到原基因组。将放在集合中，利用编程，使中的所有两两比较，将所有重合的都删除，只保留一条（具体程序见附录二）。4.2.5 基于欧拉路径的基因测序重组模型（1）传统基因组测序的概率模型传统的鸟枪法测序策略是将测序的目标通过生物工艺随机切割成长度约为的片段，再把这些片段插入适当的细菌载体中保存，建立亚克隆文库。然后

16、利用细菌载体的分裂繁殖，大量复制克隆这些片段，用于测序。因为片段拼装是通过片段之间的重叠关系还原目标序列，为了保证相邻片段之间存在重叠，测序通常采用冗余抽样的技术，即测序的对象为目标的多个版本。因此，测得的序列总长事实上为目标序列长度的数倍，并且称这个比值为测序的覆盖深度。显然，覆盖深度越高，关于目标序列的覆盖率越高，并且相邻的之间相互重叠的可能性越大，但是，被测的片段相应的也就越多，这就造成测序成本的增加。因此，如何选择一个合理的覆盖深度，在确保测序覆盖率的条件下，尽可能地降低测序的费用，是鸟枪法测序需要解决的一个重要问题。传统基因组测序的概率模型的拼装原理如下图所示：图4-6 传统基因组测

17、序拼装原理Island:一条或者多条序列经过比对后重叠在一起，形成的区域覆盖基因 组的一部分，被称为，也称为岛;Contig:由两条或者两条以上的序列构成的称为;Singleton:仅仅由一条构成的称为;Gap:目标序列中没有被序列所覆盖的部分称为，也称为缝隙;Lander-Waterman假设:每条均匀地分布在目标序列上。 传统的拼装算法，是将测序获得的每一条与其他的进行两两比对，通过它们之间的重叠来确定它们之间的相邻关系，然后将进行伸展，最后对伸展的进行多重比对获得测序的目标序列。该算法的基本拼接思想是:首先将所有的构成一个有向图，每个看成一个结点，如果两个之间存在有(重叠)，那么在相应的

18、结点之间就认为存在有一条边。然后通过寻找经过每个一次且仅一次的一条路径。这样序列拼接问题转化成哈密尔顿路径问题。这种传统的哈密尔顿路径拼装算法，可以分为如下三步:1、:这一步是对所有的进行两两比对，通常采用算法，来确定两个之间是否有重叠，当然根据实际情况，可以确定一个最小匹配值，或者是对匹配情况进行打分。以确定两条之间是否是相邻的关系，由于任何两条之间都要进行比对，所以在此将花费大量的时间。2、:这一步根据之间的重叠信息形成，即将各个重叠的合并起来。这一步把相对较小的片段拼接成了较大的片段。3、:对于每一个，按照投票或者其他原则计算出一个序列，得到一完整的序列。通过以上三步，就可以把较短的序列

19、拼接成原始的目标序列。但是，一个很重要的事实不容忽视:无论测序技术和方法如何改进，测序所得的是有错误的，国外有己有多篇文章讨论对基因组测序错误率的计算方法。无论错误率大还是小，相对于海量的数据来说，中有错的总数还是相当可观的。因此在设计一个拼装算法时，一定要把测序错误考虑进去。（2）基于欧拉路径的基因测序重组模型 基于欧拉路径的基因测序重组模型的主要思想是把序列拼接问题转化为求欧拉路问题。由于传统的“交叠排列生成一致性序列”的思维模式导致了将拼接问题抽象为路径问题，所以此模型从另一种从来没有在实际测序工程中应用但是直接导致了基因芯片工业的杂交测序方法得到启示，提出了欧拉路径算法为核心的基因测序

20、重组模型，其主要拼装原理如下图所示：图4-7 欧拉路径基因拼接原理示意图在经过筛选后的集中，对所有进行编号，在的编号下建立的碱基序列矩阵。 （1）将序列两两进行比较，当两条的序列在两端存在87个连续碱基排列相同时，即 （2）由于集中的都是正确且不重复的，所以判定这两条为相邻，将两条以重合部分为重合点进行连接，形成一条新的存入集中。 （3）以此类推，不断连接集中的，直到集再也不存在可以连接的为止。此时，称剩下的为（一代岛），然后，缩小对比碱基数量，既当两条的序列在两端存在86个连续碱基排列相同时，就将这两条以这重合部分为重合点进行连接，再将新和成的序列放回到集中，重复上述操作，直到合成最后的基因

21、组，次基因组即为所求的目标基因组的基因序列。（具体程序见附录三）5、问题二模型分析、建立与求解5.1 问题二分析本题题目提供全长约为个碱基对的细菌人工染色体，采用新一代的测序仪进行测序。附件提供了筛选好的定长数据文件。使用第一题建立的基于欧拉路径的基因测序重组模型对数据进行组装，并对结果进行误差分析。5.2 模型的计算利用将附件中的基因测序碱基数据进处理，得到部分结果如附录四，此基因组总共包含119126个碱基，与120000十分接近，初步证明了模型的可行性。5.3 基于泊松过程的检验模型5.3.1 泊松过程近似当，则的值可以用近似，因此次序统计量的联合密度可以写成： （4）显然，式（1）与泊

22、松过程内有个事件发生的到达时间序列的密度函数相同。如果将视为在时刻有一个事件发生，令表示时刻之前事件发生的总数，则是一个泊松过程，即短片段的覆盖模型可视为条件下的泊松过程。显然 （5）因此 （6）即的强度。5.3.2 覆盖深度与岛的片段数的概率分析考虑如下特殊的排队系统，顾客以强度为的泊松流到达，系统中有无限多个服务员，每个顾客接受服务的时间为固定长度，如果在时刻，有个顾客正在接受服务，称系统在时刻的负荷为。如果将视为第个顾客的到达时间，而所对应的线段长度视为第个顾客的在系统中的服务时间，则覆盖模型中点的覆盖深度与系统在时刻的负荷是等价的。 令表示第个顾客的到达时间，表示连续两个顾客的到达时间

23、间隔，则独立同分布，服从参数为的指数分布。图5-1岛长示意图 如图4-7所示，在覆盖深度为的情况下，令的最小值为，否则就不能够形成岛。假设岛的第一条片段是对应的片段，因此有令表示岛的长度，则 (8)也是一个随机变量，可以看到随着的增加而增长。对于，则（11） 上面的式(8)给出了每个岛中所含read数目的期望值的计算公式，要想求得此式，关键是求出式(7)的结果，对式(7)进行变量代换，令，则原式变为： （12）特别对，不妨假设岛的第一条片段是对应的片段，于是随机变的概率分布为： （13）是一个参数为的几何分布，覆盖岛的片段平均数为，这个结果与基于欧拉路径的基因测序重组模型的结果仍然是一致的。又

24、因为时，岛的长度为: ,因此岛的平均长度为： （14）又因为 （15）所以 （16）通过比较，我们发现岛长度的期望值与基于欧拉路径的基因测序重组模型的结果完全一致。通过前面的计算分析比较，不难发现，本文给出的模型中，对于m=1的特殊情况，所求出的覆盖率、岛的个数、岛的长度、以及岛中所含的片段的个数的期望值都是基于欧拉路径的基因测序重组模型的结果完全相同。肯定了模型的正确性，以及用泊松过程近似的合理性。8、模型评价8.1 模型的优点本文针对新一代基因组测序所得的特点，参照传统的测序工艺的概率模型，提出了基于欧拉路径的基因测序重组模型。在给出模型前，先给出了在覆盖深度为的情况下岛的概念，并通过对模

25、型的分析，求出了覆盖率，岛的个数以及岛长的期望。通过与模拟数据进行比较，发现求出的理论值与其结果吻合的很好。本文同时对模型进行了泊松近似，在覆盖深度一定的情况下，给出了一个岛中所含片段数的概率计算公式，以及岛长的期望的计算公式。对于的时候，进行了精确的计算，岛中所含的片段数以及岛的长度与基于欧拉路径的基因测序重组模型完全一致，由此可以说明模型的正确性，并能以此模型为依据进行下一步的工作。除此以外，本模型还有以下优点：1、此模型解决了传统基因组测序时出现的碱基识别错误而导致的基因组合成错误以及基因组无法合成的问题。2、以此模型拼接时，首先剔除了中大量的重复项，使重组过程更加高效。8.2 模型的不

26、足本课题所取得的结论有一定理论价值，能够为新一代测序工艺及其短序列拼装算法的设计提供一些理论上的指导。但是由于时间，资料及作者本身能力的限制，目前的研究工作还不完善，仍有许多的研究工作有待继续完成:1、对模型进行泊松近似时,只给出了一些计算公式，没有求出最终的计算结果，因此后期的首要工作就是寻找重积分的计算方法，给出覆盖深度为时的精确计算结果。2、深入挖掘与之间的关系，希望能够找出值一定的情况下，二者之间的内在联系。3、对模型的结果进行实际的应用，调用真实数据来检验模型，并对模型进行改进。参考文献1基因组学与应用生物学，2009年，第28卷，第5期2 曾培龙基于引导的基因组序列拼接；分类号：T

27、P3919/2012年12月25日。3 徐静怡基因组序列拼接算法及新基因的发现2004年6月30日4基因组组装新技术研讨,访问日期:2014年5月11日。附录附录一%提取碱基序列jianji2=;for i=2:4:length(mc2) jianji2=jianji2;cell2mat(mc2(i);end%提取质量数据zlszf2=;for i=4:4:length(mc2) zlszf2=zlszf2;cell2mat(mc2(i);end%作碱基质量统计图zlszf=;for i=4:4:length(mc1) zlszf=zlszf;cell2mat(mc1(i);endzls=ab

28、s(zlszf);m,n=size(zls);zls=reshape(zls,1,m*n);hist(zls);%删除质量数值不合格的碱基序列m,n=size(zlszf2);jianjisx2=;for i=1:mif sum(abs(zlszf2(i,:)-64)>=20)>=88jianjisx2=jianjisx2;jianji2(i,:);endendxlswrite('E:jianjisx2.xlsx',jianjisx2);%作改后碱基质量统计图zlszf=;for i=4:4:length(mc1) zlszf=zlszf;cell2mat(mc1(

29、i);endzls=abs(zlszf);m,n=size(zls);zls=reshape(zls,1,m*n);hist(zls);附录二%按照碱基互补配对原则转换的碱基序列为m,n=size(jianjimat2);for i=1:m for j=1:n if abs(jianjimat2(i,j)=65 jianjimat2(i,j)=jianjimat2(i,j)+19; elseif abs(jianjimat2(i,j)=71 jianjimat2(i,j)=jianjimat2(i,j)-4; elseif abs(jianjimat2(i,j)=67 jianjimat2(i

30、,j)=jianjimat2(i,j)+4; elseif abs(jianjimat2(i,j)=84 jianjimat2(i,j)=jianjimat2(i,j)-19; end endend%删除重复的碱基序列m,n=size(jianjimat1);num=;for i=1:m-1 for j=i+1:m if jianjimat1(i,1)=jianjimat1(j,1) continue elseif jianjimat1(i,2)=jianjimat1(j,2) continue end if sum(abs(jianjimat1(i,:)=jianjimat1(j,:)=88

31、; num=num,j; end endendnum=num'jianjizz=;for i=1:length(jianjimat1) if any(num=i)-1 jianjizz=jianjizz;jianjimat1(i,:); endend附录三% 求那些数据相同的碱基对，其中：表示有个连续碱基相同、 表示第个与第个数据相同、 假定、计算相同项的位置clc fid=fopen('jianjizz11.txt','r'); data=fread(fid);m,n=size(data);data=reshape(data,90,5091);data

32、=data'K=60; m,n=size(data); n=n-2; d1=; d2=; d3=; for j=20:K-1 for i=1:m for k=i+1:m index1=; for t=j+1:n index=abs(data(i,t)-data(k,t-j); index1=index1,index; end index2=sum(index1); if index2=0 d1=d1,j; d2=d2,i; d3=d3,k; end end end end d=d1;d2;d3;c=d'xlswrite('E:dnanew1.xlsx',c); fclose(fid);% 计算最佳路径clc c=xlsread('E:dnanew1.xlsx'); c=c' j=c(1,:); i=c(2,:); k=c(3,:); n=length(j); a=; b0=0; a0=; for t=1:n a=i(t); b=88; for p=t:n if k(t)=i(p) k(t)=k(p); i(t)=i(p); a=a,i(p); b=b+j(p); end end if b0<=b b0=b; a0=a; end end a0;b0; l=a0' xlswrite('