<environment: namespace:Rcpp>10 GO enrichment analysis
GO comprises three orthogonal ontologies, i.e. molecular function (MF), biological process (BP), and cellular component (CC).
10.1 Chapter overview
| Aspect | GO enrichment workflow |
|---|---|
| Questions | Which GO terms are over-represented in a selected gene set or enriched across a ranked profile? How can genes be classified by ontology level? |
| Input | Gene IDs compatible with an OrgDb (commonly Entrez), a named ranked vector for GSEA, or custom GO TERM2GENE annotations. |
| Methods | ORA, GSEA, GO-level classification, custom annotation enrichment, and semantic redundancy reduction. |
| Main functions | groupGO(), enrichGO(), gseGO(), buildGOmap(), enricher(), GSEA(), simplify(). |
| Output | Classification tables and enrichResult / gseaResult objects ready for enrichplot. |
| Main limitations | Coverage depends on the OrgDb release; GO is hierarchical, so ancestor propagation and ontology choice affect interpretation. Use a defensible universe and include ancestor annotations in custom mappings. |
Confirm identifier and background-set choices first. For custom GO maps, include ancestor terms where required before testing.
10.2 Supported organisms
GO analyses (groupGO(), enrichGO() and gseGO()) support organisms that have an OrgDb object available (see also session 2.2).
If a user has GO annotation data (in a data.frame format with the first column as gene ID and the second column as GO ID), they can use the enricher() and GSEA() functions to perform an over-representation test and gene set enrichment analysis. Users can also read GO annotation from GMT files using the read.gmt() function, which parses the file into a data.frame suitable for these functions.
If the genes are annotated by direct annotation, they should also be annotated by their ancestor GO nodes (indirect annotation). If a user only has direct annotation, they can pass their annotation to the buildGOmap function, which will infer indirect annotation and generate a data.frame that is suitable for both enricher() and GSEA().
10.3 Direct vs indirect GO annotation
One common issue users encounter is the discrepancy between direct gene-to-GO mapping and the results from enrichment analysis. For instance, a user might find that a specific GO term is enriched, but when they check the input gene list against the direct annotations, the term doesn’t appear to be associated with many genes.
This happens because GO structure is a Directed Acyclic Graph (DAG). If a gene is annotated to a specific term (e.g., “glucose metabolic process”), it is implicitly annotated to all its ancestor terms (e.g., “metabolic process”). enrichGO and gseGO account for this indirect annotation by propagating annotations up the graph.
If you are providing your own annotation file (e.g., using enricher or GSEA), you must ensure that it includes these indirect annotations. If your annotation file only contains direct annotations, you can use the buildGOmap() function to generate the full annotation set (direct + indirect) before performing enrichment analysis.
# Assuming 'gomap' is a data.frame with direct annotations (GeneID, GOID)
# Build indirect annotations
gomap_with_ancestors <- buildGOmap(gomap)
# Use the result for enrichment
enricher(gene, TERM2GENE = gomap_with_ancestors, ...)This step is crucial for accurate GO enrichment analysis when not using OrgDb objects.
10.4 GO classification
In clusterProfiler, the groupGO() function is designed for gene classification based on GO distribution at a specific level. Here we use the dataset geneList provided by DOSE.
library(clusterProfiler)
library(org.Hs.eg.db)
data(geneList, package="DOSE")
gene <- names(geneList)[abs(geneList) > 2]
# Entrez gene ID
head(gene)[1] "4312" "8318" "10874" "55143" "55388" "991" ggo <- groupGO(gene = gene,
OrgDb = org.Hs.eg.db,
ont = "CC",
level = 3,
readable = TRUE)
head(ggo) ID Description Count GeneRatio geneID
GO:0000133 GO:0000133 GO:0000133 0 0/207
GO:0000417 GO:0000417 HIR complex 0 0/207
GO:0000796 GO:0000796 condensin complex 2 2/207 NCAPH/NCAPG
GO:0000808 GO:0000808 origin recognition complex 0 0/207
GO:0000930 GO:0000930 gamma-tubulin complex 0 0/207
GO:0000939 GO:0000939 inner kinetochore 0 0/207 The gene parameter is a vector of gene IDs (can be any ID type that is supported by the corresponding OrgDb, see also Section 18.1). If readable is set to TRUE, the input gene IDs will be converted to gene symbols.
10.5 GO over-representation analysis
The clusterProfiler package implements enrichGO() for gene ontology over-representation test.
ego <- enrichGO(gene = gene,
universe = names(geneList),
OrgDb = org.Hs.eg.db,
ont = "CC",
pAdjustMethod = "BH",
pvalueCutoff = 0.01,
qvalueCutoff = 0.05,
readable = TRUE)
head(ego) ID Description GeneRatio
GO:0000775 GO:0000775 chromosome, centromeric region 21/201
GO:0000779 GO:0000779 condensed chromosome, centromeric region 18/201
GO:0072686 GO:0072686 mitotic spindle 17/201
GO:0000793 GO:0000793 condensed chromosome 21/201
GO:0098687 GO:0098687 chromosomal region 25/201
GO:0000776 GO:0000776 kinetochore 17/201
BgRatio RichFactor FoldEnrichment oddsRatio zScore pvalue
GO:0000775 199/11914 0.10552764 6.255006 7.560393 9.792716 1.864198e-11
GO:0000779 150/11914 0.12000000 7.112836 8.629657 9.869278 6.484361e-11
GO:0072686 137/11914 0.12408759 7.355122 8.925770 9.800345 1.316552e-10
GO:0000793 222/11914 0.09459459 5.606965 6.681924 9.076565 1.486538e-10
GO:0098687 322/11914 0.07763975 4.601990 5.459902 8.583527 1.732119e-10
GO:0000776 141/11914 0.12056738 7.146466 8.634862 9.617585 2.087156e-10
p.adjust qvalue
GO:0000775 5.443457e-09 5.025287e-09
GO:0000779 9.467167e-09 8.739892e-09
GO:0072686 1.011557e-08 9.338487e-09
GO:0000793 1.011557e-08 9.338487e-09
GO:0098687 1.011557e-08 9.338487e-09
GO:0000776 1.015749e-08 9.377188e-09
geneID
GO:0000775 SKA1/ERCC6L/TOP2A/KIF18A/TTK/NEK2/CCNB1/CDT1/AURKA/CENPE/CENPM/EZH2/CENPN/NCAPG/AURKB/HJURP/CDC20/NDC80/BIRC5/MAD2L1/CDCA8
GO:0000779 SKA1/ERCC6L/KIF18A/TTK/NEK2/CCNB1/CDT1/AURKA/CENPE/CENPM/CENPN/NCAPG/AURKB/HJURP/CDC20/NDC80/BIRC5/MAD2L1
GO:0072686 SKA1/KIF18A/DLGAP5/AURKA/CENPE/KIF18B/NUSAP1/KIF20A/CDK1/PRC1/KIFC1/KIF11/TPX2/ASPM/KIF23/AURKB/MAD2L1
GO:0000793 SKA1/ERCC6L/TOP2A/KIF18A/TTK/NEK2/CCNB1/CDT1/AURKA/CENPE/CENPM/CHEK1/CENPN/NCAPG/AURKB/HJURP/CDC20/NDC80/BIRC5/MAD2L1/NCAPH
GO:0098687 SKA1/ERCC6L/TOP2A/KIF18A/MCM5/TTK/NEK2/CCNB1/CDT1/AURKA/CENPE/CENPM/EZH2/RAD51AP1/CHEK1/CDK1/CENPN/NCAPG/AURKB/HJURP/CDC20/NDC80/BIRC5/MAD2L1/CDCA8
GO:0000776 SKA1/ERCC6L/KIF18A/TTK/NEK2/CCNB1/CDT1/AURKA/CENPE/CENPM/CENPN/AURKB/HJURP/CDC20/NDC80/BIRC5/MAD2L1
Count
GO:0000775 21
GO:0000779 18
GO:0072686 17
GO:0000793 21
GO:0098687 25
GO:0000776 17Any gene ID type that is supported in OrgDb can be directly used in GO analyses. Users need to specify the keyType parameter to specify the input gene ID type.
For enrichGO(), non-ENTREZID inputs use an Entrez-backed representation as a memory-saving internal path when possible. If source IDs cannot be mapped to Entrez, their native keyType GO annotations are added to a hybrid GSON instead of being silently dropped. When no universe is supplied, the effective annotated hybrid background is constructed and reported by the function.
The result is mapped back to the original input ID type, and the result object’s @keytype records that source type. This keeps the user-facing result aligned with the input namespace while retaining the memory protection needed for one-to-many identifiers such as ACCNUM.
For example, if we use ENSEMBL IDs as input:
gene.df <- bitr(gene, fromType = "ENTREZID",
toType = c("ENSEMBL", "SYMBOL"),
OrgDb = org.Hs.eg.db)
ego2 <- enrichGO(gene = gene.df$ENSEMBL,
OrgDb = org.Hs.eg.db,
keyType = 'ENSEMBL',
ont = "CC",
pAdjustMethod = "BH",
pvalueCutoff = 0.01,
qvalueCutoff = 0.05)
head(ego2, 3) ID Description GeneRatio
GO:0000775 GO:0000775 chromosome, centromeric region 21/202
GO:0000779 GO:0000779 condensed chromosome, centromeric region 18/202
GO:0098687 GO:0098687 chromosomal region 25/202
BgRatio RichFactor FoldEnrichment oddsRatio zScore pvalue
GO:0000775 257/19979 0.08171206 8.081808 9.606728 11.54800 1.899777e-13
GO:0000779 190/19979 0.09473684 9.370036 11.150468 11.71558 8.795611e-13
GO:0098687 417/19979 0.05995204 5.929613 6.984679 10.28124 9.528885e-13
p.adjust qvalue
GO:0000775 5.376368e-11 1.223869e-11
GO:0000779 8.988915e-11 2.046223e-11
GO:0098687 8.988915e-11 2.046223e-11
geneID
GO:0000775 ENSG00000154839/ENSG00000186871/ENSG00000131747/ENSG00000121621/ENSG00000112742/ENSG00000117650/ENSG00000134057/ENSG00000087586/ENSG00000138778/ENSG00000100162/ENSG00000106462/ENSG00000166451/ENSG00000109805/ENSG00000178999/ENSG00000123485/ENSG00000117399/ENSG00000167513/ENSG00000080986/ENSG00000089685/ENSG00000164109/ENSG00000134690
GO:0000779 ENSG00000154839/ENSG00000186871/ENSG00000121621/ENSG00000112742/ENSG00000117650/ENSG00000134057/ENSG00000087586/ENSG00000138778/ENSG00000100162/ENSG00000166451/ENSG00000109805/ENSG00000178999/ENSG00000123485/ENSG00000117399/ENSG00000167513/ENSG00000080986/ENSG00000089685/ENSG00000164109
GO:0098687 ENSG00000154839/ENSG00000186871/ENSG00000131747/ENSG00000121621/ENSG00000100297/ENSG00000112742/ENSG00000117650/ENSG00000134057/ENSG00000087586/ENSG00000138778/ENSG00000100162/ENSG00000106462/ENSG00000111247/ENSG00000149554/ENSG00000170312/ENSG00000166451/ENSG00000109805/ENSG00000178999/ENSG00000123485/ENSG00000117399/ENSG00000167513/ENSG00000080986/ENSG00000089685/ENSG00000164109/ENSG00000134690
Count
GO:0000775 21
GO:0000779 18
GO:0098687 25Gene IDs can be mapped to gene Symbols by using the parameter readable=TRUE or setReadable() function.
10.6 GO Gene Set Enrichment Analysis
The clusterProfiler package provides the gseGO() function for gene set enrichment analysis using gene ontology.
ego3 <- gseGO(geneList = geneList,
OrgDb = org.Hs.eg.db,
ont = "CC",
minGSSize = 100,
maxGSSize = 500,
pvalueCutoff = 0.05,
verbose = FALSE)The format of input data, geneList, is documented in the Section 18.1. Beware that only gene sets with size in [minGSSize, maxGSSize] will be tested.
NoteSize filtering is overlap-based, not annotation-based
minGSSize and maxGSSize operate on the overlap of each gene set with the genes you actually provide (the candidate vector for ORA, or the named geneList for GSEA), not on the raw, full annotation size stored in the ontology database. The code therefore does roughly:
gene_sets <- lapply(gene_sets, intersect, names(geneList)) # or universe for ORA
idx <- get_geneSet_index(gene_sets, minGSSize, maxGSSize)This is important for broad libraries such as MSigDB Hallmark collections, Reactome, or the upper HDO Slim. A database gene set annotated with 3 000 genes is not rejected just because 3 000 > maxGSSize = 500; it is kept if its overlap with your 12 000 measured genes has, say, 480 members. The only terms genuinely dropped are (a) very small overlaps (below minGSSize, which usually means the term has too few observed members for Fisher’s exact test or the weighted Kolmogorov-Smirnov statistic to be reliable) and (b) extremely broad overlaps (above maxGSSize, which usually correspond to root-level, biologically uninformative nodes that inflate the multiple-testing burden without adding insight).
See also Tuning minGSSize and maxGSSize.
10.7 Topology-aware GO enrichment
In classical GO enrichment, genes are treated as independent members of a gene set. When a biological network is available, clusterProfiler can now call the enrichit engine through the high-level wrappers nseGO() and mnseGO(). The former performs single-network enrichment, while the latter extends the same idea to a multi-layer setting.
To keep the example reproducible and consistent with the KEGG chapter, we reuse one shared demonstration object.
demo <- clusterprofiler_enrichit_demo()ego_nse <- nseGO(
geneList = demo$geneList_evidence,
network = demo$network,
OrgDb = org.Hs.eg.db,
keyType = "ENTREZID",
ont = "BP",
minGSSize = 5,
maxGSSize = 500,
threshold = 1e-6,
maxIter = 50,
verbose = FALSE,
pvalueCutoff = 1,
method = "sample",
nPerm = 30
)
head(as.data.frame(ego_nse)[, c("ID", "Description", "NES", "p.adjust")], 3) ID Description NES p.adjust
GO:0032496 GO:0032496 response to lipopolysaccharide 3.017036 0.3309958
GO:0009987 GO:0009987 cellular process 3.008683 0.3309958
GO:0045087 GO:0045087 innate immune response 2.875249 0.3309958The mnseGO() function keeps the same GO-aware interface, but accepts a list of layer-specific networks together with explicit inter-layer couplings.
ego_mnse <- mnseGO(
seed_list = demo$seed_list,
networks = list(RNA = demo$network, PROT = demo$network_2),
couplings = demo$couplings,
OrgDb = org.Hs.eg.db,
keyType = "ENTREZID",
ont = "BP",
collapse = "weighted_mean",
layer_weights = c(RNA = 1, PROT = 1.2),
minGSSize = 5,
maxGSSize = 500,
threshold = 1e-6,
maxIter = 50,
verbose = FALSE,
pvalueCutoff = 1,
method = "sample",
nPerm = 30
)
head(as.data.frame(ego_mnse)[, c("ID", "Description", "NES", "p.adjust")], 3) ID Description NES
GO:0051656 GO:0051656 establishment of organelle localization 2.730273
GO:0032496 GO:0032496 response to lipopolysaccharide 2.412481
GO:0009987 GO:0009987 cellular process 2.387504
p.adjust
GO:0051656 0.4033326
GO:0032496 0.4033326
GO:0009987 0.4033326These wrappers preserve the familiar GO semantics (OrgDb, ont, and keyType) while delegating the network propagation and enrichment engine to enrichit.
Both the enrichGO() and gseGO() functions require an OrgDb object as the background annotation. For organisms that don’t have OrgDb provided by Bioconductor, users can query one (if available) online via AnnotationHub. If there is no OrgDb available, users can obtain GO annotation from other sources, e.g. from biomaRt, or annotate the genes using Blast2GO or the Trinotate pipeline. Then the enricher() or GSEA() functions can be used to perform GO analysis for these organisms, similar to the examples using WikiPathways and MSigDB. Another solution is to create an OrgDb on your own using the AnnotationForge package.
Here is an example of querying GO annotation from Ensembl using biomaRt.
library(biomaRt)
ensembl <- useEnsemblGenomes(biomart = "plants_mart", dataset = "nattenuata_eg_gene")
gene2go <- getBM(attributes =c("ensembl_gene_id", "go_id"), mart=ensembl)Alternatively, you can use AnnotationHub to query and retrieve OrgDb objects for a wide range of organisms, including non-model species.
Example: GO enrichment analysis for Maize (Zea mays)
library(AnnotationHub)
hub <- AnnotationHub()
# Query for Zea mays OrgDb
query(hub, "zea")
# Retrieve the OrgDb (e.g., AH55736, assume it's the correct record)
maize <- hub[['AH55736']]
# Check keys
length(keys(maize))
columns(maize)
# Perform enrichment analysis
# Assuming 'sample_genes' is a vector of Entrez IDs
res <- enrichGO(sample_genes, OrgDb=maize, pvalueCutoff=0.05, qvalueCutoff=0.05)This approach allows you to use the standard enrichGO workflow with organisms not included in the default Bioconductor OrgDb packages.
10.8 Visualization of GO enrichment results
The clusterProfiler package provides several visualization methods to help interpret enrichment results, including barplot, dotplot, cnetplot, emapplot, goplot and plotGOgraph.
Please refer to the Chapter 26 chapter for details.
10.9 Filtering and simplifying GO terms
GO enrichment analysis often results in many redundant terms due to the hierarchical structure of GO. Parent terms and child terms often share a large proportion of genes, leading to multiple significant results that represent similar biological processes.
10.9.1 Removing specific terms or levels
The dropGO() function can be used to remove specific GO terms or GO levels from the results. The enrichGO() function tests the whole GO corpus and enriched result may contains very general terms. If users want to restrict the result to a specific GO level (e.g., level 3 or 4), they can use the gofilter() function. Both functions work with results obtained from enrichGO(), gseGO() and compareCluster().
To see the level of each term rather than filter by a level you already know, add_go_level() appends it as a level column, where level 1 is the ontology root (GO:0008150 for BP, GO:0005575 for CC, GO:0003674 for MF). This makes an arbitrary band easy to keep:
x <- enrichGO(gene, OrgDb = "org.Hs.eg.db", ont = "BP")
x <- add_go_level(x)
subset(x, level >= 3 & level <= 6)Note that GO is a directed acyclic graph, so a term can be reachable at more than one depth; the shallowest level is reported. The level is resolved from the ontology the term itself belongs to, so an ont = "ALL" result (or a compareCluster() result) is labelled term by term.
10.9.2 Simplifying enriched GO terms
To reduce redundancy, clusterProfiler provides a simplify method. It uses semantic similarity (via GOSemSim) to calculate the similarity between enriched terms and removes highly similar terms, keeping the most significant one.
The simplify method applies select_fun (which can be a user-defined function) to the feature specified by by (e.g., p.adjust) to select one representative term from redundant terms (which have similarity higher than cutoff). This results in a cleaner and more interpretable visualization.
A complementary strategy is Bayesian term selection, which uses term-gene coverage to rank terms by posterior explanatory support instead of semantic similarity alone (see Chapter 24).
In Figure 10.1(A), we can found that there are many redundant terms form a highly condense network. After removing redundant terms using the simplify() method, the result is more clear to view the whole story.
#|
library(enrichplot)
data(geneList, package="DOSE")
de <- names(geneList)[abs(geneList) > 2]
bp <- enrichGO(de, ont="BP", OrgDb = 'org.Hs.eg.db')
bp <- pairwise_termsim(bp)
bp2 <- simplify(bp, cutoff=0.7, by="p.adjust", select_fun=min)
p1 <- emapplot(bp)
p2 <- emapplot(bp2)
plot_list(p1, p2, ncol=2, tag_levels = 'A')
emapplot() function. (A) original result. (B) simplify result.
Alternatively, users can use slim version of GO and use the enricher() or gseGO() functions to analyze.
Note on enricher results: The simplify() method only needs to know which ontology the terms belong to. For enrichGO() and gseGO() results that ontology is recorded in the object, but a result built with enricher() (e.g. from a custom annotation or an MSigDB GO collection) does not carry one. In that case simplify() now derives it from the GO IDs themselves, so no manual step is required:
# 'res' is the result from enricher() using a GO gene-set collection
res_simplified <- simplify(res)The inference requires that the gene-set IDs really are GO terms (GO:0008150, …): if a single ontology is represented it is used directly, and if the terms span several ontologies they are simplified as GOALL. A result whose gene sets are not GO terms is still rejected, because simplify() reduces redundancy through GO semantic similarity. Setting the slot explicitly (res@ontology <- "BP") still works, and remains necessary only when the IDs are not GO terms.
10.10 Troubleshooting
10.10.1 No results found
It is common to receive a message saying “No gene set have size > 10 … return NULL” or finding no significant terms.
- Gene Set Size: The default
minGSSizeis 10 andmaxGSSizeis 500. Remember that filtering is based on the overlap with your provided input, not the raw database annotation size (see TuningminGSSizeandmaxGSSizeand the overlap-rule callout above). A term whose full annotation has 3 000 genes may still pass if only 200 of those genes were measured in your experiment. If no term passes both cutoffs, enlarge the window with, for example,minGSSize = 5, maxGSSize = 5000, or setmaxGSSize = Infto disable the upper cap entirely.
10.10.2 Tuning minGSSize and maxGSSize
The default maxGSSize = 500 is a conservative first-pass safety net, not a hard recommendation. It was chosen to avoid root-level GO/HDO nodes and extremely broad Hallmark-style gene sets overwhelming the multiple-testing correction in a typical short list of a few hundred differentially expressed genes. It also used to cap raw annotation sizes in older code paths, but that restriction no longer applies: size filtering now happens after intersecting with your measured genes.
The table below is a pragmatic starting point. Always verify the term list after the first run and adjust the window to the characteristic scale of your annotation library and input list.
| Scenario | Suggested minGSSize |
Suggested maxGSSize |
|---|---|---|
| Standard GO BP or HDO term-level analysis on ~200 DEGs | 5–10 | 500 (default) |
| GO Slim, DO Slim or broad candidates / Reactome | 10–20 | 2000–5000 |
| Custom GMT, MSigDB Hallmark / C2, or user-defined sets | 10 | Inf (disable cap) |
| You want to see everything and filter the result later | 1 | Inf |
Setting maxGSSize = Inf is safe and supported: get_geneSet_index() treats NA or NULL as Inf as well. If after running with wide bounds you see obviously uninformative root terms such as biological_process or broad “positive regulation of …” nodes, filter them out afterwards with dplyr::filter() or re-run with a narrower window rather than using a tight global default that may hide relevant biology.
For reproducibility, document the values you chose together with the annotation release and package versions; the Appendix reproducibility record lists this explicitly. 2. Universe: If you provide a custom universe (background gene list), ensure it is large enough. A small universe can make it statistically difficult to find enrichment. 3. P-value Adjustment: clusterProfiler uses multiple testing correction (e.g., BH) by default. Some web tools might use raw p-values or different thresholds, leading to more (but potentially false positive) results. clusterProfiler tends to be more conservative and reliable.
10.10.3 Annotation quality
Sometimes, the OrgDb provided by Bioconductor might be outdated compared to the latest online databases (e.g., TAIR for Arabidopsis).
- Check Dates: Always check the metadata of the
OrgDbobject to see the source and date of the annotation. - Custom Annotation: If the
OrgDbis too old, consider downloading the latest annotation (e.g., GAF or GOSLIM files) from the organism’s primary database and usingenricher()/GSEA()withbuildGOmap()as described in the Direct vs indirect GO annotation section.
10.11 Summary
GO semantic similarity can be calculated by GOSemSim (Yu et al. 2010). We can use it to cluster genes/proteins into different clusters based on their functional similarity and can also use it to measure the similarities among GO terms to reduce the redundancy of GO enrichment results.
10.12 Next steps
- For online annotation refreshes, downloads, and release metadata, see external-service troubleshooting.
- For custom annotations or unsupported organisms, see universal enrichment.
- To reduce redundancy, continue to GOSemSim and Bayesian term selection.
- For plots, use the visualization map.
References
Yu, Guangchuang, Fei Li, Yide Qin, Xiaochen Bo, Yibo Wu, and Shengqi Wang. 2010. “GOSemSim: An r Package for Measuring Semantic Similarity Among GO Terms and Gene Products.” Bioinformatics 26 (7): 976–78. https://doi.org/10.1093/bioinformatics/btq064.