Scientify选方法
三者在常规两组比较上结果高度重叠，选错的代价主要出现在极端样本量和数据类型上。
本页实测（R 4.5.3、DESeq2 1.50.2、edgeR 4.10.5、limma 3.68.5、clusterProfiler 4.18.4，Apple M2 8 核 16 GB）：DESeq2 得到 951 个差异基因，edgeR QL 1183 个，limma-voom 1145 个；三者交集 910 个，并集 1227 个。差异集中在效应量接近 1 的边界基因上。论文中写明所用方法和版本即可，不需要取三种方法的交集。
library(edgeR); library(limma)
design <- model.matrix(~ cell + dex, data = as.data.frame(colData(airway)))
y <- DGEList(assay(airway)); keep <- filterByExpr(y, design)
y <- normLibSizes(y[keep, , keep.lib.sizes = FALSE])
## edgeR 4 的 quasi-likelihood 流程（手册推荐）
fit <- glmQLFit(estimateDisp(y, design), design)
tt <- topTags(glmQLFTest(fit, coef = "dextrt"), n = Inf)$table
## limma-voom
v <- voom(y, design)
tl <- topTable(eBayes(lmFit(v, design)), coef = "dextrt", n = Inf)| 情况 | 推荐 | 依据与选错的后果 |
|---|---|---|
| RNA-seq counts，每组 3–6 个生物学重复 | DESeq2 或 edgeR（QL 流程） | 两者都对离散度做跨基因收缩，小样本更稳；本页实测三种方法的差异基因交集占并集的 74% |
| 每组 2 个重复 | DESeq2 或 edgeR，结论放宽为探索性 | 能跑，但 Cook 距离离群检测需要每组至少 3 个重复，2 个重复时离群样本无法被标记 |
| 没有生物学重复 | 只做描述性分析；edgeR 手设 BCV 得到候选列表 | DESeq2 报错退出；edgeR 结果随 BCV 取值变化 100 倍（见下一节） |
| 样本多（每组 20 个以上）或多因素、连续协变量 | limma-voom | 线性模型拟合快，支持 duplicateCorrelation 处理重复测量 |
| 人群队列，每组数十到数百例 | Wilcoxon 秩和检验，或 DESeq2/edgeR 加 SVA/RUV 校正技术变异 | Li 等 2022 年在 13 个人群数据集上做置换检验，DESeq2/edgeR 在目标 FDR 5% 时实际 FDR 有时超过 20% |
| 芯片（microarray）表达矩阵 | limma | 数据已是 log 尺度连续值，不能输入 DESeq2/edgeR |
| TPM、FPKM 或已标准化的矩阵 | 找回原始 counts；只有 TPM 时用 limma-trend 并在方法中说明 | DESeq2/edgeR 需要整数 counts；把 FPKM 取整后输入会让离散度估计失真 |
无重复
没有重复就无法估计组内生物学变异，任何软件给出的 p 值都依赖人为假设。
DESeq2 自 1.22 起不再支持无重复设计。本页用 airway 中同一细胞系的 1 个处理样本对 1 个对照样本运行 DESeq()，报错原文是：“The design matrix has the same number of samples and coefficients to fit, so estimation of dispersion is not possible. Treating samples as replicates was deprecated in v1.20 and no longer supported since v1.22.” 2018 年以前的中文教程里“无重复也能用 DESeq2”的做法已经失效。
edgeR 的 estimateDisp() 在无重复时把离散度设为 NA。edgeR 手册给出的做法是手动指定一个 BCV（离散度的平方根）：控制良好的人类样本取 0.4，近交系模式生物取 0.1，技术重复取 0.01。本页实测同一对样本，BCV 取 0.1、0.2、0.4 时 FDR < 0.05 的基因分别是 1645、346 和 17 个。差异基因数完全由这个假设值决定，因此只能当作候选列表。
无重复时可以做的事：按 log2FC 排序列出候选基因、画 MA 图、对排序结果做 GSEA 作为线索，并用 qPCR 或补测重复验证。不能做的事：报告 padj、宣称某基因显著差异、把这份列表当作已确认结果用于后续 ORA 和网络分析。
library(edgeR)
y <- DGEList(counts = cts[, c("ctrl_1", "trt_1")], group = c("ctrl", "trt"))
y <- normLibSizes(y)
bcv <- 0.4 # edgeR 手册：人类样本 0.4，近交模式生物 0.1，技术重复 0.01
et <- exactTest(y, dispersion = bcv^2)
topTags(et, n = 50) # 只作为候选列表，交给 qPCR 或后续实验验证DESeq2
输入必须是原始整数 counts。featureCounts 的矩阵直接读取；Salmon 的转录本定量经 tximport 汇总到基因水平后输入，不要自己把 TPM 加总。
library(DESeq2)
## 方式一：featureCounts 输出（counts.txt，前 6 列是注释）
fc <- read.delim("counts.txt", comment.char = "#", check.names = FALSE)
cts <- as.matrix(fc[, 7:ncol(fc)])
rownames(cts) <- sub("\\.\\d+$", "", fc$Geneid) # 去掉 ENSG00000000003.15 的版本号
colnames(cts) <- sub("\\.bam$", "", basename(colnames(cts)))
coldata <- read.csv("samples.csv", row.names = 1) # 行名 = 样本名，列：condition、batch 等
coldata$condition <- factor(coldata$condition, levels = c("ctrl", "trt")) # 第一个水平是对照
coldata$batch <- factor(coldata$batch)
stopifnot(identical(rownames(coldata), colnames(cts))) # 顺序不一致时结果全错且不报错
dds <- DESeqDataSetFromMatrix(cts, coldata, design = ~ batch + condition)
## 方式二：Salmon + tximport（gene 水平，传入估计 counts 与转录本长度偏移）
library(tximport)
files <- file.path("salmon", rownames(coldata), "quant.sf"); names(files) <- rownames(coldata)
tx2gene <- read.csv("tx2gene.csv") # 两列：TXNAME, GENEID，与注释版本一致
txi <- tximport(files, type = "salmon", tx2gene = tx2gene, ignoreTxVersion = TRUE)
dds <- DESeqDataSetFromTximport(txi, coldata, design = ~ batch + condition)library(DESeq2); library(airway)
data(airway)
airway$dex <- relevel(airway$dex, ref = "untrt")
dds <- DESeqDataSet(airway, design = ~ cell + dex) # cell：配对的细胞系，放在前面
## 预过滤：至少“最小组样本数”个样本 count >= 10（官方 vignette 写法）
smallestGroupSize <- 4
keep <- rowSums(counts(dds) >= 10) >= smallestGroupSize
dds <- dds[keep, ] # 63,677 -> 16,139
dds <- DESeq(dds) # 本页实测 3.7 s
resultsNames(dds) # 查看可用的 coef 名
res <- results(dds, contrast = c("dex", "trt", "untrt"), alpha = 0.05)
summary(res)
## 收缩 LFC：apeglm 只能用 coef
resLFC <- lfcShrink(dds, coef = "dex_trt_vs_untrt", type = "apeglm")
deg <- subset(as.data.frame(res), padj < 0.05 & abs(log2FoldChange) > 1)
nrow(deg) # 951（上调 490，下调 461）
write.csv(as.data.frame(res), "deseq2_all_genes.csv") # 保存全部基因，GSEA 和 universe 要用- coldata 的行顺序必须与 counts 的列顺序一致；不一致时 DESeq2 不报错，结果整体错位。
- featureCounts 的 Geneid 带版本号（ENSG00000000003.15），后续 ID 转换前要去掉版本号，否则 bitr 和 mapIds 全部匹配失败。
- tximport 的 tx2gene 要与 Salmon 索引使用的注释版本一致；转录本名带版本号时设 ignoreTxVersion = TRUE。
- 预过滤只为减少内存和计算时间：本页 63,677 行过滤到 16,139 行后，DESeq() 从 11.8 s 降到 3.7 s，差异基因数基本不变。
- results() 默认 alpha = 0.1。最终阈值用 0.05 时写 alpha = 0.05，独立过滤会按这个值优化。
- 在已经运行过 DESeq() 的对象上删样本或过滤基因后重跑，DESeq2 会提示 “using pre-existing size factors” 并沿用旧值。本页这样做得到 946 个差异基因，重新构建对象得到 951 个。删样本后从 DESeqDataSetFromMatrix 重新开始。
design 与 contrast
design 写错是结果不可信的最常见原因，而且不会产生任何报错。
本页实测：airway 的 4 个细胞系各有处理和未处理样本。design 写 ~ cell + dex 时 padj < 0.05 的基因有 4081 个，只写 ~ dex 时只有 2773 个，少了 32%（加上 |log2FC| > 1 后为 951 对 785）。配对因素不写进模型，细胞系之间的差异会计入残差，检验功效下降。
批次的正确处理是放进 design。不要先用 ComBat 或 limma::removeBatchEffect 改写 counts 再输入 DESeq2；removeBatchEffect 的输出只用于 PCA 和热图展示。批次未知时可以用 sva 估计替代变量，作为协变量加入 design。
## 两组：分子是处理组，分母是对照组
results(dds, contrast = c("condition", "trt", "ctrl"))
## 三组以上：任意两组比较
results(dds, contrast = c("condition", "drugB", "drugA"))
## 带交互的设计改用合并因子，避免手算交互项
dds$group <- factor(paste0(dds$genotype, "_", dds$condition))
design(dds) <- ~ batch + group
dds <- DESeq(dds)
results(dds, contrast = c("group", "KO_trt", "KO_ctrl"))
## 需要 apeglm 收缩的比较不是现成 coef 时：先 relevel 再重跑
dds$condition <- relevel(dds$condition, ref = "drugA")
dds <- nbinomWaldTest(dds) # 不必重估离散度
lfcShrink(dds, coef = "condition_drugB_vs_drugA", type = "apeglm")
## 或者直接用 ashr（支持 contrast）
lfcShrink(dds, contrast = c("condition", "drugB", "drugA"), type = "ashr")| 场景 | design | 要点 |
|---|---|---|
| 两组，无批次 | ~ condition | condition 的第一个水平是对照，用 relevel() 或 factor(levels=) 设置 |
| 有批次或配对（同一病人/细胞系前后） | ~ batch + condition | 感兴趣的变量放最后；batch 必须是因子，写成数字会被当成连续变量 |
| 多组比较 | ~ batch + condition | 用 contrast = c("condition", "B", "A") 取任意两组，结果是 B 相对 A |
| 两因素且关心交互 | ~ batch + group（group 为两个因素合并） | 官方推荐的写法，避免直接解读交互项系数 |
| 批次与分组完全重合（每批只有一组） | 无法校正 | 模型矩阵不满秩，DESeq2 报错；这是实验设计问题，分析阶段无法补救 |
padj 为 NA
padj 为 NA 是 DESeq2 有意设置的，删除这些行之前要知道它们属于哪一类。
预过滤后再跑，独立过滤只去掉 313 个基因。Cook 距离的离群检测要求每组至少 3 个重复；每组 2 个重复时该机制不起作用。
需要不经过滤的结果（例如做 GSEA 想保留所有基因）时，可以用 results(dds, independentFiltering = FALSE)。GSEA 的排序值用 stat 列，它不受独立过滤影响，被独立过滤的基因 stat 仍有值。
| 原因 | 表现 | airway 实测（不预过滤） | 处理 |
|---|---|---|---|
| 所有样本 count 都是 0 | baseMean = 0，log2FC、pvalue、padj 全为 NA | 30,208 / 63,677 行 | 正常，预过滤会去掉 |
| 含 Cook 距离判定的极端离群值 | pvalue 和 padj 为 NA，baseMean 不为 0 | 0 | 检查是否有坏样本；每组 ≥ 7 个重复时 DESeq2 自动替换离群值并重新拟合 |
| 被独立过滤（平均表达低） | pvalue 有值，padj 为 NA | 16,687 个（阈值 baseMean 7.16） | 这些基因检验功效太低，过滤后能提高其余基因的检出数 |
阈值
阈值决定差异基因数，进而决定后续 ORA 的结果，方法部分必须写明。
padj 是 Benjamini-Hochberg 校正后的值，padj < 0.05 表示在这批被判为显著的基因中，假发现的期望比例约为 5%。它控制的是期望比例，不保证每次分析都不超过 5%。
只用 pvalue < 0.05 的问题：本页检验了约 15,800 个基因，即使没有任何真实差异，也期望出现约 790 个 pvalue < 0.05 的基因。审稿人看到未校正 p 值会直接要求重做。pvalue 筛出的基因数太少时，改用 padj < 0.1 或放宽 |log2FC|，并如实写明。
“padj < 0.05 再加 |log2FC| > 1”和“lfcThreshold = 1”是两件事。前者检验 log2FC 是否不为 0，再按点估计筛选；点估计刚过 1 的基因，真实效应可能小于 1。需要声明“表达变化超过 2 倍”时用 lfcThreshold 检验。
| 筛选条件 | airway 实测基因数 | 说明 |
|---|---|---|
| padj < 0.1（results 默认 alpha） | 4889 | DESeq2 的默认值，大多数论文不用 |
| padj < 0.05 | 4081 | 只控制 FDR，不管效应大小 |
| padj < 0.05 且 |log2FC| > 0.585（1.5 倍） | 2025 | 效应小但样本多的数据常用 |
| padj < 0.05 且 |log2FC| > 1（2 倍） | 951 | 最常见的组合 |
| results(lfcThreshold = 1) 后 padj < 0.05 | 241 | 检验的是 |log2FC| 显著大于 1，更严格 |
| pvalue < 0.05（不做多重校正） | 5488 | 比 padj < 0.05 多 34% |
lfcShrink
lfcShrink 只改变 log2FoldChange 列，不改变 pvalue 和 padj。
用在排序、作图和报告效应量
低表达基因的原始 log2FC 噪声很大，收缩后更接近可信值。火山图、MA 图、按 LFC 排序的候选基因表、按 LFC 排序的 GSEA 都应使用收缩后的 LFC。
收缩后差异基因会变少
本页实测 |log2FC| > 1 的基因从 1091 个降到 769 个；用收缩 LFC 加 padj < 0.05 筛选得到 749 个，比原始 LFC 的 951 个少。论文中写明 log2FC 是否经过收缩。
apeglm 只接受 coef
用 contrast 调用会报错：“type='apeglm' shrinkage only for use with 'coef'”。目标比较不在 resultsNames(dds) 里时，先 relevel 参考水平再运行 nbinomWaldTest，或改用 type = "ashr"（支持 contrast，需要安装 ashr 包）。
type = "normal" 不推荐用于大 LFC
DESeq2 vignette 的对比表中，normal 不能保留大 LFC 的幅度，也不能收缩交互项；apeglm 是 1.28 起的默认推荐。本页实测 apeglm 3.6 s，ashr 1.4 s。
ID 转换
clusterProfiler 的 KEGG 和多数 GO 分析使用 Entrez ID，ENSEMBL 或 SYMBOL 都要先转换。
本页实测：16,139 个 Ensembl 基因用 bitr 转 ENTREZID，13.5% 没有对应 ID，180 个 Ensembl ID 对应多个 Entrez ID。951 个差异基因中丢失 63 个，按 airway 的 gene_biotype 统计，主要是 lincRNA（23）、pseudogene（15）和 antisense（13），蛋白编码基因只丢 5 个。丢失的基因大多本来就没有 GO/KEGG 注释，对富集结果影响很小，但要在方法部分写明丢失比例。
bitr 遇到一对多时会返回多行，直接拿去做 ORA 会重复计数。用 mapIds(multiVals = "first") 每个基因只取一个 ID，再 unique()。
- 先去掉 Ensembl 版本号（.15 之类），否则匹配率接近 0。
- SYMBOL 会随 HGNC 改名变化，跨年份数据合并时优先用 Ensembl 或 Entrez ID。
- 小鼠数据用 org.Mm.eg.db，enrichKEGG 的 organism 写 "mmu"；物种和注释包不一致时报错 “No gene can be mapped”。
- “No gene can be mapped” 的其他常见原因：输入是 SYMBOL 而 keyType 是默认的 kegg（即 Entrez）；输入是 data.frame 而不是字符向量；阈值过严导致列表为空。
ORA
ORA 检验差异基因在某个通路中是否多于随机预期，“随机”从哪个基因池里抽，由 universe 决定。
不设 universe 时，背景是数据库中所有有注释的人类基因。组织中不表达的基因永远不会成为差异基因，却被算进背景，导致在该组织高表达的通路被系统性判为富集。airway 是气道平滑肌细胞，不设 universe 时多出来的通路正是平滑肌和信号转导相关的“常见通路”。universe 用 padj 不为 NA 的基因，也就是实际参与了检验的基因。
结果表的 BgRatio 分母（本页 GO 为 11,707，KEGG 为 5,721）小于你提供的 13,776 个基因，这是正常的：clusterProfiler 会把 universe 与该数据库中有注释的基因取交集。同理，GeneRatio 的分母只计有注释的差异基因，888 个差异基因中只有 433 个有 KEGG 注释。
GO 结果冗余很多，用 simplify(cutoff = 0.7) 合并语义相近的条目，本页从 543 条降到 202 条（耗时 26.4 s）。上调和下调基因分开做 ORA 还是合并做，取决于你要回答的问题；合并时同一通路的上调和下调基因会互相抵消掉方向信息，方向用 GSEA 判断。
library(clusterProfiler); library(org.Hs.eg.db)
## 1. ID 转换：ENSEMBL -> ENTREZID，报告丢失比例
eg <- mapIds(org.Hs.eg.db, keys = rownames(res), keytype = "ENSEMBL",
column = "ENTREZID", multiVals = "first")
mean(is.na(eg)) # airway：13.5% 没有 Entrez ID
sig <- rownames(res)[which(res$padj < 0.05 & abs(res$log2FoldChange) > 1)]
tested <- rownames(res)[!is.na(res$padj)] # 真正参与检验的基因
sigE <- unique(na.omit(eg[sig])) # 888
uniE <- unique(na.omit(eg[tested])) # 13,776，作为背景
## 2. GO
ego <- enrichGO(sigE, OrgDb = org.Hs.eg.db, ont = "BP", universe = uniE,
pvalueCutoff = 0.05, qvalueCutoff = 0.2, readable = TRUE)
ego2 <- simplify(ego, cutoff = 0.7) # 去冗余：543 -> 202 条
## 3. KEGG（在线读取 rest.kegg.jp）
ekk <- enrichKEGG(sigE, organism = "hsa", universe = uniE, pvalueCutoff = 0.05)
ekk <- setReadable(ekk, OrgDb = org.Hs.eg.db, keyType = "ENTREZID")
saveRDS(ekk, "enrichKEGG_result.rds") # 保存对象，KEGG 每周更新
readLines("https://rest.kegg.jp/info/kegg")[2] # 记录 KEGG 数据版本，写进方法部分| 分析 | 设 universe（检测到的基因） | 不设 universe（全基因组） | 实测说明 |
|---|---|---|---|
| enrichKEGG | 17 条显著 | 45 条显著 | 多出的 32 条包括 MAPK、Wnt、Hippo、FoxO 等通路 |
| enrichGO BP | 543 条显著 | 555 条显著 | 两者重叠 426 条，各有 100 多条独有，排在前面的条目也不同 |
| enrichKEGG + enrichment_force_universe = TRUE | 174 条显著 | — | 把无注释基因也算进背景，显著性被夸大，不要使用 |
KEGG 联网
enrichKEGG 每次在新的 R 会话中都会实时下载 KEGG 数据，网络不通时整个分析卡住。
clusterProfiler 4.18 的 enrichKEGG 默认 use_internal_data = FALSE，会读取 https://rest.kegg.jp/link/hsa/pathway 和 https://rest.kegg.jp/list/pathway/hsa 两个地址。本页实测首次调用 9–10 s，同一会话内第二次调用使用内存缓存，约 2–3 s；重启 R 后重新下载。
KEGG 数据每周都会更新，本页查询时 rest.kegg.jp/info/kegg 显示的版本是 2026/10/09。半年后重跑同一脚本，通路基因数和显著通路可能变化。保存 enrichKEGG 结果对象，并在方法中写明 KEGG 数据日期。
网络不稳定或要求结果可复现时，用 createKEGGdb 把当前 KEGG 数据做成本地 KEGG.db 包。本页实测 createKEGGdb 0.0.5 生成人类数据用时 7–9 s，安装包约 300 KB；use_internal_data = TRUE 的结果与在线结果完全一致（17 条显著通路逐条相同），Description 列完整。Bioconductor 上的同名 KEGG.db 包自 2012 年后未更新，Bioconductor 3.12 起已弃用，不要从旧版本仓库下载它。
## 一次性：生成并安装当前 KEGG 数据的本地包（需要能访问 GitHub 和 rest.kegg.jp）
remotes::install_github("YuLab-SMU/createKEGGdb")
createKEGGdb::create_kegg_db("hsa") # 本页实测 7–9 s，生成 KEGG.db_1.0.tar.gz（约 300 KB）
install.packages("KEGG.db_1.0.tar.gz", repos = NULL, type = "source")
## 之后离线使用
library(KEGG.db)
ekk_local <- enrichKEGG(sigE, organism = "hsa", universe = uniE, use_internal_data = TRUE)
head(ekk_local$Description) # 有通路名才能正常作图GSEA
GSEA 输入全部检验过的基因的排序，不做差异基因筛选；排序指标的选择对结果的影响比任何其他参数都大。
library(clusterProfiler); library(msigdbr)
## 排序向量：全部检验过的基因，按 Wald stat 从大到小，名字是 Entrez ID
rk <- res$stat; names(rk) <- eg[rownames(res)]
rk <- rk[!is.na(rk) & !is.na(names(rk))]
rk <- sort(rk[!duplicated(names(rk))], decreasing = TRUE) # 13,952 个基因
## MSigDB hallmark（msigdbr 10 起参数是 collection，列名是 ncbi_gene）
h <- msigdbr(species = "Homo sapiens", collection = "H")
t2g <- data.frame(term = h$gs_name, gene = as.character(h$ncbi_gene))
gs <- GSEA(rk, TERM2GENE = t2g, minGSSize = 15, maxGSSize = 500,
pvalueCutoff = 0.05, eps = 0, seed = TRUE) # seed = TRUE 让结果可复现
gg <- gseGO(rk, OrgDb = org.Hs.eg.db, ont = "BP", minGSSize = 15, maxGSSize = 500,
pvalueCutoff = 0.05, eps = 0, seed = TRUE) # 本页实测 11.8 s，198 条
## 不要再写 nPerm = 1000：当前版本基于 fgsea 多层抽样，会给出不推荐 nPerm 的警告NES 与 FDR 的阈值
R 包结果看 p.adjust（BH 校正）或 qvalue，常用 < 0.05。Broad GSEA 软件的 FDR q-value 是另一种基于置换的估计，官方建议的 25% 用于产生假设；用 R 包时写“FDR < 0.25”需要说明这一出处。NES 的正负表示基因集集中在排序列表的顶端还是底端，|NES| 只用于比较，没有通用的显著阈值，|NES| > 1 的说法没有统计依据。
leading edge 的含义
leading edge 是富集分数达到峰值之前出现的那部分基因集成员，是驱动富集信号的核心基因。clusterProfiler 中基因列表在 core_enrichment 列；leading_edge 列的 tags 是核心基因占基因集的比例，list 是峰值在排序列表中的位置，signal 是二者的综合。本页 ADIPOGENESIS 的 189 个基因中有 62 个在 leading edge（tags = 33%）。挑选验证基因时从 core_enrichment 里选。
结果会随随机种子变化
不设 seed = TRUE、调用前分别 set.seed(7) 和 set.seed(8) 时，同一数据的 gseGO 得到 210 和 216 条显著条目；seed = TRUE 时在新会话中重复运行都是 198 条（本页实测）。输入的微小变化同样会改变边界结果：本页 DESeq2 沿用旧 size factor（差异基因 946 个）与重新构建对象（951 个）相比，gseGO 显著条目从 230 条变为 198 条。论文结果要固定 seed，并把 p.adjust 在 0.03–0.07 之间的条目当作不稳定结果。
msigdbr 新版参数
msigdbr 10.0 起 category 参数改为 collection，Entrez 列由 entrez_gene 改为 ncbi_gene，旧代码会给出弃用警告。当前数据版本为 2026.1.Hs；首次使用需要下载数据，本页实测 35 s。小鼠数据用 db_species = "MM" 取原生小鼠基因集。
| 排序指标 | hallmark 显著数（p.adjust < 0.05） | 特点 |
|---|---|---|
| DESeq2 的 stat（Wald 统计量） | 15（上调 11，下调 4） | 同时反映效应大小和可信度，有正负号，不受独立过滤影响；推荐默认使用 |
| lfcShrink（apeglm）后的 log2FC | 3 | 只反映效应大小；未收缩的 LFC 会把低表达噪声基因排到两端，不要使用 |
| sign(log2FC) × -log10(pvalue) | 1 | 极显著基因的值可达上百，少数基因主导富集分数；p 值为 0 时产生 Inf |
| 只取差异基因再排序 | gseGO 只剩 2 条（全基因为 198 条） | 破坏了 GSEA 的前提，等于没有背景 |
结果不一致
两者核心算法相同，结果不同几乎都来自输入和默认参数，逐项对齐即可复现。
网上有回答称 fgsea 不能计算加权统计量，这已过时：fgsea 的 gseaParam 和 clusterProfiler 的 exponent 默认都是 1，与 Broad 的 weighted 评分一致。
| 差异来源 | GSEA 桌面版（Broad） | clusterProfiler / fgsea | 对齐方法 |
|---|---|---|---|
| 输入与排序 | 标准 GSEA 输入表达矩阵，用 Signal2Noise 等指标排序 | 输入你提供的排序向量（如 Wald stat） | 桌面版改用 GSEAPreranked，导入同一个 .rnk 文件 |
| 置换方式 | 标准 GSEA 默认置换样本标签（phenotype），样本少时建议 gene_set；Preranked 只能 gene_set | 只置换基因（fgsea 多层抽样） | 比较时都用 Preranked |
| 基因集大小 | min 15，max 500 | clusterProfiler 默认 minGSSize 10；fgsea 单独使用时 minSize 默认 1 | R 中显式写 minGSSize = 15, maxGSSize = 500 |
| 基因 ID | Preranked 默认 Remap_Only，按芯片注释重映射符号 | 不做映射，名字必须与基因集一致 | 统一用 SYMBOL 或 Entrez，并去掉重复 |
| 并列值 | 要求排序值无重复，并列时顺序任意 | 给出 ties 警告，按任意顺序处理 | 排序前处理并列值，或接受微小差异 |
| 随机种子 | 默认用时间戳，每次不同 | 不设 seed 时每次不同 | 两边都固定种子 |
| 显著性指标 | FDR q-value（基于置换的 NES 分布） | BH 校正的 p.adjust | 比较 NES 方向和排名，不直接比较 FDR 数值 |
| 基因集版本 | 下载的 gmt 文件 | msigdbr 当前版本 | 使用同一版本的 gmt |
作图
作图的阈值、纵轴和颜色必须与正文的筛选条件一致。
library(ggplot2); library(ggrepel); library(enrichplot)
## 火山图：y 轴用 padj，阈值线与筛选条件一致
df <- as.data.frame(res)
df$symbol <- mapIds(org.Hs.eg.db, rownames(df), keytype = "ENSEMBL", column = "SYMBOL")
df <- df[!is.na(df$padj), ]
df$grp <- with(df, ifelse(padj < 0.05 & log2FoldChange > 1, "up",
ifelse(padj < 0.05 & log2FoldChange < -1, "down", "ns")))
lab <- head(df[df$grp != "ns", ][order(df$padj[df$grp != "ns"]), ], 15)
ggplot(df, aes(log2FoldChange, -log10(padj), colour = grp)) +
geom_point(size = 0.6, alpha = 0.6) +
scale_colour_manual(values = c(up = "#c0392b", down = "#2471a3", ns = "grey70")) +
geom_vline(xintercept = c(-1, 1), linetype = 2) +
geom_hline(yintercept = -log10(0.05), linetype = 2) +
geom_text_repel(data = lab, aes(label = symbol), colour = "black", size = 3) +
theme_bw()
## 气泡图：去冗余后的 GO
dotplot(ego2, showCategory = 15)
## GSEA 曲线图：pvalue_table 需要安装 gridExtra，宽度给到 10 英寸以上表格才不被截断
p <- gseaplot2(gs, geneSetID = 1:2, pvalue_table = TRUE)
ggsave("gsea.png", p, width = 10, height = 6, dpi = 300)火山图
纵轴用 -log10(padj)，阈值线与正文一致（本页为 padj 0.05、|log2FC| 1）。padj 为 NA 的基因要先去掉，否则 ggplot 给出缺失值警告并丢点。极显著基因的 padj 可能下溢为 0，-log10 后是 Inf，需要把 0 替换为最小非零值并在图注说明。标注基因按 padj 或预先指定的基因挑选。
气泡图
dotplot 默认横轴是 GeneRatio，点大小是基因数，颜色是 p.adjust。GeneRatio 的分母是有注释的差异基因数，不是全部差异基因数。展示去冗余后的条目，避免前 15 条全是同一个过程的不同层级。
GSEA 曲线图
上半部分的峰值是富集分数，峰值出现在左侧且为正值表示基因集在处理组上调（排序方向由你的 contrast 决定）。中间竖线是基因集成员在排序中的位置，峰值之前的成员就是 leading edge。gseaplot2 的 pvalue_table 需要安装 gridExtra，否则报错 “The package "gridExtra" is required for tableGrob2()”。
- 火山图纵轴用 pvalue，筛选却用 padj，图中“显著”的点与正文列表对不上。
- 富集图只展示上调基因的 ORA，却在正文里写通路“被激活”；方向要用 GSEA 的 NES 或分别对上下调做 ORA 后说明。
- 气泡图不写背景基因集、数据库版本和校正方法。
- 把 GSEA 中 p.adjust > 0.05 的通路曲线放进主图，并写成“显著富集”。
- 差异基因数量与补充表的行数不一致，常见原因是 SYMBOL 去重或 ID 转换丢失后没有更新数字。
国内经验
以下来自 CSDN、简书的作者实测帖，已用 clusterProfiler 4.18.4 和 createKEGGdb 0.0.5 核对。
KEGG 本地化：两位作者独立验证可用
CSDN 作者 hutong_oncology 和 qq_41547057 在 2023 年都遇到在线 enrichKEGG 连接失败或突然 “No gene can be mapped”，改用 createKEGGdb::create_kegg_db 生成本地 KEGG.db 后解决。前者还提到本地库能避免 KEGG 在线更新造成前后结果不一致。本页实测该流程在当前版本仍可用，结果与在线一致。
旧版 KEGG.db 缺 Description 导致作图报错
hutong_oncology 记录：部分人生成的 KEGG.db 缺少通路描述，富集表能出结果，但 barplot、dotplot 报错 “Error in ans[ypos] <- rep(yes, length.out = len)[ypos]” 和 “'x' is NULL so the result will be NULL”。很多帖子把它归因于阈值太严，实际原因是缺 Description。本页用 createKEGGdb 0.0.5 生成的包 Description 完整；遇到此报错先运行 head(ekk$Description) 检查。
R.utils::setOption 改下载方式的办法已失效
旧帖常用 R.utils::setOption("clusterProfiler.download.method", "auto") 解决 KEGG 下载失败。clusterProfiler 4.18.4 的 KEGG 下载改由 yulab.utils::yread 调用 readLines 完成，源码中已不读取这个选项。网络问题用系统代理或本地 KEGG.db 解决。
“padj < 0.05 保证假阳性不超过 5%”的说法不准确
简书等教程常这样解释 padj。BH 校正控制的是假发现的期望比例，单次分析中实际比例可能高于 5%。写方法时用“FDR 控制在 5%”即可，不要写“保证”。
交给 Agent
可以用一句话描述完整的差异表达与富集分析任务。
指令示例：“用 featureCounts 的 counts.txt 和 samples.csv 做 DESeq2 差异分析，design 为 ~ batch + condition，比较 trt 相对 ctrl；阈值 padj < 0.05 且 |log2FC| > 1；用 apeglm 收缩画火山图；以检测基因为背景做 GO BP 与 KEGG 富集，KEGG 数据本地化并记录版本；用 Wald stat 排序做 hallmark 与 GO BP 的 GSEA，固定种子；输出差异基因表、富集表和三张图。”
智能体在隔离云电脑中安装 R 与 Bioconductor 包，按本页流程运行：检查样本表与 counts 列顺序、预过滤、DESeq2、lfcShrink、ID 转换并报告丢失比例、带 universe 的 ORA、GSEA 和作图。工作区中保留脚本、sessionInfo、KEGG 数据日期、结果表和图，可以复现。智能体会对结果做对抗审阅，例如核对火山图阈值与筛选条件是否一致、universe 是否设置、GSEA 是否固定了种子。
你仍需要自己核对：分组和批次信息是否与实验记录一致，design 是否反映了真实的配对结构，阈值选择是否符合领域惯例，以及富集到的通路在生物学上是否说得通。
参考资料
- DESeq2 vignette（Bioconductor release） — 预过滤写法、padj 为 NA 的三类原因、lfcShrink 方法对比表、无重复 FAQ、Cook 距离与 7 个重复的替换规则
- edgeR User's Guide 2.13 What to do if you have no replicates — 无重复时的四种选择与 BCV 参考值 0.4/0.1/0.01
- Li Y, et al. Exaggerated false positives by popular differential expression methods when analyzing human population samples. Genome Biology 2022 — 大样本人群数据中 DESeq2/edgeR 的 FDR 膨胀与 Wilcoxon 建议
- clusterProfiler 书：KEGG analysis — enrichKEGG 在线数据与 KEGG.db 自 2012 年未更新的说明；书中 gseKEGG 示例仍使用已不推荐的 nPerm
- YuLab-SMU/createKEGGdb（GitHub） — 本地 KEGG.db 生成工具，本页实测 0.0.5
- Bioconductor support：BgRatio N value does not equal set gene universe size — universe 与注释基因取交集的行为
- GSEA User Guide（GSEA-MSigDB 文档） — FDR q-value 25% 阈值的含义、phenotype 与 gene_set 置换的选择、leading edge 定义
- GSEAPreranked 模块文档（v7） — Preranked 只做 gene_set 置换、默认 min 15/max 500、weighted 评分、随机种子默认时间戳、排序值不能重复
- fgsea（Bioconductor） — 多层抽样 GSEA 实现，clusterProfiler GSEA 的默认后端
- msigdbr（CRAN） — 10.0 起 collection 参数与 ncbi_gene 列
- Biostars：GSEA on preranked list with weighted enrichment statistic — 经验帖：fgsea 与 GSEA 软件结果不一致的讨论，其中“fgsea 不能加权”的说法已过时
- CSDN hutong_oncology：创建 KEGG.db 过程中的报错及解决办法 — 经验帖：createKEGGdb 流程、缺 Description 导致的作图报错原文
- CSDN qq_41547057：clusterProfiler 做 KEGG 富集时出现的错误及解决方法 — 经验帖：在线 enrichKEGG 报 No gene can be mapped 后改用本地 KEGG.db
- 简书：为什么选择 log2FC 和 padj 筛选差异基因 — 经验帖：常用阈值的中文解释，其中 FDR 的表述需修正
- KEGG REST API：info/kegg — 查询 KEGG 数据版本日期