MAX Network Deregulation Analysis in Sézary Syndrome

1 Setup and Libraries

# Data Processing & Single-Cell Analysis
library(dplyr)
library(tibble)
library(tidyr)
library(Seurat)
library(Matrix)

# Visualization
library(ggplot2)
library(patchwork)
library(ComplexHeatmap)
library(circlize)
library(grid)
library(EnhancedVolcano)

# Create output directories if they do not exist
dir.create("Output_Figures", showWarnings = FALSE)
dir.create("Output_Tables", showWarnings = FALSE)

2 Load Seurat Object and Define Conditions

# Load Seurat Object directly from cluster relative path
seurat_obj <- readRDS("../../0-Seurat_objects/Seurat_Object_With_TF_Activity.rds")
Idents(seurat_obj) <- "seurat_clusters"

# Define Malignant vs Non-Malignant (Clusters 3 & 10 = Normal CD4+ T reference)
non_malignant_clusters <- c(3, 10)
seurat_obj$Condition <- ifelse(seurat_obj$seurat_clusters %in% non_malignant_clusters, 
                               "Non-Malignant", 
                               "Malignant")
seurat_obj$Condition <- factor(seurat_obj$Condition, levels = c("Non-Malignant", "Malignant"))

print("Seurat Object Loaded successfully:")

[1] "Seurat Object Loaded successfully:"

table(seurat_obj$Condition)

Non-Malignant     Malignant 
         7873         41432 

3 Differential Statistics for MAX & Network Partners

# --- 1. Gene Expression (SCT Assay) ---
DefaultAssay(seurat_obj) <- "SCT"
Idents(seurat_obj) <- "Condition"

candidate_genes <- intersect(
  c("MAX", "MYC", "MLX", "MLXIP", "MLXIPL", "MNT", "MXD1", "MXI1", "E2F1", "CLIC1"),
  rownames(seurat_obj[["SCT"]])
)

deg_sct <- FindMarkers(
  seurat_obj,
  features = candidate_genes,
  ident.1 = "Malignant",
  ident.2 = "Non-Malignant",
  logfc.threshold = 0,
  min.pct = 0
)

  |                                                  | 0 % ~calculating  
  |+++++                                             | 10% ~00s          
  |++++++++++                                        | 20% ~00s          
  |+++++++++++++++                                   | 30% ~00s          
  |++++++++++++++++++++                              | 40% ~00s          
  |+++++++++++++++++++++++++                         | 50% ~00s          
  |++++++++++++++++++++++++++++++                    | 60% ~00s          
  |+++++++++++++++++++++++++++++++++++               | 70% ~00s          
  |++++++++++++++++++++++++++++++++++++++++          | 80% ~00s          
  |+++++++++++++++++++++++++++++++++++++++++++++     | 90% ~00s          
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=00s  

deg_sct$Gene <- rownames(deg_sct)
deg_sct$Assay <- "SCT (mRNA Expression)"

# --- 2. TF Activity (dorothea Assay) ---
DefaultAssay(seurat_obj) <- "dorothea"

candidate_tfs <- intersect(
  c("MAX", "MYC", "E2F1", "E2F4", "FOXM1", "FOXO1"),
  rownames(seurat_obj[["dorothea"]])
)

deg_tf <- FindMarkers(
  seurat_obj,
  features = candidate_tfs,
  ident.1 = "Malignant",
  ident.2 = "Non-Malignant",
  logfc.threshold = 0,
  min.pct = 0
)

  |                                                  | 0 % ~calculating  
  |+++++++++++++++++                                 | 33% ~01s          
  |++++++++++++++++++++++++++++++++++                | 67% ~00s          
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=00s  

deg_tf$Gene <- rownames(deg_tf)
deg_tf$Assay <- "dorothea (TF Activity)"

# Combine and save table for reviewer response
deg_summary <- bind_rows(deg_sct, deg_tf) %>%
  select(Assay, Gene, avg_log2FC, pct.1, pct.2, p_val, p_val_adj)

write.csv(deg_summary, "Output_Tables/MAX_Network_Differential_Summary.csv", row.names = FALSE)
print("=== Differential Statistics for MAX & Key Network Members ===")

[1] "=== Differential Statistics for MAX & Key Network Members ==="

print(deg_summary)

4 UMAP Feature Plots: MAX TF Activity vs. MAX Gene Expression

# --- Panel 1: MAX TF Activity (Red = dorothea Assay) ---
DefaultAssay(seurat_obj) <- "dorothea"
p_max_act <- FeaturePlot(seurat_obj, features = "MAX", order = TRUE, reduction = "umap") +
  scale_color_gradientn(colors = c("grey90", "firebrick")) +
  ggtitle("MAX Regulon Activity (dorothea)") +
  theme(plot.title = element_text(face = "bold", size = 12))

# --- Panel 2: MAX Gene Expression (Blue = SCT Assay) ---
DefaultAssay(seurat_obj) <- "SCT"
p_max_exp <- FeaturePlot(seurat_obj, features = "MAX", order = TRUE, reduction = "umap") +
  scale_color_gradientn(colors = c("grey90", "darkblue")) +
  ggtitle("MAX Gene Expression (SCT)") +
  theme(plot.title = element_text(face = "bold", size = 12))

# Combine side-by-side
fig_max_umap <- p_max_act | p_max_exp
fig_max_umap <- fig_max_umap + plot_annotation(
  title = "MAX Network Decoupling in Sézary Syndrome (Activity vs. mRNA Expression)",
  theme = theme(plot.title = element_text(face = "bold", size = 14))
)

# Display & Save
print(fig_max_umap)

ggsave("Output_Figures/Figure_MAX_UMAP_Activity_vs_Expression.pdf", plot = fig_max_umap, width = 11, height = 5)
ggsave("Output_Figures/Figure_MAX_UMAP_Activity_vs_Expression.png", plot = fig_max_umap, width = 11, height = 5, dpi = 300)

5 Violin Plots: Malignant vs. Non-Malignant Quantification

# Panel A: MAX Expression (SCT)
vln_exp <- VlnPlot(seurat_obj, features = "MAX", assay = "SCT", group.by = "Condition", pt.size = 0.05) +
  scale_fill_manual(values = c("Non-Malignant" = "#4DAF4A", "Malignant" = "#E41A1C")) +
  ylab("Normalized Expression") +
  ggtitle("MAX mRNA Expression (SCT)") +
  theme(legend.position = "none", plot.title = element_text(face = "bold", size = 11))

# Panel B: MAX Activity (dorothea)
vln_act <- VlnPlot(seurat_obj, features = "MAX", assay = "dorothea", group.by = "Condition", pt.size = 0.05) +
  scale_fill_manual(values = c("Non-Malignant" = "#4DAF4A", "Malignant" = "#E41A1C")) +
  ylab("TF Activity Score") +
  ggtitle("MAX TF Activity (dorothea)") +
  theme(legend.position = "none", plot.title = element_text(face = "bold", size = 11))

# Panel C: MYC Activity (dorothea) for Direct Contrast
vln_myc <- VlnPlot(seurat_obj, features = "MYC", assay = "dorothea", group.by = "Condition", pt.size = 0.05) +
  scale_fill_manual(values = c("Non-Malignant" = "#4DAF4A", "Malignant" = "#E41A1C")) +
  ylab("TF Activity Score") +
  ggtitle("MYC TF Activity (dorothea)") +
  theme(legend.position = "none", plot.title = element_text(face = "bold", size = 11))

fig_vln_max <- vln_exp | vln_act | vln_myc
print(fig_vln_max)

ggsave("Output_Figures/Figure_MAX_MYC_Violin_Comparison.pdf", plot = fig_vln_max, width = 12, height = 4)
ggsave("Output_Figures/Figure_MAX_MYC_Violin_Comparison.png", plot = fig_vln_max, width = 12, height = 4, dpi = 300)

6 Cluster-Level ComplexHeatmap (MAX vs. MYC and Core Regulators)

# Select targeted panel contrasting MYC/E2F oncogenic program with MAX and tumor suppressors
tf_panel <- c(
  # Oncogenic / Proliferative
  "MYC", "E2F1", "E2F4", "FOXM1", "SREBF1", "HIF1A",
  # MAX
  "MAX",
  # Normal Homeostasis & Tumor Suppressors
  "FOXO1", "FOXO4", "TCF7", "TBX21", "RUNX3"
)

# Extract TF Activity matrix (Seurat v4/v5 compatible)
mat_scaled <- tryCatch(
  SeuratObject::GetAssayData(seurat_obj, assay = "dorothea", layer = "scale.data"),
  error = function(e) NULL
)
mat_data <- SeuratObject::GetAssayData(seurat_obj, assay = "dorothea", layer = "data")
mat_use <- if (!is.null(mat_scaled) && nrow(mat_scaled) > 0) mat_scaled else mat_data

avail <- intersect(tf_panel, rownames(mat_use))
mat_use <- mat_use[avail, , drop = FALSE]

# Compute mean per cluster and calculate z-scores
clusters <- as.factor(seurat_obj$seurat_clusters)
avg_mat <- sapply(levels(clusters), function(cl) {
  Matrix::rowMeans(mat_use[, clusters == cl, drop = FALSE])
})
colnames(avg_mat) <- levels(clusters)

avg_mat_z <- t(scale(t(avg_mat)))
avg_mat_z[is.na(avg_mat_z)] <- 0

# Set column annotations
cluster_status <- ifelse(colnames(avg_mat_z) %in% c("3", "10"), 
                         "Normal CD4 T", 
                         "Malignant CD4 T cells")

ha_col <- HeatmapAnnotation(
  Cell_State = cluster_status,
  col = list(Cell_State = c("Normal CD4 T" = "#4DAF4A", "Malignant CD4 T cells" = "#E41A1C")),
  annotation_name_side = "left"
)

# Color scale: -3 to +3
col_fun <- circlize::colorRamp2(c(-3, 0, 3), c("#313695", "white", "#A50026"))

ht_max <- Heatmap(
  avg_mat_z,
  name = "TF activity (z)",
  col = col_fun,
  top_annotation = ha_col,
  column_split = cluster_status,
  cluster_rows = TRUE,
  cluster_columns = TRUE,
  show_row_dend = TRUE,
  show_column_dend = TRUE,
  row_names_gp = gpar(fontsize = 11, fontface = "bold"),
  column_names_gp = gpar(fontsize = 10),
  column_title = "Cluster-Level TF Activity: MYC/E2F Activation vs. MAX Repression",
  heatmap_legend_param = list(direction = "vertical")
)

# Save and Draw
pdf("Output_Figures/Figure_MAX_MYC_Cluster_Heatmap.pdf", width = 10, height = 7)
draw(ht_max, merge_legend = TRUE)
dev.off()

null device 
          1 

png("Output_Figures/Figure_MAX_MYC_Cluster_Heatmap.png", width = 10 * 300, height = 7 * 300, res = 300)

draw(ht_max, merge_legend = TRUE)
dev.off()

null device 
          1 

draw(ht_max, merge_legend = TRUE)

7 Volcano Plot Highlighting MAX & MYC Inferred Activity

# Differential TF Activity across all TFs
DefaultAssay(seurat_obj) <- "dorothea"
Idents(seurat_obj) <- "Condition"

all_tfs_diff <- FindMarkers(
  seurat_obj, 
  ident.1 = "Malignant", 
  ident.2 = "Non-Malignant", 
  logfc.threshold = 0, 
  min.pct = 0
)

  |                                                  | 0 % ~calculating  
  |+                                                 | 2 % ~10s          
  |++                                                | 3 % ~09s          
  |+++                                               | 5 % ~09s          
  |++++                                              | 7 % ~09s          
  |+++++                                             | 8 % ~09s          
  |++++++                                            | 10% ~09s          
  |++++++                                            | 12% ~09s          
  |+++++++                                           | 14% ~08s          
  |++++++++                                          | 15% ~08s          
  |+++++++++                                         | 17% ~08s          
  |++++++++++                                        | 19% ~08s          
  |+++++++++++                                       | 20% ~08s          
  |++++++++++++                                      | 22% ~08s          
  |++++++++++++                                      | 24% ~08s          
  |+++++++++++++                                     | 25% ~07s          
  |++++++++++++++                                    | 27% ~07s          
  |+++++++++++++++                                   | 29% ~07s          
  |++++++++++++++++                                  | 31% ~07s          
  |+++++++++++++++++                                 | 32% ~07s          
  |+++++++++++++++++                                 | 34% ~07s          
  |++++++++++++++++++                                | 36% ~07s          
  |+++++++++++++++++++                               | 37% ~07s          
  |++++++++++++++++++++                              | 39% ~07s          
  |+++++++++++++++++++++                             | 41% ~06s          
  |++++++++++++++++++++++                            | 42% ~06s          
  |+++++++++++++++++++++++                           | 44% ~06s          
  |+++++++++++++++++++++++                           | 46% ~06s          
  |++++++++++++++++++++++++                          | 47% ~06s          
  |+++++++++++++++++++++++++                         | 49% ~06s          
  |++++++++++++++++++++++++++                        | 51% ~06s          
  |+++++++++++++++++++++++++++                       | 53% ~06s          
  |++++++++++++++++++++++++++++                      | 54% ~05s          
  |++++++++++++++++++++++++++++                      | 56% ~05s          
  |+++++++++++++++++++++++++++++                     | 58% ~05s          
  |++++++++++++++++++++++++++++++                    | 59% ~05s          
  |+++++++++++++++++++++++++++++++                   | 61% ~05s          
  |++++++++++++++++++++++++++++++++                  | 63% ~04s          
  |+++++++++++++++++++++++++++++++++                 | 64% ~04s          
  |++++++++++++++++++++++++++++++++++                | 66% ~04s          
  |++++++++++++++++++++++++++++++++++                | 68% ~04s          
  |+++++++++++++++++++++++++++++++++++               | 69% ~04s          
  |++++++++++++++++++++++++++++++++++++              | 71% ~03s          
  |+++++++++++++++++++++++++++++++++++++             | 73% ~03s          
  |++++++++++++++++++++++++++++++++++++++            | 75% ~03s          
  |+++++++++++++++++++++++++++++++++++++++           | 76% ~03s          
  |+++++++++++++++++++++++++++++++++++++++           | 78% ~03s          
  |++++++++++++++++++++++++++++++++++++++++          | 80% ~02s          
  |+++++++++++++++++++++++++++++++++++++++++         | 81% ~02s          
  |++++++++++++++++++++++++++++++++++++++++++        | 83% ~02s          
  |+++++++++++++++++++++++++++++++++++++++++++       | 85% ~02s          
  |++++++++++++++++++++++++++++++++++++++++++++      | 86% ~02s          
  |+++++++++++++++++++++++++++++++++++++++++++++     | 88% ~01s          
  |+++++++++++++++++++++++++++++++++++++++++++++     | 90% ~01s          
  |++++++++++++++++++++++++++++++++++++++++++++++    | 92% ~01s          
  |+++++++++++++++++++++++++++++++++++++++++++++++   | 93% ~01s          
  |++++++++++++++++++++++++++++++++++++++++++++++++  | 95% ~01s          
  |+++++++++++++++++++++++++++++++++++++++++++++++++ | 97% ~00s          
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 98% ~00s          
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=12s  

# Key factors to highlight
tfs_to_label <- c("MAX", "MYC", "E2F1", "E2F4", "FOXM1", "FOXO1", "RELA", "IRF1", "STAT3", "RUNX3")

p_volcano_max <- EnhancedVolcano(
  all_tfs_diff,
  lab = rownames(all_tfs_diff),
  selectLab = tfs_to_label,
  x = 'avg_log2FC',
  y = 'p_val_adj',
  title = 'Differential TF Activity: Malignant vs. Non-Malignant CD4+ T Cells',
  subtitle = 'DecoupleR Inferred Activity Highlighting the MAX/MYC Axis',
  pCutoff = 1e-5,
  FCcutoff = 0.5,
  pointSize = 3.0,
  labSize = 4.5,
  drawConnectors = TRUE,
  col = c("grey60", "forestgreen", "royalblue", "firebrick2")
)

print(p_volcano_max)

ggsave("Output_Figures/Figure_MAX_Volcano_TF_Activity.pdf", plot = p_volcano_max, width = 10, height = 8)
ggsave("Output_Figures/Figure_MAX_Volcano_TF_Activity.png", plot = p_volcano_max, width = 10, height = 8, dpi = 300)

8 Completion Summary

print("MAX network analysis complete. All tables exported to Output_Tables/ and figures to Output_Figures/.")

[1] "MAX network analysis complete. All tables exported to Output_Tables/ and figures to Output_Figures/."