Differential Expression Analysis - Filtering and Visualization-NEWUMAP

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-robj/1-MAST_with_batch_as_Covariate_with_meanExpression.csv", row.names = 1)
markers_wilcox <- read.csv("0-robj/2-Wilcox_min.pct_logfcT-0_with_meanExpression.csv", row.names = 1)
markers_mast_no_batch <- read.csv("0-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: 74 
Number of genes with p_val_adj = 1: 19337 
Number of significant genes (p_val_adj < 0.05): 5615 
Number of upregulated genes (avg_log2FC > 1): 15048 
Number of downregulated genes (avg_log2FC < -1): 2281 

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: 4208 
Number of genes with p_val_adj = 1: 10428 
Number of significant genes (p_val_adj < 0.05): 15043 
Number of upregulated genes (avg_log2FC > 1): 15048 
Number of downregulated genes (avg_log2FC < -1): 2281 

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: 4679 
Number of genes with p_val_adj = 1: 9855 
Number of significant genes (p_val_adj < 0.05): 15579 
Number of upregulated genes (avg_log2FC > 1): 15048 
Number of downregulated genes (avg_log2FC < -1): 2281 

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-robj/1-MAST_with_batch_as_Covariate_filtered.csv", row.names = TRUE)
 write.csv(filtered_markers_wilcox, file = "0-robj/2-Wilcox_min.pct_logfcT-Default_filtered.csv", row.names = TRUE)
 write.csv(filtered_markers_mast_no_batch, file = "0-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: 35 
Number of genes with p_val_adj = 1: 0 
Number of significant genes (p_val_adj < 0.05): 2077 
Number of upregulated genes (avg_log2FC > 1): 1956 
Number of downregulated genes (avg_log2FC < -1): 121 

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: 3765 
Number of genes with p_val_adj = 1: 0 
Number of significant genes (p_val_adj < 0.05): 4639 
Number of upregulated genes (avg_log2FC > 1): 4514 
Number of downregulated genes (avg_log2FC < -1): 125 

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: 3955 
Number of genes with p_val_adj = 1: 0 
Number of significant genes (p_val_adj < 0.05): 4639 
Number of upregulated genes (avg_log2FC > 1): 4514 
Number of downregulated genes (avg_log2FC < -1): 125 

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)

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

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