WM7 Functional Annotation of Differential Expression Results

ORA, GSEA, and sample-level pathway scoring

library(dplyr)
library(readr)
library(tibble)
library(ggplot2)
library(writexl)
library(msigdbr)
library(fgsea)
library(BiocParallel)
library(DESeq2)
library(GSVA)
library(ComplexHeatmap)
library(circlize)
library(EnhancedVolcano)

register(SerialParam())

Aim

In this practical we start from the DESeq2 result table and ask which biological processes and pathways are affected by MEK inhibitor treatment in Caco-2 cells.

• manual differential gene lists for the Enrichr web interface
• one explicit Fisher exact test
• over-representation analysis (ORA) with msigdbr gene sets
• a simple pre-ranked GSEA-style analysis without a hard differential gene cutoff
• sample-level Hallmark scoring with GSVA
• extension from Hallmark gene sets to Gene Ontology Biological Process

1. Load DESeq2 Results

res <- read_tsv("DESeq_result_full.tsv", show_col_types = FALSE)

res

The important columns are:

• gene_name: gene name
• log2FoldChange: positive values are higher in MEKi, negative values are lower
• stat: Wald test statistic, useful for ranking genes
• padj: FDR-adjusted p-value

Before enrichment, remove missing or duplicated gene names. For duplicated genes, keep the row with the largest absolute test statistic.

res_clean <- res |>
  filter(!is.na(gene_name), gene_name != "", !is.na(stat)) |>
  group_by(gene_name) |>
  slice_max(order_by = abs(stat), n = 1, with_ties = FALSE) |>
  ungroup()

nrow(res)

[1] 17666

nrow(res_clean)

[1] 17666

Define the tested gene universe. This is the background set for ORA: all genes that had enough expression to enter the DESeq2 analysis and have a usable gene name.

gene_universe <- res_clean |>
  pull(gene_name) |>
  unique()

length(gene_universe)

[1] 17666

2. Manual Differential Gene Lists for Enrichr

Web tools such as Enrichr are useful for a first broad look at the biology. For this exploratory step, select a clear set of differential genes.

de_genes <- res_clean |>
  filter(!is.na(padj), padj < 0.1, abs(log2FoldChange) >= 1)

de_up <- de_genes |>
  filter(log2FoldChange > 0) |>
  arrange(padj)

de_down <- de_genes |>
  filter(log2FoldChange < 0) |>
  arrange(padj)

nrow(de_genes)

[1] 625

nrow(de_up)

[1] 432

nrow(de_down)

[1] 193

Write one-gene-per-line lists for upload to Enrichr:

write_lines(de_genes$gene_name, "enrichr_DE_genes_padj0.1_abslog2FC1.txt")
write_lines(de_up$gene_name, "enrichr_up_in_MEKi_padj0.1_abslog2FC1.txt")
write_lines(de_down$gene_name, "enrichr_down_in_MEKi_padj0.1_abslog2FC1.txt")

Task: Upload the three files to Enrichr and compare the broad biological themes. Do the up- and down-regulated lists tell the same story?

Check Different Cutoffs

The choice of cutoff changes which genes enter the list. Before interpreting results, check whether the strongest Enrichr terms are stable across a few reasonable cutoffs.

# Practice:
# Create six gene tables from res_clean:
#   padj < 0.05, up and down
#   padj < 0.05, up only
#   padj < 0.05, down only
#   padj < 0.001, up and down
#   padj < 0.001, up only
#   padj < 0.001, down only
#
# Use filter(), !is.na(padj), padj, and log2FoldChange.
# Then count the rows in each table with nrow().

# Practice:
# Use write_lines() to save one gene list for each cutoff table.
# Each file should contain one gene name per line.

Task: Upload these six lists to Enrichr. Which terms are present no matter which cutoff you use? Which terms only appear with the larger or smaller gene lists? What changes when up- and down-regulated genes are tested separately?

3. Load MSigDB Gene Sets with msigdbr

msigdbr provides MSigDB gene sets as simple R tables.

Start with Hallmark gene sets. These are compact, curated, and easier to interpret than the full Gene Ontology.

hallmark <- msigdbr(species = "Homo sapiens", collection = "H") |>
  select(gs_name, gene = 1) |>
  distinct()

hallmark

4. One Explicit Fisher Exact Test

ORA asks whether a pathway contains more differential genes than expected by chance. We can test this directly for one gene set.

Here we test HALLMARK_KRAS_SIGNALING_UP, a useful example for MEK inhibition because it is connected to the RAS/RAF/MEK/ERK axis.

example_set_name <- "HALLMARK_KRAS_SIGNALING_UP"

example_set <- hallmark |>
  filter(gs_name == example_set_name) |>
  pull(gene) |>
  intersect(gene_universe)

selected_genes <- de_genes |>
  pull(gene_name) |>
  unique()

# Practice:
# Count the four groups needed for the Fisher exact test.
#
# Use ?intersect, ?union, and ?setdiff to look up what these functions do.
#
# de_and_in_set: genes that are selected and are in the example gene set
# de_not_in_set: genes that are selected but are not in the example gene set
# not_de_in_set: genes that are not selected but are in the example gene set
# not_de_not_in_set: genes that are not selected and are not in the example gene set

fisher_table <- matrix(
  c(de_and_in_set, de_not_in_set, not_de_in_set, not_de_not_in_set),
  nrow = 2,
  byrow = TRUE,
  dimnames = list(
    differential = c("yes", "no"),
    in_gene_set = c("yes", "no")
  )
)

fisher_table

            in_gene_set
differential yes    no
         yes  11   614
         no  127 16914

fisher.test(fisher_table, alternative = "greater")

    Fisher's Exact Test for Count Data

data:  fisher_table
p-value = 0.009961
alternative hypothesis: true odds ratio is greater than 1
95 percent confidence interval:
 1.302169      Inf
sample estimates:
odds ratio 
  2.385792 

Task: Interpret the four numbers in the contingency table. Which cell contains differential genes that are members of the pathway?

5. ORA with msigdbr

Now run the same idea across all Hallmark gene sets.

hallmark_pathways <- split(hallmark$gene, hallmark$gs_name)

hallmark_ora_tab <- fora(
  pathways = hallmark_pathways,
  genes = selected_genes,
  universe = gene_universe,
  minSize = 10,
  maxSize = 500
) |>
  as_tibble() |>
  mutate(overlapGenes = vapply(overlapGenes, paste, character(1), collapse = "/")) |>
  arrange(padj)

hallmark_ora_tab

hallmark_ora_tab |>
  slice_min(padj, n = 15) |>
  ggplot(aes(x = -log10(padj), y = reorder(pathway, -log10(padj)))) +
  geom_col() +
  labs(
    x = "-log10(adjusted p-value)",
    y = NULL,
    title = "Hallmark ORA from thresholded DE genes"
  ) +
  theme_minimal()

Check whether the ORA result changes when using the same cutoffs as in the Enrichr exercise.

# Practice:
# Repeat the Hallmark ORA with the same six gene lists from the Enrichr task:
#   padj < 0.05, up and down
#   padj < 0.05, up only
#   padj < 0.05, down only
#   padj < 0.001, up and down
#   padj < 0.001, up only
#   padj < 0.001, down only
#
# For each list:
#   pull(gene_name)
#   run fora()
#   arrange by padj
#   compare the top Hallmark terms

Save the ORA table:

write_tsv(hallmark_ora_tab, "ORA_Hallmark_msigdb.tsv")
write_xlsx(hallmark_ora_tab, "ORA_Hallmark_msigdb.xlsx")

Discussion: ORA depends on the cutoff used to select genes. Which Hallmark terms are stable across the six lists? Which terms depend on testing up- and down-regulated genes separately?

6. GSEA-Style Analysis Without an Explicit Cutoff

GSEA-style analysis avoids selecting a differential-gene cutoff. Instead, it uses all genes ranked from most up-regulated in MEKi to most down-regulated in MEKi.

Here we rank genes by the DESeq2 Wald statistic.

ranked_gene_list <- res_clean |>
  arrange(desc(stat)) |>
  select(gene_name, stat)

ranks <- ranked_gene_list$stat
names(ranks) <- ranked_gene_list$gene_name

head(ranks)

   LAMA1     CYMP    PHGR1    RGPD6     KNG1  SELENOP 
17.22528 16.29956 15.93474 15.35809 15.26151 14.99842 

tail(ranks)

      ID1     CLDN2     DUSP6   CEACAM6      EGR1    ANKRD1 
-17.34017 -18.34090 -18.56878 -19.19697 -25.28795 -27.30598 

set.seed(1)

hallmark_gsea_tab <- fgsea(
  pathways = hallmark_pathways,
  stats = ranks,
  minSize = 10,
  maxSize = 500,
  eps = 0
) |>
  as_tibble() |>
  mutate(leadingEdge = vapply(leadingEdge, paste, character(1), collapse = "/")) |>
  arrange(padj)

Warning in prepareStats(stats, scoreType, gseaParam): There are ties in the preranked stats (0.03% of the list).
The order of those tied genes will be arbitrary, which may produce unexpected results.

hallmark_gsea_tab

hallmark_gsea_tab |>
  slice_min(padj, n = 15) |>
  ggplot(aes(x = NES, y = reorder(pathway, NES), fill = NES > 0)) +
  geom_col() +
  scale_fill_manual(values = c("TRUE" = "firebrick", "FALSE" = "steelblue"), guide = "none") +
  labs(
    x = "Normalized enrichment score",
    y = NULL,
    title = "Hallmark GSEA from all ranked genes"
  ) +
  theme_minimal()

Positive NES means enrichment near the top of the ranked list, i.e. genes tend to be higher in MEKi. Negative NES means enrichment near the bottom, i.e. genes tend to be lower in MEKi.

write_tsv(hallmark_gsea_tab, "GSEA_Hallmark_msigdb.tsv")
write_xlsx(hallmark_gsea_tab, "GSEA_Hallmark_msigdb.xlsx")

Make a volcano-style plot of the Hallmark GSEA result.

# Practice:
# Use EnhancedVolcano() with hallmark_gsea_tab.
#
# Put NES on the x-axis.
# Put padj on the y-axis.
# Use pathway names as labels.
# A useful cutoff for NES is 1.

Task: Compare ORA and GSEA results. Which biological themes appear in both? Which only appear when all genes are ranked?

7. Sample-Level Hallmark Scores with GSVA

ORA and GSEA use the differential expression result table. GSVA answers a different question: how active does each Hallmark gene set look in each sample?

For this we use variance-stabilized expression values, as in the PCA workflow.

dds <- readRDS("dds.rds")

sample_info <- colData(dds) |>
  as.data.frame() |>
  select(sample, condition)

sample_info

Run the same variance-stabilizing transformation used for PCA, then prepare an expression matrix with gene names.

vst_counts <- vst(dds, blind = FALSE)

expression_matrix <- assay(vst_counts) |>
  as.data.frame() |>
  rownames_to_column("gene_id") |>
  inner_join(res_clean |> select(gene_id, gene_name), by = "gene_id") |>
  select(gene_name, all_of(sample_info$sample)) |>
  group_by(gene_name) |>
  summarize(across(all_of(sample_info$sample), mean), .groups = "drop") |>
  column_to_rownames("gene_name") |>
  as.matrix()

dim(expression_matrix)

[1] 16634     6

Now calculate one GSVA score per Hallmark gene set and sample.

gsva_parameters <- gsvaParam(
  expression_matrix,
  hallmark_pathways,
  minSize = 10,
  maxSize = 500,
  kcdf = "Gaussian"
)

hallmark_gsva <- gsva(gsva_parameters)

ℹ GSVA version 2.6.2

ℹ Searching for rows with constant values

ℹ Calculating GSVA ranks

ℹ GSVA dense (classical) algorithm

ℹ Row-wise ECDF estimation with Gaussian kernels

ℹ Calculating row ECDFs

ℹ Calculating column ranks

ℹ GSVA dense (classical) algorithm

ℹ Calculating GSVA scores for 50 gene sets

✔ Calculations finished

dim(hallmark_gsva)

[1] 50  6

hallmark_gsva[1:5, ]

                                 CaCo2_1     CaCo2_2      CaCo2_3     CaCo2_4
HALLMARK_ADIPOGENESIS        -0.09536108  0.04799363 -0.004000484  0.02150758
HALLMARK_ALLOGRAFT_REJECTION  0.05645103  0.21113752 -0.064706956 -0.19217295
HALLMARK_ANDROGEN_RESPONSE   -0.24280486  0.03984928  0.008465147  0.19013019
HALLMARK_ANGIOGENESIS        -0.16997774 -0.07992886 -0.246971939 -0.04243008
HALLMARK_APICAL_JUNCTION     -0.11032190 -0.06794653 -0.172646523 -0.13008694
                                 CaCo2_5     CaCo2_6
HALLMARK_ADIPOGENESIS        -0.14806896  0.11667805
HALLMARK_ALLOGRAFT_REJECTION  0.02767541 -0.09323591
HALLMARK_ANDROGEN_RESPONSE   -0.03611282  0.02739779
HALLMARK_ANGIOGENESIS         0.25628314  0.24506663
HALLMARK_APICAL_JUNCTION      0.30505634  0.16541892

For the heatmap, we scale each Hallmark row. Red means relatively high activity for that Hallmark in that sample; blue means relatively low activity.

hallmark_gsva_scaled <- t(scale(t(hallmark_gsva)))

condition_colors <- c(control = "grey70", MEKi = "firebrick")

sample_annotation <- HeatmapAnnotation(
  condition = sample_info$condition,
  col = list(condition = condition_colors)
)

Heatmap(
  hallmark_gsva_scaled,
  name = "GSVA\nz-score",
  col = colorRamp2(c(-2, 0, 2), c("steelblue", "white", "firebrick")),
  top_annotation = sample_annotation,
  column_split = sample_info$condition,
  cluster_columns = FALSE,
  row_names_gp = grid::gpar(fontsize = 7),
  column_names_gp = grid::gpar(fontsize = 9)
)

Save the GSVA scores:

hallmark_gsva_tab <- hallmark_gsva |>
  as.data.frame() |>
  rownames_to_column("hallmark")

write_tsv(hallmark_gsva_tab, "GSVA_Hallmark_sample_scores.tsv")
write_xlsx(hallmark_gsva_tab, "GSVA_Hallmark_sample_scores.xlsx")

Task: Compare the GSVA heatmap to the Hallmark GSEA result. Do the pathways that change globally also separate the individual control and MEKi samples?

8. Extend to Gene Ontology Biological Process

GO Biological Process is much larger and more redundant than Hallmark. It is useful for depth, but interpretation should focus on themes rather than single near-duplicate term names.

go_bp <- msigdbr(
  species = "Homo sapiens",
  collection = "C5",
  subcollection = "GO:BP"
) |>
  select(gs_name, gene = 1) |>
  distinct()

go_bp

Run ORA with GO BP:

# Practice:
# Create go_bp_pathways with split(), as for the Hallmark example.
#
# Then create go_bp_ora_tab with fora().
# Use the same genes, universe, minSize, and maxSize as before.
# Convert the result to a tibble, simplify overlapGenes, and arrange by padj.

go_bp_ora_tab |>
  slice_min(padj, n = 20) |>
  ggplot(aes(x = -log10(padj), y = reorder(pathway, -log10(padj)))) +
  geom_col() +
  labs(
    x = "-log10(adjusted p-value)",
    y = NULL,
    title = "GO BP ORA from thresholded DE genes"
  ) +
  theme_minimal()

Run GSEA-style analysis with GO BP:

# Practice:
# Create go_bp_gsea_tab with fgsea().
# Use go_bp_pathways and the same ranked list used for Hallmark GSEA.
# Convert the result to a tibble, simplify leadingEdge, and arrange by padj.

go_bp_gsea_plot_tab <- bind_rows(
  go_bp_gsea_tab |>
    slice_max(NES, n = 15),
  go_bp_gsea_tab |>
    slice_min(NES, n = 15)
) |>
  distinct(pathway, .keep_all = TRUE)

go_bp_gsea_plot_tab |>
  ggplot(aes(x = NES, y = reorder(pathway, NES), fill = NES > 0)) +
  geom_col() +
  scale_fill_manual(values = c("TRUE" = "firebrick", "FALSE" = "steelblue"), guide = "none") +
  labs(
    x = "Normalized enrichment score",
    y = NULL,
    title = "GO BP GSEA from all ranked genes"
  ) +
  theme_minimal()

Make a volcano-style plot of the GO BP GSEA result, then compare it with the Hallmark GSEA volcano plot.

# Practice:
# Use EnhancedVolcano() with go_bp_gsea_tab.
#
# Put NES on the x-axis.
# Put padj on the y-axis.
# Use pathway names as labels.
# Use the same cutoffs as for the Hallmark plot.
#
# Compare the two volcano plots:
# Which broad Hallmark patterns are supported by many GO BP terms?
# Which GO BP terms add more detailed biology?
# Which GO BP terms look redundant with each other?

9. Final GSVA Task With GO BP Terms

As a final task, select the top 20 GO BP categories with positive GSEA scores and the top 20 GO BP categories with negative GSEA scores. Then calculate sample-level GSVA scores for these GO BP categories and make a heatmap.

# Practice:
# Start from go_bp_gsea_tab.
# Select the 20 categories with the largest NES.
# Select the 20 categories with the smallest NES.
# Combine these categories into one table.
#
# Use the selected pathway names to subset go_bp_pathways.
# Run GSVA with expression_matrix and the selected GO BP pathways.
# Make a Heatmap with samples split by condition.

One possible result is shown below.

Save the GO BP tables:

write_tsv(go_bp_ora_tab, "ORA_GO_BP_msigdb.tsv")
write_xlsx(go_bp_ora_tab, "ORA_GO_BP_msigdb.xlsx")

write_tsv(go_bp_gsea_tab, "GSEA_GO_BP_msigdb.tsv")
write_xlsx(go_bp_gsea_tab, "GSEA_GO_BP_msigdb.xlsx")

go_bp_gsva_tab <- go_bp_gsva |>
  as.data.frame() |>
  rownames_to_column("go_bp")

write_tsv(go_bp_gsva_tab, "GSVA_GO_BP_top_GSEA_sample_scores.tsv")
write_xlsx(go_bp_gsva_tab, "GSVA_GO_BP_top_GSEA_sample_scores.xlsx")

10. Final Interpretation

Use the result tables and plots to write a short interpretation.

Answer these questions:

1. Which broad Hallmark pathways change after MEK inhibitor treatment?
2. Which GO Biological Process themes support that interpretation?
3. Are the strongest changes mostly higher or lower in MEKi?
4. Do the GSVA sample-level scores separate control and MEKi samples?
5. Do the selected GO BP GSVA scores show the same direction as GO BP GSEA?
6. Which conclusions are robust across Enrichr, ORA, GSEA, and GSVA?
7. Which conclusions depend strongly on the differential-gene cutoff?