How to Run Pathway Analysis on Your Own RNA-Seq Data

By the end of this guide you will have a ranked table of gene sets with normalized enrichment scores and false discovery rates, a redundancy-collapsed version of that table you can read without drowning in forty near-identical GO terms, and a network figure that groups the surviving sets into themes. You need three things to start. They are a gene-level counts matrix (or salmon/kallisto quantification directories), a sample metadata table that encodes the comparison you care about, and R 4.4 or later. R needs the following packages installed.

• DESeq2
• tximport
• fgsea
• msigdbr
• clusterProfiler
• gprofiler2

Everything below assumes human data, Ensembl gene identifiers, and that you are analyzing your own longitudinal samples rather than a case-control cohort. That last assumption changes the statistics more than anything else, so we start there.

1. Build a clean gene-level counts matrix

Pathway analysis inherits every artifact of the quantification step, so it is worth getting this right before you touch a gene set database. We use salmon with selective alignment against a decoy-aware index built from the transcriptome plus the full genome, which suppresses spurious assignment of intronic and intergenic reads to transcripts.

salmon index -t gentrome.fa.gz -d decoys.txt -i salmon_index_k31 -k 31 --gencode

salmon quant -i salmon_index_k31 -l A \
  -1 sample01_R1.fastq.gz -2 sample01_R2.fastq.gz \
  --validateMappings --seqBias --gcBias --posBias \
  --numGibbsSamples 20 -p 16 -o quants/sample01

Then import to gene level with tximport, which corrects for average transcript length differences across samples. Length correction matters: if an isoform switch shortens the effective length of a gene in one condition, raw counts will move even when molar abundance does not.

library(tximport); library(DESeq2)
tx2gene <- read.csv("tx2gene_gencode_v46.csv")   # tx_id, gene_id
files <- file.path("quants", meta$sample, "quant.sf"); names(files) <- meta$sample
txi <- tximport(files, type = "salmon", tx2gene = tx2gene,
                countsFromAbundance = "lengthScaledTPM")

dds <- DESeqDataSetFromTximport(txi, colData = meta, design = ~ timepoint)
keep <- rowSums(counts(dds) >= 10) >= ceiling(0.5 * ncol(dds))
dds <- dds[keep, ]

That filter typically leaves 14,000 to 17,000 genes from a whole-blood or PBMC library sequenced to 30 to 50 million reads. Do not skip it. Genes detected in two of twelve samples contribute noise to the ranking and, worse, they inflate the apparent size of the measured universe when you later run an overrepresentation test.

One caveat specific to blood: globin and immunoglobulin transcripts can consume 20 to 70 percent of reads in whole blood without globin depletion, which compresses everything else. Check assays(dds)$counts["ENSG00000244734", ] (HBB) as a fraction of library size before interpreting anything about immune pathways.

2. Decide what comparison you are making

Pathway analysis answers a question about a contrast, never about a sample in isolation. With a single person’s longitudinal data you have three workable designs, and they are not interchangeable.

The first is a paired within-person contrast: two or more timepoints, with a covariate for anything that changed methodologically. If you have replicate draws at each timepoint, design = ~ batch + timepoint in DESeq2 works directly. With a single library per timepoint and only two timepoints, DESeq2 cannot estimate dispersion from replication, and you should not pretend otherwise. You can supply a dispersion prior from a public dataset of the same tissue and library prep, but treat the resulting p-values as heuristics for ordering genes rather than as calibrated error rates.

The second is a trend design over many timepoints, which is the strongest option if you have six or more draws. Fit ~ ns(days, df = 3) with a natural spline and use DESeq2::results() on a likelihood-ratio test against the intercept-only model. This finds genes with reproducible temporal structure and is far less sensitive to a single aberrant draw than any two-point comparison.

The third is a single-sample scoring approach, where instead of a contrast you compute a per-sample pathway activity score with GSVA::gsva() using the gsvaParam or ssgseaParam constructors, then track those scores across your own timepoints. This is how we would read a personal profile: each pathway becomes a time series, and you look for coherent drift rather than a one-off ranking. GSVA scores are relative to the sample set you pass in, so they are only interpretable within one normalized matrix, and adding new timepoints later shifts old scores. Re-run the whole matrix each time.

Whatever you choose, write the contrast down before you look at results. Pathway databases are large enough that a post hoc comparison will always produce something.

3. Build a ranking rather than a cut list

The single most common mistake in this workflow is thresholding at adjusted p < 0.05 and feeding the surviving gene symbols to a web tool. Threshold-based overrepresentation throws away the magnitude and direction of every measurement and makes the result a function of how many genes happened to clear an arbitrary cutoff, which in turn depends on your sequencing depth. Threshold-free methods that take the entire ranked gene list are the better default, and they are what published protocols recommend when a full ranking is available ¹.

Rank by the Wald statistic, or by shrunken log fold change if you want magnitude to dominate:

res <- lfcShrink(dds, coef = "timepoint_T2_vs_T1", type = "apeglm")
rk <- res$stat                       # or: res$log2FoldChange after shrinkage
names(rk) <- rownames(res)
rk <- rk[!is.na(rk)]
rk <- sort(rk, decreasing = TRUE)

Two things to avoid. Do not rank by -log10(p) * sign(lfc) unless you understand that it makes very high-expression genes dominate, since p-values scale with count depth. And do not rank by fold change alone without shrinkage, because a gene going from 3 counts to 11 will outrank a gene going from 4,000 to 9,000.

4. Choose gene set collections deliberately

The three collections people reach for answer different questions, and mixing them into one multiple-testing family is a mistake.

Gene Ontology (GO) is a controlled vocabulary of three orthogonal namespaces: biological process, molecular function, and cellular component. It is hierarchical, so a gene annotated to “positive regulation of T cell activation” is implicitly annotated to every ancestor term up to “biological process”. That structure is why GO enrichment output is so redundant: you will see the same forty genes drive fifteen nested terms. GO is broad, well-annotated, and mostly useful as a first pass or a fishing expedition. Interpret a GO result as a direction to investigate. It is not a mechanism.

KEGG PATHWAY is a set of hand-drawn metabolic and signaling maps, on the order of a few hundred for human. Maps are large and coarse, and they sometimes mix a disease name with a mechanism. So “KEGG: Pathways in cancer” enriching in your data tells you very little. KEGG’s strength is metabolism: the glycolysis, TCA, fatty acid oxidation, and one-carbon maps are curated at reaction-level detail and align well with metabolite data if you have it. Its API terms restrict bulk redistribution, which is why some tools ship Reactome and WikiPathways but require you to fetch KEGG yourself.

Reactome is a knowledgebase that is reaction-level, hierarchically organized and manually curated. Its pathways are trees of events with defined inputs, outputs, catalysts, and compartments. They span several thousand human reactions organized into thousands of pathways and subpathways ². Because it is reaction-level, you can test a specific step (“Activation of NOXA and translocation to mitochondria”) rather than an entire signaling cascade, and because it is a hierarchy you can aggregate results to whatever resolution you want. Reactome’s server-side analysis tooling also supports multi-omics comparison directly, which is convenient when you have proteomics alongside RNA ³.

Our recommendation is to run MSigDB Hallmark as the primary analysis, with its 50 sets deliberately deduplicated and curated to represent coherent biological states. Use Reactome as the secondary collection for mechanistic resolution and treat GO biological process as exploratory. Correct for multiple testing within each collection separately and report them separately.

library(msigdbr)
hallmark <- msigdbr(species = "Homo sapiens", collection = "H")
reactome <- msigdbr(species = "Homo sapiens", collection = "C2",
                    subcollection = "CP:REACTOME")
hm  <- split(hallmark$ensembl_gene, hallmark$gs_name)
rct <- split(reactome$ensembl_gene, reactome$gs_name)

5. Run fgsea

fgsea implements a multilevel Monte Carlo scheme that gives accurate p-values far below what a fixed permutation count can resolve, which matters because the strongest Hallmark set in a real dataset frequently lands below 1e-6.

library(fgsea)
set.seed(42)
fg <- fgseaMultilevel(pathways = hm, stats = rk,
                      minSize = 15, maxSize = 500, eps = 0)
fg <- fg[order(fg$padj), ]

Read four columns.

• NES (normalized enrichment score, the size-corrected effect)
• padj (Benjamini-Hochberg FDR)
• size (genes from the set present in your ranking)
• leadingEdge (the subset of genes carrying the signal)

The leading edge is where interpretation happens. A Hallmark interferon alpha response set with NES 2.4, padj 3e-8, and a leading edge of 31 genes dominated by IFIT1, IFIT3, ISG15, OAS1, and MX1 is a coherent result. The same NES with a leading edge of 4 genes out of 190 is a set whose signal comes from a handful of outliers, and you should plot those genes’ counts before believing it.

For Reactome, run the same call on rct and then collapse the hierarchy, because parent and child pathways will both appear:

fgR <- fgseaMultilevel(rct, rk, minSize = 15, maxSize = 500, eps = 0)
keep <- collapsePathways(fgR[fgR$padj < 0.05, ][order(pval), ], rct, rk)
fgR_main <- fgR[fgR$pathway %in% keep$mainPathways, ]

collapsePathways removes sets whose enrichment is explained by a more significant overlapping set. It routinely cuts a 180-row Reactome result to 25 independent themes, which is the difference between a readable table and a dump.

6. If you must run overrepresentation, fix the background

Sometimes you only have a gene list: a public supplementary table, or a set of genes carrying rare variants from your whole-genome data. Overrepresentation analysis is then appropriate, but the background set determines the answer. The correct universe is the set of genes that could have been detected in your experiment, which is exactly the filtered rows of your dds object, not all ~20,000 protein-coding genes and certainly not the tool’s default of “all genes in the annotation”. Using an inflated background systematically manufactures significance, a point the standard protocol makes explicitly ¹.

library(clusterProfiler)
universe <- rownames(dds)                       # post-filter, detected genes
sig <- rownames(res)[which(res$padj < 0.05 & abs(res$log2FoldChange) > 0.5)]

ego <- enrichGO(gene = sig, universe = universe, OrgDb = org.Hs.eg.db,
                keyType = "ENSEMBL", ont = "BP",
                pAdjustMethod = "BH", qvalueCutoff = 0.05, readable = TRUE)
ego <- simplify(ego, cutoff = 0.7, by = "p.adjust", measure = "Wang")

With g:Profiler, the equivalent is to pass your detected genes as a custom background and set the domain scope accordingly. Its g:SCS correction accounts for the dependency structure among nested GO terms, which Benjamini-Hochberg does not.

library(gprofiler2)
g <- gost(query = sig, organism = "hsapiens",
          custom_bg = universe, domain_scope = "custom_annotated",
          sources = c("GO:BP", "REAC", "WP"),
          correction_method = "gSCS", evcodes = TRUE)

7. Collapse redundancy and draw the map

A pathway table with 60 significant rows is not a result anyone can hold in their head. The fix is to build a similarity network where nodes are enriched sets and edges are gene overlap, then cluster it. EnrichmentMap in Cytoscape does this well: import the fgsea or g:Profiler output with the gene set GMT, set the edge cutoff to a combined Jaccard/overlap coefficient around 0.375, and run AutoAnnotate to label each cluster with a summary phrase ¹. The output is typically five to fifteen labeled themes, each of which you can then interrogate at gene level.

We consider this step mandatory rather than cosmetic. Pathway databases overlap heavily by construction, and reviewing a clustered map is what prevents you from reporting one biological signal fifteen times ⁴.

8. Cross-check against your other assays

A pathway result from RNA alone is a hypothesis about transcript abundance. Protein levels, metabolites, and biomarker panels frequently disagree. The disagreement is informative rather than a failure. If you have plasma proteomics from the same draw, run the same contrast through ReactomeGSA, which accepts multiple omics layers and reports pathway-level results in a common coordinate system so you can see where transcript and protein move together ³. Broader integrative tools follow the same logic across additional data types and pathway databases ⁵. The methodological tradeoffs of these integration strategies are surveyed in detail ⁶, particularly how to weight layers with very different coverage and noise. The interpretive pitfalls of multi-omic pathway output are worth reading before you draw conclusions ⁷.

Concretely, suppose the Hallmark interferon response set enriches in RNA and your proteomics shows no movement in the corresponding secreted proteins. The most likely explanations are timing (transcript precedes protein by hours to days), cell composition shift in blood rather than per-cell regulation, or noise. Deconvolution of your bulk counts with a reference signature matrix will often resolve the first two.

9. Pin everything

Pathway databases change monthly, and the same ranked list run against MSigDB 7.5 and 2024.1 will not give the same table. Record the collection version, the msigdbr and fgsea package versions, and your random seed. Record the gene filter and the exact contrast too, in a sessionInfo() dump committed alongside the results. For longitudinal self-tracking this is not optional, because your comparison next year is against this year’s output.

Common problems

Everything is significant. Usually a design problem: unmodeled batch, a sex mismatch in metadata, or degraded RNA in one sample driving a global 3’ bias. Check the PCA of vst(dds) first. If PC1 separates by extraction date rather than by your variable of interest, fix the model before touching gene sets.

Nothing is significant. With n=1 longitudinal data and a genuine but modest effect, this is the expected outcome, and it is the correct answer rather than a bug. Look at the unadjusted NES ranking for coherence across related sets, treat it as a hypothesis, and add timepoints.

Identifier loss. Mapping Ensembl to Entrez or symbol silently drops 10 to 25 percent of genes and collapses duplicates arbitrarily. Keep Ensembl IDs end to end where the tool supports it, and always print how many of your ranked genes matched the gene set universe.

Ribosomal, mitochondrial, and cell-cycle sets always on top. These reflect library composition and proliferating cell fraction more often than regulation. Inspect whether the leading edge is broad and modest (composition) or narrow and large (regulation).

Pathway names that read like a diagnosis. Many curated sets are named after diseases because they were assembled from disease literature. Enrichment of a set named for a condition is a statement about gene overlap. It is not a statement about you. Network and pathway methods were developed partly to make such coarse annotations interpretable in context ⁸, and cancer-cohort resources make clear how much of this signal is tissue and cohort specific ⁹. If a result seems clinically meaningful, that is the point at which to bring it to a physician or genetic counselor rather than to act on it yourself.

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. Erik D. Huckvale, Hunter N. B. Moseley. Predicting the Pathway Involvement of Compounds Annotated in the Reactome Knowledgebase. Metabolites, 2025. https://doi.org/10.3390/metabo15030161 ↩

3. Johannes Griss, Guilherme Viteri, Konstantinos Sidiropoulos, et al. ReactomeGSA - Efficient Multi-Omics Comparative Pathway Analysis. Molecular & Cellular Proteomics, 2020. https://doi.org/10.1074/mcp.tir120.002155 ↩ ↩²

4. 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 ↩

5. Tianyuan Liu, Pedro Salguero, Marko Petek, et al. PaintOmics 4: new tools for the integrative analysis of multi-omics datasets supported by multiple pathway databases. Nucleic Acids Research, 2022. https://doi.org/10.1093/nar/gkac352 ↩

6. Zeynab Maghsoudi, Ha Nguyen, Alireza Tavakkoli, et al. A comprehensive survey of the approaches for pathway analysis using multi-omics data integration. Briefings in Bioinformatics, 2022. https://doi.org/10.1093/bib/bbac435 ↩

7. William G. Ryan V., Smita Sahay, John Vergis, et al. Pathway Analysis Interpretation in the Multi-Omic Era. BioTech, 2025. https://doi.org/10.3390/biotech14030058 ↩

8. Gal Barel, Ralf Herwig. Network and Pathway Analysis of Toxicogenomics Data. Frontiers in Genetics, 2018. https://doi.org/10.3389/fgene.2018.00484 ↩

9. Suhas V Vasaikar, Peter Straub, Jing Wang, et al. LinkedOmics: analyzing multi-omics data within and across 32 cancer types. Nucleic Acids Research, 2017. https://doi.org/10.1093/nar/gkx1090 ↩