Figure 8D-L1 Intra-Clonal Regulatory Heterogeneity — TF State UMAP

1 load libraries

2 Load Object & Set RNA Assay

3 Define State-Defining TFs (Confirmed from Heatmap Panel C)

# orig.ident is set automatically during individual Seurat object creation
# It stores the sample name: L1, L2, L3, L4, L5, L6, L7, PBMC, PBMC-10x
seurat_L1 <- subset(seurat_obj, subset = orig.ident == "L1")

cat("Number of L1 cells:", ncol(seurat_L1), "\n")

Number of L1 cells: 5825 

# Sanity check — should show only L1
table(seurat_L1$orig.ident)

      L1       L2       L3       L4       L5       L6       L7 CD4T_lab CD4T_10x 
    5825        0        0        0        0        0        0        0        0 

4 Define 3 States (Exact match to manuscript Figure 8D)

# Inflammatory  → NF-κB/AP-1 driven (clusters 11, 12)
# Proliferative → MYC/E2F1 driven  (clusters 4, 7, 8)
# Th1_Cytotoxic → TBX21 driven     (cluster 5; stem-like via HMGA2 in text)
state_tfs <- list(
  Inflammatory  = c("STAT3", "RELA", "NFKB1"),
  Proliferative = c("MYC",   "E2F1"),
  Th1_Cytotoxic = c("TBX21")
)

5 Extract TF Activity Matrix from dorothea Assay

tf_mat <- GetAssayData(seurat_L1, assay = "dorothea", layer =  "data")

cat("TF matrix dimensions (TFs x cells):", dim(tf_mat), "\n")

TF matrix dimensions (TFs x cells): 271 5825 

# Check which state TFs are present in the dorothea assay
for (state in names(state_tfs)) {
  found   <- state_tfs[[state]][state_tfs[[state]] %in% rownames(tf_mat)]
  missing <- state_tfs[[state]][!state_tfs[[state]] %in% rownames(tf_mat)]
  cat(sprintf("%-15s Found: %-30s Missing: %s\n",
              state,
              paste(found,   collapse = ", "),
              paste(missing, collapse = ", ")))
}

Inflammatory    Found: STAT3, RELA, NFKB1             Missing: 
Proliferative   Found: MYC, E2F1                      Missing: 
Th1_Cytotoxic   Found: TBX21                          Missing: 

6 Compute Mean Activity Score per State per Cell

score_df <- as.data.frame(
  sapply(state_tfs, function(tfs) {
    tfs_found <- tfs[tfs %in% rownames(tf_mat)]
    if (length(tfs_found) == 0) {
      warning("No TFs found for this state — returning NA")
      return(rep(NA_real_, ncol(tf_mat)))
    }
    colMeans(tf_mat[tfs_found, , drop = FALSE])
  })
)

rownames(score_df) <- colnames(tf_mat)

# Preview raw scores
head(round(score_df, 3))

7 Scale Scores Across States

score_df_scaled <- as.data.frame(scale(score_df))
rownames(score_df_scaled) <- rownames(score_df)

# Preview scaled scores
head(round(score_df_scaled, 3))

8 Assign Dominant State per Cell

score_df_scaled$DominantState <- apply(
  score_df_scaled[, names(state_tfs), drop = FALSE], 1,
  function(x) {
    if (all(is.na(x))) return(NA_character_)
    names(which.max(x))
  }
)

seurat_L1$TF_State <- factor(
  score_df_scaled$DominantState,
  levels = names(state_tfs)
)

9 Also Add Raw Scores as Individual Metadata Columns

cat("\n=======================================================\n")

=======================================================

cat("     L1 TF State Distribution (Figure 8D)\n")

     L1 TF State Distribution (Figure 8D)

cat("=======================================================\n")

=======================================================

dist_table <- seurat_L1@meta.data %>%
  filter(!is.na(TF_State)) %>%
  count(TF_State) %>%
  mutate(Percentage = round(n / sum(n) * 100, 1)) %>%
  arrange(factor(TF_State, levels = names(state_tfs)))
print(dist_table)

       TF_State    n Percentage
1  Inflammatory 1917       32.9
2 Proliferative 2219       38.1
3 Th1_Cytotoxic 1689       29.0

cat("All 3 states must be > 0% to confirm intra-clonal heterogeneity\n")

All 3 states must be > 0% to confirm intra-clonal heterogeneity

cat("=======================================================\n")

=======================================================

# ── 12. Define Color Palette ──────────────────────────────────────────────────
state_colors <- c(
  Inflammatory  = "#E63946",    # red   — NF-κB/STAT3 survival program
  Proliferative = "#457B9D",    # blue  — MYC/E2F1 proliferative program
  Th1_Cytotoxic = "#2D6A4F"     # green — TBX21 Th1/cytotoxic program
)

10 Check UMAP is Present

cat("\nReductions available:", paste(Reductions(seurat_L1), collapse = ", "), "\n")

Reductions available: integrated_dr, ref.umap, pca, umap, harmony 

# Only run if UMAP is absent from the subset
if (!"umap" %in% Reductions(seurat_L1)) {
  cat("UMAP not found — recomputing from Harmony embeddings...\n")
  seurat_L1 <- RunUMAP(seurat_L1,
                        reduction = "harmony",
                        dims      = 1:20,
                        seed.use  = 42)
}

11 Panel A — UMAP Colored by TF State

p_umap <- DimPlot(
  seurat_L1,
  reduction  = "umap",
  group.by   = "TF_State",
  cols       = state_colors,
  pt.size    = 1.0,
  label      = FALSE,
  na.value   = "grey85"
) +
  ggtitle("Cell Line L1 — TF Regulatory State") +
  theme_classic(base_size = 11) +
  theme(
    plot.title      = element_text(face = "bold", size = 11, hjust = 0.5),
    legend.title    = element_text(size = 9,  face = "bold"),
    legend.text     = element_text(size = 8),
    legend.position = "right",
    axis.line       = element_line(color = "black", linewidth = 0.4),
    axis.text       = element_text(size = 7)
  ) +
  labs(color = "TF State")

p_umap

12 Panel B — Stacked Bar Chart (State Proportions)

state_counts <- seurat_L1@meta.data %>%
  filter(!is.na(TF_State)) %>%
  count(TF_State) %>%
  mutate(
    Proportion = n / sum(n) * 100,
    TF_State   = factor(TF_State, levels = names(state_colors))
  )

print(state_counts)

       TF_State    n Proportion
1  Inflammatory 1917   32.90987
2 Proliferative 2219   38.09442
3 Th1_Cytotoxic 1689   28.99571

p_bar <- ggplot(
  state_counts,
  aes(x = "L1", y = Proportion, fill = TF_State)
) +
  geom_bar(
    stat      = "identity",
    width     = 0.45,
    color     = "white",
    linewidth = 0.3
  ) +
  geom_text(
    aes(label = paste0(round(Proportion, 1), "%")),
    position = position_stack(vjust = 0.5),
    size     = 3.5,
    color    = "white",
    fontface = "bold"
  ) +
  scale_fill_manual(values = state_colors, drop = FALSE) +
  scale_y_continuous(
    limits = c(0, 100),
    breaks = seq(0, 100, 25),
    labels = function(x) paste0(x, "%")
  ) +
  labs(
    y     = "Proportion of cells (%)",
    x     = NULL,
    fill  = "TF State",
    title = "State distribution — L1"
  ) +
  theme_classic(base_size = 11) +
  theme(
    axis.text.x     = element_blank(),
    axis.ticks.x    = element_blank(),
    plot.title      = element_text(face = "bold", size = 11, hjust = 0.5),
    legend.position = "none",
    axis.line       = element_line(color = "black", linewidth = 0.4)
  )

p_bar

13 Combine Panels and Save - Figure 8D

final_fig <- p_umap + p_bar +
  plot_layout(widths = c(2.5, 1)) +
  plot_annotation(
    title    = "Intra-clonal regulatory heterogeneity — Cell Line L1",
    subtitle = "Monoclonal population simultaneously occupies Inflammatory, Proliferative, and Th1/Cytotoxic TF states",
    theme    = theme(
      plot.title    = element_text(size = 12, face = "bold",  hjust = 0.5),
      plot.subtitle = element_text(size = 9,  color = "grey40", hjust = 0.5)
    )
  )

final_fig

ggsave("Figure_8D_L1_TFstate_UMAP_Bar.pdf",
       plot = final_fig, width = 9, height = 4.5, dpi = 300, device = "pdf")

ggsave("Figure_8D_L1_TFstate_UMAP_Bar.png",
       plot = final_fig, width = 9, height = 4.5, dpi = 300)

cat("Figure 8D saved successfully.\n")

Figure 8D saved successfully.

14 Validation — Violin Plots (Supplementary Figure SX)

score_features <- paste0("Score_", names(state_tfs))

p_val <- VlnPlot(
  seurat_L1,
  features = score_features,
  group.by = "TF_State",
  pt.size  = 0,
  cols     = state_colors,
  ncol     = 3
) &
  theme_classic(base_size = 9) &
  theme(
    axis.text.x     = element_text(angle = 45, hjust = 1, size = 7),
    legend.position = "none",
    plot.title      = element_text(size = 9)
  )

p_val

ggsave("Supp_SX_TF_validation_violins.pdf",
       plot = p_val, width = 12, height = 4.5, dpi = 300)

ggsave("Supp_SX_TF_validation_violins.png",
       plot = p_val, width = 12, height = 4.5, dpi = 300, bg = "white")

cat("Validation violins saved as PDF and PNG.\n")

Validation violins saved as PDF and PNG.

15 # Session Info

sessionInfo()

R version 4.5.2 (2025-10-31)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 24.04.3 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.12.0 
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0  LAPACK version 3.12.0

locale:
 [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C               LC_TIME=fr_FR.UTF-8       
 [4] LC_COLLATE=en_GB.UTF-8     LC_MONETARY=fr_FR.UTF-8    LC_MESSAGES=en_GB.UTF-8   
 [7] LC_PAPER=fr_FR.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
[10] LC_TELEPHONE=C             LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C       

time zone: Europe/Paris
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] RColorBrewer_1.1-3 patchwork_1.3.2    dplyr_1.2.0        ggplot2_4.0.2     
[5] Seurat_5.4.0       SeuratObject_5.3.0 sp_2.2-1          

loaded via a namespace (and not attached):
  [1] deldir_2.0-4           pbapply_1.7-4          gridExtra_2.3         
  [4] rlang_1.1.7            magrittr_2.0.4         RcppAnnoy_0.0.23      
  [7] otel_0.2.0             spatstat.geom_3.7-0    matrixStats_1.5.0     
 [10] ggridges_0.5.7         compiler_4.5.2         systemfonts_1.3.1     
 [13] png_0.1-8              vctrs_0.7.1            reshape2_1.4.5        
 [16] stringr_1.6.0          pkgconfig_2.0.3        fastmap_1.2.0         
 [19] labeling_0.4.3         promises_1.5.0         rmarkdown_2.30        
 [22] ggbeeswarm_0.7.3       ragg_1.5.0             purrr_1.2.1           
 [25] xfun_0.56              cachem_1.1.0           jsonlite_2.0.0        
 [28] goftest_1.2-3          later_1.4.5            spatstat.utils_3.2-1  
 [31] irlba_2.3.7            parallel_4.5.2         cluster_2.1.8.2       
 [34] R6_2.6.1               ica_1.0-3              spatstat.data_3.1-9   
 [37] bslib_0.10.0           stringi_1.8.7          reticulate_1.44.1     
 [40] spatstat.univar_3.1-6  parallelly_1.46.1      lmtest_0.9-40         
 [43] jquerylib_0.1.4        scattermore_1.2        Rcpp_1.1.1            
 [46] knitr_1.51             tensor_1.5.1           future.apply_1.20.1   
 [49] zoo_1.8-15             sctransform_0.4.3      httpuv_1.6.16         
 [52] Matrix_1.7-4           splines_4.5.2          igraph_2.2.2          
 [55] tidyselect_1.2.1       abind_1.4-8            rstudioapi_0.18.0     
 [58] dichromat_2.0-0.1      yaml_2.3.12            spatstat.random_3.4-4 
 [61] codetools_0.2-20       miniUI_0.1.2           spatstat.explore_3.7-0
 [64] listenv_0.10.0         lattice_0.22-9         tibble_3.3.1          
 [67] plyr_1.8.9             withr_3.0.2            shiny_1.12.1          
 [70] S7_0.2.1               ROCR_1.0-12            ggrastr_1.0.2         
 [73] evaluate_1.0.5         Rtsne_0.17             future_1.69.0         
 [76] fastDummies_1.7.5      survival_3.8-3         polyclip_1.10-7       
 [79] fitdistrplus_1.2-6     pillar_1.11.1          rsconnect_1.7.0       
 [82] KernSmooth_2.23-26     plotly_4.12.0          generics_0.1.4        
 [85] RcppHNSW_0.6.0         scales_1.4.0           globals_0.19.0        
 [88] xtable_1.8-4           glue_1.8.0             lazyeval_0.2.2        
 [91] tools_4.5.2            data.table_1.18.2.1    RSpectra_0.16-2       
 [94] RANN_2.6.2             dotCall64_1.2          cowplot_1.2.0         
 [97] grid_4.5.2             tidyr_1.3.2            nlme_3.1-168          
[100] beeswarm_0.4.0         vipor_0.4.7            cli_3.6.5             
[103] spatstat.sparse_3.1-0  textshaping_1.0.4      spam_2.11-3           
[106] viridisLite_0.4.3      uwot_0.2.4             gtable_0.3.6          
[109] sass_0.4.10            digest_0.6.39          progressr_0.18.0      
[112] ggrepel_0.9.6          htmlwidgets_1.6.4      farver_2.1.2          
[115] htmltools_0.5.9        lifecycle_1.0.5        httr_1.4.7            
[118] mime_0.13              MASS_7.3-65