Differential Expression Analysis - Filtering and Visualization

1. Load Libraries

2. Load Data

2.1 Load CSV Files with Mean Expression

# Load the CSV files with mean expression data
markers_mast_batch <- read.csv("0-imp_Robj/TEST_MAST_with_batch_results1_with_meanExpression.csv", row.names = 1)
markers_wilcox <- read.csv("0-imp_Robj/TEST_Wilcox_without_batch_results2_with_meanExpression.csv", row.names = 1)
markers_mast_no_batch <- read.csv("0-imp_Robj/TEST_MAST_without_batch_results3_with_meanExpression.csv", row.names = 1)

3. Volcano Plots for All Genes

3.1 Volcano Plot for MAST with Batch Correction

EnhancedVolcano(markers_mast_batch,
                lab = rownames(markers_mast_batch),
                x = "avg_log2FC",
                y = "p_val_adj",
                title = "MAST with Batch Correction (All Genes)",
                pCutoff = 0.05,
                FCcutoff = 1.0)

3.2 Volcano Plot for Wilcoxon Test

EnhancedVolcano(markers_wilcox,
                lab = rownames(markers_wilcox),
                x = "avg_log2FC",
                y = "p_val_adj",
                title = "Wilcoxon Test (All Genes)",
                pCutoff = 0.05,
                FCcutoff = 1.0)

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

3.3 Volcano Plot for MAST without Batch Correction

EnhancedVolcano(markers_mast_no_batch,
                lab = rownames(markers_mast_no_batch),
                x = "avg_log2FC",
                y = "p_val_adj",
                title = "MAST without Batch Correction (All Genes)",
                pCutoff = 0.05,
                FCcutoff = 1.0)

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

4. Filter and Summarize Results

4.1 Apply Expression Filter First

# Apply the expression filter first
markers_mast_batch <- markers_mast_batch %>%
  filter(!(mean_expr_group1 < 0.2 & mean_expr_group2 < 0.2))

markers_wilcox <- markers_wilcox %>%
  filter(!(mean_expr_group1 < 0.2 & mean_expr_group2 < 0.2))

markers_mast_no_batch <- markers_mast_no_batch %>%
  filter(!(mean_expr_group1 < 0.2 & mean_expr_group2 < 0.2))

4.2 Apply Additional Filters for Statistical Significance and Fold Change

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

# Apply additional filters for markers_mast_batch
filtered_markers_mast_batch <- markers_mast_batch %>%
  filter(p_val_adj < p_val_threshold & abs(avg_log2FC) > logfc_threshold)

# Apply additional filters for markers_wilcox
filtered_markers_wilcox <- markers_wilcox %>%
  filter(p_val_adj < p_val_threshold & abs(avg_log2FC) > logfc_threshold)

# Apply additional filters for markers_mast_no_batch
filtered_markers_mast_no_batch <- markers_mast_no_batch %>%
  filter(p_val_adj < p_val_threshold & abs(avg_log2FC) > logfc_threshold)

# Save filtered results
write.csv(filtered_markers_mast_batch, file = "0-imp_Robj/1-MAST_with_batch_as_Covariate_filtered.csv", row.names = TRUE)
write.csv(filtered_markers_wilcox, file = "0-imp_Robj/2-Wilcox_min.pct_logfcT-0_filtered.csv", row.names = TRUE)
write.csv(filtered_markers_mast_no_batch, file = "0-imp_Robj/3-MAST_without_batch_as_Covariate_filtered.csv", row.names = TRUE)

4.3 Summary Function

summarize_markers <- function(markers) {
  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("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")
}

4.4 Summarize Markers1 (MAST with Batch Correction)

cat("Markers1 Summary (MAST with Batch Correction):\n")

Markers1 Summary (MAST with Batch Correction):

summarize_markers(filtered_markers_mast_batch)

Number of genes with p_val_adj = 0: 0 
Number of genes with p_val_adj = 1: 0 
Number of significant genes (p_val_adj < 0.05): 3466 
Number of upregulated genes (avg_log2FC > 1): 3387 
Number of downregulated genes (avg_log2FC < -1): 79 

4.5 Summarize Markers2 (Wilcoxon Test)

cat("Markers2 Summary (Wilcoxon Test):\n")

Markers2 Summary (Wilcoxon Test):

summarize_markers(filtered_markers_wilcox)

Number of genes with p_val_adj = 0: 3875 
Number of genes with p_val_adj = 1: 0 
Number of significant genes (p_val_adj < 0.05): 4535 
Number of upregulated genes (avg_log2FC > 1): 4409 
Number of downregulated genes (avg_log2FC < -1): 126 

4.6 Summarize Markers3 (MAST without Batch Correction)

cat("Markers3 Summary (MAST without Batch Correction):\n")

Markers3 Summary (MAST without Batch Correction):

summarize_markers(filtered_markers_mast_no_batch)

Number of genes with p_val_adj = 0: 4004 
Number of genes with p_val_adj = 1: 0 
Number of significant genes (p_val_adj < 0.05): 4535 
Number of upregulated genes (avg_log2FC > 1): 4409 
Number of downregulated genes (avg_log2FC < -1): 126 

5. Visualization of Filtered Results

5.1 Volcano Plot for Filtered MAST with Batch Correction

EnhancedVolcano(filtered_markers_mast_batch,
                lab = rownames(filtered_markers_mast_batch),
                x = "avg_log2FC",
                y = "p_val_adj",
                title = "Filtered MAST with Batch Correction",
                pCutoff = 0.05,
                FCcutoff = 1.0)