加载中…
个人资料
じ☆ve涛的天空
じ☆ve涛的天空
  • 博客等级:
  • 博客积分:0
  • 博客访问:1,742
  • 关注人气:0
  • 获赠金笔:0支
  • 赠出金笔:0支
  • 荣誉徽章:
相关博文
推荐博文
谁看过这篇博文
加载中…
正文 字体大小:

序列比对和数据库搜索

(2011-03-09 19:16:49)
标签:

杂谈

引言

在生物学的研究中,有一个常用的方法,就是通过比较分析获取有用的信息和知识。达尔文正是研究比较了galapagos finches同其它一些物种的形态学特征,从而提出了自然选择学说。今天,我们对基因和蛋白质序列进行比较,从本质上来讲是同达尔文一样,进行同样的分析,只不过更加精细,更加详尽。在这个意义上,我们从核酸以及氨基酸的层次去分析序列的相同点和不同点,以期能够推测它们的结构、功能以及进化上的联系。最常用的比较方法是序列比对,它为两个或更多个序列的残基之间的相互关系提供了一个非常明确的图谱。在这一章,我们只讨论一下双重比对,即只比较两个序列,至于较多的序列即多序列比对,将在 下一章介绍。

七十年代以来,DNA测序方法的飞速发展,极大地引发了序列信息量的扩增,从而使可供比较的序列数量呈现爆炸式增长。分子生物学家应该意识到,将未知序列同整个数据库中的已知序列进行比较分析已经成为他们手中一个强有力的研究手段。在过去的三十年里,即使不提及计算机的应用,序列比较的各种算法也已经发展得越来越迅速,也越来越成熟,已经能够跟上序列数据库增长的步伐。今天,我们已经拥有一些小的模式物种的基因组的全序列,还拥有人类基因序列的一些较大的样品,我们已经进入比较基因组时代,也就是说,对两个物种进行全基因组序列比较已经不再是一个梦想。

序列比对的进化基础

进行序列比对的目的之一是让人们能够判断两个序列之间是否具有足够的相似性,从而判定二者之间是否具有同源性。值得注意的是,相似性和同源性虽然在某种程度上具有一致性,但它们是完全不同的两个概念。相似性是指一种很直接的数量关系,比如部分相同或相似的百分比或其它一些合适的度量,而同源性是指从一些数据中推断出的两个基因在进化上曾具有共同祖先的结论,它是质的判断。基因之间要么同源,要么不同源,绝不象相似性那样具有多或少的数量关系。如图7.1所示,比较家鼠和小龙虾的同源的胰蛋白酶序列,发现它们具有41%的相似性。

由于受到研究进化关系这一目的的影响,大多数比对方法很自然地都希望能够在某种程度上建立起分子进化的模型。我们通常都假定同源序列是从某一共同祖先不断变化而来,但事实上,我们无法得知这个祖先序列到底是什么样子,除非能够从化石中获得它的DNA,我们所能够做到的只是从现存物种中,探求真相。从祖先序列以来所发生的变化包括取代插入以及缺失。在理想情况下,同源基因或蛋白质序列在相互比较时,残基之间相互对应,从而使取代的情况很明显地表现出来。在某些位置,一个序列中拥有某些残基而另一个序列中缺少这种残基,表明这些残基是插入到前者或是从后者中丢失的。这些空位在序列比对时用连续的短线填补。如图7.1,在序列比对中,发现了5个空位。

                 |------ S-S-----|

Mouse    IVGGYNCEENSVPYQVSLNS-----GYHFCGGSLINEQWVVSAGHCYK-------SRIQV

Crayfish IVGGTDAVLGEFPYQLSFQETFLGFSFHFCGASIYNENYAITAGHCVYGDDYENPSGLQI

                                              *

Mouse    RLGEHNIEVLEGNEQFINAAKIIRHPQYDRKTLNNDIMLIKLSSRAVINARVSTISLPTA

Crayfish VAGELDMSVNEGSEQTITVSKIILHENFDYDLLDNDISLLKLSGSLTFNNNVAPIALPAQ

                          |----- S-S--------|

Mouse    PPATGTKCLISGWGNTASSGADYPDELQCLDAPVLSQAKCEASYPG-KITSNMFCVGFLE

Crayfish GHTATGNVIVTGWG-TTSEGGNTPDVLQKVTVPLVSDAECRDDYGADEIFDSMICAGVPE

           +*|----------S-S----------|

Mouse    GGKDSCQGDSGGPVVCNG----QLQGVVSWGDGCAQKNKPGVYTKVYNYVKWIKNTIAAN

Crayfish GGKDSCQGDSGGPLAASDTGSTYLAGIVSWGYGCARPGYPGVYTEVSYHVDWIKANAV--

图7.1、保守位点通常在功能上极为重要。对老鼠的胰蛋白酶(Swiss-Prot P07146)和小龙虾的胰蛋白酶(Swiss-Prot P00765)作比对,相同的残基用下标线标出,在比对上方标出的是三个二硫键(-S-S),这些二硫键中的半胱氨酸残基极为保守,打*号的残基的侧链参与电荷传递系统,打+号的活性位点的残基负责底物的特异性。

在残基-残基比对中,很明显,某些位置的氨基酸残基相对于其它位置的残基具有较高的保守性,这个信息揭示了某些残基对于一个蛋白质的结构和功能是极为重要的。如图7.1所示,处于活性位点的残基都是极为保守的,比如形成二硫键的半胱氨酸,参与电子传递的氨基酸残基以及决定底物特异性的氨基酸残基。这些保守的残基对于保持蛋白的结构与功能非常重要,另一方面,由于历史原因,某些保守位置对蛋白功能并无太大的重要性。当我们处理非常相近的物种时必须十分小心,因为相似性在某些情况下更多地是历史的反映而不是功能的反映,比如,mouse和rat的某些序列具有高度的相似性,可能仅仅是因为没有足够的时间进行分化而已。尽管如此,系列比对仍然是从已知获得未知的一个十分有用的方法,比如通过比较一个新的蛋白同其它已经经过深入研究的蛋白,可以推断这个未知蛋白的结构与功能的某些性质。必须指出的是,不能够仅仅是通过比较分析这一判据来断定结论是否正确,结论还必须经过实验验证。

当我们发现两个基因或蛋白质具有惊人的相似性时,我们会认为他们之间具有一段共同的进化历程,从而我们判断他们会具有相似的生物学功能,但是,这个推断在成为结论之前必须经过实验的验证。例如,ζ-晶状物是脊椎动物眼睛里晶状体基质的组成部分,根据序列相似性的基础,它在E.coli中的同源物是代谢酶苯醌氧化还原酶(如图7.2),不管二者的共同祖先如何,它们的功能在进化中已经改变了(Gonzalez et al.,1994)。这就好象火车变成了铁路餐车,虽然对二者的外部结构的观察揭示了它们结构的历史,但是仅仅根据这一信息往往会得出有关其功能的错误结论。当一个基因适应了一个新的功能时,保守位置通常也会发生一些形式上的变化,比如,当蛋白具有催化功能时,活性为点的残基相当保守,而当蛋白功能改变时,这些残基将会发生漂移。

Human-ZCr MATGQKLMRAVRVFEFGGPEVLKLRSDIAVPIPKDHQVLIKVHACGVNPVETYIRSGTYS

Ecoli-QOR ------MATRIEFHKHGGPEVLQA-VEFTPADPAENEIQVENKAIGINFIDTYIRSGLYP

. . ******. . . * …. . . * *.* ..****** *

Human-ZCr RKPLLPYTPGSDVAGVIEAVGDNASAFKKGDRVFTSSTISGGYAEYALAADHTVYKLPEK

Ecoli-QOR -PPSLPSGLGTEAAGIVSKVGSGVKHIKAGDRVVYAQSALGAYSSVHNIIADKAAILPAA

* ** *.. **.. ** . * **** . . * *. **

Human-ZCr LDFKQGAAIGIPYFTAYRALIHSACVKAGESVLVHGASGGVGLAACQIARAYGLKILGTA

Ecoli-QOR ISFEQAAASFLKGLTVYYLLRKTYEIKPDEQFLFHAAAGGVGLIACQWAKALGAKLIGTV

. * * ** . * * * .. .* * * * *.***** *** *.* * *..**

Human-ZCr GTEEGQKIVLQNGAHEVFNHREVNYIDKIKKYVGEKGIDIIIEMLANVNLSKDLSLLSHG

Ecoli-QOR GTAQKAQSALKAGAWQVINYREEDLVERLKEITGGKKVRVVYDSVGRDTWERSLDCLQRR

** . . *. ** .* * ** …. * * * . .. . . . . * * .

Human-ZCr GRVIVVG-SRGTIEINPROTMAKES----SIIGVTLFSSTKEEFQQYAAALQAGMEIGWL

Ecoli-QOR GLMVSFGNSSGAVTGVNLGILNQKGSLYVTRPSLQGYITTREELTEASNELFSLIASGVI

* .. * * *.. . . . . . .*.** . . * . . * .

Human-ZCr KPVIGSQ--YPLEKVAEAHENIIHGSGATGKMILLL

Ecoli-QOR KVDVAEQQKYPLKDAQRAHE-ILESRATQGSSLLIP

* . * *** *** *. . * .*.

图7.2、最佳全局比对:对人类ζ-晶状物(Swiss-Prot Q08257)和E.coli苯醌氧化还原酶(Swiss-Prot P28304)的氨基酸序列进行比对。这是一个由CLUSTAL W程序(Higgins et al., 1996)得到的最佳全局比对结果。在比对下方,星号表示残基相同,打点表示这个残基是保守的。

早期的序列比对方法只应用于那些在全长范围内具有简单相似性的一些序列。全序列比对就是对序列进行全程扫描,进行比较。以上讨论的胰蛋白酶和ζ-晶状物之间的比较就属于全序列比对。具有简单的球形结构域的蛋白一般可以使用全序列比对的策略,以为所有的同源序列尚未经过实质上的变化

蛋白质的模块性质

许多蛋白质在全程范围内并不具有相似性,但却似乎是由众多的模块结构域搭建而成。图7.3描述了这样的一个例子,如图所示的是在血凝过程中的两种蛋白的组成结构,它们是凝血因子XII(F12)和组织型血纤蛋白溶酶原活化因子(PLAT),除了具有丝氨酸蛋白酶活性的催化结构域,这两种蛋白还具有不同数量的其它结构域单元,包括两种纤连蛋白重复,一个类似于上皮生长因子的结构域以及一个成为“kringle”域的单元。这些组分可以以不同顺序反复出现,组分形式的不同通常是由于整个外显子交换引起的。由于全程比对建立时,基因的外显子/内含子结构还没有被发现,因此全程比对并没有顾及到上述现象的重要性,这是可以理解的。在大多数情况下,使用局部比对是较为合理的,这种比对方法可能会揭示一些匹配的序列段,而本来这些序列段是被一些完全不相关联的残基所淹没的,因此,操作者应该明白,如果不恰当地使用了全程比对,很可能会掩埋一些局部的相似性。设计局部比对的另外一个很明显的原因就是在比较一个拼接后的mRNA和它的基因序列时,每个外显子都应该进行局部比对。

图7.3、血凝过程中的两中蛋白的模块结构:人类组织血纤蛋白溶酶原活化因子以及凝血因子XII的模块结构的示意图。标记为Catalytic的模块在若干种凝血蛋白中是常见的,F1和F2是较为常见的重复模块,首先在纤连蛋白中被发现。E模块同表皮生长因子极为类似。通常称为”Kringle domain”的模块被标记为K。

点阵描述方法之所以广泛流行,其部分原因就在于它能够揭示出拥有多个局部相似性的复杂关系,图7.4就是应用这种处理后的一个例子。图中F12和PLAT蛋白质序列使用DOTTER程序进行比较(软件可见本章结尾列表),其基本思路就是把两个序列分别作为一个二维坐标系中的两个坐标轴,在这个坐标系区域内,如果某一点所对应的横轴坐标和纵轴坐标所对应的两条序列的残基相同,则在这个位置上打上标记点,每个点通常都表示在一些小窗口中,序列相似性高于其它一些隔绝的区域(或者由DOTTER程序定义的隔绝区域,由不同的灰色阴影标记)。如果两个序列在一段区域内很相似,标记点将会连成一条斜线段,将这些线段的位置同图7.3中两个蛋白的已知的组成结构相比较是很有价值的,特别是要注意连续反复出现的结构域的出现方式。从PLAT的kringle结构域开始水平扫描,可以发现两条线段对应于F12序列中的两个kringle结构域,虽然现在我们已经拥有许多更复杂更精确的方法来寻求局部相似性(下面将会讨论),点阵描述方法仍然是一个很流行很有效的描述方法。

图7.4、点阵序列比较:对人类凝血因子XII(F12:Swiss-Prot P00748)和组织血纤蛋白溶酶原活化因子(PLAT:Swiss-Prot P00750)的氨基酸序列进行打点比较。这个图由DOTTER程序(Sonnhammer and durban,1996)产生。

在点阵描述方法中,某些形式的点可能会勾勒出一定的路径,但这需要操作者通过这些信息进行推理,另外一个图形描述方法即路径图提供了更直接明了的比较结果,图7.5描述了PLAT和PLAU中与EGF相似的结构域之间进行比较时的比对、点阵和路径图三种方法的关系。

c

PLAU 90 EPKKVKDHCSKHSPCQKGGTCVNMP—SGPH-CLCPQHLTGNHCQKEK---CFE 137

PLAT 23 ELHQVPSNCD----CLNGGTCVSNKYFSNIHWCNCPKKFGGQHCEIDKSKTCYE 72

图7.5、点阵、路径图和比对:所有这三种视图都表示人类尿激酶血纤蛋白溶酶原活化因子(PLAU:Swiss-Prot P00749)和组织血纤蛋白溶酶原活化因子(PLAT:Swiss-Prot P00750)中同EGF相似的模块的比对结果。a) .整个蛋白都由DOTTER程序进行比较:这里只显示了同EGF模块相似的较小区域的放大图;b)由BLASTP得到的比对的路径图;.c).用普通的字符形式显示的BLASTP空位比对。

要理解路径图,先想象一个二维格子,顶点表示序列残基之间的点(与点阵中表示残基本身相反),沿线段上连接两个顶点的边缘对应两个序列上匹配的残基,水平和竖直线段的边缘对应一个序列拥有而另一个序列上没有的残基,换句话说,这些边缘平台组成了比对中的空位,全图对应了所有可能的比对中必须审视的搜索空间,这个空间中每条可能的路径都对应于一种比对。

最佳比对方法

除了某些很不重要的问题,对于众多问题而言,比对方法多种多样,很有必要从中挑选出最好的一个或几个方法,这就是把一种比对描述成一个路径的概念所指。许多计算机科学的问题都可以简化为通过图表寻求最优路径(比如寻找从纽约打电话到旧金山的最有效的途径)。为了这一目的已经确立了许多行之有效的算法,对每一种路径都有必要对其进行某种意义上的打分,通常是对沿这一途径的每一步的增量进行加和。更精密的打分程序将在下文叙述,在这里我们只假定相同残基加正分,有插入或缺失的残基就加负分(扣分),根据这一定义,最合适的比对方法会得到最高分,也就是我们寻找的最佳路径。

今天我们所熟悉的Needleman-Wunsch算法就是针对寻求最佳序列比对这一问题所设计的动态规划寻优策略(Needleman and Wunsch,1970)。动态规划的思想是这样的,如果一条路径终止于最佳路径上的一点,那么这条路径本身就是起点到这个中间点的最佳路径,也就是说,任何一个终止于最佳路径上的一点的次级路径必然就是终止于这一点的最佳路径本身。这样,最佳路径就可以通过把各个最佳的次级路径连接而成。在基本的Needleman-Wunsch公式表达中,最佳比对必然对每个序列都由始至终,就是说从搜索空间的左上角直至右下角。换句话说,它搜索全程比对。

然而,对这种基本策略稍作修改就可以实现最佳的局部比对。这种比对的路径不需要到达搜索图的尽头,只需要在内部开始和终结。如果某种比对的打分值不会因为增加或减少比对队的数量而增加时,这种比对就是最佳的。这个过程依赖于打分系统的性质,就是说某种路径的打分会在不匹配的序列段位置减少(以下叙述的打分系统合乎这个标准)。当分值降为零时,路径的延展将会终止,一个新的路径就会应运而生。这样,我们会得到许多独立的路径,它们以不匹配的序列段为界限而不是像在全程比对中以序列的结尾作为界限。在这些路径中,拥有最高分的一个就是最佳的局部比对。

应该意识到,寻优方法总是把最佳的比对方法表达出来,而不在意它是否具有生物学意义,另一方面,寻求局部比对时可能会发现若干个重要的比对,因此,不能仅仅注意最佳的一个。改良的Smith-Waterman(Altschul and Erickson,1986;Waterman and Eggert,1987)算法把寻找K种最好的但不相互交叉的比对方式最为目标,这些思想后来都在SIM算法(Huang et al.,1990)的发展中得以体现。一个名叫LALIGN(在FASTA程序包中)的程序提供了有用的SIM工具(Pearson,1996)。对于比对多模块的蛋白质而言,寻找次优比对尤为重要。正如图7.6所示,LALIGN程序被用来获得三个最好的局部比对(比对人类凝血因子IX和因子XII)。一个标准的Smith-waterman算法只会报告出最好的一个比对,改良的算法会报告出第二和第三的比对方式,从而显示出功能结构域。

Comparison of:

  1. f9-human.aa >f9 gi|119772|sp|P00740|FA9_HUMAN COAGULATION FA -461 aa
  2. f12-hum.aa>f12 gi|119763|sp|P00748|FA12_HUMAN COAGULATION -615 aa

using protein matrix

① 35.4% identity in 254 aa overlap; score: 358

220 230 240 250 260 270

F9 QSFNDFTRVVGGEDAKPGQFPWQVVLNGKVDAFCGGSIVNEKWIVTAAHCVE---TGVKI

.:....:::::: : .:. :. ..: ..::.::... :..:::::.. . ..

F12 KSLSSMTRVVGGLVALRGAHPYIAALY-WGHSFCAGSLIAPCWVLTAAHCLQDRPAPEDL

370 380 390 400 410 420

280 290 300 310 320 330

F9 TVVAGEHNIEETEHTEQKRNVIRIIPHHNYNAAINKYNHDIALLELDEPL-----VLNSY

::: :... ... .. :. .: . :...... .:.::.::: :.: .:..:

F12 TVVLGQERRNHSCEPCQTLAVRSYRLHEAFSPV--SYQHDLALLRLQEDADGSCALLSPY

430 440 450 460 470 480

340 350 360 370 380

F9 VTPICIADKEYTNIFLKFGSGYVSGWGRVFHKGRS-ALVLQYLRVPLVDRATCLRSTKF-

: :.:... . .. :.:::. :. . . : :: .::... . : ..

F12 VQPVCLPSGAARPSETTLCQ—VAGWGHQFEGAEEYASFLQEAQVPFLSLERCSAPDVHG

490 500 510 520 530

390 400 410 420 430 440

F9 -TIYNNMFCAGFHEGGRDSCQGDSGGPHVTEVEGTS---FLTGIISWGEECAMKGKYGIY

.: .:.:::: ::: :.:::::::: : : .... : ::::::..:. ..: :.:

F12 SSILPGMLCAGFLEGGTDACQGDSGGPLVCEDQAAERRLTLQGIISWGSGCGDRNKPGVY

540 550 560 570 580 590

450

F9 TVVSRYVNWIKEKT

:.:. :..::.:.:

F12 TDVAYYLAWIREHT

600 610

------------------------------------

② 34.7% identity in 49 aa overlap; score: 120

100 110 120 130 140

F9 VDGDQCESNPCLNGGSCKDDINSYECWCPFGFEGKNCELDVTCNIKNGR

.....: .::::.::.: . . : :: :..: :..:.. . .::

F12 LASQACRTNPCLHGGRCLEVEGHRLCHCPVGYTGPFCDVDTKASCYDGR

180 190 200 210 220

-------------------------------------

③ 33.3% identity in 36 aa overlap; score: 87

100 110 120

F9 DQCESN-PCLNGGSCKDDINSYECWCPFGFECKNCE

:.:... :: .::.: . .. .: :: ..:..:.

F12 DHCSKHSPCQKGGTCVNMPSGPHCLCPQHLTGNHCQ

100 110 120 130

--------------------------------------

图7.6、最佳和次佳的局部比对:在使用LALIGN对人类凝血因子IX(F9;Swiss-Prot 900740)和凝血因子XII(F12;Swiss-Prot P00748)进行比对时发现了三个最佳的比对结果。

取代分和空位处罚

刚才描述的打分系统仅仅使用于简单的匹配/不匹配的情况,但是在比较蛋白质时,我们可以用取代矩阵来增强弱势比对的敏感性。很显然,在相关蛋白质之间,某些氨基酸可以很容易地相互取代而不用改变它们的生理生化性质,这些保守取代的例子包括异亮氨酸(isoleucine)和颉氨酸(valin)(体积小,疏水),丝氨酸(serine)和苏氨酸(threonin)(极性)。在计算比对分之时,相同的氨基酸打分会高于取代的氨基酸,而保守的取代打分高于非保守变化,换句话说,设计了一系列的分值,而且,在比对非常相近的序列(mouse和rat的同源基因)以及差异极大的序列(mouse和 yeast的基因)时会设计出不同系统的分值,考虑到这些因素,使用取代矩阵会极为有利,在这个矩阵中,任何氨基酸配对的分值会一目了然。

第一个广泛使用的最优矩阵建立在进化的点突变模型上(PAM)(Dayhoff et al.,1978)。一个PAM就是一个进化的变异单位即1%的氨基酸改变,这并不意味着经过100次PAM后,每个氨基酸都发生变化,因为其中一些位置可能会经过多次改变,甚至可能变回到原先的氨基酸,因此另外一些氨基酸可能不发生改变。如果这些变化是随机的,那么每一种可能的取代频率仅仅取决于不同氨基酸的出现的频率(称为背景频率)。然而,在相关蛋白中,已经发现的取代频率(称为目标频率)大大地倾向于那些不影响蛋白质功能的取代,换句话说,这些点突变已经被进化所接受。Dayhoff同合作者们第一次使用了log-odd处理,在这种处理中,矩阵中的取代分值同目标频率于背景频率的比值的自然对数成比例。为了评估目标频率,人们用非常相近的序列(比对时不需要取代矩阵)来收集对应于一个PAM的突变频率,然后将数据外推至250个PAM,PAM250矩阵结果如图7.7。虽然Dayhoff等人只发表了PAM250,但潜在的突变数据可以外推至其它PAM值,产生一组矩阵,在比较差异极大的序列时,通常在较高的PAM值处得到最佳结果,比如在PAM200到250之间,较低值的PAM矩阵一般使用于高度相似的序列(Altschul,1991)。

图7.7、PAM250分值矩阵。

用同样方式建立了BLOSUM取代矩阵,但在评估目标频率时,应用了不同的策略,基本数据来源于BLOCKS数据库,其中包括了局部多重比对(包含较远的相关序列,同在PAM中使用较近的相关序列相反)。虽然在这种情况下,没有进化模型,但它的优点在于可以通过直接观察获得数据而不是通过外推获得。同PAM模型一样,也有许多编号的BLOSUM矩阵,这里的编号指的是序列可能相同的最高水平,并且同模型保持独立性。举例来说,如图7.8所示的BLOSUM的矩阵,至少有62%的相同比例的序列被组合成一个序列,因此取代频率更加受到那些比空位变化还大的序列的极大影响,取代矩阵在处理高度相似序列时使用高的阈值(直至BLOSUM90),处理差异大的序列时使用低的阈值(直至BLOSUM30)。

图7.8、BLOSUM62分值矩阵。

为了补偿那些插入或缺失,可以在比对中引入一些空位,但不能太多,否则会使分子变得面目全非。每引入一个断裂,比对的分值都会有所扣除,对于这些断裂有许多罚分的规则。最常用的一个就是用一个附加的罚分比例去乘空位的长度,其中有两个参数:G(有时称为断裂开放惩罚)和L(断裂延伸惩罚),对于一个长度为n的空位,扣分总数为G+Ln,但在选择空位参数时,在很大程度上是唯经验的,所选的分值很少会有理论上的支持。通常来说,对于G会选择一个高分(在BLOSUM62中约为10-15),对于L会选择一个相对的低分(大约1-2),选择这个范围是因为插入和变异是很罕见的,但当它们一旦发生,就会影响到一系列附近的残基

比对的统计学显著性

对任何一个比队,我们都可以计算一个分值,但重要的是需要判定这个分值是否足够高,是否能够提供进化同源性的证据。在解决这一问题时,对于偶然出现的最高分,有些思想很有帮助,但是,没有一个数学理论能够描述全程比对的分值分布,其中一个能评估其重要性的方法就是将所得的比对分值和那些同样长度和组成的随机序列进行比较。

但是,对于局部比对而言,情况要好得多。正如问题总是从简单开始,人们首先注意到那些没有多少空位得局部比对,这种比对被称为高分片段配对(HSP)。HSP通常用改进得Smith-waterman算法或简单地使用大的空位罚分方法获得。Karlin-Altschul统计学为描述随机的HSP分值的分布提供了数学理论,概率密度函数形式被称为极值分布,这很值得注意,因为,更普遍更一般的分布的应用可能会夸大它的重要性,把一个已知得比对分值S同预期的分布相关联可能会计算出P值,从而给出这个分值的比对显著性的可能性。通常,P值越趋近于零,分值越有意义。

相关的变量E表示分值不低于S得可能的比对数量,而极值分布由两个参数表示,即Kλ,可以得到解析解,并且对于任何打分系统以及背景频率都是固定的。比对的显著性依赖于搜索空间的大小(就像在草堆中找针依赖于草堆的大小)。搜索空间的大小由序列长度计算出来,但由于统计的正确性,这个长度必须由局部比对的预期长度进行校正,以免出现边缘效应(Altschul and Gish,1996),需要进行这种校正还因为在搜索空间边缘开始的比对在达到一个有效分值之前就会超出序列的范围。

把比对局限于没有空位的基础之上,使问题大大简化,但是却脱离分子生物学的实际情况。实际上,要建立一个插入和缺失的精确模型需要空位,但如果空位相对较少,在这些空位之间仍然可以获得高分值区域,有代表性的是可能会获得紧密相邻的HSP,在这种情况下,从总体上去评估它的显著性是较为合理的,也许,每个片段并不显得很重要,但是几个片段同时出现就不太像是偶然事件了。Karlin-Altschul加和统计学可以计算N个HSP的统计值,这个方法的实质是把N个最佳片段的分值进行加总,从而计算事件偶然发生的可能性,其它一些论据也被用来确认这些分值只是在片段与比对一致的情况下进行加总。虽然加总的分值分布与HSP分值最大值有差异,仍然可以得到解析解。

最后,仍然有必要对局部排队的显著性进行合理评估,其中包括了模型中的空位。正如同传统的Smith-waterman比对,虽然没有先验的证据,人们仍然认为这些比对的分值也应该遵循极值分布,但是,分布参数K和λ的值不能通过计算获得,当然,通过模型获得这些值的方法已经被大大地发展了。

数据库中的相似性搜索

上述讨论主要集中于那些较为特别的匹配的序列,但是对于一个新发现的序列,我们无法得知用什么序列同它进行比对,数据库相似性搜索使我们能够从数据库中存在的数十万个序列中挑选出可能同感兴趣的序列有关联的序列,这个方法有时会导致意想不到的收获。用这种策略获得成功的第一个例子是人们因此发现病毒肿瘤基因v-sis是细胞中编码血小板派生生长因子的基因的一个变体形式(Doolittle et al., 1983; Waterfield et al., 1983)。那个时候,序列数据库还不大,因此这个发现足以另人感到万分惊奇。然而今天如果进行数据库搜索并且一无所获的话,那就更另人感到费解了。如同其它几个小的物种基因组一样,酵母saccharomyces cerevisiae的基因组全序列已经被测定出来。在脊椎动物中,大量的部分基因诸如人类和老鼠的基因都已经被测定并存入基因库(genebank)中,这也导致了表达序列标签(EST)工程。EST片段的主要用途是在数据库搜索中,用EST片段进行cDNA克隆可以分离出感兴趣的基因,包括其它模型生物中的同源基因。最近报导的多重内分泌腺肿瘤(MENI)基因就和人与老鼠的多个EST片段相匹配,其中之一在MENI发表前一年就已经入库保存了(Chandrasekharappa et al., 1997)。

在数据库搜索中,基本操作就是将查询序列和数据库中的主题序列作比对。比对结果是排列好的hit list,后面是一系列的单独的比对情况,以及不同的分值和统计值(如图7.9)。下文将会详细介绍选择不同的搜索程序、序列数据库和不同的参数都会对搜索产生影响,而且还有不同的界面,比如操作台命令、WWW形式和E-mail等。图7.10给出了一个使用Web界面进行数据库搜索的例子。这种形式的一个优点就是对任何一个感兴趣的比对,全部注解和文献应用都可以通过超文本简单方便地联接至原始的序列条目和相关的在线文献。

a

The best score are: initn initl opt z-sc E(59248)

gi|1706794|sp|P49789|FHIT_HUMAN FRAGILE HISTIDINE 996 996 996 1350.4 0

gi|1703339|sp|P49776|APH1_SCHPO BIS(5’-NUCLEOSYL) 431 395 395 536.2 2.8e-23

gi|1723425|sp|P49775|YD15_YEAST HYPOTHETICAL 24.8 290 171 316 428.1 2.9e-17

gi|1724021|sp|Q11066|YHIT_MYCTU HYPOTHETICAL 20.0 178 178 184 250.7 2.2e-07

gi|417124|sp|Q04344|HIT_YEAST HIT1 PROTEIN (ORF U 159 104 157 216.2 1.8e-05

gi|418447|sp|P32084|YHIT_SYNP7 HYPOTHETICAL 12.4 139 139 140 195.0 0.00028

gi|1351828|sp|P47378|YHIT_MYCGE HYPOTHETICAL 15.6 132 132 133 183.9 0.0012

à gi|1169826|sp|P43424|GAL7_RAT GALACTOSE-1-PHOSPHA 97 97 128 169.7 0.0072

gi|418446|sp|P32083|YHIT_MYCHR HYYPOTHETICAL 13.1 102 102 119 166.8 0.01

gi|1708543|sp|P49773|IPK1_HUMAN PROTEIN KINASE C 87 87 118 164.5 0.0014

gi|1724020|sp|P49774|YHIT_MYCLE HYPOTHETICAL 17.0 131 82 117 161.5 0.02

gi|1724019|sp|P53795|YHIT_CAEEL HYPOTHETICAL HIT- 98 98 116 161.5 0.02

gi|1170581|sp|P16436|IPK1_BOVIN PROTEIN KINASE C 86 86 115 160.4 0.023

gi|1730188|sp|Q03249|GAL7_MOUSE GALACTOSE-1-PHOSP 87 87 120 159.3 0.027

gi|1177047|sp|P42856|ZB14_MAIZE 14 KD ZINC-BIODIN 132 79 112 156.3 0.04

gi|1209081|sp|P07902|GAL7_HUMAN CALACTOSE-1-PHOSPH 78 78 117 154.8 0.048

gi|1177046|sp|P42855|ZB14_BRAJU 14 KD ZINC-BINDIN 115 76 110 154.5 0.05

gi|140775|sp|P26724|YHIT_AZOBR HYPOTHETICAL 13.2 115 65 109 152.6 0.064

gi|1169852|sp|P31764|GAL7_HAEIN GALACTOSE-1-PHOSP 62 62 104 137.9 0.42

gi|113999|sp|P16550|APA1_YEAST 5’,5’’’-P-1,P-4-TE 108 66 103 137.1 0.47

b

>>gi|1169826|sp|P43424|GAL7_RAT GALACTOSE-1-PHOSPHATE UR (379 aa)

initn: 97 init1: 97 opt: 128 z-score: 169.7 E(): 0.0072

Smith-Waterman score: 128; 30.8% identity in 107 aa overlap

10 20 30

FHIT MSFRFG-QHLIKPSVVFLKTELSFALVNRKPV

...: X.:.. . : .: ..:: :

GAL7 VWASNFLPDIAQREERSQQTYHNQHGKPLLLEYGHQELLRKERLVLTSEYWIVLVPFWAV

190 200 210 220 230 240

40 50 60 70 80

FHIT VPGHVLVCPLRPVERFHDLRPDEVADLFQTTQRVGTVVEKHFHGTSLTFSM—QDGP---

: ..:. : : :.:. .: : : :: .: ... : .. X. ::. .:: . .:

GAL7 WPFQTLLLPRRHVQRLPELTPAERDDLASTMKKLLTKYDNLFE-TSFPYSMGWHGAPMGL

250 260 270 280 290 300

90 100 110 120 130 140

FHIT EAGQTVKH--VHVHVLPRKAGDFHRNDSIYEELQKHDKEDFPASWRSEEEMAAEAAALRV

..: : : .:.: :

GAL7 KTGATCDHWQLHAHYYPPLLRSATVRKFMVGYEMLAQAQRDLTPEQAAERLRVLPEVHYC

310 320 330 340 350 360

图7.9:进行FASTA搜索的输出:(a)用人类组氨酸三联体蛋白作为(Swiss-Prot P.49789)查询序列,以Swissprot数据库为基础,进行FASTA搜索所得到的命中结果,在这个操作中,参数ktup=1;(b).以数据库中的一个条款(在命中列表中以箭头标出)为查询序列(其中包含老鼠的1-磷酸-半乳糖尿苷酸转移酶序列)所得到的最佳局部比对结果。虽然在这里,序列的相似性不太好,但是这些蛋白在结构上都显示了很好的相似性。

 

7.10:在WWW上进行数据库相似性搜索:NCBI数据库搜索的高级BLAST形式,在Web网页上容易实现。查询序列应该由剪切板中粘贴到最大的文本框中,(在本图中,框中显示的是U43746序列)。搜索中另外一些基本的元素包括搜索程序的名字以及数据库的名字,这两个元素都可以通过下拉框选择。如果需要的话,可以设定附加的选项参数。这里还有一个基本的BLAST形式,当然高级的选项参数被隐藏起来了。最后,简单地点击一下“Submit”键,提交请求后就可以开始搜索了。

 

如今的序列数据库非常之大,并且正以爆炸式的速度不断增长,在这种条件下,利用动态程序的方法直接进行数据库搜索已经变得不切实际。一个解决方法就是使用大型计算机和相关的特殊硬件,但是我们要讨论的目的是普通计算机能干些什么。当最佳方法不可行时,我们必须求助于那些启发式方法,这些方法充分利用了近似值以加快序列比较,但同时会在错过正确比对这一方面冒一点险。

有一种启发式方法建立在这样的策略之上,它将序列分解成由连续字母组成的短串(称为字串)。基于字的方法,在八十年代早期由Wilbur和Lipman提出,并且广泛使用于今天的搜索程序之中。其基本思想是这样的,一个能够揭示出正确的序列关系的比对至少包含一个两个序列都拥有的字串,把查询序列中的所有字串编成索引,并且在数据库扫描中查询这些索引,这些击中的字串就会很快被鉴定出来。

FASTA

FASTA程序是第一个广泛使用的数据库相似性搜索程序。为了达到较高的敏感程度,程序引用取代矩阵实行局部比对以获得最佳搜索。但众所周知,使用这种策略会非常耗费工作时,为了提高速度,在实施耗时的最佳搜索之前,程序使用已知的字串检索出可能的匹配。在速度和敏感度之间权衡选择依赖于ktup参数,它决定了字串的大小。增大ktup参数就会减少字串命中的数目,也就会减少所需要的最佳搜索的数目,提高搜索速度。缺省的ktup值在进行蛋白比较时选择2,但是在间距较大的情况下,将ktup值降为1较为理想。

FASTA程序并不会研究每一个遇到的字串命中,但在一开始会寻找包含若干个附近的命中的片段。使用启发式方法,这些片段会被赋予分值,最好的一个在输出时会显示为init1分值,这若干个片段会被组合起来,一个新的initn分值会从中计算出来。然后在最好的初始片段中局限于其对角线带上,会进行一次包含空位的局部比对以评估最可能的匹配。这个最佳比对的分值会在输出时显示为opt分值。对最后报导的比对来说,还要进行一次全程的Smith-Waterman比对。图7.9b显示了一个例子。对数据库中的每一个序列都只会由一个最佳的比对,但是,如果蛋白质中包含若干个模块,一些很有意义的比对就会被错过,匹配序列还必须由LALIGN程序作进一步分析。

从2.0版本开始,FASTA对每一个检索到的比对都提供一个统计学显著性的评估。程序为随机分值假定了一个极值分布,但是改写了概率密度函数的形式,其中预期的分值与数据库中的序列长度的自然对数呈线形关系,这样,可以使用简单的线形回归函数计算常规的比对的z值。最后,计算出预期的E值,从而给出那些z值不小于已知值的随机比对的预期数目。

BLAST

BLAST程序对数据库搜索进行了大量的改良,提高了搜索速度,同时把数据库搜索建立在了严格的统计学基础之上。但是,为了达到这一目的,仍然需要权衡选择,也就是说,局部比对的限制条件可能不包括空位。这个限制条件对应用Karlin-Altschul统计学极为有利,另一方面,既然空位没有明确地放在模型中,结果就不会象人们期望的那样接近于预期的比对。这并不是说插入和确实会妨碍匹配,在大多数情况下,比对仅仅会被分解为若干个明显的HSPs。无论如何,老版本的BLAST程序(1.4以前)的局限性在新版本中已经被消除了,新版本在对待空位问题上有着明确的作法(在下面讨论)。

对于一个即将被BLAST程序报告的比对,其中必然包含一个HSP,其分值不小于终止值S。这个终止值因人而异,但是使用时是很难知道其合适值的。因为程序基于Karlin-Altschul统计学,人们可以指明一个预期的终止E值,然后软件会在考虑搜索背景的性质的基础上(比如数据库的大小,取代矩阵的性质)计算出正确的S值。BLAST的一项创新就是邻近字串的思想。这个协定不需要字串确切地匹配,在引入取代矩阵的情况下,当主题序列中的字串有一个最低分值T时,BLAST就宣布找到了一个命中的字串。这个策略允许较长字串长度(W)(为了提高速度),而忽略了敏感度。于是,T值称为制衡速度和敏感度的临界参数,而W是很少会变化的。如果T值增大,可能的命中字串的数目就会下降,程序执行就会加快,减小T值会发现较远的关系。

发生一个字串命中后,程序会进行没有空位的局部寻优,比对的最低分值是S。将比对同时向左方和右方延伸并将分值加和就会得到结果。当遭遇一系列的最低分值时,加和的分值就会下降,这时,分值就不再可能反弹回S值。这个发现为附加的启发式知识提供了依据,因此,当分值的降低(与遭遇的最大值相比)超过分值下降阈值X时,命中的延伸就会终止。于是,系统回减少毫无指望的命中延伸,继续进行其它操作。

使用BLAST

可以通过e-Mail、WWW或控制台命令操作BLAST程序,无论如何,一次数据库搜索包括四种基本元素:BLAST程序的名称,数据库名称,查询序列和大量的合适的参数,很显然,当以上元素发生变化时,搜索的细节就会随之改变。为了避免混淆,我们把BLAST功能性描述为普通名词,避免提及专有工具。读者可能会要参考使用到的专有工具的有关内容。要得到关于用e-Mail执行BLAST搜索的介绍,给blast@ncbi.nlm.nih.gov发一封含有“HELP”的邮件;在WWW工具中,帮助是在线的;如果使用Unix系统,使用man blast可以获得详细的帮助信息。

表7.1、BLAST程序:

程序 数据库 查询 内容
Blastp 蛋白质 蛋白质 使用取代矩阵寻找较远的关系:可以进行SEG过滤。
Blastn 核苷酸 核苷酸 寻找较高分值的匹配,对较远关系不太适用。
Blastx 核苷酸(翻译) 蛋白质 对于新的DNA序列和ESTs的分析极为有用。
Tblastn 蛋白质 核苷酸(翻译) 对于寻找数据库中没有标注的编码区极为有用。
tblastx 核苷酸(翻译) 核苷酸(翻译) 对于分析EST极为有用。

几种不同的BLAST可以通过查询序列和数据库序列的类型来加以区分:blastp比较的是查询蛋白同蛋白质数据库;相应于核酸序列的程序是blastn;如果序列类型不同,DNA序列可以被翻译成蛋白序列(所有六种阅读框架)后同蛋白序列进行比较,blastx比较一个DNA的查询序列同一个蛋白质序列库,其结果对分析新序列和ESTs很有用;对于一个基于核酸序列库的蛋白质查询,tblastn程序对于寻找数据库中序列的新的编码区很有用;最后一个只在特殊情况下使用(在这里介绍只是出于完整的考虑),tblastx将DNA查询序列和核酸序列库中的序列全部翻译成蛋白质序列,然后进行蛋白质序列比较,这个程序主要应用于ESTs比较,尤其是当人们怀疑到其中有可能的编码区,即使并没有确切地发现这一区域。

所有这些程序使用服务器上的序列数据库,从而不需要本地的数据库,表7.2和7.3陈列了一些BLAST使用的蛋白质和核酸的序列数据库。对于常规的搜索,nr数据库拥有大量的氨基酸和核酸序列,同时合并相同的序列以减少冗余度。为了检测在过去30天里提出或更新的序列,提供了一个称为“month”的数据库。不管是nr还是month,都是日日更新。表7.2和7.3中列出的其它一些数据库在一些特别的环境里十分有用,比如在比较模型物种(酵母和大肠杆菌)的全序列时,搜索特别类型的序列(dbest或dbsts),或检测是否存在污染或问题序列(vector,alu或mito)。

表7.2、使用BLAST的蛋白序列数据库:

数据库 描述
Nr 融合了Swiss-Prot,PIR,PRF以及从GenBank序列编码区中得到的蛋白质和PDB中拥有原子坐标的蛋白质,绝非多余。
Month Nr的字集,每月(30天)更新,搜集了过去30天中的最新序列。
Swissprot Swiss-Prot数据库。
Pdb 拥有三维空间结构的原子坐标的氨基酸序列库。
Yeast 酵母基因组中基因编码的全套蛋白质。
ecoli 有大肠杆菌基因组中基因编码的全套蛋白质。

表7.3、使用BLAST的核苷酸序列数据库:

数据库 描述
Nr 极有价值的GenBank,排除了EST,STS和GSS部分。
Month Nr的字集,每月(30天)更新,搜集了过去30天中的最新序列。
Est Genbank中的EST部分(expressed sequence tags, 表达序列标签)。
Sts Genbank中的STS部分 (sequence tagged sites, 序列标签位点)。
Htgs Genbank中的HTG部分 (high throughput genomic sequences, 高容量基因组序列)。
Gss GenbankGSS(genome survey sequences,基因组测定序列)。
Yeast 酵母的全基因组序列。
Ecoli 大肠杆菌的全基因组序列。
Mito 脊椎动物线粒体的全基因组序列。
Alu 搜集了灵长类动物的Alu重复序列。
vector 搜集了流行的带菌体的克隆。

一个BLAST搜索的例子会介绍搜索输出的不同元素。如图7.11所示的例子,一种Alzheimer疾病感受性蛋白质的氨基酸序列(由GenBank中L43964翻译)作为查询序列同dbest数据库用tblastn进行搜索。进行这么一次搜索的目的是要鉴定模型生物中可能的同源物的cDNA克隆,从而为在人类中无法进行的实验打开方便之门(相应于EST序列的克隆是已经实现的)。数据库中的每一个EST序列在同alzheimer蛋白质序列比较以前,都已经按照所有的阅读框架得到翻译。图7.11a显示了此次搜索得到部分命中的列表,前两列给出了每一个显著性匹配的序列的标识和描述。尽管浏览时定义被缩短了,我们仍然可以看到老鼠和果蝇的序列都被包含进来了。下一列给出了得到最佳HSP(即使其它阅读框架翻译结果也会达到命中)的阅读框架。后面三列给出了最佳HSP的分值、p值总和及p值计算时使用到的HSP数目。

包含一种果蝇EST(由箭头标出)的比对在图7.11b中得以显示。其中包含了两个HSP,并且显示了每一个的分值,EST的概念性翻译同查询序列并排显示。相同的氨基酸残基在两个序列之间回显,+表示两个不同残基匹配的分值是正数(比如保守取代)。从不同阅读框架得到的两个HSP是显著的并且彼此相邻,这一点从序列坐标就可以看出来。这种形式表示EST序列的一种阅读框架是错误的,并且对于用相对容错性的工具进行序列单向通行数据分析时极为有效。

a

sum

Reading High Probability Y

sequence producing High-scoring Segment Pairs: Frame Score P(N) N

gb|AA056325|AA056325 zf53a03.sl Soarea retina N2b4HR H... +3 724 3.4e-102 2

gb|T03796|T03796 IBIB913 Infant brain,Bento Soares...+3 567 2.6e-78 2

gb|AA260597|AA260597 mx76g09.r1 Soares mouse NML Mus m...+2 239 4.9e-53 4

gb|H86456|H86456 yt01b06.s1 Homo sapiens cDNA clon...+2 323 4.3e-52 4

gb|N24576|N24576 yx72a04.s1 Homo sapiens cDNA clon...+1 365 5.5e-47 2

gb|AA265273|AA265273 mx91c12.r1 Soares mouse NML Mus m...+2 239 6.4e-41 2

gb|AA237206|AA237206 mx18e01.r1 Soares mouse NML Mus m...+3 159 1.5e-40 3

gb|R146001|R146001 yf34b10.r1 Homo sapiens cDNA clon...+1 278 1.5e-40 2

gb|AA200706|AA200706 mu03f12.r1 Soares mouse 3NbMs Mus...+1 343 1.9e-40 1

gb|AA045064|AA045064 zk77f12.s1 Soares pregnant ulerus...-3 269 2.3e-37 2

gb|AA087434|AA087434 mm28a04.r1 Stratagene mouse skin....+3 322 3.6e-37 1

gb|R05907|R05907 ye93h02.r1 Homo sapiens cDNA clon...+3 252 7.7e-37 2

gb|AA268820|AA268820 vb01c10.r1 Soares mouse NML Mus m...+2 234 7.7e-35 2

gb|AA162310|AA162310 mn44a07.r1 Beddington mouse embry...+1 134 8.3e-34 3

gb|N27820|N27820 yx54h10.r1 Homo sapiens cDNA clon...+3 154 7.8e-29 2

gb|AA234907|AA234907 zs38f03.r1 Soares NhHMPu S1 Homo... +2 155 1.8e-28 2

gb|AA231081|AA231081 mw11d11.r1 Soares mouse 3NME12 5... +3 134 8.8e-23 2

gb|H91652|H91652 ys80c04.s1 Homo sapiens cDNA clon... -3 215 3.7e-22 1

gb|H50532|H50532 yo30h08.s1 Homo sapiens cDNA clon... -2 211 1.2e-21 1

gb|AA150236|AA150236 zl03c01.r1 Soares pregnant uterus...+1 159 5.0e-21 2

gb|AA144382|AA144382 mr15d12.r1 Soares mouse 3NbMS Mus...+3 159 7.6e-21 2

à gb|AA390557|AA390557 LD09473.5prime LD Drosophila Embr...+3 130 1.6e-20 2

gb|AA210480|AA210480 mo86b03.r1 Beddington mouse embry...+2 128 2.0e-20 3

gb|H19021|H19021 ym44b02.r1 Homo sapeins cDNA clon...+2 134 5.9e-20 2

gb|AA283084|AA283084 zt14g09.s1 Soares NbHTGBC Homo sa...-3 175 2.3e-19 2

gb|H25759|H25795 y149d01.s1 Homo sapiens cDNA clon...-2 185 5.0e-18 1

gb|H33787|H33787 EST110123 Rattus sp.cDNA 5’ end..... +1 137 6.7e-17 2

gb|AA201988|AA201988 LD05058.5prime LD Drosophila Embr...+3 175 5.5e-15 1

gb|AA263526|AA263526 LD06652.5prime LD Drosophila Embr...+1 167 7.0e-14 1

gb|R46340|R46340 yj52c04.sl Homo sapiens cDNA clon...-1 151 5.6e-13 1

gb|AA246675|AA246675 LD05588.5prime LD Drosophila Embr...+2 117 2.8e-10 2

gb|AA282899|AA282899 zt14g09.r1 Soares NbHTGBC Homo sa...+3 118 6.1e-07 1

gb|AA247705|AA247705 csh0941.seq.F Human fetal heart,....+3 56 0.0039 2

 

b

gb|AA390557|AA390557 LD09473.5prime LD Drosophila Embryo Drosophila

melanogaster cDNA clone LD09473 5’

Length – 659

Score – 130 (60.4 bits), Expect – 1.6e-20, Sum P(2) – 1.6e-20

Identities – 25/60 (41%), Positives – 40/60 (66%), Frame - +3

Query: 105 TIKSVRFYTEKNGQLIYTTFTEDTPSVGQRLLNSVLNTLIMISVIVVMTIFLVVLYKYRC 164

+I S+ FY + L+YT F E +P + +++ ++LI++SV+VVMT L+VLYK RC

sbjct: 480 SINSISFYNSTDVYLLYTPFHEQSPEPSVKFWSALGSSLILMSVVVVMTFLLIVLYKKRC 659

Score – 117 (54.3 bits), Expect – 1.6e-20, Sum P(2) – 1.6e-20

Identities –23/30 (76%), Positives – 27/30 (90%), Frame - +1

Query: 75 LEEELTLKYGAKHVIMLFVPVTLCMIVVVA 104

+EEE LKYGA+HVI LFVPV+LCM+VVVA

sbjct: 391 MEEEQGLKYGAQHVIKLFVPVSLCMLVVVA 480

图7.11、一次TBLASTN搜索的输出:在这次TBLASTN搜索中,以dbest数据库为基础,以阿尔茨海默氏病(即进行性老年性痴呆)基因(Genbank 检索号码L43964)的蛋白质产物为查询序列,目的是为了从其它那些可能同人类基因有同源性的物种中鉴定出一些cDNA克隆。(a).命中列表的一部分显示了其中最好的25个命中。每个检索出来的序列都由它们的GenBank检索号码以及一部分定义行组成。其中包括了它们的阅读框架和最佳HSP分值,同时显示的还有一个偶然命中的可能性的加和。最后一列中的数据给出了在计算加和的可能性时所涉及到的HSP的数量。在这个命中列表中可以见到至少10条从老鼠中得到的序列和一条从果蝇中得到的序列; (b).同果蝇的EST序列(GenBank AA390557)理论上的翻译序列匹配的结果。找到了两个HSPs,每一个使用不同的阅读框架。相同的残基在两行序列中间的相应位置回显,而“+”符号标记着那些不相同但是其取代分值是正分的残基。

BLAST的最新改进

最近发布的BLAST程序的修订版提高了搜索速度、敏感度和实用性。这个完全重新写过的软件包指定为2.0版本(避免同WU-BLUST混淆,这个软件是由华盛顿大学设计的,有时称为BLAST2)。应该注意到,在发布的2.0版本中,命令行的参数有很大改变,其中一些常用的参数列在表7.4中。

一个改进来自于引发一个字串命中的延伸的标准。现在,在一个需要考虑的残基的窗口里必须找到两个字串命中。使用这种策略提高了搜索速度,因为大量随机的字串命中将会被忽略,并且很有可能得到一个显著性良好的比对。第二个改进是能够明确地而不是含蓄地处理空位。除了帮助使用者更加容易地理解产生的比对,新版本还提高了较远关系的敏感性,其中可能会包含许多插入和缺失。比较从寻找无空位的HSP这一标准策略开始,然后,这一比对中获得最高分区域的中心一列被鉴定出来,接着,从这一点向前和向后延伸,通过赋值的路径进行无空位局部比对的搜索。如同最初的HSP搜索,一个分值下降的阈值X将会促使放弃那些遭遇大量负的取代分值的路径。对剩余的HSP进行反复的这种操作,将会揭示另外的含空位的比对,并保证它们同已经报告的部分不会相交。这个系统不同于FASTA所采取的策略,FASTA只会产生一个最佳的比对。

表7.4、一些对于BLAST很有用的参数值:

参数名称 BLAST 1.4 BLAST 2.0
数据库 (database) 第一参数 -d database
查询序列文件 (query sequence file) 第二参数 -I filename
期望阈值E (expectation cutoff) E = number -e number
HSP分值阈值S (HSP score cutoff) S = number -s number
字串分值阈值T (word score cutoff) T = number -f number
多命中窗口A (multihit window) n/a -A number
打分矩阵 (score matrix) -matrix matrix -M matrix
低复杂度过滤 (low-complexity filtering) -filter seg -F
空位开放罚分 (gap opening penalty) n/a -G number
空位拓展罚分 (gap extension penalty) n/a -E number
PSI-BLAST反复 (PSI-BLAST iterations) n/a -j number

对于那些弱势的但是显著性较强的比对,进行较高敏感性的数据库搜索的一个方法就是使用诸如profile(表头轮廓)的数据结构(Gonzalez et al., 1994)。这个策略可能曾经被认为是个进行数据库搜索的比较先进的课题,但是BLAST的一个新特性简化了基于profile的搜索工作。一个profile可能会被理解为一个列表,其中列出了在一个保守的蛋白质结构域中每一个位点发现每一种氨基酸残基的频率。建立一个profile可能是很乏味的,其信息是从那些拥有我们感兴趣的蛋白质结构域的多序列比对中得到的,这些比对必须预先准备好,而且,在这里有许多技术上的问题还没有解决。

位点特性反复BLAST(PSI-BLAST)是指BLAST2.0的一个特性,其中一个profile被不断组织并且不断精练。这个过程开始于使用一个简单查询序列的一个标准的数据库搜索。在这个初始的搜索结果中,一个profile从高度显著的比对中获得,然后这个profile在第二轮的数据库搜索中使用。如果需要的话,这个过程会反复进行,并且在操作中为了精练profile,会在每一轮中加入新的序列。

为了演示PSI-BLAST方法的高敏感性,旦氨酸三联体蛋白(HIT)序列被用来作为数据库搜索中的查询序列。HIT和1-磷酸乳糖尿苷酸转移酶(GalT)基于位点重叠的三位结构相似性最近得到描述(Holm and Sander, 1997)。经过一次标准的(一轮)BLASTP搜索,没有发现一个对GalT序列有显著的命中。但是经过多次搜索,在每一次反复中都发现新的关系,正如图7.12所示。在第二次搜索中了发现老鼠的GalT蛋白质,并且在这一信息被加入profile之后,另外一些其它物种的同源物也被检测出来。

Sequences producing significant alignments: Hign E

Score Value

Pass1:

sp|P49789|FHIT_HUMAN FRAGILE HISTIDINE TRIAD PROTEIN 290 7e-79

sp|P49776|APH1_SCHPO BIS(5’ – NUCLEOSYL) – TETRAPHOSPHATASE (ASYMME... 117 8e-27

sp|P49775|YD15_YEAST HYPOTHETICAL 24.8 KD HIT – LIKE PROTEIN 88.0 6e-18

sp|Q11066|YHIT_MYCTU HYPOTHETICAL 20.0 KD HIT – LIKE PROTEIN 52.7 3e-07

sp|Q04344|HIT_YEAST HIT1 PROTEIN (ORF U) 45.3 4e-05

Pass2:

sp|P47378|YHIT_MYCGE HYPOTHETICAL 15.6 KD HIT – LIKE PROTEIN 70.5 1e-12

sp|P32083|YHIT_MYCHR HYPOTHETICAL 13.1 KD HIT – LIKE PROTEIN IN P... 59.0 3e-09

sp|P26724|YHIT_AZOBR HYPOTHETICAL 13.2 KD HIT – LIKE PROTEIN IN H... 57.6 9e-09

sp|P32084|YHIT_SYNP7 HYPOTHETICAL 12.4 KD HIT – LIKE PROTEIN IN P... 55.7 3e-08

sp|P53795|YHIT_CAEEL HYPOTHETICAL HIT – LIKE PROTEIN F21C3.3 54.3 9e-08

sp|P42856|ZB14_MAIZE 14 KD ZINC – BINDING PROTEIN (PROTEIN KINASE... 52.8 2e-07

sp|P42855|ZB14_BRAJU 14 KD ZINC – BINDING PROTEIN (PROTEIN KINASE... 50.2 1e-06

sp|P49774|YHIT_MYCLE HYPOTHETICAL 17.0 KD PROTEIN HIT – LIKE PROT... 49.5 2e-06

sp|P49773|IPK1_HUMAN PROTEIN KINASE C INHIBITOR 1 (PKCI – 1) 49.1 3e-06

sp|P16436|IPK1_BOVIN PROTEIN KINASE C INHIBITOR 1 (PKCI – 1) (17 ... 48.7 4e-06

sp|P44956|YCFF_HAEIN HYPOTHETICAL HIT – LIKE PROTEIN HI0961 47.3 1e-05

sp|P43424|GAL7_RAT GALACTOSE – 1 – PHOSPHATE URIDYLYLTRANSFERASE 41.0 8e-04

Pass3:

sp|Q03249|GAL7_MOUSE GALACTOSE – 1 – PHOSPHATE URIDYLYLTRANSFERASE 87.2 1e-17

sp|P07902|GAL7_HUMAN GALACTOSE – 1 – PHOSPHATE URIDYLYLTRANSFERASE 79.8 2e-15

sp|P31764|GAL7_HAEIN GALACTOSE – 1 – PHOSPHATE URIDYLYLTRANSFERASE 64.7 6e-11

sp|P09148|GAL7_ECOLI GALACTOSE – 1 – PHOSPHATE URIDYLYLTRANSFERASE 62.5 3e-10

sp|P22714|GAL7_SALTY GALACTOSE – 1 – PHOSPHATE URIDYLYLTRANSFERASE 58.1 6e-09

sp|P09580|GAL7_KLULA GALACTOSE – 1 – PHOSPHATE URIDYLYLTRANSFERASE 48.5 4e-06

sp|P08431|GAL7_YEAST GALACTOSE – 1 – PHOSPHATE URIDYLYLTRANSFERASE 40.8 0.001

Pass4:

sp|P40908|GAL7_CRYNE GALACTOSE – 1 – PHOSPHATE URIDYLYLTRANSFERASE 71.0 8e-13

sp|P13212|GAL7_STRLI GALACTOSE – 1 – PHOSPHATE URIDYLYLTRANSFERASE 57.0 1e-08

图7.12、使用PSI-BLAST后,敏感性提高很大:在这次BLASTP搜索中,查询序列是人类组氨酸三联体(HIT)蛋白(Swiss-Prot P49789),搜索时开启了PSI-BLAST功能。在每一次重复搜索中,新检索出来的具有统计学显著性的匹配都会显示它们的定义行,打分值以及E 数值。

 低复杂度区域

不管是蛋白还是核酸都包含一些偏颇的区域,在进行序列数据库搜索时这些区域可能会导致一些令人迷惑的结果。这些低复杂度区域(LCRs)在从明显的同性聚合顺串和短周期重复到更精细的情况(如其中某些或一些残基过多表现)的范围内变化。一个称为SEG的程序发展起来,目的是要把一个蛋白质序列分解为低复杂度和高复杂度组成的各个片段(Wootton and Federhen, 1993, 1996)。这个程序的结果表明数据库中的蛋白质有一半以上拥有至少一个LCR(Wootton and Federhen, 1993; Wootton, 1994)。LCRs的进化、功能和结构性质并没有被很好地了解。在DNA中,有许多种简单的重复,其中一些已经知道是高度多样性的,并且在作基因图谱时经常使用的。它们源起的机制可能是聚合酶滑动、偏颇核苷酸取代或者不等交换。LCRs更偏好于在结构上以非球形区域的形式存在,那些在物理化学上已经被定义为非球形的区域通常可以在使用SEG程序时获得较好的结果(Wootton, 1994)。

对于包含LCR的序列进行比对是成问题的,因为这些序列不符合残基-残基序列守恒的模型。有些时候,与功能相关的属性可能仅仅是周期性或组成结构,而不是任何特异的序列。而且,对比对作统计学显著性分析的方法是建立在一定的随机概念基础上的,LCR显然不符合这一条件,因此,对于一个包含LCR的查询序列,在进行数据库搜索的输出里会发现很多不正确的条目,因为这些匹配的显著性被过高评价了(Altschul et al., 1994)。这个问题大体上可以通过过滤(或者叫屏蔽)解决,操作是这样的,把有问题的子序列转化为不明确的字符(蛋白质用X,核酸序列用N),这样它们就不会对比对贡献正分了。

果蝇鳞甲基因产物的人类同源物就是包含LCR蛋白质的一个好例子,在用SEG分析的时候,两个低组成复杂度的序列区域被鉴定出来。图7.13a显示了缺省的树输出,其中低复杂度序列用小写字母表示在左边,高复杂度序列在右边用大写字母表示。第一个区域片段有61个残基,包含大量丙氨酸(alanine)和谷氨酸盐(glutamine)的多聚物;第二个区域片段有14个残基,偏向于精氨酸(arginine)。如果不进行过滤的话,许多包含这种偏向性序列的数据库序列都会被报告出来。使用命令行选项,SEG程序就会产生一个过滤后的查询序列版本。另外,过滤可以有BLAST程序自动完成,如果使用合适的参数。请注意在使用BLAST时,缺省情况下就可以实行过滤(比如在WWW版本)。这就解释了为什么查询序列中的不明确的字符串(在原序列中没有出现)会在比对中被偶然发现。

a

>gi|1703441|sp|P50553|ASH1_HUMAN ACHAETE – SCUTE HOMOLOG 1

1-11 MESSAKMESGG

agqqpqpqpqqpflppaacffataaaaaaa 12-72

aaaaaaqsaqqqqqqqqqqqqqqapqlrpaa
1.    DGQPSGGGHKSAPKQVKRQRSSSPELMRCKRRLNFSGFGYSLPQQQP

aavarrnerernrv 120-133
1.   KLVNLGFATLREHVPNGAANKKMSKVETLRSAVEYIRALQQLLDEHDAVSAAFQAGVLSP

TISPNYSNDLNSMAGSPVSSYSSDEGSYDPLSPEEQELLDFTBWF

b

>gi|1703441|sp|P50553|ASH1_HUMAN ACHAETE – SCUTE HOMOLOG 1

MESSAKMESGGXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

XXXXXXXXXXXXDGQPSGGGHKSAPKQVKRQRSSSPELMRCKRRLNFSGFGYSLPQQQPX

XXXXXXXXXXXXXKLVNLGFATLREHVPNGAANKKMSKVETLRSAVEYIRALQQLLDEHD

AVSAAFQAGVLSPTISPNYSNDLNSMAGSPVSSYSSDEGSYDPLSPEEQELLDFTBWF

c

>gi|540240 (U14590) achaete – scute homolog b [ Danio rerio ]

Length – 195

Score – 193 bits (512), Expect – 7e-49

Identities – 107/155 (69%), Positives – 118/155 (76%)

Gaps – 8/155 (5%)

QUERY 86 KQVKRQRSSSPELMRCKRRLNFSGFGYSLPQQQPXXXXXXXXXXXXXXKLVNLGFATLRE 145

K +KRQRSSSPEL+RCKRRL F+G GY++PQQQP K VN+GF TLR+

540240 32 KVLKRQRSSSPELLRCKRRLTFNGLGYTIPQQQPMAVARRNERERNRVKQVNMGFQTLRQ 91

QUERY 146 HVPNGAANKKMSKVETLRSAVEYIRALQQLLDEHDAVSAAFQAGVLSPTISPNYSNDLNS 205

HVPNGAANKKMSKVETLRSAVEYIRALQQLLDEHDAVSA Q GV SP++S YS

540240 92 HVPNGAANKKMSKVETLRSAVEYIRALQQLLDEHDAVSAVLQCGVPSPSVSNAYS----- 146

QUERY 206 MAG—SPVSSYSSDEGSYDPLSPEEQELLDFTNWF 238

AG SP S+YSSDEGSY+ LS EEQELLDFT WF

540240 147 -AGPESPHSAYSSDEGSYEHLSSEEQELLDFTTWF 180

图7.13、使用SEG程序检索低复杂度区域:使用SEG程序对人类achaete-scute蛋白(Swiss-Prot P50553)进行分析,发现了两段低复杂度区域。(a).以缺省的“tree”格式执行程序得到的输出结果,左边用小写字母显示了低复杂度区域,右边用大写字母显示了高复杂度区域。 (b) .开启-x命令行开关,SEG程序将会产生把低复杂度区域屏蔽掉的序列结果。 (c).为了方便使用,操作者可以使用BLAST程序来进行低复杂度区域的屏蔽。当一个低复杂度区域被屏蔽掉的序列作为查询序列被提交给数据库进行检索时,在BLASTP输出结果的比对中可能也会包括一些被屏蔽的分段序列。

重复元件

如果查询中包括一个重复元件的序列-比如说一个Alu重复-可能会出现许多错误的和令人费解的结果。虽然在蛋白质-蛋白质搜索中,这一般不会成为什么大问题,但是在包含DNA序列任何比较中,都必须对此引起必要的重视。基因组序列可能会包含大量分散的重复序列,特别是一些多基因族(例如Alus, LINEs和人的序列中的MERs),甚至mRNA序列中也可能含有重复序列,几乎都是信息的非翻译区。因此,重复元件在数据库序列中非常普遍,如果查询序列中也有这些重复,就会在比对中出现大量不正确的正分。虽然重复元件显示了大量不同成分,仍然有足够的相似性使比对具有一定的高显著性。虽然比对会跨越这些重复而不是侧面的单一序列,但是直接从数据库搜索的输出结果观察,这并不是显而易见的。

GenBank和Swiss-Prot数据库中都包含一些“暖序列(warming sequence)”,这些数据向使用者指出查询中包含重复序列(Claverie and Makalowski, 1993)。在GenBank中,这些条目表示了人类Alu重复的不同亚科的一致序列;在Swiss-Prot中的类似条目是Alu序列的六种翻译框架(一个接着一个,中间由若干X分隔)。在两种情况下,单词“WARNING”在定义行中非常显著。暖序列不必出现在命中列表的上方,而且,可以有许多包含Alu重复的数据库序列同查询序列非常相似,甚至比查询序列同暖序列还相似。这在图7.14a中有所体现,它显示了对人类转录因子CBFB(在3’UTR包含一个Alu)基于nr数据库进行一次blastn搜索的一部分命中。暖序列(用箭头标出)位于命中列表的第31位。虽然列表顶部的一些匹配显示了真正的关系(第一个是一个自命中),绝大多数只是因为具有Alu重复才会出现错误的正分。

在查询中更直接地检测Alu重复是否存在的方法就是在查询前先对alu数据库做一次搜索。如图7.14b所示,做完这个以后,包含alu的暖序列作为最高分匹配被报告出来。如果查询序列被发现包含重复元件,接下来的行动就是要对这个序列进行编辑改动,把它剔除或者屏蔽掉。在这里一个有用的工具就是CENSOR,它能够自动检测并且消除重复元件。

a

Smallest

Sum

High Probability Y

Sequences producing High – scoring Segment Pairs: Score P(N) N

gb|L20298|HUMCBFB Homo sapiens transcription factor... 8691 0.0 2

dbj|D14571|MUSPEBP2B2 Mouse mRNA for PEBP2B2 protein, co.. 2574 0.0 25

gb|L032791|MUSP215CBF Mus musculus core – binding factor m 2574 0.0 25

dbj|D14572|MUSPEBP281 Mouse mRNA for PEBP2B1 protein, co.. 2130 0.0 26

dbj|d14570|muspebp283 Mouse mRNA for PEBP2B4 protein, co.. 1701 0.0 26

gb|L03305|MUSCBFAA Mus musculus core – binding factor m 942 0.0 27

gb|L03306|MUSCBFAB Mus musculus core – binding factor m 2130 1.6e-282 10

gb|U22177|DMU22177 Drosophila melanogaster Big brothe... 382 1.5e-37 2

emb|Y10196|HSPEX H.sapins PEX gene 400 4.4e-22 1

gb|L77570|HMUDGCRCEN Homo sapiens DiGeorge syndrome cri... 409 6.7e-22 2

gb|AD00067|1010603 Homo sapiens DNA from chromosome 1... 392 2.0e-21 1

emb|Z83822|HS306D1 Human DNA sequence from PAC 306D1 ... 392 2.0e-21 1

emb|Z82097|HSF77D12 Human DNA sequence from fosmid F77... 391 2.5e-21 1

dbj|D42052|HUMKIAA000 Human cosmid Q7A10 (D21S246) inser... 391 2.5e-21 1

gb|U83511|HSUB3511 Human Xp22 cosmids U177G4,U152H5, ... 386 6.5e-21 1

gb|U52112|HSU52112 Human Xq28 genomic DNA in the regi... 386 6.5e-21 1

gb|S83170|S83170 tissue – type plasminogen activator.. 382 1.1e-20 1

emb|X9642|HSCAMF3X1 H.sapiens Y chromosome cosmid CAMF... 383 1.1e-20 1

gb|U95739|HSU95739 Human chromosome 16p11.2 – p12 BAC c. 383 1.1e-20 1

gb|95743|HSU95743 Human chromosome 16p13.1 BAC clone... 383 1.1e-20 1

gb|U91322|HSU91322 Human chromosome 16p3 BAC clone C.... 383 1.1e-20 1

gb|U82609|HSU82609 Human centromere – specific histone.. 382 1.3e-20 1

gb|AC001061|HSAC001061 Homo sapiens (subclone 2_g6 fromP.... 382 1.3e-20 1

emb|Z46940|HSPRMTNP2 H.sapiens PRM1 gene, PRM2 gene and... 382 1.4e-20 1

gb|K03021|HUMTPA Human tissue plasminogen activator... 382 1.4e-20 1

gb|U15422|HSU15422 Human protamine 1 (PRM1), protamin... 382 1.4e-20 1

gb|U91323|HSU91323 Human chromosome 16p13 BAC clone C... 382 1.4e-20 1

emb|Z54147|HSLI29H7A Human DNA sequence from cosmid L12... 381 1.7e-20 1

emb|Z82194|HSJ272J12 Human DNA sequence fom clone J272J12 374 1.7e-20 2

dbj|D0035|HIV2CAM2 Human immunodeficiency virus type-... 380 2.0e-20 1

à gb|U14567|HSU14567 ***ALU WARNING: Human Alu_J subfam... 373 2.4e-20 1

gb|L81578|HSL81578 Homo sapiens (subclone 2_b2 from P... 386 3.0e-20 2

gb|L81854|HSL81854 Homo sapiens (subclone 2_b8 from P... 377 3.4e-20 1

b

Smallest

Sum

High Probability Y

Sequences producing High – scoring Segment Pairs: Score P(N) N

à lcl|HSU14567 ***ALU WARNING: Human Alu – J subfamil... 373 4.1e-24 1

lcl|unknown gb|M94643_HSAL001949 349 1.4e-22 1

lcl|HSU14574 ***ALU WARNING: Human Alu – Sx subfami... 347 7.0e-22 1

lcl|HSU14573 ***ALU WARNING: Human Alu – Sq subfami... 347 7.0e-22 1

lcl|unknown gb|Z15026_HSAL001005 (Alu – J) 324 1.4e-21 1

lcl|unknown gb|M15657_HSAL001254 (Alu – J) 337 6.3e-21 1

lcl|unknown gb|M61839_HSAL002304 (Alu – J) 314 6.6e-21 1

lcl|unknown gb|X17354_HSAL000525 (Alu – J) 329 6.6e-21 1

lcl|HSU14572 ***ALU WARNING: Human Alu – Sp subfami... 329 2.4e-20 1

lcl|unknown gb|J03619_HSAL001939 (Alu – Sx) 329 2.8e-20 1

lcl|unknown gb|L11910_HSAL002838 (Alu – J) 307 2.8e-20 1

lcl|unknown gb|M11228_HSAL002744 (Alu – Sp) 329 2.9e-20 1

lcl|unknown gb|L18035_HSAL004322 (Alu – J) 318 9.3e-20 1

lcl|unknown gb|L05367_HSAL002551 (Alu – J) 318 1.0e-19 1

lcl|unknown gb|M58600_HSAL002004 (Alu – J) 322 1.2e-19 1

lcl|unknown gb|Z23796_HSAL005276 (Alu – J) 306 1.7e-19 1

lcl|unknown gb|M90058_HSAL002955 (Alu – J) 294 2.5e-19 1

lcl|unknown gb|D14642_HSAL003786 (Alu – J) 315 4.0e-19 1

lcl|unknown gb|M29038_HSAL002942 (Alu – J) 314 5.5e-19 1

lcl|unknown gb|M92357_HSAL001387 (Alu – J) 310 9.8e-19 1

图7.14、反复元件可能会导致令人迷惑的结果:本次blastn查询使用的查询序列是人类转录因子CBFB(GenBank L20298)的cDNA序列。(a).如果使用nr数据库,最先的一些匹配同查询序列具有真正的关联,但是也会报告许多不正确的命中结果,这些命中分布于各个人类染色体的基因组区域。在这个命中列表中,打箭头处(位于第31行)的一致的Alu-J序列被列为警告序列。 (b).如果使用alu数据库,Alu-J警告序列就成了最佳匹配序列。

为了鉴定这些潜在的搜索成果,学会怎样评估搜索结果非常重要。上述的一些策略只应用于Alu反复,它是人类以及其它一些物种中出现频率最高的,但是其它一些反复仍然存在,虽然含量较低,而且,其它物种会显示出完全不同类型的反复元件。现在有一个数据库搜索输出的附加性质,它可以指示出反复元件。例如,注意比对中与DNA序列编码区域相关的位点是非常有益的。如果非编码区域匹配而编码区域不匹配,那么反复序列就很令人怀疑;如果查询序列同大量序列匹配,但是这些序列相互之间没有什么关系,但是比对的分值都很相近,这样的结果就极为可疑。例如图7.14a中,许多匹配的相似性分值都几乎一样,而且包括了从若干不同的人类染色体上来的质粒。虽然对这个发现有很多解释,但是一个明智的看法就是至少承认这个现象可能是出于外界因素(如反复元件的存在)的影响。

小结

在世界各地科学家们每天都要执行序列比对和数据库搜索成千上万此,并且所有的分子生物学都应该熟悉这些要紧的工具。这些方法注定要不断发展,并且接受不断增长的数据库容量的挑战。特别是当可利用的信息增长时,使用者更加难以解释其结果。数据库搜索工作台致力于事后处理搜索结果并且图形显示,从而解决这一问题。这些策略的例子包括PowerBLAST(Zhang and Madden, 1997),BLIXEM(Sonnhammer and Durban, 1994)和BEAUTY(Worley et al., 1995)。

这一章描述了数据比较的一些基本概念,但是使用大量不同的程序以获得更详尽的信息非常有用。研究人员应该了解程序工作的基本操作,并且选择相应的参数。此外,他们应该了解潜在的外部影响并且知道如何避免。最重要的是,应该结合实验方法的发现和评估事物的强大威力。

本次课中所涉及到的可以在互联网上使用(获得)的软件:

CULSTAL.W ftp://ftp.ebi.ac.uk/pub/software/
DOTTER ftp://ftp.sanger.ac.uk/pub/dotter/
LALIGN.FASTA ftp://ftp.virginia.edu/pub/fasta/
BLAST ftp://ncbi.nlm.nih.gov/blast/
SEG ftp://ncbi.nlm.nih.gov/pub/seg/

----------------------------------------------------------------------

多序列比对

双序列比对是序列分析的基础。然而,对于构成基因家族的成组的序列来说,我们要建立多个序列之间的关系,这样才能
揭示整个基因家族的特征。由于可以提高序列比对的信噪比,多序列比对在阐明一组相关序列的重要生物学模式方面起着
相当重要的作用。本章中,我们将介绍一系列多序列比对的方法,从完全手动的方式到广泛应用的计算机程序,即所谓
自动比对的方法。

多序列比对有时用来区分一组序列之间的差异,但其主要用于描述一组序列之间的相似性关系,以便对一个基因家族的特征
有一个简明扼要的了解。与双序列比对一样,多序列比对的方法建立在某个数学或生物学模型之上。因此,正如我们不能对
双序列比对的结果得出“正确或错误”的简单结论一样,多序列比对的结果也没有绝对正确和绝对错误之分,而只能认为所
使用的模型在多大程度上反映了序列之间的相似性关系以及它们的生物学特征。

目前,构建多序列比对模型的方法大体可以分为两大类。第一类是基于氨基酸残基的相似性,如物化性质、残基之间的可突变
性等。另一类方法则主要利用蛋白质分子的二级结构和三级结构信息,也就是说根据序列的高级结构特征确定比对结果。显然,
这两种方法所得结果可能有很大差别。一般说来,很难断定哪种方法所得结果一定正确,应该说,它们从不同角度反映蛋白质
序列中所包含的生物学信息。

基于序列信息和基于结构信息的比对都是非常重要的比对模型,但它们都有不可避免的局限性,因为这两种方法都不能完全反映
蛋白质分子所携带的全部信息。我们知道,蛋白质序列是经过DNA序列转录翻译得到的。从信息论的角度看,它应该与DNA分子所
携带的信息更为“接近”。而蛋白质结构除了序列本身带来的信息外,还包括经过翻译后加工修饰所增加的结构信息,包括残基
的修饰,分子间的相互作用等,最终形成稳定的天然蛋白质结构。因此,这也是对完全基于序列数据比对方法批评的主要原因。
显然,如果能够利用结构数据,对于序列比对无疑有很大帮助。不幸的是,与大量的序列数据相比,实验测得的蛋白质三维结构
数据实在少得可怜。在大多数情况下,并没有结构数据可以利用,我们只能依靠序列的相似性和一些生物化学特性建立一个比较
满意的多序列比对模型。

多序列比对的定义

为了便于描述,我们对多序列比对过程给出下面的定义。把多序列比对看作一张二维表,表中每一行代表一个序列,每一列代表
一个残基的位置。将序列依照下列规则填入表中:(a)一个序列所有残基的相对位置保持不变;(b)将不同序列间相同或相似
的残基放入同一列,即尽可能将序列间相同或相似残基上下对齐(表4.1)。我们称比对前序列中残基的位置为绝对位置。如序列
Ⅰ的第3位的残基是甘氨酸G,则绝对位置Ⅰ3就是甘氨酸,而不能变成任何其它氨基酸。相应地,我们称比对后序列中残基的位置
为相对位置。显然,同一列中所有残基的相对位置相同,而每个残基的绝对位置不同,因为它们来自不同的序列。需要说明的是,
绝对位置是序列本身固有的属性,或者说是比对前的位置,而相对位置则是经过比对后的位置,也就比对过程赋予它的属性。

 

1

2

3

4

5

6

7

8

9

1

I

Y

D

G

G

A

V

-

E

A

L

II

Y

D

G

G

-

-

-

E

A

L

III

F

E

G

G

I

L

V

E

A

L

IV

F

D

-

G

I

L

V

Q

A

V

V

Y

E

G

G

A

V

V

Q

A

L

表4.1 多序列比对的定义

表示五个短序列(I-V)的比对结果。通过插入空位,使5个序列中大多数相同或相似残基放入同一列,并保持每个序列残基顺序不变。

 

算法复杂性

多序列比对的计算量相当可观,因此有必要分析以下技术的复杂性。双序列比对所需要的计算时间和内存空间与这两个序列
的长度有关,或者说正比于这两个序列长度的乘积,用O(m1m2)表示。其中m1、m2是指两条序列的长度。三序列比对则可以
理解为将双序列比对的两维空间扩展到三维,即在原有二维平面上增加一条坐标轴。这样算法复杂性就变成了O(m1m2m3),
其中m3表示第三条序列的长度。

随着序列数量的增加,算法复杂性也不断增加。我们用O(m1m2m3…mn)表示对n个序列进行比对时的算法复杂性,其中mn是最后
一条序列的长度。若序列长度相差不大,则可简化成O(mn),其中n表示序列的数目,m表示序列的长度。显然,随着序列数量
的增加,序列比对的算法复杂性按指数规律增长。

降低算法复杂性,是研究多序列比对的一个重要方面。为此,产生了不少很有实用意义的多序列比对算法。这些方法的特点是利用
启发式(heuristics)算法降低算法复杂性,以获得一个较为满意但并不一定是最优的比对结果,用来找出子序列、构建进化树、
查找保守序列或序列模板,以及进行聚类(clustering)分析等。有的算法将动态规划和启发性算法结合起来。例如,对所有的序列
进行两两比对,将所有的序列与某个特定的序列进行比对,根据某种给定的亲源树进行分组比对,等等。必须指出,上述方法求得的
结果通常不是最优解,至少需要经过n-1次双序列比对,其中n为参与比对的序列个数。

比对方法

下面介绍比对采用的几种常用方法。

手工比对方法

手工比对方法在文献中经常看到。因为难免加入一些主观因素,手工比对通常被认为有很大的随意性。其实,即使用计算机程序进行
自动比对,所得结果中的片面性也不能予以忽视。在运行经过测试并具有比较高的可信度的计算机程序基础上,结合实验结果或文献
资料,对多序列比对结果进行手工修饰,应该说是非常必要的。

多序列比对的软件已经有许多,其中一些带有编辑程序。最好的办法是将自动比对程序和编辑器整合在一起。为了便于进行交互式
手工比对,通常使用不同颜色表示具有不同特性的残基,以帮助判别序列之间的相似性。颜色的选择十分重要,如果使用不当,看
起来不很直观,就会使比对结果中一些有用的信息丢失。相反,如果选择得当,就能从序列比对结果中迅速找到某些重要的结构模式
和功能位点。例如,如果用某种颜色表示一组高度保守的残基,则某个序列的某一位点发生突变时,则由于颜色不同,就可以很快
找出。颜色的选择可以根据主观愿望和喜好,但最好和常规方法一致。用来构筑三维模型的按时氨基酸残基组件和三维分子图形软件
所用的颜色分类方法,比较容易为大家接受(表4.2)。

表4.2 氨基酸分组方法和代表性颜色

残基种类

残基特性

颜色

Asp (D), Glu (E)

酸性

红色

His (H), Arg (R), Lys (K)

碱性

兰色

Ser (S), Thr (T), Asn (N), Gln (Q)

极性

绿色

Ala (A), Val (V), Leu (L), Ile (I), Met (M)

疏水性,带支链

白色

Phe (F), Tyr (Y), Trp (W)

疏水性,带苯环

紫色

Pro (P), Gly (G)

侧链结构特殊

棕色

Cys (C)

能形成二硫键

黄色

* 表中采用的分组方法和用来区分不同组别的颜色与模型构件和三维图形软件中所用方法一致。

多序列比对程序的另一个重要用途是定量估计序列间的关系,并由此推断它们在进化中的亲缘关系。可以通过计算
完全匹配的残基数目或计算完全匹配残基和相似残基的数目得到这种定量关系。这一方法除了可以大略了解序列间
的亲缘关系外,也可用来评估比对质量。如果序列的相似性值低于预料值,那么有可能是序列间亲缘关系较远,也
可能是比对中有错误之处。

现有的软件包用的基本上是可用鼠标点击的窗口界面,其中序列编辑器位于窗口中央。这样的软件包将在以后章节中详细介绍。

使用这些软件我们将会看到,那些长度相仿且相似性程度较高的序列,采用自动比对方法将会得到相当满意的结果;而当
序列长度相差较大而相似性程度较低时,采用自动方法得出的结果则不很理想。此时,手工序列编辑器就接显得十分有用。
通过手工调整,可使结果变得接近实际。此外,采用多种不同的方法进行分析,再将结果综合,是一种行之有效的方法。
为更好地理解多序列比对的原理和规则,应该尽可能学会手工比对的方法,并把比对结果与计算机自动比对得到的结果加
以比较。

同步法

同步法实质是把给定的所有序列同时进行比对,而不是两两比对或分组进行比对。其基本思想是将一个二维的动态规划矩阵扩展
到三维或多维。矩阵的维数反映了参与比对的序列数。这类方法对于计算机的系统资源要求较高,通常是进行少量的较短的序列
的比对。

步进法

这类方法中最常用的就是Clustal,它是由Feng和Doolittle于1987年提出的(Feng和Doolittle,1987)。由于对于实际的数据
利用多维的动态规划矩阵来进行序列的比对不太现实,因此大多数实用的多序列比对程序采用启发式算法,以降低运算复杂度。
Clustal的基本思想是基于相似序列通常具有进化相关性这一假设。比对过程中,先对所有的序列进行两两比对并计算它们的相
似性分数值,然后根据相似性分数值将它们分成若干组,并在每组之间进行比对,计算相似性分数值。根据相似性分数值继续
分组比对,直到得到最终比对结果。比对过程中,相似性程度较高的序列先进行比对,而距离较远的序列添加在后面。作为程序
的一部分,Clusal可以输出用于构建进化树的数据。

Clustal程序有许多版本,ClustalW(Thompson等,1994),根据对亲缘关系较近的序列间空位情况,确定如何在亲缘关系较远的
序列之间插入空位。同样,相似性较高的序列比对结果中的残基突变信息,可用于改变某个特殊位置空位罚分值的大小,推测该
位点的序列变异性。

Clustal是免费软件,很容易从互联网上下载,和其它软件一起,广泛用于序列分析。Clustal所支持的数据格式包括EMBL/SWISSPROT、
NBRF/PIR、Pearson/FastA、GCG/MSF,以及Clustal本身定义的格式。它的输出格式可以是Clustal格式,也可以是可用于GDE、
Phylip、GCG等软件的格式。

多序列比对的数据库

多序列比对的意义在于它能够把不同种属的相关序列的比对结果按照特定的格式输出,并且在一定程度上反映它们之间的相似性。
多序列比对结果所提供的信息对于提高数据库搜索灵敏度也具有很大帮助。因此,方便实用的多序列比对数据库也就应运而生。

目前,互联网上可用的多序列比对数据库已经不少。其中一些利用计算机程序将一次数据库按家族分类;另外一些则是通过手工或
自动方法根据基因家族构建二次数据库。现在我们可以通过一些例子看看这些数据库序列比对的情况,比如说,Pfam是将一次库通过
自动比对来构建的数据库,它将大量具有结构相似性的序列归为一类,比如各种不同种类动物的转铁蛋白的基因序列具有一定的相似性,
Pfam将这些序列归为一类命名为TRANSFERRIN,我们可以在Pfam查找TRANSFERRIN来得到原始序列比对信息[下图],开头是一些注释信息,
然后给出了比对序列的名字,再下是比对结果,以“//”开始,并以“//”结束。对于一个未知的蛋白质序列在该序列库中查询,该
序列库会给出匹配的类及得分供你参考。

transferrin   MSF: 0  Type: P  TRFE_HORSE-23-349  Check: 00 ..

Name: TRFL_HUMAN/26-353  Len:    350  Check:  3574  Weight:  1.00
Name: TRFL_MOUSE/24-351  Len:    350  Check:  2044  Weight:  1.00
Name: TRFE_HORSE/23-349  Len:    350  Check:  4478  Weight:  1.00
Name: TRFE_HUMAN/25-347  Len:    350  Check:  7218  Weight:  1.00
Name: TRFE_PIG/6-332     Len:    350  Check:  5379  Weight:  1.00
Name: TRFL_BOVIN/25-352  Len:    350  Check:  5978  Weight:  1.00
Name: TRFE_CHICK/26-352  Len:    350  Check:  1647  Weight:  1.00
Name: TRF1_SALSA/25-329  Len:    350  Check:  4609  Weight:  1.00
Name: TRFE_XENLA/26-341  Len:    350  Check:  7212  Weight:  1.00
Name: TRFM_HUMAN/23-357  Len:    350  Check:  3624  Weight:  1.00

//
   TRFL_HUMAN/26-353  VQWCAVSQPE ATKCFQWQRN MRKVR...GP PVSCIKRDSP IQCIQAIAEN
   TRFL_MOUSE/24-351  VQWCAVSNSE EEKCLRWQNE MRKVG...GP PLSCVKKSST RQCIQAIVTN
   TRFE_HORSE/23-349  VRWCTVSNHE VSKCASFRDS MKSIVPA.PP LVACVKRTSY LECIKAIADN
   TRFE_HUMAN/25-347  VRWCAVSEHE ATKCQSFRDH MKSVIPSDGP SVACVKKASY LDCIRAIAAN
      TRFE_PIG/6-332  VRWCTISNQE ANKCSSFREN MSKAVKN.GP LVSCVKKSSY LDCIKAIRDK
   TRFL_BOVIN/25-352  VRWCTISQPE WFKCRRWQWR MKKLG...AP SITCVRRAFA LECIRAIAEK
   TRFE_CHICK/26-352  IRWCTISSPE EKKCNNLRDL TQQER....I SLTCVQKATY LDCIKAIANN
   TRF1_SALSA/25-329  VKWCVKSEQE LRKCHDLAAK VA........ EFSCVRKDGS FECIQAIKGG
   TRFE_XENLA/26-341  VRWCVKSNSE LKKCKDLVDT CKNKE....I KLSCVEKSNT DECSTAIQED
   TRFM_HUMAN/23-357  VRWCATSDPE QHKCGNMSEA FREAG..IQP SLLCVRGTSA DHCVQLIAAQ


   TRFL_HUMAN/26-353  RADAVTLDGG FIYEAGLAPY KLRPVAAEVY GTERQPRTHY YAVAVVKKGG
   TRFL_MOUSE/24-351  RADAMTLDGG TLFDAGKPPY KLRPVAAEVY GTKEQPRTHY YAVAVVKNSS
   TRFE_HORSE/23-349  EADAVTLDAG LVFEAGLSPY NLKPVVAEFY GSKTEPQTHY YAVAVVKKNS
   TRFE_HUMAN/25-347  EADAVTLDAG LVYDAYLAPN NLKPVVAEFY GSKEDPQTFY YAVAVVKKDS
      TRFE_PIG/6-332  EADAVTLDAG LVFEAGLAPY NLKPVVAEFY GQKDNPQTHY YAVAVVKKGS
   TRFL_BOVIN/25-352  KADAVTLDGG MVFEAGRDPY KLRPVAAEIY GTKESPQTHY YAVAVVKKGS
   TRFE_CHICK/26-352  EADAISLDGG QAFEAGLAPY KLKPIAAEVY EHTEGSTTSY YAVAVVKKGT
   TRF1_SALSA/25-329  EADAITLDGG DIYTAGLTNY GLQPIIAEDY G..EDSDTCY YAVAVAKKGT
   TRFE_XENLA/26-341  HADAICVDGG DVYKGSLQPY NLKPIMAENY GSHTETDTCY YAVAVVKKSS
   TRFM_HUMAN/23-357  EADAITLDGG AIYEAG.KEH GLKPVVGEVY D..QEVGTSY YAVAVVRRSS


   TRFL_HUMAN/26-353  SFQLNELQGL KSCHTGLRRT AGWNVPIG.. TLRPFLNWTG .PPEPIEAAV
   TRFL_MOUSE/24-351  NFHLNQLQGL RSCHTGIGRS AGWKIPIG.. TLRPYLNWN. GPPASLEEAV
   TRFE_HORSE/23-349  NFQLNQLQGK KSCHTGLGRS AGWNIPIG.. LLYWQLPE.. .PRESLQKAV
   TRFE_HUMAN/25-347  GFQMNQLRGK KSCHTGLGRS AGWNIPIG.. LLYCDLPE.. .PRKPLEKAV
      TRFE_PIG/6-332  NFQWNQLQGK RSCHTGLGRS AGWIIPMG.. LLYDQLPE.. .PRKPIEKAV
   TRFL_BOVIN/25-352  NFQLDQLQGR KSCHTGLGRS AGWIIPMG.. ILRPYLSWT. ESLEPLQGAV
   TRFE_CHICK/26-352  EFTVNDLQGK TSCHTGLGRS AGWNIPIGTL LHRGAIEWEG IESGSVEQAV
   TRF1_SALSA/25-329  AFGFKTLRGK KSCHTGLGKS AGWNIPIGTL VTESQIRWAG IEDRPVESAV
   TRFE_XENLA/26-341  KFTFDELKDK KSCHTGIGKT AGWNIIIGLL LERKLLKWAG PDSETWRNAV
   TRFM_HUMAN/23-357  HVTIDTLKGV KSCHTGINRT VGWNVPVGYL VESGRLSVMG ...CDVLKAV


   TRFL_HUMAN/26-353  ARFFSASCVP GADK.GQFPN LCRLCAGTGE NK..CAFSSQ EPYFSYSGAF
   TRFL_MOUSE/24-351  SKFFSKSCVP GAQK.DRFPN LCSSCAGTGA NK..CASSPE EPYSGYAGAL
   TRFE_HORSE/23-349  SNFFAGSCVP CADR.TAVPN LCQLCVGKGT DK..CACSNH EPYFGYSGAF
   TRFE_HUMAN/25-347  ANFFSGSCAP CADG.TDFPQ LCQLCPG... ....CGCSTL NQYFGYSGAF
      TRFE_PIG/6-332  ASFFSSSCVP CADP.VNFPK LCQQCAGKGA EK..CACSNH EPYFGYAGAF
   TRFL_BOVIN/25-352  AKFFSASCVP CIDR.QAYPN LCQLCKGEGE NQ..CACSSR EPYFGYSGAF
   TRFE_CHICK/26-352  AKFFSASCVP GAT...IEQK LCRQCKGD.. PK..TKCARN APYSGYSGAF
   TRF1_SALSA/25-329  SDFFNASCAP GATM.G..SK LCQLCKG... .D..CSRSHK EPYYDYAGAF
   TRFE_XENLA/26-341  SKFFKASCVP GAKE....PK LSQLCAGIKE HK..CSRSNN EPYYNYAGAF
   TRFM_HUMAN/23-357  SDYFGGSCVP GAGETSYSES LCRLCRGDSS GEGVCDKSPL ERYYDYSGAF


   TRFL_HUMAN/26-353  KCLRDGAGDV AFIRESTVFE DLSDEAER.. ......DEYE LLCPDNTRKP
   TRFL_MOUSE/24-351  RCLRDNAGDV AFTRGSTVFE ELPNKAER.. ......DQYK LLCPDNTWKP
   TRFE_HORSE/23-349  KCLADGAGDV AFVKHSTVLE NLPQEADR.. ......DEYQ LLCRDNTRKS
   TRFE_HUMAN/25-347  KCLKDGAGDV AFVKHSTIFE NLANKADR.. ......DQYE LLCLDNTRKP
      TRFE_PIG/6-332  NCLKEDAGDV AFVKHSTVLE NLPDKADR.. ......DQYE LLCRDNTRRP
   TRFL_BOVIN/25-352  KCLQDGAGDV AFVKETTVFE NLPEKADR.. ......DQYE LLCLNNSRAP
   TRFE_CHICK/26-352  HCLKDGKGDV AFVKHTTVNE NAPDQ..K.. ......DEYE LLCLDGSRQP
   TRF1_SALSA/25-329  QCLKDGAGDV AFIKPLAVPA AEKAS..... ........YE LLCKDGTRAS
   TRFE_XENLA/26-341  KCLQDDQGDV AFVKQSTVPE EFHKD..... ........YE LLCPDNTRKS
   TRFM_HUMAN/23-357  RCLAEGAGDV AFVKHSTVLE NTDGKTLPSW GQALLSQDFE LLCRDGSRAD


   TRFL_HUMAN/26-353  VDKFKDCHLA RVPSHAVVAR S.VNGKEDAI WNLLRQAQEK FGKDKSPKFQ
   TRFL_MOUSE/24-351  VTEYKECHLA QVPSHAVVSR S.TNDKEEAI WELLRQSQEK FGKKQASGFQ
   TRFE_HORSE/23-349  VDEYKDCYLA SIPSHAVVAR S.VDGKEDLI WGLLNQAQEH FGTEKSKDFH
   TRFE_HUMAN/25-347  VDEYKDCHLA QVPSHTVVAR S.MGGKEDLI WELLNQAQEH FGKDKSKEFQ
      TRFE_PIG/6-332  VDDYENCYLA QVPSHAVVAR S.VDGQEDSI WELLNQAQEH FGRDKSPDFQ
   TRFL_BOVIN/25-352  VDAFKECHLA QVPSHAVVAR S.VDGKEDLI WKLLSKAQEK FGKNKSRSFQ
   TRFE_CHICK/26-352  VDNYKTCNWA RVAAHAVVAR D.DN.KVEDI WSFLSKAQSD FGVDTKSDFH
   TRF1_SALSA/25-329  IDSYKTCHLA RVPAHAVVSR ..KDPE..LA NRIYNKLVA. .....VKDFN
   TRFE_XENLA/26-341  IKEYKNCNLA KVPAHAVLTR G.RDDKSKDI IEFLQEAQK. .....TQECK
   TRFM_HUMAN/23-357  VTEWRQCHLA RVPAHAVVVR ADTDGG..LI FRLLNEGQRL FSHE.GSSFQ


   TRFL_HUMAN/26-353  LFGSPSGQ.. ..KDLLFKDS AIGFSRVPPR IDSGLYLGSG YFTAIQNLRK
   TRFL_MOUSE/24-351  LFASPSGQ.. ..KDLLFKES AIGFVRVPQK VDVGLYLTFS YTTSIQNLNK
   TRFE_HORSE/23-349  LFSSPHG... ..KDLLFKDS ALGFLRIPPA MDTWLYLGYE YVTAIRNLRE
   TRFE_HUMAN/25-347  LFSSPHG... ..KDLLFKDS AHGFLKVPPR MDAKMYLGYE YVTAIRNLRE
      TRFE_PIG/6-332  LFSSSHG... ..KDLLFKDS ANGFLKIPSK MDSSLYLGYQ YVTALRNLRE
   TRFL_BOVIN/25-352  LFGSPPGQ.. ..RDLLFKDS ALGFLRIPSK VDSALYLGSR YLTTLKNLRE
   TRFE_CHICK/26-352  LFGPPGKKDP VLKDLLFKDS AIMLKRVPSL MDSQLYLGFE YYSAIQSMRK
   TRF1_SALSA/25-329  LFSSDG...Y AAKNLMFKDS AQKLVQLPTT TDSFLYLGAE YMSTIRSLKK
   TRFE_XENLA/26-341  LFSSP....V WGRDLIFKDS AVSLIPLPSS MDSFLFLGAD YFQCIQALKE
   TRFM_HUMAN/23-357  MFSSEAYG.. .QKDLLFKDS TSELVPIATQ T.YEAWLGHE YLHAMKGLLC

我们也可以看看PRINTS数据库关于TRANSFERRIN的比对信息, PRINTS数据库在自动比对的基础上进行了手工
编辑,查寻PRINTS数据库中关于TRANSFERRIN这一类的比对信息,结果可以用模体(motif)形式显示也可以
用点击链接调用JAVA APPLET进行图形显示,下图是关于TRANSFERRIN序列比对的局部图形,可见PRINTS
数据库中TRANSFERRIN一类由更多的序列比对形成。

一般来说,对于具有较高相似性的一组序列之间的比对,自动比对方法是很有效的。一旦序列的亲缘关系变
得较远,所得结果就不那么可信。若要得到比较可靠而又具有明确生物学意义的比对结果,比较有效的方法
是对比对结果进行手工编辑和调整。这对于构建二次数据库是非常重要的信息。在选择现有的序列模式或序
列模体公开数据库构建自己的数据库系统时,对这些现有数据库的可靠性必须采取谨慎的态度。

 =====================================

 

序列比对和数据库搜索

3 序列比对和数据库搜索

比较是科学研究中最常见的方法,通过将研究对象相互比较来寻找对象可能具备的特性。在生物信息学研究中,比对是最常用和最经典的研究手段。

最常见的比对是蛋白质序列之间或核酸序列之间的两两比对,通过比较两个序列之间的相似区域和保守性位点,寻找二者可能的分子进化关系。进一步的比对是将多个蛋白质或核酸同时进行比较,寻找这些有进化关系的序列之间共同的保守区域、位点和profile,从而探索导致它们产生共同功能的序列模式。此外,还可以把蛋白质序列与核酸序列相比来探索核酸序列可能的表达框架;把蛋白质序列与具有三维结构信息的蛋白质相比,从而获得蛋白质折叠类型的信息。

比对还是数据库搜索算法的基础,将查询序列与整个数据库]的所有序列进行比对,从数据库中获得与其最相似序列的已有的数据,能最快速的获得有关查询序列的大量有价值的参考信息,对于进一步分析其结构和功能都会有很大的帮助。近年来随着生物信息学数据大量积累和生物学知识的整理,通过比对方法可以有效地分析和预测一些新发现基因的功能。

3.1 序列两两比对

序列比对的理论基础是进化学说,如果两个序列之间具有足够的相似性,就推测二者可能有共同的进化祖先,经过序列内残基的替换、残基或序列片段的缺失、以及序列重组等遗传变异过程分别演化而来。序列相似和序列同源是不同的概念,序列之间的相似程度是可以量化的参数,而序列是否同源需要有进化事实的验证。在残基-残基比对中,可以明显看到序列中某些氨基酸残基比其它位置上的残基更保守,这些信息揭示了这些保守位点上的残基对蛋白质的结构和功能是至关重要的,例如它们可能是酶的活性位点残基,形成二硫键的半胱氨酸残基,与配体结合部位的残基,与金属离子结合的残基,形成特定结构motif的残基等等。但并不是所有保守的残基都一定是结构功能重要的,可能它们只是由于历史的原因被保留下来,而不是由于进化压力而保留下来。因此,如果两个序列有显著的保守性,要确定二者具有共同的进化历史,进而认为二者有近似的结构和功能还需要更多实验和信息的支持。通过大量实验和序列比对的分析,一般认为蛋白质的结构和功能比序列具有更大的保守性,因此粗略的说,如果序列之间的相似性超过30%,它们就很可能是同源的。

早期的序列比对是全局的序列比较,但由于蛋白质具有的模块性质,可能由于外显子的交换而产生新蛋白质,因此局部比对会更加合理。通常用打分矩阵描述序列两两比对,两条序列分别作为矩阵的两维,矩阵点是两维上对应两个残基的相似性分数,分数越高则说明两个残基越相似。因此,序列比对问题变成在矩阵里寻找最佳比对路径,目前最有效的方法是Needleman-Wunsch动态规划算法,在此基础上又改良产生了Smith-Waterman算法和SIM算法。在FASTA程序包中可以找到用动态规划算法进行序列比对的工具LALIGN,它能给出多个不相互交叉的最佳比对结果。

在进行序列两两比对时,有两方面问题直接影响相似性分值:取代矩阵和空位罚分。粗糙的比对方法仅仅用相同/不同来描述两个残基的关系,显然这种方法无法描述残基取代对结构和功能的不同影响效果,缬氨酸对异亮氨酸的取代与谷氨酸对异亮氨酸的取代应该给予不同的打分。因此如果用一个取代矩阵来描述氨基酸残基两两取代的分值会大大提高比对的敏感性和生物学意义。虽然针对不同的研究目标和对象应该构建适宜的取代矩阵,但国际上常用的取代矩阵有PAM和BLOSUM等,它们来源于不同的构建方法和不同的参数选择,包括PAM250、BLOSUM62、BLOSUM90、BLOSUM30等。对于不同的对象可以采用不同的取代矩阵以获得更多信息,例如对同源性较高的序列可以采用BLOSUM90矩阵,而对同源性较低的序列可采用BLOSUM30矩阵。

空位罚分是为了补偿插入和缺失对序列相似性的影响,由于没有什么合适的理论模型能很好地描述空位问题,因此空位罚分缺乏理论依据而更多的带有主观特色。一般的处理方法是用两个罚分值,一个对插入的第一个空位罚分,如10-15;另一个对空位的延伸罚分,如1-2。对于具体的比对问题,采用不同的罚分方法会取得不同的效果。

对于比对计算产生的分值,到底多大才能说明两个序列是同源的,对此有统计学方法加以说明,主要的思想是把具有相同长度的随机序列进行比对,把分值与最初的比对分值相比,看看比对结果是否具有显著性。相关的参数E代表随机比对分值不低于实际比对分值的概率。对于严格的比对,必须E值低于一定阈值才能说明比对的结果具有足够的统计学显著性,这样就排除了由于偶然的因素产生高比对得分的可能。

Genbank、SWISS-PROT等序列数据库提供的序列搜索服务都是以序列两两比对为基础的。不同之处在于为了提高搜索的速度和效率,通常的序列搜索算法都进行了一定程度的优化,如最常见的FASTA工具和BLAST工具。FASTA是第一个被广泛应用的序列比对和搜索工具包,包含若干个独立的程序。FASTA为了提供序列搜索的速度,会先建立序列片段的“字典”,查询序列先会在字典里搜索可能的匹配序列,字典中的序列长度由ktup参数控制,缺省的ktup=2。FASTA的结果报告中会给出每个搜索到的序列与查询序列的最佳比对结果,以及这个比对的统计学显著性评估E值。FASTA工具包可以在大多提供下载服务的生物信息学站点上找到。

BLAST是现在应用最广泛的序列相似性搜索工具,相比FASTA有更多改进,速度更快,并建立在严格的统计学基础之上。NCBI提供了基于Web的BLAST服务,用户可以把序列填入网页上的表单里,选择相应的参数后提交到数据服务器上进行搜索,从电子邮件中获得序列搜索的结果。BLAST包含五个程序和若干个相应的数据库,分别针对不同的查询序列和要搜索的数据库类型。其中翻译的核酸库指搜索比对时会把核酸数据按密码子按所有可能的阅读框架转换成蛋白质序列。

表1. BLAST程序:

程序

数据库

blastp

blastn

blastx

tblastn

tblastx

蛋白质

核酸

蛋白质

核苷酸(翻译)

核酸(翻译)

蛋白质

核苷酸

核酸(翻译)

蛋白质

核酸(翻译)

可能找到具有远源进化关系的匹配序列

适合寻找分值较高的匹配,不适合远源关系

适合新DNA序列和EST序列的分析

适合寻找数据库中尚未标注的编码区

适合分析EST序列

 

表2. BLAST的蛋白质数据库:

数据库

简 述

nr

month

swissprot

pdb

yeast

E.coli

Kabat

alu

汇集了SWISS-PROT,PIR,PRF以及从GenBank序列编码区中得到的

蛋白质和PDB中拥有原子坐标的蛋白质,并去除了冗余的序列

nr中过去30天内的最新序列

SWISS-PROT数据库

PDB结构数据库中的蛋白质序列

酵母基因组中编码的全部蛋白质

大肠杆菌基因组中编码的全部蛋白质

Kabat的免疫学相关蛋白质序列

由REPBASE中的Alu重复序列翻译而来,用来遮蔽查询序列中的

重复片段

  

表3. BLAST的核酸数据库:

数据库

nr

month

dbest

dbsts

htgs

yeast

E.coli

pdb

kabat

vector

mito

alu

gss

非冗余的GenBank+EMBL+DDBJ+PDB序列,除了EST、STS、

GSS和0,1,2阶段的HTGS序列

nr中过去30天的最新序列

非冗余的Genbank+EMBL+DDBJ+PDB的EST部分

非冗余的Genbank+EMBL+DDBJ+PDB的STS部分

0,1,2阶段的高产量基因组序列(3阶段完成的HTG序列在nr库里)

酵母的全基因组序列

大肠杆菌的全基因组序列

由三维结构库来的核酸序列

Kabat的免疫学相关序列库

Genbank的载体子集

线粒体核酸序列

REPBASE中Alu重复序列翻译而来,用来遮蔽查询序列中的重复片段

基因组勘测序列(Genome Survey Sequence)

BLAST对序列格式的要求是常见的FASTA格式。FASTA格式第一行是描述行,第一个字符必须是“>”字符;随后的行是序列本身,一般每行序列不要超过80个字符,回车符不会影响程序对序列连续性的看法。序列由标准的IUB/IUPAC氨基酸和核酸代码代表;小写字符会全部转换成大写;单个“-”号代表不明长度的空位;在氨基酸序列里允许出现“U”和“*”号;任何数字都应该被去掉或换成字母(如,不明核酸用“N”,不明氨基酸用“X”)。此外,对于核酸序列,除了A、C、G、T、U分别代表各种核酸之外,R代表G或A(嘌呤);Y代表T或C(嘧啶);K代表G或T(带酮基);M代表A或C(带氨基);S代表G或C(强);W代表A或T(弱);B代表G、T或C;D代表G、A或T;H代表A、C或T;V代表G、C或A;N代表A、G、C、T中任意一种。对于氨基酸序列,除了20种常见氨基酸的标准单字符标识之外,B代表Asp或Asn;U代表硒代半胱氨酸;Z代表Glu或Gln;X代表任意氨基酸;“*”代表翻译结束标志。

BLAST的当前版本是2.0,它的新发展是位点特异性反复BLAST(PSI-BLAST)。PSI-BLAST的特色是每次用profile搜索数据库后再利用搜索的结果重新构建profile,然后用新的profile再次搜索数据库,如此反复直至没有新的结果产生为止。PSI-BLAST先用带空位的BLAST搜索数据库,将获得的序列通过多序列比对来构建第一个profile。PSI-BLAST自然地拓展了BLAST方法,能寻找蛋白质序列中的隐含模式,有研究表明这种方法可以有效的找到很多序列差异较大而结构功能相似的相关蛋白,甚至可以与一些结构比对方法,如threading相媲美。PSI-BLAST服务可以在NCBI的BLAST主页上找到,还可以从NCBI的FTP服务器上下载PSI-BLAST的独立程序。

NCBI的BLUST网址是:http://www.ncbi.nlm.nih.gov/BLAST/

下载BLUST的网址是:ftp://ncbi.nlm.nih.gov/blast/

下载FASTA的网址是:ftp://ftp.virginia.edu/pub/fasta/

 

3.2 多序列比对

顾名思义,多序列比对就是把两条以上可能有系统进化关系的序列进行比对的方法。目前对多序列比对的研究还在不断前进中,现有的大多数算法都基于渐进的比对的思想,在序列两两比对的基础上逐步优化多序列比对的结果。进行多序列比对后可以对比对结果进行进一步处理,例如构建序列模式的profile,将序列聚类构建分子进化树等等。

目前使用最广泛的多序列比对程序是CLUSTALW(它的PC版本是CLUSTALX)。CLUSTALW是一种渐进的比对方法,先将多个序列两两比对构建距离矩阵,反应序列之间两两关系;然后根据距离矩阵计算产生系统进化指导树,对关系密切的序列进行加权;然后从最紧密的两条序列开始,逐步引入临近的序列并不断重新构建比对,直到所有序列都被加入为止。

CLUSTALW的程序可以自由使用,在NCBI的FTP服务器上可以找到下载的软件包。CLUSTALW程序用选项单逐步指导用户进行操作,用户可根据需要选择打分矩阵、设置空位罚分等。EBI的主页还提供了基于Web的CLUSTALW服务,用户可以把序列和各种要求通过表单提交到服务器上,服务器把计算的结果用Email返回用户。

CLUSTALW对输入序列的格式比较灵活,可以是前面介绍过的FASTA格式,还可以是PIR、SWISS-PROT、GDE、Clustal、GCG/MSF、RSF等格式。输出格式也可以选择,有ALN、GCG、PHYLIP和GDE等,用户可以根据自己的需要选择合适的输出格式。

用CLUSTALW得到的多序列比对结果中,所有序列排列在一起,并以特定的符号代表各个位点上残基的保守性,“*”号表示保守性极高的残基位点;“.”号代表保守性略低的残基位点。

EBI的CLUSTALW网址是:http://www.ebi.ac.uk/clustalw/

下载CLUSTALW的网址是:ftp://ftp.ebi.ac.uk/pub/software/

 

0

阅读 评论 收藏 转载 喜欢 打印举报/Report
  • 评论加载中,请稍候...
发评论

    发评论

    以上网友发言只代表其个人观点,不代表新浪网的观点或立场。

      

    新浪BLOG意见反馈留言板 电话:4000520066 提示音后按1键(按当地市话标准计费) 欢迎批评指正

    新浪简介 | About Sina | 广告服务 | 联系我们 | 招聘信息 | 网站律师 | SINA English | 会员注册 | 产品答疑

    新浪公司 版权所有