ScientifyChoosing a method
For a standard two-group comparison the three methods overlap heavily. A wrong choice costs you mainly at extreme sample sizes and with the wrong data type.
Measured on this page (R 4.5.3, DESeq2 1.50.2, edgeR 4.10.5, limma 3.68.5, clusterProfiler 4.18.4; Apple M2, 8 cores, 16 GB): DESeq2 found 951 DE genes, edgeR QL 1183 and limma-voom 1145; the intersection was 910 and the union 1227. The differences sit on borderline genes with effect sizes near 1. State the method and version in the paper; there is no need to intersect the three methods.
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 pipeline (recommended in the user's guide)
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)| Situation | Recommended | Basis and cost of a wrong choice |
|---|---|---|
| RNA-seq counts, 3–6 biological replicates per group | DESeq2 or edgeR (QL pipeline) | Both shrink dispersions across genes, which is more stable at small n; on this page the intersection of the three methods is 74% of their union |
| 2 replicates per group | DESeq2 or edgeR, with conclusions treated as exploratory | It runs, but Cook's distance outlier detection needs at least 3 replicates per group, so an outlier sample cannot be flagged |
| No biological replicates | Descriptive analysis only; edgeR with a fixed BCV for a candidate list | DESeq2 exits with an error; edgeR results change 100-fold with the BCV value (next section) |
| Many samples (20+ per group), multiple factors or continuous covariates | limma-voom | Fast linear-model fitting; duplicateCorrelation handles repeated measures |
| Population cohorts, dozens to hundreds per group | Wilcoxon rank-sum test, or DESeq2/edgeR with SVA/RUV for technical variation | Li et al. (2022) ran permutation tests on 13 population datasets; at a 5% target FDR, the actual FDR of DESeq2/edgeR sometimes exceeded 20% |
| Microarray expression matrix | limma | The data are continuous log-scale values and cannot go into DESeq2/edgeR |
| TPM, FPKM or other normalized matrices | Recover raw counts; if only TPM exists, use limma-trend and say so in the methods | DESeq2/edgeR need integer counts; rounding FPKM distorts dispersion estimates |
No replicates
Without replicates there is no estimate of within-group biological variation, so any p value depends on an assumption you supply.
DESeq2 has not supported designs without replicates since version 1.22. Running DESeq() on one treated and one untreated sample from the same airway cell line gives this error: “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.” Tutorials from before 2018 that run DESeq2 without replicates no longer work.
edgeR's estimateDisp() sets the dispersion to NA without replicates. The edgeR user's guide suggests fixing a BCV (the square root of the dispersion): 0.4 for well-controlled human data, 0.1 for genetically identical model organisms and 0.01 for technical replicates. On the same pair of samples, BCV values of 0.1, 0.2 and 0.4 gave 1645, 346 and 17 genes with FDR < 0.05. The number of DE genes is set entirely by this assumed value, so the result is only a candidate list.
What you can do without replicates: list candidate genes ranked by log2FC, draw an MA plot, run GSEA on the ranking as a lead, and validate by qPCR or by sequencing replicates. What you cannot do: report padj, claim that a gene is significantly differentially expressed, or feed this list into ORA and network analysis as if it were confirmed.
library(edgeR)
y <- DGEList(counts = cts[, c("ctrl_1", "trt_1")], group = c("ctrl", "trt"))
y <- normLibSizes(y)
bcv <- 0.4 # edgeR guide: 0.4 for human data, 0.1 for inbred model organisms, 0.01 for technical replicates
et <- exactTest(y, dispersion = bcv^2)
topTags(et, n = 50) # a candidate list only; validate with qPCR or follow-up experimentsDESeq2
The input must be raw integer counts. Read featureCounts matrices directly; summarize Salmon transcript quantifications to genes with tximport instead of summing TPM yourself.
library(DESeq2)
## Option 1: featureCounts output (counts.txt; the first 6 columns are annotation)
fc <- read.delim("counts.txt", comment.char = "#", check.names = FALSE)
cts <- as.matrix(fc[, 7:ncol(fc)])
rownames(cts) <- sub("\\.\\d+$", "", fc$Geneid) # strip the version suffix from ENSG00000000003.15
colnames(cts) <- sub("\\.bam$", "", basename(colnames(cts)))
coldata <- read.csv("samples.csv", row.names = 1) # row names = sample names; columns: condition, batch, ...
coldata$condition <- factor(coldata$condition, levels = c("ctrl", "trt")) # the first level is the control
coldata$batch <- factor(coldata$batch)
stopifnot(identical(rownames(coldata), colnames(cts))) # a mismatched order gives wrong results without any error
dds <- DESeqDataSetFromMatrix(cts, coldata, design = ~ batch + condition)
## Option 2: Salmon + tximport (gene level; passes estimated counts and a transcript-length offset)
library(tximport)
files <- file.path("salmon", rownames(coldata), "quant.sf"); names(files) <- rownames(coldata)
tx2gene <- read.csv("tx2gene.csv") # two columns: TXNAME, GENEID; same annotation release as the index
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: the paired cell line, placed first
## Pre-filter: count >= 10 in at least smallestGroupSize samples (as in the official vignette)
smallestGroupSize <- 4
keep <- rowSums(counts(dds) >= 10) >= smallestGroupSize
dds <- dds[keep, ] # 63,677 -> 16,139
dds <- DESeq(dds) # measured on this page: 3.7 s
resultsNames(dds) # list the available coef names
res <- results(dds, contrast = c("dex", "trt", "untrt"), alpha = 0.05)
summary(res)
## Shrink LFC: apeglm accepts coef only
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 up, 461 down)
write.csv(as.data.frame(res), "deseq2_all_genes.csv") # keep all genes; GSEA and the universe need them- The row order of coldata must match the column order of the counts; DESeq2 does not report a mismatch and the results are shifted wholesale.
- featureCounts Geneids carry a version (ENSG00000000003.15); strip it before ID conversion, otherwise bitr and mapIds match nothing.
- The tx2gene table for tximport must use the same annotation release as the Salmon index; set ignoreTxVersion = TRUE when transcript names carry versions.
- Pre-filtering only saves memory and time: filtering 63,677 rows to 16,139 cut DESeq() from 11.8 s to 3.7 s here, with essentially the same number of DE genes.
- results() defaults to alpha = 0.1. If your final cutoff is 0.05, pass alpha = 0.05 so independent filtering is optimized for it.
- If you drop samples or filter genes on an object that has already been through DESeq() and rerun it, DESeq2 prints “using pre-existing size factors” and keeps the old values. Doing this here gave 946 DE genes; rebuilding the object gave 951. After dropping samples, start again from DESeqDataSetFromMatrix.
Design and contrast
A wrong design is the most common reason for untrustworthy results, and it produces no error.
Measured on this page: each of the 4 airway cell lines has a treated and an untreated sample. With ~ cell + dex, 4081 genes have padj < 0.05; with ~ dex alone, only 2773, 32% fewer (951 vs 785 after adding |log2FC| > 1). If the pairing factor is left out, cell-line differences go into the residual and power drops.
Handle batch by putting it in the design. Do not rewrite counts with ComBat or limma::removeBatchEffect before DESeq2; use removeBatchEffect output only for PCA and heatmaps. If the batch is unknown, estimate surrogate variables with sva and add them to the design as covariates.
## Two groups: numerator is treatment, denominator is control
results(dds, contrast = c("condition", "trt", "ctrl"))
## Three or more groups: any pairwise comparison
results(dds, contrast = c("condition", "drugB", "drugA"))
## For designs with interactions, combine factors into one group factor instead of reading interaction terms
dds$group <- factor(paste0(dds$genotype, "_", dds$condition))
design(dds) <- ~ batch + group
dds <- DESeq(dds)
results(dds, contrast = c("group", "KO_trt", "KO_ctrl"))
## When the comparison you want to shrink with apeglm is not an existing coef: relevel and refit
dds$condition <- relevel(dds$condition, ref = "drugA")
dds <- nbinomWaldTest(dds) # dispersions do not need to be re-estimated
lfcShrink(dds, coef = "condition_drugB_vs_drugA", type = "apeglm")
## Or use ashr, which supports contrast
lfcShrink(dds, contrast = c("condition", "drugB", "drugA"), type = "ashr")| Scenario | design | Notes |
|---|---|---|
| Two groups, no batch | ~ condition | The first level of condition is the reference; set it with relevel() or factor(levels =) |
| Batch or pairing (same patient or cell line before/after) | ~ batch + condition | Put the variable of interest last; batch must be a factor, as a number it is treated as continuous |
| Multiple groups | ~ batch + condition | Use contrast = c("condition", "B", "A") for any pair; the result is B relative to A |
| Two factors with an interaction of interest | ~ batch + group (group combines the two factors) | The officially recommended approach; avoids interpreting interaction coefficients directly |
| Batch fully confounded with group (one group per batch) | Cannot be corrected | The model matrix is not full rank and DESeq2 errors; this is a design problem that analysis cannot fix |
NA in padj
DESeq2 sets these NAs on purpose. Know which category a row belongs to before you drop it.
After pre-filtering, independent filtering removed only 313 genes. Cook's distance outlier detection needs at least 3 replicates per group and does nothing with 2.
If you need unfiltered results, use results(dds, independentFiltering = FALSE). For GSEA, rank by the stat column, which independent filtering does not touch; filtered genes still have a stat value.
| Reason | Symptom | airway, no pre-filter | What to do |
|---|---|---|---|
| All samples have zero counts | baseMean = 0; log2FC, pvalue and padj are all NA | 30,208 of 63,677 rows | Normal; pre-filtering removes them |
| An extreme outlier flagged by Cook's distance | pvalue and padj are NA, baseMean is not 0 | 0 | Check for a bad sample; with ≥ 7 replicates per group DESeq2 replaces outliers and refits automatically |
| Removed by independent filtering (low mean count) | pvalue present, padj NA | 16,687 (threshold baseMean 7.16) | These genes have too little power; filtering them raises detections among the rest |
Thresholds
The cutoff sets the number of DE genes, which in turn sets the ORA result, so the methods section must state it.
padj is the Benjamini-Hochberg adjusted value. padj < 0.05 means the expected proportion of false discoveries among the genes called significant is about 5%. It controls an expectation and does not guarantee that any single analysis stays below 5%.
The problem with pvalue < 0.05 alone: about 15,800 genes were tested here, so even with no true differences you would expect about 790 genes below 0.05. Reviewers who see unadjusted p values will ask for a redo. If padj gives too few genes, use padj < 0.1 or relax |log2FC|, and report that.
“padj < 0.05 plus |log2FC| > 1” and “lfcThreshold = 1” are different. The first tests whether log2FC differs from 0 and then filters on the point estimate; a gene whose estimate is just above 1 may have a true effect below 1. Use the lfcThreshold test when you claim changes larger than 2-fold.
| Filter | Genes in airway | Note |
|---|---|---|
| padj < 0.1 (results default alpha) | 4889 | DESeq2's default; most papers do not use it |
| padj < 0.05 | 4081 | Controls FDR only, ignores effect size |
| padj < 0.05 and |log2FC| > 0.585 (1.5-fold) | 2025 | Common for small effects with many samples |
| padj < 0.05 and |log2FC| > 1 (2-fold) | 951 | The most common combination |
| padj < 0.05 after results(lfcThreshold = 1) | 241 | Tests whether |log2FC| is significantly above 1; stricter |
| pvalue < 0.05 (no multiple-testing correction) | 5488 | 34% more than padj < 0.05 |
lfcShrink
lfcShrink changes only the log2FoldChange column; pvalue and padj stay the same.
Use it for ranking, plots and reported effect sizes
Raw log2FC of low-count genes is very noisy, and shrinkage pulls it toward credible values. Volcano plots, MA plots, candidate tables sorted by LFC and LFC-ranked GSEA should all use shrunken LFC.
Shrinkage reduces the DE count
On this page, genes with |log2FC| > 1 fell from 1091 to 769; filtering on shrunken LFC with padj < 0.05 gives 749 genes instead of 951 with raw LFC. State in the paper whether log2FC values were shrunken.
apeglm accepts coef only
Calling it with contrast errors: “type='apeglm' shrinkage only for use with 'coef'”. If your comparison is not in resultsNames(dds), relevel the reference and run nbinomWaldTest, or use type = "ashr" (supports contrast; needs the ashr package).
type = "normal" is not recommended for large LFCs
In the DESeq2 vignette's comparison table, normal does not preserve the size of large LFCs and cannot shrink interaction terms; apeglm has been the default recommendation since 1.28. Measured here: apeglm 3.6 s, ashr 1.4 s.
ID conversion
clusterProfiler's KEGG analysis and most GO analyses use Entrez IDs, so ENSEMBL or SYMBOL inputs must be converted first.
Measured on this page: converting 16,139 Ensembl genes to ENTREZID with bitr left 13.5% without an ID, and 180 Ensembl IDs mapped to more than one Entrez ID. Of the 951 DE genes, 63 were lost; by airway's gene_biotype they were mostly lincRNA (23), pseudogene (15) and antisense (13), with only 5 protein-coding genes. Most lost genes have no GO/KEGG annotation anyway, so the effect on enrichment is small, but report the loss rate in the methods.
bitr returns several rows for one-to-many mappings, which double-counts genes in ORA. Use mapIds(multiVals = "first") to keep one ID per gene, then unique().
- Strip Ensembl version suffixes (such as .15) first, or the match rate is close to 0.
- SYMBOLs change with HGNC renaming; when merging data across years, prefer Ensembl or Entrez IDs.
- For mouse data use org.Mm.eg.db and organism = "mmu" in enrichKEGG; a species/annotation mismatch gives “No gene can be mapped”.
- Other common causes of “No gene can be mapped”: SYMBOL input with the default keyType kegg (Entrez for human); a data.frame instead of a character vector; an empty list after a strict cutoff.
ORA
ORA tests whether DE genes appear in a pathway more often than chance. The universe defines the pool that chance draws from.
Without a universe, the background is every annotated human gene in the database. Genes not expressed in the tissue can never become DE genes, yet they count toward the background, so pathways highly expressed in that tissue are systematically called enriched. airway uses airway smooth muscle cells, and the extra pathways without a universe are exactly the smooth-muscle and signaling “usual suspects”. Use the genes with non-NA padj, that is, the genes that were actually tested.
The BgRatio denominator in the results (11,707 for GO and 5,721 for KEGG here) is smaller than the 13,776 genes you supplied. This is expected: clusterProfiler intersects the universe with genes annotated in that database. Likewise, the GeneRatio denominator counts only annotated DE genes; just 433 of the 888 DE genes have KEGG annotation.
GO results are highly redundant. simplify(cutoff = 0.7) merges semantically similar terms, from 543 to 202 here (26.4 s). Whether to run ORA separately on up- and down-regulated genes depends on your question; pooling them loses direction within a pathway, so judge direction with GSEA.
library(clusterProfiler); library(org.Hs.eg.db)
## 1. ID conversion: ENSEMBL -> ENTREZID; report the loss rate
eg <- mapIds(org.Hs.eg.db, keys = rownames(res), keytype = "ENSEMBL",
column = "ENTREZID", multiVals = "first")
mean(is.na(eg)) # airway: 13.5% have no Entrez ID
sig <- rownames(res)[which(res$padj < 0.05 & abs(res$log2FoldChange) > 1)]
tested <- rownames(res)[!is.na(res$padj)] # genes that were actually tested
sigE <- unique(na.omit(eg[sig])) # 888
uniE <- unique(na.omit(eg[tested])) # 13,776, used as the background
## 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) # remove redundancy: 543 -> 202 terms
## 3. KEGG (reads rest.kegg.jp online)
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") # save the object; KEGG updates weekly
readLines("https://rest.kegg.jp/info/kegg")[2] # record the KEGG data release for the methods section| Analysis | Universe = detected genes | No universe (whole genome) | Measured note |
|---|---|---|---|
| enrichKEGG | 17 significant | 45 significant | The extra 32 include MAPK, Wnt, Hippo and FoxO pathways |
| enrichGO BP | 543 significant | 555 significant | 426 terms overlap; each has 100+ unique terms and the top terms differ |
| enrichKEGG + enrichment_force_universe = TRUE | 174 significant | — | Counts unannotated genes in the background and inflates significance; do not use |
KEGG access
enrichKEGG downloads KEGG data live in every new R session, so the whole analysis stalls when the network fails.
In clusterProfiler 4.18, enrichKEGG defaults to use_internal_data = FALSE and reads https://rest.kegg.jp/link/hsa/pathway and https://rest.kegg.jp/list/pathway/hsa. The first call took 9–10 s here; a second call in the same session uses an in-memory cache and took about 2–3 s. After restarting R it downloads again.
KEGG data are updated weekly; rest.kegg.jp/info/kegg showed release 2026/10/09 when this page was checked. Rerunning the same script six months later may change pathway gene counts and the significant pathways. Save the enrichKEGG result object and report the KEGG data date in the methods.
When the network is unreliable or you need reproducible results, use createKEGGdb to package current KEGG data as a local KEGG.db. createKEGGdb 0.0.5 built the human package in 7–9 s here, about 300 KB; results with use_internal_data = TRUE were identical to the online run (the same 17 significant pathways) with complete Description values. The KEGG.db package on Bioconductor has not been updated since 2012 and was deprecated in Bioconductor 3.12; do not download it from old repositories.
## Once: build and install a local package of current KEGG data (needs access to GitHub and rest.kegg.jp)
remotes::install_github("YuLab-SMU/createKEGGdb")
createKEGGdb::create_kegg_db("hsa") # 7–9 s on this page; produces KEGG.db_1.0.tar.gz (about 300 KB)
install.packages("KEGG.db_1.0.tar.gz", repos = NULL, type = "source")
## Offline use afterwards
library(KEGG.db)
ekk_local <- enrichKEGG(sigE, organism = "hsa", universe = uniE, use_internal_data = TRUE)
head(ekk_local$Description) # plots need pathway namesGSEA
GSEA takes a ranking of all tested genes with no DE cutoff. The choice of ranking metric affects the result more than any other parameter.
library(clusterProfiler); library(msigdbr)
## Ranked vector: all tested genes, sorted by Wald stat in decreasing order, named by 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 genes
## MSigDB hallmark (since msigdbr 10 the argument is collection and the column is 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 makes results reproducible
gg <- gseGO(rk, OrgDb = org.Hs.eg.db, ont = "BP", minGSSize = 15, maxGSSize = 500,
pvalueCutoff = 0.05, eps = 0, seed = TRUE) # 11.8 s on this page, 198 terms
## Do not pass nPerm = 1000: current versions use fgsea multilevel sampling and warn against nPermNES and FDR cutoffs
For R packages, use p.adjust (BH) or qvalue, usually < 0.05. The FDR q-value in Broad's GSEA software is a different permutation-based estimate, and the 25% cutoff it suggests is meant for hypothesis generation; if you write “FDR < 0.25” for R results, cite that origin. The sign of NES says whether the set sits at the top or bottom of the ranking. |NES| is for comparison only and has no general significance cutoff; “|NES| > 1” has no statistical basis.
What the leading edge means
The leading edge is the subset of set members that appear before the enrichment score reaches its peak; these genes drive the signal. In clusterProfiler the gene list is in core_enrichment. In the leading_edge column, tags is the fraction of the set in the core, list is where the peak falls in the ranking, and signal combines the two. Here, 62 of the 189 ADIPOGENESIS genes are in the leading edge (tags = 33%). Pick validation genes from core_enrichment.
Results change with the random seed
Without seed = TRUE, calling set.seed(7) and set.seed(8) beforehand gave 210 and 216 significant gseGO terms on the same data; with seed = TRUE, repeated runs in fresh sessions all gave 198 (measured on this page). Small input changes also move borderline results: DESeq2 reusing old size factors (946 DE genes) versus a rebuilt object (951) changed significant gseGO terms from 230 to 198. Fix the seed for published results and treat terms with p.adjust between 0.03 and 0.07 as unstable.
New msigdbr arguments
Since msigdbr 10.0, category is replaced by collection and the Entrez column entrez_gene by ncbi_gene; old code triggers a deprecation warning. The current data release is 2026.1.Hs; the first call downloads the data, 35 s here. For mouse data use db_species = "MM" to get native mouse gene sets.
| Ranking metric | Significant hallmark sets (p.adjust < 0.05) | Properties |
|---|---|---|
| DESeq2 stat (Wald statistic) | 15 (11 up, 4 down) | Reflects both effect size and confidence, signed, unaffected by independent filtering; recommended default |
| log2FC after lfcShrink (apeglm) | 3 | Reflects effect size only; unshrunken LFC pushes noisy low-count genes to the ends, do not use it |
| sign(log2FC) × -log10(pvalue) | 1 | Highly significant genes reach values above 100, so a few genes dominate the score; p values of 0 give Inf |
| Ranking DE genes only | gseGO drops to 2 terms (198 with all genes) | Breaks the premise of GSEA; there is no background |
Mismatched results
The core algorithm is the same. Differences almost always come from inputs and defaults, and aligning them one by one reproduces the result.
Some forum answers say fgsea cannot compute a weighted statistic. That is outdated: fgsea's gseaParam and clusterProfiler's exponent both default to 1, which matches Broad's weighted scoring.
| Source of difference | GSEA desktop (Broad) | clusterProfiler / fgsea | How to align |
|---|---|---|---|
| Input and ranking | Standard GSEA takes an expression matrix and ranks by Signal2Noise or similar | Takes the ranked vector you provide (such as the Wald stat) | Use GSEAPreranked in the desktop app with the same .rnk file |
| Permutation | Standard GSEA permutes phenotype labels by default, with gene_set advised for small samples; Preranked only does gene_set | Permutes genes only (fgsea multilevel sampling) | Compare using Preranked on both sides |
| Gene set size | min 15, max 500 | clusterProfiler minGSSize defaults to 10; standalone fgsea minSize defaults to 1 | Set minGSSize = 15, maxGSSize = 500 explicitly in R |
| Gene IDs | Preranked defaults to Remap_Only, remapping symbols via chip annotation | No remapping; names must match the gene sets | Use SYMBOL or Entrez consistently and remove duplicates |
| Ties | Requires no duplicate ranking values; ties are ordered arbitrarily | Warns about ties and orders them arbitrarily | Resolve ties before ranking, or accept small differences |
| Random seed | Defaults to a timestamp, different every run | Different every run without a seed | Fix the seed on both sides |
| Significance measure | FDR q-value (from the permuted NES distribution) | BH-adjusted p.adjust | Compare NES direction and rank, not FDR values |
| Gene set release | The downloaded gmt file | Current msigdbr release | Use the same gmt release |
Plots
Cutoffs, axes and colors in every plot must match the filters stated in the text.
library(ggplot2); library(ggrepel); library(enrichplot)
## Volcano plot: padj on the y axis; threshold lines match the filter
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()
## Dot plot: GO terms after simplify
dotplot(ego2, showCategory = 15)
## GSEA plot: pvalue_table needs gridExtra; use a width of at least 10 inches so the table is not clipped
p <- gseaplot2(gs, geneSetID = 1:2, pvalue_table = TRUE)
ggsave("gsea.png", p, width = 10, height = 6, dpi = 300)Volcano plot
Put -log10(padj) on the y axis with threshold lines matching the text (padj 0.05, |log2FC| 1 here). Remove genes with NA padj first, or ggplot warns about missing values and drops points. padj of extremely significant genes can underflow to 0, giving Inf after -log10; replace 0 with the smallest non-zero value and note it in the legend. Label genes chosen by padj or a predefined list.
Dot plot
By default dotplot puts GeneRatio on the x axis, gene count as point size and p.adjust as color. The GeneRatio denominator is the number of annotated DE genes, not all DE genes. Show terms after simplify so the top 15 are not all levels of the same process.
GSEA plot
The peak of the upper curve is the enrichment score. A positive peak on the left means the set is up in the treatment group (the ranking direction follows your contrast). The vertical bars mark set members in the ranking, and members before the peak form the leading edge. gseaplot2's pvalue_table needs gridExtra, otherwise it errors with “The package "gridExtra" is required for tableGrob2()”.
- The volcano plot uses pvalue on the y axis while the filter uses padj, so the “significant” points do not match the gene list in the text.
- The enrichment figure shows ORA of up-regulated genes only, but the text says pathways are “activated”; state direction from GSEA NES or from separate up/down ORA.
- The dot plot does not state the background, database release or correction method.
- GSEA curves with p.adjust > 0.05 appear in the main figure and are described as “significantly enriched”.
- The DE gene count differs from the number of rows in the supplementary table, usually because SYMBOL de-duplication or ID conversion losses were not reflected in the number.
Chinese community
These come from author-tested posts on CSDN and Jianshu and were checked against clusterProfiler 4.18.4 and createKEGGdb 0.0.5.
Local KEGG data: verified independently by two authors
CSDN authors hutong_oncology and qq_41547057 both hit online enrichKEGG failures or a sudden “No gene can be mapped” in 2023 and fixed it by building a local KEGG.db with createKEGGdb::create_kegg_db. The first also notes that a local database prevents online KEGG updates from changing results between runs. The workflow still works with current versions on this page and matches the online results.
Old KEGG.db builds missing Description break plots
hutong_oncology reports that some generated KEGG.db packages lacked pathway descriptions: the enrichment table was produced, but barplot and dotplot failed with “Error in ans[ypos] <- rep(yes, length.out = len)[ypos]” and “'x' is NULL so the result will be NULL”. Many posts blame a strict cutoff; the actual cause is the missing Description. The package built by createKEGGdb 0.0.5 here has complete descriptions; if you see this error, check head(ekk$Description) first.
The R.utils::setOption download trick no longer works
Older posts fix KEGG download failures with R.utils::setOption("clusterProfiler.download.method", "auto"). In clusterProfiler 4.18.4, KEGG downloads go through yulab.utils::yread, which calls readLines, and the source no longer reads this option. Fix network problems with a system proxy or a local KEGG.db.
“padj < 0.05 guarantees at most 5% false positives” is inaccurate
Jianshu and similar tutorials often explain padj this way. BH correction controls the expected proportion of false discoveries; in a single analysis the actual proportion can exceed 5%. Write “FDR controlled at 5%” in the methods, not “guaranteed”.
Hand it to an agent
You can describe the full differential expression and enrichment task in one instruction.
Example instruction: “Run DESeq2 on counts.txt from featureCounts and samples.csv with design ~ batch + condition, comparing trt to ctrl; cutoffs padj < 0.05 and |log2FC| > 1; draw a volcano plot with apeglm-shrunken LFC; run GO BP and KEGG enrichment with detected genes as the background, using local KEGG data and recording its release; run GSEA on hallmark and GO BP ranked by the Wald stat with a fixed seed; output the DE table, enrichment tables and three plots.”
The agent installs R and the Bioconductor packages in an isolated cloud computer and follows this page: it checks that the sample table matches the count columns, pre-filters, runs DESeq2 and lfcShrink, converts IDs and reports the loss rate, runs ORA with a universe, GSEA and the plots. The workspace keeps the scripts, sessionInfo, KEGG data date, result tables and figures, so the analysis can be reproduced. The agent reviews its results adversarially, for example checking that volcano thresholds match the filter, that a universe was set and that GSEA used a fixed seed.
You still need to check that group and batch labels match your lab records, that the design reflects the real pairing structure, that the cutoffs follow conventions in your field, and that the enriched pathways make biological sense.
References
- DESeq2 vignette (Bioconductor release) — Pre-filtering code, the three reasons padj is NA, the lfcShrink comparison table, the no-replicates FAQ, Cook's distance and the 7-replicate replacement rule
- edgeR User's Guide 2.13 What to do if you have no replicates — Four options without replicates and BCV reference values 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 — FDR inflation of DESeq2/edgeR in large population samples and the Wilcoxon recommendation
- clusterProfiler book: KEGG analysis — Online KEGG data and KEGG.db not updated since 2012; its gseKEGG example still uses the discouraged nPerm
- YuLab-SMU/createKEGGdb (GitHub) — Builds a local KEGG.db; version 0.0.5 tested on this page
- Bioconductor support: BgRatio N value does not equal set gene universe size — The universe is intersected with annotated genes
- GSEA User Guide (GSEA-MSigDB documentation) — Meaning of the 25% FDR q-value cutoff, phenotype vs gene_set permutation, definition of the leading edge
- GSEAPreranked module documentation (v7) — Preranked uses gene_set permutation only; defaults min 15/max 500, weighted scoring, timestamp seed; ranking values must not tie
- fgsea (Bioconductor) — Multilevel GSEA implementation; the default backend of clusterProfiler GSEA
- msigdbr (CRAN) — collection argument and ncbi_gene column since 10.0
- Biostars: GSEA on preranked list with weighted enrichment statistic — Forum experience: fgsea vs GSEA software discrepancies; the “fgsea cannot weight” claim is outdated
- CSDN hutong_oncology: errors when creating KEGG.db and fixes — Forum experience: createKEGGdb workflow and the plotting error text caused by missing Description
- CSDN qq_41547057: clusterProfiler KEGG enrichment errors and fixes — Forum experience: switching to a local KEGG.db after online enrichKEGG returned No gene can be mapped
- Jianshu: why filter DE genes by log2FC and padj — Forum experience: Chinese explanation of common cutoffs; its FDR wording needs correction
- KEGG REST API: info/kegg — Query the KEGG data release date