TF Activity Inference Analysis Heatmaps

1 load libraries

# Data Processing
library(dplyr)
library(Seurat)
library(tibble)
library(tidyr)
library(stringr)

# Visualization
library(ggplot2)
library(ComplexHeatmap)
library(patchwork)
library(SCpubr)

# Regulatory Network Inference
library(decoupleR)
library(dorothea)
data(dorothea_hs, package = "dorothea")
library(tictoc)

2 Load Seurat Object

# Load your Seurat Object
seurat_obj <- readRDS("../Output_Objects/Seurat_Object_With_TF_Activity.rds")

Idents(seurat_obj) <- "seurat_clusters"
print("Object Loaded.")

[1] "Object Loaded."

2.1 Run this code block to restore activities instantly:

# If 'activities' is missing but 'dorothea' assay exists, reconstruct it:
if (!exists("activities") && "dorothea" %in% names(seurat_obj@assays)) {
  
  print("Reconstructing 'activities' dataframe from Seurat object...")
  
  # Extract the matrix (Seurat v5 uses 'layer' instead of 'slot')
  # Since you ran ScaleData, we use 'scale.data'
  mat <- GetAssayData(seurat_obj, assay = "dorothea", layer = "scale.data")
  
  # Convert to long format (what SCpubr needs)
  activities <- as.data.frame(mat) %>%
    rownames_to_column("source") %>%
    pivot_longer(cols = -source, names_to = "condition", values_to = "score") %>%
    mutate(statistic = "norm_wmean") # SCpubr requires this column
    
  print("Activities dataframe restored!")
}

[1] "Reconstructing 'activities' dataframe from Seurat object..."
[1] "Activities dataframe restored!"

2.2 SCpubr Heatmap Visualization-Heatmap of averaged scores

2.3 Set the scale limits

out <- SCpubr::do_TFActivityHeatmap(sample = seurat_obj,
                                 activities = activities,
                                 min.cutoff = -1.5,
                                 max.cutoff = 1.5)

print(out)

# Save ComplexHeatmap properly
pdf("Output_Figures/SCpubr_Heatmap_Scaled.pdf", width = 10, height = 8)

print(out)
dev.off()

png 
  2 

png("Output_Figures/SCpubr_Heatmap_Scaled.png", width = 10 * 300, height = 8 * 300, res = 300)
print(out)

dev.off()

png 
  2 

2.4 Enforce Symmetry (Best for Manuscript)

out <- SCpubr::do_TFActivityHeatmap(sample = seurat_obj,
                                 activities = activities,
                                 min.cutoff = -1.5,
                                 max.cutoff = 1.5,
                                 enforce_symmetry = TRUE)

print(out)

pdf("Output_Figures/SCpubr_Heatmap_Symmetric.pdf", width = 10, height = 8)

print(out)
dev.off()

png 
  2 

png("Output_Figures/SCpubr_Heatmap_Symmetric.png", width = 10 * 300, height = 8 * 300, res = 300)
print(out)

dev.off()

png 
  2 

print(out)

2.5 Top 40 TFs

out <- SCpubr::do_TFActivityHeatmap(sample = seurat_obj,
                                 activities = activities,
                                 n_tfs = 40)

print(out)

pdf("Output_Figures/SCpubr_Heatmap_Top40.pdf", width = 14, height = 6)

print(out)
dev.off()

png 
  2 

png("Output_Figures/SCpubr_Heatmap_Top40.png", width = 14 * 300, height = 6 * 300, res = 300)
print(out)

dev.off()

png 
  2 

2.6 Top 100 TFs (Figure A for Manuscript)

out <- SCpubr::do_TFActivityHeatmap(sample = seurat_obj,
                                 activities = activities,
                                 n_tfs = 100)

print(out)

pdf("Output_Figures/Figure_3.16A_Global_TF_Heatmap_Top100.pdf", width = 32, height = 12)

print(out)
dev.off()

png 
  2 

png("Output_Figures/Figure_3.16A_Global_TF_Heatmap_Top100.png", width = 32 * 300, height = 12 * 300, res = 300)
print(out)

dev.off()

png 
  2 

2.7 Top 100 TFs (Figure A for Manuscript)

out <- SCpubr::do_TFActivityHeatmap(sample = seurat_obj,
                                 activities = activities,
                                 min.cutoff = -1.7,
                                 max.cutoff = 1.7, group.by = "seurat_clusters",
                                 n_tfs = 100)

print(out)

pdf("Output_Figures/Figure_Top100.pdf", width = 32, height = 12)

print(out)
dev.off()

png 
  2 

png("Output_Figures/Figure_Top100.png", width = 32 * 300, height = 12 * 300, res = 300)
print(out)

dev.off()

png 
  2 

3 Differential TF Activity (Malignant vs. Normal)

# Define Comparison: Clusters 3 & 10 (Normal) vs Rest (Malignant)
non_malignant_clusters <- c(3, 10)
seurat_obj$Condition <- ifelse(seurat_obj$seurat_clusters %in% non_malignant_clusters, "Non-Malignant", "Malignant")

# Perform Differential Analysis on TF Activity
DefaultAssay(seurat_obj) <- "dorothea"
Idents(seurat_obj) <- "Condition"

print("Running FindMarkers on TF Activity...")

[1] "Running FindMarkers on TF Activity..."

diff_tfs <- FindMarkers(seurat_obj, 
                        ident.1 = "Malignant", 
                        ident.2 = "Non-Malignant", 
                        logfc.threshold = 0, # Get all for volcano
                        min.pct = 0)

# Add gene column for labeling
diff_tfs$gene <- rownames(diff_tfs)

# Save Results
write.csv(diff_tfs, "Output_Tables/Differential_TF_Activity_Malignant_vs_Normal.csv")
print("Differential analysis complete.")

[1] "Differential analysis complete."

4 Figure C: Volcano Plot (Loss of Homeostasis)

# Highlight key drivers mentioned in text
highlight_tfs <- c("FOXO1", "MYC", "E2F1", "E2F4", "FOXM1", "RELA", "IRF1", "STAT1")

p_volcano <- SCpubr::do_VolcanoPlot(sample = seurat_obj,
                                    de_genes = diff_tfs
                                   )

ggsave("Output_Figures/Figure_3.16C_Volcano_TF_Activity.pdf", plot = p_volcano, width = 8, height = 6)
ggsave("Output_Figures/Figure_3.16C_Volcano_TF_Activity.png", plot = p_volcano, width = 8, height = 6, dpi = 300)
print(p_volcano)

5 Updated Figure C: EnhancedVolcano

library(EnhancedVolcano)

# Highlight key drivers mentioned in text
highlight_tfs <- c("FOXO1", "MYC", "E2F1", "E2F4", "FOXM1", "RELA", "IRF1", "STAT1", "TOX", "GATA3")

# Create the EnhancedVolcano Plot
p_volcano <- EnhancedVolcano(diff_tfs,
    lab = rownames(diff_tfs),
    x = 'avg_log2FC',
    y = 'p_val_adj',
    
    title = 'Differential TF Activity: Malignant vs. Non-Malignant',
    subtitle = 'DecoupleR Inferred Activity',
    pCutoff = 1e-5,
    FCcutoff = 0.5,
    pointSize = 3.0,
    labSize = 5.0,
    colAlpha = 0.8,
    legendPosition = 'right',
    legendLabSize = 12,
    legendIconSize = 4.0,
    drawConnectors = TRUE, # Draw lines to labels to avoid overlap
    widthConnectors = 0.5,
    colConnectors = 'grey30',
    # Custom Colors: Down (Blue), Up (Red), NS (Grey)
    col = c("grey30", "forestgreen", "royalblue", "firebrick2")
)

# Print
print(p_volcano)

# Save
ggsave("Output_Figures/Figure_3.16C_EnhancedVolcano_TF_Activity.pdf", plot = p_volcano, width = 10, height = 8)
ggsave("Output_Figures/Figure_3.16C_EnhancedVolcano_TF_Activity.png", plot = p_volcano, width = 10, height = 8, dpi = 300)

6 Figure D: Mixed Feature Plots (Activity vs Expression)

# We manually construct this to mix Assays

# Part 1: TF Activity Plots (Assay: dorothea)
DefaultAssay(seurat_obj) <- "dorothea"

p1 <- FeaturePlot(seurat_obj, features = "FOXO1", order = T, reduction = "umap") + 
      scale_color_gradientn(colors = c("grey90", "firebrick")) + ggtitle("FOXO1 Activity (Homeostasis)")
p2 <- FeaturePlot(seurat_obj, features = "RELA", order = T, reduction = "umap") + 
      scale_color_gradientn(colors = c("grey90", "firebrick")) + ggtitle("RELA Activity (Inflammatory)")
p3 <- FeaturePlot(seurat_obj, features = "IRF1", order = T, reduction = "umap") + 
      scale_color_gradientn(colors = c("grey90", "firebrick")) + ggtitle("IRF1 Activity (IFN-Response)")
p4 <- FeaturePlot(seurat_obj, features = "FOXM1", order = T, reduction = "umap") + 
      scale_color_gradientn(colors = c("grey90", "firebrick")) + ggtitle("FOXM1 Activity (Proliferation)")

# Part 2: Gene Expression Plots (Assay: SCT/RNA)
DefaultAssay(seurat_obj) <- "SCT"

p5 <- FeaturePlot(seurat_obj, features = "HMGA2", order = T, reduction = "umap") + 
      scale_color_gradientn(colors = c("grey90", "darkblue")) + ggtitle("HMGA2 Expression (Stem-like)")
p6 <- FeaturePlot(seurat_obj, features = "SOX4", order = T, reduction = "umap") + 
      scale_color_gradientn(colors = c("grey90", "darkblue")) + ggtitle("SOX4 Expression (Stem-like)")

# Combine
final_figure_D <- (p1 | p2 | p3) / (p4 | p5 | p6) + 
                  plot_annotation(title = "Figure 3.16D: Key Drivers (Red=Activity, Blue=Expression)")

ggsave("Output_Figures/Figure_3.16D_Mixed_Features.pdf", plot = final_figure_D, width = 14, height = 10)
ggsave("Output_Figures/Figure_3.16D_Mixed_Features.png", plot = final_figure_D, width = 14, height = 10, dpi = 300)
print(final_figure_D)

7 Figure E (ComplexHeatmap) chunk

library(ComplexHeatmap)
library(circlize)
library(Matrix)

# Expanded list of state-specific drivers based on your regulon analysis
literature_tfs <- c(
  "GATA3", "STAT6", "BATF", "FOXP3", "STAT3", "STAT5B", "TCF7", # Core/Memory
  "E2F1", "MYC", "FOXM1",                                       # Proliferation (Cl 7)
  "STAT1", "STAT2", "IRF1", "IRF9",                             # IFN-stimulated (Cl 13)
  "RELA", "NFKB1", "REL", "FOS",                                # Pro-inflammatory (Cl 11, 12)
  "TBX21", "RUNX3",                                             # Cytotoxic (Cl 1, 9)
  "HIF1A", "SREBF1",                                            # Metabolic shift (Cl 8)
  "RFX5", "SPI1"                                                # MHC-II High (Cl 0)
)


# Keep only TFs present in the dorothea assay
available_tfs <- intersect(literature_tfs, rownames(seurat_obj[["dorothea"]]))
if (length(available_tfs) < 5) stop("Too few TFs found in dorothea assay. Check TF naming / assay content.")

# Extract TF activity matrix (TFs x cells)
# Use scale.data if available; otherwise fall back to data layer.
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
mat_use <- mat_use[available_tfs, , drop = FALSE]

# Average per cluster (TF x cluster)
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)

# Optional: z-score across clusters (helps readability if you used raw 'data' instead of 'scale.data')
avg_mat_z <- t(scale(t(avg_mat)))
avg_mat_z[is.na(avg_mat_z)] <- 0

# Colors
col_fun <- circlize::colorRamp2(c(-2, 0, 2), c("#313695", "white", "#A50026"))

ht <- Heatmap(
  avg_mat_z,
  name = "TF activity (z)",
  col = col_fun,
  cluster_rows = TRUE,
  cluster_columns = TRUE,
  show_row_dend = TRUE,
  show_column_dend = TRUE,
  row_names_gp = grid::gpar(fontsize = 10),
  column_names_gp = grid::gpar(fontsize = 10),
  column_title = "Literature-validated Sézary TF modules (DoRothEA/decoupleR)",
  heatmap_legend_param = list(direction = "vertical")
)

# Draw to notebook
draw(ht)

# Save PDF (vector)
pdf("Output_Figures/Figure_3.16E_Literature_TF_Heatmap_ComplexHeatmap.pdf", width = 10, height = 8)

draw(ht)
dev.off()

png 
  2 

# Save PNG (raster, publication-ready)
png("Output_Figures/Figure_3.16E_Literature_TF_Heatmap_ComplexHeatmap.png",
    width = 10 * 300, height = 8 * 300, res = 300)
draw(ht)

dev.off()

png 
  2 

8 Figure F (ComplexHeatmap) chunk

library(ComplexHeatmap)
library(circlize)
library(Matrix)

# Expanded list including FOXO1 and tumor suppressors
literature_tfs <- c(
  # Top Malignant Upregulated (Oncogenic, Stress, Proliferation)
  "RFX5", "MYC", "E2F4", "HSF1", "SREBF2", "NFE2L2", 
  "RELA", "REL", "NFKB1", "IRF1", "NCOA2",
  
  # Malignant Downregulated / Normal Enriched (Tumor Suppressors & Homeostasis)
  "FOXO1", "FOXO4", "RUNX3", "TCF3", "BCL11A", "NEUROD1", "MEF2B", "PBX2",
  
  # UMAP State Drivers (Intra-tumoral heterogeneity)
  "GATA3", "STAT6", "BATF", "FOXP3", "STAT3", "STAT5B", "TCF7", # Core/Memory
  "E2F1", "FOXM1",                                              # Proliferation
  "STAT1", "STAT2", "IRF9",                                     # IFN response
  "HIF1A", "SREBF1",                                            # Metabolic
  "TBX21"                                                       # Cytotoxic
)

# Extract TF activity matrix
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
available_tfs <- intersect(literature_tfs, rownames(mat_use))
if (length(available_tfs) < 5) stop("Too few TFs found. Check assay data.")
mat_use <- mat_use[available_tfs, , drop = FALSE]

# Average per cluster and z-score
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

# Annotate Malignant vs Normal (Clusters 3, 10 = Normal)
cluster_status <- ifelse(colnames(avg_mat_z) %in% c("3", "10"), 
                         "Normal CD4 T", 
                         "Malignant CD4 T cells")

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

# Colors
col_fun <- circlize::colorRamp2(c(-3, 0, 3), c("#313695", "white", "#A50026"))

# Create heatmap with column split
ht <- Heatmap(
  avg_mat_z,
  name = "TF activity (z)",
  col = col_fun,
  top_annotation = ha,
  column_split = cluster_status, # Physically splits normal and malignant columns
  cluster_rows = TRUE,
  cluster_columns = TRUE,
  show_row_dend = TRUE,
  show_column_dend = TRUE,
  row_names_gp = grid::gpar(fontsize = 10),
  column_names_gp = grid::gpar(fontsize = 10),
  column_title = "Differential TF Modules in Sézary Heterogeneity",
  heatmap_legend_param = list(direction = "vertical")
)

# Output
pdf("Output_Figures/Figure_3.16E_Differential_TF_Heatmap.pdf", width = 11, height = 9)
draw(ht)
dev.off()

null device 
          1 

png("Output_Figures/Figure_3.16E_Differential_TF_Heatmap.png",
    width = 11 * 300, height = 9 * 300, res = 300)

draw(ht)
dev.off()

null device 
          1 

draw(ht)

9 Figure G (ComplexHeatmap) chunk

# ============================================
# LIBRARIES
# ============================================
library(ComplexHeatmap)
library(circlize)
library(Matrix)
library(grid)
library(SeuratObject)  # for GetAssayData

# ============================================
# 1. Define TF panel metadata (44-47 TFs)
# ============================================
tf_meta <- data.frame(
  TF = c(
    "MYC","E2F4","RFX5","TWIST1","JUNB","IRF4","CREB1",
    "FOS","FOSL1",
    "HSF1","NFE2L2","SREBF2",
    "RELA","REL","NFKB1","IRF1","NCOA2",
    "NFATC1","NFATC2",
    "FOXO1","FOXO4","RUNX3","ZEB1","BACH2",
    "TCF3","BCL11A","NEUROD1","MEF2B","PBX2","IRF3","BCL6",
    "GATA3","STAT6","BATF","FOXP3","STAT3","STAT5B","TCF7",
    "E2F1","FOXM1",
    "STAT1","STAT2","IRF9",
    "HIF1A","SREBF1",
    "EOMES",
    "PRDM1"
  ),
  Condition = c(
    rep("Malignant", 7),  # Oncogenic
    rep("Malignant", 2),  # AP-1
    rep("Malignant", 3),  # Stress
    rep("Malignant", 5),  # NF-kB
    rep("Malignant", 2),  # NFAT
    rep("Normal", 5),     # Tumor Suppressor
    rep("Normal", 7),     # Homeostasis
    rep("Malignant", 7),  # Th2/Memory Core
    rep("Malignant", 2),  # Proliferation
    rep("Malignant", 3),  # IFN
    rep("Malignant", 2),  # Metabolism
    rep("Malignant", 1),  # Cytotoxic
    rep("Malignant", 1)   # Terminal Effector
  ),
  Function = c(
    rep("Oncogenic", 7),
    rep("AP-1 signaling", 2),
    rep("Stress Response", 3),
    rep("Inflammatory/NF-kB", 5),
    rep("NFAT signaling", 2),
    rep("Tumor Suppressor", 5),
    rep("Normal Homeostasis", 7),
    rep("Th2/Memory Core", 7),
    rep("Proliferation", 2),
    rep("IFN Response", 3),
    rep("Metabolism", 2),
    rep("Cytotoxic", 1),
    rep("Terminal Effector", 1)
  ),
  stringsAsFactors = FALSE
)

# ============================================
# 2. Extract TF activity matrix from Seurat
# ============================================
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

# Filter for TFs present in Seurat
available_tfs <- intersect(tf_meta$TF, rownames(mat_use))
mat_use <- mat_use[available_tfs, , drop = FALSE]

# ============================================
# 3. Average per cluster & z-score
# ============================================
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

# Align tf_meta order
tf_meta_filtered <- tf_meta[match(rownames(avg_mat_z), tf_meta$TF), ]

# ============================================
# 4. Column & row annotations
# ============================================
cluster_status <- ifelse(colnames(avg_mat_z) %in% c("3","10"), "Normal CD4 T cells", "Malignant CD4 T cells")
ha_col <- HeatmapAnnotation(
  Cell_State = cluster_status,
  col = list(Cell_State = c("Normal CD4 T cells" = "#4DAF4A", "Malignant CD4 T cells" = "#E41A1C")),
  annotation_name_side = "left"
)

# Row annotation
function_colors <- c(
  "Oncogenic" = "#FF7F00",
  "AP-1 signaling" = "#FFA500",
  "Stress Response" = "#FFD700",
  "Inflammatory/NF-kB" = "#1E90FF",
  "NFAT signaling" = "#4169E1",
  "Tumor Suppressor" = "#377EB8",
  "Normal Homeostasis" = "#4DAF4A",
  "Th2/Memory Core" = "#984EA3",
  "Proliferation" = "#E41A1C",
  "IFN Response" = "#00CED1",
  "Metabolism" = "#A65628",
  "Cytotoxic" = "#F781BF",
  "Terminal Effector" = "#800080"
)

ha_row <- rowAnnotation(
  Condition = tf_meta_filtered$Condition,
  Function = tf_meta_filtered$Function,
  col = list(
    Condition = c("Normal"="#4DAF4A","Malignant"="#E41A1C"),
    Function = function_colors
  ),
  annotation_name_side = "bottom"
)

# ============================================
# 5. Heatmap colors & plotting
# ============================================
col_fun <- circlize::colorRamp2(c(-3,0,3), c("#313695","white","#A50026"))

ht <- Heatmap(
  avg_mat_z,
  name = "TF activity (z)",
  col = col_fun,
  top_annotation = ha_col,
  left_annotation = ha_row,
  column_split = cluster_status,
  row_split = tf_meta_filtered$Function,
  cluster_rows = FALSE,     # show biological split
  cluster_columns = TRUE,
  show_row_dend = FALSE,
  show_column_dend = TRUE,
  row_names_gp = gpar(fontsize = 10),
  column_names_gp = gpar(fontsize = 10),
  column_title = "Functional TF Modules in Sézary Syndrome",
  row_title_rot = 0,
  row_title_gp = gpar(fontsize = 9, fontface = "bold"),
  heatmap_legend_param = list(direction = "vertical")
)

# ============================================
# 6. Save plots
# ============================================
pdf("Output_Figures/Figure_TF_Heatmap.pdf", width=12, height=10)
draw(ht, merge_legend=TRUE)
dev.off()

null device 
          1 

png("Output_Figures/Figure_TF_Heatmap.png", width=12*300, height=10*300, res=300)

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

null device 
          1 

draw(ht, merge_legend=TRUE)

10 Figure Malignant complex heatmap chunk (with 44-47 TFs, literature-based panel)

# ============================================
# LIBRARIES
# ============================================
library(ComplexHeatmap)
library(circlize)
library(Matrix)
library(grid)
library(SeuratObject)

# ============================================
# 1. Define TF panel metadata (KEGG Aligned)
# ============================================
tf_meta <- data.frame(
  TF = c(
    # --- 1. General Malignancy & Stress ---
    "MYC","E2F4","TWIST1","IRF4",
    "HSF1","NFE2L2","SREBF2",
    
    # --- 2. TCR Signaling Triad ---
    "JUNB","FOS","FOSL1",           # AP-1
    "NFATC1","NFATC2",              # NFAT
    "RELA","REL","NFKB1","IRF1","NCOA2", # NF-kB
    
    # --- 3. Th2 / JAK-STAT Core ---
    "GATA3","STAT6","BATF","FOXP3","STAT3","STAT5B",
    
    # --- 4. Differentiation Hierarchy ---
    "TCF7","LEF1","MYB",            # Stem-like Progenitor
    "E2F1","FOXM1",                 # Proliferation (Cycling)
    "PRDM1",                        # Terminal Effector
    
    # --- 5. KEGG ALIGNED CATEGORIES (NEW) ---
    "EOMES","TBX21","RUNX3",        # NK-like Cytotoxicity (Cluster 1,9)
    "RFX5","CREB1",                 # Antigen Presentation / MHC-II (Cluster 0)
    "KLF4","ETS1","SMAD3",          # Migration / Cell Adhesion (CAMs)
    
    # --- 6. Microenvironment ---
    "STAT1","STAT2","IRF9",         # IFN Response
    "HIF1A","SREBF1",               # Metabolism
    
    # --- 7. Normal Baseline / Tumor Suppressors ---
    "FOXO1","FOXO4","ZEB1","BACH2",
    "TCF3","BCL11A","NEUROD1","MEF2B","PBX2","IRF3","BCL6"
  ),
  Condition = c(
    rep("Malignant", 4),   # Oncogenic
    rep("Malignant", 3),   # Stress
    
    rep("Malignant", 3),   # AP-1
    rep("Malignant", 2),   # NFAT
    rep("Malignant", 5),   # NF-kB
    
    rep("Malignant", 6),   # Th2 / JAK-STAT Core
    
    rep("Malignant", 3),   # Stem-like Progenitor
    rep("Malignant", 2),   # Proliferation
    rep("Malignant", 1),   # Terminal Effector
    
    rep("Malignant", 3),   # NK-like Cytotoxicity
    rep("Malignant", 2),   # Antigen Presentation
    rep("Malignant", 3),   # Migration / Adhesion
    
    rep("Malignant", 3),   # IFN
    rep("Malignant", 2),   # Metabolism
    
    rep("Normal", 4),      # Tumor Suppressor
    rep("Normal", 7)       # Homeostasis
  ),
  Function = c(
    rep("Oncogenic", 4),
    rep("Stress Response", 3),
    
    rep("AP-1 Signaling", 3),
    rep("NFAT Signaling", 2),
    rep("Inflammatory/NF-kB", 5),
    
    rep("Th2 / JAK-STAT Core", 6),
    
    rep("Stem-like Progenitor", 3),
    rep("Proliferation", 2),
    rep("Terminal Effector", 1),
    
    rep("NK-like Cytotoxicity", 3),   # KEGG aligned
    rep("Antigen Presentation", 2),   # KEGG aligned
    rep("Migration / Adhesion", 3),   # KEGG aligned
    
    rep("IFN Response", 3),
    rep("Metabolism", 2),
    
    rep("Tumor Suppressor", 4),
    rep("Normal Homeostasis", 7)
  ),
  stringsAsFactors = FALSE
)

# Lock in the precise order of the blocks from top to bottom
desired_order <- c(
  "Oncogenic", 
  "Stress Response",
  "AP-1 Signaling", 
  "NFAT Signaling", 
  "Inflammatory/NF-kB",
  "Th2 / JAK-STAT Core",
  "Stem-like Progenitor", 
  "Proliferation",        
  "Terminal Effector",
  "NK-like Cytotoxicity",   # Placed here to show effector state
  "Antigen Presentation",   # Directly links to MHC-II high cluster
  "Migration / Adhesion",   # Links to CAMs KEGG pathway
  "IFN Response",
  "Metabolism",
  "Tumor Suppressor",
  "Normal Homeostasis"
)

tf_meta$Function <- factor(tf_meta$Function, levels = desired_order)

# ============================================
# 2. Extract TF activity matrix from Seurat
# ============================================
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

# Filter for TFs present in Seurat
available_tfs <- intersect(tf_meta$TF, rownames(mat_use))
mat_use <- mat_use[available_tfs, , drop = FALSE]

# ============================================
# 3. Average per cluster & z-score
# ============================================
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

# Align tf_meta order to match mat_use precisely
row_order_idx <- match(rownames(avg_mat_z), tf_meta$TF)
tf_meta_filtered <- tf_meta[row_order_idx, ]
avg_mat_z <- avg_mat_z[order(tf_meta_filtered$Function), ]
tf_meta_filtered <- tf_meta_filtered[order(tf_meta_filtered$Function), ]

# ============================================
# 4. Column & row annotations
# ============================================
cluster_status <- ifelse(colnames(avg_mat_z) %in% c("3","10"), "Normal CD4 T cells", "Malignant CD4 T cells")
ha_col <- HeatmapAnnotation(
  Cell_State = cluster_status,
  col = list(Cell_State = c("Normal CD4 T cells" = "#4DAF4A", "Malignant CD4 T cells" = "#E41A1C")),
  annotation_name_side = "left"
)

# Row annotation colors
function_colors <- c(
  "Oncogenic" = "#808080",
  "Stress Response" = "#FFD700",
  
  "AP-1 Signaling" = "#FF8C00",
  "NFAT Signaling" = "#FF4500",
  "Inflammatory/NF-kB" = "#E31A1C",
  
  "Th2 / JAK-STAT Core" = "#984EA3",
  
  "Stem-like Progenitor" = "#FF1493",
  "Proliferation" = "#1E90FF",
  "Terminal Effector" = "#800080",
  
  # New KEGG categories
  "NK-like Cytotoxicity" = "#F781BF",   # Pink
  "Antigen Presentation" = "#66CDAA",   # Medium Aquamarine
  "Migration / Adhesion" = "#8A2BE2",   # Blue Violet
  
  "IFN Response" = "#00CED1",
  "Metabolism" = "#A65628",
  
  "Tumor Suppressor" = "#377EB8",
  "Normal Homeostasis" = "#4DAF4A"
)

ha_row <- rowAnnotation(
  Condition = tf_meta_filtered$Condition,
  Function = tf_meta_filtered$Function,
  col = list(
    Condition = c("Normal"="#4DAF4A","Malignant"="#E41A1C"),
    Function = function_colors
  ),
  annotation_name_side = "bottom"
)

# ============================================
# 5. Heatmap colors & plotting
# ============================================
col_fun <- circlize::colorRamp2(c(-3,0,3), c("#313695","white","#A50026"))

ht <- Heatmap(
  avg_mat_z,
  name = "TF activity (z)",
  col = col_fun,
  top_annotation = ha_col,
  left_annotation = ha_row,
  column_split = cluster_status,
  row_split = tf_meta_filtered$Function,
  cluster_row_slices = FALSE,    
  cluster_rows = FALSE,          
  cluster_columns = TRUE,
  show_row_dend = FALSE,
  show_column_dend = TRUE,
  row_names_gp = gpar(fontsize = 10),
  column_names_gp = gpar(fontsize = 10),
  column_title = "Functional TF Modules in Sézary Syndrome",
  row_title_rot = 0,
  row_title_gp = gpar(fontsize = 8, fontface = "bold"),
  heatmap_legend_param = list(direction = "vertical")
)

# ============================================
# 6. Save plots
# ============================================
pdf("Output_Figures/Figure_TF_Heatmap_KEGG.pdf", width=14, height=14)
draw(ht, merge_legend=TRUE)
dev.off()

null device 
          1 

draw(ht, merge_legend=TRUE)

11 Figure Malignant complex heatmap chunk (with 44-47 TFs, literature-based panel)

# ============================================
# LIBRARIES
# ============================================
library(ComplexHeatmap)
library(circlize)
library(Matrix)
library(grid)
library(SeuratObject)

# ============================================
# 1. Define TF panel metadata (KEGG Aligned)
# ============================================
tf_meta <- data.frame(
  TF = c(
    # --- 1. General Malignancy & Stress ---
    "MYC","E2F4","TWIST1","IRF4",
    "HSF1","NFE2L2","SREBF2",
    
    # --- 2. TCR Signaling Triad ---
    "JUNB","FOS","FOSL1",           # AP-1
    "NFATC1","NFATC2",              # NFAT
    "RELA","REL","NFKB1","IRF1","NCOA2", # NF-kB
    
    # --- 3. Th2 / JAK-STAT Core ---
    "GATA3","STAT6","BATF","FOXP3","STAT3","STAT5B",
    
    # --- 4. Differentiation Hierarchy ---
    "TCF7","LEF1","MYB",            # Stem-like Progenitor
    "E2F1","FOXM1",                 # Proliferation (Cycling)
    "PRDM1",                        # Terminal Effector
    
    # --- 5. KEGG ALIGNED CATEGORIES (NEW) ---
    "EOMES","TBX21","RUNX3",        # NK-like Cytotoxicity (Cluster 1,9)
    "RFX5","CREB1",                 # Antigen Presentation / MHC-II (Cluster 0)
    "KLF4","ETS1","SMAD3",          # Migration / Cell Adhesion (CAMs)
    
    # --- 6. Microenvironment ---
    "STAT1","STAT2","IRF9",         # IFN Response
    "HIF1A","SREBF1",               # Metabolism
    
    # --- 7. Normal Baseline / Tumor Suppressors ---
    "FOXO1","FOXO4","ZEB1","BACH2",
    "TCF3","BCL11A","NEUROD1","MEF2B","PBX2","IRF3","BCL6"
  ),
  Condition = c(
    rep("Malignant", 4),   # Oncogenic
    rep("Malignant", 3),   # Stress
    
    rep("Malignant", 3),   # AP-1
    rep("Malignant", 2),   # NFAT
    rep("Malignant", 5),   # NF-kB
    
    rep("Malignant", 6),   # Th2 / JAK-STAT Core
    
    rep("Malignant", 3),   # Stem-like Progenitor
    rep("Malignant", 2),   # Proliferation
    rep("Malignant", 1),   # Terminal Effector
    
    rep("Malignant", 3),   # NK-like Cytotoxicity
    rep("Malignant", 2),   # Antigen Presentation
    rep("Malignant", 3),   # Migration / Adhesion
    
    rep("Malignant", 3),   # IFN
    rep("Malignant", 2),   # Metabolism
    
    rep("Normal", 4),      # Tumor Suppressor
    rep("Normal", 7)       # Homeostasis
  ),
  Function = c(
    rep("Oncogenic", 4),
    rep("Stress Response", 3),
    
    rep("AP-1 Signaling", 3),
    rep("NFAT Signaling", 2),
    rep("Inflammatory/NF-kB", 5),
    
    rep("Th2 / JAK-STAT Core", 6),
    
    rep("Stem-like Progenitor", 3),
    rep("Proliferation", 2),
    rep("Terminal Effector", 1),
    
    rep("NK-like Cytotoxicity", 3),   # KEGG aligned
    rep("Antigen Presentation", 2),   # KEGG aligned
    rep("Migration / Adhesion", 3),   # KEGG aligned
    
    rep("IFN Response", 3),
    rep("Metabolism", 2),
    
    rep("Tumor Suppressor", 4),
    rep("Normal Homeostasis", 7)
  ),
  stringsAsFactors = FALSE
)

# Lock in the precise order of the blocks from top to bottom
desired_order <- c(
  "Oncogenic", 
  "Stress Response",
  "AP-1 Signaling", 
  "NFAT Signaling", 
  "Inflammatory/NF-kB",
  "Th2 / JAK-STAT Core",
  "Stem-like Progenitor", 
  "Proliferation",        
  "Terminal Effector",
  "NK-like Cytotoxicity",   # Placed here to show effector state
  "Antigen Presentation",   # Directly links to MHC-II high cluster
  "Migration / Adhesion",   # Links to CAMs KEGG pathway
  "IFN Response",
  "Metabolism",
  "Tumor Suppressor",
  "Normal Homeostasis"
)

tf_meta$Function <- factor(tf_meta$Function, levels = desired_order)

# ============================================
# 2. Extract TF activity matrix from Seurat
# ============================================
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

# Filter for TFs present in Seurat
available_tfs <- intersect(tf_meta$TF, rownames(mat_use))
mat_use <- mat_use[available_tfs, , drop = FALSE]

# ============================================
# 3. Average per cluster & z-score
# ============================================
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

# Align tf_meta order to match mat_use precisely
row_order_idx <- match(rownames(avg_mat_z), tf_meta$TF)
tf_meta_filtered <- tf_meta[row_order_idx, ]
avg_mat_z <- avg_mat_z[order(tf_meta_filtered$Function), ]
tf_meta_filtered <- tf_meta_filtered[order(tf_meta_filtered$Function), ]

# ============================================
# 4. Column & row annotations
# ============================================
cluster_status <- ifelse(colnames(avg_mat_z) %in% c("3","10"), "Normal CD4 T cells", "Malignant CD4 T cells")
ha_col <- HeatmapAnnotation(
  Cell_State = cluster_status,
  col = list(Cell_State = c("Normal CD4 T cells" = "#4DAF4A", "Malignant CD4 T cells" = "#E41A1C")),
  annotation_name_side = "left"
)

# Row annotation colors
function_colors <- c(
  "Oncogenic" = "#808080",
  "Stress Response" = "#FFD700",
  
  "AP-1 Signaling" = "#FF8C00",
  "NFAT Signaling" = "#FF4500",
  "Inflammatory/NF-kB" = "#E31A1C",
  
  "Th2 / JAK-STAT Core" = "#984EA3",
  
  "Stem-like Progenitor" = "#FF1493",
  "Proliferation" = "#1E90FF",
  "Terminal Effector" = "#800080",
  
  # New KEGG categories
  "NK-like Cytotoxicity" = "#F781BF",   # Pink
  "Antigen Presentation" = "#66CDAA",   # Medium Aquamarine
  "Migration / Adhesion" = "#8A2BE2",   # Blue Violet
  
  "IFN Response" = "#00CED1",
  "Metabolism" = "#A65628",
  
  "Tumor Suppressor" = "#377EB8",
  "Normal Homeostasis" = "#4DAF4A"
)

ha_row <- rowAnnotation(
  Condition = tf_meta_filtered$Condition,
  Function = tf_meta_filtered$Function,
  col = list(
    Condition = c("Normal"="#4DAF4A","Malignant"="#E41A1C"),
    Function = function_colors
  ),
  annotation_name_side = "bottom"
)

# ============================================
# 5. Heatmap colors & plotting
# ============================================
col_fun <- circlize::colorRamp2(c(-3,0,3), c("#313695","white","#A50026"))

ht <- Heatmap(
  avg_mat_z,
  name = "TF activity (z)",
  col = col_fun,
  top_annotation = ha_col,
  left_annotation = ha_row,
  column_split = cluster_status,
  row_split = tf_meta_filtered$Function,
  cluster_row_slices = FALSE,    
  cluster_rows = FALSE,          
  cluster_columns = TRUE,
  show_row_dend = FALSE,
  show_column_dend = TRUE,
  row_names_gp = gpar(fontsize = 10),
  column_names_gp = gpar(fontsize = 10),
  column_title = "Functional TF Modules in Sézary Syndrome",
  row_title_rot = 0,
  row_title_gp = gpar(fontsize = 8, fontface = "bold"),
  heatmap_legend_param = list(direction = "vertical")
)

# ============================================
# 6. Save plots
# ============================================
pdf("Output_Figures/Figure_TF_Heatmap_KEGG.pdf", width=14, height=14)
draw(ht, merge_legend=TRUE)
dev.off()

null device 
          1 

draw(ht, merge_legend=TRUE)

12 Figure Malignant complex heatmap chunk (with 44-47 TFs, literature-based panel)

# ============================================
# LIBRARIES
# ============================================
library(ComplexHeatmap)
library(circlize)
library(Matrix)
library(grid)
library(SeuratObject)

# ============================================
# 1. Define TF panel metadata (1:1 UMAP Aligned)
# ============================================
tf_meta <- data.frame(
  TF = c(
    # --- Baseline Malignancy --- 
    "MYC","E2F4","TWIST1","IRF4",
    
    # --- Clusters 2 & 6: Th2-like Core --- 
    "GATA3","STAT6","BATF","FOXP3","STAT3","STAT5B",
    
    # --- Clusters 11 & 12: Pro-inflammatory & Stress --- 
    "JUNB","FOS","FOSL1",           # AP-1
    "RELA","REL","NFKB1","IRF1","NCOA2", # NF-kB
    "NFATC1","NFATC2",              # TCR/NFAT
    "HSF1","NFE2L2","SREBF2",       # Stress
    
    # --- Cluster 4: Inflammatory-Migratory ---
    "KLF4","ETS1","SMAD3",
    
    # --- Cluster 5: Stem-like --- 
    "TCF7","LEF1","MYB",            
    
    # --- Cluster 7: Cycling (G2/M) --- 
    "E2F1","FOXM1",                 
    
    # --- Clusters 1 & 9: NK-like / Cytotoxic --- 
    "EOMES","TBX21","RUNX3","PRDM1",        
    
    # --- Cluster 0: MHC-II High --- 
    "RFX5","CREB1",                 
    
    # --- Cluster 13: IFN Stimulated --- 
    "STAT1","STAT2","IRF9",         
    
    # --- Cluster 8: Glycolytic/Metabolic --- 
    "HIF1A","SREBF1",               
    
    # --- Clusters 3 & 10: Normal CD4 T --- 
    "FOXO1","FOXO4","ZEB1","BACH2",
    "TCF3","BCL11A","NEUROD1","MEF2B","PBX2","IRF3","BCL6"
  ),
  Condition = c(
    rep("Malignant", 4),   # Oncogenic
    rep("Malignant", 6),   # Th2
    rep("Malignant", 13),  # Pro-inflammatory (AP1, NFkB, NFAT, Stress)
    rep("Malignant", 3),   # Migratory
    rep("Malignant", 3),   # Stem-like
    rep("Malignant", 2),   # Cycling
    rep("Malignant", 4),   # NK/Cytotoxic
    rep("Malignant", 2),   # MHC-II
    rep("Malignant", 3),   # IFN
    rep("Malignant", 2),   # Glycolytic
    rep("Normal", 11)      # Normal
  ),
  Function = c(
    rep("Oncogenic Core", 4),
    rep("Th2-like Core (Cl. 2, 6)", 6),
    rep("Pro-inflammatory (Cl. 11, 12)", 13),
    rep("Inflammatory-Migratory (Cl. 4)", 3),
    rep("Stem-like (Cl. 5)", 3),
    rep("Cycling G2/M (Cl. 7)", 2),
    rep("NK-like Cytotoxic (Cl. 1, 9)", 4),   
    rep("MHC-II High (Cl. 0)", 2),   
    rep("IFN Stimulated (Cl. 13)", 3),
    rep("Glycolytic/Metabolic (Cl. 8)", 2),
    rep("Normal Homeostasis (Cl. 3, 10)", 11)
  ),
  stringsAsFactors = FALSE
)

# Lock in the precise order of the blocks
desired_order <- c(
  "Oncogenic Core", 
  "Th2-like Core (Cl. 2, 6)",
  "Pro-inflammatory (Cl. 11, 12)", 
  "Inflammatory-Migratory (Cl. 4)", 
  "Stem-like (Cl. 5)", 
  "Cycling G2/M (Cl. 7)",        
  "NK-like Cytotoxic (Cl. 1, 9)",   
  "MHC-II High (Cl. 0)",   
  "IFN Stimulated (Cl. 13)",
  "Glycolytic/Metabolic (Cl. 8)",
  "Normal Homeostasis (Cl. 3, 10)"
)

tf_meta$Function <- factor(tf_meta$Function, levels = desired_order)

# ============================================
# 2. Extract TF activity matrix from Seurat
# ============================================
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

available_tfs <- intersect(tf_meta$TF, rownames(mat_use))
mat_use <- mat_use[available_tfs, , drop = FALSE]

# ============================================
# 3. Average per cluster & z-score
# ============================================
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

# Align tf_meta order to match mat_use precisely
row_order_idx <- match(rownames(avg_mat_z), tf_meta$TF)
tf_meta_filtered <- tf_meta[row_order_idx, ]
avg_mat_z <- avg_mat_z[order(tf_meta_filtered$Function), ]
tf_meta_filtered <- tf_meta_filtered[order(tf_meta_filtered$Function), ]

# ============================================
# 4. Column & row annotations
# ============================================
cluster_status <- ifelse(colnames(avg_mat_z) %in% c("3","10"), "Normal CD4 T cells", "Malignant CD4 T cells")
ha_col <- HeatmapAnnotation(Cell_State = cluster_status, col = list(Cell_State = c("Normal CD4 T cells" = "#4DAF4A", "Malignant CD4 T cells" = "#E41A1C")), annotation_name_side = "left")

# Harmonized color palette
function_colors <- c(
  "Oncogenic Core" = "#808080",
  "Th2-like Core (Cl. 2, 6)" = "#984EA3",
  "Pro-inflammatory (Cl. 11, 12)" = "#E31A1C",
  "Inflammatory-Migratory (Cl. 4)" = "#8A2BE2",
  "Stem-like (Cl. 5)" = "#FF1493",
  "Cycling G2/M (Cl. 7)" = "#1E90FF",
  "NK-like Cytotoxic (Cl. 1, 9)" = "#F781BF",   
  "MHC-II High (Cl. 0)" = "#66CDAA",   
  "IFN Stimulated (Cl. 13)" = "#00CED1",
  "Glycolytic/Metabolic (Cl. 8)" = "#A65628",
  "Normal Homeostasis (Cl. 3, 10)" = "#4DAF4A"
)

ha_row <- rowAnnotation(Condition = tf_meta_filtered$Condition, Function = tf_meta_filtered$Function, col = list(Condition = c("Normal"="#4DAF4A","Malignant"="#E41A1C"), Function = function_colors), annotation_name_side = "bottom")

# ============================================
# 5. Heatmap colors & plotting
# ============================================
col_fun <- circlize::colorRamp2(c(-3,0,3), c("#313695","white","#A50026"))

ht <- Heatmap(
  avg_mat_z, 
  name = "TF activity (z)", 
  col = col_fun, 
  top_annotation = ha_col, 
  left_annotation = ha_row, 
  column_split = cluster_status, 
  row_split = tf_meta_filtered$Function, 
  cluster_row_slices = FALSE, 
  cluster_rows = FALSE, 
  cluster_columns = TRUE, 
  show_row_dend = FALSE, 
  show_column_dend = TRUE, 
  row_names_gp = gpar(fontsize = 10), 
  column_names_gp = gpar(fontsize = 10), 
  column_title = "Regulatory Drivers of Sézary UMAP Cell States", 
  row_title_rot = 0, 
  row_title_gp = gpar(fontsize = 9, fontface = "bold"), 
  heatmap_legend_param = list(direction = "vertical")
)

# ============================================
# 6. Save plots
# ============================================
pdf("Output_Figures/Figure_TF_Heatmap_UMAP_Final.pdf", width=15, height=13)
draw(ht, merge_legend=TRUE)
dev.off()

null device 
          1 

draw(ht, merge_legend=TRUE)

13 TEST

# ============================================
# LIBRARIES
# ============================================
library(ComplexHeatmap)
library(circlize)
library(Matrix)
library(grid)
library(SeuratObject)

# ============================================
# 1. Define Literature-Validated TF panel (Exhaustive)
# ============================================
tf_meta_lit <- data.frame(
  TF = c(
    # Oncogenic 
    "MYC", "TWIST1", "IRF4",
    
    # Th2 Core
    "GATA3", "BATF", "FOXP3", "STAT3", "STAT5B",
    
    # Hyperactive TCR / Inflammatory
    "JUNB", "NFATC1", "NFATC2", "RELA", "NFKB1",
    
    # Canonical Tumor Suppressors
    "ZEB1", "BACH2", "FOXO1", "FOXO3"
  ),
  Function = c(
    rep("Oncogenic ", 3),
    rep("Th2 Core", 5),
    rep("Hyperactive TCR / Inflammatory", 5),
    rep("Canonical Tumor Suppressors", 4)
  ),
  stringsAsFactors = FALSE
)

# Lock in the order top-to-bottom
tf_meta_lit$Function <- factor(tf_meta_lit$Function, 
                               levels = c("Oncogenic ", 
                                         "Th2 Core", 
                                         "Hyperactive TCR / Inflammatory", 
                                         "Canonical Tumor Suppressors"))

# ============================================
# 2. Extract TF activity matrix
# ============================================
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

# Filter for available TFs
available_tfs <- intersect(tf_meta_lit$TF, rownames(mat_use))
mat_use <- mat_use[available_tfs, , drop = FALSE]

# ============================================
# 3. Average per cluster & z-score
# ============================================
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

# Align to our defined literature list order
row_order_idx <- match(rownames(avg_mat_z), tf_meta_lit$TF)
tf_meta_lit_filtered <- tf_meta_lit[row_order_idx, ]
avg_mat_z <- avg_mat_z[order(tf_meta_lit_filtered$Function), ]
tf_meta_lit_filtered <- tf_meta_lit_filtered[order(tf_meta_lit_filtered$Function), ]

# ============================================
# 4. Column & row annotations
# ============================================
cluster_status <- ifelse(colnames(avg_mat_z) %in% c("3","10"), "Normal CD4 T cells", "Malignant CD4 T cells")
ha_col <- HeatmapAnnotation(
  Cell_State = cluster_status, 
  col = list(Cell_State = c("Normal CD4 T cells" = "#4DAF4A", "Malignant CD4 T cells" = "#E41A1C")), 
  annotation_name_side = "left"
)

# Colors matching the literature categories
role_colors <- c(
  "Oncogenic " = "#808080",  # Gray
  "Th2 Core" = "#984EA3",                   # Purple
  "Hyperactive TCR / Inflammatory" = "#E31A1C",    # Red
  "Canonical Tumor Suppressors" = "#377EB8"        # Blue
)

ha_row <- rowAnnotation(
  Function = tf_meta_lit_filtered$Function, 
  col = list(Function = role_colors), 
  annotation_name_side = "bottom"
)

# ============================================
# 5. Heatmap colors & plotting
# ============================================
col_fun <- circlize::colorRamp2(c(-3,0,3), c("#313695","white","#A50026"))

ht_lit <- Heatmap(
  avg_mat_z, 
  name = "TF activity (z)", 
  col = col_fun, 
  top_annotation = ha_col, 
  left_annotation = ha_row, 
  column_split = cluster_status, 
  row_split = tf_meta_lit_filtered$Function, 
  cluster_row_slices = FALSE, 
  cluster_rows = FALSE, 
  cluster_columns = TRUE, 
  show_row_dend = FALSE, 
  show_column_dend = TRUE, 
  row_names_gp = gpar(fontsize = 12, fontface = "bold"), 
  column_names_gp = gpar(fontsize = 12), 
  column_title = "Literature-Validated Sézary Syndrome Regulators", 
  row_title_rot = 0, 
  row_title_gp = gpar(fontsize = 10, fontface = "bold"), 
  heatmap_legend_param = list(direction = "vertical")
)

# ============================================
# 6. Save plots
# ============================================
pdf("Output_Figures/Figure_TF_Heatmap_Literature_Validated.pdf", width=12, height=8)
draw(ht_lit, merge_legend=TRUE)
dev.off()

null device 
          1 

png("Output_Figures/Figure_TF_Heatmap_Literature_Validated.png", width=12*300, height=8*300, res=300)

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

null device 
          1 

draw(ht_lit, merge_legend=TRUE)

14 Define Th1/Th2/Th17/Th22/Treg Master Regulator Panel

# ============================================
# LIBRARIES
# ============================================
library(ComplexHeatmap)
library(circlize)
library(Matrix)
library(grid)
library(SeuratObject)

# ============================================
# 1. Define Th1/Th2/Th17/Th22/Treg Master Regulator Panel
# ============================================
tf_meta_thelper <- data.frame(
  TF = c(
    # Th1 Master Regulators
    "TBX21", "STAT1", "STAT4", "IRF1",
    
    # Th2 Master Regulators  
    "GATA3", "STAT6", "BATF", "IRF4",
    
    # Th17 / Th22 Master Regulators
    "RORC", "STAT3", "AHR", "MAF",
    
    # Treg Master Regulators
    "FOXP3", "FOXO1", "CTLA4"  # CTLA4 regulon as Treg proxy
  ),
  Function = c(
    rep("Th1", 4),
    rep("Th2", 4),
    rep("Th17/Th22", 4),
    rep("Treg", 3)
  ),
  stringsAsFactors = FALSE
)

# Lock in canonical order: Th1 → Th2 → Th17/Th22 → Treg
tf_meta_thelper$Function <- factor(tf_meta_thelper$Function, 
                                   levels = c("Th1", "Th2", "Th17/Th22", "Treg"))

# ============================================
# 2. Extract TF activity matrix
# ============================================
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

# Filter for available TFs
available_tfs <- intersect(tf_meta_thelper$TF, rownames(mat_use))
mat_use <- mat_use[available_tfs, , drop = FALSE]

# ============================================
# 3. Average per cluster & z-score
# ============================================
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

# Align to our defined helper T panel order
row_order_idx <- match(rownames(avg_mat_z), tf_meta_thelper$TF)
tf_meta_thelper_filtered <- tf_meta_thelper[row_order_idx, ]
avg_mat_z <- avg_mat_z[order(tf_meta_thelper_filtered$Function), ]
tf_meta_thelper_filtered <- tf_meta_thelper_filtered[order(tf_meta_thelper_filtered$Function), ]

# ============================================
# 4. Column & row annotations
# ============================================
cluster_status <- ifelse(colnames(avg_mat_z) %in% c("3","10"), "Normal CD4 T cells", "Malignant CD4 T cells")
ha_col <- HeatmapAnnotation(
  Cell_State = cluster_status, 
  col = list(Cell_State = c("Normal CD4 T cells" = "#4DAF4A", "Malignant CD4 T cells" = "#E41A1C")), 
  annotation_name_side = "left"
)

# Classical T helper color scheme
helper_colors <- c(
  "Th1" = "#E31A1C",      # Red
  "Th2" = "#1F78B4",      # Blue  
  "Th17/Th22" = "#FF7F00", # Orange
  "Treg" = "#33A02C"      # Green
)

ha_row <- rowAnnotation(
  Function = tf_meta_thelper_filtered$Function, 
  col = list(Function = helper_colors), 
  annotation_name_side = "bottom"
)

# ============================================
# 5. Heatmap colors & plotting
# ============================================
col_fun <- circlize::colorRamp2(c(-3,0,3), c("#313695","white","#A50026"))

ht_thelper <- Heatmap(
  avg_mat_z, 
  name = "TF activity (z)", 
  col = col_fun, 
  top_annotation = ha_col, 
  left_annotation = ha_row, 
  column_split = cluster_status, 
  row_split = tf_meta_thelper_filtered$Function, 
  cluster_row_slices = FALSE, 
  cluster_rows = FALSE, 
  cluster_columns = TRUE, 
  show_row_dend = FALSE, 
  show_column_dend = TRUE, 
  row_names_gp = gpar(fontsize = 11, fontface = "bold"), 
  column_names_gp = gpar(fontsize = 11), 
  column_title = "CD4+ T Helper Lineage Transcription Factors", 
  row_title_rot = 0, 
  row_title_gp = gpar(fontsize = 10, fontface = "bold"), 
  heatmap_legend_param = list(direction = "vertical")
)

# ============================================
# 6. Save plots
# ============================================
pdf("Output_Figures/Figure_TF_Heatmap_THelper_Lineages.pdf", width=12, height=8)
draw(ht_thelper, merge_legend=TRUE)
dev.off()

null device 
          1 

png("Output_Figures/Figure_TF_Heatmap_THelper_Lineages.png", width=12*300, height=8*300, res=300)

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

null device 
          1 

draw(ht_thelper, merge_legend=TRUE)

Final Save

print("Analysis pipeline complete. All figures and objects saved in Output_Figures folder.")

[1] "Analysis pipeline complete. All figures and objects saved in Output_Figures folder."