序列分析

一、碱基组成 DNA序列一个显而易见的特征是四种碱基类型的分布。尽管四种碱基的频率相等时对数学模型的建立可能是方便的，但几乎所有的研究都证明碱基是以不同频率分布的。

表1包含了9个完整DNA分子序列的资料，表2的数据来自两个胎儿球蛋白基因(Gr和Ar)，每个基因具有三个外显子和两个内含子(shen等1981)。这两个例子说明序列内和序列间碱基具有不同的频率。在基因每一侧的500 个任意碱基区域被称为“侧翼”，基因间区域是指两个基因间的其余序列。

表1 九种完整DNA序列的碱基组成

表2 人类胎儿球蛋白基因不同区段的碱基组成

二．碱基相邻频率 分析DNA序列的主要困难之一是碱基相邻的频率不是独立的。碱基相邻的频率一般不等于单个碱基频率的乘积 例： 鸡血红蛋白β链的mRNA编码区的438个碱基

图1 鸡β球蛋白基因编码区的DNA序列 (GenBank：CHKHBBM，记录号J00860)

表3 图1鸡β球蛋白基因序列的相邻碱基分布

　在编码区，存在某种约束来限制DNA序列编码氨基酸。在密码子水平上，这一约束与碱基相邻频率有关。 表4列出了遗传密码和图1序列中各密码子数量。尽管数目很小，难以作出有力的统计结论，但编码同一氨基酸的不同密码子(同义密码子)好像不是等同存在的。这种密码子偏倚必定与两碱基相邻频率水平有关。 表4还清楚地表明，由于密码子第3位置上碱基的改变常常不会改变氨基酸的类型，因而对第3位置上碱基的约束要比第 2位碱基小得多。

表4 64种可能的碱基三联体密码子及相应的氨基酸数（据图1序列）

　相邻碱基之间的关联将导致更远碱基之间的关联，这些关联延伸距离的估计可以从马尔科夫链(Markov chain)理论得到(Javare和Giddings，1989)

三．同向重复序列分析 　　除了分析整个序列碱基关联程度的特征外，我们常对寻找同向重复序列(direct repeats)之类的问题感兴趣。Karlin等(1983)给出了完成这一分析的有效算法。该法采用由特定的几组碱基字母组成的不同亚序列或称为字码(word)。只需要对整个序列搜索一次。给一碱基赋以值α,例如A、C、G、T的值为0、1、2、3。由X1、X2、…、Xk 共k个字母组成的每一种不同的字码按： 计算字码值。这些值的取值范围为1到4k

在本例中只有4个重复的2碱基重复序列。例如，在位置4、5、8、9、10和15均发现了字码值为1的碱基重复序列。 　例如：5字码TGACC的值为1+3×44+2×43+0×42+1×41+1×40=459。可先从低k值的字码开始搜索。记录序列中每一个位置k字码的字码值。只有在发现k字码长度重复的那些位置考虑进行长度大于k的字码搜索。 　序列TGGAAATAAAACGTAAGTAG中所有碱基2字码(k=2)的初始位置和字码值。对于完全重复、长度大于2的同向重复或亚序列的搜索可只限于2字码重复的初始位置。 　在本例中只有4个重复的2碱基重复序列。例如，在位置4、5、8、9、10和15均发现了字码值为1的碱基重复序列。 　从有重复的2碱基为起点的3字码值中发现字码值为1、45和49的序列有重复；以每一重复的3碱基为起点的4字码搜索未能发现更长的重复序列。

表5 序列TGGAAATAAAACGTAAGTAG的3字码值和位置(Karlin, 1983)

四、RNA二级结构预测 尽管现有一些RNA折叠程序可以预测RNA二级结构，但这类分析仍然是一门艺术。RNA折叠有助于找出RNA分子中可能的稳定茎区，但对给定的RNA分子来说，这一结果的生物学意义究竟有多大，还是一个未知数。即使有此局限性，二级结构的预测还是有助于找出mRNA控制区以及RNA分子中可能形成稳定折叠结构的区段。

五、从序列中寻找基因 1.基因及基因区域预测 基因按其功能可分为结构基因和调控基因：结构基因可被转录形成mRNA，并进而转译成多肽链；调控基因是指某些可调节控制结构基因表达的基因。在DNA链上，由蛋白质合成的起始密码开始，到终止密码子为止的一个连续编码序列称为一个开放阅读框(Open Reading Frame,ORF)。结构基因多含有插入序列，除了细菌和病毒的DNA中ORF是连续的，包括人类在内的真核生物的大部分结构基因为断裂基因，即其编码序列在DNA分子上是不连续的，或被插入序列隔开。断裂基因被转录成前体mRNA，经过剪切过程，切除其中非编码序列(即内含子)，再将编码序列(即外显子)连接形成成熟mRNA，并翻译成蛋白质。假基因是与功能性基因密切相关的DNA序列，但由于缺失、插入和无义突变失去阅读框而不能编码蛋白质产物。

一种典型的真核蛋白质编码基因的结构示意图。其编码序列（外显子）是不连续的，被非编码区（内含子）隔断。

所谓基因区域预测，一般是指预测DNA序列中编码蛋白质的部分，即外显子部分。 不过目前基因区域的预测已从单纯外显子预测发展到整个基因结构的预测。这些预测综合各种外显子预测的算法和人们对基因结构信号(如TATA盒等)的认识，预测出可能的完整基因

基因区域的预测是一个活跃的研究领域，先后有一大批预测算法和相应程序被提出和应用，其中有的方法对编码序列的预测准确率高达90%以上，而且在敏感性和特异性之间取得了很好的平衡 预测方法中，最早是通过序列核苷酸频率、密码子等特性进行预测(如最长ORF法等)，随着各类数据库的建立和完善，通过相似性列线比对也可以预测可能的基因。同时，一批新方法也被提了出来，如隐马尔可夫模型(Hidden Markov Model,HMM)、动态规划法(dynamic programming)、法则系统(ruled-based system)、语言学(linguistic)方法、线性判别分析(Linear Discriminant Analysis,LDA)、决策树(decision tree)、拼接列线(spliced alingment)、博利叶分析(Fourier analysis)等。 下表列出了claverie(1997)对部分程序预测基因区域能力的比较结果，表中同时列出了相应算法和程序的网址。

目前基因区域预测的各种算法均存在以下2个问题 （1）目前算法对基因中的非编码区和基因间序列不加任何区别，所以预测出的基因仍然是不完全的，对5‘和3‘非编译区（UTR，untranslated region）的预测基本上还是空白； （2）目前大多数算法都是基于已知基因序列。如相似性列线比较算法是完全依赖于已知的序列，而象HMM之类的算法都需要对已知的基因结构信号进行学习或训练，由于训练所用的序列毕竟是有限的，所以对那些与学习过的基因结构不太相似的基因，这些算法的预测效果就要大打折扣了 要解决以上两个问题，需要对基因结构进行更深入的研究，寻找隐藏在基因不同结构中的内在统计规律。

2．发现基因的一般过程 从序列中发现基因可以理解为基因区域预测和基因功能预测2个层次 第一步：获取DNA目标序列 ① 如果你已有目标序列，可直接进入第2步； ② 可通过PubMed查找你感兴趣的资料；通过GenBank或EMBL等数据库查找目标序列

第二步：查找ORF并将目标序列翻译成蛋白质序列 利用相应工具，如ORF Finder、Gene feature(Baylor College of Medicine)、GenLang(University of Pennsylvania)等，查找ORF并将DNA序列翻译成蛋白质序列 第三步：在数据库中进行序列搜索 可以利用BLAST进行ORF核苷酸序列和ORF翻译的蛋白质序列搜索 第四步：进行目标序列与搜索得到的相似序列的整体列线(global alignment) 虽然第三步已进行局部列线(local alignment)分析，但整体列线有助于进一步加深目标序列的认识

第五步：查找基因家族 第六步：查找目标序列中的特定模序 第七步：预测目标序列结构 进行多序列列线(multiple sequence alignment)和获得列线区段的可视信息。可分别在AMAS(Oxford University)和BOXSHADE(ISREC,Switzerland)等服务器上进行 第六步：查找目标序列中的特定模序 ① 分别在Procite、BLOCK、Motif数据库进行profile、模块(block)、模序(motif)检索； ② 对蛋白质序列进行统计分析和有关预测 第七步：预测目标序列结构 可以利用PredictProtein(EMBL)、NNPREDICT(University of California)等预测目标序列的蛋白质二级结构

第八步：获取相关蛋白质的功能信息 第九步：把目标序列输入“提醒”服务器 为了了解目标序列的功能，收集与目标序列和结构相似蛋白质的功能信息非常必要。可利用PubMed进行搜索 第九步：把目标序列输入“提醒”服务器 如果有与目标序列相似的新序列数据输入数据库，提醒(alert)服务会向你发出通知。可选用Sequence Alerting(EMBL)、Swiss-Shop(Switzerland)等服务器

3．解读序列(making sense of the sequence) 大致有2条途径可以发现基因： (1)基于同源性的方法，包括已知mRNA序列的应用； (2)基因家族和特殊序列间的比较。 最初的方法包括利用各种计算机手段分析外显子和其它序列信号，如酶切位点

六、基于编码区特性：最长ORF法 基因区域或蛋白质编码区的识别，特别是对高等真核生物基因组DNA序列中编码区的识别仍未能实现完全自动化。将每条链按6个读框全部翻译出来，然后找出所有可能的不间断开放阅读框(ORF)往往有助于基因的发现

预测基因组的全部编码区或称为开放阅读框的方法概括来说也可以分为三类: 一类是基于编码区所具有的独特信号，如始起密码子、终止密码子等； 二是基于编码区的碱基组成不同于非编码区，这是由于蛋白质中20种氨基酸出现的概率、每种氨基酸的密码子兼并度和同一种氨基酸的兼并密码子使用频率不同等原因造成的； 三是通过同源性比较搜寻蛋白质库或dbEST库寻找编码区。前二类方法主要是利用编码区的特性来寻找，下面对这二类方法做简单描述

最长ORF法：在细菌基因组中，蛋白质编码基因从起始密码ATG到终止密码平均有100bp，而300bp长度以上的ORF平均每36Kb才出现一次，所以只要找出序列中最长的ORF(>300bp)就能相当准确地预测出基因

利用编码区与非编码区密码子选用频率的差异进行编码区的统计学鉴别方法：由于内含子的进化不受约束，而外显子则受到选择压力，因此内含子的序列要比外显子更随机。这是目前各种预测程序中被广泛应用的一种方法，如GCG(Genetic Computer Group 研制，一种通用核酸、蛋白质分析软件包)的TestCode、美波士顿大学GeneID和Baylor Medcine College的BCM Gene Finder等程序均利用了这一方法

CpG岛：CpG岛(CpG island)一词是用来描述哺乳动物基因组DNA中的一部分序列，其特点是胞嘧啶(C)与鸟嘌呤(G)的总和超过4种碱基总和的50%，即每10个核苷酸约出现一次双核苷酸序列CG。具有这种特点的序列仅占基因组DNA总量的10%左右。从已知的DNA序列统计发现，几乎所有的管家基因(House-Keeping gene)及约占40%的组织特异性基因的5‘末端含有CpG岛，其序列可能包括基因转录的启动子及第一个外显子。因此，在大规模DNA测序计划中，每发现一个CpG岛，则预示可能在此存在基因。另外，AT含量也可以作为编码区的批示指标之一

七、序列比对 概念： 相似性和同源性 局部相似性和整体相似性 相似性分数矩阵 数据库的搜索 FastA BLAST

相似性和同源性 数据库搜索的基础是序列的相似性比对，而寻找同源序列则是数据库搜索的主要目的之一。 所谓同源序列，简单地说，是指从某一共同祖先经趋异进化而形成的不同序列。同源性可以用来描述染色体—“同源染色体”、基因—“同源基因”和基因组的一个片断—“同源片断” 必须指出，相似性(similarity)和同源性(homology)是两个完全不同的概念。

相似性是指序列比对过程中用来描述检测序列和目标序列之间相同DNA碱基或氨基酸残基顺序所占比例的高低。相似性本身的含义，并不要求与进化起源是否同一，与亲缘关系的远近、甚至于结构与功能有什么联系。 当相似程度高于50%时，比较容易推测检测序列和目标序列可能是同源序列；而当相似性程度低于20%时，就难以确定或者根本无法确定其是否具有同源性。 总之，不能把相似性和同源性混为一谈。所谓“具有50%同源性”，或“这些序列高度同源”等说法，都是不确切的，应该避免使用。

而同源又有两种不同的情况即垂直方向的(orthology)与水平方向的(paralogy)。 (1)在进化上起源于一个始祖基因并垂直传递(vertical descent)的同源基因； (2)分布于两种或两种以上物种的基因组； (3)功能高度保守乃至于近乎相同，甚至于其在近缘物种可以相互替换； (4)结构相似； (5)组织特异性与亚细胞分布相似

鉴定直系同源的实际操作标准(practical criteria)为： 如基因组Ⅰ中的A基因与基因组Ⅱ中的A‘基因被认为是直系同源，则要求： (1)A‘的产物比任何在基因组Ⅱ中所发现的其它基因产物都更相似于A产物； (2)A‘与A的相似程度比在任何一个亲缘关系较远的基因组中的任一基因都要高； (3)A编码的蛋白与A‘编码的蛋白要从头到尾都能并排比较，即含有相似以至于相同的模序(motif)

旁系同源(paralogy)基因是指同一基因组(或同系物种的基因组)中，由于始祖基因的加倍而横向(horizontal)产生的几个同源基因。 直系与旁系的共性是同源，都源于各自的始祖基因。其区别在于：在进化起源上，直系同源是强调在不同基因组中的垂直传递，旁系同源则是在同一基因组中的横向加倍；在功能上，直系同源要求功能高度相似，而旁系同源在定义上对功能上没有严格要求，可能相似，但也可能并不相似(尽管结构上具一定程度的相似)，甚至于没有功能(如基因家族中的假基因)。旁系同源的功能变异可能是横向加倍后的重排变异或进化上获得了另一功能，其功能相似也许只是机械式的相关(mechanistically related)，或非直系同源基因取代新产生的非亲缘或远缘蛋白在不同物种具有相似的功能。

序列比对的基本思想，是找出检测序列和目标序列的相似性。比对过程中需要在检测序列或目标序列中引入空位，以表示插入或删除（图2）。 局部相似性和整体相似性 序列比对的基本思想，是找出检测序列和目标序列的相似性。比对过程中需要在检测序列或目标序列中引入空位，以表示插入或删除（图2）。 图2 序列比对，图中“-”表示插入和删除，用字符表示相同的残基，“+”表示相似残基

序列比对的最终实现，必须依赖于某个数学模型。不同的模型，可以从不同角度反映序列的特性，如结构、功能、进化关系等。很难断定，一个模型一定比另一个模型好，也不能说某个比对结果一定正确或一定错误，而只能说它们从某个角度反映了序列的生物学特性。此外，模型参数的不同，也可能导致比对结果的不同。

序列比对的数学模型大体可以分为两类，一类从全长序列出发，考虑序列的整体相似性，即整体比对；第二类考虑序列部分区域的相似性，即局部比对。 局部相似性比对的生物学基础是蛋白质功能位点往往是由较短的序列片段组成的，这些部位的序列具有相当大的保守性，尽管在序列的其它部位可能有插入、删除或突变。此时，局部相似性比对往往比整体比对具有更高的灵敏度，其结果更具生物学意义。

区分这两类相似性和这两种不同的比对方法，对于正确选择比对方法是十分重要的。应该指出，在实际应用中，用整体比对方法企图找出只有局部相似性的两个序列之间的关系，显然是徒劳的；而用局部比对得到的结果也不能说明这两个序列的三维结构或折叠方式一定相同。 BLAST和FastA等常用的数据库搜索程序均采用局部相似性比对的方法，具有较快的运行速度，而基于整体相似性比对的数据库搜索程序则需要超级计算机或专用计算机才能实现。

　有2种经典方法可以计算两条序列间的最适联配。Needleman-Wunsch算法是一种整体联配(global alignment)算法，最佳联配（两条蛋白质序列具有最多匹配残基）中包括了全部的最短匹配序列。 Smith-Wateman算法是在Needleman-Wunsch算法基础上发展而来的，它是一种局部联配(Local alignment)算法。 这二种算法均可以用于核酸和蛋白质序列。在给定空位罚值和替换矩阵情况下，它们总是能给出具有最高联配值的联配。但是，这个联配并不需要达到生物学意义上的显著水平。

　　许多程序可通过匿名ftp服务用于两条序列的联配计算。GCG软件包中，BESFIT和GAP程序便是用于两对序列的联配。在一些网站可以进行两条序列的联配分析，例如：ALIGN(http://genome.eerie.fr/fasta/alignquery.html)/Align(http://www.mips.biochem.mpg.de/mips.de/mips/programs/align.html)。ALIGN允许用户提供序列进行联配，允许选择替换矩阵，但不能设置空位罚值。Align只能进行数据库中已有序列间的联配分析。

Needleman-Wunsch算法 　从整体上分析两个序列的关系，即考虑序列总长的整体比较，用类似于使整体相似 (global similarity)最大化的方式，对序列进行联配。两个不等长度序列的联配分析必需考虑在一个序列中圈掉一些碱基或在另一序列作空位(gap)处理。 Needleman 和Wunsch(1970)的法则为这些步骤提供了实例。这一算法是为氨基酸序列发展的，但也可以用于核苷酸序列。算法最初寻求的是使两条序列间的距离最小。尽管这类距离的元素是以一种特定的方式定义的，但该算法的良好特性在于它确定了最短距离。这是一个动态规划(dynamic programming)的方法。

Needleman-Wunsch算法 　将两条联配的序列沿双向表的轴放置。从任一碱基对，即表中的任一单元开始，联配可延三种可能的方式延伸：如果碱基不匹配，则每一序列加上一个碱基，并给其增加一个规定的距离权重；或在一个序列中增加一个碱基而在另一序列中增加一个空位或反之亦然。引入一个空位时也将增加一个规定的距离权重。

Needleman-Wunsch算法 　 因此，表中的一个单元可以从(至多)三个相邻的单元达到。我们把到左上角单元距离最小的方向看作相似序列延伸的方向。等距离时意味着存在两种可能的方向。将这些方向记录下来，并在研究了所有的单元之后，沿着记录的方向就有一条路径可从右下角(两个序列的末端)追踪到左上角 (两个序列的起点)。由此所产生的路径将给出具有最短距离的序列联配。

以两个短序列CTGTATC和CTATAATCCC为例： Needleman-Wunsch算法 以两个短序列CTGTATC和CTATAATCCC为例： 　设碱基错配时距离权重为1，引入一个空位时距离权重为3。该图边缘的行和列作为起始条件增加到表中。在单元5行3列，即相应较短序列(第二序列)的第2个T碱基和较长序列(第一序列)的第1个T碱基位置，有三种可能的距离增量。设在各序列中增加碱基T时 (从4行2列移动)对距离的贡献为0。从5行2列的位置作水平移动(等价于增加第二序列的碱基T而在第一序列引入一个空位)，在本例中增加一个罚值3。从3列4行向该单元作垂直移动，使第一序列增加碱基T而第二序列引入一个空位，结果也得到一个罚值3。因此从该单元(5行3列)所得到的最小距离的延伸方向是沿对角线和水平方向。在表中这两个方向用箭头表示。这两种最短方向都使从左上角到该单元的距离为6。沿箭头所指方向在表中从右下角向左上角追踪，得到6种可能的联配：

在上述6种联配中，距离均为10，即在较短序列中有6个匹配碱基、1个错配碱基和3个空位

Needleman-Wunsch算法 当两个序列被联配时，通过计算其重排序列(shuffed version)的联配距离，可以得到这两个序列间的最小距离估计。如果实际得到的联配距离小于重排序列距离的95%，则表明实际的联配距离达到了5%的显著水平，是不可能由机误造成的。

由于亲缘关系较远的蛋白质序列可能只有一些相互独立的相同片段，所以进行局部相似性分析有时可能比整体相似性分析更合理。 Smith-Waterman算法 　由于亲缘关系较远的蛋白质序列可能只有一些相互独立的相同片段，所以进行局部相似性分析有时可能比整体相似性分析更合理。 Smith和Waterman描述了一种查找具有最高相似性片段的算法。对于序列A=(a1,a2,…,am)和 B=(b1,b2,…,bn)，Hij被定义为以ai和bj 碱基对结束的片段(亚序列)的相似性值。 与Needle-Wunsch算法一样，Smith-Waterman算法也要利用递推关系来确定

相似性计算中包括2个统计量：碱基对(序列因子) 的相似性值和空位权重 (k 为空位长度)。 Smith-Waterman算法 　相似性计算中包括2个统计量：碱基对(序列因子) 的相似性值和空位权重 (k 为空位长度)。 　Smith-Waterman算法可以给出2条序列的最大相似性值。

相似性分数矩阵 在对蛋白质数据库搜索时，可采用不同的相似性分数矩阵，以提高搜索的灵敏度和准确率。常用的相似性矩阵有突变数据矩阵(Mutation Data Matrix，简称MD)和模块替换矩阵(BLOcks Substitution Matrix，简称BLOSUM)。

在序列比对中，通常希望使用能够反映一个氨基酸发生改变的概率与两个氨基酸随机出现的概率的比值的矩阵。这些比值可以用相关几率（relatedness odds）矩阵表示。这就是突变数据相似性分数矩阵产生的基础，在序列比对过程中，两个序列从头到尾逐个残基进行比对，所得几率值的乘积就是整个比对的分值。 在实际使用时，通常取几率值的对数以简化运算。因此，常用的突变数据矩阵PAM250实际上是几率值的对数矩阵（图3）。矩阵中值大于0的元素所对应的两个残基之间发生突变的可能性较大，值小于0的元素所对应的两个残基之间发生突变的可能性较小

图3 突变数据相似性分数矩阵PAM250

突变数据矩阵PAM即可接受点突变(Point Accepted Mutation，简称 PAM)。1个PAM的进化距离表示100个残基中发生一个残基突变的概率。对应于一个更大进化距离间隔的突变概率矩阵，可以通过对初始矩阵进行适当的数学处理得到[Dayhoff等，1978]，如常用的PAM250矩阵，PAM250相似性分数矩阵相当于在两个序列之间具有20%的残基匹配(图3)。

主对角线上分数值是指两个相同残基之间的相似性分数值，有些残基的分值较高，如色氨酸W为17、半胱氨酸C为12，说明它们比较保守，不易突变；有的残基的分值较低，如丝氨酸S、丙氨酸A、门冬酰氨N三种氨基酸均为2，这些氨基酸则比较容易突变。不同氨基酸之间的分数值越高，它们之间的相似性越高，进化过程中容易发生互相突变，如苯丙氨酸F和酪氨酸Y，它们之间的相似性分数值是7。而相似性分数值为负数的氨基酸之间的相似性则较低，如甘氨酸和色氨酸之间为-7，它们在进化过程中不易发生互相突变。此外，表中把理化性质相似的氨基酸按组排列在一起，如碱性氨基酸组氨酸H、精氨酸R和赖氨酸K。

突变数据矩阵的产生基于相似性较高（通常为85%以上）的序列比对，那些进化距离较远的矩阵（如PAM250）是从初始模型中推算出来而不是直接计算得到的，其准确率受到一定限制。而序列分析的关键是检测进化距离较远的序列之间是否具有同源性，因此突变数据矩阵在实际使用时存在着一定的局限性。

而模块替换矩阵BLOSUM则以序列片段为基础，它是基于蛋白质模块数据库BLOCKS，Henikoff夫妇（Henikoff和Henikoff，1992）从蛋白质模块数据库BLOCKS中找出一组替换矩阵，用于解决序列的远距离相关。在构建矩阵过程中，通过设置最小相同残基数百分比将序列片段整合在一起，以避免由于同一个残基对被重复计数而引入的任何潜在的偏差。在每一片段中，计算出每个残基位置的平均贡献，使得整个片段可以有效地被看作为单一序列。

通过设置不同的百分比，产生了不同矩阵。由此，例如高于或等于80%相同的序列组成的串可用于产生BLOSUM80矩阵（BlOcks SUbstitution Matrix 发音为blossom）；那些有62%或以上相同的串用于产生BLOSUM62矩阵，依此类推。 BLOSUM与BLOCKS对于同样的序列比对产生的结果在局部有所不同，可能是一个认为不相似不可以替换而另一个认为相似可以替换。必须说明，如果比对这两个序列高度相似，这些细微的差别对整个序列比对结果的影响不大，但在序列比对的边界区可能产生显著影响，此时增强微弱信号以探测远距离相关变得十分重要。

数据库的搜索简介

数据库查询为生物学研究提供了一个重要工具，在实际工作中经常使用。然而，在分子生物学研究中，对于新测定的碱基序列或由此翻译得到的氨基酸序列，往往需要通过数据库搜索，找出具有一定相似性的同源序列，以推测该未知序列可能属于哪个基因家族，具有哪些生物学功能。对于氨基酸序列来说，有可能找到已知三维结构的同源蛋白质而推测其可能的空间结构。因此，数据库搜索与数据库查询一样，是生物信息学研究中的一个重要工具。

数据库搜索的基础是序列的相似性比对，即双序列比对(pairwise alignment)。 新测定的、希望通过数据库搜索确定其性质或功能的序列称作检测序列(probe sequence)；通过数据库搜索得到的和检测序列具有一定相似性的序列称目标序列(subject sequence)。 为了确定检测序列和一个已知基因家族之间的进化关系，在通过数据库搜索得到某些相似序列后，还需要判断其序列相似性程度。如果检测序列和目标序列的相似性程度很低，还必须通过其它方法或实验手段才能确定其是否属于同一基因家族

比对统计学意义的评价--E值(E-Value) P值(P-Value)(概率值) BLAST程序中使用了E值而非P值，这主要是从直观和便于理解的角度考虑。比如E值等于5和10，总比P值等于0.993和0.99995更直观。但是当E<0.01时，P值与E值接近相同 参数K和λ可分别被简单地视为搜索步长(search spacesize)和计分系统(scoring system)的特征数

BLAST和FASTA数据库搜索策略 一种思路是把数据库中的所有蛋白序列与待查序列的关系都视为相同重要，也就是说对于E值均较低的短和长序列，它们是等同重要的。FASTA程序近期版本便是采用这一策略 另一种思路是把长序列视为比短序列更重要，因为长序列往往包括更多的特异功能域(domain)。如果对序列长度上进行相关优先处理，则在计算数据库序列长度为n的E值时，将乘以N/n，其中N为数据库中序列的总长度。E值的计算可简单地把整个数据库序列视为长度为N的单条序列。BLAST程序采用了这一策略 FASTA策略中E值的计算还需再乘上数据库的序列条数。如果考虑到核酸数据库的序列长度变化更大，则在DNA序列相似性搜索时，BLAST的策略可能会是合理的选择

一些数据库搜索程序，例如FASTA或其它基于Smith-Waterman算法的程序，在进行序列搜索时，会对数据库中的每条序列进行联配并给出联配值，这些值大部分与未知序列无关，但它们被用于了K和λ参数的估计。这一方法避免了随机序列模型因使用真实序列(real sequence)造成的随意性，但同时产生了使用相关序列估计参数的难题 BLAST仅通过部分而不是全部无关序列计算最适联配值，这赢得了搜索速度。因此，对于某一选定的替换矩阵和空位罚值，必须进行K和λ参数的预先估计，估计中使用真实序列，而非通过随机序列模型产生的模拟序列。这一估计的结果看来非常准确。

表6 数据库相似性搜索程序BLAST和FASTA程序清单 注：n：核酸序列或核酸序列库；p：蛋白质序列或蛋白质序列库

搜索实例

FastA和BLAST程序是目前最常用的基于局部相似性的数据库搜索程序，它们都基于查找完全匹配的短小序列片段，并将它们延伸得到较长的相似性匹配。它们的优势在于可以在普通的计算机系统上运行，而不必依赖计算机硬件系统而解决运行速度问题。

可以访问NCBI的网站在线进行BLAST和FastA的搜索 BLAST是目前常用的数据库搜索程序，它是Basic Local Alignment Search Tool的缩写，意为“基本局部相似性比对搜索工具”[Altschul, 1990, 1997]。国际著名生物信息中心都提供基于Web的BLAST服务器。BLAST程序之所以使用广泛，主要因为其运行速度比FastA等其它数据库搜索程序快，而改进后的BLAST程序允许空位的插入。 可以访问NCBI的网站在线进行BLAST和FastA的搜索

BLAST搜索 BLAST算法本身很简单，它的基本要点是序列片段对（segment pair）的概念。所谓序列片段对是指两个给定序列中的一对子序列，它们的长度相等，且可以形成无空位的完全匹配。 BLAST算法首先找出代查序列和目标序列间所有匹配程度超过一定阈值的序列片段对，然后对具有一定长度的片段对根据给定的相似性阈值延伸，得到一定长度的相似性片段，称高分值片段对（high-scoring pairs, HSPs）。这就是无空位的BLAST比对算法的基础，也是BLAST输出结果的特征。

BLAST软件包实际上是综合在一起的一组程序，不仅可用于直接对蛋白质序列数据库和核酸序列数据库进行搜索，而且可以将检测序列翻译成蛋白质或将数据库翻译成蛋白质后再进行搜索，以提高搜索结果的灵敏度(表7)。

表7 BLAST程序检测序列和数据库类型 程序名 检测序列 数据库类型 方 法 Blastp 蛋白质 用检测序列蛋白质搜索蛋白质序列数据库 Blastn 核酸 用检测序列核酸搜索核酸序列数据库 Blastx 将核酸序列按6条链翻译成蛋白质序列后搜索蛋白质序列数据库 Tblastn 用检测序列蛋白质搜索由核酸序列数据库按6条链翻译成的蛋白质序列数据库 Tblastx 将核酸序列按6条链翻译成蛋白质序列后搜索由核酸序列数据库按6条链翻译成的蛋白质序列数据库

BLAST程序是免费软件，可以从美国国家生物技术信息中心NCBI等文件下载服务器上获得，安装在本地计算机上，包括UNIX系统和WINDOWS系统的各种版本。但必须有BLAST格式的数据库，可以从NCBI下载，也可以利用该系统提供的格式转换工具由其它格式的核酸或蛋白质序列数据库经转换后得到。对核酸序列数据库而言，不论用哪种方式，都需要很大的磁盘空间；而程序运行时，需要有较大的内存和较快的运算速度，因此必须使用高性能的服务器。

对一般用户来说，目前常用的办法是通过NCBI、EBI等国际著名生物信息中心的BLAST服务器进行搜索。北京大学生物信息中心也提供了BLAST数据库搜索服务。需要说明的是，各生物信息中心BLAST用户界面有所不同，所提供的数据库也可能不完全相同，使用前最好先进行适当的选择

欧洲生物信息研究所BLAST服务器的用户界面(图4)比较简洁，提供的数据库和参数很多，用户可以根据不同要求，选择不同的数据库和各种参数。一般情况下，可以先按照系统给定的缺省参数进行初步搜索，对结果进行分析后再适当调整参数，如改变相似性矩阵、增加或减少空位罚分值、调节检测序列滑动窗口大小等。对于核酸序列数据库，一般选择重复序列屏蔽功能，而对于蛋白质序列，特别是球蛋白，通常不必选择重复序列屏蔽功能。

图4 欧洲生物信息学研究所的BLAST服务器的用户界面

图5是BLAST程序运行结果实例。这里，检测序列是与细胞凋亡有关的人自噬基因氨基酸序列，通过欧洲生物信息学研究所的BLAST服务器对包括SwissProt和TrEMBL数据库在内的蛋白质数据库进行搜索。输出结果中包括程序名称、版本号以及文献引用出处，以及检索序列的名称、数据库名称；列出相似性值较高的序列条目，以及它们在数据库中的编号和简要说明。每个条目后面给出相似性分数值Score和期望频率值E，以相似性分数值大小为序排列，分数越高，相似性越大。而E值则表示随机匹配的可能性，E值越大，随机匹配的可能性也越大。最后给出检测序列和目标序列的比对结果(限于篇幅，图中只给出检测序列和一个目标序列的比对结果)。

图5 BLAST程序运行结果实例

允许空位的 BLAST 最初的BLAST程序只能用于无空位的比对。经验表明比对结果通常会出现一些无空位但不连续的区域，不难想象，有些高分值片段对可以通过一些相似性较低且有空位的片段连接起来，组成了一些更长的或许更具实际生物学意义的比对。 基于上述思路，BLAST算法经过改进允许空位插入（Altshul等，1997）。为缩短对数据库初始搜索的时间，新的算法只找出一个最好的高分值片段，并以此为基础运用动态规划方法将这一片段向两端延伸，最终产生的比对结果可能有空位插入。由于免去了查找所有高分值片段对的步骤，新的算法比原算法快3倍。对BLAST算法的进一步扩充，可以考虑双序列比对和多序列比对的有效结合

位点特异性BLAST叠代搜索 位点特异性BLAST（Position-Specific Iterated BLAST，简称PSI-BLAST）叠代搜索（Altschul等，1997），是一种将双序列比对和多序列比对结合在一起的数据库搜索方法。

位置特异性叠代BLAST (Position-Specific Iterated BLAST，简称PSI-BLAST)则是对蛋白质序列数据库进行搜索的改进，其主要思想是通过多次叠代找出最佳结果。

尽管以下事实已经基本得到认同：基于序列模式的数据库搜索灵敏度较高、特异性较好，因而可以发现一些距离较远但却具有生物学意义的相似序列；它的不足之处也不能予以忽视。除了需要大量的计算资源这一缺点外，对于搜索结果的分析解释常常相当困难。这些制约因素限制了它的实际使用范围。 PSI-BLAST的基本思路在于根据最初的搜索结果，依照预先定义的相似性阈值将序列分成不同的组，构建一个位点特异性的序列谱，并通过多次叠代不断改进这一序列谱以提高搜索的灵敏度。

和其它叠代算法一样，PSI-BLAS方法既有不少长处，也有它的弊病。例如，如果在比对前不把胶原蛋白、同源多聚体等低复杂度的重复序列屏蔽掉，自动叠代搜索过程会因为这些重复序列的干扰而失败（Holm，1998）。假如第一轮的搜索结果出现一个错误序列，那么最终搜索结果中将会出现许多不期望的无关序列。因此，为了尽量去除大量的错误匹配，仔细分析搜索结果给出的同源关系变得非常重要。

BLAST算法 NCBI 算法：做任何事情都有一定的步骤。为解决一个问题而采取的方法和步骤，就称为算法。 BLAST算法：快速高效的保证。 将查询序列分为多个短片段及相似片段； 筛选数据库以发现具备以上片段的序列； 将匹配序列进行延伸，插入和延伸gap，根据突变矩阵（BLOSUM62）计分排序； 返回分值最高的匹配序列

BLAST结果的评价 比对好坏的评价：Bit分值 考虑了比对中相同和相似基团、gap、替代矩阵，并经过标化； Bit分值越高,比对越好 比对统计学意义的评价：E值（E-value） E值越低,则比对就更有可能具有显著性 其他：比对的长度也是一个关键因素

解读BLAST的结果 header。给出查询序列的信息和查询的数据库名称。 每一条匹配序列的描述。包括图形化方式和在线的文字描述。 每个匹配序列与查询序列的比对情况。

BLAST程序的选择 蛋白：BLASTP－tBLASTN 核酸：blastn-blastx-tblastx 数据库的选择：nr最为常用；month跟踪每个月新增数据；swissprot蛋白库注释详尽

比对结果是否有意义的判定 统计学显著性 一致性：蛋白序列>25%，核酸序列>70%（参考） 长度

FastA搜索 FastA算法是由Lipman和Pearson于1985年发表的（Lipman和Pearson，1985）。FastA的基本思路是识别与代查序列相匹配的很短的序列片段，称为k-tuple。

蛋白质序列数据库搜索时，短片段的长度一般是1-2个残基长；DNA序列数据库搜索时，通常采用稍大点的值，最多为6个碱基。通过比较两个序列中的短片段及其相对位置，可以构成一个动态规划矩阵的对角线方向上的一些匹配片段。 FastA程序采用渐进（heuristic approach）算法将位于同一对角线上相互接近的短片段连接起来。也就是说，通过不匹配的残基将这些匹配残基片段连接起来，以便得到较长的相似性片段。这就意味着，FastA输出结果中允许出现不匹配残基。这和BLAST程序中的成对片段类似。如果匹配区域很多，FastA利用动态规划算法在这些匹配区域间插入空位。

由FastA搜索产生的典型输出结果的第一行列出程序名称和版本号，以及该程序发表的杂志。接下来列出所提交的序列，然后是所用参数和运行时间，紧跟这些一般信息的是数据库搜索结果。 首先列出搜索得到的目标序列简单说明，其数目可由用户定义。所列出的目标序列的信息包括：序列所在数据库名称的缩写，目标序列的标识码、序列号和序列名等部分信息。括号中标明匹配部分的残基数。紧接着是由程序计算得到的初始化和优化后的分数值。最后一列是期望值即E值，用来判断比对结果的置信度。接近于0的E值表明两序列的匹配不大可能是由随机因素造成的。

以两条氨基酸序列的比较为例介绍算法的基本思路，算法可以分为4步： 第一步： FASTA首先找出进行比较的两条序列所有长度为K-tuple 的连续的一致序列片段。例如以下两条蛋白质序列： 设K-tuple =2，则序列2中有两个符合条件的片段(用下划线表示)，相对于序列1的偏移(offset)分别是4和1［对于一对开始位置为(x1,x2)的一致片段，偏移定义为x1-x2。在上例中有两对(x1，x2)，即（5,1）和(5,4)］。这种片段的一致性可以表示为对角线图，两条序列中的一对一致片段在图中表示为一段对角线。(图6)。

图6 序列FLWRTW和STWKTWT比较形成的对角线图

本例是两条非常短的氨基酸序列，在实际比较长的蛋白质序列或DNA序列时，对角线图如图A所示。 对于图中每一条完整的对角线(即同一偏移)上的一致片段，如果片段间距小于用户界定的界限，则将片段连接起来作为一条一致片段。.

对这些片段进行计分，每一对对应的元素，一致的加分，不一致的扣分。完成了所有一致片段的计分后，选出10条分值最高的片段进入下一轮计算，如图B

第2步： FASTA将这10对片段重新计分。这轮计分允许保守突变，对蛋白质来就，就是使用PAM250等替换矩阵。简单地说，替换矩阵就是对应于20×20种氨基酸替换(比如R替换成P)的计分规则所构成的20×20的矩阵。这种矩阵是从蛋白质进化实例中总结出来的经验矩阵，它给予进化上相对保守的氨基酸替换比非保守的替换更高的分值。在重新计算分值后，在每一条这样的片段中找出分值最高的子片段，作为“初始区域”(initial region)进入下一步。在initial region中，最高的分值计为initl。

第3步： 在这一步中，FASTA选出分值高于用户确定的界限且相互之间不重叠的初始区域，并尝试将这些初始区域连接起来。当然，由于连接而出现的缺失和插入情况要作相应的扣分。FASTA在这一步才考虑插入和缺失的情况，最终找出能够得到的最高分值的初始区域或连接起来的数个初始区域。这一步计算出的最高分计为initn。见图C

第4步： 以initl片段或(initn的片段)为中心，向前后延伸一定的长度。在这样一个区域中(见图D中虚线间的区域)，应用Smith-Waterman算法进行重新对齐，最终的得分计为opt

在实际操作中，用户可以在需要达到的灵感性程度和所需时间之间进行权衡(一般来说，要达到更高的敏感性总是需要更长的运算时间)，决定采用initn还是opt作为两条序列相似程度的分值。研究表明：使用initn与使用opt相比，前者损失的敏感性并不太大，但运算速度却快得多(Pearson WR,1991)。