RNA-seq Pathway Analysis in R, Start to Finish

Oak
rna-seq pathway-analysis fgsea deseq2 bioinformatics
A laboratory machine sorting glowing strands into thick branching bundles along a lit rail, some bundles brighter than others.

By the end of this you will have a table with one row per gene set: normalized enrichment score, BH-adjusted p-value, set size, and the leading-edge genes that drove the score, collapsed so that you are not reading the same biology fifteen times under fifteen GO labels. You need R ≥ 4.3 with DESeq2, tximport, fgsea, msigdbr, edgeR, and data.table, plus either salmon quantification directories or a gene-level count matrix someone handed you. You also need a real contrast. That last requirement is where most personal RNA-seq projects fail, so it gets its own step before any statistics.

1. Quantify, and keep the bootstraps

If you have FASTQs, quantify with salmon against a GENCODE transcriptome (the decoy-aware index, built from the transcriptome plus the full genome as decoy):

salmon index -t gentrome.fa.gz -d decoys.txt -i gencode_v44_idx -k 31 -p 16

salmon quant -i gencode_v44_idx -l A \
  -1 S01_R1.fastq.gz -2 S01_R2.fastq.gz \
  --validateMappings --seqBias --gcBias --posBias \
  --numBootstraps 30 -p 16 -o quants/S01

--gcBias matters if libraries were prepped in different batches or on different days, which is the normal case for a longitudinal series. --numBootstraps 30 costs minutes and gives you per-transcript inferential variance, which is the only way to know whether a gene-level count is dominated by ambiguous multi-isoform assignment.

Check the mapping rate in quants/S01/logs/salmon_quant.log. For a good polyA library on human, expect roughly 70–90% mapping. Below 60%, stop and look at your reads: rRNA carryover, adapter contamination, or the wrong strandedness setting will all show up here before they show up as weird pathways. Target 30–50M paired-end reads per sample for differential expression on a bulk tissue or blood library. Turnaround from a drawn tube to FASTQs is typically one to three weeks, mostly queueing.

Import to gene level:

library(tximport); library(readr)

tx2gene <- read_tsv("gencode_v44_tx2gene.tsv",
                    col_names = c("tx_id","gene_id","gene_name"))

files <- file.path("quants", samples$id, "quant.sf")
names(files) <- samples$id

txi <- tximport(files, type = "salmon", tx2gene = tx2gene[,1:2],
                countsFromAbundance = "lengthScaledTPM")

lengthScaledTPM gives you counts that are already corrected for the average transcript length per gene per sample, so downstream count models behave and you do not have to carry an offset matrix around. If you prefer the offset route, use countsFromAbundance = "no" and let DESeqDataSetFromTximport build the normalization factors.

2. Decide what you are comparing against

Transcript abundance is a snapshot of one tissue at one moment, and it is the most state-dependent layer in an omics stack: time of day, last night’s sleep, a viral exposure two days ago, and whether you ran that morning all move it 1. So the contrast is the whole experiment.

For a single person you have three honest options.

Longitudinal within-subject. You have n timepoints of the same tissue and you want genes that move with an exposure, a season, or an intervention window. This is the strongest design available to an individual, because you are your own control and the genotype term drops out. You need at least five or six timepoints before per-gene variance estimates mean anything, and ideally a few baseline draws before whatever you are tracking.

Against a reference cohort. GTEx v8 whole blood or a public PAXgene cohort. This tells you where your sample sits in a population distribution. It is contaminated by batch: different library kit, different rRNA/globin depletion, different sequencer. Expect the largest principal component to be “study”, not “you.” Use it for ranking within your own sample (a z-score per gene against the cohort mean and SD) and treat the absolute values as unreliable.

Paired conditions. Two tissues, or blood before and after a defined perturbation on the same day, run in the same library prep batch. Cleanest signal, narrowest question.

Whatever you pick, write the design formula down before you look at any results.

3. Filter, then fit

library(DESeq2); library(edgeR)

dds <- DESeqDataSetFromTximport(txi, colData = samples,
                                design = ~ batch + condition)

keep <- filterByExpr(counts(dds), group = samples$condition)
dds  <- dds[keep, ]
nrow(dds)   # expect ~13,000-16,000 for whole blood, ~15,000-18,000 for solid tissue

dds <- DESeq(dds)
res <- results(dds, contrast = c("condition", "post", "pre"))

filterByExpr is better than an arbitrary rowSums(counts) > 10: it scales the CPM threshold to library size and requires the gene to pass in at least as many samples as the smallest group. Filtering before fitting is not cosmetic. It raises power at fixed FDR by shrinking the multiple-testing burden, and it removes the near-zero genes whose log fold changes are pure noise and which will otherwise sit at the extremes of your ranking vector.

Then look at the PCA on variance-stabilized counts before you believe anything:

vsd <- vst(dds, blind = TRUE)
plotPCA(vsd, intgroup = c("condition","batch"))

If PC1 separates by batch and PC2 by nothing you recognize, fix that first. ComBat_seq from sva is the pragmatic tool when batch and condition are not confounded. If they are confounded, nothing will save you.

4. Build the ranking vector

For GSEA-style analysis you rank every gene, not just the significant ones. Use the Wald statistic from DESeq2, which already encodes both effect size and precision:

library(dplyr); library(tibble)

ranks <- res |>
  as.data.frame() |>
  rownames_to_column("gene_id") |>
  left_join(distinct(tx2gene, gene_id, gene_name), by = "gene_id") |>
  filter(!is.na(stat), !is.na(gene_name)) |>
  group_by(gene_name) |>
  slice_max(abs(stat), n = 1, with_ties = FALSE) |>   # collapse duplicate symbols
  ungroup()

stats <- setNames(ranks$stat, ranks$gene_name)
stats <- sort(stats, decreasing = TRUE)

Do not rank by raw log2 fold change. A gene with counts of 3 and 11 has a fold change of 1.87 and no information. If you want an effect-size ranking, use lfcShrink(dds, coef = ..., type = "apeglm") and rank on the shrunken estimate. Do not rank by p-value alone, because it throws away direction.

For the N=1 reference-cohort case there is no Wald statistic. Rank on a per-gene modified z-score against the cohort:

z <- (my_vst - matrixStats::rowMedians(ref_vst)) /
     (matrixStats::rowMads(ref_vst) + 0.1)

The + 0.1 floor stops low-variance genes from producing z-scores of 40. Use the median and MAD rather than mean and SD: public cohorts contain outlier samples, and a single bad one will distort a mean-based score for thousands of genes.

Sanity-check the ranking with a positive control you already know. XIST and RPS4Y1 should sit where your chromosomal sex says they should. If your sample is post-exercise blood, neutrophil genes (FCGR3B, CSF3R, S100A8) should be near the top. If they are not, your labels are swapped.

5. Get gene sets, and pick the right collections

library(msigdbr)

# msigdbr >= 10 uses `collection`/`subcollection`; earlier versions use `category`/`subcategory`
h  <- msigdbr(species = "Homo sapiens", collection = "H")
c2 <- msigdbr(species = "Homo sapiens", collection = "C2", subcollection = "CP:REACTOME")
c5 <- msigdbr(species = "Homo sapiens", collection = "C5", subcollection = "GO:BP")

to_list <- function(df) split(df$gene_symbol, df$gs_name)
sets <- c(to_list(h), to_list(c2))

Our default is Hallmark plus Reactome. Hallmark’s 50 sets are curated to be non-redundant and each one is large enough to have stable statistics, so it is the right first pass. Reactome adds mechanistic resolution with a sane hierarchy. GO Biological Process is the most commonly used and the most painful: about 7,000 testable sets, many nested, so a single real signal produces dozens of correlated “hits” and the BH correction is applied to a set of tests that are nowhere near independent. Run GO BP second, and always collapse it (step 7).

6. Run fgsea

library(fgsea)

set.seed(42)
fg <- fgseaMultilevel(pathways = sets,
                      stats    = stats,
                      minSize  = 15,
                      maxSize  = 500,
                      eps      = 0,
                      nPermSimple = 10000)

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

minSize = 15 drops sets too small for a meaningful enrichment score. maxSize = 500 drops sets so broad (“metabolic process”) that the result tells you nothing. eps = 0 turns off the lower bound on reported p-values so that strongly enriched sets get an actual estimate instead of being floored, at the cost of runtime. fgseaMultilevel is the current algorithm and is what fgsea() dispatches to by default.

Read NES, not p-value, for effect size. NES is the enrichment score normalized by set size, so it is comparable across sets. A well-powered experiment with 4 vs 4 samples typically yields Hallmark NES values in the 1.5–2.5 range for real signal. Anything above 3 with a tiny set size deserves a look at the leading edge before you believe it.

Over-representation analysis (a hypergeometric test on your padj < 0.05 gene list, via clusterProfiler::enrichGO) is the other common approach. We prefer GSEA because it does not require a significance cutoff, uses the whole ranking, and detects coordinated small shifts across a pathway that never clear per-gene FDR. ORA is faster and easier to explain, and it is the right tool when you have a gene list from somewhere else with no ranking attached. If you run ORA, your background universe must be the genes that survived filterByExpr, not all 60,000 GENCODE entries. Using the wrong universe is the single most common way to produce a confidently wrong ORA result.

7. Collapse redundancy and read the leading edge

main <- collapsePathways(fg[padj < 0.05][order(pval)],
                         pathways = sets, stats = stats)

fg_main <- fg[pathway %in% main$mainPathways]

fg_main[1:10, .(pathway, NES, padj,
                lead = sapply(leadingEdge, \(x) paste(head(x, 8), collapse = ",")))]

collapsePathways re-runs the enrichment test conditioned on the genes already explained by a stronger pathway, and keeps the parents. On a GO BP run this typically cuts 300 significant sets to 20–40 independent ones.

Now read the leading edge. A pathway hit whose leading edge is five ribosomal proteins and three mitochondrial genes is a library-composition artifact. A hit whose leading edge is the canonical members of the pathway, moving in a consistent direction, is a result. This step is manual and it is the part that cannot be automated away.

Save the table as TSV, but flatten the list column first, or fwrite will fail:

fg_out <- copy(fg_main)
fg_out[, leadingEdge := sapply(leadingEdge, paste, collapse = ";")]
data.table::fwrite(fg_out, "pathways.tsv", sep = "\t")

8. Test the confounders before you interpret

Bulk RNA-seq of blood measures a mixture. If your neutrophil-to-lymphocyte ratio shifted between draws, every neutrophil-associated pathway will light up, and the enrichment is real but it reflects cell proportions rather than regulation inside any cell. Single-cell and spatial methods exist precisely because a bulk average hides which population changed 2. You can approximate the correction with deconvolution:

library(granulator)   # or immunedeconv, or run CIBERSORTx
props <- deconvolute(m = as.matrix(assay(vsd)), sigMatrix = sigList$LM22)

Then either add the top two or three estimated proportions as covariates in the DESeq2 design, or check whether your top pathways correlate with them. If a pathway’s NES tracks neutrophil fraction across samples at r > 0.8, say so in your notes.

Globin is the other whole-blood specific problem. HBB, HBA1, HBA2 and ALAS2 can dominate an un-depleted whole-blood library and eat most of your usable depth. Check their TPM share before anything else. Different globin depletion chemistry between two of your timepoints will produce a fake pathway result every time.

Gene length bias affects ORA more than GSEA: longer genes accumulate more reads and are more likely to reach significance, and pathways enriched for long genes (extracellular matrix, synaptic signaling) therefore come up spuriously. goseq::nullp models this with a probability weighting function over median transcript length and is worth running if your headline result is length-biased.

9. Treat a pathway hit as a hypothesis for the next layer

A transcript-level enrichment says mRNA moved. It does not say protein moved, and mRNA-protein correlation across genes is modest. Studies that follow a transcriptional signal into targeted proteomics routinely find that only a subset of the pathway holds up, and the mediator that matters is sometimes a protein whose transcript barely moved 3. The useful next step is to check whether the leading-edge genes have measurable protein or metabolite correlates in your other data, and whether they move together over time. Integrating transcript and metabolite layers into a single network recovers regulators that a ranked gene list alone does not surface 4. Environmental and dietary exposures show up across layers too, which is why an unexplained pathway shift is worth cross-checking against what changed in your inputs during that window 5.

None of this output is a diagnosis. A pathway table is a description of a measurement, not a clinical finding. If something here looks medically relevant to you, take the underlying numbers to a physician and have them ordered as validated clinical assays.

Common problems

Everything is significant. Usually batch confounded with condition, or you forgot to filter low-count genes. Check the PCA colored by prep date.

Nothing is significant. With n = 3 per group, this is the expected outcome for a subtle effect. Report NES and rank rather than pretending to a p-value, and do not lower the FDR threshold to 0.25 without saying so.

fgsea warns about ties in the ranking. You have duplicate statistics, usually because genes with identical zero counts survived filtering. Re-check your filter.

Duplicate gene symbols. Ensembl IDs map many-to-one onto symbols. Collapse by maximum absolute statistic (as in step 4) or keep Ensembl IDs throughout and use an Ensembl-keyed gene set list. Do not let setNames silently overwrite.

Results change when you re-run. Set set.seed(). fgseaMultilevel is stochastic, and p-values near the threshold will move between runs.

The top pathway is “KEGG Parkinson’s disease” in your blood sample. Many disease-named KEGG sets are dominated by oxidative phosphorylation and ribosomal genes. Read the leading edge. It is a mitochondrial signal wearing a disease label.

Oak builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

  1. M. Vailati-Riboni, Valentino Palombo, Juan J. Loor. What Are Omics Sciences?. 2017. https://doi.org/10.1007/978-3-319-43033-1_1 ↩

  2. Sabrina M. Lewis, Marie-Liesse Asselin-Labat, Quan Nguyen, et al. Spatial omics and multiplexed imaging to explore cancer biology. Nature Methods, 2021. https://doi.org/10.1038/s41592-021-01203-6 ↩

  3. Xibin Tian, Wuyan Yang, Wei Jiang, et al. Multi-Omics Profiling Identifies Microglial Annexin A2 as a Key Mediator of NF-κB Pro-inflammatory Signaling in Ischemic Reperfusion Injury. Molecular & Cellular Proteomics, 2024. https://doi.org/10.1016/j.mcpro.2024.100723 ↩

  4. Stefania Savoi, Darren C. J. Wong, Asfaw Degu, et al. Multi-Omics and Integrated Network Analyses Reveal New Insights into the Systems Relationships between Metabolites, Structural Genes, and Transcriptional Regulators in Developing Grape Berries (Vitis vinifera L.) Exposed to Water Deficit. Frontiers in Plant Science, 2017. https://doi.org/10.3389/fpls.2017.01124 ↩

  5. Huan Zhong, Yuming Shi, A. S. Kozlova, et al. Omics Insights Into the Effects of Highbush Blueberry and Cranberry Crop Agroecosystems on Honey Bee Health and Physiology. PROTEOMICS, 2025. https://doi.org/10.1002/pmic.70033 ↩