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

4.7 RNA 测序

RNA Sequencing (RNA-Seq)

本节摘要　本节以短读长 cDNA 测序为主，介绍 RNA 测序（RNA-seq）的文库构建与分析框架。先说明 poly(A) 富集、rRNA 去除与链特异性文库各自保留的信息，再从理想抽样模型导出 RPKM/FPKM 和 TPM，讨论相对定量的可比性；随后介绍剪接感知比对、转录本重构与概率定量，以及用负二项模型处理过度离散的差异表达分析。最后结合早期研究比较 RNA-seq 与微阵列，并说明文库偏差、歧义比对和批次效应对结论的限制。

上一节介绍了 SAGE 与 MPSS 如何用标签计数估计转录本丰度，其通量受到当时测序能力的限制。2008 年前后，第二代测序（第 2 章）降低了大规模序列计数的成本，推动了RNA 测序(RNA sequencing, RNA-seq)的发展。Mortazavi 等（2008）定量测绘了小鼠组织转录组，Nagalakshmi 等（2008）用酵母 poly(A) RNA 改进了转录区域的注释；同期，'t Hoen 等（2008）比较了深度测序的 3′ 标签 DGE 与五种微阵列平台。后者属于标签计数技术，应与本节的随机片段 RNA-seq 区分。原书第 3 版出版于 2009 年，本节以这一时期的研究为背景，并补充 2010–2011 年的代表性分析方法（Gibson & Muse, 2009；Wang et al., 2009）。

4.7.1 文库构建与测序

本节介绍的短读长 RNA-seq 先将 RNA 转换为 cDNA 文库，再通过测序推断转录本丰度。哺乳动物细胞总 RNA 中，核糖体 RNA 通常约占 80%–95%，蛋白编码 mRNA 仅占数个百分点，具体比例随样本而变。若研究目标是 mRNA 或其他非 rRNA 转录本，直接对总 RNA 建库会使大量读段来自 rRNA，因此需要选择合适的富集策略。poly(A) 富集(poly(A) selection)以寡聚(dT)磁珠捕获带 poly(A) 尾的分子，适合研究多数成熟真核 mRNA，也会保留带 poly(A) 尾的非编码 RNA。它通常不能有效捕获不带 poly(A) 尾的 RNA，例如多数复制依赖型组蛋白 mRNA；RNA 降解还可能造成靠近 3′ 端的片段富集，使覆盖不均。

rRNA 去除(rRNA depletion)通过互补探针结合磁珠分离或酶促处理等方式降低 rRNA 含量，可扩大对非聚腺苷化 RNA 的取样范围，也常用于降解样本。但它不能保留所有 RNA：探针的物种适配性、去除效率、片段大小筛选和后续建库反应都会影响可测对象。Nagalakshmi 等（2008）的酵母实验使用的是先分离的 poly(A) RNA，再分别以 oligo(dT) 或随机六聚体引物合成双链 cDNA（GEO：GSE11209）。随机引物并不意味着省略 poly(A) 富集；该实验也不是常规意义上的链特异性建库，部分转录方向由含 poly(A) 尾的末端读段推断。

文库的后续步骤取决于具体方案。图 4.7-1 以常见流程为例：将 RNA 片段化(fragmentation)至约数百 nt，以随机引物反转录并合成双链 cDNA，随后进行末端修复、3′ 加 A、接头连接和 PCR 富集；另一些方案先合成 cDNA 再片段化，图中约 200 bp 是 cDNA 插入片段长度的示意值。片段长度需要与文库制备和测序方案相匹配。

链特异性文库(strand-specific library)保留原始 RNA 的转录方向，有助于区分反义转录本和相互重叠的基因；普通双链 cDNA 文库往往丢失这一信息。dUTP 法(dUTP method)在第二链合成时以 dUTP 取代 dTTP，再选择性处理含 U 的第二链，使后续扩增保留第一链的方向信息。尿嘧啶-DNA 糖基化酶移除 U 碱基后形成无碱基位点；选择性去除或失活第二链还需要后续处理。该类方法在 2009 年已有报道，实验中仍需检查链特异性是否达到预期（Parkhomchuk et al., 2009）。

RNA-seq 文库构建与测序流程：微阵列对照 A · RNA-seq：从 RNA 到可数读段 第 1 步：总 RNA 中 rRNA 通常占多数，非 rRNA 目标的有效覆盖受其比例影响 ① 总 RNA rRNA 占 80%–95% mRNA 仅 2%–5% 灰：rRNA　绿：其他 RNA 第 2 步：poly(A) 富集适用于多数成熟 mRNA；rRNA 去除扩大 RNA 取样范围，仍须评估样本与方案的适配性 ② 富集（二选一） poly(A) 磁珠 rRNA 去除 mRNA 为主｜范围较广 第 3 步：片段化可在 RNA 或 cDNA 水平进行；图示约 200 bp 指 cDNA 插入片段 ③ 片段化 随机打断 约 200 bp RNA 或 cDNA 水平 第 4 步：随机引物反转录并合成双链 cDNA ④ 双链 cDNA 随机引物反转录 第二链合成 第 5 步：末端修复、加 A、接头连接及 PCR 富集；dUTP 方案通过选择性处理第二链保留方向 ⑤ 接头与扩增 末端修复 · 加 A 接头连接 链特异：dUTP 法 PCR 富集文库 第 6 步：输出短读段，经比对和分配获得计数，再用于表达定量 ⑥ 测序 输出短读段 读段计数 B · 微阵列对照：RNA 起始，荧光强度输出 RNA 提取与质检 反转录为 cDNA （与 RNA-seq 相同的起点） 靶标标记 Cy3/Cy5 或生物素 双通道需评估染料偏倚 与已知探针杂交 检测依赖探针覆盖 未覆盖序列难以检测 激光扫描 连续强度（模拟量） 有背景下限与饱和上限 RNA-seq 记录读段计数，可覆盖未预设序列；微阵列测杂交荧光，受背景、饱和与探针设计限制。 计数常用负二项模型（4.7.4），微阵列常用线性模型（4.4）。链特异性文库保留转录方向，帮助区分反义与重叠转录。
图 4.7-1　短读长 RNA-seq 的示例流程及微阵列对照。上排（A）：RNA 经 poly(A) 富集或 rRNA 去除后，用于片段化、cDNA 合成、接头连接和扩增；具体步骤顺序及是否保留链方向由方案决定。图中的 RNA 比例和约 200 bp 插入片段长度仅作示意。下排（B）：微阵列同样从 RNA 提取和质检开始，随后制备带标记的杂交靶标，以荧光强度测量丰度。RNA-seq 则记录读段，经过比对和分配形成计数。悬停各步骤框可查看说明。
方法

根据目标 RNA 和样本质量选择文库。若主要关注成熟 mRNA，且 RNA 完整性较好，可选择 poly(A) 富集，以提高此类转录本的有效覆盖。若关注非聚腺苷化 RNA、微生物转录组或降解材料，可评估适配的 rRNA 去除方案；小 RNA 等目标还需要相应的片段选择和接头策略。两类文库采样的 RNA 集合不同，其定量分母也受比对、注释和过滤规则影响，不能直接合并 TPM/RPKM 进行比较。同一研究应尽量统一建库策略，并报告文库类型、RNA 质量和去除效率（Wang et al., 2009）。

4.7.2 定量模型：RPKM、FPKM 与 TPM

相对定量可以从读段随机起源的理想模型出发：假定转录本沿长度方向被均匀取样，片段独立进入测序，一条读段(read)来自某转录本的概率近似正比于其摩尔浓度与长度之积。因此，同一文库中等摩尔浓度的长转录本通常产生更多读段。为了用计数估计相对分子丰度，需要同时考虑长度和总取样量；RPKM 用转录本长度与文库大小作除数，给出每 kb、每百万映射读段的密度（Mortazavi et al., 2008）。更精细的模型还会使用有效长度，以考虑片段长度分布和可比对区域：

\[RPKM_{t} = r_{t} / ( \ell_{t} \times N_{M} )\] (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）。比较结果时须统一计数、长度和过滤规则。

实际文库未必符合均匀取样假设。随机引物有局部序列偏好，片段化位置也可能不均匀；PCR 对不同 GC 含量片段的扩增效率不同；oligo(dT) 起始和 RNA 降解还可能造成 3′ 偏置。这些因素既会改变期望计数，也会影响方差和局部覆盖，单纯除以长度不能消除它们。原始计数随长度增加则是上述理想模型本身的预期，不应与这些模型偏离混为一谈。转录本重构和低计数定量尤其需要检查覆盖是否符合模型（图 4.7-2）。

RPKM/FPKM 的总和没有预先固定的值。文库中不同长度转录本的相对丰度改变时，总读段数与长度校正后密度之和的关系也会改变，因此不同样本的 RPKM/FPKM 总和可以不同。TPM(transcripts per million)先按长度校正计数，再以所有纳入定量的转录本密度之和归一化，使结果具有明确的相对占比解释：

\[TPM_{t} = ( r_{t} / \ell_{t} ) / \sum_{u} ( r_{u} / \ell_{u} ) \times 10^{6}\] (4.7-2) 等价换算：TPMt = RPKMt / Σu RPKMu × 106。求和须覆盖该样本纳入定量的全部转录本，且分母大于 0，此时 TPM 总和为 106。在采样模型成立时，TPM 近似表示该集合中的转录本摩尔分数（Li & Dewey, 2011）。它统一了占比的标度，但仍受文库构成、技术偏差和注释差异影响，不能直接代表每个细胞的绝对 RNA 分子数。
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 者的四倍；图中计数除以长度后得到相同 RPKM。B 若计数服从泊松分布，期望计数为 λ 时至少观察到一条读段的概率为 1 − e^(−λ)；λ = 3 时约为 95%。这一概率只表示理想抽样下的观察机会，不能替代真实实验的检出限或定量精度评估。悬停图形元素可查看数值。

增加测序深度(sequencing depth)可提高低丰度转录本被取样的机会。在文库组成和取样效率不变的理想模型中，以 ρ 表示期望 RPKM 密度，则深度 N（百万）下的期望计数为 λ = ρ × ℓ(kb) × N。若某 2 kb 转录本的 ρ = 0.05，在 10 M 文库中 λ = 1，泊松模型给出的至少一次检出概率约为 63%。深度加倍使 λ 加倍；若检出概率阈值固定，达到该阈值所需的最低 ρ 与深度及长度成反比。真实检出限还受背景、比对可靠性和文库复杂度限制，不能仅由这个公式确定。泊松计数的相对标准差为 1/√λ，期望计数仅为个位数时，即使能够检出，定量仍有较大抽样误差。

习题 4.7-1

转录本 A 长 4 kb，在样本甲的 2 000 万条总映射读段中得到 1 600 条唯一比对读段。(1) 计算其 RPKM。(2) 样本乙总映射 1 000 万条读段，其中另一转录本 B 长 2 kb，得到 200 条唯一比对读段。B 的 RPKM 是多少？这两个数值能否直接说明 A 与 B 的绝对分子丰度，或证明某个转录本在两样本间差异表达？(3) 若对样本乙同一文库加深测序至 2 000 万条读段，B 的计数近似加倍至 400，RPKM 如何变化？

参考解答

(1) A 的 RPKM = 1 600 / (4 × 20) = 20。(2) B 的 RPKM = 200 / (2 × 10) = 10，因此 A 的归一化读段密度是 B 的两倍。这里比较的是不同样本中的不同转录本；仅凭这两个数值，不能判定每个细胞中的绝对分子数，也不能检验某个转录本在两样本间是否差异表达。除建库偏差外，文库构成偏置（compositional bias）也会改变相对计数：少数转录本占比增加，可使其他转录本的 RPKM 降低，即使后者的绝对分子数不变。(3) B 的 RPKM = 400 / (2 × 20) = 10。对组成不变的同一文库加深测序，RPKM 的期望保持不变；在独立泊松取样近似下，相对抽样误差随期望计数增加而减小。

交互演示 · RPKM 与 TPM 换算计算器
修改输入可查看长度和深度校正的效果。RPKM 使用样本总映射读段数；右列仅在 A、B 两条转录本组成的子集中重新归一化，不能当作全转录组 TPM。真实 TPM 的分母需要全部纳入定量的转录本。默认值在理想模型下对应 A、B 等摩尔。
转录本读段数 r长度 ℓ (kb)RPKM子集 TPM
A——
B——
A、B 的 RPKM 之和 = —。和大于 0 时，未取整的子集 TPM 总和为 1 000 000；显示值可能有 1 的舍入误差。默认参数得到两个 10 RPKM，子集 TPM 各为 500 000。若两条读段数均为 0，RPKM 为 0，子集 TPM 因分母为 0 而无定义，显示“—”。
文库深度 N 与长度 ℓ 须为正数，读段数 r 不得为负。
习题 4.7-2

在仅包含两条转录本的简化定量集合中，RPKM 分别为 30 与 10。求各自的 TPM，并解释 TPM 总和为何为 106，而 FPKM 总和没有这一性质。如果样本实际包含更多转录本，仅知道上述两个 RPKM 是否足以求出它们在全转录组中的 TPM？TPM 的固定总和是否保证跨样本的绝对丰度可比？

参考解答

在给定的两转录本集合中，TPM₁ = 30/(30+10) × 106 = 750 000，TPM₂ = 10/40 × 106 = 250 000。TPM 将长度校正后的全部密度除以其总和，再乘以 106，所以分母非零时总和固定。FPKM 以总映射片段数归一化，其总和还取决于表达构成与转录本长度。如果存在其他纳入定量的转录本，必须获得它们的密度之和，才能计算全转录组 TPM。固定总和有助于解释相对占比，但不能消除文库构成偏置，也不保证绝对丰度可比；基因级差异表达检验通常使用计数及适当的归一化和方差模型（4.7.4 节）。

4.7.3 比对与转录本重构

读段可以先比对到参考基因组，再据其位置分析转录区域。完整位于单一外显子内的读段可作连续比对；跨越外显子接头的读段则对应基因组上被内含子分隔的区域，间隔可能长达数十 kb。不支持剪接的连续比对模式通常无法正确处理这类读段，因此需要剪接感知比对(splice-aware alignment)：允许读段在基因组坐标上跨过较长间隙，在比对记录中常以 N 表示跳过的区域，并结合剪接位点序列和读段证据评估候选接头。按转录方向书写时，多数内含子使用 GT-AG 边界，也有 GC-AG、AT-AC 等类型；候选筛选不能简单排除所有非典型位点。这与 2.5 节讨论的基因结构和剪接信号相联系。早期 TopHat（2009）先定位可连续比对的读段，再构建候选外显子接头，尝试比对其余读段。

没有合适参考基因组时，可采用从头组装(de novo assembly)，根据读段之间的序列重叠重构转录本。Trinity（Grabherr et al., 2011）依次运行三个主要模块：Inchworm 以 k-mer 为基础贪婪延伸出重叠群；Chrysalis 根据序列重叠及跨接读段证据聚类，并为各组构建 de Bruijn 图；Butterfly 再结合读段和配对信息分析图中的路径，区分可能的转录本异构体（isoform）及旁系同源转录本。该方法便于研究缺乏参考基因组的物种，但低覆盖、测序错误和高度相似序列仍可能造成缺失、断裂或错误组装，结果需要结合读段支持、表达估计及同源证据评估。

表达估计需要处理读段归属的不确定性。一个基因的多个 isoform 常共享外显子，落在共享区域的读段无法仅凭序列唯一归属。概率定量将各候选转录本的丰度、有效长度和产生该读段的概率纳入模型，而不是把歧义读段一律平均分配。RSEM 用期望最大化（EM）算法迭代估计丰度：先计算每条读段来自各候选的后验权重，再据此更新参数，并可建模片段长度、测序错误和位置偏好（Li & Dewey, 2011）。Cufflinks 先利用参考基因组上的比对信息组装候选转录本，再以概率模型估计丰度；原始工作还展示了新 isoform 的重构和条件相关的 isoform 切换（Trapnell et al., 2010）。

注记

一次 EM 分配的数值示例。设有 10 条读段同时与 isoform 甲、乙兼容，另有 5 条只与甲兼容。若某轮迭代的丰度权重比为 3:1，且每条歧义读段在两候选下的条件概率相同，则 E 步将这 10 条读段分配为甲 7.5、乙 2.5；加上特异读段，期望计数为甲 12.5、乙 2.5，总数仍为 15。M 步再结合模型参数更新丰度，下一轮的比例可能改变。这个例子只说明一次分配，不能仅凭“兼容”关系确定最终丰度。EM 提高模型似然，但不能补充原始数据中缺失的区分信息；高度相似的候选可能仍不可辨识。

4.7.4 差异表达的计数统计

以明确指派到基因的读段为例，可得到“基因 × 样本”的整数计数矩阵。计数的均值与方差相关，因此分析时需要采用合适的抽样模型。在分子独立取样、取样概率稳定且目标占比较小时，泊松分布可作为技术抽样的基线。生物学重复还包含个体间或样本间的表达差异，常使条件内方差超过泊松模型的均值，即过度离散(overdispersion)；建库等技术变异也可能贡献额外方差（图 4.7-3）。若确有过度离散，却仍使用不包含额外变异的泊松检验，通常会低估不确定性、产生过小的 p 值。

\[\operatorname{Var}(K_{gj}) = \mu_{gj} + \varphi_{g} \cdot \mu_{gj}^{2}\] (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、φ = 0.1 时，两模型的方差比为 101，约一百倍；这不意味着 p 值或显著性也相差同样倍数，检验结果还取决于效应大小、重复数和检验方法。悬停点集可查看说明。

各基因的离散度可能不同，而生物学重复数往往有限，例如每组仅三至六个时，单独估计每个基因的方差仍很不稳定。DESeq 与 edgeR 都利用跨基因的信息改善估计，但实现方式不同：DESeq 的 2010 年方法拟合均值与额外方差的平滑关系；edgeR 则采用经验贝叶斯方式调整逐基因离散度，并提供共同、趋势化或逐基因等估计方式。这些方法以总体信息稳定有限重复下的方差估计，与微阵列分析中 limma 的方差收缩相联系；使用时仍需依据具体软件版本说明模型与估计规则。

差异分析还要处理文库构成偏置。总测序量固定时，少数高丰度转录本占比大幅增加，会减少分配给其他转录本的读段比例，即使后者的绝对分子数未变。因此，简单按总计数归一化未必能得到合适的比较尺度。DESeq 的中位数比率法通过基因间的稳健汇总估计样本尺度：

\[s_{j} = median_{g} ( k_{gj} / ( \prod_{j'} k_{gj'} )^{1/m} )\] (4.7-4) kgj：基因 g 在样本 j 的计数；分母为该基因在 m 个样本中的几何平均，作为参考值；对可计算的基因比值取中位数，得到尺寸因子 sj。本式不适用于几何平均为 0 的基因，应按所用方法的规则排除或处理。中位数比率法与 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，样本方差是它的 9 倍，即额外方差为泊松方差的 8 倍；由 900 ≈ 100 + φ·100² 得 φ ≈ 0.08。(3) 若两组的方差比均近似如此，且重复数等条件相同，用泊松方差计算的标准误约只有实际值的 1/3（√9 = 3），可能夸大组间差异的证据。但仅三个重复得到的方差和 φ 都很不稳定，这个例子不能给出实际检验的 p 值。分析中应采用能够描述额外变异的模型，并用跨基因信息改善方差估计。

4.7.5 与微阵列的系统比较

早期研究为测序定量提供了不同层面的实证。Mortazavi 等（2008）研究的是成年小鼠脑、肝和骨骼肌的 poly(A) RNA，并用 RNA 标准评估线性检测范围；该实验覆盖约五个数量级，也通过跨接头读段观察剪接。Nagalakshmi 等（2008）则用酵母 poly(A) RNA 改进了 UTR、内含子和重叠转录区域的注释，不能据此声称其方法完整检测了非聚腺苷化转录本。

't Hoen 等（2008）将同一批野生型和 DCLK-short 转基因小鼠海马 RNA（每组四只）的 3′ 标签 DGE 结果，与 Applied Biosystems、Affymetrix、Agilent、Illumina 和自制寡核苷酸微阵列五种平台比较。在该实验中，DGE 检出了更多低丰度及差异表达转录本。跨实验室重复性检验则使用这些 RNA 的混合样本，在另一地点按同一方案测量，得到较高一致性。这些结果支持当时深度序列计数的应用潜力，但不是所有 RNA-seq 流程在所有样本和指标上均优于芯片的证明。表 4.7-1 将测量原理与这一时期的实验表现分开列出。

表 4.7-1微阵列与 RNA 测序的比较（以 2008–2009 年研究为背景，统计方法补充至 2010 年；据本节所列原始研究与综述整理）
比较维度微阵列RNA 测序
测量原理测量探针杂交荧光，受背景、交叉杂交和饱和限制；部分早期平台的有效动态范围约 102，并非固定上限记录序列计数，早期实验的动态范围可跨 104 以上；实际范围仍受有效深度、文库复杂度及分析流程限制
先验序列需要预先设计探针，能够检测的序列受探针覆盖限制可借助参考序列识别已知转录本及新转录区域；无合适参考时可尝试从头组装
新转录本与融合基因对未设计探针的序列缺乏直接检测能力；专门设计可支持特定发现任务跨接头或嵌合读段可提示新剪接与融合转录本；编辑、突变和融合候选仍需排除比对及建库伪影
等位特异性表达SNP 可影响探针杂交效率，也可通过专门探针设计研究等位差异可利用信息性变异及相位分配读段，但需要足够覆盖，并校正参考等位基因比对偏倚
剪接定量通常需外显子或接头探针，结果受探针特异性和覆盖范围影响接头读段提供局部剪接证据；完整 isoform 丰度仍常依赖概率推断（4.7.3 节）
统计模型连续强度：线性模型（limma）+ 经验贝叶斯（见 4.4 节）基因级计数常用负二项模型（DESeq/edgeR），并处理文库尺度和构成差异
当时的成本结构平台与流程较成熟，既定任务中的预算和分析需求较易估计成本随测序深度、样本数和流程而变，当时通常需要更多计算与存储资源
成熟度（2009 年）已有 MIAME 报告规范与 GEO 数据共享实践，跨平台和跨实验室差异仍需处理建库和分析流程仍在发展；'t Hoen 等的标签 DGE 展示了特定样本及一致方案下的跨实验室重复性
案例

2008–2010 年的平台选择。已有微阵列数据需要持续随访、研究目标由既定探针充分覆盖时，沿用同一平台有利于控制测量差异，但仍须评估批次效应。若目标是发现新的转录区域、剪接形式或研究缺乏成熟探针的物种，RNA-seq 提供了更广的序列信息，代价是相应的测序和分析投入。选择应依据目标、样本质量、所需精度与验证资源，而不能只看检出特征的数量。本节介绍的 Trinity、RSEM 等 2011 年方法是对原书背景的补充，不能视为 2009 年已经普遍建立的实践。

4.7.6 偏差与陷阱

文库和比对质量会直接影响计数的含义。首先应检查 rRNA 残留：例如，若其占比由两成升至八成，在总读段数相同的情况下，留给其他 RNA 的取样量会明显减少；这只是比例示例，不是所有去除方案的残留范围。可按 rRNA、外显子、内含子和基因间区等注释类别统计读段，并明确重叠类别的计数规则。其次，旁系同源基因、假基因和重复序列可能产生多映射读段(multi-mapping reads)。丢弃它们便于使用唯一归属计数，却可能低估缺乏特异序列的成员；概率分配利用候选的比对概率和丰度估计，不必假定各成员表达相等。但在候选难以区分时，结果仍可能高度依赖模型，需报告分配规则及不确定性（Li & Dewey, 2011）。

覆盖不均也需要单独检查。序列偏好、GC 偏置和 3′ 偏置会改变沿转录本的读段分布（4.7.2 节），既可能使 isoform 丰度估计发生系统性偏移，也可能增加不确定性；如果模型未描述这些偏差，置信区间未必相应变宽。沿转录本的覆盖曲线、基因体覆盖和片段分布有助于发现这些问题。

批次与文库复杂度(library complexity)影响可用于比较的有效信息。建库日期、试剂批次、测序通道和流动槽等因素都可能造成技术差异，应尽量在批次内平衡实验分组，并随机安排样本，与 4.3 节的设计原则一致。文库复杂度指文库中独立分子的数量和多样性。起始分子少、步骤损失或部分片段过度扩增，会使测序结果被重复拷贝占据；继续加深测序未必能获得更多独立分子。独特分子标识符（UMI，原书出版后更广泛应用）在扩增前给分子加标签，用于识别和合并相应的扩增拷贝。UMI 可减少重复计数对定量的影响，但不能恢复加标签前丢失的 RNA，也不能自动消除捕获、逆转录或序列偏好；标签碰撞及测序错误同样需要处理。

警示

先确定比较对象，再选择归一化与模型。样本甲深度 40 M，基因 X 得 400 条读段；样本乙深度 10 M，基因 Y 得 200 条读段。400 大于 200 不能直接说明 X 的分子丰度更高：两个文库深度不同，而 X、Y 的长度尚未给出。跨基因比较需要考虑长度和取样偏差；同一基因跨样本的差异分析，则通常从计数出发，在模型中处理文库尺度、构成和实验设计，不能把 TPM 直接当作原始计数输入 DESeq 或 edgeR。表达热图应采用与比较目的匹配的归一化及变换；原始计数仍可用于文库大小、低计数特征等质量检查，只是解释时必须注明其用途。

RNA-seq 的结论由文库所保留的 RNA、读段归属规则和统计模型共同决定。RPKM、FPKM 与 TPM 描述相对定量，可靠的差异表达结论还需要适当的生物学重复、归一化和不确定性估计。理解这些条件后，才能将表达变化与实验问题联系起来。下一节进一步讨论功能基因组学的应用实例及其解释范围。

关键术语

RNA 测序 (RNA sequencing, RNA-seq)
通过 RNA 或其 cDNA 的序列数据研究转录本结构与丰度的技术；本节重点介绍短读长 cDNA 测序。
poly(A) 富集 (poly(A) selection)
以寡聚(dT)等方式捕获带 poly(A) 尾分子的策略，可富集多数成熟真核 mRNA 及部分非编码 RNA。
rRNA 去除 (rRNA depletion)
降低核糖体 RNA 含量、扩大对其他 RNA 取样范围的处理；保留范围仍受后续建库步骤限制。
链特异性文库 (strand-specific library)
保留转录方向信息的 RNA-seq 文库，如以 dUTP 法构建。
读段 (read)
测序产生的一段序列；双端测序的一对读段可对应同一个文库片段。
RPKM/FPKM (reads/fragments per kilobase per million)
按转录本长度与文库大小双重归一化的相对表达定量。
TPM (transcripts per million)
先按长度校正、再归一化的相对丰度；分母非零时，在全部纳入定量的转录本中总和为 106。
剪接感知比对 (splice-aware alignment)
允许读段跨越内含子形成分段比对，并结合剪接信号和读段证据评估接头的策略。
从头组装 (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.
  11. Parkhomchuk D, Borodina T, Amstislavskiy V, et al. 2009. Transcriptome analysis by strand-specific sequencing of complementary DNA. Nucleic Acids Research 37: e123.