Cell_lines vs PBMC

1. load libraries

#Differential Expression Analysis # 2. load seurat object

#Load Seurat Object L7
load("../../../0-IMP-OBJECTS/Harmony_integrated_All_samples_Merged_with_PBMC10x_with_harmony_clustering.Robj")


All_samples_Merged

An object of class Seurat 
64169 features across 59355 samples within 6 assays 
Active assay: SCT (27417 features, 3000 variable features)
 3 layers present: counts, data, scale.data
 5 other assays present: RNA, ADT, prediction.score.celltype.l1, prediction.score.celltype.l2, prediction.score.celltype.l3
 6 dimensional reductions calculated: integrated_dr, ref.umap, pca, umap, harmony, umap.harmony

DimPlot(All_samples_Merged, reduction = "umap.harmony", group.by = "cell_line",label = T, label.box = T)

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

#Differential Expression Analysis

3. L1 vs PBMC

DefaultAssay(All_samples_Merged) <- "SCT"
Idents(All_samples_Merged) <- "cell_line"

# P1 vs PBMC-Tcells
L1_vs_PBMC <- FindMarkers(All_samples_Merged, 
                        ident.1 = "L1",  
                        ident.2 = c("PBMC","PBMC_10x"),      
                        assay = "SCT")

write.csv(L1_vs_PBMC, "New_comparison_L1_vs_PBMC.csv")

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


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



# Create volcano plot for P1 vs PBMC-Tcells
volcano_L1_vs_PBMC <- EnhancedVolcano(L1_vs_PBMC, 
                                    lab = rownames(L1_vs_PBMC),
                                    x = 'avg_log2FC',
                                    y = 'p_val_adj',
                                    title = 'L1_vs_PBMC',
                                    xlab = bquote(~Log[2]~ 'fold change'),
                                    pCutoff = 1e-100,
                                    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...

print(volcano_L1_vs_PBMC)
png("volcano_L1_vs_PBMC.png", width = 12, height = 10, units = "in", res = 300)

print(volcano_L1_vs_PBMC)
dev.off()

png 
  2 

volcano2_L1_vs_PBMC <- EnhancedVolcano(L1_vs_PBMC, 
                lab = rownames(L1_vs_PBMC),
                x = "avg_log2FC", 
                y = "p_val_adj",
                selectLab = c('EPCAM', 'KIR3DL2', 'FOXM1', 'TWIST1', 'TNFSF9', 
                              'CD80', 'FOS','PTPN6','NCR1','NCR2',
                              'PCLAF', 'KIR3DL1', 'IL4','ITGA6','CCL5',
                              'IL7R', 'TCF7', 'PTTG1', 'RRM2', 'MKI67', 'CD70', 
                              'IL2RA', 'FCGR3A', 'GNLY', 'FOXP3', 'SELL',  'LEF1',
                              'CCL17', 'THY1', 'CD27', 'CD28', 'CD7',
                              # Key Sézary syndrome genes
                              'PRF1', 'GZMB', 'NCR1', 'NFATC3', 
                              'KLRK1', 'LCK', 'KLRC1', 'KLRC2', 'TNF', 
                              'KIR3DL1','KIR3DL3','KIR3DL4', 'IFNG', 'IFNGR1', 'CD244', 'FASLG'),
                title = "L1_vs_PBMC",
                subtitle = "Sézary Syndrome Cell Lines",
                xlab = bquote(~Log[2]~ 'fold change'),
                pCutoff = 0.05,
                FCcutoff = 1.5, 
                pointSize = 3.0,
                labSize = 4.0,
                labFace = 'bold',
                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...

print(volcano2_L1_vs_PBMC)
png("volcano2_L1_vs_PBMC.png", width = 12, height = 10, units = "in", res = 300)

print(volcano2_L1_vs_PBMC)
dev.off()

png 
  2 

# Display top differentially expressed genes for each comparison
head(L1_vs_PBMC)

NA
NA

4. L2_vs_PBMC

DefaultAssay(All_samples_Merged) <- "SCT"
Idents(All_samples_Merged) <- "cell_line"

# P1 vs PBMC-Tcells
L2_vs_PBMC <- FindMarkers(All_samples_Merged, 
                        ident.1 = "L2",  
                        ident.2 = c("PBMC","PBMC_10x"),      
                        assay = "SCT")

write.csv(L2_vs_PBMC, "New_comparison_L2_vs_PBMC.csv")

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


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



# Create volcano plot for P1 vs PBMC-Tcells
volcano_L2_vs_PBMC <- EnhancedVolcano(L2_vs_PBMC, 
                                    lab = rownames(L2_vs_PBMC),
                                    x = 'avg_log2FC',
                                    y = 'p_val_adj',
                                    title = 'L2_vs_PBMC',
                                    xlab = bquote(~Log[2]~ 'fold change'),
                                    pCutoff = 1e-100,
                                    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...

print(volcano_L2_vs_PBMC)
png("volcano_L2_vs_PBMC.png", width = 12, height = 10, units = "in", res = 300)

print(volcano_L2_vs_PBMC)
dev.off()

png 
  2 

volcano2_L2_vs_PBMC <- EnhancedVolcano(L2_vs_PBMC, 
                lab = rownames(L2_vs_PBMC),
                x = "avg_log2FC", 
                y = "p_val_adj",
                selectLab = c('EPCAM', 'KIR3DL2', 'FOXM1', 'TWIST1', 'TNFSF9', 
                              'CD80', 'FOS','PTPN6','NCR1','NCR2',
                              'PCLAF', 'KIR3DL1', 'IL4','ITGA6','CCL5',
                              'IL7R', 'TCF7', 'PTTG1', 'RRM2', 'MKI67', 'CD70', 
                              'IL2RA', 'FCGR3A', 'GNLY', 'FOXP3', 'SELL',  'LEF1',
                              'CCL17', 'THY1', 'CD27', 'CD28', 'CD7',
                              # Key Sézary syndrome genes
                              'PRF1', 'GZMB', 'NCR1', 'NFATC3', 
                              'KLRK1', 'LCK', 'KLRC1', 'KLRC2', 'TNF', 
                              'KIR3DL1','KIR3DL3','KIR3DL4', 'IFNG', 'IFNGR1', 'CD244', 'FASLG'),
                title = "L2_vs_PBMC",
                subtitle = "Sézary Syndrome Cell Lines",
                xlab = bquote(~Log[2]~ 'fold change'),
                pCutoff = 0.05,
                FCcutoff = 1.5, 
                pointSize = 3.0,
                labSize = 4.0,
                labFace = 'bold',
                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...

print(volcano2_L2_vs_PBMC)
png("volcano2_L2_vs_PBMC.png", width = 12, height = 10, units = "in", res = 300)

print(volcano2_L2_vs_PBMC)
dev.off()

png 
  2 

# Display top differentially expressed genes for each comparison
head(L2_vs_PBMC)

NA
NA

5. L3_vs_PBMC

DefaultAssay(All_samples_Merged) <- "SCT"
Idents(All_samples_Merged) <- "cell_line"

# P1 vs PBMC-Tcells
L3_vs_PBMC <- FindMarkers(All_samples_Merged, 
                        ident.1 = "L3",  
                        ident.2 = c("PBMC","PBMC_10x"),      
                        assay = "SCT")

write.csv(L3_vs_PBMC, "New_comparison_L3_vs_PBMC.csv")

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


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



# Create volcano plot for P1 vs PBMC-Tcells
volcano_L3_vs_PBMC <- EnhancedVolcano(L3_vs_PBMC, 
                                    lab = rownames(L3_vs_PBMC),
                                    x = 'avg_log2FC',
                                    y = 'p_val_adj',
                                    title = 'L3_vs_PBMC',
                                    xlab = bquote(~Log[2]~ 'fold change'),
                                    pCutoff = 1e-100,
                                    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...

print(volcano_L3_vs_PBMC)
png("volcano_L3_vs_PBMC.png", width = 12, height = 10, units = "in", res = 300)

print(volcano_L3_vs_PBMC)
dev.off()

png 
  2 

volcano2_L3_vs_PBMC <- EnhancedVolcano(L3_vs_PBMC, 
                lab = rownames(L3_vs_PBMC),
                x = "avg_log2FC", 
                y = "p_val_adj",
                selectLab = c('EPCAM', 'KIR3DL2', 'FOXM1', 'TWIST1', 'TNFSF9', 
                              'CD80', 'FOS','PTPN6','NCR1','NCR2',
                              'PCLAF', 'KIR3DL1', 'IL4','ITGA6','CCL5',
                              'IL7R', 'TCF7', 'PTTG1', 'RRM2', 'MKI67', 'CD70', 
                              'IL2RA', 'FCGR3A', 'GNLY', 'FOXP3', 'SELL',  'LEF1',
                              'CCL17', 'THY1', 'CD27', 'CD28', 'CD7',
                              # Key Sézary syndrome genes
                              'PRF1', 'GZMB', 'NCR1', 'NFATC3', 
                              'KLRK1', 'LCK', 'KLRC1', 'KLRC2', 'TNF', 
                              'KIR3DL1','KIR3DL3','KIR3DL4', 'IFNG', 'IFNGR1', 'CD244', 'FASLG'),
                title = "L3_vs_PBMC",
                subtitle = "Sézary Syndrome Cell Lines",
                xlab = bquote(~Log[2]~ 'fold change'),
                pCutoff = 0.05,
                FCcutoff = 1.5, 
                pointSize = 3.0,
                labSize = 4.0,
                labFace = 'bold',
                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...

print(volcano2_L3_vs_PBMC)
png("volcano2_L3_vs_PBMC.png", width = 12, height = 10, units = "in", res = 300)

print(volcano2_L3_vs_PBMC)
dev.off()

png 
  2 

# Display top differentially expressed genes for each comparison
head(L3_vs_PBMC)

NA
NA

6. L4_vs_PBMC

DefaultAssay(All_samples_Merged) <- "SCT"
Idents(All_samples_Merged) <- "cell_line"

# P1 vs PBMC-Tcells
L4_vs_PBMC <- FindMarkers(All_samples_Merged, 
                        ident.1 = "L4",  
                        ident.2 = c("PBMC","PBMC_10x"),      
                        assay = "SCT")

write.csv(L4_vs_PBMC, "New_comparison_L4_vs_PBMC.csv")

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


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



# Create volcano plot for P1 vs PBMC-Tcells
volcano_L4_vs_PBMC <- EnhancedVolcano(L4_vs_PBMC, 
                                    lab = rownames(L4_vs_PBMC),
                                    x = 'avg_log2FC',
                                    y = 'p_val_adj',
                                    title = 'L4_vs_PBMC',
                                    xlab = bquote(~Log[2]~ 'fold change'),
                                    pCutoff = 1e-100,
                                    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...

print(volcano_L4_vs_PBMC)

Registered S3 method overwritten by 'rmarkdown':
  method         from
  print.paged_df     

png("volcano_L4_vs_PBMC.png", width = 12, height = 10, units = "in", res = 300)

print(volcano_L4_vs_PBMC)
dev.off()

png 
  2 

volcano2_L4_vs_PBMC <- EnhancedVolcano(L4_vs_PBMC, 
                lab = rownames(L4_vs_PBMC),
                x = "avg_log2FC", 
                y = "p_val_adj",
                selectLab = c('EPCAM', 'KIR3DL2', 'FOXM1', 'TWIST1', 'TNFSF9', 
                              'CD80', 'FOS','PTPN6','NCR1','NCR2',
                              'PCLAF', 'KIR3DL1', 'IL4','ITGA6','CCL5',
                              'IL7R', 'TCF7', 'PTTG1', 'RRM2', 'MKI67', 'CD70', 
                              'IL2RA', 'FCGR3A', 'GNLY', 'FOXP3', 'SELL',  'LEF1',
                              'CCL17', 'THY1', 'CD27', 'CD28', 'CD7',
                              # Key Sézary syndrome genes
                              'PRF1', 'GZMB', 'NCR1', 'NFATC3', 
                              'KLRK1', 'LCK', 'KLRC1', 'KLRC2', 'TNF', 
                              'KIR3DL1','KIR3DL3','KIR3DL4', 'IFNG', 'IFNGR1', 'CD244', 'FASLG'),
                title = "L4_vs_PBMC",
                subtitle = "Sézary Syndrome Cell Lines",
                xlab = bquote(~Log[2]~ 'fold change'),
                pCutoff = 0.05,
                FCcutoff = 1.5, 
                pointSize = 3.0,
                labSize = 4.0,
                labFace = 'bold',
                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...

print(volcano2_L4_vs_PBMC)
png("volcano2_L4_vs_PBMC.png", width = 12, height = 10, units = "in", res = 300)

print(volcano2_L4_vs_PBMC)
dev.off()

png 
  2 

# Display top differentially expressed genes for each comparison
head(L4_vs_PBMC)

NA
NA

7. L5_vs_PBMC

DefaultAssay(All_samples_Merged) <- "SCT"
Idents(All_samples_Merged) <- "cell_line"

# P1 vs PBMC-Tcells
L5_vs_PBMC <- FindMarkers(All_samples_Merged, 
                        ident.1 = "L5",  
                        ident.2 = c("PBMC","PBMC_10x"),      
                        assay = "SCT")

write.csv(L5_vs_PBMC, "New_comparison_L5_vs_PBMC.csv")

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


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



# Create volcano plot for P1 vs PBMC-Tcells
volcano_L5_vs_PBMC <- EnhancedVolcano(L5_vs_PBMC, 
                                    lab = rownames(L5_vs_PBMC),
                                    x = 'avg_log2FC',
                                    y = 'p_val_adj',
                                    title = 'L5_vs_PBMC',
                                    xlab = bquote(~Log[2]~ 'fold change'),
                                    pCutoff = 1e-100,
                                    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...

print(volcano_L5_vs_PBMC)
png("volcano_L5_vs_PBMC.png", width = 12, height = 10, units = "in", res = 300)

print(volcano_L5_vs_PBMC)
dev.off()

png 
  2 

volcano2_L5_vs_PBMC <- EnhancedVolcano(L5_vs_PBMC, 
                lab = rownames(L5_vs_PBMC),
                x = "avg_log2FC", 
                y = "p_val_adj",
                selectLab = c('EPCAM', 'KIR3DL2', 'FOXM1', 'TWIST1', 'TNFSF9', 
                              'CD80', 'FOS','PTPN6','NCR1','NCR2',
                              'PCLAF', 'KIR3DL1', 'IL4','ITGA6','CCL5',
                              'IL7R', 'TCF7', 'PTTG1', 'RRM2', 'MKI67', 'CD70', 
                              'IL2RA', 'FCGR3A', 'GNLY', 'FOXP3', 'SELL',  'LEF1',
                              'CCL17', 'THY1', 'CD27', 'CD28', 'CD7',
                              # Key Sézary syndrome genes
                              'PRF1', 'GZMB', 'NCR1', 'NFATC3', 
                              'KLRK1', 'LCK', 'KLRC1', 'KLRC2', 'TNF', 
                              'KIR3DL1','KIR3DL3','KIR3DL4', 'IFNG', 'IFNGR1', 'CD244', 'FASLG'),
                title = "L5_vs_PBMC",
                subtitle = "Sézary Syndrome Cell Lines",
                xlab = bquote(~Log[2]~ 'fold change'),
                pCutoff = 0.05,
                FCcutoff = 1.5, 
                pointSize = 3.0,
                labSize = 4.0,
                labFace = 'bold',
                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...

print(volcano2_L5_vs_PBMC)
png("volcano2_L5_vs_PBMC.png", width = 12, height = 10, units = "in", res = 300)

print(volcano2_L5_vs_PBMC)
dev.off()

png 
  2 

# Display top differentially expressed genes for each comparison
head(L5_vs_PBMC)

NA
NA

8. L6_vs_PBMC

DefaultAssay(All_samples_Merged) <- "SCT"
Idents(All_samples_Merged) <- "cell_line"

# P1 vs PBMC-Tcells
L6_vs_PBMC <- FindMarkers(All_samples_Merged, 
                        ident.1 = "L6",  
                        ident.2 = c("PBMC","PBMC_10x"),      
                        assay = "SCT")

write.csv(L6_vs_PBMC, "New_comparison_L6_vs_PBMC.csv")

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


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



# Create volcano plot for P1 vs PBMC-Tcells
volcano_L6_vs_PBMC <- EnhancedVolcano(L6_vs_PBMC, 
                                    lab = rownames(L6_vs_PBMC),
                                    x = 'avg_log2FC',
                                    y = 'p_val_adj',
                                    title = 'L6_vs_PBMC',
                                    xlab = bquote(~Log[2]~ 'fold change'),
                                    pCutoff = 1e-100,
                                    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...

print(volcano_L6_vs_PBMC)
png("volcano_L6_vs_PBMC.png", width = 12, height = 10, units = "in", res = 300)

print(volcano_L6_vs_PBMC)
dev.off()

png 
  2 

volcano2_L6_vs_PBMC <- EnhancedVolcano(L6_vs_PBMC, 
                lab = rownames(L6_vs_PBMC),
                x = "avg_log2FC", 
                y = "p_val_adj",
                selectLab = c('EPCAM', 'KIR3DL2', 'FOXM1', 'TWIST1', 'TNFSF9', 
                              'CD80', 'FOS','PTPN6','NCR1','NCR2',
                              'PCLAF', 'KIR3DL1', 'IL4','ITGA6','CCL5',
                              'IL7R', 'TCF7', 'PTTG1', 'RRM2', 'MKI67', 'CD70', 
                              'IL2RA', 'FCGR3A', 'GNLY', 'FOXP3', 'SELL',  'LEF1',
                              'CCL17', 'THY1', 'CD27', 'CD28', 'CD7',
                              # Key Sézary syndrome genes
                              'PRF1', 'GZMB', 'NCR1', 'NFATC3', 
                              'KLRK1', 'LCK', 'KLRC1', 'KLRC2', 'TNF', 
                              'KIR3DL1','KIR3DL3','KIR3DL4', 'IFNG', 'IFNGR1', 'CD244', 'FASLG'),
                title = "L6_vs_PBMC",
                subtitle = "Sézary Syndrome Cell Lines",
                xlab = bquote(~Log[2]~ 'fold change'),
                pCutoff = 0.05,
                FCcutoff = 1.5, 
                pointSize = 3.0,
                labSize = 4.0,
                labFace = 'bold',
                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...

print(volcano2_L6_vs_PBMC)
png("volcano2_L6_vs_PBMC.png", width = 12, height = 10, units = "in", res = 300)

print(volcano2_L6_vs_PBMC)
dev.off()

png 
  2 

# Display top differentially expressed genes for each comparison
head(L6_vs_PBMC)

NA
NA

9. L7_vs_PBMC

DefaultAssay(All_samples_Merged) <- "SCT"
Idents(All_samples_Merged) <- "cell_line"

# P1 vs PBMC-Tcells
L7_vs_PBMC <- FindMarkers(All_samples_Merged, 
                        ident.1 = "L7",  
                        ident.2 = c("PBMC","PBMC_10x"),      
                        assay = "SCT")

write.csv(L7_vs_PBMC, "New_comparison_L7_vs_PBMC.csv")

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


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



# Create volcano plot for P1 vs PBMC-Tcells
volcano_L7_vs_PBMC <- EnhancedVolcano(L7_vs_PBMC, 
                                    lab = rownames(L7_vs_PBMC),
                                    x = 'avg_log2FC',
                                    y = 'p_val_adj',
                                    title = 'L7_vs_PBMC',
                                    xlab = bquote(~Log[2]~ 'fold change'),
                                    pCutoff = 1e-100,
                                    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...

print(volcano_L7_vs_PBMC)
png("volcano_L7_vs_PBMC.png", width = 12, height = 10, units = "in", res = 300)

print(volcano_L7_vs_PBMC)
dev.off()

png 
  2 

volcano2_L7_vs_PBMC <- EnhancedVolcano(L7_vs_PBMC, 
                lab = rownames(L7_vs_PBMC),
                x = "avg_log2FC", 
                y = "p_val_adj",
                selectLab = c('EPCAM', 'KIR3DL2', 'FOXM1', 'TWIST1', 'TNFSF9', 
                              'CD80', 'FOS','PTPN6','NCR1','NCR2',
                              'PCLAF', 'KIR3DL1', 'IL4','ITGA6','CCL5',
                              'IL7R', 'TCF7', 'PTTG1', 'RRM2', 'MKI67', 'CD70', 
                              'IL2RA', 'FCGR3A', 'GNLY', 'FOXP3', 'SELL',  'LEF1',
                              'CCL17', 'THY1', 'CD27', 'CD28', 'CD7',
                              # Key Sézary syndrome genes
                              'PRF1', 'GZMB', 'NCR1', 'NFATC3', 
                              'KLRK1', 'LCK', 'KLRC1', 'KLRC2', 'TNF', 
                              'KIR3DL1','KIR3DL3','KIR3DL4', 'IFNG', 'IFNGR1', 'CD244', 'FASLG'),
                title = "L7_vs_PBMC",
                subtitle = "Sézary Syndrome Cell Lines",
                xlab = bquote(~Log[2]~ 'fold change'),
                pCutoff = 0.05,
                FCcutoff = 1.5, 
                pointSize = 3.0,
                labSize = 4.0,
                labFace = 'bold',
                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...

print(volcano2_L7_vs_PBMC)
png("volcano2_L7_vs_PBMC.png", width = 12, height = 10, units = "in", res = 300)

print(volcano2_L7_vs_PBMC)
dev.off()

png 
  2 

# Display top differentially expressed genes for each comparison
head(L7_vs_PBMC)

NA
NA

6. Enrichment Analysis

library(clusterProfiler)
library(org.Hs.eg.db)
library(enrichplot)

perform_go_enrichment <- function(gene_list, gene_universe, title) {
  ego <- enrichGO(gene = gene_list,
                  universe = gene_universe,
                  OrgDb = org.Hs.eg.db,
                  keyType = "SYMBOL",
                  ont = "BP",
                  pAdjustMethod = "BH",
                  qvalueCutoff = 0.05,
                  readable = TRUE)
  
  if (nrow(ego@result) == 0) {
    warning(paste("No enriched GO terms found for", title))
    return(NULL)
  }
  
  p <- dotplot(ego, showCategory = 10, title = paste("GO -", title)) +
    theme(axis.text.y = element_text(size = 8))
  
  print(p)
  png(paste0("GO_enrichment_", gsub(" ", "_", title), ".png"), width = 12, height = 8, units = "in", res = 300)
  print(p)
  dev.off()
  
  return(ego)
}

perform_kegg_enrichment <- function(gene_list, gene_universe, title) {
  # Convert gene symbols to Entrez IDs
  entrez_ids <- bitr(gene_list, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)
  universe_entrez <- bitr(gene_universe, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)
  
  print(paste("Number of input genes:", length(gene_list)))
  print(paste("Number of input genes mapped to Entrez IDs:", nrow(entrez_ids)))
  print(paste("Number of universe genes:", length(gene_universe)))
  print(paste("Number of universe genes mapped to Entrez IDs:", nrow(universe_entrez)))
  
  if(nrow(entrez_ids) == 0) {
    warning(paste("No genes could be mapped for", title))
    return(NULL)
  }
  
  tryCatch({
    ekegg <- enrichKEGG(gene = entrez_ids$ENTREZID,
                        universe = universe_entrez$ENTREZID,
                        organism = 'hsa',
                        keyType = "kegg",
                        pvalueCutoff = 0.05,
                        pAdjustMethod = "BH")
    
    if(nrow(ekegg@result) == 0) {
      warning(paste("No enriched KEGG pathways found for", title))
      return(NULL)
    }
    
    p <- dotplot(ekegg, showCategory = 10, title = paste("KEGG -", title)) +
      theme(axis.text.y = element_text(size = 8))
    
    print(p)
    png(paste0("KEGG_enrichment_", gsub(" ", "_", title), ".png"), width = 12, height = 8, units = "in", res = 300)
    print(p)
    dev.off()
    
    return(ekegg)
  }, error = function(e) {
    warning(paste("Error in KEGG enrichment for", title, ":", e$message))
    return(NULL)
  })
}

gene_universe <- rownames(All_samples_Merged)

# L1_vs_PBMC comparison
upregulated_genes_L1_vs_PBMC <- rownames(L1_vs_PBMC[L1_vs_PBMC$avg_log2FC > 2.5 & L1_vs_PBMC$p_val_adj < 0.05, ])
downregulated_genes_L1_vs_PBMC <- rownames(L1_vs_PBMC[L1_vs_PBMC$avg_log2FC < -2.5 & L1_vs_PBMC$p_val_adj < 0.05, ])

go_up_L1_vs_PBMC <- perform_go_enrichment(upregulated_genes_L1_vs_PBMC, gene_universe, "Upregulated Genes in L1_vs_PBMC")

go_down_L1_vs_PBMC <- perform_go_enrichment(downregulated_genes_L1_vs_PBMC, gene_universe, "Downregulated Genes in L1_vs_PBMC")

kegg_up_L1_vs_PBMC <- perform_kegg_enrichment(upregulated_genes_L1_vs_PBMC, gene_universe, "Upregulated Genes in L1_vs_PBMC")

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

[1] "Number of input genes: 1743"
[1] "Number of input genes mapped to Entrez IDs: 1537"
[1] "Number of universe genes: 27417"
[1] "Number of universe genes mapped to Entrez IDs: 19538"

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

kegg_down_L1_vs_PBMC <- perform_kegg_enrichment(downregulated_genes_L1_vs_PBMC, gene_universe, "Downregulated Genes in L1_vs_PBMC")

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

[1] "Number of input genes: 1531"
[1] "Number of input genes mapped to Entrez IDs: 1384"
[1] "Number of universe genes: 27417"
[1] "Number of universe genes mapped to Entrez IDs: 19538"

# L2_vs_PBMC comparison
upregulated_genes_L2_vs_PBMC <- rownames(L2_vs_PBMC[L2_vs_PBMC$avg_log2FC > 2.5 & L2_vs_PBMC$p_val_adj < 0.05, ])
downregulated_genes_L2_vs_PBMC <- rownames(L2_vs_PBMC[L2_vs_PBMC$avg_log2FC < -2.5 & L2_vs_PBMC$p_val_adj < 0.05, ])

go_up_L2_vs_PBMC <- perform_go_enrichment(upregulated_genes_L2_vs_PBMC, gene_universe, "Upregulated Genes in L2_vs_PBMC")

go_down_L2_vs_PBMC <- perform_go_enrichment(downregulated_genes_L2_vs_PBMC, gene_universe, "Downregulated Genes in L2_vs_PBMC")

kegg_up_L2_vs_PBMC <- perform_kegg_enrichment(upregulated_genes_L2_vs_PBMC, gene_universe, "Upregulated Genes in L2_vs_PBMC")

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

[1] "Number of input genes: 3046"
[1] "Number of input genes mapped to Entrez IDs: 2706"
[1] "Number of universe genes: 27417"
[1] "Number of universe genes mapped to Entrez IDs: 19538"

kegg_down_L2_vs_PBMC <- perform_kegg_enrichment(downregulated_genes_L2_vs_PBMC, gene_universe, "Downregulated Genes in L2_vs_PBMC")

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

[1] "Number of input genes: 1395"
[1] "Number of input genes mapped to Entrez IDs: 1252"
[1] "Number of universe genes: 27417"
[1] "Number of universe genes mapped to Entrez IDs: 19538"

# L3_vs_PBMC comparison
upregulated_genes_L3_vs_PBMC <- rownames(L3_vs_PBMC[L3_vs_PBMC$avg_log2FC > 2.5 & L3_vs_PBMC$p_val_adj < 0.05, ])
downregulated_genes_L3_vs_PBMC <- rownames(L3_vs_PBMC[L3_vs_PBMC$avg_log2FC < -2.5 & L3_vs_PBMC$p_val_adj < 0.05, ])

go_up_L3_vs_PBMC <- perform_go_enrichment(upregulated_genes_L3_vs_PBMC, gene_universe, "Upregulated Genes in L3_vs_PBMC")

go_down_L3_vs_PBMC <- perform_go_enrichment(downregulated_genes_L3_vs_PBMC, gene_universe, "Downregulated Genes in L3_vs_PBMC")

kegg_up_L3_vs_PBMC <- perform_kegg_enrichment(upregulated_genes_L3_vs_PBMC, gene_universe, "Upregulated Genes in L3_vs_PBMC")

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

[1] "Number of input genes: 2455"
[1] "Number of input genes mapped to Entrez IDs: 2110"
[1] "Number of universe genes: 27417"
[1] "Number of universe genes mapped to Entrez IDs: 19538"

kegg_down_L3_vs_PBMC <- perform_kegg_enrichment(downregulated_genes_L3_vs_PBMC, gene_universe, "Downregulated Genes in L3_vs_PBMC")

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

[1] "Number of input genes: 1629"
[1] "Number of input genes mapped to Entrez IDs: 1438"
[1] "Number of universe genes: 27417"
[1] "Number of universe genes mapped to Entrez IDs: 19538"

# L4_vs_PBMC comparison
upregulated_genes_L4_vs_PBMC <- rownames(L4_vs_PBMC[L4_vs_PBMC$avg_log2FC > 2.5 & L4_vs_PBMC$p_val_adj < 0.05, ])
downregulated_genes_L4_vs_PBMC <- rownames(L4_vs_PBMC[L4_vs_PBMC$avg_log2FC < -2.5 & L4_vs_PBMC$p_val_adj < 0.05, ])

go_up_L4_vs_PBMC <- perform_go_enrichment(upregulated_genes_L4_vs_PBMC, gene_universe, "Upregulated Genes in L4_vs_PBMC")

go_down_L4_vs_PBMC <- perform_go_enrichment(downregulated_genes_L4_vs_PBMC, gene_universe, "Downregulated Genes in L4_vs_PBMC")

kegg_up_L4_vs_PBMC <- perform_kegg_enrichment(upregulated_genes_L4_vs_PBMC, gene_universe, "Upregulated Genes in L4_vs_PBMC")

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

[1] "Number of input genes: 3862"
[1] "Number of input genes mapped to Entrez IDs: 3328"
[1] "Number of universe genes: 27417"
[1] "Number of universe genes mapped to Entrez IDs: 19538"

kegg_down_L4_vs_PBMC <- perform_kegg_enrichment(downregulated_genes_L4_vs_PBMC, gene_universe, "Downregulated Genes in L4_vs_PBMC")

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

[1] "Number of input genes: 1155"
[1] "Number of input genes mapped to Entrez IDs: 1017"
[1] "Number of universe genes: 27417"
[1] "Number of universe genes mapped to Entrez IDs: 19538"

# L5_vs_PBMC comparison
upregulated_genes_L5_vs_PBMC <- rownames(L5_vs_PBMC[L5_vs_PBMC$avg_log2FC > 2.5 & L5_vs_PBMC$p_val_adj < 0.05, ])
downregulated_genes_L5_vs_PBMC <- rownames(L5_vs_PBMC[L5_vs_PBMC$avg_log2FC < -2.5 & L5_vs_PBMC$p_val_adj < 0.05, ])

go_up_L5_vs_PBMC <- perform_go_enrichment(upregulated_genes_L5_vs_PBMC, gene_universe, "Upregulated Genes in L5_vs_PBMC")

go_down_L5_vs_PBMC <- perform_go_enrichment(downregulated_genes_L5_vs_PBMC, gene_universe, "Downregulated Genes in L5_vs_PBMC")

kegg_up_L5_vs_PBMC <- perform_kegg_enrichment(upregulated_genes_L5_vs_PBMC, gene_universe, "Upregulated Genes in L5_vs_PBMC")

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

[1] "Number of input genes: 2974"
[1] "Number of input genes mapped to Entrez IDs: 2587"
[1] "Number of universe genes: 27417"
[1] "Number of universe genes mapped to Entrez IDs: 19538"

kegg_down_L5_vs_PBMC <- perform_kegg_enrichment(downregulated_genes_L5_vs_PBMC, gene_universe, "Downregulated Genes in L5_vs_PBMC")

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

[1] "Number of input genes: 1521"
[1] "Number of input genes mapped to Entrez IDs: 1348"
[1] "Number of universe genes: 27417"
[1] "Number of universe genes mapped to Entrez IDs: 19538"

# L6_vs_PBMC comparison
upregulated_genes_L6_vs_PBMC <- rownames(L6_vs_PBMC[L6_vs_PBMC$avg_log2FC > 2.5 & L6_vs_PBMC$p_val_adj < 0.05, ])
downregulated_genes_L6_vs_PBMC <- rownames(L6_vs_PBMC[L6_vs_PBMC$avg_log2FC < -2.5 & L6_vs_PBMC$p_val_adj < 0.05, ])

go_up_L6_vs_PBMC <- perform_go_enrichment(upregulated_genes_L6_vs_PBMC, gene_universe, "Upregulated Genes in L6_vs_PBMC")

go_down_L6_vs_PBMC <- perform_go_enrichment(downregulated_genes_L6_vs_PBMC, gene_universe, "Downregulated Genes in L6_vs_PBMC")

kegg_up_L6_vs_PBMC <- perform_kegg_enrichment(upregulated_genes_L6_vs_PBMC, gene_universe, "Upregulated Genes in L6_vs_PBMC")

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

[1] "Number of input genes: 2271"
[1] "Number of input genes mapped to Entrez IDs: 1895"
[1] "Number of universe genes: 27417"
[1] "Number of universe genes mapped to Entrez IDs: 19538"

kegg_down_L6_vs_PBMC <- perform_kegg_enrichment(downregulated_genes_L6_vs_PBMC, gene_universe, "Downregulated Genes in L6_vs_PBMC")

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

[1] "Number of input genes: 1429"
[1] "Number of input genes mapped to Entrez IDs: 1289"
[1] "Number of universe genes: 27417"
[1] "Number of universe genes mapped to Entrez IDs: 19538"

# L7_vs_PBMC comparison
upregulated_genes_L7_vs_PBMC <- rownames(L7_vs_PBMC[L7_vs_PBMC$avg_log2FC > 2.5 & L7_vs_PBMC$p_val_adj < 0.05, ])
downregulated_genes_L7_vs_PBMC <- rownames(L7_vs_PBMC[L7_vs_PBMC$avg_log2FC < -2.5 & L7_vs_PBMC$p_val_adj < 0.05, ])

go_up_L7_vs_PBMC <- perform_go_enrichment(upregulated_genes_L7_vs_PBMC, gene_universe, "Upregulated Genes in L7_vs_PBMC")

go_down_L7_vs_PBMC <- perform_go_enrichment(downregulated_genes_L7_vs_PBMC, gene_universe, "Downregulated Genes in L7_vs_PBMC")

kegg_up_L7_vs_PBMC <- perform_kegg_enrichment(upregulated_genes_L7_vs_PBMC, gene_universe, "Upregulated Genes in L7_vs_PBMC")

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

[1] "Number of input genes: 2605"
[1] "Number of input genes mapped to Entrez IDs: 2266"
[1] "Number of universe genes: 27417"
[1] "Number of universe genes mapped to Entrez IDs: 19538"

kegg_down_L7_vs_PBMC <- perform_kegg_enrichment(downregulated_genes_L7_vs_PBMC, gene_universe, "Downregulated Genes in L7_vs_PBMC")

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

[1] "Number of input genes: 1510"
[1] "Number of input genes mapped to Entrez IDs: 1348"
[1] "Number of universe genes: 27417"
[1] "Number of universe genes mapped to Entrez IDs: 19538"