GO, KEGG, and GSEA Enrichment with clusterProfiler in R

A DESeq2 results table is a list of gene names with numbers attached to them. It is not, by itself, biology. Enrichment analysis is the step that turns "these 400 genes went up" into "the interferon response and cell cycle pathways went up." clusterProfiler is the Bioconductor package most people reach for to make that jump, because it covers both flavors of enrichment testing (over-representation and gene set enrichment) against both Gene Ontology and KEGG in one consistent interface.

This tutorial takes a typical DESeq2 output data frame and walks through the full downstream path: mapping gene identifiers to Entrez IDs, running over-representation analysis with enrichGO and enrichKEGG, running ranked-list GSEA with gseGO and gseKEGG, and visualizing the results with dotplot and cnetplot.

What over-representation and GSEA actually test

These two methods answer different questions, and mixing them up is the single most common enrichment mistake.

ORA is simple and interpretable but throws away information: a gene with a log2 fold change of 1.9 and one with 0.6 count identically as long as both cleared your p-value cutoff, and everything below the cutoff is ignored entirely. GSEA keeps the full ranking, so it can detect a pathway where every member gene shifted moderately in the same direction, a pattern ORA is often too blunt to see. Good practice is to run both and treat GSEA as the more sensitive, less threshold-dependent complement to ORA rather than a replacement for it.

Setup and packages

Everything here assumes you already have a DESeqResults object or an equivalent data frame with a gene identifier column, log2FoldChange, and padj, produced by the workflow in our RNA-seq differential expression guide. Install the required Bioconductor packages once, inside an R session:

if (!requireNamespace("BiocManager", quietly = TRUE)) {
  install.packages("BiocManager")
}
BiocManager::install(c("clusterProfiler", "org.Hs.eg.db", "enrichplot", "DOSE"))

org.Hs.eg.db is the human genome-wide annotation database, it is what lets clusterProfiler translate between gene symbols, Ensembl IDs, and Entrez IDs, and it supplies the GO-to-gene mappings that enrichGO/gseGO search against. Swap it for org.Mm.eg.db (mouse) or the appropriate org.*.eg.db package if you are not working with human data.

Step 1: From DESeq2 results to a clean gene list

Start from a standard results() data frame, drop rows with NA statistics (DESeq2 assigns NA to genes filtered out by independent filtering or outlier detection), and keep the columns you need:

library(DESeq2)
library(clusterProfiler)
library(org.Hs.eg.db)
library(enrichplot)

res_df <- as.data.frame(res)              # `res` is your DESeqResults object
res_df$ensembl <- rownames(res_df)
res_df <- res_df[!is.na(res_df$padj) & !is.na(res_df$log2FoldChange), ]

If your rownames are Ensembl gene IDs with version suffixes (e.g. ENSG00000141510.19), strip the suffix before mapping, since org.Hs.eg.db keys on unversioned IDs:

res_df$ensembl <- sub("\\..*$", "", res_df$ensembl)

Step 2: Convert identifiers to Entrez IDs

enrichKEGG and gseKEGG require Entrez Gene IDs, KEGG's own gene identifiers for most organisms are just Entrez IDs, so clusterProfiler expects them directly. Use bitr (biological ID translator) to map your Ensembl IDs across:

id_map <- bitr(
  res_df$ensembl,
  fromType = "ENSEMBL",
  toType   = c("ENTREZID", "SYMBOL"),
  OrgDb    = org.Hs.eg.db
)

res_df <- merge(res_df, id_map, by.x = "ensembl", by.y = "ENSEMBL")

Duplicate Entrez IDs can appear when several Ensembl IDs map to the same gene; keep the row with the smallest padj for each ENTREZID before moving on:

res_df <- res_df[order(res_df$padj), ]
res_df <- res_df[!duplicated(res_df$ENTREZID), ]

Step 3: Over-representation analysis with enrichGO and enrichKEGG

Define your significant gene list, typically padj < 0.05 and some fold-change threshold, and use every tested, successfully-mapped gene as the background (universe). Using the whole genome as background instead of your actual tested universe is a classic source of inflated significance, so always pass universe explicitly.

sig_genes <- res_df$ENTREZID[res_df$padj < 0.05 & abs(res_df$log2FoldChange) > 1]
universe  <- res_df$ENTREZID

ego <- enrichGO(
  gene          = sig_genes,
  universe      = universe,
  OrgDb         = org.Hs.eg.db,
  keyType       = "ENTREZID",
  ont           = "BP",        # Biological Process; also "MF", "CC", or "ALL"
  pAdjustMethod = "BH",
  pvalueCutoff  = 0.05,
  qvalueCutoff  = 0.2,
  readable      = TRUE
)

ekegg <- enrichKEGG(
  gene         = sig_genes,
  universe     = universe,
  organism     = "hsa",        # KEGG organism code for human
  pvalueCutoff = 0.05
)

readable = TRUE on enrichGO swaps Entrez IDs back to gene symbols in the output for readability; enrichKEGG does not take that argument, so convert its geneID column afterward with setReadable() if you want symbols there too:

ekegg <- setReadable(ekegg, OrgDb = org.Hs.eg.db, keyType = "ENTREZID")

Inspect the top terms as a data frame:

head(as.data.frame(ego)[, c("Description", "GeneRatio", "BgRatio", "p.adjust")])

You should see something like:

                          Description GeneRatio  BgRatio  p.adjust
1 type I interferon signaling pathway    18/342 62/18670   1.2e-09
2                   response to virus    24/342 210/18670  3.8e-08
3          regulation of cell cycle     31/342 512/18670  6.1e-06

Step 4: Ranked-list GSEA with gseGO and gseKEGG

For GSEA you skip the significance cutoff entirely and instead rank every mapped gene by log2FoldChange, from most upregulated to most downregulated. The ranked vector must be named by Entrez ID and sorted in decreasing order, clusterProfiler will error or silently mis-rank otherwise.

gene_list <- res_df$log2FoldChange
names(gene_list) <- res_df$ENTREZID
gene_list <- sort(gene_list, decreasing = TRUE)

gsea_go <- gseGO(
  geneList     = gene_list,
  OrgDb        = org.Hs.eg.db,
  keyType      = "ENTREZID",
  ont          = "BP",
  minGSSize    = 15,
  maxGSSize    = 500,
  pvalueCutoff = 0.05,
  eps          = 0
)

gsea_kegg <- gseKEGG(
  geneList     = gene_list,
  organism     = "hsa",
  minGSSize    = 15,
  pvalueCutoff = 0.05,
  eps          = 0
)

eps = 0 disables the default 1e-10 boundary on the enrichment score p-value estimate, which gives more precise p-values for strongly enriched terms at some extra compute cost, worth it for anything you plan to report. minGSSize/maxGSSize filter out gene sets that are too small (noisy) or too large (uninformative) to test meaningfully; 15-500 is a reasonable default range for GO/KEGG terms.

Step 5: Visualize with dotplot and cnetplot

A dotplot is the fastest way to scan the top enriched terms across both significance and gene count:

library(enrichplot)
dotplot(ego, showCategory = 15, title = "GO Biological Process (ORA)")
dotplot(gsea_kegg, showCategory = 15, title = "KEGG (GSEA)")

Dot size typically encodes GeneRatio (ORA) or setSize (GSEA), and color encodes adjusted p-value, darker/redder dots are more significant, depending on the palette. For GSEA results specifically, a ridgeplot(gsea_kegg) is often more informative than a dotplot, since it shows the actual distribution of fold changes within each enriched gene set rather than just a summary statistic.

To see which genes drive a handful of top terms and how those terms share genes, use cnetplot (category-gene network plot):

cnetplot(
  ego,
  showCategory = 5,
  categorySize = "pvalue",
  foldChange   = gene_list
)

Passing foldChange = gene_list colors the gene nodes by their log2 fold change, so up- and down-regulated genes within an enriched term are visually distinguishable, this is usually the single most useful plot for a figure, since it connects the abstract pathway name back to the actual genes and their direction of change.

For a single GSEA term, gseaplot2(gsea_kegg, geneSetID = 1) reproduces the classic running-enrichment-score plot from the original GSEA software, with the ranked gene list shown as a barcode along the bottom.

Common gotchas

Wrap-up

The pattern is always the same: get clean Entrez IDs, decide whether you want a threshold-based hit list (ORA) or a full ranking (GSEA), run the matching clusterProfiler function against GO and/or KEGG, and visualize with dotplot for a broad scan or cnetplot for the gene-level story behind your top terms. Running both ORA and GSEA side by side, rather than picking one, is the best insurance against a real signal in the ranking getting lost because it did not quite clear an arbitrary p-value cutoff.

From here, the natural next steps are pairing this with a proper RNA-seq differential expression workflow if you have not run DESeq2 yet, or extending the same GSEA logic to single-cell cluster markers using the approach from our scRNA-seq for beginners guide.