Sézary Syndrome Cell Line Analysis_NewUMAP_Wilcox_RNA_Assay

1. load libraries

#Differential Expression Analysis

2. load seurat object

#Differential Expression Analysis # 3. Pairwise Comparisons

library(Seurat)
library(dplyr)
library(tibble)
library(EnhancedVolcano)

# Extract normalized expression values for RNA assay
expression_data_RNA <- GetAssayData(All_samples_Merged, assay = "RNA", slot = "data")  # RNA-normalized counts

# Function to calculate mean expression for the specific comparison groups
calculate_mean_expression <- function(markers, ident1_cells, ident2_cells, expression_data) {
  group1_mean <- rowMeans(expression_data[, ident1_cells, drop = FALSE], na.rm = TRUE)
  group2_mean <- rowMeans(expression_data[, ident2_cells, drop = FALSE], na.rm = TRUE)
  
  markers <- markers %>%
    rownames_to_column("gene") %>%
    mutate(mean_expr_ident1 = group1_mean[gene],
           mean_expr_ident2 = group2_mean[gene])
  
  return(markers)
}

# Function to perform DE analysis, compute mean expression, generate volcano plot, and save results
perform_comparison_and_volcano <- function(All_samples_Merged, ident1, ident2, expression_data) {
  Idents(All_samples_Merged) <- "cell_line"
  
  # Get cells belonging to each identity
  ident1_cells <- WhichCells(All_samples_Merged, idents = ident1)
  ident2_cells <- WhichCells(All_samples_Merged, idents = ident2)
  
  # Perform differential expression analysis using RNA assay
  markers <- FindMarkers(All_samples_Merged, ident.1 = ident1, ident.2 = ident2, assay = "RNA")
  
  # Compute mean expression for the two compared groups using RNA-normalized counts
  markers_with_mean <- calculate_mean_expression(markers, ident1_cells, ident2_cells, expression_data)
  
  # Save results to CSV
  output_filename <- paste0("comparison_", ident1, "_vs_", ident2, "_with_mean_expression.csv")
  write.csv(markers_with_mean, file = output_filename, row.names = FALSE)
  
  # Create volcano plot
  volcano_plot <- EnhancedVolcano(markers_with_mean,
                                  lab = markers_with_mean$gene,
                                  x = 'avg_log2FC',
                                  y = 'p_val_adj',
                                  title = paste(ident1, 'vs', ident2),
                                  pCutoff = 0.05,
                                  FCcutoff = 1)
  
  print(volcano_plot)
  png(paste0("volcano_", ident1, "_vs_", ident2, ".png"), width = 12, height = 10, units = "in", res = 300)
  print(volcano_plot)
  dev.off()
  
  return(markers_with_mean)
}

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

Warning: 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)

Warning: 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)

Warning: 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)

Warning: 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)

Warning: 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