Sézary Syndrome SCpubr Visualization

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)

3. Dim plots

sample <- All_samples_Merged

# Seurat's DimPlot.
p1 <- Seurat::DimPlot(sample, group.by = "Harmony_snn_res.0.9", reduction = "umap.harmony")

# SCpubr's DimPlot.
p2 <- SCpubr::do_DimPlot(sample = sample, group.by = "Harmony_snn_res.0.9", reduction = "umap.harmony")

p <- p1 | p2
p

SCpubr::do_DimPlot(sample = sample, group.by = "Harmony_snn_res.0.9", reduction = "umap.harmony", label = T, label.box = T)

SCpubr::do_DimPlot(sample = sample, group.by = "Harmony_snn_res.0.9", reduction = "umap.harmony", label = T, label.box = T, plot.axes = T)

SCpubr::do_DimPlot(sample = sample, group.by = "Harmony_snn_res.0.9", reduction = "umap.harmony", label = T, label.box = T, plot.axes = T, legend.position = "none")

SCpubr::do_DimPlot(sample = sample,
                   idents.keep = c("3", "8", "10", "18"), group.by = "Harmony_snn_res.0.9", reduction = "umap.harmony", label = T, label.box = T, plot.axes = T, legend.position = "none")

SCpubr::do_DimPlot(sample = sample,
                   idents.keep = c("1", "2", "13"), group.by = "Harmony_snn_res.0.9", reduction = "umap.harmony", label = T, label.box = T, plot.axes = T, legend.position = "none")

SCpubr::do_DimPlot(sample = sample,
                   idents.keep = c("4", "7", "9", "6", "16", "19"), group.by = "Harmony_snn_res.0.9", reduction = "umap.harmony", label = T, label.box = T, plot.axes = T, legend.position = "none")

NA
NA
NA
NA

4. Feature Plots

genes <- list("Naive CD4+ T" = c("IL7R", "CCR7"),
              "CD14+ Mono" = c("CD14", "LYZ"),
              "Memory CD4+" = c("S100A4"),
              "B" = c("MS4A1"),
              "CD8+ T" = c("CD8A"),
              "FCGR3A+ Mono" = c("FCGR3A", "MS4A7"),
              "NK" = c("GNLY", "NKG7"),
              "DC" = c("FCER1A", "CST3"),
              "Platelet" = c("PPBP"))

p1 <- SCpubr::do_DotPlot(sample = sample, 
                        features = genes)
p1

p2 <- SCpubr::do_DotPlot(sample = sample, 
                        features = genes, group.by = "cell_line")
p2

SCpubr::do_BarPlot(sample = sample, 
                         group.by = "Harmony_snn_res.0.9", 
                         legend.position = "none",
                         plot.title = "Number of cells per cluster", 
                         flip = TRUE)

SCpubr::do_ChordDiagramPlot(sample = sample,
                            from = "orig.ident",
                            to = "cell_line",
                            link.arr.type = "triangle")

NA
NA
NA

5. Group-wise DE analysis plots

# Set the identities correctly.
Seurat::Idents(sample) <- sample$Harmony_snn_res.0.9

# Compute DE genes and transform to a tibble.
de_genes <- tibble::tibble(Seurat::FindAllMarkers(object = sample, min.pct = 0.1, logfc.threshold = 0.1))

Calculating cluster 0
Calculating cluster 1
Calculating cluster 2
Calculating cluster 3
Calculating cluster 4
Calculating cluster 5
Calculating cluster 6
Calculating cluster 7
Calculating cluster 8
Calculating cluster 9
Calculating cluster 10
Calculating cluster 11
Calculating cluster 12
Calculating cluster 13
Calculating cluster 14
Calculating cluster 15
Calculating cluster 16
Calculating cluster 17
Calculating cluster 18
Calculating cluster 19
Calculating cluster 20
Calculating cluster 21
Calculating cluster 22
Calculating cluster 23
Calculating cluster 24

Group-wise DE analysis plots

p <- SCpubr::do_GroupwiseDEPlot(sample = sample,
                                de_genes = de_genes,
                                min.cutoff = 1)

p

# Increase the number of top DE genes by cluster.
p <- SCpubr::do_GroupwiseDEPlot(sample = sample,
                                de_genes = de_genes,
                                top_genes = 10)

p

# Define list of genes.
genes <- c("KIR3DL2", 
           "TOX", 
           "CD7", 
           "CD26", 
           "DDP4", 
           "TWIST1",
           "PLS3",
           "CD70",
           "CCR4", 
           "CCR7")

# Default parameters.
 SCpubr::do_ExpressionHeatmap(sample = sample,
                                  features = genes)

Avis : 
! The following features are not found in the row names of the specified assay (default assay if not):
CD26, DDP4

# Define the gene list
genes_to_plot <- c("SPINK2", "PRSS57", "CYTL1", "EGFL7", "GATA2", "CD34", "SMIM24", "AVP", "MYB", "LAPTM4B")

SCpubr::do_ExpressionHeatmap(sample = sample,
                                  features = genes, group.by = "Harmony_snn_res.0.9")

Avis : 
! The following features are not found in the row names of the specified assay (default assay if not):
CD26, DDP4

Group-wise DE analysis plots

p <- SCpubr::do_GroupwiseDEPlot(sample = sample,
                                de_genes = de_genes,
                                top_genes = 5)

p

# Define list of genes.
genes <- c("KIR3DL2", 
           "TOX", 
           "CD7", 
           "DPP4", 
           "TWIST1",
           "PLS3",
           "CD70",
           "CCR4", 
           "CCR7")

# Default parameters.
 SCpubr::do_ExpressionHeatmap(sample = sample,
                                  features = genes)

Group-wise DE analysis plots

library(ggplot2)
library(RColorBrewer)

# Step 1: Filter for HSPC cells
hspc_data <- All_samples_Merged@meta.data %>%
  filter(predicted.celltype.l2 == "HSPC")  # Filter for HSPC cells

# Step 2: Count HSPC cells per cell line
hspc_counts <- as.data.frame(table(hspc_data$cell_line))
colnames(hspc_counts) <- c("cell_line", "nUMI")  # Rename columns for clarity

# Step 3: Generate color palette
# Adjust the number (10) to the number of unique cell lines
cell_line_colors <- brewer.pal(min(10, nrow(hspc_counts)), "Set3")

# Step 4: Create bar plot
hspc_plot <- ggplot(hspc_counts, aes(x = cell_line, y = nUMI, fill = cell_line)) +
  geom_col() +
  theme_classic() +
  geom_text(aes(label = nUMI),
            position = position_dodge(width = 0.9),
            vjust = -0.25) +
  scale_fill_manual(values = cell_line_colors) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        plot.title = element_text(hjust = 0.5)) +  # Center the title
  ggtitle("HSPC Cells per Cell Line") +
  xlab("Cell Lines") +
  ylab("Number of HSPC Cells")

# Step 5: Print the plot
print(hspc_plot)

# Filter for HSPC cells
hspc_data <- All_samples_Merged@meta.data %>%
  filter(predicted.celltype.l2 == "HSPC")

# Count HSPC cells per cell line
hspc_counts <- as.data.frame(table(hspc_data$cell_line))
colnames(hspc_counts) <- c("cell_line", "nUMI")  # Rename columns for clarity

# Exclude "PBMC" and "PBMC 10X" cell lines
hspc_counts <- hspc_counts %>%
  filter(!cell_line %in% c("PBMC", "PBMC_10x"))

# Generate a color palette for the remaining cell lines
cell_line_colors <- brewer.pal(min(10, nrow(hspc_counts)), "Set3")

# Create the bar plot
hspc_plot <- ggplot(hspc_counts, aes(x = cell_line, y = nUMI, fill = cell_line)) +
  geom_col() +
  theme_classic() +
  geom_text(aes(label = nUMI),
            position = position_dodge(width = 0.9),
            vjust = -0.25) +
  scale_fill_manual(values = cell_line_colors) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        plot.title = element_text(hjust = 0.5)) +
  ggtitle("HSPC Cells per Cell Line (Excluding PBMC)") +
  xlab("Cell Lines") +
  ylab("Number of HSPC Cells")

# Print the plot
print(hspc_plot)

NA
NA
NA
NA
NA
NA

Azimuth Barplot

# Count cells for each cell type in predicted.celltype.l2
celltype_counts <- as.data.frame(table(All_samples_Merged$predicted.celltype.l2))
colnames(celltype_counts) <- c("cell_type", "nUMI")  # Rename columns appropriately

# Create an extended color palette
celltype_colors <- colorRampPalette(brewer.pal(12, "Set3"))(nrow(celltype_counts))

# Create bar plot
celltype_plot <- ggplot(celltype_counts, aes(x = cell_type, y = nUMI, fill = cell_type)) +
  geom_col() +
  theme_classic() +
  geom_text(aes(label = nUMI), 
            position = position_dodge(width = 0.9), 
            vjust = -0.25) +
  scale_fill_manual(values = celltype_colors) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        plot.title = element_text(hjust = 0.5)) +
  ggtitle("Cell Counts by Predicted Cell Type") +
  xlab("Cell Types") +
  ylab("Number of Cells (nUMI)")

# Print the plot
print(celltype_plot)

# List of cell types to exclude
exclude_cell_types <- c("pDC", "Plasmablast", "cDC1", "Platelet", "ILC", "ASDC", "cDC2")

# Filter out the unwanted cell types
celltype_counts_filtered <- celltype_counts[!celltype_counts$cell_type %in% exclude_cell_types, ]

# Generate the color palette
celltype_colors <- colorRampPalette(brewer.pal(12, "Set3"))(nrow(celltype_counts_filtered))

# Create the bar plot
celltype_plot <- ggplot(celltype_counts_filtered, aes(x = cell_type, y = nUMI, fill = cell_type)) +
  geom_col() +
  theme_classic() +
  geom_text(aes(label = nUMI), 
            position = position_dodge(width = 0.9), 
            vjust = -0.25) +
  scale_fill_manual(values = celltype_colors) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        plot.title = element_text(hjust = 0.5)) +
  ggtitle("Cell Counts by Predicted Cell Type") +
  xlab("Cell Types") +
  ylab("Number of Cells (nUMI)")

# Print the plot
print(celltype_plot)

NA
NA
NA