Filtering and Visualization(cell lines vs Control)

1. Load Libraries

2. Load Data

2.1 Load CSV Files with Mean Expression

# Define cell lines
cell_lines <- c("L1", "L2", "L3", "L4", "L5", "L6", "L7")

# Define filtering criteria
mean_expr_threshold <- 0.2
p_val_threshold <- 0.05
logfc_threshold <- 1

# Summary Function
summarize_markers <- function(markers, cell_line) {
  num_pval0 <- sum(markers$p_val_adj == 0)
  num_pval1 <- sum(markers$p_val_adj == 1)
  num_significant <- sum(markers$p_val_adj < p_val_threshold)
  num_upregulated <- sum(markers$avg_log2FC > logfc_threshold)
  num_downregulated <- sum(markers$avg_log2FC < -logfc_threshold)
  
  cat("\nSummary for", cell_line, ":\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")
}

# Volcano Plot Function
create_volcano_plot <- function(markers, cell_line) {
  volcano <- EnhancedVolcano(markers,
                             lab = markers$gene,
                             x = "avg_log2FC",
                             y = "p_val_adj",
                             title = paste0(cell_line, " vs Normal CD4+ T Cells"),
                             pCutoff = p_val_threshold,
                             FCcutoff = logfc_threshold,
                             pointSize = 2.5,
                             labSize = 3.0)
  
  print(volcano)
}

# Main Function to Process Each Cell Line
process_cell_line <- function(cell_line) {
  # Construct file name and read data
  file_name <- paste0("Updated_DE_Results_", cell_line, "_with_MeanExpr.csv")
  markers <- read.csv(file_name, row.names = 1)
  
  # Filter genes based on mean expression
  markers <- markers %>%
    filter(!(mean_expr_group1 < mean_expr_threshold & mean_expr_group2 < mean_expr_threshold))
  
  # Apply additional filters
  filtered_markers <- markers %>%
    filter(p_val_adj < p_val_threshold & abs(avg_log2FC) > logfc_threshold)
  
  # Summarize the filtered results
  summarize_markers(filtered_markers, cell_line)
  
  # Generate volcano plot
  create_volcano_plot(filtered_markers, cell_line)
}

# Loop through each cell line and apply the function
for (cell_line in cell_lines) {
  process_cell_line(cell_line)
}

Summary for L1 :
Number of genes with p_val_adj = 0: 1803 
Number of genes with p_val_adj = 1: 0 
Number of significant genes (p_val_adj < 0.05): 4129 
Number of upregulated genes (avg_log2FC > 1): 4013 
Number of downregulated genes (avg_log2FC < -1): 116 

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

Summary for L2 :
Number of genes with p_val_adj = 0: 2945 
Number of genes with p_val_adj = 1: 0 
Number of significant genes (p_val_adj < 0.05): 5465 
Number of upregulated genes (avg_log2FC > 1): 5324 
Number of downregulated genes (avg_log2FC < -1): 141 

Summary for L3 :
Number of genes with p_val_adj = 0: 2051 
Number of genes with p_val_adj = 1: 0 
Number of significant genes (p_val_adj < 0.05): 3725 
Number of upregulated genes (avg_log2FC > 1): 3432 
Number of downregulated genes (avg_log2FC < -1): 293 

Summary for L4 :
Number of genes with p_val_adj = 0: 3185 
Number of genes with p_val_adj = 1: 0 
Number of significant genes (p_val_adj < 0.05): 5584 
Number of upregulated genes (avg_log2FC > 1): 5419 
Number of downregulated genes (avg_log2FC < -1): 165 

Summary for L5 :
Number of genes with p_val_adj = 0: 2607 
Number of genes with p_val_adj = 1: 0 
Number of significant genes (p_val_adj < 0.05): 4688 
Number of upregulated genes (avg_log2FC > 1): 4471 
Number of downregulated genes (avg_log2FC < -1): 217 

Summary for L6 :
Number of genes with p_val_adj = 0: 1927 
Number of genes with p_val_adj = 1: 0 
Number of significant genes (p_val_adj < 0.05): 4074 
Number of upregulated genes (avg_log2FC > 1): 3915 
Number of downregulated genes (avg_log2FC < -1): 159 

Summary for L7 :
Number of genes with p_val_adj = 0: 2258 
Number of genes with p_val_adj = 1: 0 
Number of significant genes (p_val_adj < 0.05): 4455 
Number of upregulated genes (avg_log2FC > 1): 4277 
Number of downregulated genes (avg_log2FC < -1): 178 

2.2 Volcano Plot for MAST with SCT Assay

# Load necessary libraries
library(dplyr)
library(ggplot2)
library(EnhancedVolcano)

# Define cell lines
cell_lines <- c("L1", "L2", "L3", "L4", "L5", "L6", "L7")

# Define filtering criteria
mean_expr_threshold <- 0.2
p_val_threshold <- 0.05
logfc_threshold <- 1

# Define genes to highlight
highlight_genes <- c('KIR3DL2', 'TOX', 'TWIST1', 'CD52', 'CCR4', 'IL2RA', 
                     'CD7', 'DPP4', 'CCR7', 'GATA3', 'FOXP3', 'DNM3', 
                     'NCR1', 'CD70', 'AIRE', 'PDCD1 (PD-1)', 'CD274 (PD-L1)')

# Summary Function
summarize_markers <- function(markers, cell_line) {
  num_pval0 <- sum(markers$p_val_adj == 0)
  num_pval1 <- sum(markers$p_val_adj == 1)
  num_significant <- sum(markers$p_val_adj < p_val_threshold)
  num_upregulated <- sum(markers$avg_log2FC > logfc_threshold)
  num_downregulated <- sum(markers$avg_log2FC < -logfc_threshold)
  
  cat("\nSummary for", cell_line, ":\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")
}

# Volcano Plot Function
create_volcano_plot <- function(markers, cell_line) {
  volcano <- EnhancedVolcano(markers,
                             lab = markers$gene,
                             x = "avg_log2FC",
                             y = "p_val_adj",
                             title = paste0(cell_line, " vs Normal CD4+ T Cells"),
                             pCutoff = p_val_threshold,
                             FCcutoff = logfc_threshold,
                             pointSize = 2.5,
                             labSize = 3.0)
  print(volcano)
}

# Volcano Plot with Selected Genes
create_volcano_plot_highlight <- function(markers, cell_line) {
  volcano <- EnhancedVolcano(markers,
                             lab = markers$gene,
                             x = "avg_log2FC", 
                             y = "p_val_adj",
                             selectLab = highlight_genes,
                             title = paste0(cell_line, " vs Normal CD4+ T cells (Highlighted Genes)"),
                             xlab = bquote(~Log[2]~ 'fold change'),
                             pCutoff = 0.0001,
                             FCcutoff = logfc_threshold, 
                             pointSize = 3.0,
                             labSize = 3.0,
                             boxedLabels = TRUE,
                             colAlpha = 0.5,
                             legendPosition = 'right',
                             legendLabSize = 10,
                             legendIconSize = 4.0,
                             drawConnectors = TRUE,
                             widthConnectors = 0.5,
                             colConnectors = 'black',
                             arrowheads = FALSE,
                             max.overlaps = 30)
  print(volcano)
}

# Main Function to Process Each Cell Line
process_cell_line <- function(cell_line) {
  # Construct file name and read data
  file_name <- paste0("Updated_DE_Results_", cell_line, "_with_MeanExpr.csv")
  markers <- read.csv(file_name, row.names = 1)
  
  # Filter genes based on mean expression
  markers <- markers %>%
    filter(!(mean_expr_group1 < mean_expr_threshold & mean_expr_group2 < mean_expr_threshold))
  
  # Apply additional filters
  filtered_markers <- markers %>%
    filter(p_val_adj < p_val_threshold & abs(avg_log2FC) > logfc_threshold)
  
  # Summarize the filtered results
  summarize_markers(filtered_markers, cell_line)
  
  # Generate volcano plots
  create_volcano_plot(filtered_markers, cell_line)
  create_volcano_plot_highlight(filtered_markers, cell_line)
  
  # Save filtered results to CSV
  output_file <- paste0("Filtered_DE_Results_", cell_line, "_with_MeanExpr.csv")
  write.csv(filtered_markers, output_file, row.names = TRUE)
}

# Loop through each cell line and apply the function
for (cell_line in cell_lines) {
  process_cell_line(cell_line)
}

Summary for L1 :
Number of genes with p_val_adj = 0: 1803 
Number of genes with p_val_adj = 1: 0 
Number of significant genes (p_val_adj < 0.05): 4129 
Number of upregulated genes (avg_log2FC > 1): 4013 
Number of downregulated genes (avg_log2FC < -1): 116 

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

Summary for L2 :
Number of genes with p_val_adj = 0: 2945 
Number of genes with p_val_adj = 1: 0 
Number of significant genes (p_val_adj < 0.05): 5465 
Number of upregulated genes (avg_log2FC > 1): 5324 
Number of downregulated genes (avg_log2FC < -1): 141 

Summary for L3 :
Number of genes with p_val_adj = 0: 2051 
Number of genes with p_val_adj = 1: 0 
Number of significant genes (p_val_adj < 0.05): 3725 
Number of upregulated genes (avg_log2FC > 1): 3432 
Number of downregulated genes (avg_log2FC < -1): 293 

Summary for L4 :
Number of genes with p_val_adj = 0: 3185 
Number of genes with p_val_adj = 1: 0 
Number of significant genes (p_val_adj < 0.05): 5584 
Number of upregulated genes (avg_log2FC > 1): 5419 
Number of downregulated genes (avg_log2FC < -1): 165 

Summary for L5 :
Number of genes with p_val_adj = 0: 2607 
Number of genes with p_val_adj = 1: 0 
Number of significant genes (p_val_adj < 0.05): 4688 
Number of upregulated genes (avg_log2FC > 1): 4471 
Number of downregulated genes (avg_log2FC < -1): 217 

Summary for L6 :
Number of genes with p_val_adj = 0: 1927 
Number of genes with p_val_adj = 1: 0 
Number of significant genes (p_val_adj < 0.05): 4074 
Number of upregulated genes (avg_log2FC > 1): 3915 
Number of downregulated genes (avg_log2FC < -1): 159 

Summary for L7 :
Number of genes with p_val_adj = 0: 2258 
Number of genes with p_val_adj = 1: 0 
Number of significant genes (p_val_adj < 0.05): 4455 
Number of upregulated genes (avg_log2FC > 1): 4277 
Number of downregulated genes (avg_log2FC < -1): 178