Module 4: Enrichment analysis
Module 4 - Practical Exercise
In this tutorial you will perform functional enrichment analysis using:
- Over-Representation Analysis (ORA) with gProfiler and clusterProfiler
- Gene Set Enrichment analysis (GSEA) with fgsea and MSigDB
What is enrichment analysis?
Enrichment analysis tests whether specific pathways or biological functions are overrepresented in our DEG list.
Examples:
- immune response
- metabolism
- cell cycle
- apoptosis
Get DESeq2 Differentially Enriched Genes from database
library(DESeq2)
url <- "https://www.ncbi.nlm.nih.gov/geo/download/?type=rnaseq_counts&acc=GSE114360&format=file&file=GSE114360_raw_counts_GRCh38.p13_NCBI.tsv.gz"
raw_counts <- fread(url)
# convert the data.table to matrix format
raw_counts = as.matrix(raw_counts)
class(raw_counts)## [1] "matrix" "array"# set the gene ID values to be the row names for the matrix
rownames(raw_counts) = raw_counts[, "GeneID"]
# now that the gene IDs are the row names, remove the redundant column that contains them
raw_counts = raw_counts[, colnames(raw_counts) != "GeneID"]
# convert the count values from strings (with spaces) to integers, because originally the gene column contained characters, the entire matrix was set to character
class(raw_counts) = "integer"
# view the first few lines of the gene count matrix
#head(raw_counts)
# create a simple one column dataframe to start
metaData <- data.frame("Condition" = c("shRNA", "shRNA", "shRNA", "CTRL", "CTRL", "CTRL"))
# convert the "Condition" column to a factor data type
# the arbitrary order of these factors will determine the direction of log2 fold-changes for the genes (i.e. up or down regulated)
metaData$Condition = factor(metaData$Condition, levels = c("CTRL", "shRNA"))
# set the row names of the metaData dataframe to be the names of our sample replicates from the read counts matrix
rownames(metaData) = colnames(raw_counts)
# view the metadata dataframe
#head(metaData)
# check that names of htseq count columns match the names of the meta data rows
# use the "all" function which tests whether an entire logical vector is TRUE
all(rownames(metaData) == colnames(raw_counts))## [1] TRUEOver-Representation Analysis (ORA)
Enrichment analysis with gProfiler
Documentation : gProfiler
gProfiler is a popular web tool and R package used to identify biological pathways, Gene Ontology terms, and molecular functions enriched in a list of genes.
It compares our differentially expressed genes (DEGs) against known biological annotations from databases such as GO, KEGG, Reactome, and WikiPathways to help interpret the biological meaning of RNA-seq results.
gost.res = gost(
# List of differentially expressed genes
rownames(deg),
# Organism used for the enrichment analysis
# hsapiens = Homo sapiens (human)
organism = "hsapiens",
# Return only statistically significant enriched terms
significant = TRUE,
# Do not test for underrepresented pathways
# FALSE means we only look for overrepresented terms
measure_underrepresentation = FALSE,
# Significance threshold for enrichment analysis
user_threshold = 0.05,
# Multiple testing correction method
# FDR = False Discovery Rate correction
correction_method = "fdr",
# Include evidence codes for GO annotations
evcodes = TRUE
)Filter Enrichment Results
We remove:
- very small terms (term_size)
- extremely large terms (term_size)
- weakly represented terms (intersection_size)
# filter results to reduce to less specialized or general terms with at least 2 genes from our list
gostres_filt = gostres[gostres$term_size > 15 & gostres$term_size < 2000 & gostres$intersection_size > 2,]
#transform score with log10, higher is better
gostres_filt$log10=-1 * log10(gostres_filt$p_value)
#order results by decreasing score value
gostres_filt=gostres_filt[order(gostres_filt$log10,decreasing=T),]Manhattan plot for all the results
Interpretation
Each point represents:
- one biological term/pathway
Higher points indicate:
- stronger enrichment
- lower adjusted p-values
Colors correspond to different databases:
- GO
- KEGG
- Reactome
- WikiPathways
Barplot results for each database
# barplots for top 10 KEGG pathways
res_table = gostres_filt[gostres_filt$source=="KEGG",][1:10, ]
#modify margins for long terms name
par(mar = c(5, 20, 4, 2))
barplot(
res_table$log10,
horiz = TRUE,
names.arg = res_table$term_name,
col = "steelblue",
xlab = "-log10(padj)",
xlim = c(0, ceiling(max(res_table$log10))),
main = "Top 10 enriched terms",
las = 1, # horizontal labels
cex.names = 0.8 # text size
)
Enrichment analysis with clusterProfiler
clusterProfiler is a widely used Bioconductor package for statistical analysis and visualization of functional profiles from genomic data.
It allows users to perform enrichment analysis using Gene Ontology (GO), KEGG, Reactome, and Gene Set Enrichment Analysis (GSEA), while providing powerful visualization tools such as dotplots, enrichment maps, and pathway networks.
- GO Biological Process Enrichment
ego_bp = enrichGO(
# List of differentially expressed genes
gene = rownames(deg),
# Annotation database for human genes
OrgDb = org.Hs.eg.db,
# Type of gene identifiers used in the gene list
keyType = "ENTREZID",
ont = "BP",
#Method used to adjust p-values for multiple testing
#BH = Benjamini-Hochberg False Discovery Rate correction
pAdjustMethod = "BH",
# Maximum raw p-value threshold
pvalueCutoff = 0.05,
# Maximum adjusted p-value threshold (FDR/q-value)
qvalueCutoff = 0.05,
# Convert gene IDs into readable gene symbols
readable = TRUE
)
head(as.data.frame(ego_bp))Question
What does “BP” stand for?
Try changing:
or
What changes?
Visualize GO Results
GO Barplot
GO Dotplot
Interpretation of Dotplots
Dot size represents:
- number of genes
Dot color represents:
- statistical significance
KEGG Pathway Enrichment
ekegg = enrichKEGG(
gene = rownames(deg),
organism = "hsa",
pvalueCutoff = 0.05
)
head(as.data.frame(ekegg))Visualize KEGG Results
- Barplot
- Dotplot
Change the code to include only the top10 KEGG pathways.
- Pathview
The pathview package allows us to visualize gene expression changes directly on KEGG pathway diagrams.
Genes are colored according to their expression levels or fold changes, making it easier to identify activated or repressed biological pathways and interpret RNA-seq results in a biological context.
library(pathview)
# Create named vector of log2 fold changes
gene_list <- deg$log2FoldChange
# Names must be Entrez IDs
names(gene_list) <- rownames(deg)
# Visualize pathway
pathview(
gene.data = names(gene_list),
pathway.id = "hsa05205",
species = "hsa",
out.suffix = "init"
)
What is the difference with this code?
# Visualize pathway
pathview(
gene.data = gene_list,
pathway.id = "hsa05205",
species = "hsa",
out.suffix = "logFC"
)
To customize the map, we need to identify the range of expression of the DEGs involved in this pathway:
library(KEGGREST)
# Download KEGG pathway information
pathway_info = keggGet("hsa05205")
# Extract Entrez gene IDs
pathway_genes = pathway_info[[1]]$GENE
# Keep only Entrez IDs
pathway_entrez = pathway_genes[seq(1, length(pathway_genes), 2)]
head(pathway_entrez)## [1] "10000" "1026" "10451" "10818" "10855" "110117499"## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.4775 -0.2462 -0.2124 -0.1964 -0.1854 0.2699Redraw the map with an appropriate logFC range:
Show solution
# Visualize pathway
pathview(
gene.data = gene_list,
pathway.id = "hsa05205",
species = "hsa",
out.suffix = "logFCx",
limit = list(gene = 0.5) #logFC from -0.5 to 0.5
)
Now the scale color is correct, but from green to red. We usually avoid using these colors in figures. Change the scale to: yellow - blue
Show solution
# Visualize pathway
pathview(
gene.data = gene_list,
pathway.id = "hsa05205",
species = "hsa",
limit = list(gene = 0.5), #logFC from -0.5 to 0.5
out.suffix = "logFCx_YB",
# Low values (downregulated genes)
low = list(gene = "yellow"),
# High values (upregulated genes)
high = list(gene = "blue")
)
Gene Set Enrichment Analysis (GSEA)
Unlike classical enrichment analysis that uses only significant DEGs, GSEA uses ALL genes ranked by a statistic such as:
- log2FoldChange
- adjusted p-value
This approach avoids arbitrary DEG thresholds and can detect subtle but coordinated biological changes.
Prepare Ranked Gene List
We will rank genes using both:
- adjusted p-value (
padj) - direction of change (
log2FoldChange)
Genes with: - strong fold-changes - low adjusted p-values
will obtain the highest scores.
# Remove genes with missing adjusted p-values
res_gsea = as.data.frame(res)
res_gsea = res_gsea[!is.na(res_gsea$padj),]Create Ranking Score
We combine:
- statistical significance
- fold-change direction
using:
gene_rank = -log10(res_gsea$padj) * sign(res_gsea$log2FoldChange)
# Add gene names
names(gene_rank) = rownames(res_gsea)
# Sort decreasingly
gene_rank = sort(gene_rank,decreasing = TRUE)
head(gene_rank)## 9636 8638 3433 29968 16 7803
## 129.14920 61.96901 36.31050 34.86004 27.44421 25.28048Visualize Ranking Distribution
hist(
gene_rank,
breaks = 100,
col = "steelblue",
main = "Distribution of GSEA Ranking Scores",
xlab = "Ranking score"
)
Run GSEA with Gene Sets from MSigDB
MSigDB (Molecular Signatures Database) is a curated collection of annotated gene sets commonly used for pathway and enrichment analysis.
It contains multiple categories of biological signatures, including Hallmark pathways, Gene Ontology terms, immune signatures, and curated pathways from databases such as KEGG and Reactome.
We will use the Hallmark gene sets.
Hallmark pathways summarize major biological processes such as:
- inflammation
- hypoxia
- apoptosis
- MYC signaling
Create Pathway List
fgsea requires a list format:
- pathway name
- vector of genes
## $HALLMARK_ADIPOGENESIS
## [1] "19" "11194" "10449" "33" "34" "35" "47" "50"
## [9] "51" "112" "149685" "9370" "79602" "56894" "9131" "204"
## [17] "217" "226" "284" "51129" "334" "348" "369" "10124"
## [25] "64225" "483" "539" "11176" "593" "23786" "604" "718"
## [33] "847" "284119" "8436" "901" "977" "9936" "948" "1031"
## [41] "400916" "1147" "1149" "134147" "51727" "1306" "1282" "51805"
## [49] "84274" "57017" "1337" "1349" "1351" "1376" "1384" "1431"
## [57] "1537" "1580" "1629" "1652" "1666" "8694" "1717" "51635"
## [65] "25979" "1737" "1738" "4189" "29103" "128338" "1891" "1892"
## [73] "84173" "79071" "5168" "2053" "2101" "23344" "2109" "2167"
## [81] "2184" "8322" "9908" "1647" "2632" "27069" "57678" "137964"
## [89] "2820" "10243" "2878" "2879" "80273" "3033" "26275" "26353"
## [97] "3417" "3419" "3421" "3459" "10989" "3679" "80760" "6453"
## [105] "84522" "3910" "3952" "3977" "3991" "10162" "4023" "4056"
## [113] "8491" "56922" "4191" "4199" "11343" "4259" "84895" "56246"
## [121] "29088" "54996" "23788" "4638" "64859" "4698" "4706" "4713"
## [129] "4722" "28512" "4836" "4958" "5004" "27250" "10400" "5195"
## [137] "5209" "5211" "5236" "23187" "5264" "415116" "123" "5447"
## [145] "5468" "5495" "84919" "10935" "10113" "55037" "5733" "5860"
## [153] "83871" "7905" "92840" "56729" "54884" "8780" "55177" "26994"
## [161] "6239" "10313" "25813" "949" "6342" "6390" "6391" "6573"
## [169] "6510" "6576" "1468" "376497" "8884" "130814" "6623" "6647"
## [177] "10580" "65124" "8404" "58472" "8082" "6776" "2040" "8802"
## [185] "6817" "6888" "10010" "7086" "10140" "7263" "7316" "29979"
## [193] "83549" "7351" "29796" "10975" "7384" "27089" "7423" "7532"
##
## $HALLMARK_ALLOGRAFT_REJECTION
## [1] "16" "6059" "10006" "43" "92" "207" "322" "567"
## [9] "586" "8915" "602" "672" "717" "822" "9607" "6356"
## [17] "6357" "6363" "6347" "6367" "6351" "6352" "6354" "894"
## [25] "896" "1230" "729230" "1234" "912" "914" "919" "940"
## [33] "915" "916" "917" "920" "958" "959" "961" "924"
## [41] "972" "973" "941" "942" "925" "926" "10225" "1029"
## [49] "5199" "56253" "1435" "1445" "1520" "10563" "4283" "2833"
## [57] "1615" "8560" "8444" "1956" "8661" "8664" "8669" "8672"
## [65] "1984" "1991" "2000" "2069" "2113" "2147" "2149" "355"
## [73] "356" "2213" "2268" "2316" "2533" "2589" "2634" "2650"
## [81] "11146" "8477" "3001" "3002" "3059" "9734" "3091" "3105"
## [89] "3108" "3109" "3111" "3112" "3117" "3122" "3133" "3135"
## [97] "3383" "23308" "3455" "3458" "3459" "3460" "10261" "3551"
## [105] "3586" "3589" "3592" "3593" "3594" "3596" "3600" "3603"
## [113] "3606" "8807" "3553" "3558" "9466" "3559" "3560" "3561"
## [121] "3565" "3566" "3569" "3574" "3578" "3624" "3625" "3662"
## [129] "3665" "3394" "3683" "3689" "3702" "3717" "3824" "3848"
## [137] "3932" "3937" "3976" "4050" "4065" "9450" "4067" "6885"
## [145] "11184" "4153" "4318" "11222" "4528" "4689" "4690" "9437"
## [153] "114548" "4830" "4843" "4869" "5196" "5551" "5579" "5582"
## [161] "5699" "5777" "5788" "5917" "8767" "6170" "6123" "6133"
## [169] "6223" "6189" "6203" "27240" "8651" "9655" "6688" "5552"
## [177] "7903" "23166" "6772" "6775" "6890" "6891" "6892" "7040"
## [185] "7042" "7070" "7076" "7096" "7097" "7098" "10333" "7124"
## [193] "7163" "7186" "50852" "7321" "7334" "7453" "7454" "7535"
##
## $HALLMARK_ANDROGEN_RESPONSE
## [1] "10257" "11057" "2181" "87" "9510" "11047" "9590" "207"
## [9] "220" "56172" "10513" "84159" "563" "567" "2683" "658"
## [17] "10645" "595" "896" "8555" "1021" "55839" "1622" "1718"
## [25] "4189" "2005" "22936" "60481" "3992" "2289" "2773" "23171"
## [33] "2936" "2982" "3005" "8916" "3156" "3157" "9455" "3248"
## [41] "51171" "3422" "8821" "3638" "10788" "3685" "3817" "354"
## [49] "3880" "3856" "3977" "3998" "4094" "4117" "9053" "10461"
## [57] "10627" "8031" "10397" "55768" "4824" "5036" "10611" "5238"
## [65] "8554" "8611" "56937" "2185" "11099" "5867" "6197" "23223"
## [73] "6303" "6319" "9871" "6414" "6446" "1836" "54407" "6611"
## [81] "6652" "60559" "25803" "6722" "6728" "79689" "27347" "445347"
## [89] "23585" "7113" "25816" "7163" "8848" "6675" "7329" "51465"
## [97] "9218" "7520" "2547" "65986" "57178"
##
## $HALLMARK_ANGIOGENESIS
## [1] "350" "351" "894" "1281" "1290" "6372" "2260" "11167" "3685"
## [10] "182" "3714" "3764" "4023" "4043" "4060" "4487" "8829" "4973"
## [19] "5154" "5196" "8993" "10631" "5553" "5747" "6275" "5104" "6578"
## [28] "6696" "6781" "7056" "7076" "27242" "7410" "1462" "7422" "7448"
##
## $HALLMARK_APICAL_JUNCTION
## [1] "58" "60" "70" "71" "72" "87" "88" "89"
## [9] "81" "8751" "8745" "8754" "11096" "147" "208" "10000"
## [17] "247" "268" "57463" "347902" "9459" "10109" "478" "2683"
## [25] "10458" "649" "253559" "57863" "794" "10487" "30835" "29126"
## [33] "80381" "947" "942" "4267" "999" "1009" "1013" "1001"
## [41] "1002" "1004" "1006" "1024" "1041" "51148" "5010" "23562"
## [49] "24146" "51208" "149461" "1364" "7122" "9074" "1366" "9073"
## [57] "9080" "1265" "1272" "1307" "1308" "1297" "1384" "92359"
## [65] "1495" "1500" "6376" "8449" "1739" "1758" "1823" "1825"
## [73] "1956" "2037" "51466" "60412" "2200" "2318" "6624" "2533"
## [81] "2593" "2770" "2771" "2886" "2962" "3033" "3265" "3383"
## [89] "3384" "3386" "7087" "8517" "3636" "3638" "3667" "8515"
## [97] "3673" "3675" "3680" "3688" "3691" "83700" "3728" "3757"
## [105] "3881" "3909" "3914" "3918" "143903" "26119" "51474" "8174"
## [113] "51776" "5871" "5600" "5603" "1432" "4192" "4313" "4318"
## [121] "9019" "10205" "4478" "4597" "4628" "4627" "103910" "10398"
## [129] "5818" "5819" "25945" "81607" "257194" "91624" "4763" "4771"
## [137] "23114" "79849" "57555" "54413" "4892" "4902" "9379" "64398"
## [145] "84552" "55742" "5089" "5097" "5175" "5216" "5291" "8503"
## [153] "5310" "5335" "5522" "5728" "5747" "5788" "5880" "5921"
## [161] "54509" "6237" "6251" "9672" "8910" "6464" "357" "140885"
## [169] "8935" "7781" "9353" "10174" "10290" "6714" "6810" "6850"
## [177] "8189" "9344" "7045" "7059" "7070" "7073" "7082" "51754"
## [185] "4982" "7185" "7216" "7248" "7106" "7283" "7408" "7410"
## [193] "7412" "1462" "7414" "7450" "8976" "65266" "7533" "7791"
##
## $HALLMARK_APICAL_SURFACE
## [1] "102" "79602" "84632" "9465" "351" "50617" "5205" "2683"
## [9] "672" "11126" "9696" "202" "6376" "131566" "1946" "2050"
## [17] "2319" "2619" "2625" "51738" "2947" "3315" "3560" "3561"
## [25] "4067" "27076" "4118" "266727" "23054" "22854" "255738" "5314"
## [33] "5329" "51458" "146760" "80274" "357" "116085" "6517" "142680"
## [41] "8406" "55959" "7070" "51754"Run fgsea
fgsea_res = fgsea(
# Pathway list
pathways = hallmark_list,
# Ranked gene list
stats = gene_rank,
# Minimum pathway size
minSize = 15,
# Maximum pathway size
maxSize = 500,
# Number of permutations
nperm = 1000
)Sort and Explore Results
Question
What does the NES (Normalized Enrichment Score) represent?
Interpretation
Positive NES: pathway enriched among upregulated genes
Negative NES: pathway enriched among downregulated genes
Large absolute NES: stronger enrichment signal
Visualize Top Pathways
topPathwaysUp <- fgsea_res %>%
filter(NES > 0) %>%
arrange(desc(NES)) %>%
head(10)
ggplot(
topPathwaysUp,
aes(x = reorder(pathway, NES),
y = NES)
) +
geom_col(fill = "steelblue") +
coord_flip() +
theme_bw() +
xlab("Pathway") +
ylab("Normalized Enrichment Score") +
ggtitle("Top Positively Enriched Pathways")
Visualize Downregulated Pathways
topPathwaysDown <- fgsea_res %>%
filter(NES < 0) %>%
arrange(NES) %>%
head(10)
ggplot(
topPathwaysDown,
aes(x = reorder(pathway, NES),
y = NES)
) +
geom_col(fill = "tomato") +
coord_flip() +
theme_bw() +
xlab("Pathway") +
ylab("Normalized Enrichment Score") +
ggtitle("Top Negatively Enriched Pathways")
Enrichment Plot for One Pathway
plotEnrichment(
hallmark_list[["HALLMARK_INTERFERON_GAMMA_RESPONSE"]],
gene_rank
) +
labs(
title = "GSEA Enrichment Plot: Apoptosis"
)
Interpretation of Enrichment Plot
The curve shows where pathway genes appear in the ranked gene list.
Genes concentrated at the top: pathway activated
Genes concentrated at the bottom: pathway repressed
Visualize Multiple Pathways
top_pathways <- fgsea_res$pathway[1:5]
plotGseaTable(
pathways = hallmark_list[top_pathways],
stats = gene_rank,
fgseaRes = fgsea_res,
gseaParam = 1
)
Volcano-style NES Plot
ggplot(
fgsea_res,
aes(x = NES,
y = -log10(padj))
) +
geom_point(aes(color = padj < 0.05),
size = 3) +
theme_bw() +
scale_color_manual(
values = c("grey70", "red")
) +
xlab("Normalized Enrichment Score") +
ylab("-log10 adjusted p-value") +
ggtitle("fgsea Pathway Enrichment")
Biological Interpretation Questions
- Which pathways are activated?
- Which pathways are repressed?
- Are immune pathways enriched?
- Are proliferation pathways activated?
- Which biological processes dominate the dataset?
Show answer
#Which pathways are activated?
activated <- fgsea_res[fgsea_res$NES > 0 & fgsea_res$padj < 0.05, ]
dim(activated)## [1] 0 8## Empty data.table (0 rows and 8 cols): pathway,pval,padj,ES,NES,nMoreExtreme...#Which pathways are repressed?
repressed <- fgsea_res[fgsea_res$NES < 0 & fgsea_res$padj < 0.05, ]
dim(repressed)## [1] 2 8## pathway pval padj ES
## <char> <num> <num> <num>
## 1: HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION 0.001251564 0.03101266 -0.8540225
## 2: HALLMARK_MYOGENESIS 0.001265823 0.03101266 -0.7570904
## NES nMoreExtreme size leadingEdge
## <num> <num> <int> <list>
## 1: -1.845140 0 120 960,5054,11167,3673,56937,5768,...[44]
## 2: -1.608063 0 97 351,3693,3691,1266,6319,1191,...[21]#Are immune pathways enriched? Yes
immune <- fgsea_res[grep("immune|interferon|cytokine|inflammatory", fgsea_res$pathway, ignore.case = TRUE), ]
immune_sig <- immune[immune$padj < 0.05, ]
# Are proliferation pathways activated? No
prolif <- fgsea_res[grep("cell cycle|mitosis|E2F|G2M|proliferation", fgsea_res$pathway, ignore.case = TRUE), ]
prolif_sig <- prolif[prolif$padj < 0.05, ]
prolif_sig## Empty data.table (0 rows and 8 cols): pathway,pval,padj,ES,NES,nMoreExtreme...# Which biological processes dominate the dataset?
top <- fgsea_res[fgsea_res$padj < 0.05, ]
top <- top[order(-abs(top$NES)), ]
head(top, 20)## pathway pval padj ES
## <char> <num> <num> <num>
## 1: HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION 0.001251564 0.03101266 -0.8540225
## 2: HALLMARK_MYOGENESIS 0.001265823 0.03101266 -0.7570904
## NES nMoreExtreme size leadingEdge
## <num> <num> <int> <list>
## 1: -1.845140 0 120 960,5054,11167,3673,56937,5768,...[44]
## 2: -1.608063 0 97 351,3693,3691,1266,6319,1191,...[21]