PAGA using TF Activity

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

library(Seurat)
library(reticulate)
library(ggplot2)

# 1. Load Object
seurat_obj <- readRDS("temp_seurat_obj.rds")

3 SUPPLEMENTARY FIGURE: PAGA Trajectory Analysis (via Reticulate/Scanpy)-TF based

# Display in notebook
if(file.exists(output_file)){
  knitr::include_graphics(output_file)
}

4 SUPPLEMENTARY FIGURE: PAGA on Gene Expression (SCT Assay)

5 SUPPLEMENTARY FIGURE: Malignant-Only PAGA-SCT Assay

6 SUPPLEMENTARY FIGURE: Malignant-Only PAGA-TF Assay

# Display
output_file <- "Output_Figures/PAGA_Malignant_TF.png"

if (file.exists(output_file)) {
  knitr::include_graphics(output_file)
}


cat("=== TF-BASED PAGA COMPLETE ===\n")

=== TF-BASED PAGA COMPLETE ===

7 Verification of TF Availability

library(Seurat)
library(dplyr)

cat("=== CHECKING TF AVAILABILITY IN DOROTHEA ASSAY ===\n\n")

=== CHECKING TF AVAILABILITY IN DOROTHEA ASSAY ===

# 1. Identify the TF assay
available_assays <- Assays(seurat_obj)
tf_assay_name <- grep("dorothea|viper|TF", available_assays, value=TRUE)[1]

if(is.na(tf_assay_name)) {
  stop("❌ No TF assay found! Please run DoRothEA first.")
}
cat("✓ Found TF assay:", tf_assay_name, "\n")

✓ Found TF assay: dorothea 

# 2. Get all available TFs
all_tfs <- rownames(seurat_obj[[tf_assay_name]])
cat("✓ Total TFs inferred:", length(all_tfs), "\n\n")

✓ Total TFs inferred: 271 

# 3. Define our Desired Panel (Literature-Based)
desired_tfs <- list(
  "Stemness" = c("TCF7", "LEF1", "BACH2", "FOXO1", "MYC"),
  "Malignancy" = c("TOX", "STAT5A", "STAT3", "STAT5B", "GATA3", "IKZF2"),
  "Exhaustion" = c("EOMES", "TBX21", "PRDM1", "NFATC1", "BATF"),
  "Treg_Like" = c("FOXP3", "IKZF2", "SATB1") 
)

# 4. Check Presence
available_panel <- list()
missing_tfs <- c()

cat("--- PANEL VERIFICATION ---\n")

--- PANEL VERIFICATION ---

for(category in names(desired_tfs)) {
  tfs <- desired_tfs[[category]]
  present <- intersect(tfs, all_tfs)
  missing <- setdiff(tfs, all_tfs)
  
  available_panel[[category]] <- present
  missing_tfs <- c(missing_tfs, missing)
  
  cat(sprintf("%s: %d/%d found\n", category, length(present), length(tfs)))
  if(length(present) > 0) cat("  ✅ Present:", paste(present, collapse=", "), "\n")
  if(length(missing) > 0) cat("  ❌ Missing:", paste(missing, collapse=", "), "\n")
  cat("\n")
}

Stemness: 5/5 found
  ✅ Present: TCF7, LEF1, BACH2, FOXO1, MYC 

Malignancy: 4/6 found
  ✅ Present: STAT5A, STAT3, STAT5B, GATA3 
  ❌ Missing: TOX, IKZF2 

Exhaustion: 5/5 found
  ✅ Present: EOMES, TBX21, PRDM1, NFATC1, BATF 

Treg_Like: 0/3 found
  ❌ Missing: FOXP3, IKZF2, SATB1 

# 5. Summary
final_tf_list <- unlist(available_panel)
cat("=== FINAL TF LIST FOR TRAJECTORY ===\n")

=== FINAL TF LIST FOR TRAJECTORY ===

cat("Total TFs to plot:", length(final_tf_list), "\n")

Total TFs to plot: 14 

print(final_tf_list)

  Stemness1   Stemness2   Stemness3   Stemness4   Stemness5 Malignancy1 Malignancy2 Malignancy3 Malignancy4 Exhaustion1 Exhaustion2 Exhaustion3 
     "TCF7"      "LEF1"     "BACH2"     "FOXO1"       "MYC"    "STAT5A"     "STAT3"    "STAT5B"     "GATA3"     "EOMES"     "TBX21"     "PRDM1" 
Exhaustion4 Exhaustion5 
   "NFATC1"      "BATF" 

if("TOX" %in% final_tf_list) cat("\n✓ CRITICAL: TOX is available!\n")
if("TCF7" %in% final_tf_list) cat("✓ CRITICAL: TCF7 is available!\n")

✓ CRITICAL: TCF7 is available!

if("STAT5A" %in% final_tf_list) cat("✓ CRITICAL: STAT5A is available!\n")

✓ CRITICAL: STAT5A is available!

# 6. Save valid list for next step
saveRDS(final_tf_list, "valid_tf_panel.rds")

8 PAGA Gene Dynamics for Malignant Sézary Cells (TF)

library(reticulate)
library(Seurat)
library(dplyr)

# Setup Python
sc <- import("scanpy")
ad <- import("anndata")
pl <- import("matplotlib.pyplot")

cat("=== TF REGULATORY DYNAMICS (5 TRAJECTORIES) ===\n\n")

=== TF REGULATORY DYNAMICS (5 TRAJECTORIES) ===

# 1. IDENTIFY TF ASSAY & CHECK AVAILABILITY
available_assays <- Assays(seurat_obj)
tf_assay_name <- grep("dorothea|viper|TF", available_assays, value=TRUE)[1]

if(is.na(tf_assay_name)) {
  stop("❌ No TF assay found! Please run DoRothEA first.")
}
cat("✓ Using TF assay:", tf_assay_name, "\n")

✓ Using TF assay: dorothea 

# Get all available TFs
all_tfs_in_object <- rownames(seurat_obj[[tf_assay_name]])
cat("✓ Total TFs inferred:", length(all_tfs_in_object), "\n\n")

✓ Total TFs inferred: 271 

# 2. DEFINE & VALIDATE TF PANEL
desired_tfs <- list(
  "Stemness" = c("TCF7", "LEF1", "BACH2", "FOXO1", "MYC"),
  "Malignancy" = c("TOX", "STAT5A", "STAT3", "STAT5B", "GATA3", "IKZF2"),
  "Exhaustion" = c("EOMES", "TBX21", "PRDM1", "NFATC1", "BATF"),
  "Treg_Like" = c("FOXP3", "SATB1")
)

flat_valid_tfs <- c()
cat("--- TF PANEL VALIDATION ---\n")

--- TF PANEL VALIDATION ---

for(category in names(desired_tfs)) {
  tfs <- desired_tfs[[category]]
  present <- intersect(tfs, all_tfs_in_object)
  missing <- setdiff(tfs, all_tfs_in_object)
  
  if(length(present) > 0) {
    flat_valid_tfs <- c(flat_valid_tfs, present)
    cat(sprintf("✅ %s: %s\n", category, paste(present, collapse=", ")))
  }
  if(length(missing) > 0) {
    cat(sprintf("❌ %s (Missing): %s\n", category, paste(missing, collapse=", ")))
  }
}

✅ Stemness: TCF7, LEF1, BACH2, FOXO1, MYC
✅ Malignancy: STAT5A, STAT3, STAT5B, GATA3
❌ Malignancy (Missing): TOX, IKZF2
✅ Exhaustion: EOMES, TBX21, PRDM1, NFATC1, BATF
❌ Treg_Like (Missing): FOXP3, SATB1

if(length(flat_valid_tfs) == 0) stop("No TFs found!")
cat("\n✓ Using", length(flat_valid_tfs), "TFs\n\n")

✓ Using 14 TFs

# 3. PREPARE DATA (Malignant Only) - CLEAN RECREATION
malignant_clusters <- setdiff(unique(seurat_obj$seurat_clusters), c(3, 10))
seurat_malignant <- subset(seurat_obj, subset = seurat_clusters %in% malignant_clusters)

DefaultAssay(seurat_malignant) <- tf_assay_name
tf_matrix <- GetAssayData(seurat_malignant, assay = tf_assay_name, layer = "data")
tf_matrix <- tf_matrix[flat_valid_tfs, ]

# Convert to dense matrix and transpose
expr_matrix <- t(as.matrix(tf_matrix))

# Create clean obs dataframe
obs_df <- data.frame(
  seurat_clusters = as.character(seurat_malignant$seurat_clusters),
  row.names = colnames(seurat_malignant),
  stringsAsFactors = FALSE
)

# Create AnnData using Python directly to avoid reticulate issues
cat("Creating AnnData object...\n")

Creating AnnData object...

# Export data to Python environment
py$expr_matrix <- expr_matrix
py$obs_df <- obs_df
py$var_names <- flat_valid_tfs

py_run_string("
import scanpy as sc
import anndata as ad
import pandas as pd
import numpy as np

# Create AnnData cleanly in Python
adata_tf = ad.AnnData(
    X=r.expr_matrix,
    obs=r.obs_df
)

# Set variable names
adata_tf.var_names = r.var_names

print(f'✓ AnnData created: {adata_tf.n_obs} cells x {adata_tf.n_vars} TFs')
")

Error in eval(parse(text = as_r_value(code)), envir = envir) : 
  RuntimeError: object 'var_names' not found
Run ]8;;rstudio:run:reticulate::py_last_error()`reticulate::py_last_error()`]8;; for details.

9 PAGA Gene Dynamics for Malignant Sézary Cells (SCT)