Sézary Syndrome silhouette Score Test-Part2

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")


All_samples_Merged

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 predicted cell type labels and batch labels
celltype_labels <- All_samples_Merged$predicted.celltype.l2  # Using predicted.celltype.l2
batch_labels <- All_samples_Merged$orig.ident        # Adjust based on actual batch/condition column

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

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

# Compute silhouette scores for Harmony embeddings
silhouette_harmony_celltypes <- silhouette(celltype_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(celltype_labels) * 2),
  group = c(rep("Cell Type", length(celltype_labels)), rep("Batch", length(celltype_labels))),
  silhouette_width = c(
    silhouette_pca_celltypes[, "sil_width"], silhouette_pca_batches[, "sil_width"],
    silhouette_harmony_celltypes[, "sil_width"], silhouette_harmony_batches[, "sil_width"]
  )
)

# Plot silhouette scores for cell types 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 (Cell Type): ", mean(silhouette_pca_celltypes[, "sil_width"]), "\n")

PCA (Cell Type):  -0.1546204 

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

PCA (Batch):  0.3232747 

cat("Harmony (Cell Type): ", mean(silhouette_harmony_celltypes[, "sil_width"]), "\n")

Harmony (Cell Type):  -0.1344296 

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

Harmony (Batch):  0.05910068