欢迎关注：oddxix

本章主要讲序列分析与联配

序列分析是生物信息学最主要的研究内容之一，它可以分为两个主要部分：

一、是序列组成（特别是涉及到基因组层次上）分析。

二、是序列之间的比较分析。

两条序列或多条序列间的比对或联配(alignment)的目的，是对它们的序列相似性进行评估，找出这些序列中结构或功能相似性区域等。通过联配未知序列与已知序列(其功能或结构等已知)的相似程度，我们可以判断或推测未知序列的结构与功能。

序列组成及单一序列分析

在基因每一侧的500 个任意碱基区域被称为“侧翼”，基因间区域是指两个基因间的其余序列，DNA序列内和序列间碱基具有不同的频率。

分析 DNA 序列的主要困难之一是碱基相邻的频率不是独立的。碱基相邻的频率一般不等于单个碱基频率的乘积：如果 P u 是序列中碱基 u 的频率，且 P uv 为两个相邻碱基 u 和 v 的频率，则

Puv≠Pu*Pv

在编码区，存在某种约束来限制 DNA 序列编码氨基酸。在密码子水平上，这一约束与碱基相邻频率有关。编码同一氨基酸的不同密码子(同义密码子)好像不是等同存在的。这种密码子偏倚必定与两碱基相邻频率水平有关。

相邻碱基之间的关联将导致更远碱基之间的关联，这些关联延伸距离的估计可以从马尔科夫链(Markovchain)理论得到。在不援引任何生物学机制的情况下，第 k 阶马尔科夫链假定在序列中某一位置上碱基的存在只取决于前面 k 个位置上的碱基。一阶链假定一个特定碱基存在于位置 i的概率只取取决于在位置 i-1 的 4 种碱基概率。相互独立的碱基所组成的序列将与 0 阶马尔科夫链相对应。阶可以通过似然法估计。同时，马尔科夫链分析更适应于基因组水平，而非单一序列(基因)。

Needleman-Wunsch 算法是一种整体联配(global alignment)算法，最佳联配中包括了全部的最短匹配序列。Smith-Wateman算法是在 Needleman-Wunsch 算法基础上发展而来的，它是一种局部联配(Local alignment)算法。这二种算法均可以用于核酸和蛋白质序列。在给定空位罚值和替换矩阵情况下，它们总是能给出具有最高（优）联配值的联配。

从整体上分析两个序列的关系，即考虑序列总长的整体比较，用类似于使整体相似(global similarity)最大化的方式，对序列进行联配。两个不等长度序列的联配分析必需考虑在一个序列中圈掉一些碱基或在另一序列作空位(gap)处理。

Needleman-Wunsch 算法

Needleman-Wunsch 算法是为氨基酸序列发展的，但也可以用于核苷酸序列。算法最初寻求的是使两条序列间的距离最小。尽管这类距离的元素是以一种特定的方式定义的，但该算法的良好特性在于它确定了最短距离。这是一个动态规划(dynamic programming)的方法。

将两条联配的序列沿双向表的轴放置，两条序列的所有可能的联配方式都将在它们所形成的方形图中。从任一碱基对，即表中的任一单元开始，联配可延三种可能的方式延伸：如果碱基不匹配，则每一序列加上一个碱基，并给其增加一个规定的距离权重；或在一个序列中增加一个碱基而在另一序列中增加一个空位或反之亦然。引入一个空位时也将增加一个规定的距离权重。因此，表中的一个单元可以从(至多)三个相邻的单元达到。我们把达左上角单元距离最小的方向看作相似序列延伸的方向。等距离时意味着存在两种可能的方向。将这些方向记录下来，并在研究了所有的单元之后，沿着记录的方向就有一条路径可从右下角(两个序列的末端)追踪到左上角 (两个序列的起点)。由此所产生的路径将给出具有最短距离的序列联配。

Smith-Waterman 算法

由于亲缘关系较远的蛋白质序列可能只有一些相互独立的相同片段，所以进行局部相似性分析有时可能比整体相似性分析更合理。对 于序 列 A=(a 1 ,a 2 , … ,a m ) 和B=(b 1 ,b 2 , … ,b n )，H ij 被定义为以a i 和b j 碱基对结束的片段(亚序列)的相似性值。与Needle-Wunsch算法一样，Smith-Waterman算法也要利用递推关系来确定H值。

相似性计算中包括 2 个统计量：碱基对(序列因子) 的相似性值和空位权重 (k 为空位长度)。Smith-Waterman算法可以给出2 条序列的最大相似性值。以 碱基对结束的片段可以由以 和 结束片段增加碱基(因子)来获得，或者 可以删除 k 长度的碱基片段，可删除 l 长度碱基片段。

该算法可以确保具有最大H ij 值的序列片段是相似性最好的。从 为起点，向后追踪矩阵，直到到达某一负值。对于具有最大相似性片段以外部分的差异性不会影响到该片段的H值。

序列相似性的统计特性

先考虑不含有空位的局部联配问题

无空位局部联配涉及的是等长度的一对序列片段，两个片段的各部分彼此比较。一种Smith-Waterman 或 Sellers 算法的改进算法可以找到所有高比值片段对(high-scoring segment pairs,HSPs)，即这些片段对的比较分值不会因片段的延伸而进一步升高。

为了分析上述分值随机性产生的几率大小，需要建立一个随机序列模型。对于蛋白质而言，最简单的序列模型可通过从一条序列中随机地选取氨基酸残基，当然这一条序列中各种残基的频率必需一定。另外，一对随机氨基酸的联配期望值必需为负值，否则不论联配片段是否相关的，都会得到高比值。

在一定的序列长度 m 和 n 限定下，HSP 的统计值可由 2 个参数(k 和λ)确定。最简单的形式，即不小于比较值为 S的 HSP 个数，可由下列公式算得其期望值：

在给定比值的情况下，将比较序列长度加倍，则 HSP数(即 E 值)也将加倍，同样，S 值为 2X 的某个 HSP 长度必是 S 值为 X 的两倍，所以 E 值将随着 s 值的增大急剧减少。参数 K 和λ可分别被简单地视为搜索步长(search spacesize)和计分系统(scoring system)的特征数。

P 值(P-Value)(概率值)

具有大于或等于某一比值 S 的随机HSP 数可由泊松分布(Poisson distribution)确定。由此可以计算出搜索到某一比值大于或等于 S 的 HSP 的机率为

空位联配(gapped alignment)的统计问题

对于非空位联配，可用基于替换矩阵和比较序列的残基频率的办法估计统计参数；对于空位联配，参数的估计则必须根据“随机”序列的大尺度比较。

空位罚值(gappenalties)

联配中另一个重要问题是空位问题。空位处理是针对序列进化过程中可能发生的插入和缺失而设计的。插入和缺失可能只涉及 1 个或 2 个残基，也可能是整个功能域(domain)，所以，在进行空位罚值设计时必须反映这些情况。有 2 个参数应用于空位罚值设定，一个与空位设置(gap opening)有关，另一个与空位扩展(gap extension)有关。任一空位的出现均处以空位设置罚值，而任一空位的扩大必须处于空位扩展罚值。

对于一个空位长度为k的罚值W K 可用下式表示：

其中 a 是空位设置罚值，b 为空位扩展罚值。这两个参数值设置的变化对联配产生影响

替换矩阵的一般原理

我们并不能直接计算出两条序列的最佳联配，我们需要找到一个可以估计任何联配的某一统计数，使生物学关系匹配最显著的联配统计数最大。替换矩阵(substitution matrices)包括了在联配中各种匹配方式如何赋分的信息，故替换矩阵又常被称为计分矩阵(scoring matrices)。

用于蛋白质联配的替换矩阵要复杂一些，因为没有一个矩阵可以适用各种情况。构建矩阵时应考虑不同的蛋白质家族在进化过程中，一种氨基酸突变成另一种氨基酸概率的差异，根据不同的蛋白质家族和预期的相似程度构建不同的替换矩阵。2 个最有名的蛋白质替换矩阵是 PAM和 BLOSUM

PAM 氨基酸替换矩阵

在进行蛋白质序列联配时，必须通过一定的方法给联配的残基对赋予一定的分值，替换矩阵便是其中最重要的方法。已故 Dayhoff 是蛋白质列序比较的先驱，她和她的同事们通过对蛋白质进化模式的研究，建立了一组被广泛应用的替换矩阵，这些矩阵常被称为Dayhoff，MDM(Mutation Data Matrix)或 PAM(Percent Accepted Mutation)矩阵。

由于蛋白质最有可能是自然选择的目标，可以认为蛋白质序列的分析比DNA分析更具有生物学意义。蛋白质分析完全避免了几个三联体可能编码同一氨基酸的遗传密码简并问题。有必要进一步分析各种氨基酸间的同源性程度，以及在进化过程中一种氨基酸被另一种氨基酸替换的概率大小。也许把氨基酸按一定特性分成若干组更便于以上分析，例如氨基酸可分成中性疏水(G、A、V、L、I、F、P、M)、中性亲水(S、T、Y、W、N、E、C)、碱性(K、R、H)和酸性(D、E)氨基酸等。在比较许多具有相似性蛋白质序列的基础上，Dayhoff等于 1979 年构建了一个突变概率矩阵M(mutation probability matrix)。

BLOSUM 氨基酸替换矩阵

另外一种构建矩阵的方法是由 Henikoff 等于 1992 年提出的，建成的矩阵为BLOSUM(Blocks SubstitutionMatrices)。他们直接利用多序列联配(multiple alignment)分析亲缘关系较远的蛋白质，而不是用相近的序列。这方法的优点是符合实际观测结果，不足之处是它不能和进化挂起钩来。大量的试验表明，BLOSUM矩阵总体比 PAM 矩阵更适合于生物学关系的分析和局部相似性搜索。

DNA 替换矩阵

以上有关替换矩阵的讨论仅仅提及蛋白质序列的比较，但是，相关的原则同样适用于 DNA 序列的比较。在进行比较时应该意识到，用翻译而来的蛋白质序列总是好于直接用 DNA 序列。这是因为 DNA 序列的进化变化很少，在使用简单的DNA 替换矩阵比较时，获得的同源性信息远少于蛋白质序列。DNA 替换矩阵非常简单，所有 4 个碱基的匹配与不匹配的数值均设为相同，不同的只有匹配与否(0.9 和-0.1)。一个较复杂的模型是把转换(transition,两种嘧啶或两种嘌呤间的突变)频率设为高于颠换(tranversion,嘧啶与嘌呤间的突变)频率。

多序列联配

通过以上的两条序列算法，总是可以返回一个最佳匹配的联配结果。但是，当我们将两条以上的序列放在一起联配时，情况就就不一样了。现有实用的多序列联配方法还不能保证一定给出最优联配结果，只能给出一个近似值——往往人为的修正可以使联配结果更佳。

同源序列的多序列联配是生物信息学一个重要课题。通过多序列联配结果，允许你观察残基可以改变到什么程度而蛋白质仍保持功能；它也可以使你得到围绕某一残基的三级结构信息。有不少多序列联配程序可通过匿名 ftp 等服务获得，例如：ClustalW等。三条或三条以上序列的联配方法可分为几类，如用于两条序列联配的Needleman-Munsch 等算法的改进算法、等级法(hierachical method)、片段法(segment method)、一致或区段法(consensus or regions'method)等。这些方法中，等级法是目前应用最为广泛的方法。

等级法又称为树法(tree method)，是由 Feng 和 Doolittge(1987)等人发展的（ClustalW 程序）。由于两条序列的联配结果可以很容易地获得，多序列联配便可以在连续使用两条序列联配算法(如 Needleman-Wunsch 算法)基础上，通过先建“树”的思路来进行多序列联配。这一方法同样是一种动态规划方法。具体步骤如下：

• ①对所有序列进行两两联配分析，N 条序列应有 N×(N-1)/2 对；
• ②对两两联配的数据进行聚类分析，产生联配等级。该等级可用分叉树(binary tree)形式或简单的排序来表示；
• ③根据以上联配结果，首先从所有联配中相似性最好的两条序列开始，然后是剩余联配中相似性最好的两条序列⋯⋯依次类推，直至多序列联配结束。一旦两条序列的联配被列入，则序列的位置就被固定下来。例如，对于序列 A、B、C、D，如果 A 与 C、B 与 D 分别是两两联配的最佳联配结果，则 A、B、C、D 四条序列的联配则通过比对 A-C 和 B-D 两个联配(每个联配位置取平均值)来确定。这一组合方法对大量序列的多序列联配提供了实用的空位联配手段，除了最初的两序列间的联配过程，整个多序列联配过程是很快的。

数据库搜索——BLAST 和 FASTA 应用

比较和确定某一数据库中的序列与某一给定序列的相似性是生物信息学中最频繁使用和最有价值的操作。本质上这与两条序列的比较没有什么两样，只是要重复成千上万次。但是要严格地进行一次比较必定需要一定的耗时，所以必需考虑在一个合理的时间内完成搜索比较操作。目前有二个最为常用的程序服务于未知序列的数据库相似性搜索，即 BLAST 和 FASTA。FASTA 使用的是Wilbur-Lipman 算法的改进算法，进行整体联配，重点查找那些可能达到匹配显著的联配。虽然 FASTA 不会错过那些匹配极好的序列，但有时会漏过一些匹配程度不高但达显著水平的序列。BLAST(Basic Local Alignment Search Tool，基本局部联配搜索工具)是基于匹配短序列片段，用一种强有力的统计模型来确定未知序列与数据库序列的最佳局部联配。

大多数研究目前都通过国际互联网 Internet 应用 NCBI 研制的 BLAST 程序(Basic Local Alignment Search Tool)来进行 DNA 和蛋白质序列相似性搜索。用一组 BLAST 程序联配可以快速进行核酸和蛋白质序列库的相似性检索。
采用BLAST 的基本算法编成了若干各不同的程序，分别使用特定的序列库和用于特定类型的输入序列。BLASTN 是在核苷酸序列库搜索核苷酸序列。BLASTP 是在蛋白质序列库中搜索氨基酸序列。TBLASTN 则可以在核酸序列库中搜索氨基酸序列，此时序列库在搜索之前要按所有 6 种读框即时翻译。与此相反的一项分析则由BLASTX 来完成，它要将所输入的核酸序列按所有 6 种读框翻译，然后再以之搜索蛋白质序列库。PSI-BLAST 可以对数据库进行多轮循环检索，每一轮的检索速度都大约是 BLAST 的两倍，但每一轮都能提高检索的敏感性。它是目前 BLAST 程序家族中敏感性性最高的成员。

如果目的序列中有蛋白质编码区，则用翻译的蛋白质序列来搜索蛋白质序列库要比用 DNA 序列搜索核酸序列库更有价值。由于蛋白质序列的进化要比 DNA序列慢一些，在蛋白质序列水平上的远缘关系在 DNA 水平上可能被错过。如果无法确定编码区，则可利用 BLASTX 按所有 6 种读框来翻译 DNA 序列，然后用它搜索蛋白质序列库。由于蛋白质序列库仅包含已鉴定的蛋白质，所以必须采用TBLASTN 程序在现有的 GenBank、EMBL 或 DDBJ DNA 序列库中检索新确定的氨基酸或翻译过的 DNA 序列。这种检索有时可以找到一些显著相似的 DNA 序列，而原本并不知道这些序列可编码蛋白质。

另一个常用的核酸和蛋白质序列库搜索程序是 FASTA，即 FASTN 和 FASTP 程序的新版本。FASTA 首先在序列库中进行快速的初检，找出与待检序列高度相似的序列。这一快速检索局限于待检序列和序列库序列之间较短的完全相同序列区段上。

无论采用 FASTA 或 BLAST，推断相似性是否具有生物学意义都取决于研究者。要作出决断，必须充分考虑蛋白质已知的或推断的功能，与已知活性位点或模序的相似程度等等。因为BLAST和FASTA采用不同的算法，同时用这两种搜索引擎重新检索某一特定序列往往是可取的。如果用其中一种找不到显著相似序列，不妨试一试另一程序。如果BLAST和FASTA均找不到显著匹配的序列，还可以选择第 3 条比较费时的搜索策略。一些网站允许用户使用基于Smith-Waterman算法的搜索程序，如BLITZ。BLITZ( www.ebi.ac.uk/searchs/blitz.html) 被设计在大型并行计算机上运行，因此使检索更灵敏。虽然运行这样的程序比较费时，但它们有时会发现一些被BLAST和FASTA错过的勉强达到显著的联配。