Differential Expression Analysis of Malignant CD4Tcells vs Control(Normal CD4 Tcells)-GSEA

1. load libraries

2. Perform DE analysis using Malignant_CD4Tcells_vs_Normal_CD4Tcells genes

Malignant_CD4Tcells_vs_Normal_CD4Tcells <- read.csv("0-imp_Robj/1-MAST_with_batch_as_Covariate_with_meanExpression.csv", header = T)

3. Create the EnhancedVolcano plot

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

EnhancedVolcano(Malignant_CD4Tcells_vs_Normal_CD4Tcells, 
                lab = Malignant_CD4Tcells_vs_Normal_CD4Tcells$gene,
                x = "avg_log2FC", 
                y = "p_val_adj",
                selectLab = c('EPCAM', 'BCAT1', 'KIR3DL2', 'FOXM1', 'TWIST1', 'TNFSF9', 
                              'CD80',  'IL1B', 'RPS4Y1', 
                              'IL7R', 'TCF7',  'MKI67', 'CD70', 
                              'IL2RA','TRBV6-2', 'TRBV10-3', 'TRBV4-2', 'TRBV9', 'TRBV7-9', 
                              'TRAV12-1', 'CD8B', 'FCGR3A', 'GNLY', 'FOXP3', 'SELL', 
                              'GIMAP1', 'RIPOR2', 'LEF1', 'HOXC9', 'SP5',
                              'CCL17', 'ETV4', 'THY1', 'FOXA2', 'ITGAD', 'S100P', 'TBX4', 
                              'ID1', 'XCL1', 'SOX2', 'CD27', 'CD28','PLS3','CD70','RAB25' , 'TRBV27', 'TRBV2'),
                title = "Malignant CD4 T cells(cell lines) vs normal CD4 T cells",
                xlab = bquote(~Log[2]~ 'fold change'),
                pCutoff = 0.05,
                FCcutoff = 1.5, 
                pointSize = 3.0,
                labSize = 5.0,
                boxedLabels = TRUE,
                colAlpha = 0.5,
                legendPosition = 'right',
                legendLabSize = 10,
                legendIconSize = 4.0,
                drawConnectors = TRUE,
                widthConnectors = 0.5,
                colConnectors = 'grey50',
                arrowheads = FALSE,
                max.overlaps = 30)

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

library(dplyr)
library(EnhancedVolcano)

# Assuming you have a data frame named Malignant_CD4Tcells_vs_Normal_CD4Tcells
# Filter genes based on lowest p-values but include all genes
filtered_genes <- Malignant_CD4Tcells_vs_Normal_CD4Tcells %>%
  arrange(p_val_adj, desc(abs(avg_log2FC)))

# Create the EnhancedVolcano plot with the filtered data
EnhancedVolcano(
  filtered_genes, 
  lab = ifelse(filtered_genes$p_val_adj <= 0.05 & abs(filtered_genes$avg_log2FC) >= 1.0, filtered_genes$gene, NA),
  x = "avg_log2FC", 
  y = "p_val_adj",
  title = "Malignant CD4 T cells(cell lines) vs normal CD4 T cells",
  pCutoff = 0.05,
  FCcutoff = 1.0,
  legendPosition = 'right', 
  labCol = 'black',
  labFace = 'bold',
  boxedLabels = FALSE,  # Set to FALSE to remove boxed labels
  pointSize = 3.0,
  labSize = 5.0,
  col = c('grey70', 'black', 'blue', 'red'),  # Customize point colors
  selectLab = filtered_genes$gene[filtered_genes$p_val_adj <= 0.05 & abs(filtered_genes$avg_log2FC) >= 1.0]  # Only label significant genes
)

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

EnhancedVolcano(
  filtered_genes, 
  lab = ifelse(filtered_genes$p_val_adj <= 0.05 & abs(filtered_genes$avg_log2FC) >= 1.0, filtered_genes$gene, NA),
  x = "avg_log2FC", 
  y = "p_val_adj",
  title = "Malignant CD4 T cells (cell lines) vs Normal CD4 T cells",
  subtitle = "Highlighting differentially expressed genes",
  pCutoff = 0.05,
  FCcutoff = 1.0,
  legendPosition = 'right',
  colAlpha = 0.8,  # Slight transparency for non-significant points
  col = c('grey70', 'black', 'blue', 'red'),  # Custom color scheme
  gridlines.major = TRUE,
  gridlines.minor = FALSE,
  selectLab = filtered_genes$gene[filtered_genes$p_val_adj <= 0.05 & abs(filtered_genes$avg_log2FC) >= 1.0]
) 

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

4. Enrichment Analysis-1

# Step-by-Step Guide for Gene Set Enrichment Analysis (GSEA) or Over-Representation Analysis (ORA)

# Load the necessary libraries
library(clusterProfiler)
library(org.Hs.eg.db)
library(enrichplot)
library(ReactomePA)

# Get upregulated genes based on log2FC and p-value thresholds
upregulated_genes <- Malignant_CD4Tcells_vs_Normal_CD4Tcells[Malignant_CD4Tcells_vs_Normal_CD4Tcells$avg_log2FC > 2 & Malignant_CD4Tcells_vs_Normal_CD4Tcells$p_val_adj < 0.05, ]

# Get downregulated genes based on log2FC and p-value thresholds
downregulated_genes <- Malignant_CD4Tcells_vs_Normal_CD4Tcells[Malignant_CD4Tcells_vs_Normal_CD4Tcells$avg_log2FC < -1 & Malignant_CD4Tcells_vs_Normal_CD4Tcells$p_val_adj < 0.05, ]

# Gene Ontology (GO) Enrichment Analysis
# GO enrichment for upregulated genes
go_up <- enrichGO(gene = upregulated_genes$gene, 
                  OrgDb = org.Hs.eg.db, 
                  keyType = "SYMBOL", 
                  ont = "BP",   # Biological Process (BP), Molecular Function (MF), Cellular Component (CC)
                  pAdjustMethod = "BH", 
                  pvalueCutoff = 0.05)

# GO enrichment for downregulated genes
go_down <- enrichGO(gene = downregulated_genes$gene, 
                    OrgDb = org.Hs.eg.db, 
                    keyType = "SYMBOL", 
                    ont = "BP", 
                    pAdjustMethod = "BH", 
                    pvalueCutoff = 0.05)

# Visualize the top enriched GO terms
dotplot(go_up, showCategory = 20, title = "GO Enrichment for Upregulated Genes")

dotplot(go_down, showCategory = 20, title = "GO Enrichment for Downregulated Genes")

# KEGG Pathway Enrichment
# Convert gene symbols to Entrez IDs for KEGG analysis
upregulated_entrez <- bitr(upregulated_genes$gene, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)$ENTREZID

'select()' returned 1:many mapping between keys and columns
Avis : 6.03% of input gene IDs are fail to map...

downregulated_entrez <- bitr(downregulated_genes$gene, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)$ENTREZID

'select()' returned 1:1 mapping between keys and columns
Avis : 12.5% of input gene IDs are fail to map...

# KEGG pathway enrichment for upregulated genes
kegg_up <- enrichKEGG(gene = upregulated_entrez, 
                      organism = "hsa", 
                      pvalueCutoff = 0.05)

# KEGG pathway enrichment for downregulated genes
kegg_down <- enrichKEGG(gene = downregulated_entrez, 
                        organism = "hsa", 
                        pvalueCutoff = 0.05)

# Visualize KEGG pathway results
dotplot(kegg_up, showCategory = 20, title = "KEGG Pathway Enrichment for Upregulated Genes")

dotplot(kegg_down, showCategory = 20, title = "KEGG Pathway Enrichment for Downregulated Genes")

# Reactome Pathway Enrichment
# Reactome pathway enrichment for upregulated genes
reactome_up <- enrichPathway(gene = upregulated_entrez, 
                             organism = "human", 
                             pvalueCutoff = 0.05)

# Reactome pathway enrichment for downregulated genes
reactome_down <- enrichPathway(gene = downregulated_entrez, 
                               organism = "human", 
                               pvalueCutoff = 0.05)

# Visualize Reactome pathways
dotplot(reactome_up, showCategory = 20, title = "Reactome Pathway Enrichment for Upregulated Genes")

dotplot(reactome_down, showCategory = 20, title = "Reactome Pathway Enrichment for Downregulated Genes")

# Gene Set Enrichment Analysis (GSEA)
# Create a ranked list of genes (log2FC as ranking metric)
gene_list <- Malignant_CD4Tcells_vs_Normal_CD4Tcells$avg_log2FC
names(gene_list) <- Malignant_CD4Tcells_vs_Normal_CD4Tcells$gene  # Use the $gene column for gene symbols
gene_list <- sort(gene_list, decreasing = TRUE)

# Convert gene symbols to Entrez IDs for GSEA
gene_df <- bitr(names(gene_list), fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)

'select()' returned 1:many mapping between keys and columns
Avis : 27.85% of input gene IDs are fail to map...

# Ensure the gene list matches the Entrez IDs
gene_list <- gene_list[names(gene_list) %in% gene_df$SYMBOL]

# Replace gene symbols with Entrez IDs
names(gene_list) <- gene_df$ENTREZID[match(names(gene_list), gene_df$SYMBOL)]

# Run GSEA using KEGG pathways
gsea_kegg <- gseKEGG(geneList = gene_list, 
                     organism = "hsa", 
                     pvalueCutoff = 0.05)

preparing geneSet collections...
GSEA analysis...
Avis : There are ties in the preranked stats (20.52% of the list).
The order of those tied genes will be arbitrary, which may produce unexpected results.Avis : There were 1 pathways for which P-values were not calculated properly due to unbalanced (positive and negative) gene-level statistic values. For such pathways pval, padj, NES, log2err are set to NA. You can try to increase the value of the argument nPermSimple (for example set it nPermSimple = 10000)Avis : For some of the pathways the P-values were likely overestimated. For such pathways log2err is set to NA.Avis : For some pathways, in reality P-values are less than 1e-10. You can set the `eps` argument to zero for better estimation.leading edge analysis...
done...

# Plot the GSEA results
gseaplot(gsea_kegg, geneSetID = 1, title = "Top KEGG Pathway")

# Extract the name of the top KEGG pathway
top_pathway <- gsea_kegg@result[1, "Description"]

# Plot GSEA with the top pathway's name as the title
gseaplot(gsea_kegg, geneSetID = 1, title = top_pathway)

NA
NA

4.2. Enrichment Analysis-2

# Load necessary libraries
library(clusterProfiler)
library(org.Hs.eg.db)
library(msigdbr)
library(enrichplot)

# Load Hallmark gene sets from msigdbr
hallmark_sets <- msigdbr(species = "Homo sapiens", category = "H")  # "H" is for Hallmark gene sets

# Get upregulated and downregulated genes based on log2 fold change and adjusted p-value
upregulated_genes <- Malignant_CD4Tcells_vs_Normal_CD4Tcells[Malignant_CD4Tcells_vs_Normal_CD4Tcells$avg_log2FC > 2 & Malignant_CD4Tcells_vs_Normal_CD4Tcells$p_val_adj < 0.05, ]
downregulated_genes <- Malignant_CD4Tcells_vs_Normal_CD4Tcells[Malignant_CD4Tcells_vs_Normal_CD4Tcells$avg_log2FC < -1 & Malignant_CD4Tcells_vs_Normal_CD4Tcells$p_val_adj < 0.05, ]

# Convert gene symbols to uppercase for consistency
upregulated_genes$gene <- toupper(upregulated_genes$gene)
downregulated_genes$gene <- toupper(downregulated_genes$gene)

# Check for overlap between your upregulated/downregulated genes and Hallmark gene sets
upregulated_in_hallmark <- intersect(upregulated_genes$gene, hallmark_sets$gene_symbol)
downregulated_in_hallmark <- intersect(downregulated_genes$gene, hallmark_sets$gene_symbol)

# Print the number of overlapping genes for both upregulated and downregulated genes
cat("Number of upregulated genes in Hallmark gene sets:", length(upregulated_in_hallmark), "\n")

Number of upregulated genes in Hallmark gene sets: 1016 

cat("Number of downregulated genes in Hallmark gene sets:", length(downregulated_in_hallmark), "\n")

Number of downregulated genes in Hallmark gene sets: 76 

# If there are genes to analyze, proceed with enrichment analysis
if (length(upregulated_in_hallmark) > 0) {
  # Perform enrichment analysis for upregulated genes using Hallmark gene sets
  hallmark_up <- enricher(gene = upregulated_in_hallmark, 
                          TERM2GENE = hallmark_sets[, c("gs_name", "gene_symbol")],  # Ensure TERM2GENE uses correct columns
                          pvalueCutoff = 0.05)
  # Check if results exist
  if (!is.null(hallmark_up) && nrow(hallmark_up) > 0) {
    # Visualize results if available
    dotplot(hallmark_up, showCategory = 20, title = "Hallmark Pathway Enrichment for Upregulated Genes")
  } else {
    cat("No significant enrichment found for upregulated genes.\n")
  }
} else {
  cat("No upregulated genes overlap with Hallmark gene sets.\n")
}

if (length(downregulated_in_hallmark) > 0) {
  # Perform enrichment analysis for downregulated genes using Hallmark gene sets
  hallmark_down <- enricher(gene = downregulated_in_hallmark, 
                            TERM2GENE = hallmark_sets[, c("gs_name", "gene_symbol")],  # Ensure TERM2GENE uses correct columns
                            pvalueCutoff = 0.05)
  # Check if results exist
  if (!is.null(hallmark_down) && nrow(hallmark_down) > 0) {
    # Visualize results if available
    dotplot(hallmark_down, showCategory = 20, title = "Hallmark Pathway Enrichment for Downregulated Genes")
  } else {
    cat("No significant enrichment found for downregulated genes.\n")
  }
} else {
  cat("No downregulated genes overlap with Hallmark gene sets.\n")
}

NA
NA

5. Bar PLOT

# Filter for significant pathways
top_kegg_up <- kegg_up@result[kegg_up@result$p.adjust < 0.05, ]
top_kegg_down <- kegg_down@result[kegg_down@result$p.adjust < 0.05, ]

# If there are not enough pathways, consider relaxing the threshold or checking the output
top_kegg_up <- top_kegg_up[order(-top_kegg_up$p.adjust), ][1:10, ]
top_kegg_down <- top_kegg_down[order(top_kegg_down$p.adjust), ][1:10, ]

# Combine into one data frame
top_pathways <- rbind(
  data.frame(Pathway = top_kegg_up$Description, p.adjust = top_kegg_up$p.adjust, Direction = "Upregulated"),
  data.frame(Pathway = top_kegg_down$Description, p.adjust = top_kegg_down$p.adjust, Direction = "Downregulated")
)

# Convert p.adjust to -log10(p.adjust) for visualization
top_pathways$neg_log10_p <- -log10(top_pathways$p.adjust)

# Create the barplot
ggplot(top_pathways, aes(x = reorder(Pathway, neg_log10_p), y = neg_log10_p, fill = Direction)) +
  geom_bar(stat = "identity", position = position_dodge()) +
  scale_fill_manual(values = c("Upregulated" = "red", "Downregulated" = "blue")) +
  coord_flip() +  # Flip the coordinates for better readability
  labs(title = "Top Significant Pathways",
       x = "Pathways",
       y = "-Log10 Adjusted P-Value") +
  theme_minimal() +
  theme(legend.title = element_blank())

# Load necessary library
library(ggplot2)

# Create the barplot
ggplot(top_pathways, aes(x = Pathway, y = neg_log10_p, fill = Direction)) +
  geom_bar(stat = "identity", position = "identity") +  # Use position = "identity"
  scale_fill_manual(values = c("Upregulated" = "red", "Downregulated" = "blue")) +
  coord_flip() +  # Flip the coordinates for better readability
  labs(title = "Top Significant Pathways",
       x = "Pathways",
       y = "-Log10 Adjusted P-Value") +
  theme_minimal() +
  theme(legend.title = element_blank())

NA
NA
NA

5. perform gene enrichment analysis and identify pathways

# Step-by-Step Gene Enrichment Analysis and Pathway Identification

# Load necessary packages
library(clusterProfiler)
library(org.Hs.eg.db)
library(enrichplot)
library(ggplot2)

# Identify significant genes based on log2FC and p-value criteria
significant_genes <- Malignant_CD4Tcells_vs_Normal_CD4Tcells[
    (Malignant_CD4Tcells_vs_Normal_CD4Tcells$avg_log2FC > 1.5 | 
     Malignant_CD4Tcells_vs_Normal_CD4Tcells$avg_log2FC < -1.5) & 
    Malignant_CD4Tcells_vs_Normal_CD4Tcells$p_val_adj < 0.05, ]$gene

# Convert gene symbols to Entrez IDs for downstream analysis
entrez_ids <- bitr(significant_genes, fromType = "SYMBOL", 
                   toType = "ENTREZID", 
                   OrgDb = org.Hs.eg.db)

'select()' returned 1:many mapping between keys and columns
Avis : 6.31% of input gene IDs are fail to map...

# Check for any NA values in conversion
entrez_ids <- na.omit(entrez_ids)

# Perform KEGG Pathway Enrichment
kegg_results <- enrichKEGG(gene = entrez_ids$ENTREZID,
                           organism = 'hsa',
                           pvalueCutoff = 0.05)

# View KEGG enrichment results
head(kegg_results)

# Perform GO (Biological Process) Enrichment
go_results <- enrichGO(gene = entrez_ids$ENTREZID,
                       OrgDb = org.Hs.eg.db,
                       ont = "BP", # Biological Process (BP)
                       pvalueCutoff = 0.05)

# View GO enrichment results
head(go_results)

# Dot plot for KEGG results
dotplot(kegg_results, showCategory = 10) + 
  ggtitle("KEGG Pathway Enrichment")

# Dot plot for GO results
dotplot(go_results, showCategory = 10) + 
  ggtitle("GO Biological Process Enrichment")

NA
NA

5. ggplot2 for Volcano

library(ggplot2)
library(ggrepel)

# Identify top and bottom genes
top_genes <- Malignant_CD4Tcells_vs_Normal_CD4Tcells[Malignant_CD4Tcells_vs_Normal_CD4Tcells$p_val_adj < 0.05 & Malignant_CD4Tcells_vs_Normal_CD4Tcells$avg_log2FC > 1.5, ]
bottom_genes <- Malignant_CD4Tcells_vs_Normal_CD4Tcells[Malignant_CD4Tcells_vs_Normal_CD4Tcells$p_val_adj < 0.05 & Malignant_CD4Tcells_vs_Normal_CD4Tcells$avg_log2FC < -1.5, ]

# Create a new column for color based on significance
Malignant_CD4Tcells_vs_Normal_CD4Tcells$color <- ifelse(Malignant_CD4Tcells_vs_Normal_CD4Tcells$avg_log2FC > 1.5, "Upregulated genes",
                                                   ifelse(Malignant_CD4Tcells_vs_Normal_CD4Tcells$avg_log2FC < -1.5, "Downregulated genes", "Nonsignificant"))

# Create a volcano plot
ggplot(Malignant_CD4Tcells_vs_Normal_CD4Tcells, aes(x = avg_log2FC, y = -log10(p_val_adj))) +
  geom_point(aes(color = color), alpha = 0.7, size = 2) +
  
  # Add labels for top and bottom genes
  geom_text_repel(data = top_genes, aes(label = gene), color = "black", vjust = 1, fontface = "bold") +
  geom_text_repel(data = bottom_genes, aes(label = gene), color = "black", vjust = -1, fontface = "bold") +
  
  # Customize labels and title
  labs(title = "Volcano Plot",
       x = "log2 Fold Change",
       y = "-log10(p-value)") +
  
  # Add significance threshold lines
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "black") +
  geom_vline(xintercept = c(-1.5, 2), linetype = "dashed", color = "black") +
  
  # Set colors for top and bottom genes
  scale_color_manual(values = c("Upregulated genes" = "red", "Downregulated genes" = "blue", "Nonsignificant" = "darkgrey")) +
  
  # Customize theme if needed
  theme_minimal()

NA
NA
NA
NA
NA