2014年成都理工大学校内数学建模竞赛论文
二0一四年五月二十五日
摘要：本文所要研究的就是全基因组的从头测序的组装问题。
首先，本文简要介绍了测序技术及测序策略，认真分析了基因系列拼装所面临的主要挑战，比如reads数据海量、可能出现的个别碱基对识别错误、基因组中存在重复片段等复杂情况，探讨了当前基因组序列拼接所采用的主要策略，即OLC（Overlap/Layout/Consensus）方法、de Bruijn图方法，且深入探讨了de Bruijn图方法。
其次，针对题中问题，以一条reads为基本单位，分为reads拼接和contig组装两个阶段，其中contig是由reads拼接生成的长序列片段。Reads的拼接阶段主要包括数据预处理、de-Bruijn 图、contig构建等，而contig的组装阶段主要包括序列的相对位置的确定以及重叠部分overlap的检测，用序列比对的方法来提高拼接的精度。
最后，进行了算法的验证与性能的评价，并且针对问题2，进行了组装分析与验证，结果表明，得到的拼接基因组序列在小范围内与原基因组序列大致吻合。
关键词：基因组系列拼接； reads；de Bruijn图；contig组装；k-mer片段；
一.问题重述
基因组组装
快速和准确地获取生物体的遗传信息对于生命科学研究具有重要的意义。对每个生物体来说，基因组包含了整个生物体的遗传信息，这些信息通常由组成基因组的DNA或RNA分子中碱基对的排列顺序所决定。获得目标生物基因组的序列信息，进而比较全面地揭示基因组的复杂性和多样性，成为生命科学领域的重要研究内容。
确定基因组碱基对序列的过程称为测序（sequencing）。测序技术始于20世纪70年代，伴随着人类基因组计划的实施而突飞猛进。从第一代到现在普遍应用的第二代，以及近年来正在兴起的第三代，测序技术正向着高通量、低成本的方向发展。尽管如此，目前能直接读取的碱基对序列长度远小于基因组序列长度，因此需要利用一定的方法将测序得到的短片段序列组装成更长的序列。通常的做法是，将基因组复制若干份，无规律地分断成短片段后进行测序，然后寻找测得的不同短片段序列之间的重合部分，并利用这些信息进行组装。例如，若有两个短片段序列分别为
ATACCTT GCTAGCGT
GCTAGCGT AGGTCTGA
则有可能基因组序列中包含有ATACCTT GCTAGCGT AGGTCTGA这一段。当然，由于技术的限制和实际情况的复杂性，最终组装得到的序列与真实基因组序列之间仍可能存在差异，甚至只能得到若干条无法进一步连接起来的序列。对组装效果的评价主要依据组装序列的连续性、完整性和准确性。连续性要求组装得到的（多条）序列长度尽可能长；完整性要求组装序列的总长度占基因组序列长度的比例尽可能大；准确性要求组装序列与真实序列尽可能符合。
利用现有的测序技术，可按一定的测序策略获得长度约为50–100个碱基对的序列，称为读长（reads）。基因组复制份数约为50–100。基因组组装软件可根据得到的所有读长组装成基因组，这些软件的核心是某个组装算法。常用的组装算法主要基于OLC（Overlap/Layout/Consensus）方法、贪婪图方法、de Bruijn 图方法等。一个好的算法应具备组装效果好、时间短、内存小等特点。新一代测序技术在高通量、低成本的同时也带来了错误率略有增加、读长较短等缺点，现有算法的性能还有较大的改善空间。
问题一：试建立数学模型，设计算法并编制程序，将读长序列组装成基因组。你的算法和程序应能较好地解决测序中可能出现的个别碱基对识别错误、基因组中存在重复片段等复杂情况。
问题二：现有一个全长约为120,000个碱基对的细菌人工染色体（BAC），采用Hiseq2000测序仪进行测序，测序策略以及数据格式的简要说明见附录一和附录二，测得的读长数据见附录三，测序深度（sequencing depth）约为70×，即基因组每个位置平均被测到约70次。试利用你的算法和程序进行组装，并使之具有良好的组装效果。
附录一：测序策略
测序策略如下图所示。DNA分子由两条单链组成，在图中表现为两条平行直
线，两条直线上相对位置的两个碱基相互结合形成碱基对（bp），并且与碱基A 结合的碱基必为T，与碱基C结合的碱基必为G。将一个含120,000个bp的完整基因组，随机打断成500bp的片段，然后对500bp的片段进行测序。测序方法如第3步所示，分别从500bp片段的两端，对两条单链进行测序，测得的读长记为reads1，reads2。reads1，reads2的长度均为88bp，且该对reads相距500bp。
图1 测序策略示意图
附录二：数据格式
读长数据格式为fastq格式：
每4行表示一条reads
第一行：@序列ID，包含index序列及read1或read2标志；
第二行：碱基序列，大写“ACGTN”；
第三行：“+”，省略了序列ID；
第四行：质量值序列：字符的ASCII码值-64=质量值。
附录三：读长数据
测序得到的读长数据存放于两个fastq文件中（见附件一），其中McMc_BAC_1.fq.gz.clean.dup.clean和McMc_BAC_2.fq.gz.clean.dup.clean分别存放reads1和reads2的数据。
二.问题分析
正如上面问题所描述的一样，我们要解决的是要将基因小序列read组装成连续的基因大序列乃至最终的完整基因序列，而这就要将两个read1和read2
片段进行比较与拼接，比较的时候，因为相似片段的长短问题而不能确定拼接正确性，因此可以用两片段相似的权值来判断拼接的合理性，这样，若用点来代替read，用加权的边来判断到底要和哪个片段进行拼接，我们在查阅资料后，发现可以通过de bruijn图并对其进行相应的改进后来建立数学模型对问题进行求解。
设想一本杂志被复制成多份,将每份杂志均以不同的方式剪切,将多份剪切的杂志放在一起。在剪切的过程中,一些碎片丢失,一些碎片被污渍浸染,一些碎片存在着重叠现象。根据上述情况来寻找恢复原始杂志的方法。这是DNA序列拼接问题的现实模型描述。
基于de Bmijn图的序列拼接原理主要是通过构造并简化de Bmijn图结构来实现整个序列拼接的过程。
三.基于De Bruijn图的序列拼接技术分析与比较
二十世纪八十年代末，Pevzner等人提出基于de bruijn图的算法，并首次将该算法用于DNA序列拼接。基于de bruijn图的算法的核心思是将序列拼接问题转换为人们所熟悉的欧拉路径问题。Pevzner等人认为传统的
overlap-layout-consensus算法导致了将DNA序列拼接问题转换为Hamilton路径问题，他们受到杂交测序方法SBH(Sequencing by Hybridization)的启发，创造性地提出了在de Bruijn图中寻找欧拉路径的构想，尽管杂交测序方法SBH 从未在测序工程中实际应用过，但它直接引发了基因芯片工业的诞生。
构造de Bruijn图的方法如下所述：
（1）在read集合R={r1，r2，…，rn}中，首先将每一条read分割成若干k-mer(长度更短的DNA片段)，分割方法如图1-1所示。假定集合R中任意一条read的长度均为l，k-mer长度值设为k，那么集合R中的任意一条read均可被分为l–k+1条k-mer，并且这些k-mer作为de Bruijn图的顶点。
（2）对于给定的两条k-mer x和y，如果在某read ri中存在一条长度为k+1的子串，且该子串的前k个碱基与k-mer x(或y)精确匹配，同时该子串的后k个碱基与k-mer y(或x)精确匹配，那么该算法认为两条k-mer x和y之间存在一条公共边。
将采用上述方法构造的de Bruijn图记作G。对于read集合R={r1，r2，…，rn}中的任意一条read ri，若在de Bruijn图G中存在一条路径P，且该路径P 访问ri中的每一条k-mer仅一次，则欧拉路径问题便可理解为：给定某一de Bruijn图G以及G中的路径集合P，在de Bruijn图G中确定某一条欧拉路径Q，使得路径集合P中的每一个元素都是欧拉路径Q的子路径。利用欧拉路径算法进行DNA序列拼接的主要步骤如下所述：首先利用纠错软件修正read中测序错误的碱基；然后按照上述方法构建de Bruijn图；构建deBruijn图之后，应将read 集合中的所有read排列在de Bruijn图中，在deBruijn图中，每一条read均被视作一条路径；最后在de Bruijn图中寻找一条欧拉路径，使得该路径包含de Bruijn图中所有read所对应的路径。
在OLC中,在Overlap步骤中,采用了序列比对算法来寻找read之间的重叠信息,该算法的时间复杂度为0(?2),其中,《SDNA序列中read的数量。当前DNA 测序数据序列越来越短,对同一个物种进行测序,其产生的read数量大大增加,这使得OLC的计算量增加;而基于deBruijn图原理的序列拼接中,抛弃了 OLC中序列比对算法,而是采用以k-mer为图中顶点构建图,从而减少了序列比对算法所消耗的时间,提高了算法的效率与overlap-layout-consensus算法相比，基于de bruijn图的算法有更低的时间复杂度，这是因为欧拉路径问题实际上是一个线性时间的问题。利用欧拉路径思想的拼接算法有EULER-SR、ALLPATHS、Velvet
和EULER等。
四．模型建立
4.1.1模型的假设
1.假设模型中的read片段都是由一条完整的DNA经过测序而来，它们进过拼接后可以形成一个大片段。
2.模型中出现的各个序列中DNA的双链都准确
3.模型中read在拼接时合理地去掉的公共部分在误差允许的范围内。
4.由于总会在测序中出现read的碱基错误，因此，假设这少量的错误在模型求解时时在误差允许的范围内的。
5.在基因组的剪切过程中未发生基因的丢失，DNA改变，基因的重叠等4.1.2数据在拼接的预处理
Reads在拼接时，由于新一代序列数据很多，准确度较低，导致reads中含有大量错误碱基。在这种错误下，de Bruijn图的实际大小会随着reads数据量的增加呈现指数型增长，并且容易造成错误拼接。因此，在此之前需要对reads 进行预处理，修正或消除初始reads中的碱基错误。
（1）新一代测序数据错误率高，且主要分布在靠近reads3’端部分，并且越靠近3’出错率越高，而5’端比较正确，如图3-1[2]所示。为减少错误，我们的方法是：计算3’端reads长度一般的碱基的平均质量，过滤掉该区域平均质量小于15的reads。该质量值对应碱基的出错率，计算公式如下：Q=－10×lg(ε)。其中，Q为碱基的质量只值，ε为碱基的出错率。
图4-1solexa测序数据的错误特征
(2)测序数据中会有许多全A或者基本上全A的reads，这些数据很可能是solexa 测序过程中的人工数据，需要去除。方法为：该定A的含量阀值为0.9，过滤掉A含量大于0.9的reads。
（3）测序数据中含有一些未知的碱基，通常同“N”或“.”表示，代表没有被测出来的碱基，对拼接有不利的影响，因此含有未知碱基的reads需要过滤掉。
4.2算法的步骤
4.2.1de Bruijn图策略
基于de Bruijn图策略的拼接算法的做法是：
（1）构建de Bruijn图：将reads分割成一系列连续的子串k-mers(一般用k只表示kmer碱基数目的大小)，作为图中的边，相邻的两个k-mers交叠（k-1）个碱基；
（2）化简de Bruijn图：方法是合并路径出度入度唯一的节点，按照一定的规则去除图中的尖端（tips）和泡状结构（bubbles），如图3-2所示[3]；
（3）构建contigs：在de Bruijn图或其子图中寻找一条最优的欧拉路径（经过每条边一次且仅一次的路径），该路径对应的碱基序列即为contigs；
（4）生成scaffolding：利用配对数据，确定contigs之间的相对方向与位置关系，对contigs进行组装，并填充contigs之间的gaps，最终得到scaffolds序列。
图4-2 de Bruijn图的简化
4.2.2 de Bruijn图的构建
与其他基于de Bruijn图拼接算法不同的是，基于reads引导的基因序列拼接以de Bruijn图架构为基础，这里所建立的de Bruijn图只有孤立的顶点，而没有边的存在。基于de Bruijn图的数据存储方法具有快速高效的优点，然而内存开销比较大，内存是计算机的宝贵资源，直接制约着拼接的正常进行，因此，设计新的内存高效的存储方法，对拼接具有重要的影响。
（1）能够灵活操作内存，我们采用C语言作为拼接算法的实现工具。
（2）一条reads，不额外生成其反向互补，而是在拼接过程中，自动检测该reads是正向还是反向。这样就避免了数据的二倍增长，以及不必要的数据冗余和计算资源的浪费。
（3）设计高效的数据结构，de Bruijn图的主要数据结果如图4-3[4]所示。
图（a）表示的是如何在一个read上获取kmers，kmer是reads上的连续子串，其中一般用k值表示kmer的大小。从reads上的首个碱基开始，依次截取k个碱基。由于reads数据的高覆盖度及k值一般比较小，这就造成在获取kmer 的过程中会出现大量一样的kmer。因此，就没个kmer而言，需要记录出现的reads 编号和reads上几号位置，还需要对其出现的次数累加。
图（b）展现了reads1和reads2在de Bruijn图中的存在形式，每个read 在de Bruijn图中都对应着一个kmer通路。
图（c）表示的是de Bruijn图的数据结构。
图4-3 de Bruijn图
4.2.3contigs构建
综合考虑参与拼接的reads的各项信息，从reads整体角度研究拼接问题，是基于reads引导的全基因组段序列拼接算法的基本思想。其主要做法是：将正在参与拼接的reads的拼接信息记录到决策表中，综合考虑决策表中reads累计拼接信息和数据的出错特征，设计合理的评分方法，选择得分最大的kmer进行contig的扩展。根据选择的kmer更新决策表中reads的信息，contigs每次扩展一个碱基，第一轮扩展结束后，取其反向互补，进行反向扩展。
基于reads引导的contigs构建过程如下：
1.轮拼接，初始化决策表为空，在de Bruijn图中选择拼接的初始kmer作为初始拼接的contigs，将参与拼接的reads信息添加进决策表；
2.扩展候选的k-mers，根据决策表中的reads信息，对这些kmer评分，选择得分最高的k-mer；
3.候选k-mers为空，并且已经处于第二轮拼接，则拼接停止，转到第6步；如果仍处于第一轮拼接，则标记为第二轮拼接，初始化决策表为空，将contigs 反向互补，在新的3’端重新选择k-mers进行拼接；
4.选择的k-mer更新决策表中的reads的拼接信息，删除拼接成功或者失败的reads；
5.在de Bruijn图中标记拼接成功的reads的所有k-mers为删除状态，将成功的reads添加进contigs；
6.拼接结束，则保存contig及成功的reads。
该算法的拼接示意图如图4-4[5]所示。经过打分导航机制确定了一个带扩展的新kmer。该kmer在rid为3的当前read中没有出现，因此它的mismatches 增1，last pos 变为0.而该kmer在rid为6的kmer中出现了，则它的matches 增1，last pos 增1 ，变为14.如果此时rid为6的last pos不是13而是0，则last由0变为14.rid为2和8的read是新引进决策表的read。
图4-4contigs拼接示意（read长为36bp，k值取12）
4.2.4 contigs组装
本阶段以reads拼接阶段生成的contigs和配对数据文件作为输入，通过配对数据确定contigs之间的相对位置，连接contigs，填充相邻contigs间的gaps 区域等过程，最终生成长度更长的scaffolds ，如图4-5[6]所示。
图4-5 contigs组装示意图
4.3.1模型配对文库参数的校正：
在contigs组装阶段，配对数据常用语确定相邻两个contigs的相对方向和位置关系，校正contigs中的序列重排，以及填充相邻的两个contigs间gaps 区域，最终生成质量更高的超长序列scaffolds。在该阶段，reads被映射到contigs上，以获得配对reads在contigs中的位置。如果配对的两条reads分别出现在不同的contigs中，则说明这两个contigs在目标基因组中是相邻的。Contigs组装的关键是利用配对数据确定contigs的相对方向和位置关系，通常采用基于图的方法解决该问题。在图中，contigs作为顶点，contigs之间的配对连接作为边，配对俩接的数量作为边上的权重，这样就建立了任意两条contigs之间的配对关系，该问题就转化为在图中或子图中寻找遍历每个节点的最优路径。通常与给定的contigs具有配对关系的contigs数量是很少的，因此，该图并不复杂，处理起来相对简单。然而，在实际中，由于目标基因组中的重复片段以及可能测序错误的影响，contigs之间的连接有一些是不真实的（或错误的），为了保证结果的准确性，拼接算法通常忽略配对数量低于一定阀值的contigs之间的连接。
通常，配对文库的片段大小服从正态分布，并且标准差大约在均值的10%左右。在实际应用中，往往只给出配对文库的均值，而标准差是未知的。然而，配对文库的片段尺寸的分布通常与事先给定的分布并不完全相同，总会有一些差距。因此，我们对配对数据的分布参数进行了校正。
首先，随机选取若干条contigs作为参数校正的contigs样本并对数据文件的reads进行编号。一般而言，配对数据会存储在两个文件里，对数据文件1
的reads从1开始用奇数编号，对数据文件2的reads从2开始用偶数编号，这样，1、2号是配对的，3、4号是配对的，以此类推。
然后，对选取的contigs按read的长度建立索引。把两个文件里的reads
往contigs上映射。通常，如果reads数据文件太大的话，可以截取文件1和文件2相对应的各一部分，再往选取的contigs上映射。
最后，分析映射结果，计算配对数据的分布参数，提取reads数据的映射信息，并按reads编号由小到大排列，这样便于分析映射到每一条contigs上的配对数据之间的相关信息。图3-5-1[7]展现了对两对配对数据映射到contigs上
的情形。
图4-5-1配对数据映射到contigs
4.3.2 contigs相对位置的确定
在contigs组装之前，需要先无额定contigs之间的相对位置和方向。将配对数据往reads拼接阶段生成的contigs上映射，通常会有若干配对的reads
映射到不同的contigs上，如果被配对的reads映射上的两条contigs之间有交叠，则看配对的两条reads在它们上的位置之间的距离是否在合理范围内。如果距离在合理范围内且配对reads数目达到设定的阀值，则认为两条contigs是相邻的，应该连接到一起。如果这两个contigs没有交叠，则认为它们之间存在间隙gaps，需要后续的gaps填充操作。
在contigs的构建阶段，拼接错误通常发生在contigs的末端，这是因为contigs是因找不到符合田间的下一个kmer才结束拼接的，contigs靠近末端的若干碱基正确率较低。因此，我们摒弃了contigs组装通常采用的把配对reads 数据往整条contigs上映射的做法，而只取contigs的两端各Lbp长的序列片段，这样既减少了内存消耗，也有利于提高计算速度。如果contigs的长度小于L，则取其整个序列。L值的确定于配对文库的平均插入距离有关，一般取值要大于配对文库的平均插入距离，比如如果所用配对文库的平均插入距离为200bp的话，那么L就取值300。
（1）建立contigs索引；
（2）映射配对reads和contigs之间的匹配信息；
（3）确定contigs的相对位置。
4.3.3 contigs连接
contigs相对位置的确定指的是contigs排放位置的确定，并没真正的连接在一起。contigs连接时需要先计算contigs之间的交叠overlap，由于contigs
之间的相对位置已经确定，我们可以根据相邻两条contigs上配对的reads来计算contigs之间的距离distance，从而确定overlap的大小。此时的distance 大小可以大于0，也可以小于0。若计算出的distance值小于0，则表明这两条contigs之间可能有overlap存在；若计算出的distance值大于0，则表明这两条contigs之间可能有真正的gaps存在。
4.3.4 gaps的填充
contigs之间的相对位置确定后，在contigs连接阶段，有交叠的contigs 被连接到一起，合并成长度更长的contigs，没有被连接的相邻contigs之间有空隙gaps的存在，需要进行gaps填充。gaps填充就是通过局部序列拼接的方式将contigs之间缺少的那部分碱基序列构造出来。配对数据往拼接生成的contigs上映射时，有这样的一部分配对的reads：它们只有其中的一条找到了匹配位置，而与之配对的另一条没有找到，我们称配对reads中在contigs上找到匹配位置的那些reads为悬浮hanging reads。针对每两条相邻的contigs，我们将悬浮的那些reads保存在一个表中，该表记为悬浮reads表HRL。与此同时，收集那些与hang reads对应的在contigs上未找到匹配位置的reads，并用它们构建出一个De Bruijn图。有了悬浮reads表HRL和所构建出来的De Bruijn图，我们就可以很快速的对contigs之间的gaps进行填充。如图
4-5-4.1[8]所示。
图4-5-4.1gaps填充示意图
如图4-5-4.2[9]所示，选取contigs A末端READ_LEN长的序列处开始拼接，每次扩展一个碱基。待拼接kmer有效的条件：悬浮reads表中存在与kmer所在read的配对read，并且该kmer所在的read的方向与contigs的扩展方向相反以及该配对reads的距离在配对文库的误差范围之内。最优kmer的确定方法与reads拼接阶段的kmer选取方法一样，就是通过打分机制，对待选的kmer进行打分，取得分最大者。在reads拼接阶段，如果下一个kmer为空时，则停止contigs的首轮扩展，并进行contigs的反向扩展，而在填充gaps时，当下一个kmer为空时，我们把N作为扩展字符，只进行单向扩展。
图4-5-4.2局部拼接填充gaps
为了加快拼接速度，我们在实际操作时选择从gaps两端同时相向拼接，即在相邻contigs对应的两末端同时开始在悬浮reads表和De Bruijn图的引导下对gaps进行填充。在填充gaps的同时，需要时刻检测这两条contigs之间的gaps信息，一旦gaps的大小接近0时就停止填充，以防止过填充情形的出现。
4.5 算法的小结
本节详细介绍了基于reads引导的基因组序列拼接算法，该算法的实施具体包括reads拼接和contigs组装两个阶段。reads拼接阶段主要有数据的预处理、De Bruijn图的建立和contigs构建等几部分，其中重点介绍了拼接时下一个kmer的选择策略，即以配对reads数据为拼接导航、综合考虑整条reads参与拼接的累计信息，并结合reads数据的区域特征提出打分机制，对待选的所有kmer进行打分，取得分最大者。contigs组装阶段主要介绍了contigs相对位置的确定、contigs连接和gaps填充等几方面的内容，就contigs之间序列交叠overlap和间隙gaps的界定，以及contigs末端出错碱基的处理及纠正等做了深入的探究。高效的数据结构，巧妙的拼接组装策略等，为该算法的成功实施提供了重要保障。
五.模型中算法的验证
在基于reads引导的基因组序列拼接算法理论的指导下，我们结合实际应用的需要，开发出一套基因组序列拼接系统，该序列拼接系统被命名为SRGA，经数据测试取得了不错的效果，从而验证了算法的正确性。
图5-1是用相关软件做出的基因组比对相似度的一部分。
图5-1.基因组的相似度比对
5.1算法的输入与输出
SRGA序列拼接系统主要包括两个大的功能模块，reads拼接模块和contigs 组装模块。这两个模块相对独立，作用于输入数据的不同处理阶段，reads拼接模块生成的contigs是contigs组装模块的输入。在算法实现时，我们选择了C 语言作为系统开发语言，并将软件在linux环境下发布。
新一代测序技术的快速发展，产生了大量的配对reads数据，而这些配对reads的数据文件就是系统的初始输入。如图4-1所示，reads配对数据通常分布在两个数据文件，在每个单独的文件里，每4行作为一个基本单位，记录着一条reads的序列信息和碱基质量信息。其中，第1行和第3行是测序信息，第2行是reads的碱基序列，第4行是reads每个碱基的质量值，reads的碱基序列和碱基质量值是基因组序列拼接所需要的数据信息。配对reads所分布的两个文件里对应的reads相互配对，即第一个文件里的第1条reads与第二个文件的第1条reads配对，第一个文件里的第2条reads与第二个文件的第2条reads配对，依次类推。
就系统的初始输入而言，除了配对reads数据文件，还有一个重要的参数需要人为设定，即K值。K值表征kmer的大小，在不同K值条件下，基因组序列的拼接效果往往相差很大。K值的选择还没有通用的方法，一般是通过综合考虑reads长度、数据覆盖度、目标基因组长度等因素而设定。
图5-2reads拼接阶段生成的contigs示意
初始输入的配对reads经过SRGA的reads拼接功能模块处理后，会生成大量的长度更长的序列片段contigs。如图4-2所示，生产的contigs包括contigs 的编号，contigs的长度及contigs的碱基序列等。contigs只是序列拼接的中间结果，将与初始输入的配对reads一起作为系统contigs组装模块的输入。
图5-3 contigs组装阶段生成scaffolds示意
reads拼接模块生成的contigs和初始输入的配对reads数据经过contigs 组装功能模块的作用后，会生成若干条序列长达几十kbp甚至几百kbp的基因组序列片段scaffolds。如图4-3所示，生成的scaffolds包括scaffolds的编号，拼接成scaffolds的contigs数目及scaffolds碱基序列等。
scaffolds是SRGA序列拼接系统生成的最终输出结果，scaffolds数据生产后，科研人员就可以在对其进行一系列数据分析操作，比如查找功能相似的基因序列、分析单碱基变异等，为生命科学的研究提供帮助。
六评价结果及分析：
模型在建立过程中，经过严密的运算和验证后，总结出其优点如下：
⑴在reads拼接方面，从整体reads的角度入手，综合考虑拼接过程中的累计信息，同时又参考了reads的数据特征，为contigs的拼接提供了全面而客观的评分标准。并优先考虑有配对的reads数据，能够高精确性地为contigs的拼接导航，避免了传统的基于k-mers的拼接计算
⑵在contigs组装方面，提出了PEN数组的概念，用于contigs之间相对位置的确定。相邻contigs连接时，采用启发式搜索策略，进行线性扩展，降低了计算复杂度。并用序列对比的方法来检测contigs之间的overlap，对contigs 末端的碱基错误进行纠正，较好地消除了contigs末端拼接错误的影响，大大提高了contigs组装的准确性。
⑶在空间和时间复杂度上，本算法都优于其它的算法，较低些。
本文提出并实现了基于reads引导的基因组序列拼接算法，经数据验证，该算法表现出色，但仍然存在以下不足之处，有待在后续研究工作中加以改进。其不足之处如下：
⑴该算法不支持多配对文库。配对数据在序列拼接时有着重要应用，多组配对数据的使用能大大提高序列的拼接质量。但目前系统的输入只能是一个配对文库，需要加以扩展。
⑵选择下一个kmer的打分函数不够精确。当数据覆盖度较高时，由于决策表被锁定的reads会很多，造成同时多个待选kmer的得分值都很大，难以取舍，需要在区域得分设置及总分量化等方面进一步优化。
七.参考文献：
参考学位论文：
[1]韩东涛.基于概率模型的基因组从头测序算法研究. 哈尔滨工业大学，2012
[2]曾培龙.基于reads引导的基因组序列拼接. 哈尔滨工业大学,2012。
[3]杨帆 .基于 BWT 的 DNA 重叠群序列合并算法研究,哈尔滨工业大学，2012。
[4]范建蕊.基于双向de Bruiin图的序列拼接并行化研究与实现，中南大学，2012。
[5]林勇.面向下一代测序技术的 de novo 序列拼接工具综述,上海理工大学，2012.
参考书籍：
《数学模型》（第三版）姜启源谢金星叶俊编，高等教育出版社。
《全国大学生数学建模竞赛优秀论文汇编》全国大学生数学建模竞赛组委员会编，中国物价出版社。
《基因工程》梁文蔚等编，高教出版社，2011.
八.附录
另一文件夹有相关的程序和基因组装软件以及数据。
基因组学的研究内容 结构基因组学： 基因定位；基因组作图；测定核苷酸序列 功能基因组学：又称后基因组学（postgenomics基因的识别、鉴定、克隆；基因结构、功能及其相互关系；基因表达调控的研究 蛋白质组学： 鉴定蛋白质的产生过程、结构、功能和相互作用方式 遗传图谱 （genetic map）采用遗传分析的方法将基因或其它dNA序列标定在染色体上构建连锁图。 遗传标记： 有可以识别的标记，才能确定目标的方位及彼此之间的相对位置。 构建遗传图谱 就是寻找基因组不同位置上的特征标记。包括： 形态标记； 细胞学标记； 生化标记；DNA 分子标记 所有的标记都必须具有多态性！所有多态性都是基因突变的结果！ 形态标记： 形态性状：株高、颜色、白化症等，又称表型标记。 数量少，很多突变是致死的，受环境、生育期等因素的影响 控制性状的其实是基因，所以形态标记实质上就是基因标记。
细胞学标记 明确显示遗传多态性的染色体结构特征和数量特征 ：染色体的核型、染色体的带型、染色 体的结构变异、染色体的数目变异。优点：不受环境影响。缺点：数量少、费力、费时、对生物体的生长发育不利 生化标记 又称蛋白质标记 就是利用蛋白质的多态性作为遗传标记。 如：同工酶、贮藏蛋白 优点： 数量较多，受环境影响小 ?
缺点： 受发育时间的影响、有组织特异性、只反映基因编码区的信息 DNA 分子标记： 简称分子标记以 DNA 序列的多态性作为遗传标记 优点： ? 不受时间和环境的限制 ? 遍布整个基因组，数量无限 ?
不影响性状表达 ? 自然存在的变异丰富，多态性好 ? 共显性，能鉴别纯合体和杂合体 限制性片段长度多态性（restriction fragment length polymorphism ， RFLP ） DNA 序列能或不能被某一酶酶切，
DNA测序结果分析比对（实例） 关键词：dna测序结果2013-08-22 11:59来源：互联网点击次数：14423 从测序公司得到的一份DNA测序结果通常包含.seq格式的测序结果序列文本和.ab1格式的测序图两个文件，下面是一份测序结果的实例： CYP3A4-E1-1-1(E1B).ab1 CYP3A4-E1-1-1(E1B).seq .seq文件可以用系统自带的记事本程序打开，.ab1文件需要用专门的软件打开。软件名称：Chromas 软件Chromas下载 .seq文件打开后如下图： .ab1文件打开后如下图： 通常一份测序结果图由红、黑、绿和蓝色测序峰组成，代表不同的碱基序列。测序图的两端（下图原图的后半段被剪切掉了）大约50个碱
基的测序图部分通常杂质的干扰较大，无法判读，这是正常现象。这也提醒我们在做引物设计时，要避免将所研究的位点离PCR序列的两端太近（通常要大于50个碱基距离），以免测序后难以分析比对。 我的课题是研究基因多态性的，因此下面要介绍的内容也主要以判读测序图中的等位基因突变位点为主。 实际上，要在一份测序图中找到真正确实的等位基因多态位点并不是一件容易的事情。一般认为等位基因位点假如在测序图上出现像套叠的两个峰，就是杂合子位点。实际比对后才知道，情况并非那么简单，下面测序图中标出的两个套峰均不是杂合子位点，如图并说明如下：
说明： 第一组套峰，两峰的轴线并不在同一位置，左侧的T峰是干扰峰；第二组套峰，虽两峰轴线位置相同，但两峰的位置太靠近了，不是杂合子峰，蓝色的C峰是干扰峰通常的杂合子峰由一高一略低的两个轴线相同的峰组成，此处的序列被机器误判为“C”，实际的序列应为“A”，通常一个高大碱基峰的前面 1~2个位点很容易产生一个相同碱基的干扰峰，峰的高度大约是高大碱基峰的1/2，离得越近受干扰越大。 一个摸索出来的规律是：主峰通常在干扰峰的右侧，干扰峰并不一定比主峰低。最关键的一点是一定要拿疑似为杂合子峰的测序图位点与测序结果的文本序列和基因库中的比对结果相比较；一个位点的多个样本相比较；你得出的该位点的突变率与权威文献或数据库中的突变率相比较。 通常，对于一个疑似突变位点来说，即使是国际上权威组织大样本的测序结果中都没有报道的话，那么单纯通过测序结果就判定它是突变点，是并不严谨的，因一份 PCR产物中各个碱基的实际含量并不相同，很难避免不产生误差的。对于一个未知突变位点的发现，通常还需要用到更精确的酶切技术。 (责任编辑：大汉昆仑王)
2014年成都理工大学校内数学建模竞赛论文 二0一四年五月二十五日
摘要：本文所要研究的就是全基因组的从头测序的组装问题。 首先，本文简要介绍了测序技术及测序策略，认真分析了基因系列拼装所面临的主要挑战，比如reads数据海量、可能出现的个别碱基对识别错误、基因组中存在重复片段等复杂情况，探讨了当前基因组序列拼接所采用的主要策略，即OLC（Overlap/Layout/Consensus）方法、de Bruijn图方法，且深入探讨了de Bruijn图方法。 其次，针对题中问题，以一条reads为基本单位，分为reads拼接和contig组装两个阶段，其中contig是由reads拼接生成的长序列片段。Reads的拼接阶段主要包括数据预处理、de-Bruijn 图、contig构建等，而contig的组装阶段主要包括序列的相对位置的确定以及重叠部分overlap的检测，用序列比对的方法来提高拼接的精度。 最后，进行了算法的验证与性能的评价，并且针对问题2，进行了组装分析与验证，结果表明，得到的拼接基因组序列在小范围内与原基因组序列大致吻合。 关键词：基因组系列拼接； reads；de Bruijn图；contig组装；k-mer片段；
一.问题重述 基因组组装 快速和准确地获取生物体的遗传信息对于生命科学研究具有重要的意义。对每个生物体来说，基因组包含了整个生物体的遗传信息，这些信息通常由组成基因组的DNA或RNA分子中碱基对的排列顺序所决定。获得目标生物基因组的序列信息，进而比较全面地揭示基因组的复杂性和多样性，成为生命科学领域的重要研究内容。 确定基因组碱基对序列的过程称为测序（sequencing）。测序技术始于20世纪70年代，伴随着人类基因组计划的实施而突飞猛进。从第一代到现在普遍应用的第二代，以及近年来正在兴起的第三代，测序技术正向着高通量、低成本的方向发展。尽管如此，目前能直接读取的碱基对序列长度远小于基因组序列长度，因此需要利用一定的方法将测序得到的短片段序列组装成更长的序列。通常的做法是，将基因组复制若干份，无规律地分断成短片段后进行测序，然后寻找测得的不同短片段序列之间的重合部分，并利用这些信息进行组装。例如，若有两个短片段序列分别为 ATACCTT GCTAGCGT GCTAGCGT AGGTCTGA 则有可能基因组序列中包含有ATACCTT GCTAGCGT AGGTCTGA这一段。当然，由于技术的限制和实际情况的复杂性，最终组装得到的序列与真实基因组序列之间仍可能存在差异，甚至只能得到若干条无法进一步连接起来的序列。对组装效果的评价主要依据组装序列的连续性、完整性和准确性。连续性要求组装得到的（多条）序列长度尽可能长；完整性要求组装序列的总长度占基因组序列长度的比例尽可能大；准确性要求组装序列与真实序列尽可能符合。 利用现有的测序技术，可按一定的测序策略获得长度约为50–100个碱基对的序列，称为读长（reads）。基因组复制份数约为50–100。基因组组装软件可根据得到的所有读长组装成基因组，这些软件的核心是某个组装算法。常用的组装算法主要基于OLC（Overlap/Layout/Consensus）方法、贪婪图方法、de Bruijn 图方法等。一个好的算法应具备组装效果好、时间短、内存小等特点。新一代测序技术在高通量、低成本的同时也带来了错误率略有增加、读长较短等缺点，现有算法的性能还有较大的改善空间。 问题一：试建立数学模型，设计算法并编制程序，将读长序列组装成基因组。你的算法和程序应能较好地解决测序中可能出现的个别碱基对识别错误、基因组中存在重复片段等复杂情况。 问题二：现有一个全长约为120,000个碱基对的细菌人工染色体（BAC），采用Hiseq2000测序仪进行测序，测序策略以及数据格式的简要说明见附录一和附录二，测得的读长数据见附录三，测序深度（sequencing depth）约为70×，即基因组每个位置平均被测到约70次。试利用你的算法和程序进行组装，并使之具有良好的组装效果。 附录一：测序策略 测序策略如下图所示。DNA分子由两条单链组成，在图中表现为两条平行直
基因组序列的差异分析 ----mVISTA的在线使用说明 当然，除了在线版的，我们还可以在网站上填写信息申请离线的软件。但我试用了一下，需要先自己比对，然后要按照一定的格式来制作文件，当然你还必须得安装java才能运行软件；总之，我感觉没有在线版的方便。 1 将数据放入服务器中 在首页，你将被要求确定你想要分析的基因组序列的数量。输入这个数字之后，点击“提交”，将带你到主提交页面。 mVISTA服务器最多可以同时处理100条序列。 1.1主提交页面必填的内容 E-mail 地址 通过E-mail，我们可以提示你的在线处理已经得到结果。
序列 你可以用2种方式来上传你的序列： 1.使用“Browse”按钮从你的电脑上，上传纯文本的Fasta格式文件。如果是一个作为参 考的生物体的DNA序列必须作为一个contig提交(可以进行一定的定向排列将多个片段合并为一个contig)，而其他非参考序列可以在一个或多个contig中提交(draft)。 Fasta格式的示例序列(您可以在NCBI站点上找到关于该格式的更多细节)： >mouse ATCACGCTCTTTGTACACTCCGCCATCTCTCTCT … ！！！注意:序列里面我们只接受字母CAGTN和X。请确保提交序列是作为一种纯文本格式，而不是Word或HTML文件格式。 如果您以FASTA格式提交序列，我们建议您为它取一个有意义的名称（比如直接是你的物种名之类的），因为这些名称将出现在我们生成的图形中。如果您使用的是一个draft草图序列，那么结果中每个contigs的命名都将按照您在“>”符号后指示的命名进行。 2.您可以给出它的GenBank登录号，系统将自动从GenBank数据库里进行检索序列。 在这两种情况下，序列的总大小都不应超过10M，而且任何一条序列都不应超过2M。 1.2主提交页面选填的内容 这些选项允许您自定义您的VISTA分析。您可以使用独立获得的基因注释，选择合适的Repeat Masker选项，给分析的序列指定名称，并改变序列保存分析的参数。如果您没有填写这些选填选项，我们将使用它们的默认值。 比对程序 根据您分析的具体内容(参见“about”-链接中的详细信息)，您可以选择以下比对程序之一：1、AVID----全局两两比对。如果您选择使用这个程序，其中一个序列应该被完成比对，其他 所有序列可以完成或以草图draft格式完成。对于集合中所有已完成的序列，AVID生成所有相对所有成对的比对结果，可以使用任何序列作为基础(参考)来显示。如果某些序列是草图格式，AVID将生成它们与最终序列的比对，这将被用作基础(参考)。这是该服务器上唯一可以处理草图序列的比对程序。 （小知识：草图序列与完整序列DNA sequence, draft: Sequence of a DNA with less accuracy than a finished sequence. In a draft sequence, some segments are missing or are in the wrong order or are oriented incorrectly. A draft sequence is as opposed to a finished DNA sequence.）2、LAGAN----完成完整序列的全局两两比对和多重比对。如果某些序列是草图格式，您的查 询将被重定向到AVID以获得两两比对。多重比对将由VISTA可视化，它将计算并显示序列的保守区，以您指示的任何序列作为参考。这是该服务器上唯一能够产生真正的多重
学院：______ 班级:_______ 学号:_________ 姓名:__________ 成绩：______ 实验五基因结构预测分析 目的： 1、熟悉并掌握从基因组核酸序列中发现基因的方法。 内容： 1、用NCBI的ORF Finder分析原核生物核酸序列或真核生物的cDNA序列中的开放阅读框； 2、使用GENSCAN在线软件预测真核生物基因； 3、使用POL YAH在线预测转录终止信号； 4、使用PromoterScan在线预测启动子区域。 操作及问题： 随着测序技术的不断发展，越来越多的模式生物启动了全基因组测序计划，完成全基因组测序的物种也越来越多，使得基因结构和功能的预测成为可能。同时，通过基因组文库筛选也可得到目的基因所在克隆。获得克隆序列后，同样也需要对目的基因做结构预测以便指导后续功能研究。本实验介绍几种常用的基因预测分析工具，预测核酸序列的开放阅读框、转录终止信号、启动子、CpG岛等信息。 一、开放阅读框（open reading frame，ORF）的识别 ORF是指从核酸序列上5’端翻译起始密码子到终止密码子的蛋白质编码序列。原核生物与真核生物的基因结构存在很大不同，真核生物的ORF除外显子（平均150bp）外，还含有内含子，因此真核生物基因的预测远比原核生物复杂。 （一）利用NCBI ORF Finder预测原核生物核酸序列或真核生物的cDNA序列中的开放阅读框。https://www.doczj.com/doc/6213834142.html,/gorf/gorf.html 1、在NCBI上查找AC 号为AE008569 的核酸记录。（见实验五中的AE008569.mht） 问题1：这个序列的名称？ 问题2：这个序列来源物种所属的生物学大分类？
人类基因组图谱定义 1543年，比利时解剖学家A·维萨里（1514－1564）发表了划时代的著作《人体的构造》，开创了人体解剖学，使人们从宏观上了解了自己。“人类基因组计划”建立的人类基因组图，被誉为“人体的第二张解剖图”，它将从微观上或者说从根本上使人类了解自己。 人类第一个基因组草图 2000年6月26日，美国总统克林顿和英国首相布莱尔联合宣布：人类有史以来的第一个基因组草图已经完成。2001年2月12日中、美、日、德、法、英等6国科学家和美国塞莱拉公司联合公布人类基因组图谱及初步分析结果。 人类基因组计划中最实质的内容，就是人类基因组的DNA序列图，人类基因组计划起始、争论焦点、主要分歧、竞争主战场等都是围绕序列图展开的。在序列图完成之前，其他各图都是序列图的铺垫。也就是说，只有序列图的诞生才标志着整个人类基因组计划工作的完成。 2003年4月15日，在DNA双螺旋结构模型发表50周年前夕，中、美、日、英、法、德六国元首或政府首脑签署文件，六国科学家联合宣布：人类基因组序列图完成。 人类基因组图谱的绘就，是人类探索自身奥秘史上的一个重要里程碑。它被很多分析家认为是生物技术世纪诞生的标志，也就是说，21世纪是生物技术主宰世界的世纪。正如一个世纪前量子论的诞生被认为揭开了物理学主宰的20世纪一样。 全球专家拟绘癌症基因图谱 国际癌症基因组协会4月29日在英国伦敦成立。这一组织计划通过统筹各国和地区专家的合作，耗资10亿美元，历时10年，绘制较为完整的致癌基因突变图谱。目前已有英国、中国和美国等9国加入这一计划。 专家认为，图谱将为癌症预防、诊断和治疗带来一场革命，开辟癌症个案化治疗的新时代。 方法 国际癌症基因组协会计划利用更加先进、快速的基因组测序新技术，详细研究50种不同类别的癌症，希望找到所有与癌症相关的基因突变现象。 英国剑桥韦尔科姆基金会桑格研究所是这一计划的主要参与者之一。 英国《泰晤士报》援引桑格研究所专家迈克·斯特拉顿的话说：“借助更快速的脱氧核糖核酸（ＤＮＡ）测序新技术，国际癌症基因组协会雄心勃勃，志在为数以千计癌症基因组测序，制作一个涵盖所有ＤＮＡ变异的目录，绘制完整的癌变图谱。” 越来越多科研成果表明，癌症并非单一类型。而每种癌症都包括大量由不同类型基因突变导致的不同亚种类别，需要分别对症治疗。 国际癌症基因组协会计划在每种癌症的研究中，提取500名病人的细胞并测序基因组，与健康细胞作比对，以期找出导致癌细胞形成和扩散的基因突变。 意义 专家认为，图谱绘制对于癌症治疗具有革命性意义。 如果图谱绘制成功，医生就可准确掌握单个患癌病人的致癌基因突变因素，进而更为轻松地对症下药，寻找对特定病人或特定癌症种类具有针对性的治疗药物。 比如，医学界目前普遍认为，赫赛汀（Herceptin）是特定种类乳癌的治疗药物。 图谱还有助于发明新的癌症治疗药物。 韦尔科姆基金会负责人马克·沃尔波特说，识别致癌基因突变是癌症治疗领域内的一次“大跨步”进展，目的在于实现“对症治疗”。 沃尔波特披露，桑格研究所的斯特拉顿等专家已启动一项致力于研发癌症治疗新药的计划。 合作 国际癌症基因组协会成立的作用还在于促进各国专家合作。 目前已确定加入这一计划的国家为英国、中国、美国、澳大利亚、加拿大、法国、印度、日本和新加坡。 协会计划对每一种类癌症的研究资助2000万美元，因而整个研究项目将耗资约10亿美元。
实验三蛋白序列比对到基因组（GeneWise and exonerate）实验目的 1）了解基因结构，acceptor, sponsor 等概念 2）理解将蛋白序列比对到基因组的应用 3）掌握利用GeneWise 将蛋白序列定位到基因组上并得到基因结构 实验数据及软件 ftp://172.28.137.55/pub/lab_materia/biosoft/lab03/ 1、Genewise 简介 Genewise 是EBI 的Ewan Birney
实验三核酸序列分析 【实验目的】 1、掌握已知或未知序列接受号的核酸序列检索的基本步骤； 2、掌握使用BioEdit软件进行核酸序列的基本分析； 3、熟悉基于核酸序列比对分析的真核基因结构分析（内含子/外显子分析）； 4、了解基因的电子表达谱分析； 5、熟悉密码子偏好性分析。 【实验原理】 针对核酸序列的分析就是在核酸序列中寻找基因，找出基因的位置和功能位点的位置，以及标记已知的序列模式等过程。在此过程中，确认一段DNA序列是一个基因需要有多个证据的支持。一般而言，在重复片段频繁出现的区域里，基因编码区和调控区不太可能出现；如果某段DNA片段的假想产物与某个已知的蛋白质或其它基因的产物具有较高序列相似性的话，那么这个DNA片段就非常可能属于外显子片段；在一段DNA序列上出现统计上的规律性，即所谓的“密码子偏好性”，也是说明这段DNA是蛋白质编码区的有力证据；其它的证据包括与“模板”序列的模式相匹配、简单序列模式如TATA Box等相匹配等。一般而言，确定基因的位置和结构需要多个方法综合运用，而且需要遵循一定的规则：对于真核生物序列，在进行预测之前先要进行重复序列分析，把重复序列标记出来并除去；选用预测程序时要注意程序的物种特异性；要弄清程序适用的是基因组序列还是cDNA序列；很多程序对序列长度也有要求，有的程序只适用于长序列，而对EST这类残缺的序列则不适用。 1. 重复序列分析 对于真核生物的核酸序列而言，在进行基因辨识之前都应该把简单的大量的重复序列标记出来并除去，因为很多情况下重复序列会对预测程序产生很大的扰乱，尤其是涉及数据库搜索的程序。 2. 数据库搜索 把未知核酸序列作为查询序列，在数据库里搜索与之相似的已有序列是序列分析预测的有效手段。在理论课中已经专门介绍了序列比对和搜索的原理和技术。但值得注意的是，由相似性分析作出的结论可能导致错误的流传；有一定比例的序列很难在数据库里找到合适的同源伙伴。对于EST序列而言，序列搜索将是非常有效的预测手段。 3. 编码区统计特性分析 统计获得的经验说明，DNA中密码子的使用频率不是平均分布的，某些密码子会以较高的频率使用而另一些则较少出现。这样就使得编码区的序列呈现出可察觉的统计特异性，即所谓的“密码子偏好性”。利用这一特性对未知序列进行统计学分析可以发现编码区的粗略位置。这一类技术包括：双密码子计数(统计连续两个密码子的出现频率)；核苷酸周期性分析(分析同一个核苷酸在3,6,9,...位置上周期性出现的规律)；均一/复杂性分析(长同聚物的统计计数)；开放可读框架分析等。 4. 启动子分析 启动子是基因表达所必需的重要序列信号，识别出启动子对于基因辨识十分重要。有一些程序根据实验获得的转录因子结合特性来描述启动子的序列特征，
1、C值:一个单倍体基因组中DNA的总量. C值悖理 2、假基因:来源于功能基因但已失去活性的DNA 顺序. 3、遗传图 :采用遗传学分析方法将基因或其它DNA顺序标定在染色体上构建连锁图。遗传图距单位为厘摩(cM), 每单位厘摩定义为1%交换率。 4、物理图（Physical mapping）:采用分子生物学技术直接将DNA分子标记、基因或克隆标定在基因组实际位置。 5、重叠群:一群相互重叠的克隆或DNA顺序,可以是草图顺序或精确顺序(finished), 包括连续的(内部无间隙)或不连续的(内部含间隙)DNA顺序,未锚定到染色体上. 6、序列间隙：指测序时遗漏的序列，这些序列仍然保留在尚未挑选到的克隆中。 7、物理间隙：指构建基因组文库时被丢失的DNA序列，已从已有的克隆群体中永久性消失 8、全基因组鸟枪法测序：将基因组打成小片段后将其克隆到质粒载体中，然后随机挑取克隆对插入片段测序，并以获得的测序序列构建重叠群。在此基础上进一步搭建序列支架，最后以分子标记为向导将序列支架锚定到基因组整合图上。 9、支架(scaffold):一组已锚定在染色体上的重叠群, 内部含间隙或不含间隙. 10、作图测序: 按照大分子DNA克隆绘制的物理图分别在单个大分子DNA 克隆内部进行测序与序列组装,然后将彼此相连的大分子克隆按排列次序搭建支架,最后以分子标记为向导将搭建好的支架逐个锚定到基因组整合图上. 11、开放阅读框 ORF：指由一系列指令氨基酸的密码子组成，包括一个起始密码子（ATG），还有一个终止密码子（TAA，TAG，TGA） 12、基因敲除：将一段无关的DNA片段用来取代某一特定的基因。 13、同源性(homology):基因(序列同源性) 指起源于同一祖先但顺序已经发生变异的基因成员, 分布在不同物种间的同源基因又称直系基因. 同一物种的同源基因则称水平基因, 水平基因由重复后趋异产生. 14、一致性(identity):指同源DNA顺序的同一碱基位置的相同的碱基成员, 或者蛋白质的同一氨基酸位置的相同氨基酸成员, 可用百分比表示. 15、相似性(similarity):指同源蛋白质的氨基酸顺序中一致性氨基酸和可取代氨基酸所占的比例. 可取代氨基酸系指具有相同性质如极性氨基酸或非极性氨基酸的成员, 它们之间的代换不影响蛋白质(或酶)的生物学功能. 16、异染色质:深色区分布在细胞核的周缘，称为异染色质
实验一生物序列统计分析 一．实验目的 一般情况下，真核细胞中的线粒体是主要的能量生产中心。人类线粒体基因组在GenBank中的编号为“NC_001807”。以这条序列为例，学习有关DNA序列和蛋白质序列的统计分析方法。 1．学习和掌握在MATLAB平台上应用Bioinformatics工具包访问GenBank，并读取DNA序列。 2．学习和掌握在MATLAB平台上应用Bioinformatics工具包统计DNA序列的组成成分及含量，分析DNA序列的性质。 3．学习和掌握在MATLAB平台上应用Bioinformatics工具包搜索DNA序列的开放阅读框ORFs。 4．学习和掌握在MATLAB平台上应用Bioinformatics工具包，根据已定位的ORFs，实现DNA序列向蛋白质序列的转换。 5．学习和掌握在MATLAB平台上应用Bioinformatics工具包统计蛋白质序列中各种氨基酸含量。 二．实验内容 1．在MATLAB平台上应用Bioinformatics工具包访问GenBank，读取DNA序列。 ①用“web”命令在MATLAB平台上打开NCBI网页。 web('https://www.doczj.com/doc/6213834142.html,/') web('https://www.doczj.com/doc/6213834142.html,/genomes/framik.cgi?db=Genome&gi=12188') ②用“getgenbank”功能从GenBank中读序列信息到MARLAB mitochondria = getgenbank('NC_001807','SequenceOnly',true); 选项“SequenceOnly”使我们从GenBank中只读取“NC_001807”的序列信息。 “Mitochondria”是我们定义的变量，存在MATLAB的Workspace中。 ③查看变量mitochondria whos mitochondria 2．在MATLAB平台上应用Bioinformatics工具包统计DNA序列的组成成分及含量，分析DNA序列的性质。 ①查看DNA序列的性质 ntdensity(mitochondria)
什么是高通量测序？ 高通量测序技术（High-throughput sequencing，HTS）是对传统Sanger测序（称为一代测序技术）**性的改变, 一次对几十万到几百万条核酸分子进行序列测定, 因此在有些文献中称其为下一代测序技术(next generation sequencing，NGS )足见其划时代的改变, 同时高通量测序使得对一个物种的转录组和基因组进行细致全貌的分析成为可能, 所以又被称为深度测序(Deep sequencing)。 什么是Sanger法测序（一代测序） Sanger 法测序利用一种DNA聚合酶来延伸结合在待定序列模板上的引物。直到掺入一种链终止核苷酸为止。每一次序列测定由一套四个单独的反应构成，每个反应含有所有四种脱氧核苷酸三磷酸(dNTP)，并混入限量的一种不同的双脱氧核苷三磷酸(ddNTP)。由于ddNTP 缺乏延伸所需要的3-OH基团，使延长的寡聚核苷酸选择性地在G、A、T或C处终止。终止点由反应中相应的双脱氧而定。每一种dNTPs和ddNTPs的相对浓度可以调整，使反应得到一组长几百至几千碱基的链终止产物。它们具有共同的起始点，但终止在不同的的核苷酸上，可通过高分辨率变性凝胶电泳分离大小不同的片段，凝胶处理后可用X-光胶片放射自显影或非同位素标记进行检测。 什么是基因组重测序（Genome Re-sequencing） 全基因组重测序是对基因组序列已知的个体进行基因组测序，并在个体或群体水平上进行差异性分析的方法。随着基因组测序成本的不断降低，人类疾病的致病突变研究由外显子区域扩大到全基因组范围。通过构建不同长度的插入片段文库和短序列、双末端测序相结合的策略进行高通量测序，实现在全基因组水平上检测疾病关联的常见、低频、甚至是罕见的突变位点，以及结构变异等，具有重大的科研和产业价值。 什么是de novo测序 de novo测序也称为从头测序：其不需要任何现有的序列资料就可以对某个物种进行测序，利用生物信息学分析手段对序列进行拼接，组装，从而获得该物种的基因组图谱。获得一个物种的全基因组序列是加快对此物种了解的重要捷径。随着新一代测序技术的飞速发展，基因组测序所需的成本和时间较传统技术都大大降低，大规模基因组测序渐入佳境，基因组学研究也迎来新的发展契机和**性突破。利用新一代高通量、高效率测序技术以及强大的生物信息分析能力，可以高效、低成本地测定并分析所有生物的基因组序列。 什么是外显子测序（whole exon sequencing） 外显子组测序是指利用序列捕获技术将全基因组外显子区域DNA捕捉并富集后进行高通量测序的基因组分析方法。外显子测序相对于基因组重测序成本较低，对研究已知基因的SNP、
应用3'Full-RACE技术克隆PCDHB1基因3'非翻译区 廖柔霞：汕头大学医学院2008级本科一班 指导教师：许丽艳，方王楷，黄俏，吴健谊，吴炳礼，李恩民 【摘要】目的：克隆PCDHB1基因3'非翻译区，为深入研究该基因在食管癌等肿瘤或内分泌性疾病中的功能提供基本实验材料。方法：3'Full-RACE技术。结果：未能成功获得PCDHB1基因3'非翻译区克隆。结论：1）提示PCDHB1基因的结构组织形式可能十分复杂；2）今后，拟在具体方法学上改进，争取获得该基因3'非翻译区克隆。 PCDHB1（别名PCDH-BETA1），位于5q31。该基因共含有2534个碱基。这个基因是原钙黏连蛋白基因簇的一员，是5号染色体上的三个串连连接的成员之一。该基因组显示了与B的细胞和T细胞受体基因簇的一个不寻常的基因相似。其具体职能是未知的，但它很有可能在神经细胞间的连接和建立中发挥关键的作用。β钙黏连蛋白在特定的细胞粘连的事件扮演一个重要的角色。见诸于人类恶性肿瘤中的表达和调控功能,在很多情况下,导致肿瘤细胞浸润和转移的加剧。由于该基因的3＇端还未知，所以必须用3＇Full RACE的技术将PCDHB1的3＇端测出来。对此实验设计如下：应用3’RACE试剂盒及设计好的基因特异性引物扩增得到的基因片段，并将其连接到pEASYTM-Blunt Simple coloning vector 并测序。该基因片段由于出现“双峰”现象导致测序失败。 关键词：基因克隆，PCDHB1基因，3'非翻译区，3'Full-RACE The clone of 3' untranslational region of PCDHB1 gene using 3'Full-RACE RX Liao: 2008 Undergraduates of Medical College of Shantou University Supervisor: LY Xu, WK Fang, Q Huang, JY Wu, BL Wu, EM Li Abtract:Objective Methods Results Conclusion PCDHB1(PCDHB-BETA1),is located in 5q31. This gene is a member of the protocadherin beta gene cluster, one of three related gene clusters tandemly linked on chromosome five. The gene clusters demonstrate an unusual genomic or ganization similar to that of B-cell and T-cell receptor gene clusters. Their specific functions are unknown but they most likely play a critical role in the
DNA关键词： WG-BSA (全基因组重测序BSA) 对已有参考基因组序列的物种的所有作图群体（F1、F2、RIL、DH 和BC1等），对亲本进行个体重测序，对某个极端性状材料混池测序，检测SNP，获得与性状紧密关联的分子标记和精细定位区域，是目前最高效的基因定位方法。通过选取某个极端性状，利用高效率低成本的混池测序技术，勿需开发分子标记进行遗传图的构建，快速定位与性状相关的候选QTL。 MP-Reseq (多混池全基因组重测序) 针对特有的优良地方品种中的不同品种/品系，通过群体内pooling 建库的方法，进行全基因组重测序，采用生物信息学方法全基因组范围内扫描变异位点，能快速的定位不同混池样品基因组中明显经过人工或自然选择的区域，检测与性状相关的基因区域及其功能基因。 全基因组个体重测序 基于全基因组重测序的变异图谱通过测序手段结合生物信息分析研究同一物种不同个体之间的变异情况，获得大量的变异信息，如SNP、Indel、SV 等。主要可以快速地获得大量的分子标记以及不同个体在基因组水平上的差异。 全基因组关联分析-GWAS 通过重测序对动植物重要种质资源进行全基因组基因型鉴定，与关注的表型数据进行全基因组关联分析，找出与关注表型相关的SNP位点，定位数量性状基因，与数量性状相关的基因紧密连锁的SNP标记，后续可用于分子标记辅助育种，助力育种进程。 全基因组重测序-遗传进化 通过对来自全国各地、具有代表性的XX 份XX 材料进行全基因组重测序，检测SNP、Indel、SV，并利用获得的SNP 与SV 数据进行群体多样性分析，包括连锁不平衡分析、群体进化分析、群体结构分析、群体主成分分析等。 全基因组重测序-遗传图谱 基于全基因组重测序技术对已有参考基因组序列的物种进行个体或群体的全基因组测序,利用高性能计算平台和生物信息学方法，检测单核苷酸多态性位点（SNP），并计算多态性标记间的遗传连锁距离，绘制高密度的遗传图谱。通过与表型性状进行关联分析，利用获得的强关联性标记进行下游基因的精细定位。遗传图可用于分子标记辅助育种，重要性状候选基因克隆，辅助基因组组装，比较基因组学等研究。 细菌基因组de novo 测序 细菌是生物的主要类群之一，是所有生物中数量最多的一类。细菌广泛分布于土壤和水中，或者与其他生物共生，也有部分种类分布在极端环境中，例如温泉，甚至是放射性废弃物中。由于细菌自身的营
实验七核酸序列分析（附加部分） 1、发现核酸序列中的蛋白质编码区域。 1）利用NCBI ORF Finder。https://www.doczj.com/doc/6213834142.html,/gorf/gorf.html A、在NCBI上查找AC号为AE008569的核酸记录，思考：1、这个序列的名称？2、 这个序列所属的生物学分类？ B、进入OFR Finder，首先在页面下方的Genetic codes 下拉菜单中浏览现有的22种遗 传密码选择项（这里我们只使用默认的standard code），利用AC号或其裸序列（想 一想怎么能得到）进行ORF finding。 C、在结果显示页面中，按照序列的正向+1、+2、+3以及反向的-1、-2、-3进行的六框 翻译结果以图形的方式显示在页面中。利用默认的100bp阈值所发现的各框内的 ORF以绿色条状显示。同时，按照六框内所有发现的ORF的大小顺序，在页面的 右侧有一个列表，分别显示了ORF的翻译框、在基因组上的位置以及ORF的长度。 你可以改变ORF鉴别中的长度阈值（50，100，300），点击Redraw重新进行计算。 D、点击图形上的绿色条框，就可以对这个ORF进行检查（当然也可以点击右侧的ORF 列表），页面上会显示预测的氨基酸序列，同时页面上还嵌入了BLAST程序以及 NCBI的有关序列数据库以便于发现与此ORF相似的库记录。非常方便！ E、SixFrames是以另外一种方法计算并显示结果，点击SixFrames，结果中各框上边拉 下的绿色短线表示为一个起始密码子，而各框下方的粉色短线表示为一个终止密码 子。 F、如果你拥有一个高等生物的cDNA时，可以利用ORF finder这个简单的工具来找到 你的蛋白编码区域。因为cDNA不含有intron，因此可拥有与微生物相似的ORF结 构。 G、ORF finder可以正确地鉴定85%左右的蛋白编码区，但要发现一些很短的蛋白序 列，shadow gene或使用了非常用遗传密码子的基因，则需要使用那些包含了密码 子使用频率及使用偏好等统计学特性的程序，如GeneMark。这里给出两个GeneMark 网址：https://www.doczj.com/doc/6213834142.html,/GeneMark/ , https://www.doczj.com/doc/6213834142.html,/genemark/。2）发现真核生物基因组（如脊椎动物）序列中的蛋白质编码区域。 A、剪切位点（splice site）的预测。 脊椎动物的外显子很小（平均150bp），它们的剪切位点还有一定的变化。因此发现外显子要比利用ORF finder或GeneMark发现ORF困难得多。下面是一种外显子预测程序：MZEF。点击https://www.doczj.com/doc/6213834142.html,/，这是位于冷泉港实验室Michae Q. Zhang’s的主页，点击左侧的databases and Software Tools，进入的页面中包含了多个物种的启动子数据库、外显子发现工具等，点击页面中间的Gene –Finding (public)连接，则进入了MZEF页面（https://www.doczj.com/doc/6213834142.html,/tools/genefinder/）。程序的相关说明文件在页面下方的For more information about MZEF行的here链接中，事先阅读一下此文件，有助于程序的使用以及对输出结果的理解(https://www.doczj.com/doc/6213834142.html,/tools/genefinder/readme.htm )，你也可以阅读实验数据-实验七中的MZEFexample.PDF文件，这一文件也可以从Michae Q. Zhang’s的数据库及软件工具页面上找到(https://www.doczj.com/doc/6213834142.html,/reprints/mzefexample.pdf）。回到MZEF主页，点击Human 链接（https://www.doczj.com/doc/6213834142.html,/tools/genefinder/human.htm），进入由先前统计数据校准的人类编码外显子预测MZEF程序页面。 在NCBI上找到一条AC号为AF018429的人类核酸记录，这是一个包含了外显子1和外显子2的dUTPase基因（注意一下这两个外显子在基因上的位置）。将FASTA格式的序列粘贴到人类MZEF程序页面的检索框中，点击submit。程序很快给你返回结果。它发现
基因组信息学实验课 序列分析 第一部分：课堂知识复习 形式：题目问答，请同学单独回答以下问题，并陈述答题理由。 1．判断对错：两条序列的同源程度为60%。 2．判断对错：两条序列的相似性很高，所以它们一定是同源序列。错误3．计算下面两条序列的海明距离： AGCAAACACACTA ACATAAGCACACA 4．通过字符编辑操作将序列s转换成t s：AG-CA t：ACAC- 5．分析两条序列的关系时，（）方法可以通过观察矩阵对角线迅速发现可能的序列比对。 A Dot-plot B Pairwise-Alignment C BLAST D FASTA E Score Matrix 6．下列哪些是核酸序列数据库（）。 A GenBank B PDB 蛋白质结构数据库 C Entrez D EMBL-Bank E DDBJ 7．下列哪些是蛋白质序列数据库（）。 A PIR蛋白质信息资源 B SWISS-PROT C TrEMBL D EPD 真核生物启动子 E PDB大分子结构图 8．PDB文件的显示序列信息中，关键字（）作为显式序列标记，以该关键字打头的每一行都是关于序列的信息。 A HEADER B REMARK C SEQRES D EXPDTA 9．下列哪些是可以用来显示分子结构的软件（）。 A GCG-DS Visualizer B RasMol C ChemView D DSSP 10．Entrez数据库集成系统中集成了NCBI中哪些数据库中的信息（）。 A 核酸序列 B 蛋白质序列 C 生物大分子结构 D 基因组数据 E 生物 分类数据库 F 孟德尔人类遗传学数据（OMIM）G Pubmed
全基 1. 简 通过变（d 的功况，dise 比较 实验 （1）（2） 基因组重测序简介(Introduc 过高通量测序识deletioin, du 功能性进行综合杂合性缺失ease （cance 较基因组学，群验设计与样本 Case-Contr ）家庭成员组序数据分析 ction) 识别发现de plication 以及合分析；我们（LOH ）以及r ）genome 中群体遗传学综ol 对照组设计 组设计：父母novo 的som 及copy numb 们将分析基因及进化选择与中的mutation 综合层面上深计 ； -子女组（4 人matic 和germ ber variation 因功能（包括与mutation 之n 产生对应的深入探索疾病基人、3 人组或m line 突变，）以及SNP miRNA ），重之间的关系；以的易感机制和基因组和癌症多人）； 结构变异-SN 的座位；针对重组率（Rec 以及这些关系功能。我们将症基因组。 NV ，包括重排对重排突变和combination ）系将怎样使得 将在基因组学排突 SNP ）情在 学以及
初级数据分析 1．数据量产出：总碱基数量、Total Mapping Reads、Uniquely Mapping Reads统计，测序深度分析。 2．一致性序列组装：与参考基因组序列（Reference genome sequence）的比对分析，利用贝叶斯统计模型检测出每个碱基位点的最大可能性基因型，并组装出该个体基因组的一致序列。3．SNP检测及在基因组中的分布：提取全基因组中所有多态性位点，结合质量值、测序深度、重复性等因素作进一步的过滤筛选，最终得到可信度高的SNP数据集。并根据参考基因组信息对检测到的变异进行注释。 4．InDel检测及在基因组的分布: 在进行mapping的过程中，进行容gap的比对并检测可信的short InDel。在检测过程中，gap的长度为1~5个碱基。对于每个InDel的检测，至少需要3个Paired-End序列的支持。 5．Structure Variation检测及在基因组中的分布: 能够检测到的结构变异类型主要有：插入、缺失、复制、倒位、易位等。根据测序个体序列与参考基因组序列比对分析结果，检测全基因组水平的结构变异并对检测到的变异进行注释。 高级数据分析 1.测序短序列匹配（Read Mapping） （1）屏蔽掉Y染色体上假体染色体区域（pseudo-autosomal region）, 将Read与参考序列NCBI36进行匹配（包括所有染色体，未定位的contig，以及线粒体序列mtDNA（将用校正的剑桥参考序列做替代）)。采用标准序列匹配处理对原始序列文件进行基因组匹配， 将Read与参考基因组进行初始匹配；给出匹配的平均质量得分分布； （2）碱基质量得分的校准。我们采用碱基质量校准算法对每个Read中每个碱基的质量进行评分，并校准一些显著性误差，包括来自测序循环和双核苷酸结构导致的误差。 （3）测序误差率估计。 pseudoautosomal contigs，short repeat regions（包括segmental duplication，simple repeat sequence-通过tandem repeat识别算法识别）将被过滤； 2. SNP Calling 计算（SNP Calling） 我们可以采用整合多种SNP探测算法的结果，综合地，更准确地识别出SNP。通过对多种算法各自识别的SNP进行一致性分析，保留具有高度一致性的SNP作为最终SNP结果。这些具有高度一致性的SNP同时具有非常高的可信度。在分析中使用到的SNP识别算法包括基于贝叶斯和基因型似然值计算的方法，以及使用连锁不平衡LD或推断技术用于优化SNP识别检出的准确性。 统计SNV的等位基因频率在全基因组上的分布