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!")
}

2.2 SCpubr Heatmap Visualization-Heatmap of averaged scores

library(SCpubr)
# General heatmap (Top Variable TFs)
out <- SCpubr::do_TFActivityHeatmap(sample = seurat_obj,
                                 activities = activities)

print(out)

# 1. Save as PDF
pdf("Output_Figures/SCpubr_Heatmap_Default.pdf", width = 10, height = 8)
print(out) # ComplexHeatmap requires explicit print() inside pdf()
dev.off()

# 2. Save as PNG
png("Output_Figures/SCpubr_Heatmap_Default.png", width = 10 * 300, height = 8 * 300, res = 300)
print(out)
dev.off()


print(out)

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("Output_Figures/SCpubr_Heatmap_Scaled.png", width = 10 * 300, height = 8 * 300, res = 300)
print(out)
dev.off()

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("Output_Figures/SCpubr_Heatmap_Symmetric.png", width = 10 * 300, height = 8 * 300, res = 300)
print(out)
dev.off()


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("Output_Figures/SCpubr_Heatmap_Top40.png", width = 14 * 300, height = 6 * 300, res = 300)
print(out)
dev.off()

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("Output_Figures/Figure_3.16A_Global_TF_Heatmap_Top100.png", width = 32 * 300, height = 12 * 300, res = 300)
print(out)
dev.off()

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("Output_Figures/Figure_Top100.png", width = 32 * 300, height = 12 * 300, res = 300)
print(out)
dev.off()

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...")
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.")

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()

# 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()

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()

png("Output_Figures/Figure_3.16E_Differential_TF_Heatmap.png",
    width = 11 * 300, height = 9 * 300, res = 300)
draw(ht)
dev.off()
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()

png("Output_Figures/Figure_TF_Heatmap.png", width=12*300, height=10*300, res=300)
draw(ht, merge_legend=TRUE)
dev.off()

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()
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()

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()

draw(ht, merge_legend=TRUE)

13 TEST

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

# ============================================
# 1. Define TF panel metadata
# ============================================
tf_meta <- data.frame(
  TF = c(
    # --- 1. General Malignancy ---
    "MYC","E2F4","TWIST1","IRF4",

    # --- 2. TCR Signaling Triad ---
    "JUNB","FOS","FOSL1",                # AP-1
    "NFATC1","NFATC2",                   # NFAT
    "RELA","REL","NFKB1","IRF1","NCOA2", # NF-kB

    # --- 3. Stress Response ---
    "HSF1","NFE2L2","SREBF2",

    # --- 4. Th2 / JAK-STAT Core ---
    "GATA3","STAT6","BATF","FOXP3","STAT3","STAT5B",

    # --- 5. Differentiation Hierarchy ---
    "TCF7","LEF1","MYB",                 # Stem-like
    "E2F1","FOXM1",                      # Proliferation
    "PRDM1",                             # Terminal Effector

    # --- 6. KEGG Aligned ---
    "EOMES","TBX21","RUNX3",             # NK-like Cytotoxicity
    "RFX5","CREB1",                      # Antigen Presentation
    "KLF4","ETS1","SMAD3",              # Migration / Adhesion

    # --- 7. Microenvironment ---
    "STAT1","STAT2","IRF9",              # IFN Response
    "HIF1A","SREBF1",                    # Metabolism

    # --- 8. Tumor Suppressors / Normal ---
    "FOXO1","FOXO4","ZEB1","BACH2",
    "TCF3","BCL11A","NEUROD1","MEF2B","PBX2","IRF3","BCL6"
  ),
  Condition = c(
    rep("Malignant", 4),   # Oncogenic
    rep("Malignant", 3),   # AP-1
    rep("Malignant", 2),   # NFAT
    rep("Malignant", 5),   # NF-kB
    rep("Malignant", 3),   # Stress Response
    rep("Malignant", 6),   # Th2/JAK-STAT
    rep("Malignant", 3),   # Stem-like
    rep("Malignant", 2),   # Proliferation
    rep("Malignant", 1),   # Terminal Effector
    rep("Malignant", 3),   # NK-like
    rep("Malignant", 2),   # Antigen Presentation
    rep("Malignant", 3),   # Migration
    rep("Malignant", 3),   # IFN Response
    rep("Malignant", 2),   # Metabolism
    rep("Normal", 4),      # Tumor Suppressor
    rep("Normal", 7)       # Normal Homeostasis
  ),
  Function = c(
    rep("Oncogenic", 4),
    rep("AP-1 Signaling", 3),
    rep("NFAT Signaling", 2),
    rep("NF-\u03baB Signaling", 5),  # RENAMED from Inflammatory/NF-kB
    rep("Stress Response", 3),
    rep("Th2 / JAK-STAT Core", 6),
    rep("Stem-like Progenitor", 3),
    rep("Proliferation", 2),
    rep("Terminal Effector", 1),
    rep("NK-like Cytotoxicity", 3),
    rep("Antigen Presentation", 2),
    rep("Migration / Adhesion", 3),
    rep("IFN Response", 3),
    rep("Metabolism", 2),
    rep("Tumor Suppressor", 4),
    rep("Normal Homeostasis", 7)
  ),
  stringsAsFactors = FALSE
)

# Lock in order
desired_order <- c(
  "Oncogenic",
  "AP-1 Signaling",
  "NFAT Signaling",
  "NF-\u03baB Signaling",
  "Stress Response",
  "Th2 / JAK-STAT Core",
  "Stem-like Progenitor",
  "Proliferation",
  "Terminal Effector",
  "NK-like Cytotoxicity",
  "Antigen Presentation",
  "Migration / Adhesion",
  "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. Cluster status & labels
# ============================================
cluster_status <- ifelse(
  colnames(avg_mat_z) %in% c("3", "10"),
  "Normal CD4 T cells",
  "Malignant CD4 T cells"
)

# Cluster biological annotations
cluster_labels <- c(
  "0"  = "MHC-II High",
  "1"  = "NK-like",
  "2"  = "Th2-like",
  "3"  = "Naive/TCM",
  "4"  = "Inflammatory",
  "5"  = "Stem-like",
  "6"  = "Th2 Activated",
  "7"  = "Cycling",
  "8"  = "Glycolytic",
  "9"  = "Cytotoxic",
  "10" = "Central Memory",
  "11" = "Pro-inflammatory",
  "12" = "GZMB-high",
  "13" = "IFN Stimulated"
)
cluster_annotation <- cluster_labels[colnames(avg_mat_z)]

# ============================================
# 5. Column annotation (top)
# ============================================
ha_col <- HeatmapAnnotation(
  Cell_State   = cluster_status,
  Cluster_Type = cluster_annotation,
  col = list(
    Cell_State = c(
      "Normal CD4 T cells"    = "#4DAF4A",
      "Malignant CD4 T cells" = "#C00000"
    ),
    Cluster_Type = c(
      "MHC-II High"      = "#66CDAA",
      "NK-like"          = "#F781BF",
      "Th2-like"         = "#984EA3",
      "Naive/TCM"        = "#4DAF4A",
      "Inflammatory"     = "#8A2BE2",
      "Stem-like"        = "#FF1493",
      "Th2 Activated"    = "#DA70D6",
      "Cycling"          = "#1E90FF",
      "Glycolytic"       = "#A65628",
      "Cytotoxic"        = "#E41A1C",
      "Central Memory"   = "#32CD32",
      "Pro-inflammatory" = "#FF4500",
      "GZMB-high"        = "#B22222",
      "IFN Stimulated"   = "#00CED1"
    )
  ),
  annotation_name_side = "left",
  annotation_name_gp   = gpar(fontsize  = 8,
                               fontface = "bold"),
  simple_anno_size     = unit(4, "mm"),
  border               = TRUE
)

# ============================================
# 6. Row annotation colors
# ============================================
function_colors <- c(
  "Oncogenic"              = "#808080",
  # TCR Signaling Triad — related but distinct
  "AP-1 Signaling"         = "#FF8C00",
  "NFAT Signaling"         = "#FF4500",
  "NF-\u03baB Signaling"  = "#E31A1C",
  "Stress Response"        = "#FFD700",
  "Th2 / JAK-STAT Core"   = "#984EA3",
  "Stem-like Progenitor"   = "#FF1493",
  "Proliferation"          = "#1E90FF",
  "Terminal Effector"      = "#800080",
  "NK-like Cytotoxicity"   = "#F781BF",
  "Antigen Presentation"   = "#66CDAA",
  "Migration / Adhesion"   = "#8A2BE2",
  "IFN Response"           = "#00CED1",
  "Metabolism"             = "#A65628",
  "Tumor Suppressor"       = "#377EB8",
  "Normal Homeostasis"     = "#4DAF4A"
)

# ============================================
# 7. Row annotation (left)
# ============================================
ha_row <- rowAnnotation(
  Condition = tf_meta_filtered$Condition,
  Function  = tf_meta_filtered$Function,
  col = list(
    Condition = c(
      "Normal"    = "#4DAF4A",
      "Malignant" = "#C00000"
    ),
    Function = function_colors
  ),
  annotation_name_side = "bottom",
  annotation_name_gp   = gpar(fontsize  = 8,
                               fontface = "bold"),
  simple_anno_size     = unit(5, "mm"),
  border               = TRUE,
  gap                  = unit(1, "mm")
)

# ============================================
# 8. Heatmap color scale — 5 point gradient
# ============================================
col_fun <- circlize::colorRamp2(
  c(-3, -1.5, 0, 1.5, 3),
  c("#2166AC",   # Deep blue
    "#92C5DE",   # Light blue
    "#F7F7F7",   # Off-white
    "#F4A582",   # Light red
    "#B2182B")   # Deep red
)

# ============================================
# 9. Build heatmap
# ============================================
ht <- Heatmap(
  avg_mat_z,
  name = "TF activity\n(z-score)",
  col  = col_fun,

  # Annotations
  top_annotation  = ha_col,
  left_annotation = ha_row,

  # Splitting
  column_split = cluster_status,
  row_split    = tf_meta_filtered$Function,

  # Clustering
  cluster_row_slices = FALSE,
  cluster_rows       = FALSE,
  cluster_columns    = TRUE,
  show_row_dend      = FALSE,
  show_column_dend   = TRUE,
  column_dend_height = unit(15, "mm"),

  # Cell appearance
  rect_gp = gpar(col = "white", lwd = 0.8),

  # Row names
  row_names_gp   = gpar(fontsize  = 8.5,
                          fontface = "italic",
                          col      = "grey20"),
  row_names_side = "right",

  # Column names
  column_names_gp  = gpar(fontsize  = 10,
                            fontface = "bold",
                            col      = "grey10"),
  column_names_rot = 0,

  # Row title (functional block labels on left)
  row_title_rot  = 0,
  row_title_gp   = gpar(fontsize  = 8,
                          fontface = "bold",
                          col      = "grey10"),
  row_title_side = "left",

  # Gaps between blocks
  row_gap    = unit(2.5, "mm"),
  column_gap = unit(4,   "mm"),

  # Main title
  column_title    = "Functional TF Modules in Sézary Syndrome",
  column_title_gp = gpar(fontsize  = 13,
                           fontface = "bold",
                           col      = "grey5"),

  # Legend
  heatmap_legend_param = list(
    direction     = "vertical",
    title         = "TF activity\n(z-score)",
    title_gp      = gpar(fontsize  = 9,
                           fontface = "bold"),
    labels_gp     = gpar(fontsize  = 8),
    legend_height = unit(4, "cm"),
    border        = "grey50",
    at            = c(-3, -1.5, 0, 1.5, 3),
    labels        = c("-3", "-1.5", "0", "1.5", "3")
  )
)

# ============================================
# 10. Save
# ============================================
# PDF — vector quality for publication
pdf("Output_Figures/Figure_TF_Heatmap_Final.pdf",
    width       = 16,
    height      = 15,
    useDingbats = FALSE)
draw(ht,
     merge_legend           = TRUE,
     heatmap_legend_side    = "right",
     annotation_legend_side = "right",
     padding = unit(c(5, 5, 5, 5), "mm"))
dev.off()

null device 
          1 

# PNG — high resolution for presentations
png("Output_Figures/Figure_TF_Heatmap_Final.png",
    width  = 16 * 300,
    height = 15 * 300,
    res    = 300)

draw(ht,
     merge_legend           = TRUE,
     heatmap_legend_side    = "right",
     annotation_legend_side = "right",
     padding = unit(c(5, 5, 5, 5), "mm"))
dev.off()

null device 
          1 

# Draw in notebook
draw(ht,
     merge_legend           = TRUE,
     heatmap_legend_side    = "right",
     annotation_legend_side = "right",
     padding = unit(c(5, 5, 5, 5), "mm"))

14 Test2

15 Test3

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

# ============================================
# 1. Define TF panel metadata
# ============================================
tf_meta <- data.frame(
  TF = c(
    # --- 1. General Malignancy ---
    "MYC","E2F4","TWIST1","IRF4",

    # --- 2. TCR Signaling Triad ---
    "JUNB","FOS","FOSL1",                # AP-1
    "NFATC1","NFATC2",                   # NFAT
    "RELA","REL","NFKB1","IRF1","NCOA2", # NF-kB

    # --- 3. Stress Response ---
    "HSF1","NFE2L2","SREBF2",

    # --- 4. Th2 / JAK-STAT Core ---
    "GATA3","STAT6","BATF","FOXP3","STAT3","STAT5B",

    # --- 5. Differentiation Hierarchy ---
    "TCF7","LEF1","MYB",                 # Stem-like
    "E2F1","FOXM1",                      # Proliferation
    "PRDM1",                             # Terminal Effector

    # --- 6. KEGG Aligned ---
    "EOMES","TBX21","RUNX3",             # NK-like Cytotoxicity
    "RFX5","CREB1",                      # Antigen Presentation
    "KLF4","ETS1","SMAD3",               # Migration / Adhesion

    # --- 7. Microenvironment ---
    "STAT1","STAT2","IRF9",              # IFN Response
    "HIF1A","SREBF1",                    # Metabolism

    # --- 8. Tumor Suppressors / Normal ---
    "FOXO1","FOXO4","ZEB1","BACH2",
    "TCF3","BCL11A","NEUROD1","MEF2B","PBX2","IRF3","BCL6"
  ),
  Condition = c(
    rep("Malignant", 4),   # Oncogenic
    rep("Malignant", 3),   # AP-1
    rep("Malignant", 2),   # NFAT
    rep("Malignant", 5),   # NF-kB
    rep("Malignant", 3),   # Stress Response
    rep("Malignant", 6),   # Th2/JAK-STAT
    rep("Malignant", 3),   # Stem-like
    rep("Malignant", 2),   # Proliferation
    rep("Malignant", 1),   # Terminal Effector
    rep("Malignant", 3),   # NK-like
    rep("Malignant", 2),   # Antigen Presentation
    rep("Malignant", 3),   # Migration
    rep("Malignant", 3),   # IFN Response
    rep("Malignant", 2),   # Metabolism
    rep("Normal", 4),      # Tumor Suppressor
    rep("Normal", 7)       # Normal Homeostasis
  ),
  Function = c(
    rep("Oncogenic", 4),
    rep("AP-1 Signaling", 3),
    rep("NFAT Signaling", 2),
    rep("NF-\u03baB Signaling", 5),
    rep("Stress Response", 3),
    rep("Th2 / JAK-STAT Core", 6),
    rep("Stem-like Progenitor", 3),
    rep("Proliferation", 2),
    rep("Terminal Effector", 1),
    rep("NK-like Cytotoxicity", 3),
    rep("Antigen Presentation", 2),
    rep("Migration / Adhesion", 3),
    rep("IFN Response", 3),
    rep("Metabolism", 2),
    rep("Tumor Suppressor", 4),
    rep("Normal Homeostasis", 7)
  ),
  stringsAsFactors = FALSE
)

# Lock in order
desired_order <- c(
  "Oncogenic",
  "AP-1 Signaling",
  "NFAT Signaling",
  "NF-\u03baB Signaling",
  "Stress Response",
  "Th2 / JAK-STAT Core",
  "Stem-like Progenitor",
  "Proliferation",
  "Terminal Effector",
  "NK-like Cytotoxicity",
  "Antigen Presentation",
  "Migration / Adhesion",
  "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 orig.ident & z-score
# ============================================
orig_ident <- as.factor(seurat_obj$orig.ident)

avg_mat <- sapply(levels(orig_ident), function(id) {
  Matrix::rowMeans(mat_use[, orig_ident == id,
                            drop = FALSE])
})
colnames(avg_mat) <- levels(orig_ident)

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. Define sample metadata
# FIXED: names match actual orig.ident levels
# ============================================
sample_to_patient <- c(
  "L1"       = "Patient 1",
  "L2"       = "Patient 1",
  "L3"       = "Patient 2",
  "L4"       = "Patient 2",
  "L5"       = "Patient 3",
  "L6"       = "Patient 3",
  "L7"       = "Patient 3",
  "CD4T_lab" = "Normal",
  "CD4T_10x" = "Normal"
)

sample_to_type <- c(
  "L1"       = "Malignant",
  "L2"       = "Malignant",
  "L3"       = "Malignant",
  "L4"       = "Malignant",
  "L5"       = "Malignant",
  "L6"       = "Malignant",
  "L7"       = "Malignant",
  "CD4T_lab" = "Normal",
  "CD4T_10x" = "Normal"
)

# Extract annotations
patient_annotation  <- sample_to_patient[colnames(avg_mat_z)]
celltype_annotation <- sample_to_type[colnames(avg_mat_z)]

# Verify no NAs — should print character(0) for both
cat("Patient NAs:", colnames(avg_mat_z)[is.na(patient_annotation)], "\n")

Patient NAs:  

cat("Type NAs:",    colnames(avg_mat_z)[is.na(celltype_annotation)], "\n")

Type NAs:  

# ============================================
# 5. Column annotation (top)
# ============================================
ha_col <- HeatmapAnnotation(
  Cell_Type = celltype_annotation,
  Patient   = patient_annotation,
  col = list(
    Cell_Type = c(
      "Normal"    = "#4DAF4A",
      "Malignant" = "#C00000"
    ),
    Patient = c(
      "Patient 1" = "#E41A1C",
      "Patient 2" = "#377EB8",
      "Patient 3" = "#FF7F00",
      "Normal"    = "#4DAF4A"
    )
  ),
  annotation_name_side = "left",
  annotation_name_gp   = gpar(fontsize  = 9,
                               fontface = "bold",
                               col      = "black"),
  simple_anno_size     = unit(5, "mm"),
  border               = TRUE
)

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

# ============================================
# 7. Row annotation (left)
# ============================================
ha_row <- rowAnnotation(
  Condition = tf_meta_filtered$Condition,
  Function  = tf_meta_filtered$Function,
  col = list(
    Condition = c(
      "Normal"    = "#4DAF4A",
      "Malignant" = "#C00000"
    ),
    Function = function_colors
  ),
  annotation_name_side = "bottom",
  annotation_name_gp   = gpar(fontsize  = 9,
                               fontface = "bold",
                               col      = "black"),
  simple_anno_size     = unit(5, "mm"),
  border               = TRUE,
  gap                  = unit(1, "mm")
)

# ============================================
# 8. Heatmap color scale — 5 point gradient
# ============================================
col_fun <- circlize::colorRamp2(
  c(-3, -1.5, 0, 1.5, 3),
  c("#2166AC",   # Deep blue
    "#92C5DE",   # Light blue
    "#F7F7F7",   # Off-white
    "#F4A582",   # Light red
    "#B2182B")   # Deep red
)

# ============================================
# 9. Build heatmap
# ============================================
ht <- Heatmap(
  avg_mat_z,
  name = "TF activity\n(z-score)",
  col  = col_fun,

  # Annotations
  top_annotation  = ha_col,
  left_annotation = ha_row,

  # Splitting — columns by patient, rows by function
  column_split = patient_annotation,
  row_split    = tf_meta_filtered$Function,

  # Clustering
  cluster_row_slices = FALSE,
  cluster_rows       = FALSE,
  cluster_columns    = TRUE,
  show_row_dend      = FALSE,
  show_column_dend   = TRUE,
  column_dend_height = unit(15, "mm"),

  # Cell appearance
  rect_gp = gpar(col = "white", lwd = 0.8),

  # Row names — BOLD for slide visibility
  row_names_gp   = gpar(fontsize  = 12,
                          fontface = "bold",
                          col      = "black"),
  row_names_side = "right",

  # Column names — BOLD for slide visibility
  column_names_gp  = gpar(fontsize  = 12,
                            fontface = "bold",
                            col      = "black"),
  column_names_rot = 45,

  # Row title
  row_title_rot  = 0,
  row_title_gp   = gpar(fontsize  = 11,
                          fontface = "bold",
                          col      = "black"),
  row_title_side = "left",

  # Gaps
  row_gap    = unit(2.5, "mm"),
  column_gap = unit(6,   "mm"),

  # Main title — defined only once
  column_title    = "Functional TF Modules per Cell Line — Sézary Syndrome",
  column_title_gp = gpar(fontsize  = 14,
                           fontface = "bold",
                           col      = "black"),

  # Legend — at/labels removed to avoid version errors
  heatmap_legend_param = list(
    direction     = "vertical",
    title         = "TF activity\n(z-score)",
    title_gp      = gpar(fontsize  = 10,
                           fontface = "bold",
                           col      = "black"),
    labels_gp     = gpar(fontsize  = 9,
                           col      = "black"),
    legend_height = unit(4, "cm"),
    border        = "grey50"
  )
)

# ============================================
# 10. Save
# ============================================
# PDF — vector quality for publication
pdf("Output_Figures/Figure_TF_Heatmap_OrigIdent_Final.pdf",
    width       = 18,
    height      = 16,
    useDingbats = FALSE)
draw(ht,
     merge_legend           = TRUE,
     heatmap_legend_side    = "right",
     annotation_legend_side = "right",
     padding = unit(c(5, 5, 5, 5), "mm"))
dev.off()

null device 
          1 

# PNG — high resolution for presentations
png("Output_Figures/Figure_TF_Heatmap_OrigIdent_Final.png",
    width  = 18 * 300,
    height = 16 * 300,
    res    = 300)

draw(ht,
     merge_legend           = TRUE,
     heatmap_legend_side    = "right",
     annotation_legend_side = "right",
     padding = unit(c(5, 5, 5, 5), "mm"))
dev.off()

null device 
          1 

# Draw in notebook
draw(ht,
     merge_legend           = TRUE,
     heatmap_legend_side    = "right",
     annotation_legend_side = "right",
     padding = unit(c(5, 5, 5, 5), "mm"))

16 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()

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()

draw(ht_thelper, merge_legend=TRUE)

Final Save

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