Running GSEA on Your Own RNA-seq Data

By the end of this guide you will have a table of gene sets ranked by normalized enrichment score with false discovery rates, a leading-edge gene list for each set that survives filtering, and a set of single-sample pathway scores you can track across timepoints. You need a gene-level count matrix from your RNA-seq run, and the quant.sf files from salmon or ReadsPerGene.out.tab from STAR are both fine. You also need a sample table describing what you are comparing, and R 4.3 or later. You need these packages as well.

• DESeq2
• tximport
• fgsea
• msigdbr
• GSVA
• clusterProfiler

Budget about an hour of compute on a laptop. Running fgsea on 20,000 genes and 1,500 gene sets takes under a minute, and the slow part is everything upstream of it.

1. Decide what you are contrasting, before you touch the data

GSEA answers a question about an ordered list of genes, so the first job is to produce an ordering that means something. The method itself is agnostic about where the ordering comes from, which is exactly why it is easy to misuse: a ranked list built from a confounded comparison produces confidently significant pathways that describe the confound.

For a single person, there are three contrasts that hold up. The first is longitudinal within-person, with several draws in one state versus several draws in another. Examples are before and after a sustained dietary change, a training block or a documented infection. You treat replicate draws within a state as the replicates. The second is within-person paired across a manipulation, modeling draw as a blocking factor. The third is you versus an external cohort, which is the weakest of the three because library prep, globin depletion, sequencing platform, time of day, and cell composition all differ between your sample and a public dataset, and every one of those differences will land on top of your biology.

If your RNA comes from whole blood, one more thing matters more than any pathway database. Whole-blood expression is dominated by cell composition, so a shift in neutrophil-to-lymphocyte ratio will light up sets named for inflammation, interferon signaling, and neutrophil degranulation without any per-cell change in regulation. Run a deconvolution first, and carry the estimated fractions as covariates in the model or at least as an interpretation caveat. This is the single most common way a clean-looking personal GSEA result turns out to be a blood count in disguise.

2. Build the count matrix and fit the model

Import transcript quantifications and collapse to gene level, then fit a model with the covariates you decided on above. The Wald statistic that DESeq2 reports is the ranking metric we prefer, because it already combines effect size and precision, and unlike the raw log2 fold change it does not put a gene with three counts and a tenfold swing at the top of the list.

library(tximport); library(DESeq2)

tx2gene <- read.csv("tx2gene.gencode.v44.csv")            # tx_id, gene_symbol
files   <- file.path("quants", samples$id, "quant.sf")
names(files) <- samples$id

txi <- tximport(files, type = "salmon", tx2gene = tx2gene,
                countsFromAbundance = "lengthScaledTPM")

dds <- DESeqDataSetFromTximport(txi, colData = samples,
                                design = ~ neut_frac + state)

# keep genes with a real chance of being informative
keep <- rowSums(counts(dds) >= 10) >= min(table(samples$state))
dds  <- dds[keep, ]

dds <- DESeq(dds)
res <- results(dds, name = "state_post_vs_pre")

Filtering low-count genes before ranking is not cosmetic. Genes with near-zero counts produce statistics clustered tightly around zero and, worse, exact ties, and ties degrade the enrichment score computation because the order within a tied block is arbitrary. Dropping genes that are not detected in at least the size of your smallest group typically removes 30 to 45 percent of annotated genes from a blood sample and sharpens everything downstream.

One deliberate omission: do not apply lfcShrink and then rank by the shrunken log2 fold change. Shrinkage is the right tool for reporting effect sizes and for volcano plots, but it compresses exactly the low-information genes that GSEA already handles correctly through the statistic, and it introduces its own ties at zero.

3. Turn results into a clean ranked vector

GSEA consumes a named numeric vector, one value per gene, sorted descending. Two problems appear here reliably: duplicate gene symbols (from multi-mapping identifiers, readthrough transcripts, or an Ensembl-to-symbol join) and NA statistics from genes that DESeq2 filtered independently.

library(dplyr); library(tibble)

ranks <- as.data.frame(res) |>
  rownames_to_column("symbol") |>
  filter(!is.na(stat), symbol != "") |>
  group_by(symbol) |>
  slice_max(abs(stat), n = 1, with_ties = FALSE) |>   # one row per symbol
  ungroup() |>
  arrange(desc(stat))

rv <- setNames(ranks$stat, ranks$symbol)
length(rv); sum(duplicated(names(rv)))   # expect 14k-18k and 0
summary(rv)

Collapsing duplicates by maximum absolute statistic is a judgment call. Taking the mean is the other defensible option and will slightly flatten the list. What you must not do is leave duplicates in, because fgsea will treat the same gene as two independent hits and inflate the enrichment score for any set containing it.

Look at summary(rv). If the median is far from zero, your normalization or your design is off, and no pathway result from that list is interpretable.

4. Choose gene sets, and choose few

Gene set collections differ enormously in redundancy and in how much interpretation they support. We start with MSigDB Hallmark (50 sets, hand-curated, deliberately non-overlapping) because a Hallmark result is short enough to read in full and coherent enough to reason about. Reactome and the GO Biological Process branch come second, and they are where redundancy becomes a real problem: hundreds of nested parent and child terms will report the same underlying signal many times over.

library(msigdbr)

hall <- msigdbr(species = "Homo sapiens", collection = "H")
gs_h <- split(hall$gene_symbol, hall$gs_name)

reac <- msigdbr(species = "Homo sapiens", collection = "C2", subcollection = "CP:REACTOME")
gs_r <- split(reac$gene_symbol, reac$gs_name)

This is also the place to settle the question of whether GSEA and “pathway analysis” are the same thing, and how either relates to GO. Pathway analysis is the general activity of summarizing a gene-level result at the level of annotated groups, and GSEA is one method for doing it. GO is a database of annotations, not a method, and you can run GO terms through GSEA or through overrepresentation analysis, which are different statistical procedures on the same annotations. The full workflow, from ranked list through enrichment testing to network visualization of the redundant results, is laid out well by Reimand and colleagues ¹.

On the choice between overrepresentation analysis (ORA) and GSEA: ORA takes a list of genes that passed a threshold and asks whether a set is enriched among them relative to a background, while GSEA uses every gene in your ranked list and asks whether a set’s members cluster toward either extreme. For a personal dataset with few replicates, where the differential expression list at FDR 0.05 might contain twelve genes, GSEA is clearly better because it does not need a threshold and it detects coordinated small shifts across a set. ORA remains useful when you have a genuinely discrete gene list, for example the genes carrying rare variants in your genome. The background must then be the set of genes you could have detected. It must not be all 20,000 in the annotation. Using the wrong background is the most common error in ORA and it makes p-values look far better than they are ^{1 2}.

5. Run fgsea

Use fgsea, which implements the multilevel adaptive approach and gives accurate p-values far into the tail without running a million permutations.

library(fgsea)
set.seed(42)

fg <- fgsea(pathways = gs_h,
            stats    = rv,
            minSize  = 15,
            maxSize  = 500,
            eps      = 0,
            nproc    = 4)

fg <- fg[order(padj), ]
fg[padj < 0.05, .(pathway, NES, size, padj, log2err)]

The size filters matter. Sets below about 15 genes give unstable scores that are dominated by one or two genes, and sets above 500 are usually so broad (“metabolic process”, “signal transduction”) that enrichment tells you nothing actionable. Setting eps = 0 lets the multilevel estimator report p-values below the default floor, which you will want for the strongest sets. Check log2err: if it is large or NA for a pathway you care about, the p-value estimate has not converged and you should rerun with higher nPermSimple.

For Reactome or GO, add a redundancy collapse so you are reading distinct signals rather than a term hierarchy:

fg_r <- fgsea(gs_r, rv, minSize = 15, maxSize = 500, eps = 0)
main <- collapsePathways(fg_r[order(pval)][padj < 0.05, ], gs_r, rv)
fg_r[pathway %in% main$mainPathways, .(pathway, NES, size, padj)]

clusterProfiler is the alternative front end, and it is the better choice if you want the ID translation, the GO semantic-similarity simplification, and the plotting to come from one package, or if you plan to extend the same analysis to other data types later ^{3 4}.

6. Read the enrichment score, and then the plot

The enrichment score (ES) is the maximum deviation from zero of a running sum that walks down your ranked list, incrementing when it hits a set member and decrementing when it does not. A large positive ES means the set’s genes pile up at the top of the list, and a large negative ES means they pile up at the bottom. The raw ES depends on set size, so always interpret the normalized enrichment score (NES), which divides by the mean of the permutation null for that set and is therefore comparable across sets of different sizes. In practice, NES above roughly 1.5 in absolute value with an FDR under 0.05 is where we start paying attention. The sign tells you direction relative to how you defined the contrast, so confirm which level was the numerator.

Plot the running score for your top set and read three features.

plotEnrichment(gs_h[["HALLMARK_INTERFERON_ALPHA_RESPONSE"]], rv) +
  ggplot2::labs(title = "IFN-alpha response")

The curve should rise steeply and peak early if the set is genuinely up. The black tick marks below show where set members sit in the ranked list, and you want them concentrated near the left edge rather than scattered evenly, because an even scatter with a modest peak often reflects a handful of extreme genes plus noise. The leading edge is the set of genes from the start of the list up to the peak, and it is the part worth reading gene by gene:

fg[pathway == "HALLMARK_INTERFERON_ALPHA_RESPONSE", leadingEdge][[1]]

If the leading edge for an interferon set is six genes and four of them are also in your top-ten differential list, the pathway result is essentially restating those four genes. If it is forty genes each moving modestly in the same direction, that is the coordinated shift GSEA was designed to find. This distinction does not appear anywhere in the FDR column, and it is the main thing that separates a result you can build on from one you cannot.

7. Score single samples with GSVA when there is no contrast

Sooner or later you will have one RNA-seq sample and no comparison group, and GSEA has nothing to rank. GSVA solves a different problem: it converts an expression matrix into a gene-set-by-sample matrix, where each cell is a score for how coordinately high or low that set’s genes are in that sample relative to the other samples in the matrix.

library(GSVA)
vst_mat <- assay(vst(dds, blind = FALSE))

par <- gsvaParam(exprData = vst_mat, geneSets = gs_h,
                 kcdf = "Gaussian", minSize = 15, maxSize = 500)
sc <- gsva(par)      # 50 x n_samples

Use kcdf = "Gaussian" for log-scale continuous input such as VST or log-CPM, and "Poisson" only if you are feeding integer counts. The interpretation caveat is the whole story here: a GSVA score is relative to the sample set you passed in, so adding or removing samples changes every score. Scoring one sample alone is meaningless, and scoring your sample inside a public cohort mostly measures batch. What does work for a personal profile is scoring your own timepoints together, which makes the score a within-person trajectory. Then you can regress a set score on time, exercise load, or sleep, and the reference is your own distribution.

So: GSEA tests one ranked list for enrichment and returns a p-value per set, while GSVA returns a score per set per sample with no inferential claim attached. They answer different questions and are not substitutes.

8. Run the negative controls

Before you act on anything, break the analysis on purpose. Permute your condition labels, refit, rerun fgsea, and confirm that essentially nothing survives FDR 0.05. Then run fgsea with size-matched random gene sets sampled from your detected genes and confirm the NES distribution is centered at zero.

rand <- lapply(1:500, function(i) sample(names(rv), 100))
names(rand) <- paste0("RAND", 1:500)
summary(fgsea(rand, rv, minSize = 15, maxSize = 500)$NES)

If permuted labels still yield significant Hallmark sets, you have structure in the data that your design does not model, usually batch or cell composition. The other check worth doing is external: if the biology you think you are seeing has been characterized in published cohorts, ask whether the direction and the specific sets agree with those studies. Cross-study GSEA comparisons of this kind have been used to establish which transcriptional signals reproduce across independent datasets and which are study-specific ⁵.

Finally, interpretation has a hard boundary. A significant NES for a coagulation or interferon set is a statement about coordinated transcript abundance in one tissue at one moment. It is not a finding about your health, it does not name a cause, and it cannot be mapped onto a condition or a treatment. If a result seems clinically meaningful to you, the next step is a conversation with a physician who can order validated diagnostic testing. The result is not a reason to change anything you are doing.

Common problems

Ties in the ranked list. fgsea warns when it sees them, and the usual cause is ranking by p-value with many exact zeros or by a statistic computed on unfiltered low-count genes. Filter first, rank by the Wald statistic, and the warning usually disappears.

Everything is significant. Hundreds of GO terms at FDR 0.001 almost always means the ranked list is driven by one dominant axis, such as cell composition in blood or a large batch effect. Add the covariate, or accept that you are reading that axis.

Nothing is significant. With two or three replicates per state this is the normal outcome for a modest effect. The right response is to report the top sets by nominal p-value as hypotheses and to add timepoints. Dropping the size filters until something crosses the line is not the right response.

Symbol drift. MSigDB set membership is keyed to a specific HGNC vintage, and gene symbols change. If a set of 200 genes only matches 140 in your list, check whether the gap is real filtering or failed symbol matching, and consider joining on Ensembl gene IDs throughout.

Direction confusion. results(dds, name = ...) and results(dds, contrast = ...) can differ in sign depending on factor level order. Set your reference level explicitly with relevel and verify the sign on a gene whose direction you already know.

Reading a pathway name as a mechanism. Set names are curator shorthand, and overlapping sets with different names often carry the same leading edge. Always compare leading-edge lists across your top results before concluding you have found several independent signals, and use network visualization when the count of significant sets exceeds what you can hold in your head ¹. If you later want to bring proteomic or metabolomic measurements into the same framework, the enrichment machinery generalizes across data types, with the main work being identifier mapping between layers ^{6 4}.

Footnotes

1. Jüri Reimand, Ruth Isserlin, Veronique Voisin, et al. Pathway enrichment analysis and visualization of omics data using g:Profiler, GSEA, Cytoscape and EnrichmentMap. Nature Protocols, 2019. https://doi.org/10.1038/s41596-018-0103-9 ↩ ↩² ↩³

2. Kangmei Zhao, Seung Yon Rhee. Interpreting omics data with pathway enrichment analysis. Trends in Genetics, 2023. https://doi.org/10.1016/j.tig.2023.01.003 ↩

3. Tianzhi Wu, Erqiang Hu, Shuangbin Xu, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. The Innovation, 2021. https://doi.org/10.1016/j.xinn.2021.100141 ↩

4. Shuangbin Xu, Erqiang Hu, Yantong Cai, et al. Using clusterProfiler to characterize multiomics data. Nature Protocols, 2024. https://doi.org/10.1038/s41596-024-01020-z ↩ ↩²

5. Mayte Suárez-Fariñas, Michelle A. Lowes, Lisa C. Zaba, et al. Evaluation of the Psoriasis Transcriptome across Different Studies by Gene Set Enrichment Analysis (GSEA). PLoS ONE, 2010. https://doi.org/10.1371/journal.pone.0010247 ↩

6. Sebastian Canzler, Jörg Hackermüller. multiGSEA: a GSEA-based pathway enrichment analysis for multi-omics data. BMC Bioinformatics, 2020. https://doi.org/10.1186/s12859-020-03910-x ↩