Functional Enrichment: GO, KEGG, and GSEA

From a gene list to biology

Differential expression often hands you hundreds of changed genes. Reading them one by one is hopeless, and your eye will only recognise the few you already know. Functional enrichment solves this by grouping genes by what they do, then asking a statistical question: is any biological function or pathway showing up in your list far more than you would expect by chance? If "immune response" genes are all over your hits, that is a strong clue about what your experiment did.

Two ways to ask the question

There are two common approaches, and they suit different situations.

• Over-representation analysis (ORA). You take your list of significant genes (your hits) and ask: is this list unusually full of genes from a particular function or pathway, compared to all genes? It is simple and intuitive, but it throws away everything below your significance cutoff.
• Gene Set Enrichment Analysis (GSEA). You keep all genes, ranked from most up to most down, and ask: does any group of related genes sit together near the top or the bottom of the ranking? It catches coordinated, subtle shifts that ORA can miss.

1Set up and get a ranked gene list

Enrichment needs real gene IDs (so the tools can look up what each gene does), which our small synthetic dataset does not have. So we'll use a ready-made example list of real human genes that ships with the analysis package. It is a named vector: each gene's ID paired with its fold change, sorted from most up to most down. That is exactly the shape your own DESeq2 results take once you map them (we show how at the end).

# clusterProfiler is the standard R toolkit for enrichment
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install(c("clusterProfiler", "org.Hs.eg.db", "DOSE", "enrichplot"), update = FALSE, ask = FALSE)
library(clusterProfiler); library(org.Hs.eg.db); library(enrichplot)

# a ranked list of real human genes (named vector: gene ID -> fold change)
data(geneList, package = "DOSE")
head(geneList)

The names are Entrez IDs, a standard numeric code for each gene. org.Hs.eg.db is the human annotation database that lets the tools translate those IDs into functions and pathways. (The "Hs" means Homo sapiens; there are equivalent packages for mouse, and other species.)

2Over-representation: GO terms

First the simple approach. We take the clearly significant genes and ask which Gene Ontology terms are over-represented.

# treat genes with a large fold change as our "hits"
sig_genes <- names(geneList)[abs(geneList) > 2]

ego <- enrichGO(gene          = sig_genes,
                OrgDb         = org.Hs.eg.db,
                ont           = "BP",        # Biological Process
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.05,
                readable      = TRUE)       # show gene symbols, not IDs

head(ego)
dotplot(ego, showCategory = 12) + ggtitle("Over-represented biological processes")

The dot plot is the payoff: each row is a biological process, the dot size is how many of your genes fall in it, and the colour is the adjusted p-value. Processes near the top are the ones your changed genes are concentrated in.

3GSEA: the whole ranked list

Now the more powerful approach. GSEA uses every gene, in ranked order, and finds gene sets that trend up or down as a group, even if no single gene was a dramatic hit.

gse <- gseGO(geneList     = geneList,    # the full ranked list
             OrgDb        = org.Hs.eg.db,
             ont          = "BP",
             pvalueCutoff = 0.05,
             verbose      = FALSE)

head(gse[, c("Description", "NES", "p.adjust")])

# the classic GSEA running-score plot for the top gene set
gseaplot2(gse, geneSetID = 1, title = gse$Description[1])

4Reading the results

Whichever approach you use, the result is a table of biological terms with a few key columns:

• Description: the name of the GO term or pathway.
• p.adjust: the significance, corrected for testing thousands of gene sets. Trust terms below 0.05, just like padj for genes.
• GeneRatio / setSize: how many of your genes are in that term, and how big the term is.
• NES (GSEA only): the direction and strength of the shift.
• core_enrichment (GSEA only): the specific genes driving the result, your shortlist for follow-up.

The goal is not a giant table. It is a sentence: "the genes that changed are enriched for X, which makes sense because our treatment does Y." That sentence is the biological result of the whole analysis.

Using your own DESeq2 results

To run this on your own data, you need two small steps: translate your gene symbols into Entrez IDs, and rank them. Here is the bridge from a DESeq2 results table:

# res is your DESeq2 results data frame, with gene symbols as row names
ids <- bitr(rownames(res), fromType = "SYMBOL",
            toType = "ENTREZID", OrgDb = org.Hs.eg.db)

# build a ranked named vector (rank by the test statistic, high to low)
ranks <- res$stat
names(ranks) <- rownames(res)
ranks <- ranks[ids$SYMBOL]
names(ranks) <- ids$ENTREZID
ranks <- sort(ranks, decreasing = TRUE)
# now feed 'ranks' into gseGO() exactly as above

Check your understanding

Track 3

See the full transcriptomics roadmap →