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/1-MAST_with_batch_as_Covariate_with_meanExpression.csv", row.names = 1)
markers_wilcox <- read.csv("0-imp_Robj/2-Wilcox_min.pct_logfcT-0_with_meanExpression.csv", row.names = 1)
markers_mast_no_batch <- read.csv("0-imp_Robj/3-MAST_without_batch_as_Covariate_with_meanExpression.csv", row.names = 1)

3. Summarize Results

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

3.2 Summarize Markers1 (MAST with Batch Correction)

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

Markers1 Summary (MAST with Batch Correction):

summarize_markers(markers_mast_batch)

Number of genes with p_val_adj = 0: 1 
Number of genes with p_val_adj = 1: 16799 
Number of significant genes (p_val_adj < 0.05): 8601 
Number of upregulated genes (avg_log2FC > 1): 14796 
Number of downregulated genes (avg_log2FC < -1): 2494 

3.3 Summary for Markers2 (Wilcoxon Test)

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

Markers2 Summary (Wilcoxon Test):

summarize_markers(markers_wilcox)

Number of genes with p_val_adj = 0: 4502 
Number of genes with p_val_adj = 1: 10049 
Number of significant genes (p_val_adj < 0.05): 15357 
Number of upregulated genes (avg_log2FC > 1): 14796 
Number of downregulated genes (avg_log2FC < -1): 2494 

3.4 Summary for Markers3 (MAST without Batch Correction)

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

Markers3 Summary (MAST without Batch Correction):

summarize_markers(markers_mast_no_batch)

Number of genes with p_val_adj = 0: 4957 
Number of genes with p_val_adj = 1: 9570 
Number of significant genes (p_val_adj < 0.05): 15882 
Number of upregulated genes (avg_log2FC > 1): 14796 
Number of downregulated genes (avg_log2FC < -1): 2494 

4. Volcano Plots for All Genes

4.1 Volcano Plot for MAST with Batch Correction

EnhancedVolcano(markers_mast_batch,
                lab = markers_mast_batch$gene,
                x = "avg_log2FC",
                y = "p_val_adj",
                title = "MAST with 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.2 Volcano Plot for Wilcoxon Test

EnhancedVolcano(markers_wilcox,
                lab = markers_wilcox$gene,
                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...

4.3 Volcano Plot for MAST without Batch Correction

EnhancedVolcano(markers_mast_no_batch,
                lab = markers_mast_no_batch$gene,
                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...

5. Filter and Summarize Results

5.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))

5.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-Default_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)

5.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")
}

5.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): 3133 
Number of upregulated genes (avg_log2FC > 1): 3042 
Number of downregulated genes (avg_log2FC < -1): 91 

5.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: 3933 
Number of genes with p_val_adj = 1: 0 
Number of significant genes (p_val_adj < 0.05): 4630 
Number of upregulated genes (avg_log2FC > 1): 4503 
Number of downregulated genes (avg_log2FC < -1): 127 

5.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: 4090 
Number of genes with p_val_adj = 1: 0 
Number of significant genes (p_val_adj < 0.05): 4630 
Number of upregulated genes (avg_log2FC > 1): 4503 
Number of downregulated genes (avg_log2FC < -1): 127 

6. Visualization of Filtered Results

6.1 Volcano Plot for Filtered MAST with Batch Correction

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

6.2 Volcano Plot for Wilcoxon Test

EnhancedVolcano(filtered_markers_wilcox,
                lab = filtered_markers_wilcox$gene,
                x = "avg_log2FC",
                y = "p_val_adj",
                title = "Filtered Wilcoxon Test",
                pCutoff = 0.05,
                FCcutoff = 1.0)

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

6.3 Volcano Plot for MAST without Batch Correction

EnhancedVolcano(filtered_markers_mast_no_batch,
                lab = filtered_markers_mast_no_batch$gene,
                x = "avg_log2FC",
                y = "p_val_adj",
                title = "Filtered MAST withrout Batch Correction",
                pCutoff = 0.05,
                FCcutoff = 1.0)

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