4. 估计大小因子（sizeFactors）：方法与影响

做差异表达分析，第一步就是归一化。这一步要是错了，后面再漂亮的统计检验也救不回来。DESeq2 估计大小因子的方法，说白了就是帮你把不同样本的测序深度拉齐。我刚开始接触 RNA-seq 时觉得这步就是个「库大小校正」，直到踩过几次坑才明白——细节远比想象中多。

4.1 为什么要单独估计大小因子？

你可能会想：直接用总 reads 数归一化不就好了？嗯，这恰恰是新手最容易犯的错。总 reads 数对样本里高表达的少数基因特别敏感。假设样本 A 里有一个基因表达量爆表，总 reads 数就跟着虚高，其他基因的归一化表达量反而被压低了。差异分析的结果？一塌糊涂。

DESeq2 的 sizeFactors 只依赖所有基因的表达量，对高表达基因做了「中位数」处理，稳得很。我个人的经验是：只要文库构建没问题，DESeq2 默认估计出来的 sizeFactors 极少翻车。

核心思想：每个样本的 sizeFactor 反映其整体测序深度与表达谱的「相对关系」。不是简单的文库大小比例，而是所有基因中位数比值的代表值。

4.2 默认方法：比值中位数法（Median of Ratios）

DESeq2 的 estimateSizeFactors() 函数，默认就是比值中位数法。流程分三步：

  1. 构造伪参考样本：对每个基因，计算所有样本的几何均值。
  2. 计算比值：对于每个样本、每个基因，用原始 count 除以该基因的几何均值。
  3. 取中位数：每个样本的 sizeFactor = 该样本所有基因比值的中位数。

为什么用几何均值而不是算术均值？因为 counts 是正偏态分布，几何均值能压制极高表达基因的影响。记得有一次我在项目里看到某个样本的 sizeFactor 异常大，排查后发现是文库预处理时混入了接头二聚体——中位数比值直接被稀释了。

实际用法：

library(DESeq2)

# 假设你已经创建好了 dds 对象
dds <- estimateSizeFactors(dds)

# 查看估计结果
sizeFactors(dds)

这一行代码背后就在跑上面的三步。如果想偷懒，DESeq() 函数会自动调用。

4.3 其他估计方法：你有哪些选择？

DESeq2 允许你手动指定归一化策略。下面这些是我在实际项目里用过的变体：

方法 实现方式 适用场景
默认比值中位数 estimateSizeFactors(dds) 标准 Bulk RNA-seq，推荐使用
自定义几何均值基因集 geoMeans 参数 剔除某些基因（如线粒体、核糖体）后做分母
指定定位函数 locfunc 参数 样本量极小（如 n=2）时改用 shorth 函数
手动设置 sizeFactors sizeFactors(dds) <- values 需要统一归一到外部标准时

举个具体例子：当样本只有 2 个重复时，默认的中位数容易飘。我会把 locfunc 改成 shorth——这是一个更鲁棒的位置估计量，能稍微提高稳定性。不过这招属于「微调」，不是万金油。

我曾经踩过的坑：有一次我手动把 sizeFactors 设为文库大小的比例（Total Counts），结果差异基因数目暴增 30%。仔细一查，那些多出来的基因全是高表达的—因为总 reads 归一化本身就没有剔除高表达基因的「自我膨胀」效应。从此以后，除非有极其特殊的理由，我绝不手动绕过 DESeq2 的默认方法。

4.4 如何评估 sizeFactors 是否合理？

光跑完函数还不够，我建议做两个快速检查。

  • 检查 sizeFactors 的范围：正常应在 0.3 ~ 3 之间。如果某个样本大于 3 或小于 0.3，说明该样本测序深度或文库质量与其他样本差异巨大，需要回看 QC。
  • 画箱线图看比值分布：对每个样本，画出 log2(count / 伪参考) 的箱线图。所有箱体的中位数应该落在 0 附近。如果某个样本的箱子整体偏移，sizeFactor 估计可能不准。

我自己常用这段代码快速诊断：

# 导出比值分布
normFactors <- estimateSizeFactors(dds)
ratios <- counts(dds, normalized = FALSE) / 
          rowMeans(counts(dds, normalized = FALSE), na.rm = TRUE)
boxplot(log2(ratios), las = 2, main = "Ratio distribution per sample")
abline(h = 0, col = "red", lty = 2)

你会看到什么？如果所有样本的中线都在红线附近，OK。如果某个样本的盒子整体下沉或上浮，说明它的 sizeFactor 可能没被正确捕获。我曾经在一个批次效应严重的项目里发现，其中一个处理组的样品因为文库 DNA 降解导致比值整体偏低，手动剔除后差异分析才正常。

4.5 核心流程一览

下面这张图把 sizeFactors 估计的整个数据流画了出来。从原始 count 矩阵到最终的归一化系数，每一步的逻辑都在这了。

大小因子估计流程 原始 count 矩阵 每基因几何均值（伪参考） 每个样本每个基因：count / 几何均值 每个样本取比值中位数 → sizeFactor 原始 count / sizeFactor = 归一化 count 注意：中位数替代均值是为了抵抗高表达基因的杠杆效应；几何均值构造伪参考比算术均值更稳健

4.6 sizeFactors 对差异结果的实际影响

差异表达的检验统计量（Wald 统计量或 LRT）直接依赖归一化之后的 counts。如果 sizeFactors 偏了，后果有：

  • 假阳性增加：一个样本的 sizeFactor 被低估，它自身的归一化 counts 就会虚高，导致它和其他样本对比时产生虚假的 fold change。
  • 假阴性增加：真实的差异被错误归一化掩盖掉。有时候你跑完差异分析发现基因全不显著，先别急着怀疑实验，看看 sizeFactors 是不是被某个异常样本拉偏了。
  • p值分布畸形：正常情况下 p-value 直方图应该接近均匀分布（零假设部分）。如果 sizeFactors 不对，p-value 会在 0~0.05 区间出现异常峰或完全平坦。

一句话总结：DESeq2 默认的比值中位数法是久经考验的。别轻易改它，除非你明确知道自己在做什么。我见过太多人因为「想优化」而手动调整，最后改出一个更差的模型。

4.7 特殊情况处理：直接指定 sizeFactors

有时候实验设计特殊，比如 spike-in 对照（外源 ERCC）或者你需要对齐到某个外部参考。DESeq2 允许你直接赋值：

# 假设你通过 ERCC 或 qPCR 算出了样本间的校正系数
manualSF <- c(0.8, 1.0, 1.2, 0.9, 1.1)
sizeFactors(dds) <- manualSF
dds <- estimateDispersions(dds)  # 重新估计离散度
dds <- nbinomWaldTest(dds)        # 重新检验

注意：一旦手动设置 sizeFactors，必须重新估计离散度（estimateDispersions）和统计检验（nbinomWaldTest），因为归一化后的 counts 变了，离散度估计也会跟着变。这个顺序不要搞反——我有个同事直接把 dds 跑 DESeq() 覆盖了手动设置，结果等于没改。

另外，手动设置的系数不能是负数或 0，否则 DESeq2 会报错。每个样本最好在 0.1 ~ 10 之间，太极端会破坏模型假设。

避坑指南：如果你同时用 ERCC 做了滴定实验，别直接拿 ERCC 的 counts 除以总 reads 作为 sizeFactor。ERCC 的拷贝数在不同样本间往往是恒定的（理论上），但实际加样量可能有误差。我建议用 ERCC 推导的系数作为参考，再和默认 sizeFactors 对比，如果差异超过 2 倍，优先排查实验操作。

总的来说，sizeFactors 估计是 DESeq2 里最不起眼但影响最深远的步骤。花 5 分钟检查一下，可能省掉你后面一周的返工。嗯，今天就聊到这儿。

← 上一章 📖 返回目录 下一章 →