📚 课程目录第4章差异表达分析

4、差异表达分析：edgeR、DESeq2与特征初筛

拿到RNA-seq的count矩阵后，很多人第一反应就是：跑差异表达。但说实话，这一步才是整个分析中最容易翻车的环节。你想想看，不同方法背后的统计模型不一样，数据假设也不一样，用错了工具，结果就差得离谱。

我做肿瘤亚型项目的时候，遇到过好几次这样的场景：两种算法跑出来的差异基因列表只有一半重叠。当时我就意识到，理解每个工具的脾性，比盲目跑脚本重要得多。

4.1 差异表达分析的核心问题

我们先明确一下：差异表达分析到底在做什么？

说白了，就是比较两组样本（比如肿瘤A vs 肿瘤B，或者肿瘤 vs 正常）中，每个基因的表达量有没有统计学显著差异。但RNA-seq数据有个特点——它基于短读长计数，服从的是负二项分布，不是正态分布。如果你拿t检验硬上，p值会严重偏小，假阳性一大堆。

所以业内主流工具都是基于负二项模型，加上各种均值和方差的估计矫正。最常见的两个就是edgeR和DESeq2。我个人的习惯是：小样本量（每组少于5）首选edgeR，大样本量（每组10以上）两个都跑，然后取交集。

特性	edgeR	DESeq2
统计方法	精确检验 + 经验贝叶斯	负二项GLM + 收缩估计
离散度估计	加权经验贝叶斯	分步式收缩（先全局，再基因级）
小样本表现	较好（对离群值鲁棒）	保守，但p值校准更准
数据预处理要求	推荐先TMM归一化	内置大小因子归一化
运行速度	快（约5-10分钟10k基因）	慢一些（约20-30分钟）

技巧： 我一般在项目初期先用edgeR快速筛选一轮（因为快），留下候选基因后，再用DESeq2做精细验证。这样可以平衡速度和准确性。

4.2 edgeR实战：从count矩阵到基因列表

我们先走一遍edgeR的完整流程。假设你已经有了一个count矩阵，行是基因，列是样本。这里给一个标准的R脚本片段：

library(edgeR)

# 读取count矩阵（假设行为基因，列为样本）
counts <- read.table("count_matrix.txt", header=TRUE, row.names=1)

# 构建分组信息
group <- factor(c(rep("TumorA",3), rep("TumorB",3)))

# 创建DGEList对象
dge <- DGEList(counts=counts, group=group)

# 过滤低表达基因（关键步骤）
keep <- filterByExpr(dge, min.count=10, min.total.count=15)
dge <- dge[keep, , keep.lib.sizes=FALSE]

# TMM归一化
dge <- calcNormFactors(dge, method="TMM")

# 估计离散度
dge <- estimateDisp(dge)

# 精确检验
et <- exactTest(dge)

# 提取结果
results <- topTags(et, n=Inf)

这里面有几处坑，我踩过不止一次：

过滤低表达基因：你别把这一步省了。不过滤的话，低表达的基因会贡献大量噪声，影响离散度估计。我曾经在一个项目里忘了过滤，结果跑出来2000个“显著基因”，后来发现全是测序背景噪音。
TMM归一化：edgeR默认的归一化方法。它假设大部分基因不是差异表达的，用这个假设来矫正文库大小。如果你的数据里有大量差异基因（比如肿瘤vs正常），TMM可能会出问题。我建议先用密度图看看整体分布，如果偏差太大，可以考虑改用RLE或手动设置归一化因子。
精确检验 vs 似然比检验：如果你的分组多于2个，或者有批次效应等协变量，建议用glmFit + glmLRT代替exactTest。

注意： edgeR的精确检验exactTest只适用于两组比较（如TumorA vs TumorB）。如果你要做配对设计或多因素分析，必须切换至GLM工作流。

4.3 DESeq2实战：稳健的替代方案

DESeq2的逻辑和edgeR不太一样。它不依赖TMM归一化，而是自己算一个“大小因子”（size factor），然后用负二项广义线性模型拟合每个基因。

library(DESeq2)

# 准备Count矩阵，注意DESeq2要求整数
counts <- read.table("count_matrix.txt", header=TRUE, row.names=1)

# 构建colData（包含样本分组信息）
colData <- data.frame(condition = c("TumorA","TumorA","TumorA",
                                     "TumorB","TumorB","TumorB"),
                       row.names = colnames(counts))

# 创建DESeqDataSet对象
dds <- DESeqDataSetFromMatrix(countData = counts,
                               colData = colData,
                               design = ~ condition)

# 过滤低表达基因（DESeq2推荐最少3个样本总和>=10）
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]

# 运行差异分析
dds <- DESeq(dds)

# 提取结果（对比TumorB vs TumorA）
res <- results(dds, contrast=c("condition", "TumorB", "TumorA"))

我个人的经验是：DESeq2的p值校准更保守，假阳性更低。但代价是灵敏度不如edgeR。如果你的项目需要高召回率（比如做差异基因筛选后还要做功能富集），我建议用edgeR的结果。如果是为了后续机器学习建模，希望特征更干净，那就用DESeq2。

不过有一点要注意——DESeq2的DESeq()函数默认包含了离群值检测和自动替换。如果你发现某些基因的p值异常小（比如< 1e-300），可能是离群样本导致的。可以检查一下assay(dds, "cooks")矩阵，看看Cook距离是否过大。

4.4 特征初筛：从差异基因到候选特征

跑完差异分析，你会得到一个表，每行一个基因，包括log2FC、p值、FDR（矫正后的p值）。通常我们按以下阈值初筛：

FDR < 0.05 —— 控制假发现率
|log2FC| > 1 —— 相当于2倍表达量变化
baseMean > 20 —— 排除极低表达基因（可选）

但实际项目中，肿瘤亚型分类往往不只需要差异显著。我还见过很多情况：某基因虽然不显著，但它的表达量在两组间分离得很好（比如A组全部大于10，B组全部小于1）。这种基因对分类器也很重要。所以我在做特征初筛时，会额外加一个指标——AUC值：

# 计算每个基因在两组间的AUC（区分能力）
library(pROC)

calc_auc <- function(gene_counts, group) {
  roc_obj <- roc(response=group, predictor=gene_counts)
  auc(roc_obj)
}

# 假设dge对象已过滤归一化
expr_norm <- cpm(dge, normalized.lib.sizes=TRUE)
group_binary <- ifelse(dge$samples$group == "TumorA", 0, 1)
auc_vec <- apply(expr_norm, 1, function(x) calc_auc(x, group_binary))

我会筛选出AUC > 0.7的基因，和差异基因取并集。这样既能保证统计学显著，又能捕捉那些区分度好但不显著的特征。在肿瘤亚型项目中，我经常用这种策略拿到更鲁棒的特征集。

  核心策略： 差异表达分析是做减法，特征初筛是做加法。前者控制假阳性，后者扩大候选池。两者结合才能为后续机器学习提供足够优秀的特征。

4.5 知识体系总览

下面这张图把本章的逻辑串起来了：

4.6 避坑指南

最后总结几个我踩过的坑，你们遇到类似情况可以少走弯路：

坑1： 记得检查count矩阵是否有NA或重复的基因名。我碰到过一次因为Excel保存时把基因名自动转化成了日期格式，结果读进来都是乱码。
坑2： edgeR的exactTest对低表达基因的假阳性控制不够好。如果你发现差异基因太多（比如超过总基因的50%），先检查过滤阈值是不是太松了。
坑3： 做归一化前一定要确认文库大小有没有极端的离群样本（比如一个样本的total count是其他样本的10倍）。我遇到过技术重复跑偏的情况，归一化因子会变得离谱。解决办法：手动删除离群样本或用上限截断法。
坑4： 特征初筛时别只盯着p值。我见过很多p值显著但表达量极低的基因，在分类器里全变成了噪声。加上baseMean过滤或者基因表达变异系数过滤（CV > 0.3）会好很多。

4.7 本章小结

差异表达分析是肿瘤亚型分类特征工程的起点。edgeR和DESeq2就像两个老朋友——各有长短，但都能帮你拿到核心差异基因。我建议你上手时先跑通一个端到端的流程（count→过滤→归一化→检验→结果输出），再去调参数。这个基础打牢了，后面的机器学习建模才站得稳。