How to Run GO Enrichment on Your Own RNA-Seq Data

Oak
rna-seq gene-ontology gsea clusterprofiler transcriptomics
A lab machine lifts glowing strands from a sample tray into a branching tree of light inside its glass column.

By the end of this you will have two complementary outputs from one RNA-seq experiment: a table of Gene Ontology (GO) terms overrepresented among your differentially expressed genes, with Benjamini-Hochberg adjusted p-values and a bias-corrected version of the same test, and a gene set enrichment result computed on the full ranked gene list so that you are not dependent on a significance cutoff. You will also have a redundancy-collapsed summary, because raw GO output is dominated by parent and child terms that say the same thing. What you need going in: a gene-level count matrix (genes by samples, integers, from featureCounts, HTSeq, or the salmon/tximport route), sample metadata with at least one contrast of interest, R 4.4 or later, and the packages DESeq2, clusterProfiler, org.Hs.eg.db, GO.db, goseq, fgsea, msigdbr, and rrvgo. Two biological replicates per condition is the practical floor for differential expression; with serial samples from one person, the “replicates” are timepoints and the design needs to reflect that, which we cover in step 2.

1. What GO tests, and how it differs from KEGG and GSEA

Gene Ontology is a structured vocabulary of gene function organized as a directed acyclic graph. Each node is a GO term with a stable identifier (for example GO:0006954, inflammatory response) and terms are connected by relations such as is_a and part_of. The graph is split into three aspects: biological process (BP), the broader program a gene contributes to; molecular function (MF), the biochemical activity of the gene product; and cellular component (CC), where it acts. Annotations propagate up the graph, so a gene annotated to a specific child term is implicitly annotated to all of its ancestors. This is the single most important structural fact about GO, and it is why your results will contain nested, highly correlated terms.

GO enrichment answers a narrow statistical question: among the genes in my list, are genes annotated to term T more frequent than you would expect by chance given the background set of genes I could have detected? The standard test is a one-sided Fisher exact test on a two-by-two table, repeated across thousands of terms, with a multiple-testing correction. KEGG is a different resource, a curated collection of pathway maps with directed molecular interactions and metabolic reactions, so KEGG enrichment tests membership in defined pathways rather than position in a functional hierarchy. KEGG terms are fewer, less redundant, and more mechanistically specific, while GO covers far more genes and far more of the transcriptome. Studies that need functional coverage usually run both and read them together, treating agreement between the two as the meaningful signal rather than relying on either alone 1. Multi-omics work routinely uses exactly this pairing to interpret differentially expressed gene sets 2.

Gene set enrichment analysis (GSEA) differs along a different axis. Overrepresentation tests require you to split genes into “significant” and “not,” which discards the magnitude of every effect and makes the result sensitive to your threshold. GSEA takes the entire gene list ranked by some statistic and asks whether the members of a gene set are distributed non-randomly toward either end of the ranking, using a running-sum enrichment score with a permutation-based null. GO and GSEA are not alternatives in the sense people often assume: GSEA is a method, GO is a source of gene sets, and gseGO runs GSEA over GO terms. The real choice is between thresholded and threshold-free testing, and we run both.

2. Build a gene list you are willing to defend

Everything downstream inherits the quality of this step. Run DESeq2 on raw counts, not on TPM or FPKM, and give it the actual design.

library(DESeq2)
dds <- DESeqDataSetFromMatrix(counts, coldata, design = ~ condition)
keep <- rowSums(counts(dds) >= 10) >= min(table(coldata$condition))
dds <- dds[keep, ]
dds <- DESeq(dds)
res <- lfcShrink(dds, coef = "condition_B_vs_A", type = "apeglm")

Three decisions matter. First, prefilter by requiring a count of at least 10 in as many samples as the smallest group, which removes genes with no power and shrinks the multiple-testing burden without cherry-picking. Second, use shrunken log fold changes (apeglm) for ranking and for any fold-change filter, because unshrunken estimates for low-count genes are enormous and will drag irrelevant terms into your results. Third, if you have repeated samples from the same person over time, include the subject or batch as a blocking factor, or use ~ subject + condition for a paired within-person contrast. A within-person design is the strongest one available in an n-of-1 setting because it removes the genotype and most of the stable regulatory variation from the comparison.

For the thresholded test, take padj < 0.05 and, if the gene list is very large, add abs(log2FoldChange) > 0.5. Separate the up- and down-regulated genes and test them independently. Merging directions is common and it is a mistake: a term with half its genes up and half down will be called enriched while the underlying biology is mixed, and you will have no idea which way to read it.

The background universe is the other half of the test and it is where most published GO analyses go wrong. The correct universe is the set of genes that could have entered your DE list, meaning the genes surviving your prefilter, not all ~20,000 protein-coding genes and certainly not the full annotation. Using an inflated universe systematically inflates significance for any term enriched in expressed genes generally, which in most tissues means translation, mitochondrion, and RNA processing appear in every analysis you ever run.

3. Map identifiers cleanly

GO annotation packages are keyed on Entrez gene IDs and symbols, while your counts are almost certainly keyed on Ensembl gene IDs, often with version suffixes. Strip versions before mapping and keep the mapping explicit so you can audit losses.

library(org.Hs.eg.db)
library(AnnotationDbi)

strip <- function(x) sub("\\..*$", "", x)
de_ens  <- strip(rownames(res)[which(res$padj < 0.05 & res$log2FoldChange > 0)])
uni_ens <- strip(rownames(res))

map <- AnnotationDbi::select(org.Hs.eg.db, keys = uni_ens,
                             keytype = "ENSEMBL",
                             columns = c("ENTREZID", "SYMBOL"))
map <- map[!is.na(map$ENTREZID), ]
map <- map[!duplicated(map$ENSEMBL), ]   # keep one Entrez per Ensembl

Expect to lose 20 to 40 percent of rows, mostly non-coding genes, pseudogenes, and retired identifiers. That loss is not neutral. Long non-coding RNAs are nearly unannotated in GO, so a real lncRNA-driven signal will vanish silently. Check how many of your DE genes survived mapping and report that number alongside the enrichment table.

4. Run the overrepresentation test

clusterProfiler::enrichGO handles the Fisher test, the propagation through the GO graph, and the correction.

library(clusterProfiler)

ego_up <- enrichGO(gene          = map$ENTREZID[map$ENSEMBL %in% de_ens],
                   universe      = map$ENTREZID,
                   OrgDb         = org.Hs.eg.db,
                   keyType       = "ENTREZID",
                   ont           = "BP",
                   pAdjustMethod = "BH",
                   pvalueCutoff  = 0.05,
                   qvalueCutoff  = 0.10,
                   minGSSize     = 10,
                   maxGSSize     = 500,
                   readable      = TRUE)

Set maxGSSize = 500. Terms with thousands of members (GO:0008152, metabolic process, has tens of thousands) are statistically significant in almost any large gene list and carry no information. minGSSize = 10 removes terms where one or two genes produce an enormous apparent fold enrichment. Read the output by FoldEnrichment and Count together, not by p-value alone: a term with 4 of 11 genes and padj = 0.03 is a much weaker claim than one with 60 of 300 genes at the same padj.

5. Correct for length and expression bias

RNA-seq detects differential expression in long, highly expressed genes more readily than in short ones, because read count scales with transcript length and count-based power scales with count. GO terms are not length-neutral: extracellular matrix, cell adhesion, and synaptic terms are enriched for very long genes, and ribosomal and mitochondrial terms for very short ones. A plain Fisher test cannot tell length-driven detection bias from biology. The goseq package fits a probability weighting function of DE likelihood against gene length and uses a Wallenius non-central hypergeometric approximation in place of the hypergeometric.

library(goseq)

all_genes <- as.integer(uni_ens %in% de_ens)
names(all_genes) <- uni_ens

pwf <- nullp(all_genes, genome = "hg38", id = "ensGene", plot.fit = TRUE)
go_wall <- goseq(pwf, genome = "hg38", id = "ensGene", test.cats = "GO:BP")
go_wall$padj <- p.adjust(go_wall$over_represented_pvalue, method = "BH")

Look at the nullp plot before anything else. If the fitted curve is close to flat, length bias is negligible in your data and the corrected and uncorrected results will agree. If it rises steeply, treat any term that is significant in enrichGO but not in goseq as suspect. We use the intersection of the two as the reportable set and mention the difference explicitly. The same weighting logic applies to expression level rather than length if you supply mean normalized counts as the bias vector, which is often the more relevant covariate when you are using transcript-level quantification that already accounts for length.

6. Run GSEA on the full ranked list

Now drop the threshold entirely. Rank every mapped gene by the Wald statistic from DESeq2, which incorporates both effect size and its uncertainty and is better behaved than ranking on log fold change or on a p-value.

ranks <- res$stat
names(ranks) <- strip(rownames(res))
ranks <- ranks[!is.na(ranks)]
ranks <- ranks[names(ranks) %in% map$ENSEMBL]
names(ranks) <- map$ENTREZID[match(names(ranks), map$ENSEMBL)]
ranks <- sort(ranks, decreasing = TRUE)

gse <- gseGO(geneList = ranks, OrgDb = org.Hs.eg.db, keyType = "ENTREZID",
             ont = "BP", minGSSize = 15, maxGSSize = 500,
             pvalueCutoff = 0.05, pAdjustMethod = "BH",
             eps = 0, nPermSimple = 10000, seed = TRUE)

Set eps = 0 so that very small p-values are not floored, and raise nPermSimple to 10,000 if you see the multilevel-estimation warning. Ties in the ranking break GSEA’s assumptions, so if many genes share a statistic of exactly zero, remove them rather than leaving a block of ties in the middle of the list. The sign of NES (normalized enrichment score) gives direction, and core_enrichment lists the genes driving the leading edge, which is where the interpretable content lives. For a signal that appears in GSEA but not in the thresholded test, the usual explanation is a coordinated but individually modest shift across many genes in a set, which is precisely the case GSEA exists to catch.

If you also want curated pathway sets, pull Hallmark and Reactome collections through msigdbr and run fgseaMultilevel on the same ranking. Hallmark’s 50 sets are the most readable starting point in our experience because they were explicitly built to reduce redundancy.

7. Collapse redundancy before you read the results

A typical BP result contains 200 significant terms describing perhaps eight distinct processes. Use clusterProfiler::simplify, which drops terms whose Wang semantic similarity to a more significant term exceeds a cutoff.

ego_s <- simplify(ego_up, cutoff = 0.7, by = "p.adjust",
                  select_fun = min, measure = "Wang")

For a coarser view, rrvgo computes a semantic similarity matrix, clusters it, and gives you a treemap where block size reflects term size and block grouping reflects semantic relatedness. Report the collapsed set and keep the full table as a supplement. Do not hand-pick the three terms that match your hypothesis out of a list of 200, since with a redundant hierarchy you can find support for almost any narrative that way.

8. Interpret with the design in mind

GO enrichment describes the gene list you gave it. It does not establish causation, and it does not transfer a group-level finding to an individual. If your RNA-seq is whole blood or PBMC, the dominant driver of transcriptome variation is usually cell type composition, not per-cell regulation, so a result loaded with neutrophil degranulation or lymphocyte activation terms may reflect a shift in white cell proportions. Estimate composition with a deconvolution method (CIBERSORTx, xCell) and include the estimated fractions as covariates in the DESeq2 design, then rerun the enrichment and see what survives.

The strongest version of this analysis at the individual level is longitudinal. With serial samples, you can contrast timepoints within yourself, use each earlier timepoint as its own control, and require that a GO term move consistently across multiple independent timepoints before you take it seriously. Pairing transcriptome changes with proteomic and metabolomic readouts on the same samples raises confidence further, since concordance across layers is much harder to produce by chance than enrichment in one layer 3. Integrated designs of this kind are how multi-omics studies localize a process to a specific tissue and pathway rather than stopping at a list of terms 4.

None of this output is a diagnosis. A GO term is a label on a statistical pattern in one sample set, and translating any of it into a statement about your health, or into any action involving a medication, supplement, or dose, requires a physician who can see your clinical picture. Bring the tables and the design description, not the interpretation.

9. Make the run reproducible

GO is revised continuously and annotations are added and retracted, so the same gene list will give different results six months later. Record the versions.

sessionInfo()
packageVersion("org.Hs.eg.db")
GO.db::GO_dbInfo()   # includes the GO release date

Save the exact gene list, the universe, and the full unfiltered result table to disk alongside the script. When you compare a new timepoint to an old one, rerun both through the same annotation snapshot rather than comparing a fresh result to an archived one.

Common problems

Everything is significant and the top terms are all translation, ribosome, and mitochondrion. Your universe is wrong, almost certainly the whole genome instead of the expressed genes. Set universe to your prefiltered, successfully mapped gene set.

Nothing is significant. Check the size of your DE list first. Below roughly 50 genes, overrepresentation testing has little power, and GSEA on the ranked list is the better tool. Also check whether your contrast is confounded by batch, which inflates dispersion and suppresses everything.

The GO result contradicts the KEGG result. This is common and usually not an error. GO covers more genes with looser functional definitions while KEGG covers fewer genes in defined pathways, and the two can disagree in emphasis on the same gene list. Look at the gene-level overlap directly rather than trying to reconcile term names 1.

Long-gene terms dominate. Compare enrichGO with goseq and inspect the nullp fit. If the bias curve is steep, report only terms significant under the Wallenius correction.

Results change when you switch the fold-change cutoff. That instability is the main weakness of thresholded testing. Prefer GSEA for the headline result and use the thresholded test as a confirmation.

Duplicate or missing identifiers throw errors in gseGO. The ranking vector must have unique names, be numeric, and be sorted in decreasing order. Collapse duplicated Entrez IDs by keeping the one with the largest absolute statistic before sorting.

An enriched term looks biologically exciting but rests on four genes. Check Count and FoldEnrichment and look at the individual genes’ fold changes and adjusted p-values in the DESeq2 output. Single-gene-driven terms are the most common source of overinterpretation in published GO analyses, and they are equally common in personal data 5.

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. Yu-Hang Zhang, Tao Zeng, Lei Chen, et al. Determining protein–protein functional associations by functional rules based on gene ontology and KEGG pathway. Biochimica et Biophysica Acta (BBA) - Proteins and Proteomics, 2021. https://doi.org/10.1016/j.bbapap.2021.140621 ↩ ↩2

  2. Huiyan Zhao, Guoxia Shang, Nengwen Yin, et al. Multi-omics analysis reveals the mechanism of seed coat color formation in Brassica rapa L. Theoretical and Applied Genetics, 2022. https://doi.org/10.1007/s00122-022-04099-8 ↩

  3. Yun Fan, Shiqi Li, Xiancheng Yang, et al. Multi-omics approach characterizes the role of Bisphenol F in disrupting hepatic lipid metabolism. Environment International, 2024. https://doi.org/10.1016/j.envint.2024.108690 ↩

  4. Eva Bang Harvald, Richard R. Sprenger, Kathrine Brændgaard Dall, et al. Multi-omics Analyses of Starvation Responses Reveal a Central Role for Lipoprotein Metabolism in Acute Starvation Survival in C. elegans. Cell Systems, 2017. https://doi.org/10.1016/j.cels.2017.06.004 ↩

  5. Zekai Nian, Yicheng Mao, Zexia Xu, et al. Multi-omics analysis uncovered systemic lupus erythematosus and COVID-19 crosstalk. Molecular Medicine, 2024. https://doi.org/10.1186/s10020-024-00851-6 ↩