2.6 序列比对与比较基因组学

2.5 节的注释流程在几乎每一个环节都重复同一个动作：把一条新序列与一条已知序列并排放在一起，逐位指出“这两个字母对应”。这个动作就是序列比对(sequence alignment)，它是基因组科学中使用频率最高的计算操作，也是本章从“读出序列”走向“读懂序列”的枢纽。比对之所以需要认真的理论，是因为两个问题远比看上去困难：其一，什么样的比对才算“好”，需要一套客观的打分体系；其二，如何在可以接受的时间内，在数以亿计的候选位置中找到最好的比对。本节先把比对问题形式化（2.6.1），再讨论工业上最重要的启发式工具 BLAST（2.6.2–2.6.3），最后把两两比较放大到整个基因组，进入比较基因组学（2.6.4–2.6.5）。

2.6.1 比对问题的形式化：打分与最优解

给定长度为 m 与 n 的两条序列 x、y，一个比对是指：在任一序列的任意位置插入空位(gap)“−”，使两行长度相等，从而把同一列的两个字符置于相互对应的关系。空位并不存在于真实序列之中，它是历史上一次插入/缺失(insertion/deletion, indel)事件留下的痕迹；比对因此不是排版游戏，而是关于两条序列共同历史的一个假设。既然是假设，就有优劣之分，而优劣必须可以度量。

打分体系由两部分组成：列内残基对的替换得分，以及对空位的处罚。替换得分由替换矩阵(substitution matrix)给出：矩阵元素 s(a, b) 表示把 a 与 b 放在同一列应得的分数。核苷酸序列只有四种字符且多数应用只关心是否配对，常用“匹配 +1、错配 −m”的简单方案（blastn 缺省即属此类），更精细的方案可对转换（嘌呤间或嘧啶间）与颠换（嘌呤对嘧啶）设置不同罚分，因为转换在进化上远为常见。氨基酸的替换倾向则悬殊得多：理化性质相近的替换（如异亮氨酸与缬氨酸）频繁且常被选择容忍，而色氨酸突变为半胱氨酸一类的事件罕见且多受净化选择排斥。因此蛋白比对必须使用经验矩阵：PAM 系列由 Dayhoff 等自密切同源蛋白中观察到的接受点突变(accepted point mutation)计数外推而来；BLOSUM 系列由 Henikoff 与 Henikoff 从未经比对的保守序列块(blocks)中直接统计替换频率而得（Henikoff & Henikoff, 1992）。两个家族共享同一个核心思想——以“取代的可观测性”计分：矩阵元素本质上是观测替换频率对随机背景频率的比值取对数（对数几率，log-odds），常见替换得正分，罕见替换得负分，从而使总分与“两序列同源与否”的证据强度挂钩。其中 BLOSUM62 由块内相似度不低于 62% 的序列先合并成簇再统计，是蛋白数据库检索的缺省矩阵。

空位的处罚采用仿射空位罚分(affine gap penalty)：开设一个新空位付出较重的开放罚分(gap-open penalty) d，此后空位每延长一位只付出较轻的延伸罚分(gap-extension penalty) e，故长度为 L 的一个空位共罚 d + e·L。这一设计的生物学动机十分明确：一次真实的插入/缺失事件——复制滑移、不等交换或转座子插入——往往一次牵连多个碱基；如果按“每位重罚”的线性罚分，把一个长空位拆散成许多单碱基空位反而“更便宜”，最优比对将系统性地背离进化事实。综合两部分，比对 A 的总分定义为：

求解式 (2.6-1) 最优比对的经典框架是动态规划(dynamic programming)。全局比对(global alignment)要求两条序列从头到尾全部参与，由 Needleman 与 Wunsch 于 1970 年提出：以 F(i, j) 记“x 的前 i 位与 y 的前 j 位的最优比对得分”，则 F(i, j) 只依赖三个更小的子问题 F(i−1, j−1)、F(i−1, j)、F(i, j−1)，按行列填满整张表后回溯即得最优比对（Needleman & Wunsch, 1970）。局部比对(local alignment)由 Smith 与 Waterman 于 1981 年给出：在递推式中增加取 0 的一支，使比对可以在任意位置开始与结束，最终只保留整张表中得分最高的一条路径（Smith & Waterman, 1981）。局部比对对应数据库检索的常态——查询与库序列往往只有一段同源，例如只共享一个结构域，强迫首尾对齐的全局比对反而会被无关区段淹没。仿射空位罚分需要三条并行的递推线，复杂度量级不变（Gotoh, 1982）。无论全局还是局部，填表需要约 m×n 个格子、每格常数次比较，故时间与空间均为 O(mn)。

把这笔账具体化：一条 500 aa 的查询对总计 10⁹ aa 残基的数据库做逐条 Smith–Waterman，需要约 5×10¹¹ 次格子更新。精确算法在离线比对两条指定序列时毫无问题，却撑不起“任何实验室随时提交一条序列、数秒内得到结果”的在线检索。理论与工程之间的这道缺口，由下一小节的启发式填补。

2.6.2 BLAST：种子—扩展的启发式加速

BLAST（基本局部比对搜索工具）(Basic Local Alignment Search Tool)（Altschul et al., 1990）建立在一条生物学直觉之上：一段足够显著的同源片段，几乎总是含有若干短而无空位的精确匹配。据此它把检索分成两步，即种子—扩展(seed-and-extend)策略。

第一步播种：把查询序列切成所有长度为 k 的字(word)（核酸缺省 k = 11，蛋白缺省 k = 3），在数据库序列中快速定位完全相同的字；对蛋白还进一步考虑邻字(neighborhood words)——凡与查询字比对得分不低于阈值 T 的所有 k 字（按替换矩阵计）都算作潜在种子，例如一个三肽字通常能展开出数十个邻字。T 越小，种子越多、敏感性越高而速度越慢。第二步扩展：从每个种子出发沿两个方向作无空位延伸并逐列累加分数，一旦累计分从已经达到的峰值跌落超过阈值 X（X-drop）即停止——分数一旦跌落过深便几乎不可能回升，继续延伸只是浪费。扩展后仍达标的局部比对称为高分片段对(high-scoring segment pair, HSP)，即 BLAST 的基本产出（图 2.6-1）。1997 年的 gapped BLAST 改为先以“同一链上相距有限的两个种子”触发、再作带空位的最终比对，并把逐条序列的统计折算得更精细，在提高速度约三倍的同时降低了假命中率（Altschul et al., 1997）。启发式的代价是明确的：BLAST 不保证找到最优比对，它以放弃理论最优换取数千倍的提速；在绝大多数检索场景中，这一交换是划算的。

BLAST 的第二重贡献是给每个 HSP 配上严格的统计。Karlin 与 Altschul 证明：在随机序列背景下，无空位局部比对的最高得分服从极值分布，由此可以算出得分不低于 S 的命中按期望会出现多少个（Karlin & Altschul, 1990）：

式 (2.6-2) 的解读必须紧扣“期望计数”这一含义：E = 0.001 表示纯靠巧合也能在千次检索中期望撞出一次这么高分的命中，通常可以认真对待；E = 10 则表示这个分数在随机背景下本来就该出现约十次，没有意义。常用的工作阈值：严格的功能推断取 E < 10⁻⁵，一般注释取 E < 10⁻³，粗筛可放宽到 E < 0.05，但必须与一致性、覆盖率合读。还应注意 E 值是相对于检索规模的：同一命中在扩大 100 倍的数据库中 E 值约增大 100 倍——数据库越大，巧合越多，这正是公式的直观内容。

式 (2.6-2) 中得分 S 与参数 λ、K 耦合，使不同打分体系之间、不同规模检索之间的原始得分不可直接比较。为此把 λ 与 K 吸收进得分定义，得到位得分(bit score)：

2.6.3 数据库搜索的实践：程序、参数与规范

BLAST 不是一个程序而是一个家族，五个成员按“查询与数据库各是核酸还是蛋白”划分；核酸库的六个阅读框可以在检索时动态翻译，从而使核酸查询直接接受蛋白层面的证据（表 2.6-1）。

实践中的另一个要点是对低复杂度区域(low-complexity region)的处理。胶原的甘氨酸—脯氨酸重复、卫星序列、富 Alu 区段等不携带多少特异性信息，却极易制造虚高得分，把毫不相干的序列“粘”在一起。BLAST 缺省以 SEG（蛋白）与 DUST（核酸）算法把这类区段遮蔽为 X 或 N 再行检索；解读结果时应留意比对图中被遮蔽的片段——一段看似漂亮的命中若主要落在遮蔽区，其证据价值应大打折扣。

最后是学术写作中的引用与报告规范。BLAST 是一套参数体系：矩阵与空位罚分、字长与邻字阈值、E 值阈值、过滤选项乃至数据库的版本，都会影响结果。报告检索结果时应给出程序名称与版本、数据库及版本（或访问日期）、关键参数设置和命中的具体位置（ accession 与坐标区间），以保证可重复。引用时应当引算法原文——用 BLAST 引 Altschul 等（1990、1997），提到全局或局部最优比对引 Needleman 与 Wunsch（1970）和 Smith 与 Waterman（1981），用 BLOSUM 矩阵引 Henikoff 与 Henikoff（1992）——而不是只给出网站地址。这既是学术溯源的要求，也使方法链条可以被后来的读者完整检查。

2.6.4 比较基因组学：直系、旁系与同线性

当比较的对象从两条序列扩展到两个基因组，必须先把“同源”这个词拆开，因为它对应两种不同的事件。

直系同源的推断在算法上并不容易，因为复制与基因丢失会把“最像”与“对应”搅在一起。最简单的方法是最佳命中互惠(reciprocal best hit, RBH)：若物种 A 的基因 a 在物种 B 的蛋白组中的最佳命中是 b，而 b 反过来在 A 中的最佳命中恰是 a，则推断 a 与 b 为直系同源。RBH 简单快速，是快速注释的常用手段，但一旦基因家族在某一谱系发生过扩张，旁系拷贝可能恰好比真正的直系拷贝更相似而被误判。全基因组规模的方法转而把所有种内、种间相似关系连成一张图再行聚类：COG 数据库以多物种间的最佳命中传递关系聚类原核基因（Tatusov et al., 2001），OrthoMCL 以种内种间同源图配合 MCL 图聚类处理真核基因家族（Li et al., 2003）。它们的共同思想是把直系/旁系判定从孤立的“两两比较”提升为跨物种基因家族的系统关系问题，从而部分抵消家族扩张带来的误导。

在基因次序的层面，比较还有更大尺度的结构信号。相关物种的基因组往往整段整段地保留着基因次序与方向的保守性，称为同线性(synteny)；一段内部次序基本一致的区段称为共线性区块(collinear block)。人与小鼠约 7 500 万年前分化，小鼠基因组计划的首篇比较分析仍能将两个基因组对应地切成 300 余个（约 340 个）大尺度保守同线区块，覆盖基因组的大部分（Waterston et al., 2002）。考察同线性的标准可视化是点阵图(dot plot)：把两个基因组分别放在横轴与纵轴上，在每一对同源位置打一个点（图 2.6-2）。连续的对角线就是同线区块；对角线的中断处（断点）对应分化以来发生的染色体重排——倒位、易位或断裂；方向相反的反向斜线提示一次倒位；主对角线之外的零散灰点则多是旁系同源或重复序列的相似。点阵图因此是全基因组比较的第一张图，读图即读重排历史。

同线性的价值在于它把功能推断从“序列相似”升级为“位置同源”。一个功能未知的基因若恰好落在两个物种互相保守的同线区块中、且其上下游邻居也彼此对应，则其注释的可信度显著高于单凭序列相似；反过来，两个物种同线断裂的分布本身就是一份重排事件的编年记录，可用来重构祖先染色体的结构。比较基因组学由此成为注释流程的放大器：序列层面的 BLAST 找到候选同源基因，结构层面的同线性为它们提供位置证据。

2.6.5 保守性与功能元件的发现

比较基因组学还有一条威力强大的推理路径，出发点是一个群体遗传学事实：不承担功能的序列在中性漂变下以大致恒定的速率自由积累替换，而承担功能的序列因净化选择而进化受限。于是在适当选择的物种距离上，“比中性背景更保守”本身就成了功能的证据。物种距离的选择是一门折中艺术：太近（如人–黑猩猩）中性背景尚未发散，几乎处处“保守”，无从分辨；太远则连功能性同源也已饱和。人与鼠是经典组合——中性位点自分化以来平均已置换约 0.5 次，中性背景已充分发散，而功能元件依然醒目。

在基因调控层面，这一思想发展为系统发育足迹法(phylogenetic footprinting)：对一组物种中同一基因的旁侧非编码区做多重比对，跨物种保守出现的短 motif 多为转录因子结合位点等调控元件。编码区有密码子结构可循，而裸露的非编码背景上，保守性几乎是我们赖以辨认调控信号的唯一通用线索。全基因组尺度的同类比较揭示了大量保守非编码序列(conserved noncoding element, CNE)：它们不编码蛋白、却保守得异乎寻常，且在基因组中富集于发育调控基因附近，其中相当一部分经报告基因实验验证为增强子。

保守性逻辑的极端情形是所谓超保守元件：人、小鼠与大鼠三个物种间长达 200 bp 以上完全一致的片段有 480 余段，多位于调控区与发育基因附近（Bejerano et al., 2004）。如此长度的严格不变远超蛋白编码的需要，提示尚待完全阐明的强功能约束，至今仍是比较基因组学的著名谜题。本节的三层工具——序列层面的同源检索、结构层面的同线性、功能层面的保守性证据——共同把“比较”变成了基因组注释与发现的通用方法；下一章（第 3 章）将把比较的对象从物种之间转向同一物种的个体之间，进入 SNP 与变异的讨论。