TF Activity Inference Analysis (DecoupleR + DoRothEA)

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)

# Create output directories to keep project clean
dir.create("Output_Figures", showWarnings = FALSE)
dir.create("Output_Tables", showWarnings = FALSE)
dir.create("Output_Objects", showWarnings = FALSE)

2 Load Seurat Object

# Load your Seurat Object
seurat_obj <- readRDS("../0-Seurat_RDS_OBJECT_FINAL/All_samples_Merged_with_Renamed_Clusters_Cell_state-03-12-2025.rds.rds")

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

3 Infer TF Activity (Run Weighted Mean)

# Select Assay: SCT or RNA
assay_to_use <- "SCT" 

# 1. Load Local DoRothEA Network
library(dorothea)
data(dorothea_hs, package = "dorothea")

# Filter for High/Medium Confidence (A, B, C)
network <- dorothea_hs %>%
  dplyr::filter(confidence %in% c("A", "B", "C"))

print(paste("Running TF inference using", assay_to_use, "assay..."))

# --- OPTIMIZATION: Filter for Shared Genes to Save Memory ---
# Identify genes present in both the Seurat object and the Network
network_genes <- unique(c(network$tf, network$target))
genes_to_keep <- intersect(rownames(seurat_obj@assays[[assay_to_use]]@data), network_genes)

print(paste("Original features:", nrow(seurat_obj)))
print(paste("Features optimized for inference:", length(genes_to_keep)))

# Create smaller matrix (only relevant genes)
expression_matrix <- as.matrix(seurat_obj@assays[[assay_to_use]]@data[genes_to_keep, ])

# Load Parallel Library
library(BiocParallel)

# Set up Parallel Backend
# If on Windows, use SnowParam(). If Linux/Mac, MulticoreParam() is better.
# Let's use SnowParam as it is safer for stability.
n_cores <- 6  # Start conservative. Don't use all cores yet.
bpparam <- SnowParam(workers = n_cores, progressbar = TRUE)

print(paste("Running DecoupleR on", n_cores, "cores..."))

gc()  # Clean memory before running inference)

# 2. Run Weighted Mean (VIPER-like method)
tic("Run DecoupleR")
activities <- decoupleR::run_wmean(mat = expression_matrix,
                                   network = network,
                                   .source = "tf",       # Must match DoRothEA column name
                                   .target = "target",
                                   .mor = "mor",
                                   times = 100,
                                   minsize = 5)
gc()  # Clean memory before running inference)
toc()
gc()  # Clean memory before running inference)
# 3. Reshape and Add as a New Assay ('dorothea')
# Transform long-format to matrix (Features x Cells)
tf_activity_matrix <- activities %>%
  filter(statistic == "norm_wmean") %>%
  pivot_wider(id_cols = "source", names_from = "condition", values_from = "score") %>%
  column_to_rownames("source") %>%
  as.matrix()

# Ensure strict column alignment with Seurat object cells
# (Note: Using expression_matrix colnames ensures perfect match)
common_cells <- colnames(expression_matrix)
tf_activity_matrix <- tf_activity_matrix[, common_cells]

# Add to Seurat Object
seurat_obj[["dorothea"]] <- CreateAssayObject(data = tf_activity_matrix)

# Scale the data (Required for Heatmaps)
seurat_obj <- ScaleData(seurat_obj, assay = "dorothea")

# 4. SAVE OBJECT WITH TF ACTIVITY
saveRDS(seurat_obj, file = "Output_Objects/Seurat_Object_With_TF_Activity.rds")
print("TF Assay added and Object Saved.")