8  Functional Enrichment

The Human Genome Project produced the first reference sequence for most of the human genome (Lander et al., 2001). Modern reference assemblies and annotations have changed substantially since that draft, so exact gene counts and genome fractions must always be tied to an assembly and annotation release. Many proteomic, genetic, and metabolic analyses now produce a ranked list of interesting biomolecules — genes ranked by GWAS or differential-expression statistics, microarray probesets, or small molecules from LC-MS metabolomics.

Such ranked lists lack structure and context. Typical questions to answer from a list of biomarkers are (Tipney & Hunter, 2010): what does this gene and its protein product do; how and where; does it make sense on this list; does it interact with other genes/proteins; does its behaviour change during disease or therapy? Answering these gene-by-gene from the literature is cumbersome and biased. Functional enrichment analysis addresses this systematically.

By the end of this laboratory you will be able to:

• annotate a ranked gene list with symbols and names using org.Hs.eg.db
• perform GO and KEGG over-representation analysis with clusterProfiler
• visualise enrichment results (dotplots, concept networks)
• run network-based metabolomics enrichment by diffusion with FELLA

In this lab we perform two functional analyses: one with gene-expression data and one with metabolomics (LC-MS) data.

NoteRequirements

The gene-annotation step runs with org.Hs.eg.db (using resultTable.Rda bundled in the repo under data/). The GO/KEGG enrichment (clusterProfiler) and the metabolomics diffusion (FELLA) need packages and live KEGG downloads that are not available at build time, so those chunks are shown but not executed here. Install the packages with BiocManager::install("clusterProfiler") / BiocManager::install("FELLA") to run them yourself.

NoteLab at a glance

Goal: turn a ranked gene list into biological meaning (GO, KEGG, GSEA). Prerequisites: Differential Expression; Searching NCBI & Ensembl. Time: ~2 h. Packages: org.Hs.eg.db (runs); clusterProfiler, enrichplot, fgsea, FELLA (install to run). Deliverable: a lab report.

8.1 Functional analysis in differential expression

8.1.1 Gene set functional enrichment

Load a sample results table (an .Rda from the course campus), similar to the results from the differential-expression lab:

library(org.Hs.eg.db)
load("data/resultTable.Rda")
resultTable

This small teaching table includes only 10 results ranked by p-value, with the log fold change (\(\log_2 \mathrm{FC} = \log_2 B - \log_2 A\)) and Entrez Gene IDs as row names. Ten selected hits are enough to practise identifier mapping, but usually too few for a well-powered enrichment analysis.

Tip✅ Exercise 8.1 — Annotation

Use select() and merge() to create a new data frame with SYMBOL and gene-name information.

resultTable$ENTREZID <- rownames(resultTable)
annots <- AnnotationDbi::select(org.Hs.eg.db, keys = resultTable$ENTREZID,
                                columns = c("SYMBOL", "GENENAME"), keytype = "ENTREZID")
resultTable <- merge(resultTable, annots, by = "ENTREZID", all.x = TRUE)
gene_ids <- unique(resultTable$ENTREZID)
head(resultTable)

8.1.1.1 GO enrichment analysis

Install and load clusterProfiler, which tests a selected gene list against a defined background universe. In a real analysis, background_ids must contain all genes that could have been selected—for example, all genes retained after expression filtering—not every annotated human gene. The broader fallback below is included only because this teaching file does not contain the original tested universe:

library(clusterProfiler)
background_ids <- keys(org.Hs.eg.db, keytype = "ENTREZID") # replace when available
egoBP <- enrichGO(gene         = gene_ids,
                  universe     = background_ids,
                  OrgDb        = org.Hs.eg.db,
                  pAdjustMethod = "BH",
                  pvalueCutoff = 0.05,
                  ont          = "BP",
                  keyType      = "ENTREZID",
                  readable     = TRUE)
head(as.data.frame(egoBP))

Dotplots are a good graphical summary of the over-representation data:

dotplot(egoBP)

Tip✅ Exercise 8.2 — GO

1. Draw a dot plot or enrichment map and interpret the gene ratio, adjusted p-value, and gene count.
2. Obtain over-represented GO terms for the Cellular Component and Molecular Function ontologies. Interpret the results.

8.1.1.2 KEGG enrichment analysis

Repeat the analysis for over-represented KEGG pathways. The KEGG.db annotation package has not been updated since 2012 (when KEGG went partly paywalled); enrichKEGG downloads the latest online KEGG data for the analysis:

ekegg <- enrichKEGG(gene         = gene_ids,
                    organism     = "hsa",
                    pvalueCutoff = 0.05)
head(as.data.frame(ekegg))

The relationship between KEGG pathways and the input gene list can be plotted with cnetplot() (a concept network); enrichplot also gives dot plots and enrichment maps:

library(enrichplot)
dotplot(ekegg)                       # top pathways by gene ratio
cnetplot(ekegg, categorySize = "pvalue")  # pathway–gene network
emapplot(pairwise_termsim(egoBP))    # term–term similarity map

8.1.2 ORA vs GSEA

The analyses above are over-representation analysis (ORA): they ask whether a cut list of top genes is enriched for a term. Gene-set enrichment analysis (GSEA) instead uses the whole ranked list (e.g. by log fold change), so it can detect coordinated shifts of many small-effect genes that never pass a cutoff. GSEA requires the complete, uniquely named ranking—not this file of 10 preselected hits. Given a full differential-expression table called full_results, the workflow is:

# rank every tested gene by a signed statistic; names must be unique Entrez IDs
stopifnot(all(c("ENTREZID", "statistic") %in% names(full_results)))
full_results <- full_results[!duplicated(full_results$ENTREZID) &
                             !is.na(full_results$ENTREZID), ]
geneList <- sort(setNames(full_results$statistic, full_results$ENTREZID),
                 decreasing = TRUE)
gse <- gseGO(geneList, OrgDb = org.Hs.eg.db, ont = "BP", keyType = "ENTREZID")
ridgeplot(gse)      # distribution of ranks per enriched set
gseaplot2(gse, geneSetID = 1)

Tip✅ Exercise 8.3 — ORA vs GSEA

Using a complete ranked result table, run both ORA and GSEA. Where do they agree, and where does GSEA recover pathways ORA misses? Explain why GSEA cannot be validly run on the provided top-10 table.

TipSolution 8.3

GSEA needs the full, uniquely-named ranked vector of every tested gene. With only a top-10 (already-thresholded) list the ranked distribution is incomplete, the running-sum statistic is meaningless, and duplicated IDs break the ranking — so GSEA cannot be validly run on it. Agreement appears for the strongest pathways; GSEA also recovers pathways driven by many small-effect genes that never survive the ORA cutoff.

8.2 Enrichment in metabolomics with FELLA

This case study examines the effects of environmental contamination on gilt-head bream (Sparus aurata) (Ziarrusta et al., 2018).¹ Fish were exposed to oxybenzone over 14 days, and changes were sought in brain, liver, and plasma using untargeted metabolomics (UHPLC-qOrbitrap MS, positive and negative modes, C18 and HILIC separation). Mortality and brain metabolites were unchanged, but liver and plasma showed perturbations — evidence of adverse effects beyond the well-studied hormonal activity.

We use network-based enrichment by diffusion with Bioconductor’s FELLA (Picart-Armada et al., 2017). FELLA starts from a list of affected metabolites and introduces a comprehensive output relating the top-ranked pathways to other biological entities, through a diffusion process on a holistic graph of the KEGG database. FELLA needs the KEGG graph and complementary data (a FELLA.DATA S4 object) and a list of affected metabolites as KEGG compounds (a FELLA.USER object).

8.2.1 Building a KEGG object

KEGG has no pathway annotations for Sparus aurata (Kanehisa et al., 2017), so it is common to use zebrafish (Danio rerio, organism code dre) as an approximation:

library(FELLA)
library(igraph)
library(magrittr)
tmpdir <- paste0("./kegg_database")

# Filter the dre01100 overview pathway, as in the article
graph <- buildGraphFromKEGGREST(organism = "dre", filter.path = c("01100"))

# Ensure the database does not already exist (FELLA will not overwrite it)
unlink(tmpdir, recursive = TRUE)
buildDataFromGraph(keggdata.graph = graph,
                   databaseDir = tmpdir, internalDir = FALSE,
                   matrices = "none", normality = "diffusion", niter = 100)

Load the FELLA.DATA object, which represents the relationships among KEGG entities (Pathways, Modules, Enzymes, Reactions, Compounds):

fella.data <- loadKEGGdata(databaseDir = tmpdir, internalDir = FALSE, loadMatrix = "none")
fella.data

8.2.2 Enrichment analysis on liver tissue

Find the KEGG identifiers associated with metabolic changes in liver tissue from (Ziarrusta et al., 2018), build a vector cpd.liver of KEGG compound IDs, and enrich it against fella.data:

cpd.liver <- c("C12623", "C01179", "C05350", "C05598", "C01586")

analysis.liver <- enrich(compounds = cpd.liver, data = fella.data,
                         method = "diffusion", approx = "normality")

Check whether all input metabolites were mapped:

analysis.liver %>% getInput %>% getName(data = fella.data)

Plot the resulting network of entities found by FELLA:

plot(analysis.liver, method = "diffusion", data = fella.data,
     nlimit = 250, plotLegend = FALSE)

Retrieve the sub-network as an igraph object, or as a table:

g.liver <- generateResultsGraph(object = analysis.liver, data = fella.data,
                                method = "diffusion")
tab.liver <- generateResultsTable(object = analysis.liver, data = fella.data,
                                  method = "diffusion")
head(tab.liver)

Tip✅ Exercise 8.4 — Plasma metabolites

1. Repeat the analysis with the plasma metabolites reported in (Ziarrusta et al., 2018).
2. Discuss the enriched entities in relation to the paper’s goals.

8.2.3 FELLA graphical user interface

FELLA also ships a graphical user interface (based on shiny, launched with launchApp()) with tabs for Compounds upload (submit the metabolite list and see mapped vs mismatching compounds), Advanced options (database, method and approximation choice, parameter tweaking, GO label for semantic similarity), Results (an interactive sub-graph of the top-k KEGG entries with a linked table), and Export (several tabular and network export options).

8.3 Further reading

• Yu et al. (2012), OMICS — clusterProfiler for ORA/GSEA.
• Subramanian et al. (2005), PNAS — the original GSEA method.
• The FELLA vignette for the metabolomics diffusion analysis.

Worked solutions for the exercises are in the collapsed callouts below each exercise.

8.4 Report

Write a laboratory report following the structure of this text. It should contain an introduction, the goals of each section, and a description of the results of each part. Include any graphics relevant to making your work understandable. The lab should be delivered within two weeks.

De Castro-Català, N., Muñoz, I., Riera, J. L., & Ford, A. T. (2017). Evidence of low dose effects of the antidepressant fluoxetine and the fungicide prochloraz on the behavior of the keystone freshwater invertebrate gammarus pulex. Environmental Pollution (Barking, Essex : 1987), 231, 406–414. https://doi.org/10.1016/j.envpol.2017.07.088

Kanehisa, M., Furumichi, M., Tanabe, M., Sato, Y., & Morishima, K. (2017). KEGG: New perspectives on genomes, pathways, diseases and drugs. Nucleic Acids Research, 45, D353–D361. https://doi.org/10.1093/nar/gkw1092

Lander, E. S., Linton, L. M., Birren, B., Nusbaum, C., Zody, M. C., Baldwin, J., Devon, K., Dewar, K., Doyle, M., FitzHugh, W., Funke, R., Gage, D., Harris, K., Heaford, A., Howland, J., Kann, L., Lehoczky, J., LeVine, R., McEwan, P., … Consortium, I. H. G. S. (2001). Initial sequencing and analysis of the human genome. Nature, 409, 860–921. https://doi.org/10.1038/35057062

Picart-Armada, S., Fernández-Albert, F., Vinaixa, M., Rodríguez, M. A., Aivio, S., Stracker, T. H., Yanes, O., & Perera-Lluna, A. (2017). Null diffusion-based enrichment for metabolomics data. PloS One, 12, e0189012. https://doi.org/10.1371/journal.pone.0189012

Tipney, H., & Hunter, L. (2010). An introduction to effective use of enrichment analysis software. Human Genomics, 4, 202–206.

Ziarrusta, H., Mijangos, L., Picart-Armada, S., Irazola, M., Perera-Lluna, A., Usobiaga, A., Prieto, A., Etxebarria, N., Olivares, M., & Zuloaga, O. (2018). Non-targeted metabolomics reveals alterations in liver and plasma of gilt-head bream exposed to oxybenzone. Chemosphere, 211, 624–631. https://doi.org/10.1016/j.chemosphere.2018.08.013

---

1. Low concentrations of toxic products such as antidepressants in rivers can change swimming and feeding behaviour in some creatures. Such products reach the water environment — they are not fully removed in wastewater treatment — and can accumulate in fish brains (De Castro-Català et al., 2017).