Benchmark cell Prediction Methods-5-012-2025 on Clusters

load libraries————————————

1. Load and Subset Normal CD4 T Cells

Idents(All_samples_Merged) <- "seurat_clusters"
table(All_samples_Merged@active.ident)

   0    1    2    3    4    5    6    7    8    9   10   11   12   13 
6789 5275 4663 4661 4086 3634 3536 3409 3338 3273 3212 1675 1063  691 

make sure your labels are present

colnames(scATOMIC_obj@meta.data)

 [1] "orig.ident"                         "nCount_RNA"                        
 [3] "nFeature_RNA"                       "nCount_ADT"                        
 [5] "nFeature_ADT"                       "nUMI"                              
 [7] "ngene"                              "cell_line"                         
 [9] "Patient_origin"                     "Patient_Immunophenotype"           
[11] "condition_of_amplification_in_vivo" "culture_medium"                    
[13] "Stromal_cells"                      "Cell_line_Immunophenotype"         
[15] "TP53_mutation"                      "age_at_diagnosis"                  
[17] "stage_diagnosis"                    "stage_analysis"                    
[19] "Treatments_analysis"                "TCRVB2"                            
[21] "CD3_M"                              "CD30_M"                            
[23] "CCR4_M"                             "CD162_BL"                          
[25] "CD26_BL"                            "CD7_M"                             
[27] "CLA_BD"                             "CD4_BD"                            
[29] "CCR7_M"                             "CD45RO_BD"                         
[31] "CD45RA_BD"                          "CellName"                          
[33] "percent.mt"                         "CD26_BD"                           
[35] "CD45RA_M"                           "predicted.celltype.l1.score"       
[37] "predicted.celltype.l1"              "predicted.celltype.l2.score"       
[39] "predicted.celltype.l2"              "predicted.celltype.l3.score"       
[41] "predicted.celltype.l3"              "mapping.score"                     
[43] "percent.rb"                         "nCount_SCT"                        
[45] "nFeature_SCT"                       "S.Score"                           
[47] "G2M.Score"                          "Phase"                             
[49] "old.ident"                          "CC.Difference"                     
[51] "SCT_snn_res.0.4"                    "SCT_snn_res.0.5"                   
[53] "SCT_snn_res.0.6"                    "SCT_snn_res.0.7"                   
[55] "SCT_snn_res.0.8"                    "seurat_clusters"                   
[57] "harmony_res_0.1"                    "harmony_res_0.2"                   
[59] "harmony_res_0.3"                    "harmony_res_0.4"                   
[61] "harmony_res_0.5"                    "harmony_res_0.6"                   
[63] "harmony_res_0.7"                    "harmony_res_0.8"                   
[65] "harmony_res_0.9"                    "harmony_res_1"                     
[67] "harmony_res_1.2"                    "predicted.celltype.l1_backup"      
[69] "predicted.celltype.l2_backup"       "predicted.celltype.l3_backup"      
[71] "Prediction"                         "Cluster"                           
[73] "Uncertainty_score"                  "Condition"                         
[75] "scATOMIC_pred"                     

Transfer scATOMIC labels into your main object

All_samples_Merged$scATOMIC_annotation <- scATOMIC_obj@meta.data[colnames(All_samples_Merged), "scATOMIC_pred"]

2. Create the annotation summary table

library(dplyr)
library(tidyr)
library(purrr)

methods <- c(
  "predicted.id",             #scPred
  "predicted.celltype.l2",    # Azimuth (l2 prediction)
  "singler.immune",           # Example of other annotation methods like SingleR
  "Prediction",               #STCAT
  "scATOMIC_annotation"      #scATOMIC
  )

annotation_summary <- map_dfr(methods, function(m) {
  df <- All_samples_Merged@meta.data
  df %>%
    group_by(seurat_clusters, label = .data[[m]]) %>%
    summarise(n = n(), .groups = "drop") %>%
    group_by(seurat_clusters) %>%
    slice_max(n, n = 1, with_ties = FALSE) %>%
    mutate(method = m)
})



# Set the order of methods
annotation_summary$method <- factor(annotation_summary$method, levels = methods)

head(annotation_summary)

# Should show columns: seurat_clusters | label | n | method

3. Convert to matrix format for plotting

library(tidyr)
library(ggplot2)

annotation_matrix <- annotation_summary %>%
  select(method, seurat_clusters, label) %>%
  pivot_wider(names_from = seurat_clusters, values_from = label)

# Rename methods for the plot
annotation_summary <- annotation_summary %>%
  mutate(method = recode(method,
                         "predicted.id" = "scPred",
                         "predicted.celltype.l2" = "Azimuth.l2",
                         "singler.immune" = "SingleR(Immune)",
                         "Prediction" = "STCAT",
                         "scATOMIC_annotation" = "scATOMIC"))

#4. Creating a Table for All Annotations

# Summarize the number of cells assigned to each label by method and cell line
annotation_table <- annotation_summary %>%
  group_by(seurat_clusters, method, label) %>%
  summarise(n_cells = n(), .groups = 'drop') %>%
  spread(key = label, value = n_cells, fill = 0)  # Pivot the table to have labels as columns

# Display the summary table
annotation_table

#  Save the table to a CSV file
write.csv(annotation_table, "annotation_summary_table_clusters-5-12-2025.csv", row.names = FALSE)

#5. Plot categorical heatmap (tile plot)

library(ggplot2)

ggplot(annotation_summary, aes(x = seurat_clusters, y = method, fill = label)) +
  geom_tile(color = "white", linewidth = 0.5) +
  scale_fill_discrete() +   # You can assign your own palette later
  labs(
    x = "Sézary cell line",
    y = "Annotation method",
    fill = "Assigned state",
    title = "Cross-method comparison of T-cell state annotations"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    axis.title = element_text(face = "bold")
  )

#6. Plot categorical heatmap (tile plot)

library(RColorBrewer)

num_colors <- length(unique(annotation_summary$label))
my_colors <- colorRampPalette(brewer.pal(8, "Set2"))(num_colors)

ggplot(annotation_summary, aes(x = seurat_clusters, y = method, fill = label)) +
  geom_tile(color = "white", linewidth = 0.5) +
  scale_fill_manual(values = my_colors) +   # 22 or more colors here
  labs(
    x = "Sézary cell line clusters",
    y = "Annotation method",
    fill = "Assigned state",
    title = "Cross-method comparison of T-cell state annotations"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    axis.title = element_text(face = "bold")
  )

library(RColorBrewer)

num_colors <- length(unique(annotation_summary$label))
my_colors <- colorRampPalette(brewer.pal(8, "Set2"))(num_colors)

p <- ggplot(annotation_summary, aes(x = seurat_clusters, y = method, fill = label)) +
  geom_tile(color = "white", linewidth = 0.5) +
  scale_fill_manual(values = my_colors) +
  labs(
    x = "Sézary cell line clusters",
    y = "Annotation method",
    fill = "Assigned state",
    title = "Cross-method comparison of T-cell state annotations"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    axis.title = element_text(face = "bold")
  )
ggsave("annotation_heatmap.png", plot = p, width = 14, height = 6, units = "in", dpi = 300)
ggsave("annotation_heatmap.pdf", plot = p, width = 14, height = 6, units = "in")