# How to Perform Genomic Annotation Enrichment Analysis in R By the end of this guide you will be able to take a list of genes, for example candidate genes from a GWAS, a set of differentially expressed genes, or genes near QTL intervals, and test whether their functional annotations are overrepresented or collectively shifted, using R. You will run over-representation analysis (ORA) and gene set enrichment analysis (GSEA) against Gene Ontology (GO) and KEGG terms, visualize the results, export them, and adapt the workflow to a non-model organism. Before you start, have the following on hand: - R 4.3 or newer and RStudio, or any environment where you can install Bioconductor packages - The `clusterProfiler`, `enrichplot`, `BiocManager`, and `AnnotationHub` packages - An organism annotation resource, either an `OrgDb` package such as `org.Hs.eg.db` or a custom annotation you built yourself - A gene list with a known identifier type (SYMBOL, ENTREZID, ENSEMBL, or similar) - A decision about your background universe, which I explain in the walkthrough Genomic annotation enrichment analysis is one of the most frequently requested tasks in R for bioinformatics, and it is also one of the most frequently done wrong. The steps below follow the conventions used in published studies, from soybean meta-QTL candidate genes [1] to shrimp functional annotation platforms [2] and fungal secretome analyses [3], so the workflow you build here matches what reviewers expect to see. ## What Genomic Annotation Enrichment Analysis Does Functional annotation assigns biological meaning to genes or genomic regions. Gene Ontology describes genes across three ontologies: biological process (BP), molecular function (MF), and cellular component (CC). KEGG maps genes to reference pathways. Other common vocabularies include Reactome, Disease Ontology, and InterPro protein domains. Enrichment analysis asks whether some annotation category appears more often in your gene list than expected by chance. There are two main formulations. **Over-representation analysis (ORA)** takes a discrete list of genes, counts how many fall into each term, and tests the count against a hypergeometric or Fisher's exact null model built on a background universe. This is the approach behind GO and KEGG over-representation analyses of predicted secretomes [3] and behind the GO term testing that identified response to salicylic acid as the only significant biological process after FDR correction in a soybean meta-QTL study [1]. **Gene set enrichment analysis (GSEA)** takes all genes ranked by a continuous statistic, such as log2 fold change or a GWAS signal, and asks whether genes belonging to a term cluster toward one end of the ranking. ChromTag, an R Shiny platform for CUT&Tag and ChIP-seq data, supports both ORA-based GO and KEGG analysis and preranked GSEA using gene-level summaries derived from peak annotation [4]. Both approaches produce a p-value per term, which you then correct for multiple testing. The standard correction in `clusterProfiler` is the Benjamini-Hochberg false discovery rate (FDR) procedure, which controls the expected proportion of false discoveries among your called-significant terms. A term with `p.adjust = 0.03` is expected to contribute false positives at a 3 percent rate across your entire significant set. ## Choosing the Right Analysis for Your Input ![Decision path from input data type to ORA, preranked GSEA, comparative ORA or region-based chromatin enrichment.](/article-images/illustrations/genomic-annotation-enrichment-analysis-r-method-choice.webp) Your input data type determines the method. Most enrichment mistakes begin here, with an analysis format that does not match the data. | Input data type | Appropriate method | Typical R function | |--|--|--| | Discrete gene list (DE genes, candidate genes under QTL) | ORA | `enrichGO()`, `enrichKEGG()` | | All genes ranked by fold change or statistic | GSEA (preranked) | `gseGO()`, `gseKEGG()` | | Multiple gene lists to compare | Comparative ORA | `compareCluster()` | | ChIP-seq or CUT&Tag peaks assigned to nearby genes | ORA plus preranked GSEA on gene summaries | `enrichGO()`, `GSEA()` [4] | | GWAS lead variants | Variant annotation first, then gene-list ORA | `SNPannotator`, then `enrichGO()` [5] | | Predefined genomic regions with epigenomic signal | Region-based chromatin enrichment | `ChromCall` package [6] | Two recent studies illustrate the pattern. In a maize drought GWAS, candidate genes identified under 45 associated loci were functionally annotated and then interpreted for drought-response biology [7]. In a yak GWAS for average daily gain, candidate genes from significant loci were subjected to functional enrichment, which highlighted the cAMP signaling pathway and extracellular matrix remodeling terms [8]. In both cases, the enrichment step ran on a candidate gene list derived from association signals, not on the full SNP table. That distinction matters, and it leads directly to the question of the background universe. ## Step-by-Step Walkthrough: Genomic Annotation Enrichment Analysis in R This walkthrough uses a human example with `org.Hs.eg.db`. The same logic applies to any organism with an `OrgDb`, and the non-model organism section covers organisms without one. ### Step 1: Install and Load the Packages ```r install.packages("BiocManager") BiocManager::install(c( "clusterProfiler", "enrichplot", "org.Hs.eg.db", "AnnotationHub", "ReactomePA" )) ``` Then load them: ```r library(clusterProfiler) library(enrichplot) library(org.Hs.eg.db) ``` `clusterProfiler` is the workhorse package for this task in R and is also the engine behind web platforms built for specific organism groups, including the ShEnrich platform for penaeid shrimp, which constructed GMT-format libraries for KEGG and GO enrichment and implemented them with `clusterProfiler` [2]. If a paper reports GO or KEGG enrichment in R, there is a good chance this package produced it. ### Step 2: Prepare Your Gene List For this example, imagine you ran a GWAS or a differential expression study and extracted 20 candidate genes by symbol. In practice, a GWAS-derived list would come from genes near significant loci, exactly as in the maize and yak studies cited above [7], [8]. ```r gene_symbols <- c( "TP53", "BRCA1", "PTEN", "EGFR", "KRAS", "MYC", "CDKN2A", "RB1", "ATM", "CHEK2", "ATR", "BRCA2", "MLH1", "MSH2", "APC", "PIK3CA", "AKT1", "MTOR", "BRAF", "NRAS" ) ``` Record the source of every gene before this step. A gene list with an unclear provenance cannot be defended later, and the background universe question in the next step depends on knowing where the list came from. ### Step 3: Define the Background Universe The universe is the set of genes that had a chance to appear in your list. Choices include: - All genes in the genome annotation - All genes detected in your RNA-seq experiment (genes with nonzero counts or above an expression filter) - All genes represented on your SNP array or capture probe set - All genes in the genomic regions you searched If your candidate genes came from genes near significant SNPs on an array, the correct universe is all genes near SNPs on that array, not the whole genome. Using the whole genome as background when your detectable set was much smaller inflates apparent enrichment. The variant annotation workflow in the insulin and osteoarthritis study mapped GWAS variants through Ensembl VEP and other annotation tools before enrichment, which implicitly defines the universe as genes mappable from the variant set [9]. For the toy example, define the universe as all protein-coding genes in the `OrgDb`: ```r universe_keys <- keys(org.Hs.eg.db, keytype = "ENTREZID") ``` ### Step 4: Map Your Identifiers `clusterProfiler` works best with Entrez Gene IDs for KEGG and with Entrez or Ensembl for GO. Convert symbols first and check the mapping rate: ```r gene_map <- bitr( gene_symbols, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db ) ``` Expected output is a data frame with columns `SYMBOL` and `ENTREZID`, plus console messages such as `0.95% of input gene IDs are fail to map...`. That message is informative, not an error. If a large fraction fails to map, for example more than 20 percent, stop and check your identifier type. A frequent cause is passing Ensembl IDs with version suffixes such as `ENSG00000141510.16`. Strip the version first: ```r gene_symbols_clean <- sub("\\..*$", "", gene_symbols) ``` Do the same for the universe so both sets share an identifier space: ```r universe_map <- bitr( universe_keys, fromType = "ENTREZID", toType = "ENTREZID", OrgDb = org.Hs.eg.db ) universe <- universe_map$ENTREZID gene_list <- gene_map$ENTREZID ``` A common runtime error at this stage is `Error in bitr(): No gene can be mapped`. This means the identifier type does not match the `OrgDb`. Print the first few entries of your input and compare them against `head(keys(org.Hs.eg.db, keytype = "SYMBOL"))` to confirm the format. ### Step 5: Run GO Over-Representation Analysis ```r ego <- enrichGO( gene = gene_list, universe = universe, OrgDb = org.Hs.eg.db, keyType = "ENTREZID", ont = "BP", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.05, readable = TRUE ) ``` Key parameters: - `ont`: the ontology, one of `BP`, `MF`, `CC`, or `ALL` - `pAdjustMethod`: `BH` (Benjamini-Hochberg FDR) is the default and the field standard. Use `bonferroni` only if you want maximum conservatism. - `pvalueCutoff` and `qvalueCutoff`: thresholds applied to raw p and adjusted q values. Set them to 1 during exploration and filter afterward, because cutoffs applied inside the function also trim the output object. - `universe`: the background defined in Step 3. Omitting it silently defaults to all genes in the `OrgDb`, which is only correct if that truly was your detectable set. - `readable`: converts output Entrez IDs back to symbols for readability. Inspect the results: ```r head(as.data.frame(ego)) ``` Expected output structure (your values will differ, this shows the columns): | ONTOLOGY | ID | Description | GeneRatio | BgRatio | pvalue | p.adjust | qvalue | geneID | |--|--|--|--|--|--|--|--|--| | BP | GO:0006915 | apoptotic process | 12/19 | 1452/19556 | small | small | small | TP53, BRCA1, ... | `GeneRatio` is the fraction of your genes in the term. `BgRatio` is the fraction of the universe in the term. If `BgRatio` looks wrong, your universe is wrong. ### Step 6: Run KEGG Over-Representation Analysis ```r ekk <- enrichKEGG( gene = gene_list, universe = universe, organism = "hsa", keyType = "ncbi-geneid", pAdjustMethod = "BH", pvalueCutoff = 0.05 ) ekk <- setReadable(ekk, OrgDb = org.Hs.eg.db, keyType = "ENTREZID") ``` The `organism` argument takes a three-letter KEGG code such as `hsa` for human, `bta` for cattle, or `osa` for rice. If KEGG lookups fail with connection errors, the package fetches the reference pathways online, so a stable internet connection is required during the run. ### Step 7: Run Preranked GSEA When You Have a Ranking If you have fold changes or another continuous statistic for all measured genes, GSEA uses more information than ORA because it does not require an arbitrary significance cutoff. ```r log2fc <- rnorm(length(universe), 0, 1) names(log2fc) <- universe ranked <- sort(log2fc, decreasing = TRUE) eg <- gseGO( geneList = ranked, ont = "BP", OrgDb = org.Hs.eg.db, minGSSize = 10, maxGSSize = 500, pAdjustMethod = "BH", eps = 1e-50 ) ``` Notes on the parameters: - `geneList` must be a named, sorted numeric vector, highest value first. Sorting is not optional. An unsorted vector produces meaningless enrichment scores. - `minGSSize` and `maxGSSize` exclude terms that are too small to test reliably or too broad to interpret. Defaults of 10 and 500 are reasonable starting points. - The ranking metric choice matters. Log2 fold change alone ignores significance, so many analysts use the Wald statistic or a signed combination such as `-log10(p) * sign(log2fc)`. ChromTag implements the same logic for chromatin data, separating upregulated and downregulated regions before running GSEA so that activated and repressed regulatory programs are not mixed into one signature [4]. ### Step 8: Visualize the Results ```r dotplot(ego, showCategory = 15, split = "ONTOLOGY") + theme(axis.text.y = element_text(size = 9)) emapplot(ego, showCategory = 20) cnetplot(ego, showCategory = 5, foldChange = log2fc) ``` - `dotplot()` shows GeneRatio on the x axis, term on the y axis, and a color or size channel for gene count or p.adjust. This is the plot most journals expect. - `emapplot()` draws terms as nodes connected by shared genes, which reveals redundant terms clustering together. - `cnetplot()` shows which genes drive each term, which is the fastest way to sanity-check an enrichment result. For multi-condition comparisons, `compareCluster()` runs ORA per group in one call and produces a single dot plot with condition on a facet, which is the cleanest way to present contrasts. ### Step 9: Export Results ```r write.csv(as.data.frame(ego), "go_BP_enrichment.csv", row.names = FALSE) write.csv(as.data.frame(eg), "gsea_BP_results.csv", row.names = FALSE) ``` Keep the full unfiltered table in your exported file. Reviewers and collaborators frequently want to see terms just above your cutoff, and a filtered file destroys that option. ## Working With Non-Model Organisms Model organism `OrgDb` packages cover only a fraction of species. Everything else requires a custom annotation, and the R ecosystem supports two practical routes. The first route is building a GMT file (a text format with one term per line, the term name, then the gene members) from your own annotation pipeline and feeding it to `clusterProfiler` through `enricher()` and a TERM2GENE mapping. ShEnrich took exactly this approach for five penaeid shrimp species, building annotation from BLASTP against NCBI nr, InterProScan domain prediction, and eggNOG-mapper orthology, then constructing GMT libraries for KEGG and GO that `clusterProfiler` consumes [2]. Similar multi-tool annotation pipelines, combining eggNOG-mapper, GO, KEGG, and specialized databases such as PHI-base for pathogenicity, are standard practice in fungal genomics [3], and pan-genome studies use the same annotation stack to characterize core and accessory gene functions [10]. The second route is packaging the annotation as an installable R library. A published protocol walks through building custom annotation R packages from NCBI genomic data that import directly into R and work with `clusterProfiler`, demonstrated with Mycobacterium tuberculosis H37Rv [11]. This route keeps your annotation versioned and reproducible, which matters when collaborators rerun the analysis months later. The central caveat for non-model organisms is annotation transfer. Homology-based transfer from a close relative can introduce errors and misses species-specific functions, which the custom annotation protocol explicitly flags as the motivation for building species-specific libraries [11]. When you cite an enrichment result from a non-model organism, state where the annotation came from and how recent it is. A GO term assigned by a 2018 BLAST transfer and one assigned by current eggNOG-mapper orthology are not equivalent evidence. ## Beyond Gene Lists: Variants and Genomic Regions Two extensions of the basic workflow deserve their own mention because they are increasingly common in post-GWAS pipelines. For variant-centered work, annotate the variants first and derive the gene list afterward. SNPannotator is a CRAN package written in R that automates this for GWAS lead variants, finding LD proxies, querying Ensembl, GTEx, the eQTL Catalog, and STRING DB, and compiling functional enrichment results into reports [5]. Integrative studies commonly pair this kind of variant annotation with exploratory GO biological process and KEGG enrichment on the mapped genes [9]. For region-centered work, such as asking whether promoters or enhancers carry specific chromatin marks, gene-list ORA is the wrong tool. ChromCall is an R package for region-based chromatin enrichment analysis that assigns chromatin status to predefined genomic regions, and it was applied to glioblastoma ChIP-seq to show that promoters of treatment-resistance genes carry a specific histone mark combination [6]. Use this family of methods when your question is about the genomic regions themselves rather than the genes nearby. ## Practical Tips These tips address the recurring operational questions in genomic annotation enrichment analysis in R. **Run all three GO ontologies separately.** BP, MF, and CC terms answer different questions and mixing them in one plot confuses readers. `enrichGO(ont = "ALL")` is available for exploration but split the ontologies for figures. **Set cutoffs to 1 during exploration.** `pvalueCutoff = 1, qvalueCutoff = 1` returns everything, so you can see the shape of the result before deciding where to cut. Filter in the table, not inside the call, when you are still deciding. **Use `simplify()` on GO BP results.** GO has hierarchical redundancy, so you often get five versions of the same process. `ego_simplified <- simplify(ego, cutoff = 0.7)` removes semantically redundant terms based on GO similarity. **Check the gene count per term before interpreting anything.** A term supported by 2 of your 19 genes is a weak lead even with a small p-value. The dotplot gene count channel exists for exactly this reason. **Record session information for reproducibility.** ```r sessionInfo() ``` Annotation databases update over time, and a result that reproduces today may not reproduce next year against a refreshed KEGG release. Save `sessionInfo()` output alongside your results and freeze annotation versions for publications. **Keep the mapping table.** The output of `bitr()` tells you which input genes failed to map. Those dropped genes are part of your analysis story and belong in supplementary material. ## Common Mistakes and Limitations **Using the wrong background universe.** This is the single most damaging error. If your gene list came from a probe array, a capture panel, or expressed genes only, the whole-genome background produces inflated significance. Always pass `universe =` explicitly and state the universe in your methods section. **Mismatched identifier types.** Mixing Ensembl gene IDs, Entrez IDs, and symbols in one vector silently drops the un-matched ones. Run `bitr()` first, check the mapping rate, and verify counts at every step: `length(gene_symbols)` versus `nrow(gene_map)` versus the total genes inside the enrichment result. **Treating GSEA and ORA as interchangeable.** ORA on a significance-thresholded list discards ranking information and depends on your cutoff. If reviewers ask whether results are cutoff-sensitive, the answer comes from a preranked GSEA on the full ranking [4]. **Over-interpreting annotation transfer in non-model species.** As noted in the non-model organism protocol, homology-based annotation transfer introduces errors and misses species-specific biology [11]. Enrichment on transferred annotations is enrichment on transferred annotations, not on curated species knowledge. Report the annotation pipeline and version. **Ignoring multiple testing correction or citing raw p-values.** With thousands of GO terms tested, dozens of terms will appear significant at raw p less than 0.05 by chance alone. Report `p.adjust`, name the correction method, and match the practice in the published literature, where FDR correction is what separated real findings from noise in the soybean meta-QTL enrichment [1]. **Reading GO terms too literally.** GO enrichment of "response to salicylic acid" in soybean SDS resistance MQTLs supported a defense-signaling interpretation [1], but the enrichment tells you the term is overrepresented, not that the mechanism is confirmed. Enrichment is hypothesis-generating. Follow-up validation, whether expression assays or functional studies, is a separate project. **Running enrichment on very small gene lists.** With fewer than about 10 genes, ORA has almost no power and GSEA is unstable. If your candidate list is 6 genes, consider presenting the individual gene annotations in a table instead of p-values that cannot be robust. **Circularity between list and universe.** If you selected your gene list using the same annotation you are testing enrichment against, the result is partly tautological. Define the list from one data source and test against an independent annotation layer where possible. **Forgetting that KEGG is an online resource.** `enrichKEGG()` downloads current pathway definitions at run time, so results can shift between runs when KEGG updates. For archival reproducibility, save the returned object with `saveRDS()`. A final limitation applies to every step in this guide: enrichment analysis describes associations between your gene set and curated annotations. It does not establish causation, and candidate genes emerging from these analyses require experimental validation before biological conclusions are drawn. ## Frequently Asked Questions ### What is the difference between ORA and GSEA in clusterProfiler? ORA tests whether a discrete gene list contains more members of a term than expected from the background, while GSEA tests whether genes from a term cluster at the extremes of a full ranked list. Use ORA for gene lists and preranked GSEA when you have fold changes or statistics for all genes. ### What background universe should I use for enrichment analysis? Use the set of genes that actually had a chance to enter your list, such as all genes on your array, all expressed genes in your RNA-seq run, or all genes mappable from your variants. Passing no universe defaults to all genes in the annotation database, which is only correct for a truly genome-wide list. ### Can I do genomic annotation enrichment analysis in R for a non-model organism? Yes. Build a GMT file and TERM2GENE table from your own annotation pipeline, for example BLASTP plus InterProScan plus eggNOG-mapper, and use `enricher()`, or install a custom annotation package that works with `clusterProfiler`. Published workflows demonstrate both routes for shrimp and Mycobacterium. ### Why do I get no significant enrichment terms? The most common causes are a gene list that is too small, a background universe that matches the list too closely, a cutoff that is too strict, or an annotation database with poor coverage for your organism. Relax cutoffs during exploration, check the mapping rate, and confirm your organism has annotation in the resource you queried. ### Which gene identifier type should I use? Entrez Gene IDs are the safest default because KEGG and most `OrgDb` resources key on them directly. Whatever you choose, convert all inputs to one type with `bitr()` before running the analysis and report the mapping rate. ### How do I handle hundreds of significant GO terms? Apply `simplify()` to collapse semantically redundant terms, filter to terms with at least 5 to 10 supporting genes, and group related terms using `emapplot()`. Present a focused table of 10 to 20 terms rather than the full result. ### Is clusterProfiler the only tool for enrichment analysis in R? No. ChromTag provides ORA and GSEA through an R Shiny interface for chromatin data, SNPannotator handles variant-level annotation and enrichment for GWAS results, and specialized platforms such as ShEnrich wrap `clusterProfiler` for specific organism groups. `gprofiler2` and `limma::goana` are additional options, but `clusterProfiler` remains the most widely used engine. ### How do I report enrichment results in a paper? Report the analysis type, the annotation database and version, the background universe, the identifier type, the multiple testing correction method, and the significance threshold. Provide the full unfiltered results table as supplementary data and cite adjusted p-values in figures and text. ## Sources 1. [Meta-QTL Analysis Reveals Consensus Genomic Regions and Candidate Genes for Resistance to Sudden Death Syndrome in Soybean.](https://pubmed.ncbi.nlm.nih.gov/42739458/) 2. [ShEnrich: a functional annotation database and enrichment platform for penaeid shrimp.](https://pubmed.ncbi.nlm.nih.gov/42963765/) 3. [Genomic insights into M. phaseolina: Multifactorial virulence via secreted effectors, carbohydrate enzymes and secondary metabolites.](https://pubmed.ncbi.nlm.nih.gov/42801983/) 4. [ChromTag: an interactive R-shiny platform for analysis and visualization of CUT&Tag and ChIP-seq peak profiling data.](https://pubmed.ncbi.nlm.nih.gov/42428152/) 5. [SNPannotator: automated functional annotation of genetic variants and linked proxies.](https://pubmed.ncbi.nlm.nih.gov/42635531/) 6. [ChromCall: assigning chromatin status to defined genomic regions using epigenomic profiling data.](https://pubmed.ncbi.nlm.nih.gov/42178345/) 7. [Genome-wide association identifies and validates genomic region controlling grain yield and agronomic traits in extra-early orange maize inbred lines under drought.](https://pubmed.ncbi.nlm.nih.gov/42448753/) 8. [Genome-Wide Association Study and Genomic Selection for Average Daily Gain in Ashidan Yak.](https://pubmed.ncbi.nlm.nih.gov/42651857/) 9. [Integrative Multi-Omics-Informed Analysis Identifies Insulin Therapy-Associated Loci Implicated in Osteoarthritis Susceptibility.](https://pubmed.ncbi.nlm.nih.gov/42644749/) 10. [Pan-genomic analysis reveals ecological adaptation and biocontrol of Bacillus pumilus.](https://pubmed.ncbi.nlm.nih.gov/42435517/) 11. [Enhanced RNA-Seq Expression Profiling and Functional Enrichment in Non-model Organisms Using Custom Annotations.](https://pubmed.ncbi.nlm.nih.gov/42367215/)