Differential Expression Analysis - Filtering and Visualization(MAST on SCT and RNA

1. Load Libraries

2. Load Data

2.1 Load CSV Files with Mean Expression

# Load the CSV files with mean expression data
markers_mast_SCT <- read.csv("../0-robj/1-MAST_with_SCT_batch_patient_cellline_as_Covariate_with_meanExpression.csv", row.names = 1)

markers_mast_RNA <- read.csv("../0-robj/2-MAST_with_RNA_batch_patient_cellline_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 SCT Assay)

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

Markers1 Summary (MAST with Batch Correction):

summarize_markers(markers_mast_SCT)

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 Markers3 (MAST with RNA Assay)

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

Markers3 Summary (MAST without Batch Correction):

summarize_markers(markers_mast_RNA)

Number of genes with p_val_adj = 0: 83 
Number of genes with p_val_adj = 1: 26400 
Number of significant genes (p_val_adj < 0.05): 9221 
Number of upregulated genes (avg_log2FC > 1): 7334 
Number of downregulated genes (avg_log2FC < -1): 15426 

4. Volcano Plots for All Genes

4.1 Volcano Plot for MAST with SCT Assay

EnhancedVolcano(markers_mast_SCT,
                lab = markers_mast_SCT$gene,
                x = "avg_log2FC",
                y = "p_val_adj",
                title = "MAST with SCT assay (All Genes)",
                pCutoff = 0.05,
                FCcutoff = 1)

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

4.2 Volcano Plot for MAST with SCT Assay

EnhancedVolcano(markers_mast_SCT, 
                lab = markers_mast_SCT$gene,
                x = "avg_log2FC", 
                y = "p_val_adj",
                selectLab = c('KIR3DL2', 'TOX', 'TWIST1', 'CD52', 'CCR4', 'IL2RA', 'CD7', 'DPP4', 'CCR7', 'GATA3', 'FOXP3', 'DNM3', 'NCR1', 'CD70', 'AIRE', 'PDCD1 (PD-1)', 'CD274 (PD-L1)'),
                title = "Sézary CD4 T cells vs Normal CD4 T cells",
                xlab = bquote(~Log[2]~ 'fold change'),
                pCutoff = 0.0001,
                FCcutoff = 1, 
                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)

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

4.3 Volcano Plot for MAST with RNA Assay

EnhancedVolcano(markers_mast_RNA,
                lab = markers_mast_RNA$gene,
                x = "avg_log2FC",
                y = "p_val_adj",
                title = "MAST with RNA Assay (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.4 Volcano Plot for MAST with RNA Assay

EnhancedVolcano(markers_mast_RNA, 
                lab = markers_mast_RNA$gene,
                x = "avg_log2FC", 
                y = "p_val_adj",
                selectLab = c('KIR3DL2', 'TOX', 'TWIST1', 'CD52', 'CCR4', 'IL2RA', 'CD7', 'DPP4', 'CCR7', 'GATA3', 'FOXP3', 'DNM3', 'NCR1', 'CD70', 'AIRE', 'PDCD1 (PD-1)', 'CD274 (PD-L1)'),
                title = "Sézary CD4 T cells vs Normal CD4 T cells",
                xlab = bquote(~Log[2]~ 'fold change'),
                pCutoff = 0.0001,
                FCcutoff = 1, 
                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)

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_SCT <- markers_mast_SCT %>%
  filter(!(mean_expr_group1 < 0.2 & mean_expr_group2 < 0.2))


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

5.1.1 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_SCT <- markers_mast_SCT %>%
  filter(p_val_adj < p_val_threshold & abs(avg_log2FC) > logfc_threshold)



# Apply additional filters for markers_mast_no_batch
filtered_markers_mast_RNA <- markers_mast_RNA %>%
  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.2 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.3 Summarize Markers1 (MAST with Batch Correction)

cat("Markers1 Summary (markers_mast_SCT):\n")

Markers1 Summary (markers_mast_SCT):

summarize_markers(filtered_markers_mast_SCT)

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.4 Summarize Markers2 (markers_mast_RNA)

cat("Markers2 Summary (markers_mast_RNA):\n")

Markers2 Summary (markers_mast_RNA):

summarize_markers(filtered_markers_mast_RNA)

Number of genes with p_val_adj = 0: 16 
Number of genes with p_val_adj = 1: 0 
Number of significant genes (p_val_adj < 0.05): 1488 
Number of upregulated genes (avg_log2FC > 1): 791 
Number of downregulated genes (avg_log2FC < -1): 697 

6. Visualization of Filtered Results

6.1 Volcano Plot for markers_mast_SCT

EnhancedVolcano(filtered_markers_mast_SCT,
                lab = filtered_markers_mast_SCT$gene,
                x = "avg_log2FC",
                y = "p_val_adj",
                title = "MAST with SCT assay (All Genes)",
                pCutoff = 0.05,
                FCcutoff = 1)

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

6.2 Volcano Plot for MAST with SCT Assay

EnhancedVolcano(filtered_markers_mast_SCT, 
                lab = filtered_markers_mast_SCT$gene,
                x = "avg_log2FC", 
                y = "p_val_adj",
                selectLab = c('KIR3DL2', 'TOX', 'TWIST1', 'CD52', 'CCR4', 'IL2RA', 'CD7', 'DPP4', 'CCR7', 'GATA3', 'FOXP3', 'DNM3', 'NCR1', 'CD70', 'AIRE', 'PDCD1 (PD-1)', 'CD274 (PD-L1)'),
                title = "Sézary CD4 T cells vs Normal CD4 T cells",
                xlab = bquote(~Log[2]~ 'fold change'),
                pCutoff = 0.0001,
                FCcutoff = 1, 
                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)

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

6.3 Volcano Plot for MAST with RNA Assay

EnhancedVolcano(filtered_markers_mast_RNA,
                lab = filtered_markers_mast_RNA$gene,
                x = "avg_log2FC",
                y = "p_val_adj",
                title = "MAST with RNA Assay (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...

6.4 Volcano Plot for MAST with RNA Assay

EnhancedVolcano(filtered_markers_mast_RNA, 
                lab = filtered_markers_mast_RNA$gene,
                x = "avg_log2FC", 
                y = "p_val_adj",
                selectLab = c('KIR3DL2', 'TOX', 'TWIST1', 'CD52', 'CCR4', 'IL2RA', 'CD7', 'DPP4', 'CCR7', 'GATA3', 'FOXP3', 'DNM3', 'NCR1', 'CD70', 'AIRE', 'PDCD1 (PD-1)', 'CD274 (PD-L1)'),
                title = "Sézary CD4 T cells vs Normal CD4 T cells",
                xlab = bquote(~Log[2]~ 'fold change'),
                pCutoff = 0.0001,
                FCcutoff = 1, 
                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)

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