Sézary Syndrome silhouette Score Test

1. load libraries

#Differential Expression Analysis

2. load seurat object

#Load Seurat Object L7
load("../0-IMP-OBJECTS/Harmony_integrated_All_samples_Merged_with_PBMC10x_with_harmony_clustering.Robj")

Le chargement a nécessité le package : SeuratObject
Le chargement a nécessité le package : sp

Attachement du package : ‘SeuratObject’

Les objets suivants sont masqués depuis ‘package:base’:

    intersect, t

All_samples_Merged

Le chargement a nécessité le package : Seurat
Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

An object of class Seurat 
64169 features across 59355 samples within 6 assays 
Active assay: SCT (27417 features, 3000 variable features)
 3 layers present: counts, data, scale.data
 5 other assays present: RNA, ADT, prediction.score.celltype.l1, prediction.score.celltype.l2, prediction.score.celltype.l3
 6 dimensional reductions calculated: integrated_dr, ref.umap, pca, umap, harmony, umap.harmony

DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "cell_line",label = T, label.box = T)

DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "Harmony_snn_res.0.9",label = T, label.box = T)

3. silhouette Score Test

library(Seurat)
library(cluster)
library(ggplot2)

# Assuming `All_samples_Merged` is already created and contains both PCA and Harmony embeddings
# Define the range of PCs used in SCTransform
pc_range <- 1:13

# Extract PCA and Harmony embeddings
pca_embedding <- Embeddings(All_samples_Merged, reduction = "pca")[, pc_range]
harmony_embedding <- Embeddings(All_samples_Merged, reduction = "harmony")[, pc_range]

# Retrieve cluster labels and batch labels
cluster_labels <- All_samples_Merged$seurat_clusters  # Adjust based on actual cluster column
batch_labels <- All_samples_Merged$orig.ident        # Adjust based on actual batch/condition column

# Ensure labels are numeric (for silhouette function)
cluster_labels <- as.numeric(as.factor(cluster_labels))
batch_labels <- as.numeric(as.factor(batch_labels))

# Compute silhouette scores for PCA embeddings
silhouette_pca_clusters <- silhouette(cluster_labels, dist(pca_embedding))
silhouette_pca_batches <- silhouette(batch_labels, dist(pca_embedding))

# Compute silhouette scores for Harmony embeddings
silhouette_harmony_clusters <- silhouette(cluster_labels, dist(harmony_embedding))
silhouette_harmony_batches <- silhouette(batch_labels, dist(harmony_embedding))

# Convert silhouette objects to data frames for visualization
silhouette_df <- data.frame(
  method = rep(c("PCA", "Harmony"), each = length(cluster_labels) * 2),
  group = c(rep("Cluster", length(cluster_labels)), rep("Batch", length(cluster_labels))),
  silhouette_width = c(
    silhouette_pca_clusters[, "sil_width"], silhouette_pca_batches[, "sil_width"],
    silhouette_harmony_clusters[, "sil_width"], silhouette_harmony_batches[, "sil_width"]
  )
)

# Plot silhouette scores for clusters and batches across PCA and Harmony
ggplot(silhouette_df, aes(x = method, y = silhouette_width, fill = group)) +
  geom_boxplot() +
  theme_minimal() +
  labs(
    title = "Comparison of Silhouette Scores for PCA and Harmony",
    x = "Embedding Method",
    y = "Silhouette Width"
  ) +
  facet_wrap(~ group, scales = "free") +
  theme(legend.position = "none")

# Print average silhouette scores
cat("Average silhouette scores:\n")

Average silhouette scores:

cat("PCA (Cluster): ", mean(silhouette_pca_clusters[, "sil_width"]), "\n")

PCA (Cluster):  0.05326484 

cat("PCA (Batch): ", mean(silhouette_pca_batches[, "sil_width"]), "\n")

PCA (Batch):  0.3232747 

cat("Harmony (Cluster): ", mean(silhouette_harmony_clusters[, "sil_width"]), "\n")

Harmony (Cluster):  0.07464966 

cat("Harmony (Batch): ", mean(silhouette_harmony_batches[, "sil_width"]), "\n")

Harmony (Batch):  0.05910068