Visualization of TF activity

1 Load necessary libraries

# IMPORTANT: Replace this with the actual path and command to load your object
SS <- readRDS("../../1-Seurat_RDS_OBJECT_FINAL/Seurat_Object_With_TF_Activity.rds")
# Placeholder:
if (!exists("SS")) {
  stop("The 'All_samples_Merged' Seurat object is not loaded. Please load your data.")
}

# Check the object structure
print(SS)

An object of class Seurat 
63171 features across 49305 samples within 7 assays 
Active assay: RNA (36601 features, 0 variable features)
 2 layers present: data, counts
 6 other assays present: ADT, prediction.score.celltype.l1, prediction.score.celltype.l2, prediction.score.celltype.l3, SCT, dorothea
 5 dimensional reductions calculated: integrated_dr, ref.umap, pca, umap, harmony

2 Adapting Your Dorothea Assay for SeuratExtend Visualizations

2.1 Rename or Copy Dorothea Assay to “TF”

library(SeuratExtend)
library(Seurat)

# Option 1: Copy dorothea assay to "TF" (recommended - keeps both)
SS[["TF"]] <- SS[["dorothea"]]
DefaultAssay(SS) <- "TF"

3 WaterfallPlot - Top 200 TFs (Main Interest)

library(SeuratExtend)
library(Seurat)

# Copy dorothea to TF assay
SS[["TF"]] <- SS[["dorothea"]]
DefaultAssay(SS) <- "TF"

# Correct way to get data in Seurat v5
tf_data <- GetAssayData(SS, layer = "data", assay = "TF")  # Changed slot to layer

# Calculate variance for top 200 TFs
tf_var <- apply(tf_data, 1, var)
top200_tfs <- names(sort(tf_var, decreasing = TRUE)[1:200])

cat("Top 200 TFs selected:\n")

Top 200 TFs selected:

print(head(top200_tfs, 20))

 [1] "RFX5"   "MYC"    "E2F4"   "FOXM1"  "MYCN"   "E2F1"   "STAT4"  "TBX21"  "E2F2"   "STAT6"  "NFKB1"  "HIF1A"  "ATF6"   "TFDP1"  "LYL1"   "SREBF2"
[17] "TCF7L2" "MEF2C"  "THAP11" "NFYA"  

# WaterfallPlot: Malignant vs Normal
Idents(SS) <- SS$Condition  # Set your condition column

WaterfallPlot(
  SS,
  features = top200_tfs,  # Top 200 TFs
  ident.1 = "MalignantCD4T",
  ident.2 = "NormalCD4T",
  exp.transform = FALSE,  # Important: dorothea scores are already normalized
  top.n = 50  # Show top 50 in the plot
)

NA
NA

3.1 Rename or Copy Dorothea Assay to “TF”

library(SeuratExtend)

DefaultAssay(SS) <- "TF"

# This is the CORRECT syntax - no transpose parameter exists
tf_zscore1 <- CalcStats(
  SS, 
  features = rownames(SS[["TF"]]),
  group.by = "orig.ident",
  method = "zscore",
  order = "p",
  n = 4  # Top 4 TFs per cell line
)

# This is the CORRECT syntax - no transpose parameter exists
tf_zscore2 <- CalcStats(
  SS, 
  features = rownames(SS[["TF"]]),
  group.by = "seurat_clusters",
  method = "zscore",
  order = "p",
  n = 4  # Top 4 TFs per cell line
)

# Create heatmap - the output is already in the right format
Heatmap(tf_zscore1, lab_fill = "zscore", color_scheme = "RdBu")

# Create heatmap - the output is already in the right format
Heatmap(tf_zscore2, lab_fill = "zscore", color_scheme = "RdBu")

4 DimPlot2 - Show Key TFs on UMAP

# Select top 20 most differential TFs for visualization
top20_tfs <- top200_tfs[1:20]

DimPlot2(
  SS,
  features = top20_tfs[1:20],reduction = "umap",  # Show first 6 TFs
  cols = "RdYlBu",
  theme = NoAxes(),
  ncol = 3
)

# Compare gene expression vs TF activity
DimPlot2(
  SS,
  features = c("STAT3", "tf_STAT3", "NFKB1", "tf_NFKB1"), reduction = "umap",  # If you have gene expression too
  cols = list("tf_STAT3" = "D", "tf_NFKB1" = "D"),
  theme = NoAxes()
)

5 VlnPlot2 - Distribution of Top TFs

# Top 10 differential TFs
VlnPlot2(
  SS,
  features = top20_tfs[1:20],
  group.by = "orig.ident",
  ncol = 5
)

# Compare Malignant vs Normal
Idents(SS) <- SS$Condition
VlnPlot2(
  SS,
  features = top20_tfs[1:20],
  group.by = "Condition",
  ncol = 3
)

6 DotPlot2 - Top TFs Activity

DotPlot2(
  SS,
  features = top20_tfs,
  group.by = "orig.ident"
)

# Split by condition
DotPlot2(
  SS,
  features = top20_tfs,
  group.by = "orig.ident",
  split.by = "Condition"
)

7 Custom Waterfall Plot with ggplot2

library(Seurat)
library(ggplot2)
library(dplyr)
library(ggrepel)

# Step 1: Find differential TF activity
DefaultAssay(SS) <- "dorothea"
Idents(SS) <- SS$Condition

tf_markers <- FindMarkers(
  SS,
  ident.1 = "MalignantCD4T",
  ident.2 = "NormalCD4T",
  assay = "dorothea",
  logfc.threshold = 0,
  min.pct = 0,
  test.use = "wilcox"
)

# Step 2: Prepare data for waterfall plot
tf_markers$TF <- rownames(tf_markers)
tf_markers <- tf_markers %>%
  arrange(avg_log2FC) %>%
  mutate(rank = row_number(),
         significance = case_when(
           p_val_adj < 0.05 & avg_log2FC > 0 ~ "Up in Malignant",
           p_val_adj < 0.05 & avg_log2FC < 0 ~ "Up in Normal",
           TRUE ~ "Not Significant"
         ))

# Step 3: Select top N TFs for labeling
top_n <- 20
top_tfs <- tf_markers %>%
  arrange(desc(abs(avg_log2FC))) %>%
  head(top_n)

# Step 4: Create waterfall plot
ggplot(tf_markers, aes(x = rank, y = avg_log2FC, fill = significance)) +
  geom_bar(stat = "identity", width = 1) +
  scale_fill_manual(
    values = c("Up in Malignant" = "#D62728", 
               "Up in Normal" = "#2166AC", 
               "Not Significant" = "grey70"),
    name = "TF Activity"
  ) +
  geom_text_repel(
    data = top_tfs,
    aes(x = rank, y = avg_log2FC, label = TF),
    size = 3,
    max.overlaps = 20,
    box.padding = 0.5,
    segment.color = "grey50"
  ) +
  geom_hline(yintercept = 0, linetype = "solid", color = "black", linewidth = 0.5) +
  labs(
    title = "Differential TF Activity: Malignant vs Normal",
    subtitle = paste0("Top ", top_n, " TFs labeled"),
    x = "Ranked Transcription Factors",
    y = "log2 Fold Change (TF Activity)"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 14, face = "bold"),
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    legend.position = "right"
  )

NA
NA

8 Enhanced Version - Top 50 TFs Only

library(Seurat)
library(ggplot2)
library(dplyr)
library(ggrepel)

DefaultAssay(SS) <- "dorothea"
Idents(SS) <- SS$Condition

# Find differential TFs
tf_markers <- FindMarkers(
  SS,
  ident.1 = "MalignantCD4T",
  ident.2 = "NormalCD4T",
  assay = "dorothea",
  logfc.threshold = 0,
  min.pct = 0
)

# Select top 50 by absolute log2FC
tf_markers$TF <- rownames(tf_markers)
top50_tfs <- tf_markers %>%
  arrange(desc(abs(avg_log2FC))) %>%
  head(50)

# Prepare for waterfall plot
top50_tfs <- top50_tfs %>%
  arrange(avg_log2FC) %>%
  mutate(
    rank = row_number(),
    color_cat = case_when(
      p_val_adj < 0.05 & avg_log2FC > 0.5 ~ "Highly Up (Malignant)",
      p_val_adj < 0.05 & avg_log2FC > 0 ~ "Up (Malignant)",
      p_val_adj < 0.05 & avg_log2FC < -0.5 ~ "Highly Down (Normal)",
      p_val_adj < 0.05 & avg_log2FC < 0 ~ "Down (Normal)",
      TRUE ~ "Not Significant"
    )
  )

# Label top 20
top20_labels <- top50_tfs %>%
  arrange(desc(abs(avg_log2FC))) %>%
  head(20)

# Waterfall plot
ggplot(top50_tfs, aes(x = rank, y = avg_log2FC, fill = color_cat)) +
  geom_bar(stat = "identity", width = 1, color = "black", linewidth = 0.1) +
  scale_fill_manual(
    values = c(
      "Highly Up (Malignant)" = "#B2182B",
      "Up (Malignant)" = "#EF8A62",
      "Not Significant" = "grey70",
      "Down (Normal)" = "#67A9CF",
      "Highly Down (Normal)" = "#2166AC"
    ),
    name = "TF Activity Change"
  ) +
  geom_text_repel(
    data = top20_labels,
    aes(label = TF),
    size = 3.5,
    fontface = "bold",
    max.overlaps = 25,
    box.padding = 0.5,
    point.padding = 0.3,
    segment.size = 0.3,
    segment.color = "grey40"
  ) +
  geom_hline(yintercept = 0, linetype = "solid", color = "black", linewidth = 0.8) +
  geom_hline(yintercept = c(-0.5, 0.5), linetype = "dashed", color = "grey30", alpha = 0.5) +
  annotate("text", x = 5, y = 0.5, label = "logFC = 0.5", size = 3, color = "grey30", vjust = -0.5) +
  annotate("text", x = 5, y = -0.5, label = "logFC = -0.5", size = 3, color = "grey30", vjust = 1.5) +
  labs(
    title = "Differential Transcription Factor Activity",
    subtitle = "Malignant Sézary CD4+ T cells vs Normal CD4+ T cells (Top 50 TFs)",
    x = "Ranked Transcription Factors (by log2FC)",
    y = "log2 Fold Change (TF Activity)",
    caption = paste0("p-adjusted < 0.05, Top 20 labeled | Total TFs analyzed: ", nrow(tf_markers))
  ) +
  theme_classic() +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
    plot.subtitle = element_text(size = 12, hjust = 0.5),
    axis.title = element_text(size = 12, face = "bold"),
    axis.text.y = element_text(size = 10),
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    legend.position = "right",
    legend.title = element_text(size = 10, face = "bold"),
    plot.caption = element_text(size = 9, hjust = 0.5, color = "grey40")
  )

# Save the plot
ggsave("TF_waterfall_plot.pdf", width = 12, height = 8, dpi = 300)

9 Waterfall Plot for Each Cell Line

library(Seurat)
library(ggplot2)
library(dplyr)
library(cowplot)
library(ggrepel)

DefaultAssay(SS) <- "dorothea"

# Simpler function using group.by
create_waterfall_cellline_v2 <- function(seurat_obj, cell_line, top_n = 30) {
  
  # Find markers using group.by
  markers <- FindMarkers(
    seurat_obj,
    ident.1 = cell_line,
    ident.2 = c("CD4T_lab", "CD4T_10x"),  # Your normal controls from orig.ident
    group.by = "orig.ident",
    assay = "dorothea",
    logfc.threshold = 0,
    min.pct = 0,
    verbose = FALSE
  )
  
  if (nrow(markers) == 0) {
    warning(paste("No markers found for", cell_line))
    return(NULL)
  }
  
  # Prepare data
  markers$TF <- rownames(markers)
  markers_top <- markers %>%
    arrange(desc(abs(avg_log2FC))) %>%
    head(top_n) %>%
    arrange(avg_log2FC) %>%
    mutate(
      rank = row_number(),
      significant = ifelse(p_val_adj < 0.05, "Significant", "Not Significant")
    )
  
  # Top 15 for labels
  top_labels <- markers_top %>%
    arrange(desc(abs(avg_log2FC))) %>%
    head(15)
  
  # Plot
  p <- ggplot(markers_top, aes(x = rank, y = avg_log2FC, fill = significant)) +
    geom_bar(stat = "identity", width = 1, color = "black", linewidth = 0.1) +
    scale_fill_manual(
      values = c("Significant" = "#D62728", "Not Significant" = "grey70")
    ) +
    geom_text_repel(
      data = top_labels,
      aes(label = TF),
      size = 2.5,
      max.overlaps = 15,
      box.padding = 0.3
    ) +
    geom_hline(yintercept = 0, color = "black", linewidth = 0.5) +
    labs(
      title = paste0(cell_line, " vs Normal Controls"),
      x = "Ranked TFs",
      y = "log2FC (TF Activity)"
    ) +
    theme_minimal() +
    theme(
      axis.text.x = element_blank(),
      axis.ticks.x = element_blank(),
      legend.position = "none",
      plot.title = element_text(face = "bold", hjust = 0.5, size = 12)
    )
  
  return(p)
}

# Create waterfall plots
cell_lines <- c("L1", "L2", "L3", "L4", "L5", "L6", "L7")
plots <- lapply(cell_lines, function(x) create_waterfall_cellline_v2(SS, x, top_n = 30))

# Remove NULL plots
plots <- plots[!sapply(plots, is.null)]

# Combine
combined_plot <- plot_grid(plotlist = plots, ncol = 3, nrow = 3)

# Save
ggsave("TF_waterfall_celllines.pdf", 
       plot = combined_plot, 
       width = 16, height = 16, dpi = 300)

combined_plot

10 Export Results Table

# Export differential TF activity results
DefaultAssay(SS) <- "dorothea"
Idents(SS) <- SS$Condition

tf_markers <- FindMarkers(
  SS,
  ident.1 = "MalignantCD4T",
  ident.2 = "NormalCD4T",
  assay = "dorothea",
  logfc.threshold = 0,
  min.pct = 0
)

tf_markers$TF <- rownames(tf_markers)

# Export full results
write.csv(tf_markers, 
          "TF_differential_activity_MalignantVsNormal.csv", 
          row.names = FALSE)

# Export top 50
top50 <- tf_markers %>%
  arrange(desc(abs(avg_log2FC))) %>%
  head(50)

write.csv(top50, 
          "TF_differential_activity_top50.csv", 
          row.names = FALSE)

cat("Analysis complete!\n")

Analysis complete!

cat("Total TFs analyzed:", nrow(tf_markers), "\n")

Total TFs analyzed: 114 

cat("Significant TFs (p-adj < 0.05):", sum(tf_markers$p_val_adj < 0.05), "\n")

Significant TFs (p-adj < 0.05): 106