基因组科学概论A PRIMER OF GENOME SCIENCE
第四章 · 4.7

4.7 RNA 测序

RNA Sequencing (RNA-Seq)

本节摘要　本节讨论 RNA 测序（RNA-seq）的原理与分析框架：自 poly(A) 富集与 rRNA 去除两类文库策略起，经片段化、双链 cDNA 与链特异性接头连接构建测序文库；继而建立读段随机起源的抽样模型，导出 RPKM/FPKM 与 TPM 定量并辨析其可比性；讨论剪接感知比对、转录本重构与多 isoform 表达估计的似然框架；引入负二项分布处理过度离散的计数差异分析；最后与微阵列系统比较，并归纳主要偏差与陷阱。

上一节说明 SAGE 与 MPSS 确立了"表达量即可数标签"的思想，其通量受制于当时的测序能力。2008 年前后，第二代测序（第 2 章）把每碱基成本压低到使"把整个转录组直接测序"成为常规操作——RNA 测序(RNA sequencing, RNA-seq)在一年之内由三篇标志性工作确立：对哺乳动物转录组的定量测绘（Mortazavi et al., 2008）、对酵母全转录组的链特异性图谱（Nagalakshmi et al., 2008），以及与五种微阵列平台的正面对比（'t Hoen et al., 2008）。原书第 3 版出版于 2009 年，恰是这一技术起飞的时刻。本节按数据产生的顺序展开：文库构建（4.7.1）、定量模型（4.7.2）、比对与转录本重构（4.7.3）、差异分析的计数统计（4.7.4），继而与微阵列作系统比较（4.7.5），并归纳偏差与陷阱（4.7.6）（Gibson & Muse, 2009；Wang et al., 2009）。

4.7.1 文库构建与测序

RNA-seq 的测量对象不是 RNA 本身，而是由 RNA 制成的测序文库。哺乳动物细胞的总 RNA 中，核糖体 RNA 通常占 80%–95%，蛋白编码 mRNA 仅占数个百分点；若不加处理直接建库，绝大部分测序通量将被核糖体序列吞噬。因此建库的第一个决定就是如何处置 rRNA，两条路线各有明确的适用范围。其一是 poly(A) 富集(poly(A) selection)：以寡聚(dT)磁珠捕获带 poly(A) 尾的分子，简便有效，每条 mRNA 分得更多读段，适合以成熟 mRNA 为主的"表达浏览"；代价是丢弃了不带 poly(A) 尾的转录本（多数组蛋白 mRNA、许多非编码 RNA 与微生物转录本），且对降解样本（如福尔马林固定石蜡包埋材料）只余 3′ 端可用，偏置严重。其二是 rRNA 去除(rRNA depletion)：以与 rRNA 互补的寡核苷酸探针结合磁珠清除 rRNA，保留其余全部 RNA，适合全转录组研究（含非编码 RNA）与降解样本；代价是流程更繁琐，且此类商品化试剂在原书出版的 2009 年时尚未标准化。Nagalakshmi 等（2008）在酵母中不经 poly(A) 选择而以随机引物建库，首次得到覆盖非聚腺苷化转录本的链特异性图谱，即依赖这一思路（图 4.7-1）。

富集之后的步骤与普通基因组文库一致：将 RNA（或其反转录产物）随机片段化(fragmentation)为数百万至二百 bp 的片段，与测序仪的读长及簇生成要求匹配；以随机引物反转录并合成双链 cDNA；末端修复、3′ 加 A、连接测序接头，最后以 PCR 富集文库。其中一步值得特别强调——链特异性文库(strand-specific library)：常规双链 cDNA 文库在接头连接后丢失了"转录来自哪条链"的信息，而重叠基因、反义转录本的分析都以方向信息为前提。dUTP 法(dUTP method)的解决方式是：第二链合成时以 dUTP 取代 dTTP，PCR 扩增前用尿嘧啶-DNA 糖基化酶特异性消化含 U 的第二链，只保留第一链进入扩增，从而完整保留转录方向（该法在 2009 年前后逐渐定型）。

RNA-seq 文库构建与测序流程——与微阵列的并行对照 A · RNA-seq：从 RNA 到可数读段 第 1 步：总 RNA 中 rRNA 占绝对多数，直接建库将浪费通量 ① 总 RNA rRNA 占 80%–95% mRNA 仅 2%–5% 灰：rRNA　绿：其他 RNA 第 2 步：poly(A) 富集适合 mRNA 浏览；rRNA 去除适合全转录组与降解样本 ② 富集（二选一） poly(A) 磁珠 rRNA 去除 mRNA 浏览｜全转录组 第 3 步：随机打断为约 200 bp 片段，可在 RNA 或 cDNA 水平进行 ③ 片段化 随机打断 约 200 bp RNA 或 cDNA 水平 第 4 步：随机引物反转录并合成双链 cDNA ④ 双链 cDNA 随机引物反转录 第二链合成 第 5 步：末端修复、加 A、接头连接；dUTP 法实现链特异性；PCR 富集 ⑤ 接头与扩增 末端修复 · 加 A 接头连接 链特异：dUTP 法 PCR 富集文库 第 6 步：输出短读段，表达量以整数计数呈现 ⑥ 测序 输出短读段 整数计数（数字量） B · 微阵列并行对照：同一起点，模拟量输出 RNA 提取与质检 反转录为 cDNA （与 RNA-seq 相同的起点） 荧光标记 Cy3/Cy5 或生物素 染料偏倚需交换设计 与已知探针杂交 只能测阵列上的基因 假设空间封闭 激光扫描 连续强度（模拟量） 有背景下限与饱和上限 分野：RNA-seq 数"可数读段"（整数、动态范围跨数量级、对未知序列开放）；微阵列测"杂交荧光"（连续、受背景与饱和限制、限于已知探针）。 统计后果：前者以负二项计数模型分析（4.7.4），后者以线性模型分析（4.4 节）。链特异性文库（⑤）保留转录方向，是反义转录与重叠基因分析的前提。
图 4.7-1　RNA-seq 文库构建与测序流程，及其与微阵列的并行对照。上排（A）：总 RNA 经 poly(A) 富集或 rRNA 去除后片段化，反转录为双链 cDNA，末端修复、加 A 并连接接头（dUTP 法保留链方向），PCR 富集后上机测序，产出可数的短读段。下排（B）：微阵列路线与 RNA-seq 共享"RNA → cDNA"的起点，分歧发生在测量环节——杂交荧光的连续强度取代了读段计数。悬停各步骤框可查看说明。
方法

选择富集策略的决策清单。目标限于成熟 mRNA 且 RNA 完整（RIN 高）：poly(A) 富集，经济且每条 mRNA 覆盖更深；关注非编码 RNA、微生物转录组、全转录组，或样本已降解（如 FFPE 材料）：rRNA 去除。两条路线产出的定量分母不同——一个是"聚腺苷化 RNA 池"，一个是"除 rRNA 外的全部 RNA"——两类文库的 TPM/RPKM 不可直接混编比较；同一研究内应统一策略，并在报告文库类型（Gibson & Muse, 2009；Wang et al., 2009）。

4.7.2 定量模型：RPKM、FPKM 与 TPM

计数的统计解释建立在一个理想化假设上：读段随机起源——测序是对文库中分子池的均匀随机抽样，一条读段(read)落入某转录本的概率正比于该转录本的摩尔浓度与长度之积。在该假设下，期望读段数同时编码了"多少"（浓度）与"多大"（长度）：等摩尔浓度的两条转录本，长的那条吸引成比例更多的读段。长度是技术混杂而非生物学信号，必须除掉；文库大小（总映射读段数）同样是技术变量，也必须除掉。两个偏倚、两个除法，这就是 RPKM 的全部逻辑（Mortazavi et al., 2008）：

RPKMt = rt / ( ℓt × NM ) (4.7-1) rt：唯一比对到转录本 t 的读段数；ℓt：转录本长度（kb）；NM：该样本总映射读段数（百万计）。数值例：rt = 200，ℓt = 2 kb，NM = 10，则 RPKM = 200 / (2 × 10) = 10。双端测序以"片段"（fragment，一对读段计一次）为计数单位时称 FPKM（Trapnell et al., 2010）。RPKM 可近似理解为"每百万读段中落在该转录本每 kb 上的读段数"。

随机起源假设在三个方面被系统性违背。其一，序列偏好：随机引物对模板局部序列有结合偏好，片段化也存在热点，使读段沿转录本的分布并不均匀。其二，长度偏置：如上所述，原始计数与长度成正比，这是定义式 (4.7-1) 显式校正的部分。其三，GC 偏置：PCR 扩增对高（或低）GC 片段效率不同，造成覆盖随 GC 含量漂移。此外，以 oligo-dT 起始或 RNA 降解会产生 3′ 偏置。这些偏好通常不摧毁大样本下的期望结构，但会改变方差并使局部覆盖失真，在转录本重构与低计数区间须格外小心（图 4.7-2）。

RPKM/FPKM 校正了长度与文库大小，却仍有一个次级缺陷：分母里的总读段数在"除长度"之前就已固定，各转录本长度分布的变化会使一个样本内的 FPKM 之和不守恒，不同样本的 FPKM 之和也不同，比例解释因此漂移。TPM(transcripts per million)把归一化顺序颠倒——先除长度、再对总量归一：

TPMt = ( rt / ℓt ) / Σu ( ru / ℓu ) × 106 (4.7-2) 等价换算：TPMt = RPKMt / Σu RPKMu × 106，即对同一样本内各 RPKM 按比例重新标定。构造上每个样本的 TPM 之和恒等于 106，TPM 因此是"转录本摩尔分数"的直接类比——不同样本间比较占比或分布时，TPM 的可比性优于列和随样本漂移的 RPKM/FPKM（Li & Dewey, 2011）。
A · 等摩尔浓度不等于等计数：长度偏置 转录本 X：1 kb，5 条读段。RPKM = 5/(1×N) 转录本 X：1 kb 读段 5 条 → RPKM = 5/(1×N) 转录本 Y：4 kb，20 条读段。RPKM = 20/(4×N)，与 X 相等 转录本 Y：4 kb 读段 20 条 → RPKM = 20/(4×N) 两者 RPKM 相等：长度的技术混杂已被除法抹平 原始计数 ∝ 浓度 × 长度——跨基因比较必须先除长度 示意：同一文库，总映射读段 N 百万 B · 深度与检出下限（泊松直觉） P(检出 ≥ 1 条读段) = 1 − e^(−λ) λ=0.5：检出概率约 39% 39%0.5 λ=1：检出概率约 63% 63%1 λ=2：检出概率约 86% 86%2 λ=3：检出概率约 95% 95%3 λ=5：检出概率约 99% 99%5 期望读段数 λ = RPKM × ℓ(kb) × N(M) 深度加倍 → λ 加倍 → 可检出的最低 RPKM 减半：检出下限与测序深度、转录本长度均成反比； 且"检出"不等于"定量可靠"——计数区间内相对误差约为 1/√λ，计数个位数时仍很大。
图 4.7-2　读段抽样模型的两个直觉。A 长度偏置：同一文库中等摩尔浓度的两条转录本，4 kb 者吸引的读段是 1 kb 者的四倍；除以长度（式 4.7-1）后两者 RPKM 相等，说明原始计数不可跨基因直接比较。B 深度与检出下限：读段到达近似泊松过程，一个期望计数为 λ 的转录本至少被读到一次的概率为 1 − e^(−λ)；λ 达到 3 时约 95%。检测低丰度转录本的手段是加深测序或（在允许时）富集目标 RNA。悬停图形元素可查看数值。

测序深度(sequencing depth)决定检出下限。给定 RPKM 定量 ρ，某转录本在深度为 N（百万）文库中的期望计数为 λ = ρ × ℓ(kb) × N；泊松抽样下它至少被读到一次的概率为 1 − e^(−λ)。例如 2 kb 转录本在 10 M 文库中定量为 0.05 RPKM 时 λ = 1，检出概率约 63%；深度加倍则同一转录本 λ 加倍，可检出的最低 RPKM 相应减半——检出下限与深度成反比，也与长度成反比（图 4.7-2B）。但"检出"与"定量可靠"是两回事：计数为 λ 的基因，其相对标准差约为 1/√λ，计数个位数时定量极不稳定。

习题 4.7-1

某转录本长 4 kb，在一个总映射读段为 2 000 万的样本中得到 1 600 条唯一比对读段。(1) 计算其 RPKM。(2) 另一样本总映射 1 000 万，其中 2 kb 的转录本得到 200 条读段，其 RPKM 是多少？两个数值可以直接说明谁的"表达水平"更高吗？(3) 若第二个样本加深测序至 2 000 万读段，该转录本读段数近似加倍至 400，RPKM 如何变化？这体现了 RPKM 的什么设计意图？

参考解答

(1) RPKM = 1 600 / (4 × 20) = 20。(2) RPKM = 200 / (2 × 10) = 10；两者都经长度与文库大小双重归一，"每 kb 每百万读段的密度"是第一个的两倍，可以比较。(3) RPKM = 400 / (2 × 20) = 10，保持不变——RPKM 的意图正是让同一转录本在不同深度下的定量可比；但计数 λ 加倍使其相对误差按 1/√λ 下降，即定量更可靠。注意结论的前提是两样本用同一套转录本注释计数，且比较的是同一转录本。

交互演示 · RPKM 与 TPM 换算计算器
修改任一输入框即实时重算：先按式 (4.7-1) 求两条转录本的 RPKM，再按式 (4.7-2) 在样本内重新标定为 TPM。观察默认值（等摩尔的两条转录本）与扰动后的结果。
转录本读段数 r长度 ℓ (kb)RPKMTPM
A——
B——
样本内 Σ RPKM = —；TPM 之和恒为 1 000 000（构造性性质，见式 4.7-2）。默认参数即习题 4.7-1 之外的一个等摩尔例子：A 得 10 RPKM，B 得 10 RPKM，TPM 各为 500 000。
输入须为正数（读段数 r 允许为 0，此时 TPM 无定义）。
习题 4.7-2

同一样本中两转录本的 RPKM 分别为 30 与 10。求各自的 TPM，并回答：为什么任何样本的 TPM 之和都恒等于 106，而 FPKM 之和没有这一性质？据此说明跨样本比较"转录本占比"时为何应优先使用 TPM。

参考解答

TPM₁ = 30/(30+10) × 106 = 750 000；TPM₂ = 10/40 × 106 = 250 000。TPM 先做长度归一得到"每 kb 读段密度"，再把这些密度按比例缩放到总和 106，分母是样本内全部转录本密度之和，因此列和被构造性地固定。FPKM 的分母是总映射读段数，各转录本长度分布与表达构成变化时，FPKM 之和随样本漂移，比例解释不一致。比较占比（摩尔分数式的量）时，TPM 在两个方向上都归一到同一标尺，故跨样本可比性更好；但同一基因在两样本间的差异检验仍应回到原始计数与负二项模型（4.7.4 节）。

4.7.3 比对与转录本重构

测序产出数以千万计的短读段后，计算有两条路。基因组比对把读段放回参考基因组：读段大多完整落在单一外显子内，可与普通比对器对齐；但跨越外显子接头的读段在基因组上不连续——前半段属于一个外显子、后半段属于数十 kb 之外的另一个外显子，中间隔着被剪接删除的内含子。普通比对器会把这类读段判为错配或直接丢弃，因此需要剪接感知比对(splice-aware alignment)：比对器允许读段中间出现大的间隙（N-gap），并把间隙两端约束在符合剪接位点拓扑的位置上——供体位点几乎总是 GT、受体位点几乎总是 AG（另有 GC-AG、AT-AC 等少数类型），这与 2.5 节基因结构模型中供体/受体状态的发射规则一致。以 TopHat（2009 年发表）为代表的工具先用全基因组比对定位"体"读段，再借助候选剪接位点索引处理"接头"读段，从而同时获得表达计数与剪接证据。

没有参考基因组（或参考质量差）时走从头组装(de novo assembly)：不依赖任何外部序列，由读段之间的重叠直接重构转录本。Trinity（Grabherr et al., 2011）是其代表：以高频 k-mer 为种子贪婪延伸出转录本骨架，再把共线的 k-mer 打包为组件，最后在图上用路径区分 isoform 与旁系同源序列。从头组装对发现新物种的转录组不可替代，但高表达基因的深覆盖与低表达基因的稀疏覆盖都给图结构带来困难，重构结果通常须与表达估计、同源证据联合评估。

比对之后是表达估计。困难在于一个基因常有多个 isoform，共享大部分外显子：落在公共外显子上的读段在序列层面无法判定来自哪个 isoform。Cufflinks 与 RSEM 给出的答案是最大似然框架下的概率分配：对每条读段列出全部"兼容"的转录本，按当前的丰度估计把读段按比例分摊给各候选，用分摊结果更新丰度，迭代至收敛——即期望最大化（EM）算法的思想。Cufflinks 在此框架下同时完成转录本组装与定量，能报告未经注释的新 isoform 及其随条件的切换（Trapnell et al., 2010）；RSEM 则显式建模比对概率与测序错误，输出转录本与基因两个层次的定量（Li & Dewey, 2011）。

注记

EM 分配的直观。把歧义读段想象成"只能看出捐给了某几个候选之一的款项"：第一轮按均分或初始比例入账；第二轮起，账面更丰厚的候选按比例分得更多；反复"按比例入账—更新账面"，最终收敛到一组自洽的丰度，使观测到的读段模式出现概率最大。一个微型数字例：10 条读段均与 isoform 甲、乙兼容，另有 5 条读段只与甲兼容；收敛结果是甲 8、乙 2——只与甲兼容的 5 条全归甲，歧义的 10 条按 3:1 的证据倾向分配。EM 不创造信息，只是把"部分信息"按自洽原则用尽。

4.7.4 差异表达的计数统计

把读段指派到基因（或转录本）后，得到"基因 × 样本"的整数计数矩阵——形式上与微阵列矩阵相同，但取值是计数而非强度，这决定了统计模型必须更换。计数的抽样噪声天然由泊松分布刻画：早期研究表明技术重复的计数散布确实接近泊松预期；但生物学重复之间存在真实的个体差异——同一基因在不同个体里的实际表达量本就波动——使观测方差系统性超过泊松预测，即过度离散(overdispersion)（图 4.7-3）。对过度离散的计数强行使用泊松检验或 t 检验，会低估方差、夸大显著性，重复越少失真越重。

Var(Kgj) = μgj + φg · μgj2 (4.7-3) Kgj：基因 g 在样本 j 的计数；μgj：与样本条件相应的期望计数；φg：基因 g 的离散度，刻画超出泊松的额外散布。φ 趋于 0 时退化为泊松。该式是负二项分布(negative binomial distribution)的方差结构，等价于"基因真实表达水平在样本间服从伽马分布波动、给定水平后计数为泊松抽样"的混合模型——φμ² 项吸收的正是生物学变异（Robinson et al., 2010；Anders & Huber, 2010）。
计数方差随均值：泊松（技术重复）与负二项（生物学重复） 1101001000 基因平均计数 μ（对数刻度） 10⁰10¹10²10³10⁴10⁵ 跨重复方差（对数刻度） 生物学重复的点：围绕负二项曲线 Var = μ + 0.1μ² 散布 技术重复的点：贴近泊松线 Var = μ 负二项：生物学重复 泊松：技术重复 Var = μ + φμ²（φ ≈ 0.1） 高表达端偏离泊松可达约百倍 示意数据：生物学重复围绕负二项曲线散布，技术重复贴近泊松线；高表达基因的绝对偏离最大，这正是需要离散度参数的原因。
图 4.7-3　均值-方差关系中的过度离散（示意数据，双对数坐标）。技术重复只含泊松抽样噪声，散布贴近 Var = μ（灰虚线）；生物学重复叠加了个体间真实的表达波动，散布围绕负二项曲线 Var = μ + φμ²（绿实线，此处取 φ ≈ 0.1）。在 μ = 1000 处负二项方差约为泊松的一百倍——用泊松模型检验生物学重复，会按同样倍数高估显著性。悬停点集可查看说明。

式 (4.7-3) 中的离散度 φ 逐基因而异，而每个基因通常只有三至六个重复——直接由样本估计 φ 极不稳定。DESeq 与 edgeR 的共同对策是跨基因"借用信息"：假定离散度随表达水平平滑变化，先以全部基因拟合一条趋势，再把个别基因的估计向趋势收缩（经验贝叶斯思想，与微阵列时代 limma 收缩方差的做法一脉相承）；edgeR 进一步分出共同、趋势化、逐基因三层估计供不同重复规模选用。标准化方面，除文库大小外还须处理文库构成偏置：少数高表达基因的大幅上调会按比例吸收更多读段，压低其余所有基因的相对计数，即便它们真实不变。DESeq 的中位数比率法为此设计：

sj = mediang ( kgj / ( ∏j′ kgj′ )1/m ) (4.7-4) kgj：基因 g 在样本 j 的计数；分母为该基因在 m 个样本中的几何平均，充当"虚拟参考"；对每个样本取全体基因比值的中位数即尺寸因子 sj，用于把各样本拉到同一标尺。以中位数聚合，使少数强差异基因或极端计数不至牵动整体；edgeR 的 TMM（截尾 M 值平均）以截尾均值达到同一目标。归一化校正的是技术标尺，不应吞掉生物学差异。
习题 4.7-3

某基因在对照组三个生物学重复中的计数为 70、130、100。(1) 估计均值与样本方差。(2) 泊松模型预测的方差是多少？用式 (4.7-3) 反推离散度 φ 的量级。(3) 若无视过度离散、按泊松方差做两组比较，结论会发生什么系统性偏差？

参考解答

(1) 均值 = 100；样本方差 = [(70−100)² + (130−100)² + (100−100)²]/(3−1) = 900。(2) 泊松预测 Var = μ = 100，实测 900，超出约八倍；由 900 ≈ 100 + φ·100² 得 φ ≈ 0.08。(3) 泊松把标准误低估约三倍（√9），检验统计量随之虚高，p 值系统性偏小——重复越少、基因越高表达，假阳性越严重。正确做法是采用负二项模型，并借助跨基因趋势收缩获得稳定的 φ 估计（4.7.4 节）。

4.7.5 与微阵列的系统比较

三篇早期工作提供了比较的实证基础。Mortazavi 等（2008）在小鼠组织与人类细胞中给出 RPKM 定量框架，显示 RNA-seq 动态范围跨越数个数量级、并能以接头读段直接观察剪接；Nagalakshmi 等（2008）在酵母中以链特异性文库修正了大量注释错误并揭示非聚腺苷化转录本；'t Hoen 等（2008）将 RNA-seq 与五种微阵列平台同台对比，结论是 RNA-seq 在重复性、分辨率与实验室间可移植性上均占优、动态范围更宽、低丰度检出更好，但当时单位成本与计算门槛显著更高。综合这些证据与原书的判断，整理为表 4.7-1。

表 4.7-1微阵列与 RNA 测序的系统比较（2009 年时点；据 Mortazavi et al., 2008；'t Hoen et al., 2008；Wang et al., 2009；Gibson & Muse, 2009 整理）
比较维度微阵列RNA 测序
测量原理与固定探针杂交的荧光强度：连续模拟量，受背景下限与饱和上限约束，动态范围约 102读段计数：整数数字量，无饱和上限，动态范围可跨 104 以上
先验序列必须预先设计探针，假设空间封闭有参考时可注释已知并对新转录本开放；无参考可从头组装
新转录本与融合基因探针之外的目标不可见接头读段可定位新剪接与融合转录本；读段可携带编辑或突变信息
等位特异性表达探针杂交受 SNP 影响产生"基因型假差异"可按 SNP 相位把读段归属到等位基因（需足够覆盖）
剪接定量需专门的外显子/接头探针阵列，交叉杂交干扰明显接头读段直接计数；isoform 级定量属概率推断（4.7.3 节）
统计模型连续强度：线性模型（limma）+ 经验贝叶斯（见 4.4 节）整数计数：负二项（DESeq/edgeR）+ 尺寸因子归一化
成本结构仪器成熟、单价低、数据分析负担轻测序深度即成本变量；计算与存储负担显著更高
成熟度（2009 年）MIAME/GEO 体系完备，跨实验室可比性久经检验流程未标准化；'t Hoen 等（2008）显示其实验室间一致性已领先各芯片平台
案例

RNA-seq 与微阵列的交界年代（2008–2010）：互补而非立即替代。这一时期装备完善的实验室常并行运转两套平台：大样本临床队列、既定基因 panel 与历史数据的衔接用芯片——便宜、标准化、结果与十余年积累直接可比；发现导向的研究（新转录本、可变剪接、融合基因）、无参考物种与宽动态范围需求用 RNA-seq。't Hoen 等（2008）的对比研究正是这一格局的注脚：新平台在多项指标上领先，但作者同时报告了成本、通量分配与数据分析门槛的现实约束。此后数年，随着去 rRNA 试剂、UMI 与工具链（如后续版本的比对器与 DESeq2）成熟，测序逐渐成为默认平台——这一后见之明超出了原书 2009 年的时点，彼时作者给出的判断是务实的"互补"（Gibson & Muse, 2009）。

4.7.6 偏差与陷阱

RNA-seq 的常见陷阱多数发生在计数矩阵生成之前。第一，rRNA 残留：去 rRNA 试剂失效或操作失当，可使两成到八成读段落在核糖体序列上，有效深度骤降；惯例的质控是把读段比对结果按注释类别（rRNA、外显子、内含子、基因间区）分层统计。第二，多映射读段(multi-mapping reads)：旁系同源基因家族与假基因使部分读段无法唯一归属。丢弃多映射读段简单保守，但会系统性压低这些家族成员的表观表达；按候选位置概率分摊能保留信息，却引入"成员间均匀表达"等难以检验的假设。两种策略必须在研究内一致并如实报告。第三，覆盖不均一：序列偏好、GC 偏置与 3′ 偏置使读段沿转录本分布疏密不均（4.7.2 节），直接后果是 isoform 定量与短转录本的置信区间变宽，绘图检查沿转录本的覆盖曲线应成为例行步骤。

第四，批次与文库复杂度(library complexity)。RNA-seq 同样存在批次效应——不同测序通道、流动槽、建库日期的差异会沉淀进数据，分组必须跨批次随机平衡，与 4.3 节的设计原则一致。文库复杂度指文库中独立分子（而非 PCR 拷贝）的数目：起始量过低或扩增过度时复杂度下降，重复读段占比上升，有效深度缩水；独特分子标识符（UMI，原书出版后普及）通过给每个原始分子加标签，使重复拷贝可被识别与合并，从根源上缓解 PCR 偏置。

警示

覆盖深度不等于表达量：未归一化的原始计数不可直接比较。样本甲深度 40 M，基因 X 得 400 条读段；样本乙深度 10 M，基因 Y 得 200 条读段。直接比较 400 对 200 会得出"X 更高表达"，但 X 与 Y 长度不同、两文库总读段数相差四倍——换成 TPM 或经尺寸因子校正的比较完全可能反转。同理，同一样本内"基因 X 读段多于基因 Y"不能推出"X 表达更高"，因为 X 可能只是更长。原始计数只在与模型相连（负二项 + 尺寸因子）或经式 (4.7-1)/(4.7-2) 归一化之后才有解释力；热图、箱线图等任何跨样本图形都应以归一化量绘制。

总结本节：RNA-seq 把 SAGE 的"计数"思想与第二代测序的通量结合——文库策略决定测什么（poly(A) 或全转录组、是否保留链方向），抽样模型决定怎么算（RPKM/FPKM/TPM），比对与似然估计决定把读段放到哪（剪接感知比对、EM 定量），负二项模型决定怎么比（离散度收缩 + 尺寸因子）。它与微阵列的交接及其在具体研究中的应用，见下一节。

关键术语

RNA 测序 (RNA sequencing, RNA-seq)
对转录组片段文库进行大规模平行测序、以读段计数测量表达的技术。
poly(A) 富集 (poly(A) selection)
以寡聚(dT)磁珠捕获带 poly(A) 尾分子的 mRNA 富集策略。
rRNA 去除 (rRNA depletion)
以互补探针清除核糖体 RNA、保留其余 RNA 的全转录组建库策略。
链特异性文库 (strand-specific library)
保留转录方向信息的 RNA-seq 文库，如以 dUTP 法构建。
读段 (read)
测序仪自单个克隆化片段读出的短序列，RNA-seq 的基本计数单位。
RPKM/FPKM (reads/fragments per kilobase per million)
按转录本长度与文库大小双重归一化的相对表达定量。
TPM (transcripts per million)
先长度归一、再总量归一的摩尔分数式定量，样本内之和恒为 106。
剪接感知比对 (splice-aware alignment)
允许读段跨外显子接头、以 GT-AG 等剪接拓扑约束间隙的比对策略。
从头组装 (de novo assembly)
不依赖参考基因组、由读段重叠直接重构转录本序列（如 Trinity）。
过度离散 (overdispersion)
生物学重复的计数方差超过泊松预测（Var > μ）的现象。
负二项分布 (negative binomial distribution)
方差为 μ + φμ² 的计数分布，RNA-seq 差异分析的标准模型。
多映射读段 (multi-mapping reads)
可比对到多个基因组位置的读段，归属须丢弃或按概率分摊。

参考文献与延伸阅读

  1. Gibson G, Muse SV. 2009. A Primer of Genome Science, 3rd ed. Sunderland (MA): Sinauer Associates. (Chapter 4)
  2. Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. 2008. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nature Methods 5: 621–628.
  3. Nagalakshmi U, Wang Z, Waern K, Shou C, Raha D, Gerstein M, Snyder M. 2008. The transcriptional landscape of the yeast genome defined by RNA sequencing. Science 320: 1344–1349.
  4. Wang Z, Gerstein M, Snyder M. 2009. RNA-Seq: a revolutionary tool for transcriptomics. Nature Reviews Genetics 10: 57–63.
  5. 't Hoen PAC, Ariyurek Y, Thygesen HH, et al. 2008. Deep sequencing-based expression analysis shows major advances in robustness, resolution and inter-lab portability over five microarray platforms. Nucleic Acids Research 36: e141.
  6. Trapnell C, Williams BA, Pertea G, et al. 2010. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nature Biotechnology 28: 511–515.
  7. Robinson MD, McCarthy DJ, Smyth GK. 2010. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26: 139–140.
  8. Anders S, Huber W. 2010. Differential expression analysis for sequence count data. Genome Biology 11: R106.
  9. Li B, Dewey CN. 2011. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics 12: 323.
  10. Grabherr MG, Haas BJ, Yassour M, et al. 2011. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nature Biotechnology 29: 644–652.