2.5 基因组注释

经过 2.3 节的拼装与 2.4 节的完成工作，测序仪产出的海量读段变成了以碱基为单位的连续序列——对人类基因组而言，约 31 亿个字符的文本。然而对计算机而言，这仍是一串没有句读、没有标题的字符串：哪些区段是基因、基因从何处开始到何处结束、其产物行使什么功能，序列本身并未言明。基因组注释(genome annotation)就是在基因组序列上系统地标定生物学元件并解释其生物学意义的过程（Stein, 2001）。注释之于基因组，犹如索引与笺注之于一部无人读过的巨著：它把"字符流"变成可检索、可比较的"基因目录"，是此后一切变异、表达与调控分析的基础。

注释也不是一次性的终点。草图阶段的注释以计算预测为主，此后随转录本、蛋白质与功能实验证据的累积不断修订——本节 2.5.3 将以人类基因数目的变迁为例说明这一点。本节按注释流程的自然顺序展开：先界定注释的两重含义（2.5.1），再讨论从头预测的统计原理（2.5.2）与多源证据的整合（2.5.3），继而把视野扩展到非编码元件（2.5.4）与重复序列（2.5.5），最后讨论功能标注体系及其证据可信层级（2.5.6）。

2.5.1 注释的两重含义：结构与功能

结构注释(structural annotation)回答"元件在哪里、边界如何"：蛋白编码基因的坐标（外显子—内含子结构、编码区起止）、非编码 RNA 基因、调控区域、重复序列与假基因的位置。功能注释(functional annotation)回答"编码什么、做什么"：为基因产物指定名称、结构域组成、分子功能、所参与的生物过程与细胞定位。二者的输出形式不同：前者是一组坐标（现代数据库常以 GFF/GTF 等格式记录），后者是基因与受控词汇之间的关联表。实践中二者相互依赖——功能注释以结构注释的结果为对象，而结构注释又常以功能信息为证据（如与已知蛋白的比对可用于锚定外显子边界）。

基因模型的复杂程度决定了结构注释的难度，而原核与真核恰好位于两个极端。

图 2.5-1 以一个典型的真核蛋白编码基因座为例，列出结构注释需要确定的全部坐标：启动子元件（TATA 框、GC 框等，多数位于转录起始位点上游数百碱基内，增强子则可远达数十 kb 之外）、转录起始位点(transcription start site, TSS)、5′ 非翻译区（5′UTR）、旁侧符合 Kozak 共识的起始密码子 ATG、外显子—内含子边界（供体与受体剪接位点）、终止密码子、3′UTR 内的加尾信号 AATAAA 与 polyA 位点。转录后内含子被剪接体去除，外显子连同两端 UTR 组装为成熟 mRNA（图 2.5-1B）。

需要强调，图中的"典型"结构只是统计意义上的众数：首、末外显子既含 UTR 又含编码区；部分基因（如多数组蛋白基因）没有内含子；另一些内含子位于 UTR 之内。结构注释的输出因此总带有不确定性，坐标应当附有来源与置信度——这也是主流数据库（如 Ensembl/GENCODE）把基因集按证据强度分级的原因。

2.5.2 从头预测：统计特征与隐马尔可夫模型

不依赖任何转录本或同源蛋白证据、仅凭序列本身的统计规律识别基因，称为从头预测(ab initio prediction)。其信息来源可归为两类：一是信号——短小而位置关键的离散 motif，如剪接位点、起始与终止密码子、分支点、Kozak 共识与加尾信号；二是含量——编码区特有的统计签名。前者定位边界，后者区分"像编码区"与"不像编码区"的区段。

编码区的统计签名主要有三种表现形式。其一为密码子三周期性(period-3 periodicity)：编码区每三个碱基一个密码子，碱基组成随其在密码子中的位置（第 1、2、3 位）系统变化，对序列做频谱分析会在 1/3 频率处出现峰值，而内含子与基因间区无此规律。其二为密码子使用偏倚(codon usage bias)：同义密码子的使用频率在各物种中显著偏离均匀且物种间各异；把一段序列按某个阅读框分成三联体后，其频率分布偏离均匀期望的程度（如 Karlin 等提出的按密码子位置分组的碱基组成签名统计量 Ks 一类偏差度量）随该区段是否编码、以及阅读框是否正确而显著不同。其三为六聚体统计(hexamer statistics)：把窗口内全部相邻六碱基在编码与非编码训练集中的频率比取对数累加，作为该窗口的编码得分——六聚体比单碱基或二核苷酸捕捉到更强的相邻依赖，比单个密码子更能区分弱编码区与背景，是 GENSCAN 等程序中编码状态的核心发射得分。

剪接信号则提供了另一类强约束。供体位点（5′ 端，外显子—内含子交界）的一致性序列可写作外显子末尾的 CAG 接内含子开头的 gtaagt：内含子前两位 GT 在经典内含子中几乎不变；受体位点（3′ 端）则是内含子末端 AG 之前一段 10–20 个嘧啶组成的富嘧啶区，再上游 18–40 nt 处有分支点 A（YNYURAY 一类共识）。需要注意"几乎"二字：绝大多数内含子为 GT–AG 型，另有极少数为 AT–AC 型（由 U12 型剪接体识别）及 GC–AG 等罕见形式——例外虽然稀少，却提示剪接信号应当以权重矩阵（概率）而非硬规则来建模。

把信号与含量统一起来的主流框架是隐马尔可夫模型(hidden Markov model, HMM)：把基因结构拆成一串隐藏状态（基因间区 → 起始外显子 → 供体位点 → 内含子 → 受体位点 → 内部外显子 → … → 终末外显子 → 加尾），每个状态以特定概率分布"发射"碱基，状态之间以转移概率相连；观测到序列后，用 Viterbi 算法找出最能解释该序列的状态路径，路径即预测的基因结构（图 2.5-2）。

代表性工具方面，GENSCAN（Burge & Karlin, 1997）以广义 HMM 首次在脊椎动物基因组上达到实用精度：核苷酸水平的敏感度与特异度约 0.9，外显子完全正确的比例约四分之三（其原始测试集）；GeneMark 系列以马尔可夫链编码测度见长，其自训练变体可在毫无注释的新基因组上迭代估计参数；AUGUSTUS（2003 年发表）延续了广义 HMM 框架，并在 2006 年以后逐步把 RNA-seq 与同源证据纳入同一概率模型，成为多证据整合时代的代表（见 2.5.3）。

如何评价预测好坏，需要一对互补指标：敏感度(sensitivity, Sn)与特异度(specificity, Sp)（后者在预测语境下亦称精确率）。

使用这对指标时有两点必须注意。其一，层次不同结论迥异：核苷子水平的数字容易偏高（边界大体正确即计入），而"整条基因的全部外显子与边界完全正确"的比例通常远低——对人类基因组，单靠从头预测能完整复原的基因至今只有一半上下的量级（Brent, 2008）。其二，特异度存在基数陷阱：真实外显子仅占基因组约 1%–2%，把全部序列都判为非基因也能获得 Sp≈0.99 的核苷酸水平特异度，因此必须看外显子与基因层次、并结合假发现率一起报告。Sn 与 Sp 之间存在此消彼长的权衡：放宽判定阈值可提高敏感度，代价是更多假阳性，习题 2.5-1 即演练这一权衡。

2.5.3 证据整合：转录本、同源蛋白与共识注释

从头预测的精度上限受制于编码信号本身的微弱，实际注释总是多源证据的融合。第一类证据是转录本。表达序列标签(expressed sequence tag, EST)是从某组织 cDNA 文库随机挑取克隆、单端测序 400–700 bp 所得的片段；dbEST 累计收录的 EST 以百万计，构成最早的规模化转录证据，全长 cDNA（FLcDNA）则更进一步。把转录本比对回基因组：比对块即外显子、比对块之间的间隙即内含子，剪接位点由此被直接"锚定"——只要转录本跨越某个内含子，该内含子两端的边界就不再依赖统计推断。

第二类证据是同源蛋白。别的物种中已获实验证实的基因，其编码序列在进化中被相当程度地保守，于是可把已知蛋白与基因组做蛋白对核酸的比对：tBLASTn 一类的搜索把基因组动态翻译成六个阅读框再与蛋白查询比对，能跨物种锚定外显子；Exonerate 一类的模型化比对程序更进一步，把比对打分与剪接信号语法（如内含子必须符合供体—受体结构）联合建模，直接输出符合基因语法的外显子链。比对与同源关系的原理详见 2.6 节。（RNA-seq 括注：原书出版于 2009 年，彼时大规模平行测序刚开始改造转录组学；此后 RNA-seq 拼接转录本迅速取代 EST/cDNA 成为主要转录证据，其与基因预测器的融合正是 AUGUSTUS 等工具的现代形态。）

把这些证据组织起来的是共识注释流水线。以 Ensembl 的自动注释系统为例（Curwen et al., 2004）：先对基因组做重复掩蔽（见 2.5.5）；并行运行多种从头预测器与蛋白同源搜索、转录本比对，各自产生候选基因模型；随后按证据优先级合并——全长 cDNA 支持的模型优先，其次蛋白同源，再次纯预测——冲突处消解后输出共识基因集，并随证据更新定期重注释。整个过程强调两点：证据等级与可追溯性，即每一条注释都能回答"依据是什么"。

2.5.4 非编码元件：tRNA、rRNA、其他 ncRNA 与假基因

蛋白编码基因之外，结构注释还须清点非编码 RNA。tRNA 基因数量不多却高度规律：其三叶草二级结构与关键位点的序列在进化中严格保守。tRNAscan-SE（Lowe & Eddy, 1997）把快速初筛与协方差模型打分相结合，既保持速度又把随机序列上的假阳性压到极低水平，长期是该领域的标准工具；人类基因组可检出数百个正常 tRNA 基因，另有数量更多的 tRNA 衍生假基因。rRNA 则以串联簇组织：18S–5.8S–28S rRNA 作为同一个前体转录单位，以数百拷贝串联排列于近端着丝粒染色体的核仁组织区（NOR），5S rRNA 独立成簇；rRNA 序列高度保守，同源比对即可定位，但其高拷贝串联结构在拼装中常发生坍缩（与 2.3 节重复序列问题同源）。

miRNA 的识别依赖两条线索的交集：前体发夹结构（约 60–110 nt 的茎环）与成熟序列尤其是"种子"区在近缘物种间的保守性。更一般地，snRNA、snoRNA 等结构化非编码 RNA 由 Rfam 统一收录——Rfam 为每个家族构建协方差模型，同时建模序列与二级结构的保守性（对照 Pfam 只建模序列的 profile HMM，见 2.5.6），是结构 ncRNA 注释的参照库。（lncRNA 括注：长链非编码 RNA 的系统注释是原书出版之后才展开的领域，多依赖 RNA-seq 覆盖与染色质标记等证据，此处从略。）

2.5.5 重复序列注释：类别、占比与掩蔽

注释流程通常以重复序列的处理为第一步，原因有二。其一，占比：人类基因组约一半序列可归为重复（IHGSC, 2001）；其二，干扰：未加处理的重复会在预测层面制造大量假开放阅读框与假外显子，在拼装层面造成坍缩与错连（2.3 节已详述拼装一侧，此处补注释一侧）。重复的主体是转座子，按转座机制分为两大类（图 2.5-3）。

第一类为逆转座子（Class I），走"复制—粘贴"路线：以 RNA 为中间体，元件先转录、再逆转录为 DNA 插入新位点，原件保留，拷贝数净增。其中 LINE-1（L1）是自主元件：全长约 6.1 kb，编码 RNA 结合蛋白（ORF1）与兼具逆转录酶和内切酶活性的 ORF2，人类基因组中约 90 万拷贝（绝大多数 5′ 截断，具完整转座能力的全长拷贝仅百余个的量级）；SINE 则是非自主的短元件，灵长类特有的 Alu（约 300 bp，由 7SL RNA 衍生）自身不编码蛋白，借用 LINE-1 的机器实现扩增，拷贝约 110 万。LTR 逆转座子（如内源性逆转病毒 HERV 家族）同属此类的化石遗迹。第二类为 DNA 转座子（Class II），走"剪切—粘贴"路线：转座酶识别末端反向重复，把元件切离后插入新位点，拷贝数基本不变；此类在人类基因组中约占 3% 且无活跃家族。

注：散在重复（转座子来源）合计约 45%，加上卫星、简单重复与区段重复，人类基因组中可归为重复的序列约占一半（IHGSC, 2001；原书出版于 2009 年，此处沿用其引用的量级数据）。

注释重复的主力工具是 RepeatMasker，其策略是"已知库匹配"：以 Repbase（Jurka, 2000）等重复文库收录的各类元件一致性序列为数据库，用 cross_match 一类敏感搜索把基因组序列与库中元件逐一比对，命中的片段按库中的分类标注；输出既包括全基因组的重复构成谱（表 2.5-1 一类统计的来源），也包括"掩蔽"后的序列——硬掩蔽把命中区段替换为 N，软掩蔽仅改为小写字母而保留序列信息，下游基因预测通常在软掩蔽序列上进行，以便保留的非重复区仍可被完整搜索。对新测序物种，库中未必收录其特异性重复，需先以 RepeatModeler 一类工具从头构建重复文库再行掩蔽（原书出版后发展起来的常规做法）。

重复对基因注释的干扰具体而实在：Alu 与 LINE 片段中随处可见无意义的开放阅读框，未掩蔽时会被预测器当作外显子，逆转录酶片段更可能被完整地"注释"成基因；同时，重复也是拼装断裂的主因（2.3 节），scaffold 的断点多落在重复内部，因此重复注释、拼装与基因注释三者实际上交替迭代、互相校准。

2.5.6 功能标注体系：GO、结构域与证据可信层级

结构注释给出基因目录之后，功能注释回答"它们做什么"。最普遍的路径是同源转移：将新基因与已获实验表征的直系同源基因比对，把后者的功能描述与术语转移到新基因（比对方法与直系/旁系之辨见 2.6 节）。第二条路径是结构域分析：蛋白质常由保守的结构域模块组成，Pfam 等数据库以概型隐马尔可夫模型(profile HMM)为每个家族建模（注意与基因预测 HMM 的分工不同：此处 HMM 用于给"序列—家族"的比对打分），InterPro 则整合 Pfam、PROSITE、SMART 等多个资源给出统一的家族、结构域与位点注释。"若干锌指＋一个激酶域"这样的组合往往直接提示蛋白的功能类别。此外，酶学委员会的 EC 编号以四段层级为酶反应赋予唯一代码，KEGG Orthology（KO）则在通路语境中界定直系同源群——各体系侧重不同，注释中常交叉引用。

为使不同数据库、不同物种的功能描述可以互相比较，基因本体(Gene Ontology, GO)项目建立了受控词汇体系，把功能描述组织为三个互相独立的本体（Ashburner et al., 2000）：分子功能（molecular function，产物做什么，如"ATP 结合"）、生物过程（biological process，参与什么事件，如"DNA 修复"）与细胞组分（cellular component，位于何处，如"核糖体"）。

每条 GO 注释还携带证据代码(evidence code)，标明知识的来源：IDA（直接测定）、IMP（突变表型）、IGI（遗传互作）等表示实验证据；ISS（序列相似）与 IEA（inferred from electronic annotation，电子注释）表示计算推断。其中 IEA 由流水线自动转移、未经人工审核，历史上曾占 GO 注释的绝大多数；其可信层级低于任何实验代码，且可能源自另一条 IEA（错误会沿传递链放大）。多数严格的富集分析与注释比较默认剔除 IEA，使用时务必查看注释的来源（with/from 字段）并追溯原始证据。

至此，本节完成了从"序列"到"目录"再到"功能"的两级跳跃：结构注释给出元件坐标，功能注释赋予生物学含义，而二者都建立在证据与可信度的显式表达之上。下一节（2.6 节）将讨论注释中使用频率最高的计算工具——序列比对，并由两两比对延伸到全基因组的比较基因组学。