Gene set enrichment analysis (GSEA) helps translate gene-level expression changes into broader biological processes. It starts with a ranked list of all tested genes, usually based on a signed statistic from differential gene expression analysis, and asks whether genes from the same pathway tend to cluster near the top or bottom of that ranking.
The quality of the result depends heavily on the experimental design and the gene ranking. In both bulk and single-cell studies, the biological replicate may be a donor, animal, or independently treated sample. Bulk RNA-seq measures the average signal across all cells in that sample, whereas single-cell RNA-seq can resolve changes within specific cell types or populations. For condition comparisons, cells are often aggregated into donor- and cell-type-specific pseudobulk profiles before differential gene expression analysis. Tools such as clusterProfiler can then run GSEA against GO, KEGG, or custom gene sets and return the normalized enrichment score, adjusted p-value, and leading-edge genes, together with useful visualizations.
What Is Gene Set Enrichment Analysis (GSEA)?
Gene set enrichment analysis (GSEA) is a method for identifying biological pathways in gene expression data. It ranks genes according to their association with a comparison and then asks whether genes from the same pathway cluster near the top or bottom of that ranking.
Unlike methods that use only genes passing a significance threshold, such as gene ontology (GO) enrichment analysis, GSEA uses the full ranked gene list. This helps detect small but coordinated changes across many genes belonging to the same biological pathway.
What goes into GSEA?
Ranked list of genes
GSEA needs two main inputs. The first is a ranked list of genes.
Before discussing gene ranking, it is important to define the contrast of interest. For example, we may want to determine whether immune-response pathways are more active in treated samples than controls, or whether DNA-repair pathways are enriched in drug-resistant cells.
Because GSEA evaluates pathway enrichment relative to a specific comparison, genes must first be ranked according to a statistic that captures both the direction and strength of that contrast.
For bulk RNA-seq, genes are typically ranked using the DESeq2 Wald statistic. In single-cell RNA-seq, counts are often aggregated into pseudobulk profiles by donor and cell type, after which differential expression can be performed with DESeq2 and genes are ranked using the same Wald statistic for GSEA.
Gene sets
The second input is the gene sets. Each gene set represents a pathway, biological process, or regulatory program.
Common sources include:
- Gene Ontology (GO)
- KEGG
- Reactome
- MSigDB
You can also provide custom gene sets.
Gene sets are typically restricted to about 15–500 genes, although some tools, including clusterProfiler, use a default minimum of 10.
Gene-set size affects the raw enrichment score (ES), which measures how strongly the genes in a set are concentrated near the top or bottom of the ranked gene list. Each matching gene has a larger influence in a small gene set than in a large one.
For example, a small gene set may produce an ES of 0.60, while a larger set with a comparable enrichment pattern produces an ES of 0.30. This occurs because hits in a smaller list make up a larger proportion of that set and therefore contribute more strongly to the score. However, random small sets may also achieve higher scores by chance.
To account for this, GSEA calculates the normalized enrichment score (NES) by dividing the observed ES by the typical null ES obtained from repeated permutations for that gene set.
For example:
- Small gene set:
- ES = 0.60
- Mean null ES = 0.30
- NES = 2.0
- Large gene set:
- ES = 0.30
- Mean null ES = 0.15
- NES = 2.0
Both gene sets are therefore enriched about twice as strongly as expected under their respective null distributions.
NES makes scores from differently sized gene sets more comparable, but very small or very large gene sets are still excluded. Very small gene sets are sensitive to individual genes, whereas very large sets may dilute specific biological signals. Size filtering therefore improves both statistical robustness and biological interpretability.
What Comes Out of GSEA?
The main output is typically a table with one row per gene set.
Important elements include:
- Enrichment score (ES)
- Normalized enrichment score (NES)
- Nominal p-value
- Adjusted p-value
- Gene-set size
- Leading-edge genes
The normalized enrichment score (NES) is the primary statistic used to compare enrichment across gene sets.
- A positive NES indicates that the gene set is concentrated near the top of the ranked list.
- A negative NES indicates that the gene set is concentrated near the bottom.
Example output:
| ID | Description | setSize | NES | p.adjust |
|---|---|---|---|---|
| GO:0006955 | Immune response | 142 | 2.15 | 0.003 |
| GO:0007049 | Cell cycle | 96 | -1.82 | 0.018 |
Results can be visualized using:
- Dot plots
- Enrichment curves
- Ridge plots
- Gene–pathway networks
A dot plot provides a compact overview of several pathways, while an enrichment curve shows where genes from one pathway occur in the ranked list.
Figure 1: GSEA of a simulated RNA-seq experiment (Treated vs Control, n = 12 per group).
(A) Running enrichment score across 13,494 genes ranked by the DESeq2 Wald statistic, from most upregulated in Treated (left) to most upregulated in Control (right). Tick marks indicate genes belonging to each gene set, while the bottom strip shows the Wald statistic, ranging from −20.3 to +26.6 and crossing zero at rank 6,526. Curves peak early for pathways enriched in Treated (E2F_TARGETS, G2M_CHECKPOINT) and trough late for pathways enriched in Control (OXIDATIVE_PHOSPHORYLATION, FATTY_ACID_METABOLISM).
(B) Variance-stabilized counts (row z-scores) for the leading-edge genes of E2F_TARGETS (125 of 185 genes), representing the genes contributing most strongly to pathway enrichment.
(C) Normalized enrichment scores (NES) for all 18 gene sets. Color indicates −log₁₀ FDR q-value and dot size represents the number of leading-edge genes. Pathways range from OXIDATIVE_PHOSPHORYLATION (NES −3.66) to E2F_TARGETS (NES +3.31), while APICAL_JUNCTION and KRAS_SIGNALING_DN are not significantly enriched. Counts were simulated using a negative-binomial model, analyzed with DESeq2, and enrichment was calculated using the weighted Kolmogorov–Smirnov statistic with 1,000 gene-set permutations.
Scaling Up Single-Cell: How to Analyze Large-Scale Datasets
Watch how our bioinformatics experts tackle large-scale RNA-seq data analysis from start to finish.
Experimental Design: Bulk or Single-Cell RNA-seq?
GSEA is only as reliable as the differential gene expression analysis used to generate the ranked list. This requires a clearly defined experimental design, appropriate contrasts, and sufficient biological replication, meaning multiple independent experimental units, or independent observations, per condition.
In bulk RNA-seq, the biological replicate is typically an individual donor, animal, or tissue sample, with each sample contributing a single expression profile. However, because bulk RNA-seq averages expression across all cells, it may obscure cell-type-specific pathway changes or confound changes in gene expression with differences in cell-type abundance.
Single-cell RNA-seq preserves cellular resolution, allowing differential expression analysis and GSEA to be performed separately within defined cell types. This makes it possible to identify which cell populations drive pathway-level changes and distinguish cell-type-specific responses that may be diluted or undetectable in bulk data.
Although thousands of cells may be measured, the experimental units remain the biological donors or samples from which those cells were collected. To respect this structure, counts are often aggregated across cells from the same donor and cell type to create pseudobulk profiles. DESeq2 is then applied to these pseudobulk samples, treating each donor as an independent biological replicate while individual cells contribute to the expression estimate rather than serving as independent replicates.
Experimental conditions should also be balanced across batches and sequencing runs. Relevant variables such as age, sex, batch, time point, or paired measurements should be recorded and included in the statistical model where appropriate. As with differential gene expression analysis, a confounded experimental design can produce a misleading gene ranking and, consequently, misleading pathway enrichment results.
How Does GSEA Work?
Ranking the Genes
Suppose we compare treated samples with controls. A differential gene expression model produces one statistic for every tested gene.
- A large positive statistic means the gene is more strongly associated with the treated condition.
- A large negative statistic means the gene is more strongly associated with the control condition.
- Values close to zero indicate little evidence of differential expression.
A simplified ranking might look like this:
| Gene | Test statistic | Position |
|---|---|---|
| ISG15 | 8.2 | Top |
| IFIT1 | 7.6 | Top |
| STAT1 | 6.9 | Top |
| ACTB | 0.1 | Middle |
| MCM5 | -5.8 | Bottom |
| MKI67 | -7.1 | Bottom |
The ranking must preserve direction. Genes should therefore not be ranked using:
abs(log2FoldChange)- Unsigned p-values
These measures mix genes changing in opposite directions.
A signed test statistic is generally the preferred choice:
- DESeq2: use the
statcolumn. - limma: use the moderated t-statistic.
- Other methods: use a signed statistic that combines effect direction with statistical evidence.
Ranking by log2 fold change is possible, but genes with low counts or greater uncertainty may receive disproportionately extreme ranks.
Calculating the Enrichment Score
For each gene set, GSEA walks from the top of the ranked list to the bottom. The running score increases when a gene belongs to the set and decreases when it does not. Genes with more extreme ranking statistics usually have a stronger effect on the score.
The enrichment score (ES) is the largest deviation of this running score from zero.
Consider an interferon-response gene set. If many of its genes, including ISG15, IFIT1, IRF7 and STAT1, occur near the top of the ranked list, the running score rises quickly, and the pathway receives a positive enrichment score.
The signal does not rely on one highly significant gene. Instead, it reflects a coordinated shift across several related genes.
Because the raw enrichment score depends partly on gene-set size, GSEA calculates a normalized enrichment score (NES). Statistical significance is evaluated against a null distribution and corrected for the number of gene sets tested.
The genes that contribute to the score up to its maximum deviation form the leading edge. These genes often provide the clearest link between a broad pathway name and the underlying expression data.
clusterProfiler: Your Go-To R Package for GSEA
The clusterProfiler R package provides a convenient workflow for running, interpreting, and visualizing GSEA. One of its main strengths is its unified interface.
The GSEA() function accepts custom gene sets through a two-column term-to-gene table, where each row links a biological term or pathway to one of its genes. Functions such as gseGO() and gseKEGG() connect directly to widely used pathway databases.
Assume that res is a DESeq2 results object. Its stat column contains the Wald statistic used to rank genes.
library(clusterProfiler)
library(enrichplot)
library(org.Hs.eg.db)
gene_list <- res$stat
names(gene_list) <- rownames(res)
gene_list <- gene_list[!is.na(gene_list)]
gene_list <- sort(gene_list, decreasing = TRUE)The gene_list must be a named numeric vector containing the full ranked gene list, not only statistically significant genes.
go_result <- gseGO(
geneList = gene_list,
OrgDb = org.Hs.eg.db,
keyType = "ENSEMBL",
ont = "BP",
minGSSize = 10,
maxGSSize = 500,
pvalueCutoff = 0.05,
pAdjustMethod = "BH",
verbose = FALSE
)The gseGO() function tests Gene Ontology Biological Process terms. The output contains one row for each enriched biological process.
An illustrative result might look like this:
| ID | Description | setSize | NES | p.adjust |
|---|---|---|---|---|
| GO:0006955 | Immune response | 142 | 2.15 | 0.003 |
| GO:0007049 | Cell cycle | 96 | -1.82 | 0.018 |
A positive NES indicates enrichment near the top of the ranked gene list, whereas a negative NES indicates enrichment near the bottom. The setSize column shows the number of genes from the pathway included in the analysis, and p.adjust reports the multiple-testing-adjusted p-value.
Results can be inspected and visualized directly:
as.data.frame(go_result)
dotplot(go_result, showCategory = 15)
gseaplot2(go_result, geneSetID = 1)The dot plot summarises several enriched pathways, while gseaplot2() displays the running enrichment score together with the positions of pathway genes within the ranked list. Before running the analysis, ensure that gene identifiers match the selected keyType and that all identifiers are unique.
GSEA Data Analysis and Bioinformatics Support
GSEA is only as powerful as the data and analysis behind it. From experimental design and differential gene expression analysis to pathway enrichment and biological interpretation, every step influences the quality of the final results. Single Cell Discoveries provides comprehensive Data Consulting services to help researchers extract the maximum value from their transcriptomics data. Whether you need support with bulk RNA-seq, single-cell RNA-seq, or downstream analyses such as GSEA, our bioinformatics experts work alongside your team to deliver robust, reproducible, and publication-ready results.
Scaling Up Single-Cell
Learn about practical solutions for working with large datasets
You will see how expert-driven data support can turn computational hurdles into scientific discoveries, and why partnering with the right team makes all the difference.
