暑期数模培训作业

基因组组装摘要快速和准确地或其生物提的遗传信息对生命科学研究具有重要的意义。测序技术从第一代到现在普遍应用的第二代以及正在兴起的第三代，能直接读取的碱基对序列长度远小于基因组长度。所以测序之前DNA分子要经过复制若干份、随机打断成短片段。要获得整个DNA片段，需要把这些片段利用重合部分信息组装连接。如何在保证组装序列的连续性、完整性和准确性的同时设计耗时短、内存小的组装算法是本题的关键。对于第一问，我们根据题目所给的一些方法再结合了相应的资料，最终我们依据基于Hamilton路径的拼接算法利用用优化后的Phrap算法并根据Phrap算法的内容选取了Smith-Waterman算法以及Tarjan算法等算法对问题建立了数学模型。对于第二问，我们根据第一问所建立的模型根据算法过程对题目所给的数据进行了预处理，然后编写了Java程序，把处理后的数据输入程序得到相应结果。关键词Hamilton图；Phrap算法；DNA拼接一、问题重述快速和准确地获取生物体的遗传信息对于生命科学研究具有重要的意义。对每个生物体来说，基因组包含了整个生物体的遗传信息，这些信息通常由组成基因组的DNA或RNA分子中碱基对的排列顺序所决定。获得目标生物基因组的序列信息，进而比较全面地揭示基因组的复杂性和多样性，成为生命科学领域的重要研究内容。测序技术始于20世纪70年代，伴随着人类基因组计划的实施而突飞猛进。从第一代到现在普遍应用的第二代，以及近年来正在兴起的第三代，测序技术正向着高通量、低成本的方向发展。利用现有的测序技术，可按一定的测序策略获得长度约为50100个碱基对的序列，称为读长（reads）。基因组复制份数约为50100。基因组组装软件可根据得到的所有读长组装成基因组，这些软件的核心是某个组装算法。常用的组装算法主要基于OLC（Overlap/Layout/Consensus）方法、贪婪图方法、de Bruijn图方法等。一个好的算法应具备组装效果好、时间短、内存小等特点。新一代测序技术在高通量、低成本的同时也带来了错误率略有增加、读长较短等缺点，现有算法的性能还有较大的改善空间。问题一：试建立数学模型，设计算法并编制程序，将读长序列组装成基因组。你的算法和程序应能较好地解决测序中可能出现的个别碱基对识别错误、基因组中存在重复片段等复杂情况。问题二：现有一个全长约为120,000个碱基对的细菌人工染色体（BAC）， 采用Hiseq2000测序仪进行测序，测序策略以及数据格式的简要说明见附录一和附录二，测得的读长数据见附录三，测序深度（sequencing depth）约为70，即基因组每个位置平均被测到约70次。试利用你的算法和程序进行组装，并使之具有良好的组装效果。二、问题分析伴随着人类基因组计划的实施和突飞猛进，基因组测序组装的研究意义越来越重要。如果不能快速、准确的获取信息，基因组的研究将会是举步维艰，人类基因组计划也难以进行下去。所以，建立合理有效的基因组组装模型迫在眉睫。对于问题一，首先，我们通过DNA (RNA)分子中碱基对的互补特征，确定基因拼接方式，即将DNA分子中的一条链的互补链与其对应的另一条链比对，找到阈值允许范围内个数的相同碱基序列将单个基因片段记录并相连。其次，根据Hamiltom图论原理，将拼接上的特征子串看做是一个个结点，并将这些结点根据其中信息的联系通过有向线段连接，然后在这些所有的有向线段中找到一条最长、通过结点最多的路径，次路径就是我们要找到的基因链。最后，基于对拼接中的 de Bruijn图结构的研究，为使repeat结构拼接准确性更高，讨论的 de Bruijn图中错位数对拼接结果的影响。将之前的模型进一步优化三、模型的假设与符号的说明四、模型的建立与求解对于问题一的求解 本题是基于新一代测序技术的基因组装算法问题，要求设计算法针对性的解决新一代测序技术带来的一些弊端。1.先是通过DNA (RNA)分子中碱基对的互补特征，确定基因拼接方式DNA（RNA）分子中碱基对的互补特征1：在DNA分子结构中，由于碱基之间的氢键具有固定的数目和DNA两条链之间的距离保持不变，使得碱基配对必须遵循一定的规律，这就是Adenine（A，腺嘌呤）一定与Thymine（T，胸腺嘧啶）配对，Guanine（G，鸟嘌呤）一定与Cytosine（C，胞嘧啶）配对，反之亦然。碱基间的这种一一对应的关系叫做碱基互补配对原则。 腺嘌呤与胸腺嘧啶之间有两个氢键，鸟嘌呤与胞嘧啶之间有三个氢键，即AT, GC。基因拼接方式2：由于碱基间的这种一一对应的关系 碱基互补配对原则，而当前测序仪的准度不能 超过500bp(base pair，碱基对)的序列具有很高的准确性。所以对于长的基因（如人类基因组DNA长达30亿bp），在生物学上普遍采用逆向构造的方法，将较长的序列做多个克隆，用酶将这些克隆随机割成小片段，然后又用测序反应测出小片段上各个位置上的碱基，直观的读出A,C,G,T的序列。通过一定的处理将互相重叠的序列拼接起来，还原出片段克隆的完整序列。而当前广泛采用的拼接方法是Celer公司的鸟枪法。基因组DNA 随机打碎克隆序列序列拼接图1 全基因组鸟枪图法示意图全基因鸟枪法的主要工作包括两部分：shotgun 测序和序列拼接。其中shotgun测序主要流程如下，用途表示如图2.1. DNA的提取和纯化。2. 载体预备：和DNA片段结合，从而能够在细菌中扩增。3. DNA片段的制备：将DNA用超声波切成能够测序的小片段。4. 转化培养：小片段和载体结合，植入细菌中进行扩增。5. 提质粒：从细菌中提取出繁殖好的质粒。6. 电泳检测：检测质量的好坏。7. 测序：用测序仪侧序。DNA整体切成小段小段和载体结合结合进行测序图2 Shotgun 测序流程示意图全基因测序大致包括以下几个步骤。序列测定：使用Sanger的双脱氧链终止法，即通过化学作用是DNA链在每一个碱基处终止，并在终止处染色（是终端带放射性物质或荧光物质），通过变性凝胶的电泳和激光照射可以得到终止于不同碱基的DNA 片段的谱带。在激光照射到不同的碱基对上时，会反射不同的颜色，所以得到的谱带一般有四条，它记录的一般是不同颜色的放射光的亮度，测定美国片段的序列。碱基读出：运用计算机独处软件进行图像处理，以确定谱带所蕴含的DNA 序列。四条谱带在测序时被同步保留，处理软件试图按照一定的间隔读出一个碱基。序列拼接：在得到每个片段的序列后，利用这些片段见的重叠，将它们的拼接成原来的序列，或尽可能接近原来序列的几个较长的序列。因为片段数据中含有序列的夺个克隆且片段的打断都是随机而独立的，这给拼接提供了足够的可能性。拼接的关键问题是得到每个片段在一个长序列中的位置信息，这种组合的集合称为contig (contiggoussegment).序列拼接是测序工作的关键所在，序列拼接方法和实现技术的好坏，直接影响到拼接结果的优劣，同时它与其他领域的联系广泛，也成为各学科研究的热点问题。4.1一问模型Hamilton路径基于Hamilton路径的拼接算法总共有三个过程：找出序列片段间的的重复信息将存在有重合的片段组合起来形成一个contig结构根据片段中每个碱基的质量值在contig寻找中寻找一条“Consensus”最终序列根据上面的三个步骤，我们打算运用Phrap算法来对这个问题进行求解，Phrap算法大致也分为三步：QverlapLayoutConsensus1) 寻找Overlap区域使用实验所得到的基因片段reads作为数据源，片段间两两对比，找到所有的重叠区域，并对其进行筛选得到的所有的符合条件的重叠区域。2) Layout阶段利用在Overlap阶段得到的所有重叠区域，建立序列间的组合关系，称为“contig”。3) Consensus阶段使用Layout阶段得到的所有contig结构，在contig结构的基础上建立加权有向图，遍历图得到权重最重的路径，获取与路径相对应的的序列称为“Consensus”。一、 获得片段间的重要信息（Overlap）Phrap算法假定两个片段之间如果发生拼接，则这两个片段之间必定有一部分是重叠部分。在它们重叠的部分上必然有一个最小长度为L的精确匹配区域。Phrap算法中将这个最小长度为L的子序列称为字（Word），没一个字的长度为。在Phrap算法中，当一个Word的长度小于min-word时，该算法认为这两个片段不能够进行拼接，而当一个word的长度大于max-word时，Phrap算法认为这两条片段完全相同。即，Phrap算法设置了一个字长的区间用于对重叠区域进行筛选。对所有的reads进行两两对比，在碱基序列集合F中随机找出一个没有被处理过的碱基序列fi，在中寻找所有可能的重叠区域并对区域进行筛选输出符合条件的重叠区域Overlap，直到F中所有序列处理完毕。重叠区域判定条件计算法：1）判定条件满足条件的重叠区域就叫Overlap,其意义在于如果，则表示这两个序列由于重叠区域太短不能进行拼接。如果则表示重叠区域过长，此时这两个序列视为完全相同的重复序列。由于序列的测量存在误差，因此在进行重叠区域的寻找时应才用模糊匹配，在score(fi,fj)大于阈值的情况下才视两个模糊序列匹配成功，否则认为这两个序列的重叠区域误差过大，不视为一个Overlap。Phrap算法描述为Input: n reads sequence,F=f1,f2,fnOutput: all the Overlaps between any two readsBegin1. Let the minlen,maxlen, be the minimum length,maximum length and a threshold of Overlap respectively;2. Trim of the low qualitys reads;3. Use the Smith-Waterman algorithm to calculate the score of any two reads fi and fj (fiF,fjF,ij and I,j1,n)4. If score (fi,fj and then compute the Largest Likely Ratio (LLR) of fi and fj,and output the Overlap(fi,fj).Else repeat the step 3;Iterate the step3 and 4 until each read of F has been handled.Phrap算法中运用Smith-Waterman算法来精确计算两条片段之间可能存在的重叠并对每一个可能的重叠进行打分，同时，设立一个阈值参数W，当重叠的分值大于阈值参数时我们就认为该重叠为一个Word。输入：具有n个顶点的有向图G，指定顶点v和v 第一步：在G中随机生成大量的各种各样的路径 第二步：去掉所有不以顶点v为起点和不以顶点v为终点的路 第三步：去掉所有没有恰含n个顶点的路 第四步：对于这n个顶点中的每一个顶点v，去掉所有不包含顶点v的路 设G=（V,E,W）其中V,E,W分别图G的顶点，边，边上权值的集合。设图G是有P个顶点的图，即V(G)=(V0，V1，V2Vp-1)；如果存在以V0为始点以VP-1为终点一条包含图G每个顶点一次且仅一次的路P,则称P是V0为始点VP-1为终点的一条Hamilton路，显然当图G是赋权V图且存在以V0为始点VP-1为终点一条包含G每个顶点一次且仅一次的路径P,则必存在一条权值最低的Hamilton路。这条权值最低的路径也就是我们所要找的。我们通过此方法即可得出读长序列组装成基因组。Hamiltom图对于基因的重复问题我们用全局比对，全局对比就是对两条序列进行编辑操作，通过字符匹配（或错配）、插入及删除使得序列达到一样长度并使得相似性得分最大的过程。相似性得分的定义如下：设两条序列S1和S2进行全局比对，其中=m,=n,比对后的序列和,长度为，序列和经过编辑操作得到，每种编辑操作都赋予一定的得分，比对后所得到的序列和的相似性得分定义为所有操作得分总和（此处假设所有操作都是独立的）。匹配（或错配）操作定义为两序列中定义的字符相等或不等，相等时赋予一个正的得分m，若不等，赋予一个负的得分d，用函数 插入操作就是在中插入空位，此时得分按照空位罚分的规则计算;删除操作就是在中删除字符，这等价于在中的相应的位置插入空位，所以删除操作的得分也按照空位罚分的规则计算。空位罚分是为了避免进行过多插入和删除操作使得比对失去意义。空位罚分通常有两种罚分规则：权值恒定模型仿射空位罚分模型。权值恒定模型中，每个空位的罚分相同，都是-p，因此k个空位的罚分同空位的长度k成比例penalty=-pk。而仿射空位模型中，罚分包括两部分，开始一个空位的罚分和按照空位长度的罚分。罚分公式：penalty= 其中q为开始一个空位的罚分，k为空位的长度，r为每一个空位的罚分。由于放射空位罚分模型中引入一个开始罚分q,并且q的取值通常比r大的多，这样避免了轻易引入空位，也既避免了过多的插入删除操作，另外，一旦引入一个空位后，紧接着出现几个连续空位的可能性要比出现几个单一空位的可能性大得多。可见放射空位罚分模型虽然比权值恒定模型复杂一些，但它更合理更具有生物学意义。计算全局比对有许多种方法，其中最经典的就是Needleman-Wunsch【3】算法。它实际上使用了动态规划的思想。动态规划就是把较大的问题按照时间步骤或空间分布分割处理，从局部最优逐步寻找全局最优结果的一种方法【4】。下面将要介绍的动态规划全局序列比对算法就是在Needleman-Wunsch算法的基础上得来的。列罚分ISIIIDSSSDDD行行详细图解列罚分罚分罚分罚分罚分罚分罚分注：图中上部分是矩阵的全局示意图，下部分是邻近节点的详细图解。每一个节点都有三种可能的状态：插入、删除和匹配（错配），分别用I、D、S来表示。图中各个节点到相邻节点的罚分均已标注。按箭头所指方向可以得到矩阵的计算方法，反之，若按箭头所指的方向就可以的到最佳对比。Smith-Waterman算法：设两个序列，。和（其中，）皆属于某个字符集。对中的任何元素和空符号，它们两两之间进行都有一个计分值，用计分函数表示。表示序列的前缀与序列的前缀之间最优相似性比较的得分。然后可得如下式子：根据上面的式子，我们可以得到一个计分矩阵通过得分矩阵进行动态规划回溯获得相似性比较的算法如下：其中函数表示在一个序列b的第c个位置插入一个字符a。二、 建立片段间的组合关系（Layout）我们通过第一阶段（Overlap）得到了相应的Overlap和LLR分值，在Layout阶段对所有的Overlap区域所在的reads进行组合，形成contig结构，即用contig结构记录reads之间的组合关系。Contig结构包含F，S,Offset :reads集合，S：F中的reads表达的模糊序列，Offest：F中的reads相对S的偏移量。然后对所有的Overlap按LLR值降序排序，并按照此顺序依次处理所有的Overlap。这样做的好处是可以在一定程度上解决重复序列repeat对拼接的影响重复序列repeat是指生物基因组中的重复片段，有两种类型：一、 简单repeat，基因序列中存在单一的repeat区域RR二、 复杂repeat，除了单一repeat区域外，在repeat区间内还有子重复区间。R1R2R3R1R2R3由于上面的原因，可能导致拼接错误，例如：s1RRs2原始DNA序列s1Rs2错误的拼接Phrap算法对Overlap按LLR得分进行排序，按排序后的顺序进行处理，其本质还是使用“贪心算法”对DNA片段进行处理，这种处理方法在一定程度上减少了重复序列对拼接的影响，但对于重复序列过长的，或者repeat区域过多的基因则无法避免rpeat对拼接的影响。Layout阶段的目的是对所有的reads进行组合，形成contig结构。首先要对每个reads都初始化为contig结构，之后按照Overlap排序的顺序一次处理每个Overlap。在处理Overlap（fi,fj）时的步骤：1. 找到fi,fj所属的contig(i)和contig(j)。2. 如果contig(i)和contig(j)不是同一个contig，则对两个contig进行合并。3. 如果contig(i)和contig(j)是同一个contig，则不处理Overlap。4. 按照排序处理下一个Overlap，直到处理完毕。在这些步骤中，合并contig(i)与contig(j)为contig（i,j）的步骤为：1. 更新F集合为FFj，将Si和Sj合并为Sij2. 重新计算F合集中所有的reads相对于Sij的偏移量Offset，同时删除contig(i)与contig(j)。此时，Fi和Fj中的reads所属的contig就变成了contig(i,j)。这样Phrap算法将所有存在Overlap的reads组合起来形成了一个contig结构。由于数据源reads存在实验误差，所以Layout阶段生成一个contig集合最好的情况是集合中只含有唯一一个contig。算法描述：Input: Input：Overlapset，RLLset，nreadssequencesF=fl，21，fn1Output：contigsetBegin1. Create the contig structure(contig(fi) for each read in F,fiF,i1,n;2. Sort the Overlap by the LLR score in decreasing order;3. Under the LLR decreasing order,if the Overlap (fi,fj) is unhandled,then merge the contig(fi) and contig(fj) into contig(fi,fj);4. Iterate the step 3 until all the Overlap have been handled.三、形成Consensus序列（Consensus）在Layout阶段对每一个contig建立一个模糊序列S并计算contig中的所有reads对于S的相对偏移量Offset，是为构造每个contig的加权有向图做准备。由于contig结构中S是使用模糊序列拼接而成（经由Smith-Waterman序列匹配算法得到）故形成的S是存在很大误差的。正是由于S存在很大误差，因此，S并不作为曾哥cintig拼接结果的输出而是在reads之间建立加权有向图考虑每个碱基的碱基品质值，输出contig中由高品质的序列片段拼接成的结果。在contig结构的基础上，Phrap使用“Mosaic”方法，将每一个contig结构中高品质的序列片段连接起来形成一条Consensus序列作为此contig最终的拼接结果。对于每个contig，在此结构的基础上构成一个加权有向图G，求解Consensus序列转换为在图G中寻找一条从起始点到终点权重最大的路径。图G的构造方法如下1. 节点。在一个contig的序列集合F（F=f1,f2,f3）中，将偏移量Offset最小的reads的最左端的K个碱基定义为起始节点（K=1），偏移量Offset的最右端K个碱基定义为终止节点。2. 边。同一个read中两个相邻的节点之间定义一条单向边，其权重为两个节点间所有的碱基的品质值之和；在两个reads的Overlap区域对应的节点间定义一条权重为0的双向边。Contig的有向加权图Phrap算法利用contig结构，在一个contig中的所有reads之间建立了一个加权有向图，由两个不同的contig(i)和contig(j)建立的加权有向图Gi与Gj是完全相互独立的。Phrap算法使用Tarjan算法来寻找最大的路径。4.2二问模型4.2.1. 问题分析本题题目提供全长约为120,000个碱基对的细菌人工染色体，采用新一代的Hiseq2000测序仪进行测序。附件提供了筛选好的定长read数据文件。使用第一题设计的基于de bruijn图的组装算法对read数据进行组装，并对结果进行误差分析。4.2.2.数据分析由测序策略可知，read1和read2为相互补的两条单链，故选用read1的数据带入算法进行组装，read2可以作为校准链备用。分析read1文件可知，本题数据已满足如下条件：（1）序列片段被切成固定长度L=88bp；（2）经过复制，原基因至少有6个副本；（3）所有片段上的碱基都已经被识别出来，不存在未知碱基；（4）由于技术限制，本文不对质量数进行讨论，假设read1中的所有片段满足正确率要求。4.2.3.带入模型求解4.2.4建立de bruijn图（1）将k值定为4。把上述read1文件中的序列存入库中，开始建立read条目的数据结构和k-mer条目的数据结构。预读数据，逐条读取read数据，每条readID进行升序保存；生成该read上所有k-mer（共85个k-mer），统计这些k-mer出现的次数，填写k-mer结构中的num字段。如图所示，为相关代码片段。相关数据录入程序源代码见附录。（2）遍历de bruijn图，根据上一步统计的k-mer数量，申请readID_pos数组所需要的内存空间。依次读取每一个read，填写readID_pos数组中的第cur行，填好之后把cur值加1。（3）将碱基替换成2位二进制数。A=00，C=01，G=10，T=11。4.2.5模型求解由于数据非常庞大，演算拼接过程不能完整的展示，接下来将列举一段算法拼接的过程：（1）初始k-mer定为 AAAC即 00000001，该k-mer出现在4条read上，且k-mer出现在每条read上的pos为1.这四条read开始参与拼接。如图为k-mer比对拼接相关代码（2）此时num 为4，addr为1，cur为-1；read100000001011111010100001011101111000100001001110100001000read200000001000011100000101010010000011000111000110011010100read300000001110001001010101110111111110101101101000010010010read400000001111110001110101000000011010001011111111111000010初始kmer 00000001后继kmercontig 00000001（3）当前候选后继k-mer情况如下图：初始kmer 00 00 00 01候选后继kmer100 00 01 01候选后继kmer200 00 01 00候选后继kmer300 00 01 11候选后继kmer400 00 01 11（4）选定后继k-mer为 AACT,即00000111；进行下一段拼接，此时num为5，前驱结点cur 为1，addr为2；此时contig增加了一个碱基；read100000001011111010100001011101111read200000001000011100000101010010000read300000001110001001010101110111111read400000001111110001110101000000011read500000100111111001111110000000010初始kmer00000001后继kmer00000111contig0000000111（5）重复（4）步骤，直到无法找到符合条件的后继k-mer则该条contig拼接结束；若该条contig的长度大于100bp则保存下来，否则删除；程序继续运行以上案列，得到一个长度为360bp的contig如下： 110111110010010100001000100011000111111000111010010110010011100011011010111011000000010000111000010111110100001011110100010000111011000011001101111000000111111111010000001010010100100010010011011001011001111011001110100010010110001001110010010110010100111110110101111101010001010001011100001110100000110111101011011011011011011011010101111000101111010111110001001000111101101001101000110010000000110010000101001010001011100010011000111011101100010000010011100110101111100110101110000010100100000101111011001100110110001010101010110010001101001111010111001011111011001101001111111111101111110100111011110111111111001100001111001110111111001101001001001101110010011100101111010011101101110011110111001111001101001101111011 （6）到此一条contig拼接结束，开始重复（1）步骤；最终得到若干条长度大于100bp 的contig。 101100010011100100101100101001111101101011111010100010100010111000011101000001101111010110110110110110110110101011110001011110101111100010010001111011010011010001100100000001100100001010010100010111000100110001110111011000100000100111001100 4.2.6 模型二总结通过基于de bruijn图的基因组装算法，成功拼接2条长度在100bp以上congtig；中途算法自动放弃长度不够的contig 364条。contig拼接成功率为0.54945 %。理想情况下，S中的read 也能完 整覆盖原基因组，即使是这样，也不能保证 g与原基因组完全一样，主要原因是基因组中存在大量重复片段。另外，测序过程中的错误也不容忽视。设ra，rb是S中的两条完全一样的read ，则有如下几种可能： （1）rb 是ra 的复本（或ra 是rb 的复本）； （2 ）rb 不是ra 的复本，他们是来自重复片段的两条read ； （3 ）rb 中包含测序错误，ra 中无错误，导致rb 碱基与ra 相同；（4 ）ra 中包含测序错误，rb 中无错误，导致rb 碱基与ra 相同； （5 ）ra和rb中都有错误，它们的碱基相同纯属巧合。所以，现在能做的是使碱基序列 g 尽可能与原始基因组一致，存在误差在所难免，本文要使该误差尽可能的小。模型评价模型的优点本模型针对新一代测序技术出现的问题逐一进行了解决。尽可能地拉长了对所能读取组装得到的碱基对序列长度，而组装序列的总长度站基因序列长度的比例尽可能的大，组装序列与真实序列尽可能符合。Needleman-Wunsch的优点：Needleman-Wunsch.算法定义为一个抽象类，该类继承在动态编程框架抽象类未进行优化，查找全局和局部比对的时间，这个算法是动态编程的应用。动态编程还可用于矩阵链相乘、装配线规划和计算机象棋程序。使用动态编程能够解决的一些问题，这些问题都有共同的特征：每个问题的解都能用递归关系表示;用递归方法对这些递归关系的直接实现会造成解决方案效率低下，因为其中包含了对子问题的多次计算;一个问题的最优解能够用原始问题的子问题的最优解构造得到.模型的缺点 由于技术的限制和实际情况的复杂性，最终组装得到的序列与真实基因组序列之间仍可能存在差异，甚至于还存在若干条无法进一步连接起来的序列。Hamilton图论缺点：无论是存储开销还是时间开销上都是非常大的。另外，Hamilton图论还需要保留大量临时和过渡数据；而且Hamilton图论由于其在进行片段比对前没有进行相应的处理，导致重复子序列信息进入了片段相互重叠阶段，干扰了拼接，使得它的结果受重复子序列的影响较大模型的改进在某些情况下采用合理的方法或技术，可以起到意想不到的效果。针对生物数量大、消耗内存较多的具体特性，我们可以使用到并行和分布式技术和Gamma模型来降低算法的运行时间和提高结果的精确度。关于Phrap算法也有许多改进的方法，例如其中包含的Smith-Waterman算发。介于时间问题，我们不再进行验算。参考文献1薛金星，高中生物必修2 遗传与进化，北京，人民教育出版社，2004年。2王磊，中南大学硕士学位论文，基因序列拼接算法的研究 ，2006年。3S.B.Needleman,G.E.Wunsch，A,general method applicade to the search for similarities in the amino acid sequences of two proteins. J.Mol.Biol.1970;4R.Bellman,S.E.Dreyfus,Applied Dynamic Programming, 1962.5胡嘉瑜，Phrap算法在DNA序列拼接中的应用与改进，现代计算机（ModemComputer）下半月版，2012（10），1.附录：源代码文件（1）DataInput.javapackage 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(B length:+b.length );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 sql = ;sql =insert into datasave(dataId,data,quality) values(?,?,?);PreparedStatement pstmt = conn.prepareStatement(sql, Statement.RETURN_GENERATED_KEYS);for(int i = 0;ia.length;i+)/sql = insert into datasave(dataId,data,quality) values(+ai+,+bi+,+ci+);pstmt.setString(1,ai);pstmt.setString(2,bi);pstmt.setString(3,ci);pstmt.executeUpdate();/stmt.executeUpdate(sql);/stmt.close();pstmt.close();conn.close(); catch (ClassNotFoundException e1) / TODO Auto-generated catch blocke1.printStackTrace();System.out.println(a1:+a1);public String StringToByte(String str)String reString=;char c = str.toCharArray();for(int i=0;ic.length;i+)if(ci=A)reString +=00;else if(ci=C)reString +=01;else if(ci=G)reString +=10;else if(ci=T)reString +=11;return reString;(2) DataToKmer.javapackage cn.data;import java.util.*;import java.sql.*;public class DataToKmer /* * param args */public static void main(String args) throws SQLException,ClassNotFoundException/ TODO Auto-generated method stub/*HashMap hs =new HashMap();*/* * 建立两个集合一个村存储read 一个存储Kmer * * */*String sql1 = ;*/String contig = ;readBank readBank =new readBank 46847;Read read = new Read46847;HashMap hsToRead = new HashMap();List bankList = new ArrayList();String getRead = ;/初始Kmertry Class.forName(com.mysql.jdbc.Driver);String url = jdbc:mysql:/localhost/genmdata?user=cyh&password=1011325;Connection conn =DriverManager.getConnection(url);ResultSet rs1;Statement stmt = conn.createStatement();PreparedStatement pstmt;String sql = select count(*) from datasave;int datalength = 0;rs1 = stmt.executeQuery(sql);if(rs1.next()datalength = Integer.parseInt(rs1.getString(count(*);System.out.println(datalength=+datalength);sql = select * f