序列比对和数据库搜索

Gregory D.Schuler

National Center for Biotechnology Information

National Library of Medicine. National Institutes of Health

Bethesda. Maryland

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

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

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

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

列中缺少这种残基，表明这些残基是插入到前者或是从后者中丢失的。这些空位在序列比对时用连续的短线填补。如图7.1，在序列比对中，发现了5个空位。

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

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

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

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

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

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

图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描述了ＰＬＡＴ和ＰＬＡＵ中与ＥＧＦ相似的结构域之间进行比较时的比对、点阵和路径图三种方法的关系。

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

图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算法只会报告出最好的一个比对，改良的算法会报告出第二和第三的比对方式，从而显示出功能结构域。

F12 VQPVCLPSGAARPSETTLCQ—VAGWGHQFEGAEEYASFLQEAQVPFLSLERCSAPDVHG

图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

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|113999|sp|P16550|APA1_YEAST 5’,5’’’-P-1,P-4-TE 108 66 103 137.1 0.47

FHIT VPGHVLVCPLRPVERFHDLRPDEVADLFQTTQRVGTVVEKHFHGTSLTFSM—QDGP---

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

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

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

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

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