Pathway Analysis on Your Own RNA-Seq Data

By the end of this you will have a ranked pathway table for your own transcriptome, collapsed to a non-redundant set. Its columns are gene set name, normalized enrichment score, p-value, BH-adjusted p-value and leading-edge genes. You will also have the intermediate ranking vector you can rerun against any gene set collection. You need a counts or TPM matrix with at least four samples (more on that in step 2) and R 4.3 or newer. You also need these packages.

• tximport
• DESeq2
• limma
• edgeR
• fgsea
• msigdbr
• org.Hs.eg.db
• data.table
• ggplot2

Budget an afternoon. Most of the work is getting a defensible ranking to feed the enrichment call. The call itself takes seconds and is not where the effort goes.

1. Start from transcript-level quantification and not from a gene table someone handed you

If you have FASTQs, quantify with salmon against a decoy-aware index. Selective alignment with a genome decoy reduces spurious mappings to transcripts that are really genomic background:

# index once
grep "^>" GRCh38.primary_assembly.genome.fa | cut -d " " -f 1 | tr -d ">" > decoys.txt
cat gencode.v45.transcripts.fa GRCh38.primary_assembly.genome.fa > gentrome.fa
salmon index -t gentrome.fa -d decoys.txt -p 16 -i salmon_idx_v45 --gencode -k 31

# per sample
salmon quant -i salmon_idx_v45 -l A \
  -1 S01_R1.fastq.gz -2 S01_R2.fastq.gz \
  --validateMappings --seqBias --gcBias --posBias \
  --numGibbsSamples 20 -p 16 -o quant/S01

--gcBias matters if libraries were prepped in different batches, which they will be in a longitudinal profile. Check the mapping rate in quant/S01/logs/salmon_quant.log. For whole blood poly-A libraries expect 75-90%. Below 60%, stop and look at adapter content and rRNA/globin fraction before you do anything downstream.

Then collapse to gene level with tximport, using countsFromAbundance = "lengthScaledTPM" so that transcript-length changes between timepoints do not masquerade as expression changes:

library(tximport); library(AnnotationDbi); library(org.Hs.eg.db)
t2g <- read.delim("gencode_v45_tx2gene.tsv")           # tx_id, gene_id, gene_symbol
txi <- tximport(file.path("quant", samples$id, "quant.sf"),
                type = "salmon", tx2gene = t2g[,1:2],
                countsFromAbundance = "lengthScaledTPM")

If instead you were given a gene-level matrix, find out three things before proceeding: the annotation version, whether values are counts or TPM, and whether globin and rRNA were removed. Pathway results shift visibly with annotation version, because gene set memberships are keyed on symbols that get renamed.

2. Decide what comparison you are making

This is where personal RNA-seq diverges from the textbook. A case/control differential expression analysis of the kind used to rank pathways in tumor versus normal tissue does not exist for one person.¹ You have three workable designs.

Within-person over time. Two or more conditions you imposed on yourself, each with repeated sampling. Examples are fasted versus fed, pre- versus post-exercise, and before versus during a sleep restriction block. The samples are the replicates and you are the block. With limma-voom this is a paired design:

library(edgeR); library(limma)
y <- DGEList(counts = round(txi$counts), genes = gene_meta)
keep <- filterByExpr(y, group = meta$condition)          # ~13-16k genes survive in blood
y <- y[keep, , keep.lib.sizes = FALSE]
y <- calcNormFactors(y, method = "TMM")
design <- model.matrix(~ condition + draw_batch, data = meta)
v <- voom(y, design, plot = TRUE)
corfit <- duplicateCorrelation(v, design, block = meta$visit)
fit <- lmFit(v, design, block = meta$visit, correlation = corfit$consensus)
fit <- eBayes(fit)
res <- topTable(fit, coef = "conditionfed", number = Inf, sort.by = "none")

Trend over a continuous covariate. If you have 8+ timepoints spread over a year, model expression against something you measured continuously rather than against a categorical label. ~ hba1c + neut_frac + draw_batch gives you a per-gene slope, and the moderated t-statistic on that slope is a perfectly good ranking metric.

Deviation from a reference distribution. One sample, no comparison. You can still compute a z-score per gene against a matched external cohort (GTEx whole blood, or the Human Protein Atlas blood atlas), but the batch effect between your library prep and theirs is usually larger than the biology. We do this only for genes with very large effects (say |z| > 4) and we treat the resulting pathway list as a hypothesis generator, nothing more.

Pick one and write it down. The single most common failure in personal omics is running enrichment on a gene list whose provenance nobody can reconstruct.

3. Choose the method class deliberately

Three families, and they answer different questions.²³

Over-representation analysis (ORA) takes a thresholded gene list and asks whether a pathway is over-represented relative to a background. Fast, interpretable, and sensitive to your cutoff and your choice of background universe. Use the filtered expressed genes as the universe, never all ~20,000 annotated genes. ORA on a personal dataset with modest effect sizes tends to return nothing, because your list is short.

Functional class scoring (GSEA-style) uses the full ranked gene vector and needs no threshold. This is what we run first, and for n=1 longitudinal data it is the right default, because coordinated 15% shifts across a 60-gene set are detectable while no individual gene clears FDR 0.05.

Topology-based methods use the pathway graph, not just membership, so that a change in an upstream kinase counts differently from a change in a terminal effector. ToPASeq implements several of these (SPIA, TopologyGSA, clipper, DEGraph) over a common interface and handles RNA-seq count input directly.⁴ Coverage is the limit: topology exists for KEGG and Reactome-derived graphs. It does not exist for the MSigDB hallmark or immunologic collections. We run topology as a second pass on the 10-20 pathways that FCS flags, to ask where in the pathway the signal sits.

4. Build the ranking vector, then run fgsea

Rank by the moderated t-statistic (limma) or the shrunken Wald statistic (DESeq2 with lfcShrink(type="apeglm")). Do not rank by raw log fold change: low-count genes generate enormous unreliable fold changes and they will dominate the tails of your ranking. Do not rank by p-value alone, since it is unsigned.

library(fgsea); library(msigdbr); library(data.table)

stats <- res$t
names(stats) <- res$gene_symbol
stats <- stats[!is.na(names(stats)) & !duplicated(names(stats))]
stats <- sort(stats, decreasing = TRUE)

h  <- msigdbr(species = "Homo sapiens", collection = "H")
r  <- msigdbr(species = "Homo sapiens", collection = "C2", subcollection = "CP:REACTOME")
sets <- split(c(h$gene_symbol, r$gene_symbol),
              c(h$gs_name, r$gs_name))

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

minSize = 15 drops sets too small to give a stable enrichment score. maxSize = 500 drops sets so broad they are unfalsifiable (“metabolic process”). eps = 0 lets fgseaMultilevel estimate p-values below the multilevel floor, which matters when your top hits are strong. Start with hallmark (50 sets, curated to be non-redundant) before Reactome (~1,600 sets, deeply nested). If hallmark returns nothing at FDR 0.10, Reactome will return a long list of overlapping parent and child terms that all reflect the same handful of genes.

5. Collapse redundancy and read the leading edge

A raw Reactome result will show “Immune System”, “Cytokine Signaling in Immune System”, “Interferon Signaling”, and “Interferon alpha/beta signaling” as four hits. They are one finding. collapsePathways fits the parent sets and keeps only those with independent signal:

main <- collapsePathways(fg[padj < 0.1][order(pval)], sets, stats)
fg_indep <- fg[pathway %in% main$mainPathways]

# what genes are driving it
fg_indep[1, leadingEdge][[1]]
plotEnrichment(sets[[fg_indep$pathway[1]]], stats) 

Always read the leading edge before you believe a pathway. Two checks. First, how many genes: if a 300-gene set is enriched on the strength of 9 leading-edge genes, the label on the pathway is doing more work than the data. Second, what those genes are: in whole blood, a leading edge of HBB, HBA1, HBA2, ALAS2, SLC4A1 means you measured reticulocyte content. A leading edge of DEFA3, DEFA4, LCN2, CAMP, MPO means you measured neutrophil degranulation from a different draw time. Pathway enrichment inherits every bias in the upstream gene ranking, and annotation bias toward well-studied genes is real and directional.³

6. Regress out cell composition, because in blood it is the dominant signal

Whole blood bulk RNA-seq is a weighted average of cell types. Neutrophil fraction varies 30-70% across draws in a healthy person and drives the largest principal component in most personal datasets. If you do not model it, half your pathway hits will be immune cell proportion.

Estimate fractions from the expression matrix itself, then include them as covariates:

# xCell or CIBERSORTx on log2(TPM+1); or use the CBC differential from the same draw
frac <- read.delim("cibersortx_abs.txt", row.names = 1)
meta$neut  <- frac[meta$id, "Neutrophils"]
meta$lymph <- frac[meta$id, "T.cells.CD8"] + frac[meta$id, "T.cells.CD4.naive"]

design <- model.matrix(~ condition + neut + lymph + draw_batch, data = meta)

If you have a same-draw CBC with differential, use the measured fractions instead of deconvolved ones. They are better. Run the analysis both ways and look at what changes. Pathways that survive adjustment for cell composition are the interesting ones. Pathways that vanish were composition.

7. Cross-check with a second tool, and use web tools for that

Method choice changes the answer. Different enrichment algorithms on the same data return overlapping but non-identical lists, so a hit that appears under FCS, ORA, and a topology method is worth more than a hit that appears under one. The web tools are good for this cross-check.

The iDEP tool will take a raw counts matrix and run the whole chain against a large collection of annotation databases. The chain is filtering, normalization and PCA, then DESeq2 or limma, then GSEA/PGSEA and pathway visualization. That makes it a fast independent replication of your R pipeline.⁵ The ReactomeGSA tool runs the analysis server-side over Reactome’s curated pathways and is built for comparative multi-omics, so the same call handles RNA-seq and proteomics.⁶ The GeneAnalytics tool scores gene sets against tissue and cell-type contexts as well as pathways, which is useful when you suspect a composition artifact.⁷ The ShinyGO tool is the fastest ORA sanity check if you just want to paste 200 symbols.

Our practice: R pipeline is the record, web tools are the second opinion. Web tools hide their filtering and background-universe choices, and those choices change ORA results substantially.

8. Bring in the other layers

If you have proteomics and metabolomics on the same timepoints, joint pathway analysis is stronger than three separate analyses, mostly because concordance across layers is a much harder thing to produce by chance than a single enriched list.⁸⁹ Two tools do this well without asking you to write integration code. The first is PaintOmics 3. It accepts transcriptomics, proteomics, and metabolomics matrices and maps them onto KEGG and Reactome pathway diagrams. It colors each node by layer, so you see where a transcript change is and is not matched by a protein change.¹⁰ The second is ReactomeGSA, which handles multiple omics as separate datasets in one comparative analysis.⁶

Expect transcript-protein discordance. Plasma proteomics measures a secreted compartment, not the blood cell transcriptome, and correlation between the two at the single-analyte level is generally weak. A pathway that moves in both layers is a real finding. A pathway that moves only in RNA is a finding about RNA.

9. Freeze the analysis

Record the gene set collection version (MSigDB releases change set membership), the annotation GTF version, the R package versions, and your random seed. Write the ranking vector to disk as a two-column TSV. When a new gene set collection ships in six months, you rerun enrichment in thirty seconds without re-deriving anything upstream.

saveRDS(list(stats = stats, res = res, fgsea = fg, sessionInfo = sessionInfo()),
        "pathway_run_2026-09-15.rds")

Common problems

Nothing is significant. Usually correct. Personal longitudinal data with 4-8 samples and modest perturbations often has no pathway at FDR 0.05. Report the top 10 with nominal p and NES and treat them as hypotheses for the next draw. Lowering the FDR threshold until something appears is how people end up believing in noise.

Everything is significant. Almost always cell composition, library batch, or both. Check whether PC1 of the log-CPM matrix correlates with neutrophil fraction or with prep date. If it does, step 6.

The same pathway appears under six names. Run collapsePathways, or restrict to hallmark. Nested collections like Reactome and GO produce parent/child duplicates by construction.

Gene symbols do not match. Roughly 5-10% of symbols in an older matrix will be aliases or withdrawn. Map through Ensembl gene IDs with org.Hs.eg.db or biomaRt, keep the ID as the key, and attach symbols only at the end. Unmapped genes silently vanish from every gene set and they bias the result in whatever direction the missing genes were going.

Duplicate symbols in the ranking vector. fgsea will take the vector but the result is undefined. Aggregate first, by keeping the entry with the highest mean expression. Do not average statistics.

You want to know what a pathway result means for you. Pathway enrichment on one person’s transcriptome is a measurement and a hypothesis. It is not a clinical finding, it does not establish causation, and it cannot tell you whether anything should change. An enriched interferon signature has dozens of possible sources, several of them benign and several not, and distinguishing them needs a clinician with your history and the ability to order confirmatory tests. Take the table to them, ask what would confirm or rule it out, and do not act on a gene set name.

Footnotes

1. A. Siavoshi, M. Taghizadeh, E. Dookhe, et al. Gene expression profiles and pathway enrichment analysis to identification of differentially expressed gene and signaling pathways in epithelial ovarian cancer based on high-throughput RNA-seq data. Genomics, 2022. https://doi.org/10.1016/j.ygeno.2021.11.031 ↩

2. Jiung-Wen Chen, Lisa Shrestha, George Green, et al. The hitchhikers’ guide to RNA sequencing and functional analysis. Briefings in Bioinformatics, 2023. https://doi.org/10.1093/bib/bbac529 ↩

3. 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 ↩ ↩²

4. Ivana Ihnatova, Eva Budinska. ToPASeq: an R package for topology-based pathway analysis of microarray and RNA-Seq data. BMC Bioinformatics, 2015. https://doi.org/10.1186/s12859-015-0763-1 ↩

5. Steven Xijin Ge, Eun Wo Son, Runan Yao. iDEP: an integrated web application for differential expression and pathway analysis of RNA-Seq data. BMC Bioinformatics, 2018. https://doi.org/10.1186/s12859-018-2486-6 ↩

6. 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 ↩ ↩²

7. Shani Ben-Ari Fuchs, Iris Lieder, Gil Stelzer, et al. GeneAnalytics: An Integrative Gene Set Analysis Tool for Next Generation Sequencing, RNAseq and Microarray Data. OMICS: A Journal of Integrative Biology, 2016. https://doi.org/10.1089/omi.2015.0168 ↩

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

9. Upasna Srivastava, Swarna Kanchan, Minu Kesheri, et al. Integrative omics approaches for identification of biomarkers. Integrative Omics, 2024. https://doi.org/10.1016/b978-0-443-16092-9.00010-2 ↩

10. Rafael Hernández-de-Diego, Sonia Tarazona, Carlos Martínez-Mira, et al. PaintOmics 3: a web resource for the pathway analysis and visualization of multi-omics data. Nucleic Acids Research, 2018. https://doi.org/10.1093/nar/gky466 ↩