Differential Expression Analysis of SS vs PBMC10X+PBMC

1. load libraries

2. Load Seurat Object

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 = "cell_line", 
        reduction = "umap.harmony",
        label.size = 3,
        repel = T,
        label = T)


DimPlot(All_samples_Merged,
        group.by = "seurat_clusters", 
        reduction = "umap.harmony",
        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.harmony",
        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.harmony",
        label.size = 3,
        repel = T,
        label = T, label.box = T)

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

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

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

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

table(All_samples_Merged$predicted.celltype.l2, All_samples_Merged$SCT_snn_res.0.5)

                   
                        0     1     2     3     4     5     6     7     8     9    10    11    12    13    14    15    16    17
  ASDC                  0     0     0     0     0     0     3     0     0     0     0     0     0     0     0     0     0     0
  B intermediate        0     1     0     0     0     0     3     0   664     1    17     4     0     2     2     2     0     0
  B memory              7     3     0     7    14     0     3   205   261     0     2     4     0    13     1     2     1     0
  B naive               0     1     0     0     0     0     4     0  1170     1     2    13     0     0     1     0     0     0
  CD14 Mono             3     1     8     0     0    20  2595    39     0     0   814   327     6     0     0     7     3     0
  CD16 Mono             1     0     0     0     0     0     4     0     0     0   120     1     0     0     0     0     0     0
  CD4 CTL               0     0     0     0     0     0     0     0     0    16     0     1     0     0     0     0     0     0
  CD4 Naive             0  1989     0     0     0     0     0     0     7     0     0     6     0     0    40     0     0     1
  CD4 Proliferating 10573     8  4839  6254  5489   412     1   738     2     1     0    23    48   490   132     1     0     0
  CD4 TCM             909  5922   335   115   212  2649    18  1293    65    38    32   398   579    43   256     3     0    11
  CD4 TEM               0    80     0     0     0     1     0     0     0    10     0     3     0     0     0     0     0     0
  CD8 Naive             0  1330     0     0     0     0     0     0     2     1     1    20     1     0    17     0     1     0
  CD8 Proliferating     0     0     0     0     0     0     0     2     0     0     0     0     0     0     0     0     0     0
  CD8 TCM               1   370     9     0     0     0     0     0     2    80     0     6     7     0     1     0     0     1
  CD8 TEM               0    64     6     0     1     0     1     6     1   305     0     0     0     0     7     0     0     0
  cDC1                  0     0     0     0     0     0    13     5     0     0     0     0     0     2     1    21     0     0
  cDC2                  0     0     1     1     2     0   123    45     0     0     1     0     0     3     0    53     0     0
  dnT                   0    31     0     0     0     1     0     7     2     0     3    17     0     0    20     0     0     1
  gdT                   0    15     0     0     0     0     0     0     0    78     0     0     0     0     0     0     0     0
  HSPC                 55     0     0  1429   258     3     6    11     1     1     0     3     6    38     7     4     7     5
  ILC                   0     0     0     0     0     0     0     0     1     4     0     1     0     0     1     0     0     0
  MAIT                  0     9     0     0     0     0     0     0     0   227     0     3     0     0     3     0     0     0
  NK                    0     0     0     0     0     0     0     0     1   523     1     7     0     0     2     0     0     0
  NK Proliferating     24     0  2818    15   254     0     0     6     0     2     0     4     0    11    32     0     1     0
  NK_CD56bright         0     0     0     0     0     0     0     0     0    14     0     0     0     0     2     0     0     0
  pDC                   0     0     0     0     0     0    56     0     0     0     0     0     0     0     0     0     0     0
  Plasmablast           0     0     0     0     0     0     5     0    13     0     0     0     0     0     1     0     0     0
  Platelet              0     0     0     0     0     0     1     0     0     0     1     0     0     0     0     0    30     0
  Treg                  0   258     1     0     0     0     0     9     5     0     0    40     0     0    39     1     0     0

3. Perform DE analysis using the FindMarkers or FindAllMarkers function

# Find markers for all clusters
All_markers_default <- FindAllMarkers(All_samples_Merged)

Calculating cluster L1
Calculating cluster L2
Calculating cluster L3
Calculating cluster L4
Calculating cluster L5
Calculating cluster L6
Calculating cluster L7
Calculating cluster PBMC
Calculating cluster PBMC_10x

All_markers_Default_min_diff  <- FindAllMarkers(All_samples_Merged, 
                              min.pct = 0.25,
                              logfc.threshold = 0.25)

Calculating cluster L1
Calculating cluster L2
Calculating cluster L3
Calculating cluster L4
Calculating cluster L5
Calculating cluster L6
Calculating cluster L7
Calculating cluster PBMC
Calculating cluster PBMC_10x

All_markers_0.5_min_diff <- FindAllMarkers(All_samples_Merged, 
                              min.pct = 0.25,
                              logfc.threshold = 0.25,
                              min.pct.diff= 0.5)

Calculating cluster L1
Calculating cluster L2
Calculating cluster L3
Calculating cluster L4
Calculating cluster L5
Calculating cluster L6
Calculating cluster L7
Calculating cluster PBMC
Calculating cluster PBMC_10x

Patient_cell_lines_vs_PBMC_Tcells <- FindMarkers(All_samples_Merged, 
                           ident.1 = c("2", "5", "12", "0", "3", "4", "7", "13"),
                           ident.2 = c("1", "9", "11", "14", "17")
                           )

Error in WhichCells.Seurat(object = object, idents = ident.1) : 
  Cannot find the following identities in the object: 2512034713

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)

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

Please cite:

S Xu, E Hu, Y Cai, Z Xie, X Luo, L Zhan, W Tang, Q Wang, B Liu, R Wang, W Xie, T Wu, L Xie, G Yu. Using
clusterProfiler to characterize multiomics data. Nature Protocols. 2024, doi:10.1038/s41596-024-01020-z

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.25.3 Learn more at https://yulab-smu.top/contribution-knowledge-mining/

Please cite:

Guangchuang Yu, Fei Li, Yide Qin, Xiaochen Bo, Yibo Wu and Shengqi Wang. GOSemSim: an R package for measuring
semantic similarity among GO terms and gene products. Bioinformatics. 2010, 26(7):976-978

library(ReactomePA)

ReactomePA v1.49.1 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.64% 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: 14.88% 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.13% 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.2% of the list).
The order of those tied genes will be arbitrary, which may produce unexpected results.Warning: There were 2 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: 11.08% 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