Differential Expression Analysis of SS vs PBMC10X+PBMC

1. load libraries

2. Load Seurat Object

load ("/home/bioinfo/0-imp_Robj/Harmony_integrated_All_samples_Merged_with_PBMC10x_with_harmony_clustering.Robj")
 

2. Initial Visualization

All_samples_Merged <- SetIdent(All_samples_Merged, value = "Harmony_snn_res.0.9")
  
DimPlot(All_samples_Merged,group.by = "cell_line", 
        reduction = "umap.harmony",
        label.size = 3,
        repel = T,
        label = T)

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

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

3. Perform DE analysis using the FindMarkers

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

Patient_cell_lines_vs_PBMC_Tcells <- FindMarkers(All_samples_Merged, 
                           ident.1 = c("3", "8", "10", "18", "1", "2", "13", "6", "16", "19", "4", "7", "9"),
                           ident.2 = c("0", "5", "14", "24", "20")
                           )



# Convert to data frame and add gene names as a new column
Patient_cell_lines_vs_PBMC_Tcells <- as.data.frame(Patient_cell_lines_vs_PBMC_Tcells)
Patient_cell_lines_vs_PBMC_Tcells$gene <- rownames(Patient_cell_lines_vs_PBMC_Tcells)

# Rearranging the columns for better readability (optional)
Patient_cell_lines_vs_PBMC_Tcells <- Patient_cell_lines_vs_PBMC_Tcells[, c("gene", "p_val", "avg_log2FC", "pct.1", "pct.2", "p_val_adj")]

#write.csv(Patient_cell_lines_vs_PBMC_Tcells, "Patient_cell_lines_vs_PBMC_Tcells.csv", row.names = FALSE)

Warning messages:
1: ggrepel: 25 unlabeled data points (too many overlaps). Consider increasing max.overlaps 
2: ggrepel: 25 unlabeled data points (too many overlaps). Consider increasing max.overlaps 
3: ggrepel: 15 unlabeled data points (too many overlaps). Consider increasing max.overlaps 
4: ggrepel: 13 unlabeled data points (too many overlaps). Consider increasing max.overlaps 
5: ggrepel: 8 unlabeled data points (too many overlaps). Consider increasing max.overlaps 
6: ggrepel: 7 unlabeled data points (too many overlaps). Consider increasing max.overlaps 
7: ggrepel: 7 unlabeled data points (too many overlaps). Consider increasing max.overlaps 

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.5, 
                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.5, 
                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 = Patient_cell_lines_vs_PBMC_Tcells$gene,
                x = "avg_log2FC", 
                y = "p_val_adj",
                selectLab = c('EPCAM', 'BCAT1', 'KIR3DL2', 'FOXM1', 'TWIST1', 'TNFSF9', 
                              'CD80',  'IL1B', 'TRBV7.6', 'TRBV5.4', 'TRBV12.4',
                              'UBE2C', 'PCLAF', 'TYMS', 'CDC20', 'RPS4Y1', 
                              'IL7R', 'TCF7', 'PTTG1', 'RRM2', '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 = "Sézary Cell Lines vs PBMC 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)

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)

clusterProfiler v4.14.0 Learn more at https://yulab-smu.top/contribution-knowledge-mining/

Please cite:

T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu.
clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. The Innovation. 2021,
2(3):100141

Attaching package: 'clusterProfiler'

The following object is masked from 'package:stats':

    filter

library(org.Hs.eg.db)

Loading required package: AnnotationDbi
Loading required package: stats4
Loading required package: BiocGenerics

Attaching package: 'BiocGenerics'

The following object is masked from 'package:gridExtra':

    combine

The following objects are masked from 'package:dplyr':

    combine, intersect, setdiff, union

The following object is masked from 'package:SeuratObject':

    intersect

The following objects are masked from 'package:stats':

    IQR, mad, sd, var, xtabs

The following objects are masked from 'package:base':

    anyDuplicated, aperm, append, as.data.frame, basename, cbind, colnames, dirname, do.call, duplicated,
    eval, evalq, Filter, Find, get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
    order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank, rbind, Reduce, rownames, sapply, saveRDS,
    setdiff, table, tapply, union, unique, unsplit, which.max, which.min

Loading required package: Biobase
Welcome to Bioconductor

    Vignettes contain introductory material; view with 'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages 'citation("pkgname")'.

Loading required package: IRanges
Loading required package: S4Vectors

Attaching package: 'S4Vectors'

The following object is masked from 'package:clusterProfiler':

    rename

The following objects are masked from 'package:dplyr':

    first, rename

The following object is masked from 'package:utils':

    findMatches

The following objects are masked from 'package:base':

    expand.grid, I, unname


Attaching package: 'IRanges'

The following object is masked from 'package:clusterProfiler':

    slice

The following objects are masked from 'package:dplyr':

    collapse, desc, slice

The following object is masked from 'package:sp':

    %over%


Attaching package: 'AnnotationDbi'

The following object is masked from 'package:clusterProfiler':

    select

The following object is masked from 'package:dplyr':

    select

library(enrichplot)

enrichplot v1.26.0 Learn more at https://yulab-smu.top/contribution-knowledge-mining/

Please cite:

Guangchuang Yu, Li-Gen Wang, and Qing-Yu He. ChIPseeker: an R/Bioconductor package for ChIP peak annotation,
comparison and visualization. Bioinformatics. 2015, 31(14):2382-2383

library(ReactomePA)

ReactomePA v1.50.0 Learn more at https://yulab-smu.top/contribution-knowledge-mining/

Please cite:

Guangchuang Yu, Qing-Yu He. ReactomePA: an R/Bioconductor package for reactome pathway analysis and
visualization. Molecular BioSystems. 2016, 12(2):477-479

# 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: 10.4% 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: 16.3% 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)

Reading KEGG annotation online: "https://rest.kegg.jp/link/hsa/pathway"...
Reading KEGG annotation online: "https://rest.kegg.jp/list/pathway/hsa"...

# 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.19% 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.27% of the list).
The order of those tied genes will be arbitrary, which may produce unexpected results.Warning: 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)Warning: For some of the pathways the P-values were likely overestimated. For such pathways log2err is set to NA.Warning: 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

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: 10.94% 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 > 1.5, ]
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 > 1.5, "Upregulated genes",
                                                   ifelse(Patient_cell_lines_vs_PBMC_Tcells$avg_log2FC < -1.5, "Downregulated genes", "Nonsignificant"))

# 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("Upregulated genes" = "red", "Downregulated genes" = "blue", "Nonsignificant" = "darkgrey")) +
  
  # Customize theme if needed
  theme_minimal()

NA
NA
NA
NA
NA