KEGG Pathway Analysis of Your Own RNA-Seq Data

Oak
rna-seq kegg pathway-analysis clusterprofiler deseq2
A small feathered reptile on a steel perch, its translucent plumage revealing a branching network of glowing amber filaments, lit against black.

By the end of this you will have a ranked table of KEGG pathways for a contrast you define (two timepoints, two states, or a trajectory across many draws), a set of annotated pathway maps with your fold changes painted onto them, and a clear sense of which of those results survive scrutiny. You need: FASTQ files from a bulk RNA-seq run (ideally 30M+ paired-end reads per sample, PAXgene or similar whole-blood collection), R 4.3+ with DESeq2, tximport, clusterProfiler, org.Hs.eg.db, pathview, and fgsea, plus salmon or STAR+featureCounts on a machine with 32 GB RAM. Everything below assumes human (hsa in KEGG’s organism codes). KEGG covers thousands of organisms with complete genomes, but the curated map content is deepest for human, mouse, and a handful of model species.

1. Quantify to transcript level, then collapse to gene

KEGG works on genes, not transcripts, but transcript-level quantification with a downstream collapse gives better gene counts than counting reads against gene models directly, because it handles reads that map to multiple isoforms.

# index once, with decoys (full genome as decoy sequence)
salmon index -t gentrome.fa.gz -d decoys.txt -i salmon_idx -k 31 -p 16

salmon quant -i salmon_idx -l A \
  -1 sample_R1.fastq.gz -2 sample_R2.fastq.gz \
  --gcBias --seqBias --posBias \
  --numGibbsSamples 20 \
  -p 16 -o quant/sample01

--gcBias matters more than people expect on library preps separated in time, which is exactly the situation in a longitudinal personal profile: a draw from March and a draw from September went through different reagent lots. --numGibbsSamples 20 gives you inferential replicates so you can see which genes have quantification uncertainty rather than biological signal. If a “differentially expressed” gene turns out to be a two-isoform paralog pair with wide Gibbs spread, drop it before it drags a pathway with it.

Collapse to gene:

library(tximport); library(DESeq2)
tx2gene <- read.csv("tx2gene_gencode_v45.csv")  # tx_id, gene_id (Ensembl)
txi <- tximport(files, type = "salmon", tx2gene = tx2gene,
                countsFromAbundance = "lengthScaledTPM")
dds <- DESeqDataSetFromTximport(txi, colData, design = ~ state)
keep <- rowSums(counts(dds) >= 10) >= min(table(colData$state))
dds <- dds[keep, ]

That filter is the first thing that changes your KEGG results, and nobody reports it. Keeping genes with near-zero counts inflates the universe you test against and shifts every over-representation p-value.

2. Decide what the contrast is

This is where personal RNA-seq breaks from the textbook. You do not have a treatment group and a control group. You have one person sampled repeatedly. Two designs work:

Repeated states. You have four or more draws in state A and four or more in state B (before and after a sustained change in sleep, training load, altitude, season). Model it as ~ state and treat the replicate-to-replicate variance as your error term. Four per state is the floor for DESeq2’s dispersion estimation to be worth anything at n=1 subject.

Trajectory. You have 8-20 draws over a year with no clean grouping. Do not force a two-group contrast. Score every sample against KEGG gene sets directly (step 6) and look at the time series.

If you have whole blood, add one more step before anything else: estimate cell composition. Neutrophil fraction in whole blood swings widely with acute illness, exercise, and cortisol rhythm, and neutrophil-heavy samples will light up KEGG’s neutrophil extracellular trap formation (hsa04613), phagosome (hsa04145), and cytokine-cytokine receptor interaction (hsa04060) maps every time. Run a deconvolution (CIBERSORTx signature matrix LM22, or a marker-gene score), then include the top one or two fractions as covariates:

design(dds) <- ~ neutrophil_frac + state

Without this, the description of most whole-blood KEGG results is “a differential white blood cell count with extra steps.”

3. Map identifiers to Entrez, and check what you lost

KEGG’s human gene IDs are Entrez Gene IDs. clusterProfiler will do the conversion, but silently drops what it cannot map.

library(clusterProfiler); library(org.Hs.eg.db)

res <- lfcShrink(dds, coef = "state_B_vs_A", type = "apeglm")
res <- as.data.frame(res)
res$ensembl <- sub("\\..*$", "", rownames(res))   # strip version suffix

map <- bitr(res$ensembl, fromType = "ENSEMBL", toType = "ENTREZID",
            OrgDb = org.Hs.eg.db)
nrow(map) / nrow(res)   # expect ~0.65-0.75 starting from all of GENCODE

Two things to know. First, the version suffix on GENCODE IDs (ENSG00000141510.17) will fail every lookup if you leave it on, and bitr reports it as a percentage of unmapped input rather than an error. Second, ENSEMBL-to-ENTREZ is many-to-many. Roughly 2-4% of IDs map to multiple Entrez IDs. Deduplicate deliberately, keeping the row with the largest absolute test statistic, rather than letting merge silently duplicate genes into your ranking.

The unmapped 25-35% is mostly lncRNA, pseudogenes, and novel transcripts that KEGG has no annotation for anyway. That is fine, but it means KEGG is testing a protein-coding subset of your data. Note the number and move on.

4. Run GSEA over KEGG, not over-representation

Over-representation analysis (ORA) takes your significant gene list, intersects with each pathway, and runs a hypergeometric test. It throws away effect size and it depends entirely on where you set the significance cutoff. With n=1 longitudinal data your cutoff is arbitrary and your list is short. We use ranked GSEA as the primary method and keep ORA as a cross-check.

Rank by the Wald statistic from DESeq2, not by fold change. Fold change alone puts noisy low-expression genes at the extremes. The Wald statistic already incorporates uncertainty.

library(fgsea)

rk <- res[!is.na(res$stat), ]
rk <- merge(rk, map, by.x = "ensembl", by.y = "ENSEMBL")
rk <- rk[order(-abs(rk$stat)), ]
rk <- rk[!duplicated(rk$ENTREZID), ]
ranks <- setNames(rk$stat, rk$ENTREZID)
ranks <- sort(ranks, decreasing = TRUE)

kegg <- clusterProfiler::download_KEGG("hsa")
pathways <- split(kegg$KEGGPATHID2EXTID$to, kegg$KEGGPATHID2EXTID$from)

set.seed(42)
fg <- fgseaMultilevel(pathways, ranks, minSize = 12, maxSize = 400, eps = 0)
fg <- fg[order(fg$padj), ]

minSize = 12 drops the tiny maps where three genes drive everything. maxSize = 400 drops “Metabolic pathways” (hsa01100), which contains over a thousand genes and is never interpretable. eps = 0 lets fgseaMultilevel estimate p-values below the permutation floor instead of clamping them at 1e-10, which matters when you are ranking dozens of pathways by significance.

Note the KEGG release you pulled, because download_KEGG hits the live REST API and your results are not reproducible without it:

KEGGREST::keggInfo("kegg")   # prints release number and date

Then the ORA cross-check, with an explicit universe:

sig <- rk$ENTREZID[rk$padj < 0.05 & !is.na(rk$padj)]
ora <- enrichKEGG(gene = sig, organism = "hsa",
                  universe = rk$ENTREZID,      # every tested gene, not all of KEGG
                  pvalueCutoff = 0.1, qvalueCutoff = 0.2)

The universe argument is the single most common mistake in published KEGG analyses. Leave it out and clusterProfiler uses all annotated genes, including the thousands you filtered for low expression, which makes almost everything look enriched. Pathways that appear in both the GSEA and ORA output are worth your attention. Pathways that appear only in ORA usually reflect the cutoff.

5. Read the maps instead of the table

A KEGG pathway map is a hand-drawn diagram of a biochemical process: boxes are gene products (specifically KEGG Orthology groups, so one box often covers a whole gene family), arrows are reactions or regulatory relations, circles are small-molecule compounds. The map has a defined topology, unlike a GO term, which is only a list. That topology is the reason to use KEGG at all.

Paint your data on it:

library(pathview)
fc <- setNames(rk$log2FoldChange, rk$ENTREZID)
pathview(gene.data = fc, pathway.id = "hsa04151",   # PI3K-Akt
         species = "hsa", limit = list(gene = 2),
         low = "blue", mid = "grey", high = "red",
         kegg.native = TRUE, out.suffix = "stateB_vs_A")

Now look at what colored. Three patterns to distinguish:

  • Coherent directional shift along a branch of the map (a receptor, its adaptor, and two downstream effectors all up). This is the interesting case.
  • Scattered coloring with no topological relationship. Usually the pathway contains a housekeeping module or a set of large gene families that appear on 30 different maps.
  • Coloring concentrated in genes that are also annotated to five other significant pathways. KEGG maps overlap heavily. Check with setdiff on the leading-edge gene sets before treating two hits as two findings.

The leading edge is the part of fgsea output that carries information:

fg$leadingEdge[[1]]
mapIds(org.Hs.eg.db, fg$leadingEdge[[1]], "SYMBOL", "ENTREZID")

If a pathway’s leading edge is six genes out of 180, the enrichment is a statement about those six genes, not about the pathway.

6. For trajectories, score each sample

With no clean two-group split, per-sample pathway scoring is a better readout than any contrast. GSVA computes a pathway-level enrichment score per sample directly from the expression matrix.

library(GSVA)
vsd <- assay(vst(dds, blind = TRUE))
rownames(vsd) <- rk$ENTREZID[match(sub("\\..*$", "", rownames(vsd)), rk$ensembl)]
vsd <- vsd[!is.na(rownames(vsd)), ]

par <- gsvaParam(vsd, pathways, minSize = 12, maxSize = 400, kcdf = "Gaussian")
scores <- gsva(par)    # pathways x samples

Now you have a time series per pathway. Plot them against draw date. The thing to establish first is your noise floor: if you have a technical replicate (same blood draw, two library preps), the spread between them tells you how large a score change has to be before it means anything. In our experience that floor is larger than most people assume, and a lot of apparent seasonal pathway drift sits under it.

GSVA scores are relative to the sample set, so adding new draws changes old scores. Re-score the whole matrix every time, and never compare scores across two runs with different sample sets.

7. Cross-check against another molecular layer

A KEGG hit from RNA alone is a hypothesis about one layer. Transcript changes do not reliably propagate to protein, and KEGG maps mix enzymes, receptors, and metabolites that you can measure independently. The strongest personal-data workflow is to check whether the same map lights up from a second assay: does the proteomics panel show the same direction for the leading-edge genes that have protein coverage, do the relevant metabolites or clinical biomarkers move, does CGM-derived glycemic variability track a hit on insulin signaling (hsa04910).

This is how the large multi-omics studies reach conclusions worth acting on. The NASA twin and rodent spaceflight analysis converged on mitochondrial stress as a shared hub precisely because the signal appeared across transcriptomic, proteomic, and metabolomic layers and across tissues, rather than in one assay.1 Integrated mitochondrial analyses in osteoarthritis followed the same structure, combining expression data with additional layers before treating a pathway as a finding.2 Tooling exists for the overlay itself: metaKEGG renders multiple omics layers onto a single KEGG map so you can see which nodes are supported by more than one measurement.3 Pathway-level frameworks that join genotype to expression take a similar view, using pathway membership as the unit that links variants to transcriptional signal.4

Nothing in this workflow is diagnostic. KEGG enrichment describes coordinated movement in a gene set relative to your own earlier samples. It does not establish disease, and pathway names that sound clinical (hsa05010, hsa04930) are named after the research context in which the map was curated, not after your state. If something here concerns you, take the underlying biomarker values, not the pathway plot, to a physician.

Common problems

Everything is significant. Almost always a missing universe in enrichKEGG, or no low-count filter. Second most common cause: cell composition in whole blood. Re-run with the neutrophil fraction in the design.

Nothing is significant, but the volcano plot looks alive. Check your ID mapping rate. If bitr mapped 40% of your genes, you probably left version suffixes on or passed symbols where Ensembl IDs were expected. Also check that your ranking vector is named with Entrez IDs, not symbols, since fgsea will silently return an empty result rather than erroring on a total ID mismatch.

Results changed between last month and today. KEGG is versioned and updated continuously. download_KEGG fetches current content. Cache the gene sets to disk with a release stamp and reuse that file for any comparison across time.

Globin dominates the library. Whole blood without globin depletion routinely gives half or more of reads to HBB, HBA1, HBA2. Your effective depth collapses and low-expression genes vanish. Either depletion at prep time or globin read removal before quantification, plus honest accounting of the remaining depth.

A pathway hit is driven by one gene family. Look at leading-edge symbols. If they are all HLA-, all COL, or all RPL*/RPS*, you found a module, and it will appear on many maps.

Directionality is wrong on the map. Pathview colors by whatever vector you pass. If you passed unshrunken log2 fold changes, low-count genes will show extreme colors. Pass apeglm-shrunken values, and set limit = list(gene = 2) so the color scale is not blown out by one outlier.

You cannot license the data. KEGG is free to browse and query over the REST API for academic use. Bulk FTP download and commercial use require a paid license from Pathway Solutions. If you are building anything productized, read the terms before you cache the whole database.

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. Willian A. da Silveira, Hossein Fazelinia, Sara Brin Rosenthal, et al. Comprehensive Multi-omics Analysis Reveals Mitochondrial Stress as a Central Biological Hub for Spaceflight Impact. Cell, 2020. https://doi.org/10.1016/j.cell.2020.11.002 ↩

  2. Yinteng Wu, Haifeng Hu, Tao Wang, et al. Characterizing mitochondrial features in osteoarthritis through integrative multi-omics and machine learning analysis. Frontiers in Immunology, 2024. https://doi.org/10.3389/fimmu.2024.1414301 ↩

  3. Michail Lazaratos, Neele Haacke, Jasmin Gaugel, et al. metaKEGG: A comprehensive algorithm package to visualize multi-omics pathway enrichment. Biochemistry and Biophysics Reports, 2026. https://doi.org/10.1016/j.bbrep.2026.102723 ↩

  4. Yuvraj Singh, Nikhil Sharma, Kabir Raghuvanshi, et al. CardioWAS Pathway-Based RNA-Seq and Genome Integration Platform for CVD Omics-Wide Associations. bioRxiv (Cold Spring Harbor Laboratory), 2025. https://doi.org/10.1101/2025.10.06.679779 ↩