第四章：负对照基因的选择与构建

4.1 为什么我们需要负对照基因？

RUVseq 的核心思路其实挺朴素的——它想从数据里偷偷学出批次效应的「形状」，然后把它减掉。但问题来了：你怎么知道哪些变化是批次带来的，哪些是生物学信号？

答案就在负对照基因（Negative Control Genes）上。这些基因有个特点：它们不受你关心的生物学条件影响。说白了，它们的表达差异只来源于批次效应。我刚开始接触这个思路时觉得有点绕，但想通了就很简单——你拿一批已知不受处理影响的基因，它们在不同批次间的波动，就是批次效应的「代言人」。

那怎么定义「不受生物学条件影响」？通常有两种思路：

• 要么你事先知道哪些基因是「管家基因」，在不同条件下表达稳定
• 要么你从实验设计里找——比如用空白对照、或者用很多样本跑一次差异分析，挑出那些p值很大、没明显变化的基因

我个人习惯第二种方式，因为更贴合实际数据。你给我一个真实的RNA-seq数据集，我很难百分之百保证某个基因在所有条件下都不变，但用统计方法筛出来的负对照基因更靠谱。

4.2 构建负对照基因的两种主流策略

这里我整理了一张图，你可以看到整个选择过程是怎么迭代的：

上图展示了两条主路径。我先说说具体的实现细节。

4.3 策略一：基于先验知识的低技术门槛法

如果你手头有成熟的参考基因集（比如 HRT Atlas 里列出的那些管家基因），直接用就行。我在做小鼠组织样本时常用这个方案，很省事。

典型的做法是构建一个基因ID的向量：

# 假设你已知一组管家基因的ensembl ID
housekeeping_genes <- c("ENSG00000111640", 
                        "ENSG00000184009",
                        "ENSG00000113575",
                        "ENSG00000139618")

# 在你的表达矩阵中提取这些基因的表达
ctrl_genes <- rownames(se) %in% housekeeping_genes

# 传给RUVseq
seq <- RUVg(se, ctrl_genes, k=5)

不过要小心一点：这些管家基因不一定在所有批次里都稳定。我以前遇到过一个案例，某个所谓的「万能管家基因」在某种应激条件下表达量翻了一倍，结果把批次校正搞砸了。所以我的建议是——先验集可以当种子，但最好再用数据验证一下。

4.4 策略二：数据驱动的自筛选法（推荐）

这个方法是正路。你想想看，与其依赖外部知识，不如直接从你自己的数据里找那些在组间没有差异的基因。具体怎么做？

1. 先跑一轮差异分析——用 DESeq2 或 edgeR 都可以，给出每个基因的 p-value 和 padj。
2. 设定一个宽松的过滤阈值——比如 padj > 0.9，或者 p-value > 0.5。这里不要卡太严，否则样本量不够。
3. 再做一些质量过滤——去掉表达量太低的基因（它们在哪个组都不变，但也没信息），去掉极高离散度的基因。

我自己的经验：用 padj > 0.95 起步，如果筛出来的基因少于 200 个，就放宽到 0.9；如果太多（超过 5000 个），可以收紧到 0.99。下面是一个具体实现：

library(DESeq2)
# 假设你已经有了 dds 对象
dds <- DESeq(dds)
res <- results(dds, alpha = 0.05)

# 选择负对照基因：padj > 0.95 且 baseMean > 10
ctrl_idx <- which(res$padj > 0.95 & 
                  res$baseMean > 10 & 
                  !is.na(res$padj))
ctrl_genes <- rownames(res)[ctrl_idx]

cat("筛选到", length(ctrl_genes), "个负对照基因\n")

4.5 避坑指南：我踩过的三个坑

坑一：把 housekeeping 基因当万能钥匙。 我曾经在一个细胞分化实验里直接用了文献里的管家基因列表，结果校正后反而把分化信号都抹掉了。后来发现那些管家基因在分化早期确实有波动。用自己的数据再筛一遍，稳妥很多。

坑二：负对照基因里混入了「假稳定基因」。 比如某个基因在所有样本里表达量都是 0，它当然没有差异，但也没带来任何信息。这种基因会让 RUVseq 估计出的批次效应偏小，校正效果打折。所以我通常会加一条 baseMean > 10 的过滤。

坑三：忽略了批次内重复。 如果你只有每个条件一个样本，那就没法区分批次效应和生物学信号。负对照基因再多，RUVseq 也分不清。这点在实验设计阶段就要想好——每个批次至少两个生物学重复。

4.6 如何验证你选的负对照基因是否合格？

选完之后别急着往下跑，先做一次可视化快检：

• 热图检查：用 pheatmap 把负对照基因的表达画出来，看看是不是按批次聚集而不是按组别聚集。如果看到明显的组别分块，说明混入了生物学差异基因，得重新筛选。
• PCA 检查：只用负对照基因做 PCA，第一批主成分应该揭示批次结构，而不是实验条件。

这里给一个检查代码：

library(pheatmap)
# 提取负对照基因的表达矩阵（已标准化）
expr_ctrl <- assay(se)[ctrl_genes, ]
# 去掉行均值和列均值后画热图
pheatmap(expr_ctrl, scale = "row",
         annotation_col = colData(se)[, "batch", drop=FALSE])