序列比对（精选5篇）
序列比对 篇1
0 引言
1977年, Sanger测序法[1]的诞生是DNA测序技术的一个里程碑性质的大事件。在其后的三十多年中, 几乎所有的测序技术都只是Sanger测序法的改进, 而后研究人员又将Sanger测序法的研发推进到了自动化的层面, 从而大大提高了DNA序列的测定速度。在2004年, 454、SOLi D, Illumina等测序技术的兴起, 给序列的测定带来了飞跃式的变化, 但随着形态的多样化和应用的复杂化, 由于Sanger测序法的某些缺陷, 使得测序的通量和技术已经迟滞于该领域的发展需求, 而相对于Sanger测序法的新一代测序技术因其具有的高通量, 低能耗的优点, 使得新一代测序技术代替Sanger测序法, 而获得广泛的使用已成为势所必然。而由于新一代测序技术的产生, RNA序列的研究也随之发生了重大的改变。在此之前, RNA序列的测定主要是通过外显子芯片技术。外显子芯片可以用来测定RNA的序列信息, 也可以用来分析外显子表达量, 同时也能发现外显子的可变剪接等信息, 但是外显子组芯片的制作却需要丰富的先验知识, 并且与新一代测序技术相比, 外显子组芯片的花费是巨大的, 同时只能测定小范围内的序列, 不能对整个基因组实施全方位的分析。随着新一代测序技术的发展, 两种技术之间的差异也会越来越大。目前, 对于RNA序列的大部分研究已经转向了新一代测序技术产生的序列数据, 但即便如此, 外显子组芯片也依然在其中发挥着独特的重要作用。新一代测序技术的产生给RNA序列的分析技术也带来了重大的改变。新一代测序技术产生的数据序列短、覆盖度高, 但数据量大, 这给传统的RNA序列分析工具设置了难题, 因而应运而生地出现了多种基于新一代测序技术的RNA序列分析工具。相对于传统Sanger测序法, 新一代测序数据产生的序列较短, 通常称为短序列 (reads) , 但是新一代测序数据产生的数据量却要远远大于Sanger测序法。必须正视这一问题的积极解决, 才能确保新一代测序技术的先进性和有效性得以充分的发挥。
1 DNA序列比对工具现状
对于当前的RNA序列比对工具的研究, 首先就要研究DNA序列比对工具, 因为当前的RNA序列比对工具都是以DNA序列比对工具为基础发展得来的。
新一代测序技术产生后, 曾经应用于外显子芯片技术的RNA序列分析方法已经不再适用, 但是这些方法却可留下许多有益的启发。新一代测序技术的产生, 给序列比对也带来了很大的挑战, 人们都致力于研发更为有效的DNA序列比对软件。众所周知, 只有找到新的、性能更佳的DNA序列比对方法, 才能使高通量数据问题获得理想的解决。而RNA序列比对工具就是根据DNA序列比对软件工具, 在其基础之上并根据RNA的不同性质和各种分析需求, 构造可用于RNA序列的分析工具和分析策略。
基于新一代测序技术设计了很多DNA序列比对工具。由于建立索引的不同, 目前DNA序列比对工具主要分为两类, 一类是用Hash表来构建索引, 另一类是用BWT (Burrows-Wheeler Transform) 来建立索引结构, 该索引结构由于占用空间小、搜索速度快等优点正被广泛地关注和使用[2]。
在高通量的序列比对中, 索引是一个非常有效的机制。通过构建索引可以提高检索速度, 从而提高了整体比对速度。基于Hash表的索引构建可以分为两种。
一种是将参考序列 (reference sequence) 构建成Hash表索引, 建立索引时根据所需的短序列特性, 例如长度等信息, 将原始的参考序列分成连续重叠的短序列, 根据不同的Hash算法将这些短序列存储起来, 然后将实验得到的短序列与参考序列生成的Hash表进行比较, 从而确定短序列的比对位置。基于Hash表的全部索引结构比对工具都可以比对有插入删除的序列, 但是时间和空间的开销却很可观。
另一种是将短序列 (reads) 数据构建成Hash表索引, 这种序列比对工具却较少。
还有的软件两种方式都采用以提高比对速度。基于Hash表的索引软件主要有Blast、Eland、MAQ、Bfast等。其中, Blast是出现最早的基于Hash表的索引软件, 目前有很多学者正致力于减少基于Hash表的索引比对算法所需花费的时间和占用的空间。
基于BWT索引结构的DNA序列比对软件在目前的学术界较为流行。BWT变化方式比较复杂, 在这里就不多做介绍了, 但需要知道的是, 该方式占用空间小, 比对速度快。基于BWT索引结构的DNA序列比对软件也自然会有其无法忽视的弱点, 即在处理插入删除上显然没有基于Hash表的DNA比对软件有效。当基于BWT索引结构的DNA序列比对每增加一个插入删除位点, 就会大大增加比对负担, 并且截至目前为止, 也没有找到这个问题的合理解决方式。但是在不允许插入删除的比对中, 人们还是更为倾向于选择基于BWT索引结构的DNA序列比对软件。基于BWT索引结构的DNA序列比对软件中, 最具有代表性的是Bowtie和BWA。其中, Bowtie不接受插入删除, 只处理失配位点, 所以速度更快一些;而BWA却可允许少量的插入删除, 速度相对来说就会慢一些, 这主要是由处理插入删除时消耗较多资源而引起的。
2 RNA序列比对工具分析
对于RNA序列的研究并不能完全等同于DNA序列, 主要是由于RNA序列是由不连续的片段组合而成, 这种不连续的片段就叫做外显子 (exon) 。RNA序列虽然是以D-NA序列为模板转录而来, 但是与DNA序列又有很大的不同, 因为不是所有的DNA序列都会出现在成熟的RNA中, 并且最后翻译成蛋白质, 也只有外显子才能获得这种表达。RNA序列的转录及翻译过程如图1所示。初始转录成的RNA会经过一系列的生物活动, 剪接掉内含子, 保留外显子, 并将外显子连接在一起, 同时在5’端加上一个帽子, 3’端加上一个多聚腺苷的尾巴, 最后还要经过一系列的修饰, 才能转运到细胞核外, 翻译成蛋白质。可想而知, 将RNA序列直接比对到DNA参考基因组上, 将会产生很大的问题, 所以需要设计适用于RNA的序列比对策略。
在DNA序列比对软件基础上, 根据不同的需求, 产生了很多RNA序列分析工具。其中, 以Top Hat[3]的应用最为广泛, Top Hat是建立在Bowtie发展之上的, 速度快, 占用空间小, 但是同样也具有不允许插入删除的缺陷。Top Hat首先利用Bowtie将所有的短序列比对到参考基因组上, 然后将比对上的短序列连接成外显子区域, 再将外显子区域外延几个bp的长度, 并参考已知的外显子剪接组合, 试用外显子区域上的不同组合, 将Bowtie在第一轮没有获得比对成功的短序列继续比对至组合而成的参考序列上, 如果确有短序列实现了这种有效比对, 就认定这种组合是正确的。
由于Bowtie是基于BWT索引的DNA比对工具中最早研发成功的, 所以后续研究开展得较为充分, 配套工具又很丰富, 知名度也相对较高, 所以使用选择者也就较多。同样, Top Hat的开发时间也是目前较为有效的几种RNA序列分析工具中位居首位的, 因而也成为当前流传甚广的分析工具。即使后面推出了更多的RNA序列分析工具, 研究学者们也依然重点关注Bowtie和Top Hat。Top Hat的开发带动了RNA序列比对软件在新一代测序技术上的策略改变, 由原来的主要依靠分析来解决RNA的比对, 转变为依靠序列本身的信息来解决RNA序列的比对。在前文提到了RNA在转录后需要经过修饰, 在很大程度上与原始的DNA序列已经有所不同, 又由于Top Hat不能处理插入删除, 可想而知在RNA序列比对上, Top Hat还是存在着一些问题。而且在一定限度上, Top Hat还需要依靠RNA序列的先验知识, 所以在寻找未发现的外显子上面, 效果不是很好。
在Top Hat之后, 相继又产生了其他的RNA序列比对工具, 例如Map Splice[4]、Splice Map[5]等, 这些工具也是建立在DNA比对工具基础之上, 构造出的适用于m RNA的序列比对软件。这些软件中, Map Splice是利用Bowtie来进行短序列比对, 但是Map Splice的比对策略却与Top Hat存在着不同。首先, Map Splice将实验测得的短序列分成连续不重叠的小片段, 将小片段比对到参考基因组上, 再利用小片段之间的联系, 找出外显子所在位置;然后, 利用统计学特性最终确定外显子位置以及外显子边缘。Splice Map则主要是应用Eland。首先, Splice Map将短序列分成重叠的50bp长的小片段, 将小片段的两端25bp长的序列比对到参考基因组上, 而后根据两端序列比对情况再分析外显子区域, Splice Map在时间和空间上与Top Hat和Map Splice都要偏长、偏大, 并且准确性还较低。另外, 也还有很多其他的比对工具, 但应用却较少, 诸如Split Seek[6], ABMapper[7]等。
此外, 还有一些RNA序列分析策略, 虽然没有产生新的算法来解决新一代测序技术之下的RNA序列比对的问题, 但是通过组合现有的DNA序列比对方法, 产生了一个有效的RNA序列比对流程, 使得RNA序列的比对结果更为精确。例如RUM[8]和RNA-MATE[9]。RUM不但包括序列比对流程, 同时还包括一个RNA序列模拟生成器。RUM首先利用Bowtie将序列比对到参考基因组和转录组上, 将剩余没有得到比对的序列运用BLAT再次进行比对。但是RUM却需要依靠现有的转录库来分析序列, 在RNA序列分析上表现了很大的局限性, 如不能发现新的外显子以及可变剪接组合信息, 在新功能的发现上应用空间也不大。RNA-MATE允许使用任何比对软件。首先将所有的短序列比对到参考基因组上, 再将剩余的序列分割成较短的序列, 进行比对, 如此循环反复, 直至达到一个设定的限度停止。
还有一类m RNA的分析策略是, 首先, 将测序得到的所有短序列 (reads) 利用新一代测序技术的组装工具拼装到一起, 形成长的contigs, 再利用简单的DNA序列比对工具, 就可以将RNA序列比对到参考基因组上, 而不是只能应用基于新一代测序技术的DNA序列比对工具, 才可以解决RNA序列中的比对不连续问题。照此举例即如Trans-ABy SS[10]。
RNA的种类繁多, 其生物学特性也为数众多, 所以基于不同性质的各类分析工具也一定会有很多。目前研究更多地集中在基因融合方面, 代表性的有shortfuse、Fusion Map和Tophat Fusion等。这些研究都是利用pair-end序列数据相对位置的改变以分析得出基因位置的相对变化, 因而发现癌细胞中的基因融合现象。在小RNA的序列比对方面还有Micro Razer S等, 小RNA由于序列较短, 比对将更加困难。
综上所述, 对现有的基于新一代测序数据的RNA序列比对算法做以总结, 结果如表1所示。
3 RNA序列分析工具分析
除了RNA序列比对软件外, 还有一些比较著名的RNA序列分析软件, 例如:Cufflinks[11]和Scripture[12]。这两种工具都是首先利用Top Hat进行RNA序列比对, 然后通过各自的分析策略, 来推断isoform的工具。通过将RNA的可变剪接清楚地呈现在人们面前, 使得m RNA序列分析在整体上具备了完备性。这也是Top Hat之所以受到欢迎的另一个原因。
首先, Cufflinks可使用任何版本的Top Hat, 将所有的pair-end序列数据比对到参考基因组上, 然后利用组装算法, 将互有交叠的pair-end序列组装到一起, 同时依据pair-end序列的交叠信息发现不同的组装路径。而后, 再根据每个位置上的序列覆盖度, 运用统计学的方法分析出每种isoform的比例。对全基因组中的每个可变位置都计算该比例, 最后分析得出整体的isoform分布。
Scripture和Cufflink都是利用Top Hat将pair-end序列比对到参考基因组上, 根据pair-end序列数据的相对位置信息, 将可能的外显子组合寻找出来。两者不同的是, Cufflink以pair-end序列为节点构造出连通图, 而Scripture则是以每个碱基为节点构造连通图。Scripture首先列举出参考序列上的碱基, 在参考序列上相邻的碱基之间有一条边, 在比对序列上相邻的碱基之间也有一条边, 最后形成连通图。Scripture同时利用这样的方法在统计学上排除了错误剪接位点, 重新确定了外显子边界, 从而根据序列的覆盖信息来确定isoform的组份。
Cufflinks和Scripture的分析结果都可以利用基因组浏览器进行观看, 更直观地反映转录组信息。除了Cufflinks和Scripture之外, 又新近涌现了一些RNA序列分析工具, 例如:FDM[13]。
4 结束语
RNA序列的分子在遗传上具有重要的应用, 在疾病的发现和治疗上也表现出了非同寻常的意义。RNA生物学性质的多种多样又给RNA序列的比对分析带来了巨大的困难, 目前的RNA序列比对软件依然无法满足已有需求, 因而在RNA序列的研究和分析上, 依然任重而道远。
生物信息学中的序列比对算法 篇2
生物信息学是80年代末随着人类基因组计划的启动而兴起的一门新的交叉学科,最初常被称为基因组信息学。生物信息学是在生命科学的研究中,以计算机为工具对生物信息进行储存、检索和分析的科学。它是当今生命科学和自然科学的重大前沿领域之一,同时也将是21世纪自然科学的核心领域之一。其研究重点主要体现在基因组学和蛋白组学两方面,具体说,是从核酸和蛋白质序列出发,分析序列中表达结构与功能的生物信息。
生物信息学的研究重点主要体现在基因组学和蛋白质学两方面,具体地说就是从核酸和蛋白质序列出发,分析序列中表达结构和功能的生物信息。生物信息学的基本任务是对各种生物分析序列进行分析,也就是研究新的计算机方法,从大量的序列信息中获取基因结构、功能和进化等知识。在从事分子生物学研究的几乎所有实验室中,对所获得的生物序列进行生物信息学分析已经成为下一步实验之前的一个标准操作。而在序列分析中,将未知序列同已知序列进行相似性比较是一种强有力的研究手段,从序列的片段测定,拼接,基因的表达分析,到RNA和蛋白质的结构功能预测,物种亲缘树的构建都需要进行生物分子序列的相似性比较。例如,有关病毒癌基因与细胞癌基因关系的研究,免疫分子相互识别与作用机制的研究,就大量采用了这类比较分析方法。这种相似性比较分析方法就称为系列比对(Sequence Alignment)。目前,国际互联网上提供了众多的序列比对分析软件。然而,不同的分析软件会得到不同的结果,同时所使用的参数在很大程度上影响到分析的结果。有时常常会由于采用了不合适的参数而丢失了弱的但却具有统计学显著性意义的主要信息,导致随后的实验研究走弯路。因此,生物信息学中的序列比对算法的研究具有非常重要的理论与实践意义。
序列比对问题根据同时进行比对的序列数目分为双序列比对和多序列比对。双序列比对有比较成熟的动态规划算法,而多序列比对目前还没有快速而又十分有效的方法。一般来说,评价生物序列比对算法的标准有两个:一为算法的运算速度,二为获得最佳比对结果的敏感性或准确性。人们虽已提出众多的多序列比对算法,但由于问题自身的计算复杂性,它还尚未得到彻底解决,是生物信息学中一个非常重要且具有挑战性的研究课题。
2 序列比对
比较是科学研究中最常见的方法,通过将研究对象相互比较来寻找对象可能具备的特性。在生物信息学研究中,比对是最常用的研究手段。
最常见的比对是蛋白质序列之间或核酸序列之间的两两比对,通过比较两个序列之间的相似区域和保守性位点,寻找二者可能的分子进化关系。进一步的比对是将多个蛋白质或核酸同时进行比较,寻找这些有进化关系的序列之间共同的保守区域、位点和profile,从而探索导致它们产生共同功能的序列模式。此外,还可以把蛋白质序列与核酸序列相比来探索核酸序列可能的表达框架;把蛋白质序列与具有三维结构信息的蛋白质相比较,从而获得蛋白质折叠类型的信息。
序列比对的理论基础是进化学说,如果两个序列之间具有足够的相似性,就推测二者可能有共同的进化祖先,经过序列内残基的替换、残基或序列片段的缺失、以及序列重组等遗传变异过程分别演化而来。
早期的序列比对是全局的序列比较,但由于蛋白质具有的模块性质,可能由于外显子的交换而产生新蛋白质,因此局部比对会更加合理。通常用打分矩阵描述序列两两比对,两条序列分别作为矩阵的两维,矩阵点记录两个维上对应的两个残基的相似性分数,分数越高则说明两个残基越相似。因此,序列比对问题变成在矩阵里寻找最佳比对路径,目前最有效的方法是NeedlemanWunsch动态规划算法,在此基础上又改良产生了Smith-Waterman算法。
在进行序列两两比对时,有两方面问题直接影响相似性分值:取代矩阵和空位罚分。粗糙的比对方法仅仅用相同/不同来描述两个残基的关系,显然这种方法无法描述残基取代对结构和功能的不同影响效果。用一个取代矩阵来描述氨基酸残基两两取代的分值会大大提高比对的敏感性和生物学意义。虽然针对不同的研究目标和对象应该构建适宜的替换矩阵,但国际上常用的替换矩阵有PAM和BLOSUM等。它们来源于不同的构建方法和不同的参数选择。对于不同的对象可以采用不同的替换矩阵以获得更多信息。
多序列比对就是把两条以上可能有系统进化关系的序列进行比对的方法。目前对多序列比对的研究还在不断前进中,现有的大多数算法都基于渐进的比对思想,在两两比对的基础上逐步得到多序列比对的结果。
多序列比对算法是生物信息学中的最基本算法,是生物体的进化分析、蛋白质的分析和预测等生物体研究的基础,具有重要的理论意义和使用价值。
3 序列同源性与序列相似性
序列相似和序列同源是不同的概念,序列之间的相似程度是可以量化的参数,而序列是否同源需要有进化事实的验证。序列同源(homology)指的是序列来自相同的祖先,意味着这些序列具有相同的进化历史,而序列的相似性(similarity)指的是两序列在某参数条件下的相像,它可以用相同残基的百分比或是其他的方法来表示。序列之间的相似度是可以量化的参数,而序列是否同源需要有进化事实的验证,显著的相似性通常意味着同源。
序列比对是运用某种特定的数学模型或算法,找出两个或多个序列之间的最大匹配碱基或残基数,比对算法的结果在很大程度上反映了序列之间的相似性程度以及它们的生物学特征。序列比对根据同时进行比对的序列数目多少可分为双序列比对(pairwise sequence alignment)和多序列比对(multiple sequence alinment)。序列比对从比对范围考虑也可分为全局比对(global alignment)和局部比对(local alignment),全局比对考虑序列的全局相似性,局部比对考虑序列片断之间的相似性。如下所示。
全局比对:
在实际应用中,用全局比对方法企图找出只有局部相似性的两个序列之间的关系显然是徒劳的;而用局部比对得到的局部相似性结果则同样不能说明这两个序列的三维结构或折叠方式是否相同。
4 序列比对算法
在生物分子信息处理过程中,将生物分子序列抽象为字符串,其中的字符取自特定的字母表。字母表是一组符号或字符,字母表中的元素组成序列。如DNA序列由四种核苷酸组成,用“A”,“T”,“C”,“G”代表四种碱基,其复杂度为4,“CCATGCTAGAT”可代表一个简单的DNA序列。蛋白质序列由20中氨基酸组成,由{ABCDEFGHIKLMNPQRSTV WXYZ}代表不同的残基。“X”表示某个不确定的残基。“B”表示天冬胺或天冬胺酸,用三个字符表示“Asx”。“Z”表示谷氨酰胺或谷氨酸,用三个字符表示为“Glx”,其复杂度为23,“BEGSSTTNMABNNMA”可代表一个简单的蛋白质序列。因此生物序列比对可以看作字符串的比对。对字符串的编辑操作有以下三种:插入———在序列中插入一个或多个字符;删除———在序列中删除一个或多个字符;替换———用另一个字符替代某个字符。
4.1 序列比对基本定义
定义1序列是有限长度的字符串,序列中的字符由某个有限字符集合Ω确定。对于DNA,Ω={A,C,T,G}。对于蛋白质,Ω由20种代表氨基酸的字符组成。
定义2对于序列S,|S|表示S中字符个数。S[i]表示序列的第i个字符。S[1…i]表示序列的前i个字符组成的子序列。
定义3我们用“-”来表示插入和删除所产生的空位,则:
(1)(a,a)表示匹配(从序列S到序列T没有发生变化);
(2)(a,-)表示从S中删除字符a,或是在T中插入空位;
(3)(a,b)表示用T中的字符b替代S中的a,(a≠b);
(4)(-,b)表示在S中插入空位,是从T中删除字符b。
定义4对于x,y∈Ω∪{-},定义σ(x,y)为计分函数,表示x,y比较时的得分。以下是最简单的一种定义公式:
定义5 S和T的一个比对A用序列S和T中字符的一一对应表示,其中
(1)|S'|=|T'|;
(2)S',T'去掉空格就是S和T。
定义6序列比对A的得分为M,得分M越高表示序列的相似程度越高。
4.2 序列比对算法
Needleman-Wunsch算法是双序列比对的经典算法,其使用的是动态规划的基本思想。对于长度分别为m和n的两个序列S和T,构造矩阵T,矩阵T中的最后一个元素T[m][n]即对应于最优比对的得分,而最优比对本身则可以通过回溯算法得到。该算法的时间和空间复杂度均为O(mn)。
Smith-Waterman对Needleman-Wunsch算法稍加改动,使其可以计算局部最优比对,其所需的时间和空间复杂度仍是O(mn)。
Mayers和Miller使用Hirschberg提出的技巧在时间复杂度不变的前提下将空间约减到O(m+n)。
M.Crochemore等人对经典算法加以改进,提出了一个可以在O(n2/log n)时间内实现的双序列比对算法。其主要思路是对序列进行压缩编码,从而将序列分为若干段,从而将比对所构造的矩阵分为若干块来计算。后面的块的计算可以利用前面的块的结果在常数时间内计算得出。
除了利用矩阵来计算序列比对外,还有两种常用于序列分析的后缀队列Suffix Array和后缀树。
AVID是一个双序列全局比对算法,首先,用后缀树找出所有的最大匹配子序列,并在其中选择所有不重叠,不交叉的序列作为锚点。然后用锚点作为最后比对的一部分,在锚点之间的序列部分则递归的用此算法进行比对。
生物数据的信息量极大,序列比对的计算需要耗费大量的时间。由于进行算法可以大大地加快问题求解速度,近年来对该问题并行化的研究也引起研究者的注意。
在CREW-PARM模型上,Aggarwal和Apostolico等人独立地提出了一个O(log m log n)时间,使用mn/log m个处理机的并行算法;Mi Lu等人设计了两个并行算法;一个使用mn/logm台处理机,时间复杂度为O(log2m+log m);另一个使用mn/log2mlog log m台处理机,时间复杂度为O(log2m log log m)。
对于多序列比对问题,传统方法所采用的表示模型是行一列模型,即对于输入的多个序列插入空位并排列比对,使其达到相同的长度。对于N个序列S1,…,Sn,其多序列比对是一个新的序列集S'=(S1',…,SN'),S'的所有序列长度相同,并且每一个序列Si'由Si插入空位‘-’得到。如果将各序列在垂直方向排列起来,则可以根据每一列观察各序列中字符的对应关系。图1是6个蛋白质序列片断基于行-列模型的多序列比对。
多序列比对问题实际上是两条序列比对问题的一般化推广。但是由于DNA或蛋白质数据库容量的指数级增长,当比对的序列大大超过两个时,基于基本动态规划法的多序列比对算法的计算量是非常惊人的,这使得多序列比对这一NP难题变得更加复杂。因此,为了解决这一问题,许多近似算法和启发式算法被提出。以下介绍几种典型的多序列比对算法。
动态规划方法:给定k条长度均为n的序列,根据在两条序列比对中的动态规划算法的思想,需要计算一个K维的超级立方体,该立方体的尺寸为(n+1)k。在双序列比对的动态规划解决方案中,每一项(i,j)要由(i-1,j-1)、(i-1,j)和(i,j-1)这三项来决定,在这个超级立方体中的每一项要有2k-1个相邻的项来决定。这样该问题的时间复杂度是O((2n)k),空间复杂度是O(2nk)。
渐进比对算法:渐进比对算法是最常用的、简单而又有效的启发式多序列比对方法,它所需要的时间较短、所占内存较少。这个算法首先是Hogeweg和Hesper给出的,随后Feng和Doolittle对此做了进一步研究和改进。基于渐进比对算法并被广泛使用且成为多序列比对标准方法的软件有:Clustal W和T-Coffee等。渐进比对算法的基本思想是迭代地利用两序列动态规划比对算法,先由两个序列的比对开始,逐渐添加新序列,直到所有序列都加入为止。但是不同的添加顺序会产生不同的比对结果。因此,确定合适的比对顺序是渐进比对算法的一个关键问题。而两个序列越相似人们对它们的比对就越有信心。因此,整个序列的比对应该从最相似的两个序列开始,由近至远逐步完成。作为全局多序列比对的渐进比对算法有个基本的前提假设:所有要比对的序列是同源的,即由共同的祖先序列经过一系列的突变积累,并经自然选择遗传下来的。分化越晚的序列之间相似程度就越高。因此,在渐进比对过程中,应该对近期的进化事件比远期的进化事件给予更大的关注。由于同源序列是进化相关的,因此可以按着序列的进化顺序,即沿着系统发育树(指导树)的分支,由近至远将序列或已比对序列按双重比对算法逐步进行比对,重复这一过程直到所有序列都已添加到这个比对中为止。
渐进比对算法主要由三个步骤组成:计算距离矩阵;构建指导树;依据指导树进行渐进比对。
这类算法的主要优点是:简单、快速,但存在两个主要问题:比对参数选择问题和局域最小化问题。
迭代比对方法:这种方法是使用比对记分函数反复添加一附加的序列到已比对的比对序列中,首先在所有的两条序列比对中找出距离值最小的一组,组成最优比对,然后反复地找出与最优比对距离值最小的序列。与最优比对的表头文件进行匹配,并且根据所得的结果相应的修改最优比对和表头文件。
Clustal W算法:比对过程中,先对所有的序列进行两两比对并计算它们的相似性分数值,然后根据相似性分数值将他们分成若干个组,并在每组之间进行比对,计算相似性分数值。根据相似性分数值继续分组比对,直到得到最终比对结果。
当得到多序列比对后,需要对比对的质量进行评价,SP模型是评价比对优劣的最常用模型。设得分函数具有可加性,多序列比对的得分是各列得分之和,对于某一列字符的得分可用公式(3)进行计算,即某一列字符的SP得分为一列中所有字符对得分之和:
其中ci表示该列中的第i个字符,f(ci,cj)表示字符ci和字符cj比较所得分值。具体计算时,可以先对多序列比对的每一列进行计算,然后将各列得分相加,也可以先计算所有两两序列比对的得分,然后再将得分相加。这两种计算在f('-','-')=0这一条件成立下等价。
多序列比对的目标是:在计分机制确定的情况下,寻找使得比对得分最高的多序列之间的最优比对。可以证明,利用SP模型寻找最优多重序列比对是一个NP完全问题。
要获得给定的多个基因或蛋白质序列之间的一个正确的比对是一个困难的计算问题,其困难在于两个方面:一是如何根据包括结构信息在内的生物学意义对给定比对打分,即如何获得一个完美的目标函数(Obj ective Function简称OF);二是在目标函数确定的情况下,如何求得分值最高的最优比对。前者要依据生物学的知识和实际问题的需要来决定。假设已经求得的目标函数相当完美且简单,后者也将是一个非常困难的计算问题。
5 结束语
随着生物学数据的大量积累,对序列比对算法的敏感性和运算速度提出了更高的要求,对计算的挑战就令人生畏,序列比对中的主要困难就是如何研究和设计同时具备高敏感性和高速度的算法,序列比对算法研究仍然是生物信息学中一个非常重要且具有挑战性的研究课题,对序列比对算法研究具有非常重要的意义。
参考文献
[1]Katoh.K,Kuma.K,Toh.H.,and Miyata.T.MAFFT version5:improvement in accuracy of multiple sequence alignment[J].Nucleic Acids Research.2005,33(2):511-518.
[2]Morgenstern,B.Werner,N.,Prohaska,S.J.,Steinkamp,R.,Schneider,I.,Subramanian,A.R.,Stadler,P.F.,and Weyer-Menkhoff,J.Multiple sequence alignment with user_defined constraints[J].Bioinformatics.2004.
[3]Simossis.V.A,Kleinjung.J and Heringa.J.Hommology-extended sequence alignment[J].Nucleic Acids Research,2005,33(3):816-824.
[4]Zhang,M.,Fang,W.W.,Zhang,J,H.,and Chi,Z.X.MSAID:Multiple Sequence Alignment Based on a Measure of Information Discrepancy[J].Computational Biology and Chemistry,2005,29(2):175-181.
[5]Edgar,R.C.MUSCLE:multiple sequence alignment with high accuracy and high throughput[J].Nucleic Acids Research,2004,32(5):1792-1797.
[6]Edbert,R.C.,and Sjolander,K.COACH:profile-profile alignment of protein families using hidden Markov models[J].Bioinformatics,2004,20(8):1309-1318.
[7]T K Attwood,D J Parry-Smith著.罗静初,等译.生物信息学概论[M].北京:北京大学出版社,2001:141-145.
[8]张敏.生物序列比对算法研究现状与展望[J].大连大学学报,2004,25(4):75-78.
[9]张永,李其申,江泽涛,蔡虹.基于序列结构信息的多序列比对算法[J].微计算机信息,2007,23(21):240-242.
序列比对 篇3
Needleman算法[2]是经典的全局动态规化比对算法，Smith-Waterman算法[3]是对Needleman算法进行改进而形成的局部动态规划比对算法，适用范围很广。BLAST算法[4]是启发式的局部相似性搜索算法，采用搜索种子-扩展种子的策略，可以快速给出数据库中与查询序列有较高相似度的序列集合，但该算法的序列比对精度稍差。BLAT算法[5]与BLAST相似，不同之处在于BLAT对数据库建立索引而非对查询序列建立索引;BLAT将高得分片段对接合在一起作为比对结果返回，而BLAST将每个高得分片段对都返回。国防科大的张阳等用FP-GA构成的脉动阵列对Smith-Waterman算法进行了优化[6]，取得了不错的进展。
现场可编程门阵列(FPGA)内部支持程序并行执行，时钟周期在ns级，很适合用来进行序列比对。遗憾的是，FPGA的存储容量一般不大。目前，Altera公司的新一代FPGA产品Arria II GX，其Memory大小不到1 MB。作者所在小组设计的布尔逻辑DNA序列比对系统采用PC内存存储-FPGA比对的策略，弥补FPGA容量不足的问题。这种策略使得PC和FPGA之间的数据传输成为比对系统的关键。USB 2.0最高速度达480 Mb/s，可以满足FP-GA对数据的处理速率，作为比对系统的数据传输方式在理论上是可行的。
本文论述了布尔逻辑DNA序列比对系统软件平台的设计和实现。
1 布尔逻辑比对系统组成
布尔逻辑比对系统由硬件(协处理器)和软件(BooleanAlign)两部分组成。BooleanAlign配合协处理器完成整个系统的功能，二者之间通过USB系统互连。协处理器是一个外接USB设备，通过主机上的USB接口与主机相连。系统的组成框图如图1所示。
当开始比对的时候，BooleanAlign首先向协处理器发送ASCII码表示的DNA序列。由于DNA序列往往很长，BooleanAlign与协处理器协商好分段发送，每次发送512 B。每发送一段数据，协处理器就处理一段数据，然后返回处理结果。接着，BooleanAlign继续向协处理器发送下一段数据。这样循环，直到序列全部发送完毕。序列发送完毕的同时，比对结果也全部返回到BooleanAlign。最后，BooleanAlign依据这些数据构建两条序列的最佳匹配和系统发生树。比对系统的工作流程图如图2所示。
2 BooleanAlign平台设计
BooleanAlign是布尔逻辑比对系统的总控平台。DNA序列的输入、比对的开启、数据的发送与接收、控制命令的发送以及最佳匹配的构建和系统发生树[7]的重建，都由BooleanAlign控制和完成。BooleanAlign的运行环境为Windows，采用Microsoft公司的高效编程语言C++加以实现。BooleanAlign的程序调用关系如图3所示。
应用程序层(D2XX API Application)通过调用动态链接库(FTD2XX.DLL)函数实现与USB设备(FTD2XX.SYS)的通信。FTD2XX.SYS和FTD2XX.DLL分别是FTDI(Future Technology Devices International Ltd)公司为其USB接口转换芯片FT245BM提供的驱动程序和编程接口函数库。BooleanAlign平台的设计包括人机界面设计、数据通信、比对结果处理和基于距离的系统发生树重建。
2.1 人机界面设计
利用MFC(Microsoft Foundation Class)可以快速构建Windows应用程序框架[8]。VB也可以快速构建应用程序框架，但执行效率不如VC++.Java源代码经编译后变成字节码(.class)文件，程序需要边解释边执行，不适应BooleanAlign程序特性。基于以上情况，BooleanAlign选择MFC作为人机界面的实现工具。BooleanAlign人机界面采用MFC多文档结构，以满足用户同时打开多个窗口查看比对数据的要求。采用“多视窗”技术，创建动态分割窗口，将界面分割成独立的多个功能区域。菜单命令的响应函数，用来实现用户的各种操作。实现过程中的难点在于:视图窗格(pane)中显示的四种碱基字符(A、T、C、G)的字符宽度不相等，导致鼠标移动到相应字符上时，字符位置(posi)显示错误。解决方法是:调用CDC的文本绘图函数ExtTextOut()，将每一个字符的实际显示宽度调整成一致。这样，程序就可以直接利用鼠标的逻辑坐标值来计算字符在碱基序列中的位置了。字符位置计算公式如式(1)。
式(1)中，posi是字符在碱基序列中的位置，x_LP-Point是鼠标的逻辑坐标，widthCharacterRect是每个字符的实际显示宽度。人机界面如图4所示。
2.2 数据通信
2.2.1 数据通信格式
上层应用程序与底层硬件的数据交换必须按照一定的协议进行。BooleanAlign发送给协处理器的数据和协处理器返回的数据必须具有一定的格式，双方才能识别所接收到的数据的意义。
BooleanAlign发送给协处理器的数据单元的格式如图5。
“数据块”是长度为k个字节的DNA序列，即k个碱基。k在一般情况下都取512，只有在最后一个数据单元k的值可能会小于512。由于插入和缺失是一对相对的概念，为了区别这两个概念是相对于哪条DNA序列而言，有必要对发送的数据加个标识。“标识”给数据块加上了一个标志，用来说明该数据块来自两条序列中的哪一条，占一个字节，取值范围是三个ASCII字符“E”、“N”、“D”，如表1所示。当标识是“D”时，表示参与比对的两条序列已全部发送到协处理器，即数据发送完毕。
协处理器返回给BooleanAlign的数据单元的格式如图6所示。
cnt代表发生变异的碱基位点(碱基字符在DNA序列中的位置)个数，占一个字节，最大可表示256个变异位点。每个变异位点信息用3个字节描述，position描述变异发生在DNA序列的什么位置，占两个字节;type描述变异是四种类型中的哪一种，占一个字节，分别用00000001B表示颠换，00000010B表示插入，00000100B表示缺失，00001000B表示转换。
2.2.2 数据通信过程
(1)过程描述
在BooleanAlign与协处理器通信过程中，数据要分段多次发送接收，每次通信过程包括一次发送一次接收。要保证整个通信过程不出差错，需要一种机制来控制通信过程的进行。BooleanAlign采用MFC的线程机制，通过控制线程的阻塞与唤醒，实现数据通信过程中发送与接收的同步。
如图7(a)所示，首次发送数据时，send线程将first Wr置为1，然后进入WriteDev状态。在WriteDev状态，send线程截取DNA序列的一段并发送到协处理器。发送完毕且发送正确后，send线程等待50 ms(等待协处理器将本次序列比对完毕)，然后将信号eventRd置为有效，eventWr置为无效，send线程进入Block(阻塞)状态;否则，writeDevError有效，endRd被置为1,send线程进入writeDev Error状态，接着立即终止。在Block状态receive线程将接收并存储协处理器返回的数据，这将在下面叙述。当eventWr被receive线程置为有效时，send线程又进入WriteDev状态。此时，如果还有DNA序列要发送，send线程继续截取DNA序列的一段并发送到协处理器;如果DNA序列已经全部发送完毕，那么send线程进入终止态。
如图7(b)所示，当eventRd被send线程设置为有效后，receive进程被启动并进入ReadDev状态。在ReadDev状态，receive线程接收并存储协处理器返回的数据。这里有三种可能。第一，如果数据接收正确，则eventWr置为有效，eventRd置为无效，receive线程进入Block(阻塞)状态。第二，如果这是最后一段数据(lastWr=1)，则receive线程进入ReadDev End状态。第三，如果接收数据出错，则readDevError有效，endWr置为1,receive线程进入ReadDev Error状态，接着立即终止。在Block状态send线程发送下一段DNA序列。当eventRd被send线程置为有效时，receive线程又进入ReadDev状态。
(2)通信的实现细节
协处理器是一个USB2.0设备，FTDI公司提供了访问USB2.0设备的编程接口函数(如表2)，使应用程序可以直接调用接口函数实现与协处理器的通信。Windows操作系统本身带有USB设备的驱动程序，不过这是个通用的USB驱动程序，不能实现针对协处理器的最大传输率。FTD2XX.SYS可以针对协处理器的特点，对USB通信过程进行简化，从而实现通信速率的提高。如图3所示，在上层应用程序，WriteDev()函数调用FTD2XX.DLL里的FT_Write()库函数，实现向协处理器发送数据;ReadDev()函数调用FTD2XX.DLL里的FT_Read()库函数，实现接收数据。BooleanAlign使用隐式连接调用动态连接库文件中的函数。程序在运行时动态调用所需要的代码和数据，节省了内存占用。FTD2XX.SYS向上提供上层软件访问协处理器的接口，向下直接操纵协处理器，允许动态链接库函数通过驱动程序直接访问协处理器。FTD2XX.SYS采用批量传输方式，有一个接口，两个传输管道，每个管道支持最大64字节的数据包。
2.3 比对结果的处理
主机接收到的比对结果包含position和type，分别表示突变发生的位点和突变类型。应用程序根据这两种数据决定如何处理该位点的这对碱基。如果是转换和颠换，将该位点碱基以黑色显示;如果是插入，A序列上该位点的碱基与B序列上的“_”配对，B序列该位点以后的碱基依次后移一位;如果是缺失，A序列以“_”与B序列该位点的碱基配对，A序列该位点以后的碱基依次后移一位。假设，有序列S:TACAGTTGGATC，序列T::TTTGGA，其比对结果的集合R={<2，插入>，<3，插入>，<4，插入>，<5，插入>，<11，插入>，<12，插入>}。经应用程序处理以后，序列呈现为:
2.4 基于距离的系统发生树重建
非加权组平均法(UPGMA)[9]是一种简单的基于距离的系统发生树够建方法。Kimura模型[10]是常用的一种距离计算模型，其距离计算公式为:
式(2)中
d1表示比对排列后转换的位点数目与序列总长度的比值;d2表示比对排列后颠换的位点数目与序列总长度的比值，K表示两条DNA序列的距离[11]。构建基于距离的系统发生树，首先要建立一个距离矩阵，用来选择一对距离最小的DNA序列。程序中用一个二维数组(GDVector)存储距离矩阵，矩阵元素是一个结构体，名为Element，如图8所示。
Element存储了行序列与列序列的距离、行序列和列序列各自的名称以及本节点的名称。例如，图9(a)是a、b、c、d四条DNA序列的距离矩阵。a列b行处的矩阵元素(命名为E)，其结构体内容为:E.GD=0.000 161,E.up=“a”,E.down=“b”,E.name=“ab”。
聚合节点是BooleanAlign绘制系统发生树的基础，它代表了一个如图10所示的形状。每个聚合节点有四个顶点，图10示意了聚合节点u和v的8个顶点。BooleanAlign借助距离矩阵找到构建系统发生树所需的全部聚合节点，然后用MFC画图函数LineTo()绘制系统发生树。
以来自GeneBank加德纳菌属的四条序列a(GI号:297532249)、b(GI号:297532250)、c(GI号:297532251)、d(GI号:297532252)为例，构建系统发生树的步骤为:
Step 1:建立初始化距离矩阵(GDVector)。以a、b为例，经协处理器比对后得到的转换位点个数为12，颠换位点个数为26，序列总长度为195。根据式(2)，得K=0.098 008，将K填入距离矩阵的a列b行。依次计算序列ac、ad、bc、bd、cd之间的距离，得到初始距离矩阵如图9(a)所示。
Step 2:如果GDVector.size()=0，转到Step 6;否则，找出距离矩阵中GD最小的元素min Element。在本例中为元素ac，它的GD为0.049 519。
Step 3:计算聚合节点坐标。a、c首先聚合出一个节点u。根据BooleanAlign给定的绘图初始点P(x,y)和u的分支长度du=min Element.GD/2，计算u的四个顶点point1(x,y)、point2(x,y)、point3(x,y)、point4(x,y)坐标。
Step 4:删减距离矩阵。已经参加过聚合的序列不再参加下次聚合，要从距离矩阵中删掉。在本例中，删掉序列a、c。然后，计算u与其它序列之间的距离。u与序列m的距离dum=(dam+dcm)/2,m=b,d。删减后的距离矩阵如图9(b)所示。
Step 5:重复Step2。
Step 6:根据聚合节点，绘制系统发生树。
按照上述的方法，来自狼狗的五条DNA序列的系统发生树如图11所示。
分支a,b,c,d,e,f的长度如表3所示。
3 测试结果和分析
在Windows xp操作系统，1.61 GHz处理器主频，512 MB内存条件下，进行了两项测试。一项测试BooleanAlign在序列相似度较高情况下的性能，另一项测试BooleanAlign在序列相似度较低情况下的性能。每项测试包括五组测试数据，序列长度从300到16 195不等，每组测试都是双序全局列比对。测试数据如表4、表5所示，给出了在不同长度情况下，BooleanAlign和Needleman各自的性能。
注:max Len表示序列最大长度
令Tb表示BooleanAlign比对耗费时间，Tn表示Needleman算法比对耗费时间。在第一项高相似度测试中，五组数据的Tb/Tn的均值为0.06。在第二项低相似度测试中，五组数据的Tb/Tn的均值为0.42。
4 结论
测试结果显示，在两条序列相似度很高的情况下，BooleanAlign的比对耗费时间是Needleman算法的6%;在两条序列相似度较低情况下，BooleanAlign的比对耗费时间也仅是Needleman算法的42%。可以得出，BooleanAlign结合FPGA的速度优势和PC内存容量的优势，提高了双序列全局比对的速率。
参考文献
[1]孙啸,陆祖宏,谢建明.生物信息学基础.北京:清华大学出版社,2005:81
[2] Needlman S B,Wunsch C D.A genera method application to thesearch for similarities in the amino acid sequence of two proteins.Journal of Molecular Biology,1970,48(3):443—453
[3] Smith T F,Waterman M S.Identification of common molecular se-quence.Journall of Molecular Biology,1981;147:195—197
[4] Altschul S F,Gish W,Miller W,et al.Basic local alignment searchtool.Journall of Molecular Biology,1990;215(3):403—410
[5] Kent W J.BLAT—the BLAST-like alignment tool.Genome Res,2002;10:656—664
[6]张阳,窦勇,夏飞.生物信息学双序列比对算法加速器设计与实现.计算机科学与探索,2008;2(5):519—528
[7]郭静,王超,陈崚.基于遗传算法的系统发生树构建方法.计算机工程与应用,2009;45(16):72—76
[8][美]Shephend G,Wingo S.深入解析MFC.赵剑云,卿瑾,译.北京:中国电力出版社,2003:5—17
[9] Sokal R R,Michener C D.A statistical method for evaluating system-atic relationships.University of Kansas Scientific Bulletin,1958;28:1409—1438
[10] Kimura M.A simple method for estimating evolutionary rates of basesubstitutions through comparative studies of nucleotide sequences.Journal of Molecular Evolution,1980;16(2):111—120
编码理论在双序列比对中的应用 篇4
随着人类基因组计划的实施和完成以及分子生物学相关学科的迅猛发展, 越来越多的动物、微生物基因组序列得以测定。基因序列数据正在以前所未有的速度迅速增长。然而, 怎样去研究如此众多的基因在生命过程中所担负的功能就成了全世界生命科学工作者共同的课题。DNA (脱氧核糖核酸) 科学的发展就是顺应这一科学发展要求的产物。遗传物质的主要载体是染色体, 而控制生物性状遗传的主要物质是DNA。组成一个链的四个核苷酸是腺嘌呤 (A) 、腺嘧啶 (T) 、鸟嘌呤 (G) 和胞核嘧啶 (C) , 他们通常称为碱基。A和T配对, G和C配对。在分子生物学中, 一个DNA序列是由A, T, G, C组成的序列。DNA序列不产生一个马上有用的或提供消息的刻画。即使被考虑的DNA是短的片段, 比较他们也是困难的[1]。
生物体包含两类核酸, DNA和RNA, DNA是遗传信息的载体, 但DNA并不直接决定蛋白质的合成, 需要通过RNA才能把DNA上控制蛋白质合成的遗传信息传递给核糖体。同时RNA还参与蛋白质的合成, 与细胞分化、代谢、记忆的储存等有重要关系。人们发现RNA (而非DNA) 包含了一些病毒比如HIV中的基因信息, 而这些基因信息与这些病毒功能有关系。由于它的催化功能, RNA目前获得了高度的关注, 并吸引人们研究它的结构信息。为了更好地获取生物序列或结构信息, 近几年科学研究工作者发展了一种便于可视化的新的序列或结构表达方式—图形表示, 并基于图形表示给出了序列相似性分析、突变分析和进化分析方法[2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21]。
基于DNA序列四个碱基的物理化学性质分类, 给出了二联体一种四位二进制字符串表示方法, 并基于编码序列和异或操作给出了一种简单快速地点突变分析方法, 现将基于核苷酸二联体编码序列, 给出一种快速双序列比对算法。
1 基于序列编码的双序列比对
根据DNA序列四个碱基的物理化学性质分类, 给出了二联体一种四位二进制字符串表示方法:φ=a1a2a3a4, ai∈{0, 1}, i=1, 2, 3, 4, φ∈Ω, 其中Ω={AA, CC, GG, TT, AC, AG, AT, CA, CG, CT, GA, GC, GT, TA, TC, TG}, 这样将每条DNA序列转换成0, 1序列, 并利用异或操作进行序列比较, 很好地解决了序列比较的四个基本问题: (1) 对于两条长度相近的序列相似, 找出序列的差别; (2) 判断一条序列的前缀与另一条序列的后缀相似; (3) 判断一条序列是否是另一条序列的子序列; (4) 判断两条序列中是否有非常相似的子序列。
生物信息学计算的核心是序列比较, 而序列比较往往通过序列比对来实现, 序列比对目的就是找到尽可能多的匹配字符。现将利用异或操作, 通过搜寻最长公共连续0子序列来寻找最优比对。
定理1[22]:α⊕β=0000⇔α=β。
推论1[22]:如果α⊕β≠0000, 则α与β不匹配。
结合定理1和推论1可知, 如果两序列有公共字符, 则它们异或所得结果为0000, 通过搜寻连续0000序列给出上序列比对算法如下:设比对两序列分别为a, b则有:
例, 寻找序列a= TAGGCCTCTGCCTAATCACA-CAG, b= CGGCCTCTGCCTTATTACACAA, 首先定义16个二联体编码规则如下。
AA:0000, AC:0001, AG:0010, AT: 0011, GA: 1000, GC: 1001, GG: 1010, GT: 1011, CA: 0100, CC: 0101, CG: 0110, CT: 0111, TA: 1100, TC: 1101, TG: 1110, TT: 1111
根据定义的编码规则将序列a和b转换成0, 1序列
利用比对算法, 找到第一个最长公共0序列l=0000000000000000000000000000000000000000, 其对应的字符序列为S=GGCCTCTGCCT, 这样将序列a和b分别表示成a=A1SA3, b=B1SB3, 利用同样的方法寻找A1和B1, A3和B3的比对, 序列a和b的最终比对结果如下:
TA-GGCCTCTGCCT-AATC-ACACA-G-CGGCCTCTGCCTT-AT-TACACAA-
2 结束语
基于核苷酸二联体的一种编码规则, 利用异或操作的结果, 通过搜寻最长公共0序列来寻找最优比对。虽然这种二值化的思想使得序列变得更长了, 但是利用异或操作, 每次考虑的是四位二进制数的异或结果, 并未影响序列比对。这种二值化的思想既提高了电子计算机的可操作性, 也提高了计算的实时性, 计算简单, 也能体现各核苷酸的分布规律。
序列比对 篇5
核基因组序列中的内转录间隔区 (Internal transcribed spacer, ITS) 主要编码植物的核糖体RNA, 属于中度保守序列, 因此可用于较低分类阶元的系统发育研究[6,7]。本文以Ziziphus Mill.的部分ITS序列为对象, 介绍了Clustal X2.1[8,9]、Boxshade、DNAman和MEGA6[10]等软件在序列比对、着色美化、序列分析和植物系统发育树构建中的应用。
1 序列获取
从NCBI网站 (http://www.ncbi.nlm.nih.gov/) 下载Ziziphus Mill.的4条ITS核酸序列, 以Fasta格式存储于桌面文档中, 命名为four sequences.txt。
2 使用Clustal X进行多序列比对
2.1 打开Clustal X软件
找开Clustal X软件 (Version 2.1, 2010年发布, http://www.clustal.org/) , 在“Multiple alignment”模式下点击“File”菜单→“Load sequece”, 打开存储于桌面上的four sequences.txt文件 (图1) 。
2.2 进行序列比对
打开“Alignment”菜单→“Do complete alignment”, 系统提示“Do complete alignment”结果的输出文件类型及保存位置 (图2) 。其中, dnd格式为输出向导树文件, 可以用Treeview软件打开, aln为输出比对文件格式, 打开后对序列进行手动调整。
2.3 aln文件的着色美化
利用Boxshade在线工具 (http://www.ch.embnet.org/software/BOX_form.htm L) 或Espript在线工具 (http://espript.ibcp.fr/ESPript/cgi-bin/ESPript.cgi) 对aln文件进行着色美化。将aln文件序列粘贴到Boxshade序列框中, 点击“Run boxshade”, 结果以.ps格式存储, 用Ghostview和Ghostspript软件打开, 如图3所示。
3 DNAman多序列比对和序列一致性
DNAman软件除了可对序列进行比对之外, 还可进行引物的设计、限制性内切酶和质粒图谱的绘制, 功能较为强大。双击DNAman图标, 打开“File”菜单→“Open”, 选择four sequences.aln文件, 再打开“Sequence”菜单→“Alignment”→“Multiple sequence alignment”, 从显示结果可看出, 这4条ITS序列的序列一致性为90.10% (图4) , 结果可以以Phylip、Clustal和GCG等多种格式输出。
4 MEGA进行系统发育分析
4.1 打开分子进化遗传分析 (Molecular evolution genetic analysis, MEGA) 软件
由于MEGA只能识别.meg文件格式, 故需要对.aln格式的文件进行转化。点击“File”→“Convert file format to MEGA”, 将.aln格式文件转化为.meg格式, 存储于桌面上, 命名为four sequences.meg。或进行在线文件格式转换 (http://sing.ei.uvigo.es/ALTER/) 。
4.2 主窗口中打开“Phylogeny”菜单
可以看到, MEGA提供了UPGMA法、NJ法、MP法和ML法等多种算法, 如前所述, NJ法在处理相似度较高、亲缘关系较近的序列时是最可靠的一种算法。选择“Construct/Test neighbor-joined tree”, 在弹出的对话框“Options summary”→“Test phylogeny”中, 选择“Bootstrap method”, 重复次数为1000次, 模式框中选择Kimura-2-parameter model对进化树进行评估 (图5) , 本次所构建的系统进化树如图6所示。
从构建的NJ树可知, 基因登录号为EU075099和EU075097的序列聚为一类, bootstrap支持率为97%, 表明该进化枝的可信度较高。
5 结语
本文以Ziziphus Mill.的部分ITS序列为对象, 介绍了Clustal X、Boxshade、DNAman和MEGA等软件在序列比对、着色美化、序列分析和植物系统发育树的构建中的应用。在此需要特别指出, 在基于序列构建进化树的过程中, 由于不同的算法, 不同的重复次数, 选用不同的模型所构建的树是不完全相同的[11], 因此, 这种单纯的依靠某段核酸序列建立的进化树只能提供物种进化的部分信息, 而不能完全代表该物种进化的全过程。
摘要：通过Clustal X2.1、Boxshade、DNAman和MEGA6等软件的使用, 介绍了生物学软件在核酸序列比对、着色美化、序列分析和植物系统发育树构建中的应用。并结合具体实例, 采用邻接法对核基因组中的内转录间隔区序列进行了构树, 进化枝的可信度较高, 表明该方法适用于相似度较高、亲缘关系较近序列的系统发育树的构建。
关键词：序列比对,邻接法,系统发育树,构建
参考文献
[1]王禄山, 高培基.生物信息学应用技术[M].北京:化学工业出版社, 2007.
[2]谢强, 卜文俊.核苷酸序列比对在生物系统发育研究中的应用[J].动物分类学报, 2005 (2) :281-286.
[3]杨学森.基于汉明距离的DNA短序列比对算法研究[D].哈尔滨:哈尔滨工业大学, 2013.
[4]路明.利用进化树研究基因序列的进化[D].南宁:广西大学, 2014.
[5]高凯.NJ进化树构建方法的改进及其应用[D].北京:北京工业大学, 2008.
[6]樊杰, 白妍, 束明月.远志属7种药用植物ITS1和ITS2序列分析[J].中草药, 2015 (4) :562-565.
[7]BRUCE G, BALDWIN.Phylogenetic utility of the internal transcribed spacers of nuclear ribosomal DNA in plants:An example from the compositae[J].Molecular Phylogenetics and Evolution, 1992, 1 (1) :3-16.
[8]LARKIN M A, BLACKSHIELDS G, BROWN N P, et al.Clustal W and Clustal X version 2.0[J].Bioinformatics, 2007, 23:2947-2948.
[9]李彬彬, 黄培春, 钟复光.生物学软件在线粒体DNA序列多态性分析中的应用[J].生物信息学, 2010 (2) :153-155.
[10]TAMURA K, STECHER G, PETERSON D, et al.MEGA6:Molecular Evolutionary Genetics Analysis Version 6.0[J].Molecular Biology and Evolution, 2013, 30:2725-2729.