数学建模“基因重组问题”

基因组组装摘要快速和准确地获取生物体的遗传信息获得目标生物基因组的序列信息，进而比较全面地揭示基因组的复杂性和多样性，是生命科学领域的重要研究内容。本文主要研究的是利用一定的方法将测序得到的短片段序列组装成更长的序列。本文根据目前的测序新策略，首先按照策略的思想，确定主要算法，再根据所确定的算法编写对应的程序，最后导入数据，得出最长的组装序列。对于问题一：本文采取分部解决问题的方法，具体步骤如下：步骤一：我们运用Matlab软件，对数据进行了预处理，将读长序列完整的提取出来；步骤二：我们采取了基因组组装算法优化模型，使得碱基数目尽可能大及组装序列的总长度的比例尽可能大；步骤三：在步骤二的基础上，我们建立了deBruijn图方法的模型，该模型能较好地解决测序中可能出现的个别碱基对识别错误、基因组中存在重复片段等复杂情况，形象地描述了我们前面所建立的模型。对于问题二：这个问题刚好成为问题一我们所建立的模型与编写的程序的检验，根据问题一所建立的模型和编写的程序，我们将数据导入程序之中，利用matlab软件运行并得出结果，最后在所得结果进行连续性、完整性、准确性的检验。关键词：基因组组装，Matlab编程，deBruijn图，读长（read）序列一、问题重述快速和准确地获取生物体的遗传信息对于生命科学研究具有重要的意义。获得目标生物基因组的序列信息，进而比较全面地揭示基因组的复杂性和多样性，成为生命科学领域的重要研究内容。确定基因组碱基对序列的过程称为测序(sequencing)。当然，由于技术的限制和实际情况的复杂性，最终组装得到的序列与真实基因组序列之间仍可能存在差异，甚至只能得到若干条无法进一步连接起来的序列。对组装效果的评价主要依据组装序列的连续性、完整性和准确性。连续性要求组装得到的(多条)序列长度尽可能长；完整性要求组装序列的总长度占基因组序列长度的比例尽可能大；准确性要求组装序列与真实序列尽可能符合。基因组组装软件可根据得到的所有读长组装成基因组，这些软件的核心是某个组装算法。常用的组装算法主要基于OLC(Overlap/Layout/Consensus)方法、贪婪图方法、deBruijn图方法等。问题一：建立数学模型，设计算法并编制程序，将读长序列组装成基因组。算法和程序应能较好地解决测序中可能出现的个别碱基对识别错误、基因组中存在重复片段等复杂情况。问题二：现有一个全长约为120,000个碱基对的细菌人工染色体(BAC)，采用Hiseq2000测序仪进行测序，测序策略以及数据格式的简要说明见附录一和附录二，测得的读长数据见附录三，测序深度(sequencingdepth)约为70×，即基因组每个位置平均被测到约70次。利用你的算法和程序进行组装，并使之具有良好的组装效果。二、问题分析问题一分析：首先，我们把问题一细分成几个小问题，需要我们逐步攻破。(1)根据测序策略，基于读长(reads)建立数学模型，并设计相关算法和程序将读长序列组装成基因组；(2)尽量解决组装过程可能出现的若干问题，比如个别碱基对错误、基因组存在重复片段等；(3)建立的模型要尽可能符合连续性要求、完整性要求、准确性要求；问题一中我们根据最新的测序策略，确定出自己的算法，根据我们确定下来的算法编写相应程序，具体算法步骤如下：步骤一：将原始数据构造成矩阵D；步骤二：计算矩阵维度m,n（m表示矩阵宽度，n表示矩阵长度）；步骤三：赋值，j=1,i=1,k=2。其中j表示忽略的字符数，i表示第i组数据，k表示第k组数据；步骤四：将矩阵中两组数据作差（一组数据从第j项到第n项，另一组数据从第1项到第，并且j依次递增）。判断其值是否为0，若真，则记录此时的i，j，k。若假，则继续步骤五；步骤五：1kk，mk?若真，则返回步骤四，若假，则继续步骤六；步骤六：mi若真，则1ii，2ik。返回步骤四，若假，则结束。问题二分析：问题二是在问题一的基础上，根据测序获取的两组读长序列，基于上述建立的的算法和程序，组装一个全长约为120,000个碱基对的细菌人工染色体(BAC)。首先我们对数据进行预处理，用Matlab软件，从两组reads(读长)筛选出所需碱基序列，我们只需由所给数据中提取出其中的碱基序列，即按从第二行起，将每隔四行的数据提取，保存为新的文本以供后续程序运行调用。结合问题二的已知数据，在MATLAB软件中输入问题一所编写的相关模型求解代码，最后对所得结果进行连续性、完整性和准确性的分析，验证我们的算法具有良好的组装效果。三、问题假设1.假设所提供的数据与真实数据误差相对较小；2.假设测序及组装过程操作中碱基的属性和对应原则不受影响；3.读长序列匹配按照首尾衔接原则，且碱基对质量的筛选问题不考虑；4.假设衔接过程中碱基重复个数低于4个均为巧合，超过4个读长序列按原则连接。四、符号说明符号说明ijH表示第i个读长第j个碱基sB序列匹配中第s个碱基iA表示第i个读长碱基A的个数iC表示第i个读长碱基C的个数iG表示第i个读长碱基G的个数iT表示第i个读长碱基T的个数Z表示最终某一条组装序列的长度Y表示匹配序列中的碱基个数X表示匹配组装过后的最后长度K表示最少实行匹配的相同个数M表示最多实行匹配的相同个数R表示读长的碱基个数（常数88）L表示序列匹配中与读长匹配的碱基个数五、模型建立及求解5.1问题一：5.1.1模型一：基因组组装算法优化模型根据已知模型组装要求，即要求碱基数目尽可能大及组装序列的总长度的比例尽可能大，可建立目标函数：)(XMaxMaxZ(1)因为Z是最后得到的所有碱基序列中的部分碱基序列X组装之后的最长碱基序列，此时Z定然满足与其它的读长序列不能再匹配的条件。Y为有重叠部分的碱基链，R为每一个读长碱基个数，L为相匹配的序列中相同的碱基个数，因此，我们易得X约束条件为：LRYX(2)查阅资料并进行相关推导得K(最少实行匹配的相同个数)取4，由于事先处理了全部相同的读长，故M取87。读长记为reads1，reads2，并且将reads2中的读长序列全部转换(碱基互补配对原则)进reads1，同时我们得到其它相关的约束条件：RGTCAiiii(3)KLM(4)11LYiBH(5)22LYiBH(6)⋮YiLBH(7)在以上条件确立的情况下，我们根据测序策略确定了如下关于所建立模型的核心算法：步骤一：将原始数据构造成矩阵D；步骤二：计算矩阵维度m,n（m表示矩阵宽度，n表示矩阵长度）；步骤三：赋值，j=1,i=1,k=2。其中j表示忽略的字符数，i表示第i组数据，k表示第k组数据；步骤四：将矩阵中两组数据作差（一组数据从第j项到第n项，另一组数据从第1项到第，并且j依次递增）。判断其值是否为0，若真，则记录此时的i，j，k。若假，则继续步骤五；步骤五：1kk，mk?若真，则返回步骤四，若假，则继续步骤六；步骤六：mi若真，则1ii，2ik。返回步骤四，若假，则结束。为了方便理解我们所建立的模型的算法与相应程序，我们做出了程序框图如下（具体程序见附录程序二）：Yj,i,kP=0P=i-k开始构造原始数据矩阵D矩阵维度m,nj=1,i=1,k=2k=k+1mkmi结束Yi=i+1k=i+2NYNN5.1.2建立deBruijn图方法模型基于deBruijn图数据结构的read之间对比拼接算法可概括一下几个步骤：(1)把筛选过得序列集合S作为参与比对的read库；(2)确定k值，建立debruijn图。这是需要扫描所有read数据，将每一个长为L的read拆分成L-k+1个k-mer,并用所有read的所有k-mer来累加，建立节点和边都加权的debruijn图；(3)化简debruijn图，连续线性延伸节点合并成为单一节点，产生一些碱基序列更长的节点；(4)遍历debruijn图产生contig。以上是对基于debruijn图的算法做了简单介绍。debruijn图是一种特殊的加权图，不仅图的及节点上有权值，而且边上也有权值。使用debruijn图只能将较短的read拼接成较长的contig，要得到全基因组，还需要contig的组装过程。contig的组装，即可以将read转化成定长的k-mer，并将这些k-mer存入bebruijn图中，以备之后查找使用。此时要设定的一个重要参数是k-mer的长度。选定k值之后，要将长度L的read拆成L-K+1个k-mer。根据一定策略，选定一个初始k-mer，接下来就可以在该k-mer为结点开始搜索后继的k-mer。搜索时采用贪婪图策略，每一步选择在当时看来最优的后继k-mer，直到满足事先设定的终止条件，结束一条contig的拼接，接着开始下一条contig的拼接。直到没有合适的初始k-mer可供选择，整个拼接过程结束。deBruijn图的拼接算法的核心是K值的选择方法，本文假设出现4个碱基相同的视为巧合情况。故将K=4进行研究，这样使得deBruijn图有相对好的连通性。下面简单用一个例子表达以上模型的含义:将基因组复制若干份，无规律地分断成短片段后进行测序，然后寻找测得的不同短片段序列之间的重合部分，并利用这些信息进行组装。例如，若有三个短片段序列分别为：ATACCTTGCTAGCGTGCTAGCGTAGGTCTGAGCTTAACTGCTTAACTGTTACCAGATC则有可能基因组序列中包含有ATACCTTGCTAGCGTAGGTCTGAGCTTAACTGTTACCAGATC这一段。我们的程序就是按此方法，一步一步找出最长的组合完成之后的序列。5.2问题二：对于问题二，我们在问题一的基础上，将所给的读长数据进行提取，导入matlab软件中，结合问题一解答所编写之程序，运行后得出结果。现在我们就可以用问题一的程序来解决问题二，并让显示的结果来验证问题一所建立的模型的组装序列的连续性、完整性和准确性。首先我们利用编写的程序（见附录程序一）将原始数据进行提取，得出一个只包含碱基序列的文本，并将其导入matlab软件，结合问题一模型所编写的程序，得出问题二所要求的组装后的DNA序列。（由于数据量较大，假若将完整DNA打乱之后的全部片段导入并运行程序，耗时相当长，而且后续结果也难以处理和展示，所以在此我们只提取碱基序列文本的前10000行碱基序列进行组装），组装结果如下：（因为基因的起始序列不同，所以程序运行完毕会得出多个结果，在此我们选取较长的几个结果序列表示出来）结果一：CGGTTAGGGTTAGCAATAAAACAAAAATCTAAAAAGAAAAGAAAAATTAGGGTTAATCAAGACCTAGATTAGATCGATTTTTTCTTTTGACGAGTCTCGGTAAGACCATGAGAAGCTTGTCTTGTTTTGACAACATCTCATTTATAATATTTTGGTAGGACTAAACTTTCCCGCGCTCGAAAACATGACCGGATGCTCGAGTTTGAGCTCCTCCGGCGAAGCACCGCCGGGTGCCGCCGTCTCCGCCGCGCGCCGCACAACACATCGGATACATCCCGAAAGGGTCGTCATCGCCAGGGTCTCCGCCGTCATCACCATCTTCATTGTCGTCGCTGCCGTCTGGCAGCGCTGTCGCTGGCGATTGTCATCAGTGCGCTGCTGGTGAACGCTGTGGGGATTATCGGCTTTATCGGGTTGTTCGCGC结果二：AAATATGACTCGTGCAGAAGCATGGGCTGATAGACTAGTTTTCACAAATTTAGATTCATAGTCGTTATTATATATCTACTAGACACAACGTTTACATGCCCTTATGGAACATTTCCATATAGAAGGATGCCCTTCGGTCTATGCAACGCTCCAGCCTCTTTTCAACGTTGCAATACACCACTTTGATTTAATAGTGGGATTGATATTCTTTCATG