58 生物信息学 第一篇生物信息学基础 第 1-3 章两条序列联配算法及序列搜索 序列分析是生物信息学最主要的研究内容之一, 它可以分为两个主要部分 : 一是序列组成分析 ( 包括基因和基因组层次 ), 二是序列之间的比较分析 两条序列或多条序列间的联配或比对, 目的是对它们的序列相似性进行评估, 找出这些序列中结构或功能相似性区域等 通过联配未知序列与已知序列 ( 其功能或结构等已知 ) 的相似程度, 我们可以判断或推测未知序列的结构与功能 第一节序列联配基本概念 序列联配 (sequence alignment) 也叫序列对比, 是生物信息学中的重要内容之一, 许多生物信息学分析均涉及序列联配方法 如下两条 DNA 序列 : >seq1 TGCGGAGC >seq2 TCGGAGC 我们简单地把它们联配如下, 仅有两个碱基匹配 : TGCGGAGC TCGGAGC 如果我们在一条序列中引入一个空位或空格 (gap), 即一条 DNA 序列在进化过程中经常发生的碱基删除事件, 它们的联配结果会显著改善 : TGCGGAGC T-CGGAGC 序列联配我们可以定义如下 : 根据特定的计分规则, 通过一定的算法对两条或多条 DNA 或蛋白质序列进行比较, 找出它们之间最优匹配或最大相似度匹配 序列联配的意义在于, 通过序列比对可以获得一个序列相似性度比对值, 通过这个值的统计学特征分析 ( 见本章第四小节 ), 可以获得两条序列相似度或同源性的一个统计学判断 如果达到显著水平, 说明两条序列相似, 进化上有亲缘关系 在序列分析中, 大量问题涉及这样的相似度判断, 如基因家族 系统发育等分析 根据序列联配的目的不同, 序列联配可以分为全局联配和局部联配两种方式 全局联配 (global alignment) 目的是对两条序列的全长进行比对, 目标是基于它们的全长序列获得最优匹配结果 同样上例两条序列, 以一定计分规则 ( 如碱基匹配得 1 分, 错配罚 1 分, 一个空位罚 3 分 ), 它们的全局最优联配结果如下 :
第 1-3 章两条序列联配算法及序列搜索 59 TGCGGAGC T-CGGAGC (7-3 = 4 分 ) 局部联配 (local alignment) 的目的是获得两条序列比对中得分最高的匹配片段 上例的局部最优联配结果为 : CGGAGC CGGAGC (6 分 ) 上述案例的两条序列很短, 略加比较就可以获得最佳匹配结果 而实际应用中, 序列往往很长, 可能的联配方式很多, 这样就需要一个算法找到两条序列的最佳匹配结果 ( 详见第三节 ) 由此可见, 进行序列匹配, 需要一个计分规则 / 系统 ( 生物信息学专业名称叫计分矩阵 ) 和一个确定最优联配的算法, 同时, 还需要一个统计方法确定序列间的相似程度 以下各节就分别对计分矩阵 联配算法和统计方法进行讲解 第二节计分矩阵 计分矩阵 (scoring matrix) 是序列联配过程中使用的计分规则, 是序列比对的重要组成部分, 它给出序列联配中碱基或氨基酸匹配或错配值, 故又称替换矩阵 ( substitution matrix) DNA 序列相对比较简单, 只有 4 种碱基, 而蛋白质序列有 20 种氨基酸, 如何给出这些氨基酸匹配和错配的一个科学准确评价值, 即准确反映它们的生物学特征, 是生物信息学发展之初就面临的问题, 也是最早被解决的序列联配关键问题 一 计分矩阵的一般原理构建计分矩阵, 我们需要找到一个可以估计任何联配的某一统计数, 使生物学关系最显著的联配统计数最大 先看以下 2 条氨基酸序列的联配情况 如果我们将各残基按相同率处理, 则 2 种联配方式 (a 和 b) 的得分是相等的 (9 个残基中 5 个匹配 ): ( a) TTYGAPPWCS ( b) TTYGAPPWCS TGYAPPPWS TGYAPPPWS 但是联配 (a) 中是一些相对常见的残基 (A P S 和 T) 保持一致, 而联配 (b) 则是有一些相对稀有残基 (W- 色氨酸和 Y- 酪氨酸 ) 相一致 我们需要一个更科学的赋分方法来反映匹配氨基酸间生物学和化学关系 实际联配中,C-C 匹配相对比 S-S 匹配更重要些, 因为半胱氨酸 (C) 是具有非常特殊性质的相对稀有氨基酸, 而丝氨酸 ( S) 则相对常见或普通 同样, D-E 错配值应取正值, 因为这两个残基具有相同的化学性质, 在两条联配的蛋白质序列中能起到相同的功用 但是,V-K 匹配结果则应被罚分, 因为这两个残基毫无相似, 不可能在两条序列中起到一样的作用 用于 DNA 序列联配的替换矩阵相对比较直观 以下是一个常被使用的 DNA 替换
60 生物信息学 第一篇生物信息学基础矩阵 : A C G T A 0.9-0.1-0.1-0.1 C -0.1 0.9-0.1-0.1 G -0.1-0.1 0.9-0.1 T -0.1-0.1-0.1 0.9 矩阵中每个匹配的碱基对均计为 0.9 分, 每个不匹配 ( 错配 ) 的碱基对被罚 0.1 分, 这样, 下面一个联配的得分应为 5 0.9+2 ( -0.1)= 4.3: GCGCCTC GCGGGTC 用于蛋白质联配的替换矩阵相对复杂, 因为没有一个矩阵可以适用各种情况 构建矩阵时, 应考虑不同的蛋白质家族在进化过程中一种氨基酸突变成另一种氨基酸概率的差异, 根据不同的蛋白质家族和预期的相似程度构建不同的替换矩阵 两个最著名的氨基酸替换矩阵分别是 PAM 和 BLOSUM, 它们分别是在 1979 年和 1992 年被提出的 一个重要的概念必须明确 同源性 (homology) 和相似性 (similarity) 是不同的 2 个概念, 不能混淆和混用 2 条序列具有同源性, 意味着这两条序列存在进化方面的关系, 它们从一条共同的祖先序列进化而来 ; 而相似性, 只是表明两条序列间具有一定的相似程度 序列联配计分中另一个重要问题是空位问题 空位处理是针对序列进化过程中可能发生的插入和缺失而设计的 插入和缺失可能只涉及 1 个或多个碱基或残基, 也可能是整个功能域 (domain), 所以, 在进行空位罚值设计时必须反映这些情况 一般用两个参数应用于空位罚值 (gap penalties) 设定, 一个与空位设置 (gap opening) 有关, 另一个与空位扩展 (gap extension) 有关 任一空位的出现均处以空位设置罚值, 而任一空位的扩大则处于空位扩展罚值 对于一个空位长度为 k 的罚值 w K 可用下式表示 : w k = a+bk 其中 a 是空位设置罚值,b 为空位扩展罚值 这两个参数值设置的变化会对联配产生明显影响 ( 表 1-3.1) 表 1-3.1 空位设置和空位扩展罚值对联配的影响空位设置罚值 (a) 空位扩展罚值 (b) 联配效果大大极少插入或缺失, 适用于非常相关蛋白质间的联配 ; 大小少量大块插入, 用于整个功能域可能插入的情况小大大量小块插入, 适用于亲缘关系较远的蛋白质同源性分析如何设定罚值并无明确的理论可循, 大的空位设置罚值配以很小的空位扩展罚值, 被普遍证实是最佳的设定思路 经过多年的试验, 一个合适的空位罚值方式已经被确定下来, 大多数联配程序均对特定的替换矩阵设定了空位罚值的缺略值 (default) 如果使用者希望使用不同的替换矩阵, 则可以根据特殊问题设置合适的空位罚值标准
二 氨基酸替换矩阵 第 1-3 章两条序列联配算法及序列搜索 61 1. 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) 确定了进化过程中一种氨基酸被另一种氨基酸替换的经验数据 他们 收集了大量蛋白质家族序列, 通过比较, 她们共观测到 1572 次替换 事件 以此为基础, 她 们建立了表 1-3.2 的 可观测或可接受点突变矩阵 A (accepted point mutation matrix)( 由于 舍入误差使表中的数值相加不完全等于 1 572) 氨基酸 i 被氨基酸 j 替换的经验次数 ( 记作 A ij ) 可从上表中找到 矩阵 A 可被称为原始 PAM 矩阵 表 1-3.2 氨基酸替换次数表 (Dayhoff 等,1979) A R N D C Q E G H I L K M F P S T W Y R 30 N 109 17 D 154 0 532 C 33 10 0 0 Q 93 120 50 76 0 E 266 0 94 831 0 422 G 579 10 156 162 10 30 112 H 21 103 226 43 10 243 23 10 I 66 30 36 13 17 8 35 0 3 L 95 17 37 0 0 75 15 17 40 253 K 57 477 322 85 0 147 104 60 23 43 39 M 29 17 0 0 0 20 7 7 0 57 207 90 F 20 7 7 0 0 0 0 17 20 90 167 0 17 P 345 67 27 10 10 93 40 49 50 7 43 43 4 7 S 772 137 432 98 117 47 86 450 26 20 32 168 20 40 269 T 590 20 169 57 10 37 31 50 14 129 52 200 28 10 73 696 W 0 27 3 0 0 0 0 0 3 0 13 0 0 10 0 17 0 Y 20 3 36 0 30 0 10 0 40 13 23 10 0 260 0 22 23 6 V 365 20 13 17 33 27 37 97 30 661 303 17 77 10 50 43 186 0 17 注 : 总计观测到 1572 次替换 ; 表中次数均已乘 10; 祖先序列不明时, 次数以平分处理
62 生物信息学 第一篇生物信息学基础 由矩阵 A 可以进一步获得 突变概率矩阵 M ( mutation probability matrix) 矩阵 M 的 元素 M ij 表示经过一定的进化时期氨基酸 j 被氨基酸 i 所替换的经验频率 Dayhoff 等进而把 可观测突变百分率 (percent accepted mutation 或 point accepted mutation per 100 residues), 即 PAM 作为一种时间度量单位 假设同一位点不会发生二次以上的突变, 则 1PAM 等于 100 个氨基酸多肽链中预期发生一次替换所需的时间 Schwarts 和 Dayhoff(1979) 发现将突变概 率矩阵 M 250 次方处理获得的 250PAM 矩阵, 对于研究远缘蛋白质之间进化关系是一个合 适的时间单位 Dayhoff 等 ( 1979) 进一步定义了一个相对概率矩阵 R ( relatedness odds matrix), 表 1-3.3 中各元素已经对数处理, 为对数概率矩阵 (log-odds matrix), 并将最有可能 发生相互替换的氨基酸归类排列 表 1-3.3 PAM250 的对数概率矩阵 (Dayhoff 等,1979) C S T P A G N D E Q H R K M I L V F Y W C 12 S 0 2 T -2 1 3 P -3 1 0 6 A -2 1 1 1 2 G -3 1 0-1 1 5 N -4 1 0-1 0 0 2 D -5 0 0-1 0 1 2 4 E -5 0 0-1 0 0 1 3 4 Q -5-1 -1 0 0-1 1 2 2 4 H -3-1 -1 0-1 -2 2 1 1 3 6 R -4 0-1 0-2 -3 0-1 -1 1 2 6 K -5 0 0-1 -1-2 1 0 0 1 0 3 5 M -5-2 -1-2 -1-3 -2-3 -2-1 -2 0 0 6 I -2-1 0-2 -1-3 -2-2 -2-2 -2-2 -2 2 5 L -6-3 -2-3 -2-4 -3-4 -3-2 -2-3 -3 4 2 6 V -2-1 0-1 0-1 -2-2 -2-2 -2-2 -2 2 4 2 4 F -4-3 -3-5 -4-5 -4-6 -5-5 -2-4 -5 0 1 2-1 9 Y 0-3 -3-5 -3-5 -2-4 -4-4 0-4 -4-2 -1-1 -2 7 10 W -8-2 -5-6 -6-7 -4-7 -7-5 -3 2-3 -4-5 -2-6 0 0 17 表中数值均乘以 10 1PAM 相当于所有的氨基酸平均有 1% 发生了变化, 经过 100PAM 的进化, 并非每个氨 基酸的残基均发生变化 : 有一些可能突变多次, 甚至又变成原来的氨基酸, 而另一些氨基酸 可能根本没有发生过变化 总体上说, 利用大于 100PAM 的时间间隔可能区分进化关系较 远的同源性蛋白质 应该注意,PAM 与进化时间之间没有大致对应关系, 因为不同的蛋白 质家族的进化速率是不同的 当 2 条序列进行相似性比较时, 事先不知道怎样的进化时间
(PAM) 是恰当的 对于相近的序列, 比较容易选择, 即使不太合适的矩阵也无妨 在很多年 里,PAM250( 矩阵后面数字表示一种进化的距离, 数字越大, 进化距离越远 ) 是应用最广的替 换矩阵, 因为该矩阵是唯一由 Dayhoff 最初发表的矩阵 后来一些学者利用大量新出现的蛋 白质序列数据更新 Dayhoff 矩阵的频率数值, 由此构建新的 PAM 矩阵, 但它们与最初的 PAM 矩阵没有太大的差异 2. BLOSUM 替换矩阵 另外一种构建矩阵的方法是由 Henikoff 等于 1992 年提出的, 建成的矩阵称为 BLOSUM (blocks substitution matrix) 他们直接利用多序列联配 (multiple alignment) 分析亲缘关系较 远的蛋白质, 而不是用近缘的序列 这方法的优点是符合实际观测结果, 不足之处是它不能 和进化挂起钩来 大量的试验表明,BLOSUM 矩阵总体比 PAM 矩阵更适合于生物学关系的 分析和局部相似性搜索 假设 f ij 为序列联配中氨基酸 i 和 j 对 ( 忽略顺序 ), 则 i,j 对氨基酸所占比例为 : q ij = f ij / f ij i,j 在完全独立的状况下, 该比例的期望值为 : e ij = p2 i i = j { 2pi p j i j 则 BLOSUM 矩阵元素 (i,j) 定义为 : p i = q ii + 1 2 j iq ij s ij = 2log 2 (q ij / e ij ) 蛋白序列的高度保守区 ( highly conserved regions) 或称为模块 ( block) 数据被用于构建 BLOSUM 矩阵 BLOSUM 矩阵后的数字表示用于构建矩阵模块的最小相似比例, 例如 BLOSUM62 为用于构建矩阵的模块序列数据中, 序列片段的各联配点上至少 62% 是相同的 矩阵后的数字越大, 则表示关系越近 三 位置特异性计分矩阵 (PSSM) PSSM(position-specific scoring matrix) 是由一个简单对数变换而来的矩阵, 它给出不同 来源的一小段保守序列 ( 基序 ) 各个特定位置氨基酸的频率 PSSM 可以用于一条序列的保 守序列的搜索 一条序列中, 与 PSSM 最相似的位置即为 PSSM 代表的基序位置 PSSM 在 数据库搜索, 特别是保守短序列 ( 功能域 ) 搜索方面有很好的应用 例如 PSI - BLAST, DELTA-BLAST 等均使用 PSSM 进行搜索等 ( 详见第四节 ) 下面举例说明 PSSM 构建过程 : 例如如下一个保守序列联配结果 : 由 5 条序列, 每条 6 个碱基构成, 即从左到右 6 个位 置上均有 5 个碱基构成 第 1-3 章两条序列联配算法及序列搜索 AGGCTT AAGCTA AAACTT TAACTA AGACTT 63
64 我们如何将其转换成一个 PSSM? 步骤 1: 联配 6 列位置 上四种碱基数量统计 ; 合计 ( 即背景 ) 统计四种碱基 步骤 2: 各个位置四种碱基频率计算 ; 合计四种碱基频率计算 步骤 3: 标准化 ( 去除背景碱基的影响 ), 即各个位置碱基频率除以背景相应碱基频率 步骤 4: 进行以 2 为底对数 log2 数据转换, 完成 PSSM 构建 #1 #2 #3 #4 #5 #6 合计 A 4 3 3 0 0 2 12 G 0 2 2 0 0 0 4 C 0 0 0 5 0 0 5 T 1 0 0 0 5 3 9 #1 #2 #3 #4 #5 #6 背景频率 A 0.8 0.6 0.6 0 0 0.4 0.40 G 0 0.4 0.4 0 0 0 0.13 C 0 0 0 1.0 0 0 0.17 T 0.2 0 0 0 1.0 0.6 0.30 #1 #2 #3 #4 #5 #6 背景频率 A 2.0 1.5 1.5 0 0 1.0 0.40 G 0 0.3 0.3 0 0 0 0.13 C 0 0 0 5.9 0 0 0.17 T 0.7 0 0 0 3.3 2.0 0.30 #1 #2 #3 #4 #5 #6 A 1.0 0.6 0.6-3.3-3.3 1.0 G -3.3-1.7-1.7-3.3-3.3-3.3 C -3.3-3.3-3.3 2.6-3.3-3.3 T 0.7-3.3-3.3-3.3 1.7 1.0 PSSM 可以直接用于序列搜索, 确定未知序列中与 PSSM 序列相似的保守区段 对于一 条未知序列, 可以用 PSSM 对该序列进行扫描 ( 即从左到右以 PSSM 宽度为窗口宽度, 逐个碱 基位置进行比对 ), 这样就可以获得未知序列中各个位置片段与 PSSM 相似性概率 例如某 一未知序列 : 未知序列 :...A G A C T A... #1 #2 #3 #4 #5 #6 A 1.0 0.6 0.6-3.3-3.3 1.0 G -3.3-1.7-1.7-3.3-3.3-3.3 C -3.3-3.3-3.3 2.6-3.3-3.3 T 0.7-3.3-3.3-3.3 1.7 1.0 利用上述 PSSM(6 个碱基窗口 ) 对该条未知序列进行逐个碱基位点扫描 在上述未知 序列的某一个特定位点, 基于该 PSSM 的几率对数值 (log odds score) 为 (1.0-1.7+0.6+2.6+ 1.7+1.0) = 5.2 生物信息学 第一篇生物信息学基础 5.2 意味着 37 / 1 比率支持未知序列中该 6 个碱基连续片段与 PSSM 并非随机匹配 也就是说, 该未知序列的确包含与该保守区域高度同源的序列片段 以上有关替换矩阵的讨论仅仅提及蛋白质序列的比较, 相关的原则同样适用于 DNA 序
第 1-3 章两条序列联配算法及序列搜索 65 列的比较 如前所述,DNA 替换矩阵非常简单, 所有四个碱基的匹配与不匹配的替换数值均设为相同, 不同的只有匹配与否 ( 匹配 0.9 分和错配罚分 0.1) 一个较复杂的模型可以把转换 ( 两种嘧啶或两种嘌呤间的突变 ) 频率设为高于颠换 ( 嘧啶与嘌呤间的突变 ) 频率 第二节两条序列联配算法一 Needleman-Wunsch 算法 Needleman-Wunsch 算法是一种全局联配算法, 它从整体上分析两个序列的关系, 即考虑序列总长的整体比较, 用类似于使整体相似 (global similarity) 最大化的方式, 对序列进行联配 两个不等长度序列的联配分析必需考虑在一个序列中一些碱基的删除即在另一序列做空位 (gap) 处理 Needleman 和 Wunsch(1970) 最初提出的算法是寻求使两条序列间的距离最小, 或最短距离, 它使用的是一个动态规划 (dynamic programming) 的方法 该算法可以用于核酸和蛋白质序列, 是生物信息学最经典算法之一 在给定计分规则 ( 替换矩阵和空位罚值 ) 情况下, 它们总是能给出具有最高 ( 优 ) 联配值的联配结果 但是, 这个联配结果并不一定具有生物学意义, 因为它可能达不到生物学意义上的显著水平 如果将两条联配的序列沿双向表的上轴和左侧轴放置, 两条序列的所有可能的联配方式都将在它们所形成的方形图中 ( 例举见图 1-3.1a) 图中标出了一条序列所有碱基与另外一条序列碱基所有可能的联配方式 ( 碱基匹配 不匹配和删除即空位 ), 这样所有可能的联配方式都在这个方形图中 从最上角出发, 到右下角结束, 任何一个联配方式均可以画出一条联配路径, 或反过来, 任何一条路径也对应一种联配方式 ( 如图中标注的一个路径及其对应的联配结果 ) 这样一来, 我们确定任意 2 条序列最优联配结果的问题就转化为寻找最优路径问题了 所有可能的路径 ( 联配方式 ) 如果都在这个方形图中, 那么我们如何找到最短路径? 对于任一联配位点, 即图中的任一单元格, 仅有三种可能的方式延伸联配过来 ( 图 1-3.1b):(1) 碱基匹配或不匹配, 即每一序列均加上一个碱基 ( x 路径 ), 并给其增加一个规定的距离权重 ( 匹配加分, 错配罚分 );(2) 在一个序列中增加一个碱基而在另一序列中增加一个空位或反之亦然 (y 和 z 路径 ) 这三种延伸方式的权重值分别加上到达上一个位点的累计得分 (x,y,z), 就可以得到三种可能联配方式的得分, 然后得分最高 ( H) 的路径作为到达本位点的最佳路径 引入一个空位时也将增加一个规定的距离权重 ( 空位罚分 ) 因此, 图中的一个单元可以从 ( 最多 ) 三个相邻的单元达到 为了获得最优路线, 我们必须保证从一开始就每步最优, 即把到达单元格距离最小的方向作为序列延伸的方向 将这些方向记录下来, 并在计算了所有的单元之后, 沿着记录的方向就有一条路径可从方形图右最下角 ( 两个序列的末端 ) 追踪到左最上角 ( 两个序列的起点 ) 由此所产生的路径将给出具有最短距离的序列联配 ( 即最优连配结果 ) 如两条路径获得等距离 ( 相同得分 ), 意味着存在两条最短路径
６６ 生物信息学 第一篇 图 １ ３ １ 生物信息学基础 两条序列联配方式与路径 以两个短序列为例 将上述过程说明如下 ｓｅｑ１ ＣＴＧＴＡＴＣ ｓｅｑ２ ＣＴＡＴＡＡＴＣＣＣ 设定计分规则为 碱基匹配不罚分 ０ 碱基错配时距离权重 罚分 为 １ 引入一个空 位时距离权重为 ３ 两条序列分列于图两侧 图 １ ３ ２ 初始值设为 ０ 然后依次从左上角向右下角计算最 佳路径 根据图 １ ３ １ｂ 每个联配位点均可以计算其 ｘ ｙ ｚ Ｈ 值 例如本例中 第一位点 浅 色框 为 ｘ ｙ ｚ Ｈ ０ ３ ３ ０ 意思是进行碱基联配路径 ｘ 得 ０ 分 不罚分 而在 其中任何一条序列中加入一个空位 ｙ ｚ 均罚 ３ 分 最后到达这个位点的最高得分是进行碱 基联配 即最佳路径 Ｈ 得 ０ 分 初始位点值 ０ 碱基匹配 ０ 而其他两个路径分别得分 ６ 上一个位点得分 ３ 该路径得分 ３ 再看一个位点 深色框 ｘ ｙ ｚ Ｈ １ ３ ３ ３ 其最佳路径是在序列 １ 中插入一个空位 ｚ 其得分为 ３ 上一个位点 ０ ３ 而序列 ２ 中插入空位的路径 ｙ 达本位点最后得分 ９ 上一个位点 ６ ３ 碱 基联配路径 ｘ 得分 ４ 上一个位点 ３ １ 后两个路径值均不及第一个 ｚ 路径 因此 Ｈ ３ 联配过程中可以用箭头标出最佳路径 如此依次类推 就可以最终确定这两条序列最 后一个联配位点的得分 反推得到这个得分的路径 就可以得到其联配方式 而该联配方式 就是这两条序列的最佳联配结果 本例中 沿箭头所指方向在表中从右最下角向左最上角追踪 可以发现 ６ 条路径可以从 初始位点抵达最后一个位点 说明存在 ６ 种最佳联配方式 ＺＤ３９ 生物信息学 兴邦二校 成品 １８５ｍｍ ２６０ｍｍ ３９ 行 ３９ 字 １０ ５ＳＳ
第 １ ３ 章 ６７ 两条序列联配算法及序列搜索 ＣＴＡＴＡＡＴＣＣＣ ＣＴＧＴＡ ＴＣ ＣＴＡＴＡＡＴＣＣＣ ＣＴＧＴＡ Ｔ Ｃ ＣＴＡＴＡＡＴＣＣＣ ＣＴＧＴＡ Ｔ Ｃ ＣＴＡＴＡＡＴＣＣＣ ＣＴＧＴ ＡＴＣ ＣＴＡＴＡＡＴＣＣＣ ＣＴＧＴ ＡＴ Ｃ ＣＴＡＴＡＡＴＣＣＣ ＣＴＧＴ ＡＴ Ｃ 在上述 ６ 种联配中 罚分均为 １０ 即在较短序列中有 ６ 个匹配碱基 １ 个错配碱基和 ３ 个空位 Ｃ Ｔ Ａ Ｔ Ａ Ａ Ｔ Ｃ Ｃ Ｃ ０ ３ ３ ３ ６ ３ ９ ３ １２ ３ １５ ３ １８ ３ ２１ ３ ２４ ３ ２７ ３ ３０ Ｃ Ｔ Ｇ Ｔ Ａ Ｔ Ｃ ３ ３ ３ ９ ３ ３ ３ ６ ３ １２ ３ １５ ３ １８ ３ ２１ ３ ０ ３ ３ ３ ６ ３ ９ ３ １２ ３ １５ ３ １８ ３ ２１ ３ ２４ ３ ２７ ３ ６ ３ ９ ３ １２ ３ ３ ３ ６ ３ ９ ３ ４ ３ ６ ３ １８ ３ １２ ３ １５ ３ １ ３ ４ ３ ７ ３ １０ ３ １３ ３ １６ ３ １９ ３ ２２ ３ ０ ３ ３ ３ ６ ３ ９ ３ １２ ３ １５ ３ １８ ３ ２１ ３ ２４ ３ １５ ３ ９ ３ １２ ３ １ ３ ４ ３ ７ ３ １０ ３ １３ ３ １６ ３ １９ ３ ４ ３ ６ ３ ９ ３ １ ３ ４ ３ ７ ３ １０ ３ １３ ３ １６ ３ ４ ３ ７ ３ ２ ３ ５ ３ ４ ３ ７ ３ １０ ３ １３ ３ ３ ３ ４ ３ ７ ３ １０ 图 １ ３ ２ Ｎｅｅｄｌｅｍａｎ Ｗｕｓｃｈ 算法实例 计分方式 设定碱基错配罚 １ 分 单个碱基缺失或插入时罚 ３ 分 碱基匹配不罚分即得 ０ 分 图中箭头 为每个位点最佳联配方向 来源 每个单元格 ３ 个特别标出举例 内有四个数字分别为 ｘ ｙ ｚ Ｈ 值 Ｈ 值 下划线 图 １ ３ １ 两条序列的全局最优联配 Ｈ 值 １０ 即罚 １０ 分 其联配路径见反向箭头 该算法可以用代数形式来描述 设具有碱基 ａ ｉ 和 ｂ ｊ 的两个序列 ａ 和 ｂ 这两个序列间 距离为 ｄ ａ ｂ 通过评价序列 ａ 中前 ｉ 个位置和序列 ｂ 前 ｊ 位置的距离 ｄ ａ ｉ ｂ ｊ 递归地得 到距离 ｄ ａ ｂ 如果 ａ 和 ｂ 的长度为 ｍ 和 ｎ 则其期望距离为 ｄ ａ ｍ ｂ ｎ 上表中引入的第 ＺＤ３９ 生物信息学 兴邦二校 成品 １８５ｍｍ ２６０ｍｍ ３９ 行 ３９ 字 １０ ５ＳＳ
68 生物信息学 第一篇生物信息学基础 1 行 1 列单元的距离为 0( 相当于空序列或初始位点 ), 在单元 (i,j) 内, 使到达该单元距离增加的三种可能事件为 : 1. 从单元 (i-1,j) 向 (i,j) 的垂直移动, 相当于在 b 序列中插入一个空位使相似序列延伸 换言之,b 序列由 a 序列中 a i 的缺失所产生, 这一事件的权重记作 w_(a i ) 2. 从单元 (i-1,j-1) 向 (i,j) 的对角线移动, 相当于增加碱基 a i 和 b j 使相似序列延伸 换言之,b 序列由 a 序列中的 a i 被 b j 取代所产生, 这一事件的权重记为 w_(a i,b j ) 3. 从单元 (i,j-1) 向 (i,j) 的水平移动, 相当于在序列 b 中插入一个空位使相似序列延伸 换言之,b 序列由 b j 插入 a 序列所产生, 这一事件的权重记为 w + (b j ) 因此, 单元 (i,j) 的距离 d( a i,b j ) 可看成三个相邻单元的距离加上相应权重后的最小者, 即 ìd(a i -1,b j ) + w - (a i ) ï d(a i,b j ) = miníd(a i -1,b j-1 ) + w(a i,b j ) ï îd(a i,b j-1 ) + w + (b i ) 且初始条件为 d(a 0,b 0 )= 0 d(a 0,b j )= j k = 1 d(a i,b 0 )= i k = 1 w + (b k ) w - (a k ) 在图 1-3.2 的实例中 w - (a i )= -3 ( 对于每一个 i) 0 (i = j) w(a i,b j )= { -1 (i j) w + (b j )= -3 ( 对于每一个 j) 当联配两个序列时, 通过计算其随机重排序列的联配距离 ( 如 1 000 次 ), 可以得到这两个序列间的最小距离估计 ( 距离分布 ) 如果实际得到的联配距离小于 95% 的重排序列距离 (1 000 个距离值 ), 则表明实际获得的联配距离达到了 5% 的显著水平, 是不可能由机误造成的 二 Smith-Waterman 算法 Smith-Wateman 算法是在 Needleman-Wunsch 算法基础上发展而来的, 它是一种局部联配算法 由于亲缘关系较远的蛋白质序列可能只有一些相互独立的保守片段, 所以进行局部相似性分析有时可能比整体相似性分析更合理 Smith 和 Waterman(1981) 提出了一种查找具有最高相似性片段的算法 对于序列 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 值,H 的初始值为 : H i0 = 0, 0 i n, H 0j = 0, 0 j m 相似性计算中包括 2 个统计量 : 碱基对 ( 或氨基酸等序列因子对 ) a i,b j 的相似性值 S (a i,b j ) 和空位权重 ( 罚分 )w k = v+uk(k 为空位长度 ) Smith-Waterman 算法可以给出两条序
第 1-3 章两条序列联配算法及序列搜索 69 列的最大相似性值 以 a i,b j 碱基对结束的片段可以由以 a i-1 和 b j-1 结束片段增加碱基 ( 因子 ) 来获得, 或者 a i 可以删除 k 长度的碱基片段,b j 可删除 l 长度碱基片段 具体算法如下 : P ij = max(h i-1,j -w 1,P i-1,j -u) Q ij = max(h i,j-1 -w 1,P i,j-1 -u) ìï H i-1,j-1 +S(a i,b j ) ïp ij = max 1 k i (H i-k,j -w k ) 则 H ij = maxí,(1 i m,1 j n) ïq ij = max 1 l j (H i,j-l -w l ) ï î0 其中 P 0,0 = P 0,j = Q 0,0 = Q i,0 = 0 该算法可以确保具有最大 H ij 值的序列片段是最优联配结果或相似性最好的 从 ( a i, b j ) 为起点, 向后追踪 H ij 矩阵, 直到到达 0 值 对于具有最大相似性片段以外部分的差异性, 不会影响到该片段的 H 值 举例说明这一算法 我们同样以 Needleman-Wunsch 算法中的两条短序列为例 将两条序列 (CTGTATC 和 CTATAATCCC) 排于表 1-3.4 的两侧, 相应的 H ij,p ij 和 Q ij 值分别列入表 中 本例的权重等根据 Smith 和 Waterman(1981) 的例子设定为 : { S(a i,b j )= 1-1 / 3 a i = b j a i b j w k = -(1+k / 3) (1 k) 对于 4 个碱基具有相同频率的随机长序列,S( a i,b j ) 值的平均值为零 w k 值应至少不 小于匹配与不匹配权重的差值 表 1-3.4 的最大 H ij 为 4.33(8 行与 7 列相交处 ), 说明该位点是局部可以达到最高序列 相似度的联配位点, 反推其到达路径, 就得到一条最优路径 ( 星号 表示 ), 其对应具有最大 相似性的片段联配方式为 : CTGTA-TC CTATAATC 表 1-3.4 Smith-Waterman 算法例举 j = 0 j = 1 j = 2 j = 3 j = 4 j = 5 j = 6 j = 7 C T G T A T G i = 0 H ij 0 0 0 0 0 0 0 0 P ij 0 0 0 0 0 0 0 0 Q ij 0 0 0 0 0 0 0 0 i = 1 C H ij 0 1.00 0 0 0 0 0 1.00 P ij 0-0.33-0.33-0.33-0.33-0.33-0.33-0.33 Q ij 0-0.33-0.33-0.67-1.00-1.33-1.33-1.33 i = 2 T H ij 0 0 2.00 0.67 1.00 0 1.00 0 P ij 0-0.33-0.67-0.67-0.67-0.67-0.67-0.33 Q ij 0-0.33-0.67 0.67 0.33 0-0.33-0.33
70 生物信息学 第一篇生物信息学基础 续表 j = 0 j = 1 j = 2 j = 3 j = 4 j = 5 j = 6 j = 7 C T G T A T G i = 3 A H ij 0 0 0.67 1.67 0.33 2.00 0.67 0.67 P ij 0-0.67 0.67-0.67-0.33-1.00-0.33-0.67 Q ij 0-0.33-0.67-0.67 0.33 0 0.67 0.33 i = 4 T H ij 0 0 1.00 0.33 2.67 1.33 3.00 1.67 P ij 0-1.00 0.33 0.33-0.67 0.67-0.67-0.67 Q ij 0-0.33-0.67-0.33-0.67 1.33 1.00 1.67 i = 5 A H ij 0 0 0 0.67 1.33 3.67 2.33 2.67 P ij 0-1.33 0 0 1.33 0 1.67 0.33 Q ij 0-0.33-0.67-1.00-0.67 0.00 2.33 2.00 i = 6 A H ij 0 0 0 0 1.00 2.33 3.33 2.00 P ij 0-1.33-0.33-0.33 1.00 2.33 1.33 1.33 Q ij 0-0.33-0.67-1.00-1.33-0.33 1.00 2.00 i = 7 T H ij 0 0 1.00 0 1.00 2.00 3.33 3.00 P ij 0-1.33-0.67-0.67 0.67 2.00 2.00 1.00 Q ij 0-0.33-0.67-0.33-0.67-0.33 0.67 2.00 i = 8 C H ij 0 1.00 0 0.67 0.33 1.67 2.00 4.33 P ij 0 1.33-0.33-1.00 0.33 1.67 2.00 1.67 Q ij 0-0.33-0.33-0.67-0.67 1.00 0.33 0.67 i = 9 C H ij 0 1.00 0.67 0 0.33 1.33 1.67 3.00 P ij 0-0.33-0.67-0.67 0 1.33 1.67 3.00 Q ij 0-0.33-0.33-0.67-1.00-1.00 0 0.33 i = 10 C H ij 0 1.00 0.67 0.33 0 1.00 1.33 2.67 P ij 0-0.33-0.67-1.00-0.33 1.00 1.33 2.67 Q ij 0-0.33-0.33-0.67-1.00-1.33-0.33 0 第四节 BLAST 算法及数据库搜索 当你面对大量两条序列间的比对时, 运算时间变得非常重要 例如数据库序列搜索就是这样一个问题, 未知序列或递交序列 (query sequence) 需要跟数据库 ( 如 GenBank) 中大量已有序列进行比对, 确定与递交序列的相似序列, 这时就需要在短时间内 ( 如 1 分钟 ) 返回搜索比对结果 直接利用上述两种算法, 将需要大量运算时间, 难以达到时间要求, 必须提出新算法来解决这一问题 Altschul 等人 1990 年提出的用于数据库搜索的 BLAST ( Basic Local Alignment Search Tool) 算法和 1988 年 Pearson 和 Lipman 提出的 FASTA(Fast All) 算法
很好地解决了这一问题 这两个算法 ( 尤其 BLAST) 是目前应用最广泛的方法, 美国生物信 息技术中心 (NCBI) 和欧洲生物信息学研究所 ( EBI) 等分别应用它们来进行 DNA 和蛋白质 序列相似性搜索 这两个算法都是一种基于局部联配搜索工具, 给出的是递交序列与数据 库序列的最佳局部联配结果 两者算法相似, 本节以 BLAST 为例进行讲解 一 BLAST 算法 BLAST 算法同样是利用动态规划算法, 与 Smith-Waterman 算法类似, 其不同之处是引 入了所谓 字 或 字符串 (Word 或 K-tuple,K-mer 等 ) 的检索技术 所有序列其实都是由 若干字符串组成, 例如我们以 3 个碱基长度的字符串为例, 下列 DNA 序列包括了 6 个字长 为 3 的字符串, 其中 2 个字符串 (GCG 和 CGG) 各出现了 2 次 : >seq1 TGCGGAGCGG 其包含的字符串 :TGC,GCG,CGG,GGA,GAG,AGC 为了降低比对时间,BLAST 算法中一个重要手段, 是建立序列数据库的 字 检索系统, 即将数据库中所有序列所包含的不同长度字符串进行扫描, 并建立索引 这样, 数据库中存 贮的序列包含哪些特定长度字符串就已知了, 如常见的 11 个碱基长度 DNA 和 3 个氨基酸 长度蛋白质字符串等 当对递交序列进行数据库搜索比对时, 首先对递交序列进行扫描, 确 定其包含的所有特定长度字符串, 然后仅对包含相同字符串的数据库序列 ( 一般仅占整个数 据库序列很小的比例 ) 进行进一步比对 这样就节省了大量无关序列比对的耗时 对于包 含相同字符串的数据库序列, 进行进一步序列比对时, 基于匹配上的字符串, 分别向两端延 伸, 序列延伸规则基于动态规划算法 下面用图解方式进行具体说明 : 第一步 : 扫描递交序列的每个位点, 发现特定位点 (p) 起始 字长为 w 匹配值超过 T 的所有字并列表 第二步 : 根据上述字列表, 确定目标数据库中所有包含与列表中字完全一致的序列 第 1-3 章两条序列联配算法及序列搜索 71
７２ 生物信息学 第一篇 生物信息学基础 续表 第三步 对于每个匹配的字得 到的联 配 所 谓 ｈｉｔ 向 两 端以动态规划算法向外延伸 延伸过程中联配值 Ｓ 不断 变 动 当联配过程中联配值 Ｓ 降 低超过某一临界值 Ｘ 延伸结 束 这样 我 们 可 以 获 得 所 谓 高 得 分 联 配 对 ＨＳＰ ｈｉｇｈ ｓｃｏｒｉｎｇ ｓｅｇｍｅｎｔ ｐａｉｒ 然后列 出所有 超 过 某 一 设 定 临 界 联 配值或 Ｅ 值 的 ＨＳＰ Ｅ 值 是 一个统计参数 详见下节 表 明随机情况下 可以获得等于 或超过联配 Ｓ 值的 ＨＳＰ 数量 该值越小 说明在统计学意义 上递交序列与得到的 ｈｉｔ 序列 相似性越高 上述对于每个匹配字得到的联配 以动态规划算法向两端延伸 延伸方式是以一个特定 数值 Ｘ 为限定 如果超过该值就终止联配 这种动态规划联配方式是一种数值限定类型 ｓｃｏｒｅ ｌｉｍｉｔｅｄ ＤＰ 图 １ ３ ３ 当然我们也可以用其他条件进行联配限定 如不允许插入 空格 ｕｎｇａｐｐｅｄ ＤＰ 或插入空格数量进行限定 ｂａｎｄｅｄ ＤＰ 如果没有限定 就是我们熟知 的全局联配方式 所谓允许空格的全联配 这些限定 ＤＰ 往往针对特殊问题或目的进行选 用 以提高搜索速度和高效获得特定靶序列 二 利用 ＢＬＡＳＴ 进行数据库序列搜索 采用 ＢＬＡＳＴ 的基本算法目前形成了若干不同的工具 分别用于特定序列数据库和特定 目的的序列搜索 以 ＮＣＢＩ 提供的在线序列数据库搜索工具 ＢＬＡＳＴ ２０１４ 年 １２ 月 为例 图 １ ３ ４ ＢＬＡＳＴＮ 是对核苷酸递交序列库搜索核苷酸序列数据库 ＢＬＡＳＴＰ 是在蛋白质序列 库中搜索蛋白质序列 ＴＢＬＡＳＴＮ 则可以在核酸序列库中搜索蛋白质序列 此时序列库在搜 索之前要按所有 ６ 种读框即时翻译 与此相反的一项分析则由 ＢＬＡＳＴＸ 来完成 它要将所 递交的核酸序列按所有 ６ 种读框翻译 然后再用它搜索蛋白质序列库 同时 特定目的或方 式搜索工具不断被开发出来 图 １ ３ ５ 例如 Ａｌｔｓｃｈｕｌ 等人 １９９７ 提出了一个寻找蛋白质家 ＺＤ３９ 生物信息学 兴邦二校 成品 １８５ｍｍ ２６０ｍｍ ３９ 行 ３９ 字 １０ ５ＳＳ
第 １ ３ 章 图 １ ３ ３ 两条序列联配算法及序列搜索 ７３ 序列动态规划联配方式 可以利用空格 是否允许 空格数量和联配值等进行限定联配过程 族保守序列新算法 ＰＳＩ ＢＬＡＳＴ Ｐｏｓｉｔｉｏｎ Ｓｐｅｃｉｆｉｃ Ｉｔｅｒａｔｅｄ ＢＬＡＳＴ 算法 并开发了相应的 软件 ＰＳＩ ＢＬＡＳＴ 可以对数据库进行多轮循环检索 每一轮的检索速度都大约是 ＢＬＡＳＴ 的 两倍 每一轮都能提高检索的敏感性 图 １ ３ ４ 美国国家生物信息中心 ＮＣＢＩ 提供的在线 ＢＬＡＳＴ 序列搜索工具 ２０１４ 年 １２ 月 ＺＤ３９ 生物信息学 兴邦二校 成品 １８５ｍｍ ２６０ｍｍ ３９ 行 ３９ 字 １０ ５ＳＳ
７４ 生物信息学 第一篇 生物信息学基础 图 １ ３ ５ 数据库 ＢＬＡＳＴ 搜索工具 基于局部动态规划算法 除了标准 ＢＬＡＳＴ 工具 如 ＢＬＡＳＴＰ 等 不断发展了 新的用于蛋白质和核苷酸序列数据的 ＢＬＡＳＴ 工具 下面框中 ＰＳＩ ＢＬＡＳＴ 等 利用 ＢＡＬＳＴ 工具搜索数据库案例 对一条未知 ＤＮＡ 序列 利用 ＢＬＡＳＴＮ 工具搜索核苷酸数据库 获得如下搜索结果 ＺＤ３９ 生物信息学 兴邦二校 成品 １８５ｍｍ ２６０ｍｍ ３９ 行 ３９ 字 １０ ５ＳＳ
第 1-3 章两条序列联配算法及序列搜索 75 上述结果中可见未知 DNA 序列与搜索获得的数据库序列联配得分 (score) 和统计测验结果 (E-value) 等 其中一条数据库记录 (KP751445) 与未知序列具体联配结果如下 :
７６ 生物信息学 第一篇 生物信息学基础 该搜索相关参数列表 该参数列表中 可以看到我们用的搜索参数 默认搜索字长 Ｗｏｒｄ ｓｉｚｅ 为 ２８ｎｔ 结果输 出的默认联配值 Ｅｘｐｅｃｔ ｖａｌｕｅ １０ 调低该值 可以减少低相似度序列的输出 序列联配计 分系统 Ｍａｔｃｈ Ｍｉｓｍａｔｃｈ ｓｃｏｒｅｓ Ｇａｐｃｏｓｔｓ 为匹配 １ 分 错配罚 ２ 分 不考虑空位罚分 下面 ＢＬＡＳＴｘ 为新开空位罚 １１ 分 空位每延伸一个罚 １ 分 数据库中 ３ ４６６ ６ 万条序列进行了比 对 统计测验参数 Ｋａｒｌｉｎ Ａｌｔｓｃｈｕｌ ｓｔａｔｉｓｔｉｃｓ ｋ λ 和 Ｈ 取值情况 详见下节 数据库搜索的 有效序列数据统计结果等列在最后 同样 可以对该未知 ＤＮＡ 序列用 ＢＬＡＳＴＸ 工具进行蛋白质序列数据库搜索 该工具首先按照 ６ 个 ＯＲＦ 阅读框对未知 ＤＮＡ 序列进行翻译 然后比对蛋白质数据库序 列数据 搜索结果如下 ＺＤ３９ 生物信息学 兴邦二校 成品 １８５ｍｍ ２６０ｍｍ ３９ 行 ３９ 字 １０ ５ＳＳ
第 1-3 章两条序列联配算法及序列搜索 77
78 生物信息学 第一篇生物信息学基础可见 BLASTX 与 BLASTN 分别在蛋白质和核苷酸序列水平上进行搜索, 所使用的搜索参数完全不同 前文提及 PSSM, 该计分矩阵在生物信息学领域应用非常广泛, 一个最常见的应用就是 BLAST 搜索中的应用 在蛋白质序列搜索的几个算法中,PSI-BLAST 和 DELTA-BLAST 均使用了 PSSM( 见下图 :NCBI 的 BLAST 主页 ) 这两种算法均利用了一种循环搜索策略, 即利用初次 BLASTP 搜索结果, 临时构建一个新的计分矩阵 ( 即 PSSM), 然后利用该 PSSM 作为下次搜索的计分矩阵, 搜索获得结果后再构建新的 PSSM, 如此循环往复, 直到搜索结果不再变化为止 这样的搜索算法可以最大限度地找到与递交序列具有同源性的序列 ( 有时序列相似性会很低 )
第 １ ３ 章 ７９ 两条序列联配算法及序列搜索 三 序列相似性的统计推断 ＢＬＡＳＴ 搜索返回的结果中 提供了递交序列与数据库中序列比对结果的得分 ｓｃｏｒｅ 和 一个统计测验结果 Ｅ ｖａｌｕｅ 到目前为止 对局部联配的统计学问题已基本搞清楚 特别 是那些不含有空位的局部联配更是如此 我们不妨首先考虑不含有空位的局部联配问题 ＢＬＡＳＴ 最初的搜索程序便是以此为基础的 无空位局部联配涉及的是等长度的一对序列片段 两个片段的各部分彼此比较 Ｓｍｉｔｈ Ｗａｔｅｒｍａｎ 算法方法可以找到所有最高比值片段对 ＨＳＰ 即这些片段对的比值 Ｓ 不会因 片段的延伸而进一步升高 为了分析上述分值随机性产生的几率大小 需要建立一个随机 序列模型 对于蛋白质而言 最简单的序列模型可通过从一条序列中随机地选取氨基酸残 基 当然这一条序列中各种残基的频率必须一定 另外 一对随机氨基酸的联配期望值必 需为负值 否则不论联配片段是否相关 都会得到高比值 统计理论也将派不上用场 就像独 立 随 机 变 量 之 和 总 是 倾 向 于 正 态 分 布 ｎｏｒｍａｌ ｄｉｓｔｒｉｂｕｔｉｏｎ 一样 独立随机变量的最大值倾 向于极值分布 ｅｘｔｒｅｍｅ ｖａｌｕｅ ｄｉｓｔｒｉｂｕｔｉｏｎ 在进行两 条序列最佳局部联配结果统计测验时 主要涉及的是 后一种情况 对于两条序列 在一定的序列长度 ｍ 和 ｎ 限定 下 ＨＳＰ 的统计值可由 ２ 个参数 ｋ 和 λ 确定 最简 单的形式 即不小于比分值为 Ｓ 的 ＨＳＰ 个数 可由下 列公式估计其期望值 Ｅ ｋｍｎｅ λｓ １ ３ １ 我们称该期望值为比值 Ｓ 的 Ｅ 值 Ｅ ｖａｌｕｅ 上述公式非常灵敏 在给定比值的情况下 将两 图 １ ３ ６ 概率密度函数正态分布 虚线 和极值分布 实线 比较 条比较序列长度加倍 则 ＨＳＰ 数 即 Ｅ 值 也将加倍 同样 Ｓ 值为 ２Ｘ 的某个 ＨＳＰ 其长度必是 Ｓ 值为 Ｘ 的 ＨＳＰ 的两倍 所以 Ｅ 值将随着 Ｓ 值的 增大急剧减少 参数 ｋ 和 λ 可分别被简单地视为搜索空间 ｓｅａｒｃｈ ｓｐａｃｅ ｓｉｚｅ 即数据库序列 数据量 和计分系统的特征数 最初获得的比值 Ｓ 在没有计分系统 或统计量 ｋ 和 λ 的辅助下 没有什么意义 单独 的比值就如同没有单位 米或者光年 的距离 可将比值按下式标准化 λｓ ｌｎｋ Ｓ ｌｎ２ １ ３ ２ 获得 Ｓ 值就如同得到了具有标准单位的数值 Ｅ 值因此可简化为 Ｅ ｍｎ２ Ｓ １ ３ ３ 二进制值 ｂｉｔ ｓｃｏｒｅ 使所用的计分系统赋予了统计学意义 除了可以确定搜索空间外 同样可以计算相应的显著水平 具有大于或等于某一比值 Ｓ 的随机 ＨＳＰ 数可由泊松分布确定 由此可以计算出搜索 到某一比值大于或等于 Ｓ 的 ＨＳＰ 的概率为 ＺＤ３９ 生物信息学 兴邦二校 成品 １８５ｍｍ ２６０ｍｍ ３９ 行 ３９ 字 １０ ５ＳＳ
80 e -E E X X! 的概率为 式中 E 由 (1-3.1) 式确定 (1-3.4) 作为一个特例, 搜索不到比值 S 的 HSP 概率为 e -E,, 所以至少发现一个 HSP 比值 S P = 1-e -E = 1-exp( -kmne -λx ) (1-3.5) 这是与比值 S 相关的 P 值 ( 概率值 ) 例如, 在可能搜索到 3 个 HSP 比值 S 的情况下, 至少发现一个 HSP 的机率为 0.95[ 可由 (1-3.5) 式算得 ] BLAST 程序中使用了 E 值而非 P 值, 这主要是从直观和便于理解的角度考虑 比如 E 值等于 5 和 10, 总比 P 值等于 0.993 和 0.999 95 更直观 但是当 E<0.01 时,P 值与 E 值趋于相同 E 值计算公式 [ 公式 (1-3.1)] 可以直接应用于 2 个蛋白质序列长度分别为 m 和 n 的比 较, 但是对于某一序列长度为 m 的蛋白序列, 如何在那些长短不一的数据库序列中找到与之 匹配良好的序列呢? 一种思路是把数据库中的所有蛋白序列与递交序列的关系都视为相同 重要, 也就是说对于 E 值均较低的短或长序列, 它们是等同重要的 FASTA 程序便是采用这 一策略 另一种思路是把长序列视为比短序列更重要, 因为长序列往往包括更多的特异功 能域 (domain) 如果对序列长度上进行相关优先处理, 则在计算数据库序列长度为 n 的 E 值时, 将乘以 N / n, 其中 N 为数据库中序列的总长度 根据公式 (1-3.1),E 值的计算可简单 地把整个数据库序列视为长度为 N 的单条序列 BLAST 程序采用了这一策略 FASTA 策 略中 E 值的计算还需再乘上数据库的序列条数 如果考虑到核酸数据库的序列长度变化更 大, 则在 DNA 序列相似性搜索时,BLAST 的策略可能会是合理的选择 一些数据库搜索程序 ( 例如 FASTA 或其它基于 Smith-Waterman 算法的程序 ) 在进行序 列搜索时, 会对数据库中的每条序列进行联配并给出联配值, 这些值大部分与递交序列无 关, 但它们被用于了 k 和 λ 参数的估计 这一方法避免了因使用真实序列 (real sequence) 导 致随机序列模型的偏向性, 但同时产生了使用相关序列估计参数的难题 BLAST 仅通过部 分而不是全部无关序列计算最适联配值, 这赢得了搜索速度 因此, 对于某一选定的替换矩 阵和空位罚值, 必须进行 k 和 λ 参数的预先估计, 估计中使用真实序列, 而非通过随机序列 模型产生的模拟序列 根据统计理论, 上述统计方法只适用于不含有空位的局部联配 ( 非空位联配 ) 但是, 许 多计算试验和分析结果充分证明, 上述统计方法同样适用于包括空位的联配结果 对于非 空位联配, 可用基于替换矩阵和比较序列的残基频率的办法, 估计统计参数 ; 对于空位联配, 参数的估计则必须根据大量无关序列的比较 以上统计学方法对于短序列来说有些偏差 这些统计方法的基础理论是一个渐近理 论, 该理论假设局部联配可以适用于任何规模的联配 但是, 一个高比值联配必须有一定的 长度, 不能从接近二条序列末端的地方开始 这种边际效应 (edge effect) 可以通过计算序列 的 有效长度 (effective length) 来修正 BLAST 程序中包含了这一修正过程 对于长于 200 残基的序列可以不进行边际效应的修正 生物信息学 第一篇生物信息学基础 局部联配的结果与所选用的替换矩阵紧密相关 没有任何一个计分系统 ( 即替换矩阵 + 罚分办法 ) 可以适用于所有研究目标, 对于局部联配的计分基础理论的正确理解可以极大促 进序列分析准确性 ( 相关内容详见本章第二节 )
第 1-3 章两条序列联配算法及序列搜索 81 如前所述,BLAST 算法中空位罚值应用 2 个参数, 一个与空位设置 ( gap opening) 有关, 另一个与空位扩展 (gap extension) 有关 任一空位的出现均处以空位设置罚值, 而任一空位的扩大必须加入空位扩展罚值 经过多年的试验, 一个合适的空位罚值已经被确定下来 大多数联配程序均对特定的替换矩阵设定了空位罚值的默认值 ( default), 如果使用者希望使用不同的替换矩阵, 则原来的空位罚值设定不一定合适 如何设定罚值并无明确的理论可遁, 较大的空位设置罚值配以很小的空位扩展罚值 ( 如 11 / 1) 被普遍证实是最佳的设定思路 在上节 BLAST 搜索结果中, 有一个统计量 H 它用于估计 BLAST 搜索所用计分矩阵相对信息量或熵值 ( 具体见第 1-4 章具体介绍 ),H 值越大, 说明所用替换矩阵或计分系统的特异性越高, 区分相关与不相关序列的能力就越强 (Mount,2004) 习题 1. 目前计分矩阵主要有哪些? 比较它们的异同 2. 请利用动态规划 Needleman-Wunsch 算法对下列两条蛋白质序列进行全局联配, 获 得最优联配结果 : P1 = AGWGAHEA P2 = PAWHEAEAG 计分系统 : 计分矩阵 BLOSUM50, 空位罚 8 分 BLOSUM50 ( 部分 ) A E G H P W A 5-1 0-2 -1-3 E 6-3 0-1 -3 G 8-2 -2-3 H 10-2 -3 P 10-4 W 15 3. 数据库搜索同源性或相似性较远序列时, 为什么用蛋白质序列搜索比 DNA 序列 要好? 4. 请解释 BLAST 搜索结果中 score 和 E value 含义 5. BLASTN 和 BLASTX 两者搜索方式有何不同?