Sézary Syndrome Cell Line derived from patient Analysis_NewUMAP_Wilcox_RNA_Assay

1. load libraries

#Differential Expression Analysis

2. load seurat object

All_samples_Merged <- readRDS("/home/nabbasi/isilon/PHD_3rd_YEAR_Analysis/0-Seurat_RDS_OBJECT_FINAL/All_samples_Merged_Harmony_integrated_Cell_line_renamed_03-07-2025.rds")

Set Up Identifiers for Clustering

# Assign cluster identities to the Seurat object
Idents(All_samples_Merged) <- "seurat_clusters"

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

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

#Differential Expression Analysis # 3. Pairwise Comparisons

# Perform comparisons dynamically for each cell line pair
p1_comparison <- perform_comparison_and_volcano(All_samples_Merged, "L1", "L2", expression_data_RNA)

Avis : One or more p-values is 0. Converting to 10^-1 * current lowest non-zero p-value...

p2_comparison <- perform_comparison_and_volcano(All_samples_Merged, "L3", "L4", expression_data_RNA)

Avis : One or more p-values is 0. Converting to 10^-1 * current lowest non-zero p-value...

p3_comparison_L5L6 <- perform_comparison_and_volcano(All_samples_Merged, "L5", "L6", expression_data_RNA)

Avis : One or more p-values is 0. Converting to 10^-1 * current lowest non-zero p-value...

p3_comparison_L5L7 <- perform_comparison_and_volcano(All_samples_Merged, "L5", "L7", expression_data_RNA)

Avis : One or more p-values is 0. Converting to 10^-1 * current lowest non-zero p-value...

p3_comparison_L6L7 <- perform_comparison_and_volcano(All_samples_Merged, "L6", "L7", expression_data_RNA)

Avis : One or more p-values is 0. Converting to 10^-1 * current lowest non-zero p-value...

# Print the first few rows of each result
head(p1_comparison)

head(p2_comparison)

head(p3_comparison_L5L6)

head(p3_comparison_L5L7)

head(p3_comparison_L6L7)

4. Pairwise Comparisons Filtering

library(dplyr)

# Function to summarize DE markers
summarize_markers <- function(markers, comparison_name) {
  num_pval0 <- sum(markers$p_val_adj == 0)
  num_pval1 <- sum(markers$p_val_adj == 1)
  num_significant <- sum(markers$p_val_adj < 0.05)
  num_upregulated <- sum(markers$avg_log2FC > 1)
  num_downregulated <- sum(markers$avg_log2FC < -1)
  
  cat("\nSummary for", comparison_name, ":\n")
  cat("Number of genes with p_val_adj = 0:", num_pval0, "\n")
  cat("Number of genes with p_val_adj = 1:", num_pval1, "\n")
  cat("Number of significant genes (p_val_adj < 0.05):", num_significant, "\n")
  cat("Number of upregulated genes (avg_log2FC > 1):", num_upregulated, "\n")
  cat("Number of downregulated genes (avg_log2FC < -1):", num_downregulated, "\n")
}

# Function to filter DE markers and save results
filter_and_save <- function(markers, ident1, ident2) {
  comparison_name <- paste(ident1, "vs", ident2)
  
  # Summary before filtering
  summarize_markers(markers, comparison_name)
  
  # Apply the expression filter
  markers_filtered <- markers %>%
    filter(!(mean_expr_ident1 < 0.2 & mean_expr_ident2 < 0.2))

  # Save filtered results
  output_filename <- paste0("comparison_", ident1, "_vs_", ident2, "_with_mean_expression_filtered.csv")
  write.csv(markers_filtered, file = output_filename, row.names = FALSE)
  
  cat("Filtered results saved to:", output_filename, "\n")
  
  return(markers_filtered)
}

# Apply the function to all comparisons
p1_filtered <- filter_and_save(p1_comparison, "L1", "L2")

Summary for L1 vs L2 :
Number of genes with p_val_adj = 0: 1301 
Number of genes with p_val_adj = 1: 2209 
Number of significant genes (p_val_adj < 0.05): 10169 
Number of upregulated genes (avg_log2FC > 1): 1297 
Number of downregulated genes (avg_log2FC < -1): 3216 
Filtered results saved to: comparison_L1_vs_L2_with_mean_expression_filtered.csv 

p2_filtered <- filter_and_save(p2_comparison, "L3", "L4")

Summary for L3 vs L4 :
Number of genes with p_val_adj = 0: 389 
Number of genes with p_val_adj = 1: 2998 
Number of significant genes (p_val_adj < 0.05): 8391 
Number of upregulated genes (avg_log2FC > 1): 943 
Number of downregulated genes (avg_log2FC < -1): 1537 
Filtered results saved to: comparison_L3_vs_L4_with_mean_expression_filtered.csv 

p3_filtered_L5L6 <- filter_and_save(p3_comparison_L5L6, "L5", "L6")

Summary for L5 vs L6 :
Number of genes with p_val_adj = 0: 375 
Number of genes with p_val_adj = 1: 3837 
Number of significant genes (p_val_adj < 0.05): 7656 
Number of upregulated genes (avg_log2FC > 1): 1350 
Number of downregulated genes (avg_log2FC < -1): 1700 
Filtered results saved to: comparison_L5_vs_L6_with_mean_expression_filtered.csv 

p3_filtered_L5L7 <- filter_and_save(p3_comparison_L5L7, "L5", "L7")

Summary for L5 vs L7 :
Number of genes with p_val_adj = 0: 189 
Number of genes with p_val_adj = 1: 4450 
Number of significant genes (p_val_adj < 0.05): 5982 
Number of upregulated genes (avg_log2FC > 1): 1063 
Number of downregulated genes (avg_log2FC < -1): 1129 
Filtered results saved to: comparison_L5_vs_L7_with_mean_expression_filtered.csv 

p3_filtered_L6L7 <- filter_and_save(p3_comparison_L6L7, "L6", "L7")

Summary for L6 vs L7 :
Number of genes with p_val_adj = 0: 172 
Number of genes with p_val_adj = 1: 4312 
Number of significant genes (p_val_adj < 0.05): 6485 
Number of upregulated genes (avg_log2FC > 1): 1113 
Number of downregulated genes (avg_log2FC < -1): 929 
Filtered results saved to: comparison_L6_vs_L7_with_mean_expression_filtered.csv