Differential Expression Analysis at 0.9

1. load libraries

2. Load Seurat Object

 load("/home/bioinfo/Cluster_to_Computer_Transfer_files_folder/All_Normal-PBMC_Abnormal-cellLines_T_cells_Merged_Annotated_UMAP_on_Clusters_to_USE.Robj")
 
All_samples_Merged <- SetIdent(All_samples_Merged, value = "cell_line")
  
DimPlot(All_samples_Merged,group.by = "cell_line", 
        reduction = "umap",
        label.size = 3,
        repel = T,
        label = T)

DimPlot(All_samples_Merged,
        group.by = "SCT_snn_res.0.9", 
        reduction = "umap",
        label.size = 3,
        repel = T,
        label = T)

#cluster_table <- table(Idents(All_samples_Merged))


DimPlot(All_samples_Merged, group.by = "predicted.celltype.l2", 
        reduction = "umap",
        label.size = 3,
        repel = T,
        label = T)

# library(clustree)
# clustree(All_samples_Merged, prefix = "SCT_snn_res.")


DimPlot(All_samples_Merged, group.by = "predicted.celltype.l1", 
        reduction = "umap",
        label.size = 3,
        repel = T,
        label = T, label.box = T)

DimPlot(All_samples_Merged, group.by = "predicted.celltype.l1", 
        reduction = "umap",
        label.size = 3,
        repel = T,
        label = F)

DimPlot(All_samples_Merged, group.by = "predicted.celltype.l2", 
        reduction = "umap",
        label.size = 3,
        repel = T,
        label = T, label.box = T)

DimPlot(All_samples_Merged, group.by = "predicted.celltype.l2", 
        reduction = "umap",
        label.size = 3,
        repel = T,
        label = F)

DimPlot(All_samples_Merged, group.by = "predicted.celltype.l2", 
        reduction = "umap",
        label.size = 3,
        repel = T,
        label = T, label.box = T)

# table(All_samples_Merged$predicted.celltype.l2, All_samples_Merged$SCT_snn_res.0.9)

3. Perform DE analysis using the FindMarkers or FindAllMarkers function

EnhancedVolcano(Patient_cell_lines_vs_PBMC_Tcells , 
                lab=rownames(Patient_cell_lines_vs_PBMC_Tcells),
                x ="avg_log2FC", 
                y ="p_val_adj",
                title = "Sézary Cell Lines vs PBMC T cells",
                pCutoff = 0.05,
                FCcutoff = 1, 
                legendPosition = 'right', 
                labCol = 'black',
                labFace = 'bold',
                boxedLabels = TRUE,
                pointSize = 3.0,
                labSize = 5.0, 
                drawConnectors = TRUE,
                widthConnectors = 0.25)

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

EnhancedVolcano(Patient_cell_lines_vs_PBMC_Tcells , 
                lab=rownames(Patient_cell_lines_vs_PBMC_Tcells),
                x ="avg_log2FC", 
                y ="p_val_adj",
                selectLab = c('EPCAM','BCAT1','KIR3DL2',
      'FOXM1','TWIST1','TNFSF9','CD80','CD7','IL1B', 'TRBV7.6','TRBV5.4','TRBV12.4'),
                title = "Sézary Cell Lines vs PBMC T cells",
                 xlab = bquote(~Log[2]~ 'fold change'),
                pCutoff = 0.05,
                FCcutoff = 1, 
                legendPosition = 'right', 
                legendLabSize = 14,
                legendIconSize = 4.0,
                labCol = 'black',
                labFace = 'bold',
                boxedLabels = TRUE,
                pointSize = 3.0,
                labSize = 5.0, 
                drawConnectors = TRUE,
                widthConnectors = 0.75,
                colConnectors = 'black')

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

EnhancedVolcano(Patient_cell_lines_vs_PBMC_Tcells, 
                lab = ifelse(Patient_cell_lines_vs_PBMC_Tcells$avg_log2FC > 1 & Patient_cell_lines_vs_PBMC_Tcells$p_val_adj < 0.05, 
                             rownames(Patient_cell_lines_vs_PBMC_Tcells), 
                             ""),  # Label only significant genes
                x = "avg_log2FC", 
                y = "p_val_adj",
                title = "Sézary Cell Lines vs PBMC T cells",
                pCutoff = 0.05,
                FCcutoff = 1, 
                legendPosition = 'right', 
                labCol = 'black',
                labFace = 'bold',
                boxedLabels = TRUE,
                pointSize = 3.0,
                labSize = 5.0, 
                drawConnectors = TRUE,
                widthConnectors = 0.25)

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

EnhancedVolcano(Patient_cell_lines_vs_PBMC_Tcells, 
                lab = ifelse((Patient_cell_lines_vs_PBMC_Tcells$avg_log2FC > 1.5 | Patient_cell_lines_vs_PBMC_Tcells$avg_log2FC < -1.5) & 
                             Patient_cell_lines_vs_PBMC_Tcells$p_val_adj < 0.05, 
                             rownames(Patient_cell_lines_vs_PBMC_Tcells), 
                             ""),  # Label only significant genes
                x = "avg_log2FC", 
                y = "p_val_adj",
                title = "Sézary Cell Lines vs PBMC T cells",
                pCutoff = 0.05,
                FCcutoff = 1, 
                legendPosition = 'right', 
                labCol = 'black',
                labFace = 'bold',
                boxedLabels = TRUE,
                pointSize = 3.0,
                labSize = 5.0, 
                drawConnectors = TRUE,
                widthConnectors = 0.25)

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

All_samples_Merged <- SetIdent(All_samples_Merged, value = "SCT_snn_res.0.9")

L2_Thesholds <- FindMarkers(All_samples_Merged, ident.1 = "4", ident.2 = "9", min.pct = 0.10, thresh.use = 0.25)

EnhancedVolcano(L2_Thesholds , 
                lab=rownames(L2_Thesholds),
                x ="avg_log2FC", 
                y ="p_val_adj",
                title = "4_vs_9",
                pCutoff = 0.05,
                FCcutoff = 1, 
                legendPosition = 'right', 
                labCol = 'black',
                labFace = 'bold',
                boxedLabels = TRUE,
                pointSize = 3.0,
                labSize = 3.0, 
                drawConnectors = FALSE,
                widthConnectors = 0.75)

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

4. Enrichment Analysis

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

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


# Get upregulated genes
upregulated_genes <- rownames(Patient_cell_lines_vs_PBMC_Tcells[Patient_cell_lines_vs_PBMC_Tcells$avg_log2FC > 1.5 & Patient_cell_lines_vs_PBMC_Tcells$p_val_adj < 0.05, ])

# Get downregulated genes
downregulated_genes <- rownames(Patient_cell_lines_vs_PBMC_Tcells[Patient_cell_lines_vs_PBMC_Tcells$avg_log2FC < -1.5 & Patient_cell_lines_vs_PBMC_Tcells$p_val_adj < 0.05, ])

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

# GO enrichment for downregulated genes
go_down <- enrichGO(gene = downregulated_genes, 
                    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, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)$ENTREZID

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

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

'select()' returned 1:1 mapping between keys and columns
Warning: 8.54% 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) (Optional)
# Create a ranked list of genes
gene_list <- Patient_cell_lines_vs_PBMC_Tcells$avg_log2FC
names(gene_list) <- rownames(Patient_cell_lines_vs_PBMC_Tcells)
gene_list <- sort(gene_list, decreasing = TRUE)

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

'select()' returned 1:many mapping between keys and columns
Warning: 11.26% 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
gsea_kegg <- gseKEGG(geneList = gene_list, 
                     organism = "hsa", 
                     pvalueCutoff = 0.05)

using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).

preparing geneSet collections...
GSEA analysis...
Warning: There are ties in the preranked stats (0.82% of the list).
The order of those tied genes will be arbitrary, which may produce unexpected results.Warning: There were 7 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)Warning: For some of the pathways the P-values were likely overestimated. For such pathways log2err is set to NA.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

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

# Load the packages
library(clusterProfiler)
library(org.Hs.eg.db)

# Assuming `significant_genes` is your list of significant gene symbols
significant_genes <- rownames(Patient_cell_lines_vs_PBMC_Tcells[
    (Patient_cell_lines_vs_PBMC_Tcells$avg_log2FC > 1.5 | 
     Patient_cell_lines_vs_PBMC_Tcells$avg_log2FC < -1.5) & 
    Patient_cell_lines_vs_PBMC_Tcells$p_val_adj < 0.05, ])

entrez_ids <- bitr(significant_genes, fromType = "SYMBOL", 
                   toType = "ENTREZID", 
                   OrgDb = org.Hs.eg.db)

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

kegg_results <- enrichKEGG(gene = entrez_ids$ENTREZID,
                           organism = 'hsa',
                           pvalueCutoff = 0.05)

# View results
head(kegg_results)

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

# View 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 <- Patient_cell_lines_vs_PBMC_Tcells[Patient_cell_lines_vs_PBMC_Tcells$p_val_adj < 0.05 & Patient_cell_lines_vs_PBMC_Tcells$avg_log2FC > 2, ]
bottom_genes <- Patient_cell_lines_vs_PBMC_Tcells[Patient_cell_lines_vs_PBMC_Tcells$p_val_adj < 0.05 & Patient_cell_lines_vs_PBMC_Tcells$avg_log2FC < -1.5, ]

# Create a new column for color based on significance
Patient_cell_lines_vs_PBMC_Tcells$color <- ifelse(Patient_cell_lines_vs_PBMC_Tcells$avg_log2FC > 2, "Top Gene",
                                                   ifelse(Patient_cell_lines_vs_PBMC_Tcells$avg_log2FC < -1.5, "Bottom Gene", "Neutral"))

# Create a volcano plot
ggplot(Patient_cell_lines_vs_PBMC_Tcells, 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("Top Gene" = "red", "Bottom Gene" = "blue", "Neutral" = "darkgrey")) +
  
  # Customize theme if needed
  theme_minimal()

NA
NA
NA
NA
NA

library(ggplot2)
library(ggrepel)

# Identify top and bottom genes based on criteria
top_genes <- Patient_cell_lines_vs_PBMC_Tcells[Patient_cell_lines_vs_PBMC_Tcells$p_val_adj < 0.05 & Patient_cell_lines_vs_PBMC_Tcells$avg_log2FC > 2, ]
bottom_genes <- Patient_cell_lines_vs_PBMC_Tcells[Patient_cell_lines_vs_PBMC_Tcells$p_val_adj < 0.05 & Patient_cell_lines_vs_PBMC_Tcells$avg_log2FC < -1.5, ]

# Sort and select the most significant top and bottom genes
top_genes <- top_genes[order(top_genes$p_val_adj), ][1:100, ]  # Top 100 significant upregulated genes
bottom_genes <- bottom_genes[order(bottom_genes$p_val_adj), ][1:100, ]  # Top 100 significant downregulated genes

# Create a new column for color based on significance
Patient_cell_lines_vs_PBMC_Tcells$color <- ifelse(Patient_cell_lines_vs_PBMC_Tcells$avg_log2FC > 2, "Top Gene",
                                                   ifelse(Patient_cell_lines_vs_PBMC_Tcells$avg_log2FC < -1.5, "Bottom Gene", "Neutral"))

# Create a volcano plot
ggplot(Patient_cell_lines_vs_PBMC_Tcells, aes(x = avg_log2FC, y = -log10(p_val_adj))) +
  geom_point(aes(color = color), alpha = 0.7, size = 2) +
  
  # Add labels only for top 100 and bottom 100 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("Top Gene" = "red", "Bottom Gene" = "blue", "Neutral" = "darkgrey")) +
  
  # Customize theme if needed
  theme_minimal()

NA
NA
NA
NA
NA