BioUnfold.in
RNA-Seq Data Analysis

Lesson 12 of 13 · 12 min

Functional and GO Enrichment

A ranked DE gene list tells you what changed, but not why it matters. Over-representation analysis (ORA) asks whether a biological theme appears among your significant genes more often than expected by chance. That expectation is set by a background gene universe, so choosing the universe correctly matters as much as the DE list itself.

Convert to Entrez IDs

clusterProfiler links genes to GO terms and KEGG pathways through Entrez gene IDs, so symbols or Ensembl IDs must be converted first. The bitr helper maps IDs against the org.Hs.eg.db annotation package; it drops anything it cannot map and warns you what fraction of the input failed. Convert your background the same way, since both lists must live in the same ID space.

r
library(clusterProfiler)
library(org.Hs.eg.db)

# de_genes  : SYMBOLs of significant genes (e.g. padj < 0.05 & |log2FC| > 1)
# all_genes : every gene that passed DESeq2 independent filtering (the universe)

de_map <- bitr(de_genes,  fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)
bg_map <- bitr(all_genes, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)

de_entrez <- unique(de_map$ENTREZID)
bg_entrez <- unique(bg_map$ENTREZID)

GO and KEGG enrichment

r
ego <- enrichGO(gene          = de_entrez,
                universe      = bg_entrez,
                OrgDb         = org.Hs.eg.db,
                keyType       = "ENTREZID",
                ont           = "BP",          # biological process
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.05,
                qvalueCutoff  = 0.10,
                minGSSize     = 10,            # drop tiny gene sets
                maxGSSize     = 500,           # drop over-generic sets
                readable      = TRUE)          # map Entrez back to symbols

ekegg <- enrichKEGG(gene         = de_entrez,
                    universe      = bg_entrez,
                    organism      = "hsa",
                    keyType       = "kegg",    # KEGG hsa IDs == Entrez IDs
                    pAdjustMethod = "BH",
                    pvalueCutoff  = 0.05,
                    minGSSize     = 10,
                    maxGSSize     = 500)
ekegg <- setReadable(ekegg, OrgDb = org.Hs.eg.db, keyType = "ENTREZID")
> head(as.data.frame(ego), 5)[, c("ID", "Description", "GeneRatio", "p.adjust", "Count")]
                   ID                            Description GeneRatio p.adjust Count
GO:0006954 GO:0006954                 inflammatory response    45/612 3.10e-12    45
GO:0002250 GO:0002250              adaptive immune response    38/612 8.24e-10    38
GO:0009611 GO:0009611                  response to wounding    41/612 2.71e-08    41
GO:0030198 GO:0030198       extracellular matrix organization    29/612 6.05e-07    29
GO:0071345 GO:0071345  cellular response to cytokine stimulus    52/612 1.18e-06    52

clusterProfiler adjusts each term's raw hypergeometric p-value for multiple testing with Benjamini-Hochberg, reported as p.adjust. The minGSSize and maxGSSize filters remove tiny sets that are statistically unstable and huge ones too generic to interpret. GeneRatio (k/n) versus BgRatio tells you the fold enrichment behind each adjusted p-value.

Warning: a proper background universe is every gene that had a real chance of being called DE, typically the genes that passed DESeq2 independent filtering, not the whole genome. Omitting universe makes all annotated genes the background, which overstates significance and produces misleadingly long term lists.

GSEA threshold-free alternative

r
# Rank ALL tested genes by a signed statistic (log2FC or the Wald stat log2FC/lfcSE).
# res must already carry an Entrez column (mapped with bitr/mapIds beforehand).
# Names must be Entrez IDs; the vector must be sorted decreasing with unique names.
gene_list <- res$log2FoldChange
names(gene_list) <- res$entrez
gene_list <- gene_list[!is.na(names(gene_list)) & !duplicated(names(gene_list))]
gene_list <- sort(gene_list, decreasing = TRUE)

gse <- gseGO(geneList      = gene_list,
             OrgDb         = org.Hs.eg.db,
             keyType       = "ENTREZID",
             ont           = "BP",
             minGSSize     = 10,
             maxGSSize     = 500,
             pAdjustMethod = "BH",
             pvalueCutoff  = 0.05,
             eps           = 0)             # allow very small p-values (fgsea)

# Unlike ORA, GSEA needs no significance cutoff: it uses the full ranking and
# reports a signed NES, catching coordinated shifts ORA would miss.

Try it yourself: from your own DESeq2 results build de_entrez and bg_entrez as above, run enrichGO(ont = "BP"), then summarize the top themes with dotplot(ego, showCategory = 15). Re-run once with universe = bg_entrez and once with the universe omitted, and note how many more terms clear p.adjust < 0.05 with the wrong background.

Visualizing Results: Volcano and MA PlotsWrap-up and Where to Go Next