WM7 Regulator Activity

Pathway and transcription factor activity with decoupleR

library(DESeq2)
library(dplyr)
library(readr)
library(tibble)
library(tidyr)
library(ggplot2)
library(writexl)
library(decoupleR)
library(ComplexHeatmap)
library(circlize)

Aim

In Day 1 we asked which gene sets were enriched among the differential expression results. Today we ask a different question: which upstream pathways and transcription factors look active?

• infer pathway activity with PROGENy and decoupleR
• infer transcription factor activity with DoRothEA and decoupleR
• compare contrast-level activity scores with sample-level activity scores
• connect activity scores back to the ORA, GSEA, and GSVA results from Day 1

1. Load Differential Expression Results

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

res

Remove rows that cannot be used for activity inference. If a gene appears more than once, keep the row with the strongest DESeq2 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

2. Make A Contrast-Level Input Matrix

For decoupleR, rows are genes and columns are samples or contrasts. First use the DESeq2 Wald statistic as a one-column contrast matrix.

# Practice:
# Start from res_clean.
# Keep gene_name and stat.
# Move gene_name into the row names with column_to_rownames().
# Convert the result to a matrix with as.matrix().
# Give the column the name MEKi_vs_control.

dim(contrast_matrix)

[1] 17666     1

head(contrast_matrix)

          MEKi_vs_control
5_8S_rRNA      1.36968445
7SK            0.08579961
A1BG           0.43044299
A1BG-AS1       0.09589767
A1CF           0.97165659
A2M            0.70371021

Positive values mean higher expression in MEKi. Negative values mean lower expression in MEKi.

3. Make A Sample-Level Input Matrix

For sample-level activity scores, use variance-stabilized counts, as in the PCA and GSVA workflows.

dds <- readRDS("dds.rds")

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

sample_info

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

4. Load Activity Models

PROGENy connects pathways to responsive genes. DoRothEA connects transcription factors to regulated target genes.

progeny_model <- read_tsv("progeny_human_top500.tsv", show_col_types = FALSE)
dorothea_network <- read_tsv("dorothea_human_ABC.tsv", show_col_types = FALSE)

progeny_model

dorothea_network

Check how many targets each pathway or transcription factor has.

progeny_model |>
  count(source, sort = TRUE)

dorothea_network |>
  count(source, sort = TRUE)

5. PROGENy On The DE Contrast

Run decoupleR on the contrast matrix. The score is an inferred pathway activity score for the MEKi versus control comparison.

# Practice:
# Use run_mlm().
# mat should be contrast_matrix.
# network should be progeny_model.
# Use weight as the direction and strength column.
# Arrange the result by score.

progeny_contrast

progeny_contrast |>
  ggplot(aes(x = score, y = reorder(source, score), fill = score > 0)) +
  geom_col() +
  scale_fill_manual(values = c("TRUE" = "firebrick", "FALSE" = "steelblue"), guide = "none") +
  labs(
    x = "PROGENy activity score",
    y = NULL,
    title = "Pathway activity from the DESeq2 statistic"
  ) +
  theme_minimal()

Task: Which pathways are lower in MEKi? Which pathways are higher? Does the MAPK result match the treatment?

6. PROGENy At Sample Level

Run the same model on the variance-stabilized expression matrix. This gives one pathway activity score per sample.

# Practice:
# Use run_mlm() again.
# This time mat should be expression_matrix.
# Convert the long result into a matrix with pivot_wider().
# Rows should be pathways and columns should be samples.

progeny_sample

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

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

Heatmap(
  progeny_sample_matrix,
  name = "activity",
  col = colorRamp2(c(-4, 0, 4), c("steelblue", "white", "firebrick")),
  top_annotation = sample_annotation,
  column_split = sample_info$condition,
  cluster_columns = FALSE,
  row_names_gp = grid::gpar(fontsize = 9),
  column_names_gp = grid::gpar(fontsize = 9)
)

Task: Which pathways separate control and MEKi samples? Are these the same pathways that were strongest in the contrast-level plot?

7. DoRothEA On The DE Contrast

Now infer transcription factor activity from the same DESeq2 statistic matrix.

# Practice:
# Use run_ulm().
# mat should be contrast_matrix.
# network should be dorothea_network.
# Arrange the result by the strongest absolute scores.

tf_contrast

tf_contrast_top |>
  ggplot(aes(x = score, y = reorder(source, score), fill = score > 0)) +
  geom_col() +
  scale_fill_manual(values = c("TRUE" = "firebrick", "FALSE" = "steelblue"), guide = "none") +
  labs(
    x = "TF activity score",
    y = NULL,
    title = "Top transcription factor activities from the DESeq2 statistic"
  ) +
  theme_minimal()

Positive scores mean that the target-gene pattern is higher in MEKi. Negative scores mean that the target-gene pattern is lower in MEKi.

8. Inspect One Transcription Factor

Pick the strongest transcription factor and look at its target genes in the differential expression table.

# Practice:
# Choose one transcription factor from tf_contrast.
# Filter dorothea_network to this transcription factor.
# Join the targets to res_clean.
# Arrange by padj.

selected_tf

[1] "STAT2"

STAT2

selected_tf_targets

Task: Are the target genes mostly moving in the same direction as the transcription factor activity score suggests?

9. DoRothEA At Sample Level

Run the transcription factor model on the variance-stabilized expression matrix. For the heatmap, show the transcription factors with the strongest contrast-level scores.

# Practice:
# Use run_ulm() again.
# This time mat should be expression_matrix.
# Convert the long result into a matrix with pivot_wider().
# Keep the strongest transcription factors from tf_contrast for the heatmap.

tf_sample

Heatmap(
  tf_heatmap_scaled,
  name = "TF\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 = 8),
  column_names_gp = grid::gpar(fontsize = 9)
)

Task: Which transcription factors separate control and MEKi samples? Which ones are also strong in the contrast-level result?

10. Save Result Tables

write_tsv(progeny_contrast, "PROGENy_contrast_activity.tsv")
write_xlsx(progeny_contrast, "PROGENy_contrast_activity.xlsx")

write_tsv(progeny_sample, "PROGENy_sample_activity.tsv")
write_xlsx(progeny_sample, "PROGENy_sample_activity.xlsx")

write_tsv(tf_contrast, "DoRothEA_contrast_activity.tsv")
write_xlsx(tf_contrast, "DoRothEA_contrast_activity.xlsx")

write_tsv(tf_sample, "DoRothEA_sample_activity.tsv")
write_xlsx(tf_sample, "DoRothEA_sample_activity.xlsx")

11. Final Interpretation

Use the pathway and transcription factor activity results together with the Day 1 ORA, GSEA, and GSVA results.

Answer these questions:

1. Which PROGENy pathways best explain the MEKi effect?
2. Which transcription factors support the same interpretation?
3. Which activity scores are visible both in the contrast-level results and at sample level?
4. Which findings agree with the Hallmark and GO results from Day 1?
5. Which findings are new compared with enrichment analysis?
6. Which conclusions should be treated cautiously because they depend on one model or one transcription factor?