Libraries & functions

First, I’m loading the required libraries & functions

suppressPackageStartupMessages(library(Seurat, lib.loc = "/software/Seuratv4/lib/")) # Seurat v4.4.0, library for single-cell analysis
suppressPackageStartupMessages(library(SeuratObject, lib.loc = "/software/Seuratv4/lib/")) # For compatibility purpose
suppressPackageStartupMessages(library(data.table)) # For writing DE gene file
suppressPackageStartupMessages(library(readxl)) # For writing DE gene file
suppressPackageStartupMessages(library(loomR)) # For building a file for ASAP
suppressPackageStartupMessages(library(Matrix)) # For building a file for ASAP
suppressPackageStartupMessages(library(ggplot2)) # For plotting
suppressPackageStartupMessages(library(ggpubr)) # For grid plotting
suppressPackageStartupMessages(library(plotly)) # For interactive plots
suppressPackageStartupMessages(library(crayon)) # Just for bolding the console output :D

cat(bold("Seurat"), "version", as.character(packageVersion("Seurat")), "\n")
## Seurat version 4.4.0
cat(bold("SeuratObject"), "version", as.character(packageVersion("SeuratObject")), "\n")
## SeuratObject version 4.1.4
cat(bold("data.table"), "version", as.character(packageVersion("data.table")), "\n")
## data.table version 1.16.4
cat(bold("readxl"), "version", as.character(packageVersion("readxl")), "\n")
## readxl version 1.4.3
cat(bold("loomR"), "version", as.character(packageVersion("loomR")), "\n")
## loomR version 0.2.0
cat(bold("Matrix"), "version", as.character(packageVersion("Matrix")), "\n")
## Matrix version 1.6.5
cat(bold("ggplot2"), "version", as.character(packageVersion("ggplot2")), "\n")
## ggplot2 version 3.5.1
cat(bold("ggpubr"), "version", as.character(packageVersion("ggpubr")), "\n")
## ggpubr version 0.6.0
cat(bold("plotly"), "version", as.character(packageVersion("plotly")), "\n")
## plotly version 4.10.4
# Random seed
set.seed(42)

Parameters

Paper: Unravelling cell type-specific responses to Parkinson’s Disease at single cell resolution Martirosyan et al. Mol Neurodegeneration 19, 7 (2024) https://molecularneurodegeneration.biomedcentral.com/articles/10.1186/s13024-023-00699-0

Data repo: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE243639

Loom file: NA

## Parameters
input_matrix <- "./data/Martirosyan_2024/GSE243639_Filtered_count_table.csv.gz"
input_metadata <- "./data/Martirosyan_2024/GSE243639_UMAP_coordinates.xlsx"
input_Oxphos_135_AUCell_values_path <- "./data/Martirosyan_135_Oxphos_AUCell_auc.tsv"
input_UPR_93_AUCell_values_path <- "./data/Martirosyan_93_UPR_AUCell_auc.tsv"
output_figure_prefix <- "./figures/Figure_5A/"

I. Loading dataset

I.1. Reading count matrix from csv provided by authors on GEO

data.matrix <- fread(input_matrix, header = T, stringsAsFactors = F, showProgress = F, sep = ",", data.table = F, check.names = F)
rownames(data.matrix) <- data.matrix$V1
data.matrix <- data.matrix[,-1]

I.2. Build Seurat object

# Build a Seurat object with this data
data.seurat <- CreateSeuratObject(data.matrix, project = "Martirosyan")

# Restore hyphens in cell barcodes
data.seurat <- RenameCells(data.seurat, new.names = gsub(pattern = ".1$", replacement = "-1", x = colnames(data.seurat)))

# Change identity
data.seurat$orig.ident.parsed <- data.seurat$orig.ident
data.seurat$orig.ident <- "Martirosyan"
data.seurat <- SetIdent(data.seurat, value = "orig.ident")

# Add MT content
data.seurat[["percent.mt"]] <- PercentageFeatureSet(data.seurat, pattern = "^MT-")

# Run default Normalization
data.seurat <- NormalizeData(data.seurat, normalization.method = "LogNormalize", scale.factor = 10000)
# Visualize QC metrics as a violin plot
VlnPlot(data.seurat, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

I.3. Reading global metadata from authors

data.metadata <- as.data.frame(read_xlsx(path = input_metadata, sheet = 1))
rownames(data.metadata) <- data.metadata$CELL_ID
data.metadata$UMAP_1 <- as.numeric(data.metadata$UMAP_1)
data.metadata$UMAP_2 <- as.numeric(data.metadata$UMAP_2)
data.metadata

I.4 Reading OXPHOS AUC values from AUCell (using the 133 OXPHOS genes from KEGG Oxphos pathway https://www.genome.jp/entry/mmu00190)

data.auc.oxphos <- fread(input_Oxphos_135_AUCell_values_path, data.table = F, check.names = F)
rownames(data.auc.oxphos) <- data.auc.oxphos$Cell
data.auc.oxphos

I.5 Reading UPR AUC values from AUCell (using the 93 UPR genes from REACTOME

data.auc.upr <- fread(input_UPR_93_AUCell_values_path, data.table = F, check.names = F)
rownames(data.auc.upr) <- data.auc.upr$Cell
data.auc.upr

I.6 Add to Seurat Object

# Add metadata
data.seurat$AUC_OXPHOS_135_KEGG <- data.auc.oxphos[colnames(data.seurat), "AUC"]
data.seurat$AUC_UPR_93_REACTOME <- data.auc.upr[colnames(data.seurat), "AUC"]
data.seurat$CLUSTER <- data.metadata[colnames(data.seurat), "CLUSTER"]
data.seurat$IDENT <- data.metadata[colnames(data.seurat), "IDENT"]

# Add UMAP
data.seurat[["UMAP"]] <- CreateDimReducObject(embeddings = as.matrix(data.metadata[,c("UMAP_1", "UMAP_2")]), key = "UMAP_", assay = "RNA")

II Plot interesting stuff

II.1 Plot author’s UMAP

DimPlot(data.seurat, reduction = "UMAP", group.by = "IDENT")

DimPlot(data.seurat, reduction = "UMAP", group.by = "CLUSTER")

II.2.a Plot Distribution OXPHOS

hist(data.seurat$AUC_OXPHOS_135_KEGG, breaks = 100)

II.2.b Plot Boxplot OXPHOS

boxplot_data <- data.seurat@meta.data

# Filter out <5 cells per annotation
boxplot_data <- boxplot_data[boxplot_data$IDENT %in% names(which(table(boxplot_data$IDENT) >= 5)),]

# Calculating median values of "AUC_OXPHOS_135_KEGG" for each group
median_values <- tapply(boxplot_data$AUC_OXPHOS_135_KEGG, boxplot_data$IDENT, median)

# Sorting the metadata based on median values
boxplot_data$IDENT <- factor(boxplot_data$IDENT, levels = names(sort(median_values, decreasing = F)))

# Plot
p <- ggplot(boxplot_data, aes(x = IDENT, y = AUC_OXPHOS_135_KEGG, fill = IDENT)) +
  geom_boxplot() +
  labs(x = NULL, y = "AUC_OXPHOS_135_KEGG") +
  theme_minimal() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  NoLegend()

p

#ggsave(p, filename = paste0(output_figure_prefix, "Martirosyan_135_Oxphos_AUCell_ordered_boxplot_IDENT.pdf"), width = 5, height = 5, bg = "white")
#ggsave(p, filename = paste0(output_figure_prefix, "Martirosyan_135_Oxphos_AUCell_ordered_boxplot_IDENT.png"), width = 5, height = 5, dpi = 1000, bg = 'white')

II.2.c Plot UMAP OXPHOS

FeaturePlot(data.seurat, reduction = "UMAP", features = "AUC_OXPHOS_135_KEGG")

II.3.a Plot Distribution UPR

hist(data.seurat$AUC_UPR_93_REACTOME, breaks = 100)

II.3.b Plot Boxplot UPR

boxplot_data <- data.seurat@meta.data

# Filter out <5 cells per annotation
boxplot_data <- boxplot_data[boxplot_data$IDENT %in% names(which(table(boxplot_data$IDENT) >= 5)),]

# Calculating median values of "AUC_UPR_93_REACTOME" for each group
median_values <- tapply(boxplot_data$AUC_UPR_93_REACTOME, boxplot_data$IDENT, median)

# Sorting the metadata based on median values
boxplot_data$IDENT <- factor(boxplot_data$IDENT, levels = names(sort(median_values, decreasing = F)))

# Plot
p <- ggplot(boxplot_data, aes(x = IDENT, y = AUC_UPR_93_REACTOME, fill = IDENT)) +
  geom_boxplot() +
  labs(x = NULL, y = "AUC_UPR_93_REACTOME") +
  theme_minimal() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  NoLegend()

p

#ggsave(p, filename = paste0(output_figure_prefix, "Martirosyan_93_UPR_AUCell_ordered_boxplot_IDENT.pdf"), width = 5, height = 5, bg = "white")
#ggsave(p, filename = paste0(output_figure_prefix, "Martirosyan_93_UPR_AUCell_ordered_boxplot_IDENT.png"), width = 5, height = 5, dpi = 1000, bg = 'white')

II.3.c Plot UMAP UPR

FeaturePlot(data.seurat, reduction = "UMAP", features = "AUC_UPR_93_REACTOME")

II.4.a Plot Boxplot HSP90AA1

boxplot_data <- data.seurat@meta.data
boxplot_data$HSP90AA1 <- data.seurat@assays$RNA@data["HSP90AA1",]

# Filter out <5 cells per annotation
boxplot_data <- boxplot_data[boxplot_data$IDENT %in% names(which(table(boxplot_data$IDENT) >= 5)),]

# Calculating median values of "HSP90AA1" for each group
median_values <- tapply(boxplot_data$HSP90AA1, boxplot_data$IDENT, median)

# Sorting the metadata based on median values
boxplot_data$IDENT <- factor(boxplot_data$IDENT, levels = names(sort(median_values, decreasing = F)))

# Plot
p <- ggplot(boxplot_data, aes(x = IDENT, y = HSP90AA1, fill = IDENT)) +
  geom_boxplot() +
  labs(x = NULL, y = "HSP90AA1 Normalized expression") +
  theme_minimal() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  NoLegend()

p

#ggsave(p, filename = paste0(output_figure_prefix, "Martirosyan_HSP90AA1_ordered_boxplot_IDENT.pdf"), width = 5, height = 5, bg = "white")
#ggsave(p, filename = paste0(output_figure_prefix, "Martirosyan_HSP90AA1_ordered_boxplot_IDENT.png"), width = 5, height = 5, dpi = 1000, bg = 'white')

II.4.b Plot UMAP HSP90AA1

FeaturePlot(data.seurat, reduction = "UMAP", features = "HSP90AA1")

II.5.a Plot Boxplot HSPA8

boxplot_data <- data.seurat@meta.data
boxplot_data$HSPA8 <- data.seurat@assays$RNA@data["HSPA8",]

# Filter out <5 cells per annotation
boxplot_data <- boxplot_data[boxplot_data$IDENT %in% names(which(table(boxplot_data$IDENT) >= 5)),]

# Calculating median values of "HSPA8" for each group
median_values <- tapply(boxplot_data$HSPA8, boxplot_data$IDENT, median)

# Sorting the metadata based on median values
boxplot_data$IDENT <- factor(boxplot_data$IDENT, levels = names(sort(median_values, decreasing = F)))

# Plot
p <- ggplot(boxplot_data, aes(x = IDENT, y = HSPA8, fill = IDENT)) +
  geom_boxplot() +
  labs(x = NULL, y = "HSPA8 Normalized expression") +
  theme_minimal() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  NoLegend()

p

#ggsave(p, filename = paste0(output_figure_prefix, "Martirosyan_HSPA8_ordered_boxplot_IDENT.pdf"), width = 5, height = 5, bg = "white")
#ggsave(p, filename = paste0(output_figure_prefix, "Martirosyan_HSPA8_ordered_boxplot_IDENT.png"), width = 5, height = 5, dpi = 1000, bg = 'white')

II.5.b Plot UMAP HSPA8

FeaturePlot(data.seurat, reduction = "UMAP", features = "HSPA8")

II.6.a Plot Boxplot ATF4

boxplot_data <- data.seurat@meta.data
boxplot_data$ATF4 <- data.seurat@assays$RNA@data["ATF4",]

# Filter out <5 cells per annotation
boxplot_data <- boxplot_data[boxplot_data$IDENT %in% names(which(table(boxplot_data$IDENT) >= 5)),]

# Calculating median values of "ATF4" for each group
median_values <- tapply(boxplot_data$ATF4, boxplot_data$IDENT, median)

# Sorting the metadata based on median values
boxplot_data$IDENT <- factor(boxplot_data$IDENT, levels = names(sort(median_values, decreasing = F)))

# Plot
p <- ggplot(boxplot_data, aes(x = IDENT, y = ATF4, fill = IDENT)) +
  geom_boxplot() +
  labs(x = NULL, y = "ATF4 Normalized expression") +
  theme_minimal() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  NoLegend()

p

#ggsave(p, filename = paste0(output_figure_prefix, "Martirosyan_ATF4_ordered_boxplot_IDENT.pdf"), width = 5, height = 5, bg = "white")
#ggsave(p, filename = paste0(output_figure_prefix, "Martirosyan_ATF4_ordered_boxplot_IDENT.png"), width = 5, height = 5, dpi = 1000, bg = 'white')

II.6.b Plot UMAP ATF4

FeaturePlot(data.seurat, reduction = "UMAP", features = "ATF4")

III Subsetting to only neurons

III.1 Read metadata from authors

data.neurons <- as.data.frame(read_xlsx(path = input_metadata, sheet = "Neurons"))
rownames(data.neurons) <- data.neurons$CELL_ID
data.neurons$UMAP_1 <- as.numeric(data.neurons$UMAP_1)
data.neurons$UMAP_2 <- as.numeric(data.neurons$UMAP_2)
data.neurons

Authors have subset their dataset and redo a UMAP and clustering on the subset dataset.

III.1 Subset the Seurat object

# Subset to neuron cells from metadata
neuron.cells <- rownames(data.neurons)[rownames(data.neurons) %in% colnames(data.seurat)]
data.seurat <- subset(data.seurat, cells = neuron.cells)

# Add/replace new metadata
data.seurat$CLUSTER <- data.neurons[colnames(data.seurat), "CLUSTER"]
data.seurat$IDENT <- data.neurons[colnames(data.seurat), "IDENT"]

# Add UMAP
data.seurat[["UMAP"]] <- CreateDimReducObject(embeddings = as.matrix(data.neurons[,c("UMAP_1", "UMAP_2")]), key = "UMAP_", assay = "RNA")

IV Plot interesting stuff

IV.1 Plot author’s UMAP

DimPlot(data.seurat, reduction = "UMAP", group.by = "IDENT")

DimPlot(data.seurat, reduction = "UMAP", group.by = "CLUSTER")

IV.2.a Plot Distribution OXPHOS

hist(data.seurat$AUC_OXPHOS_135_KEGG, breaks = 100)

IV.2.b Plot Boxplot OXPHOS

boxplot_data <- data.seurat@meta.data

# Filter out <5 cells per annotation
boxplot_data <- boxplot_data[boxplot_data$IDENT %in% names(which(table(boxplot_data$IDENT) >= 5)),]

# Calculating median values of "AUC_OXPHOS_135_KEGG" for each group
median_values <- tapply(boxplot_data$AUC_OXPHOS_135_KEGG, boxplot_data$IDENT, median)

# Sorting the metadata based on median values
boxplot_data$IDENT <- factor(boxplot_data$IDENT, levels = names(sort(median_values, decreasing = F)))

# Plot
p <- ggplot(boxplot_data, aes(x = IDENT, y = AUC_OXPHOS_135_KEGG, fill = IDENT)) +
  geom_boxplot() +
  labs(x = NULL, y = "AUC_OXPHOS_135_KEGG") +
  theme_minimal() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  NoLegend()

p

#ggsave(p, filename = paste0(output_figure_prefix, "Martirosyan_135_Oxphos_AUCell_ordered_boxplot_IDENT_Neurons.pdf"), width = 5, height = 5, bg = "white")
#ggsave(p, filename = paste0(output_figure_prefix, "Martirosyan_135_Oxphos_AUCell_ordered_boxplot_IDENT_Neurons.png"), width = 5, height = 5, dpi = 1000, bg = 'white')

IV.2.c Plot UMAP OXPHOS

FeaturePlot(data.seurat, reduction = "UMAP", features = "AUC_OXPHOS_135_KEGG")

IV.3.a Plot Distribution UPR

hist(data.seurat$AUC_UPR_93_REACTOME, breaks = 100)

IV.3.b Plot Boxplot UPR

boxplot_data <- data.seurat@meta.data

# Filter out <5 cells per annotation
boxplot_data <- boxplot_data[boxplot_data$IDENT %in% names(which(table(boxplot_data$IDENT) >= 5)),]

# Calculating median values of "AUC_UPR_93_REACTOME" for each group
median_values <- tapply(boxplot_data$AUC_UPR_93_REACTOME, boxplot_data$IDENT, median)

# Sorting the metadata based on median values
boxplot_data$IDENT <- factor(boxplot_data$IDENT, levels = names(sort(median_values, decreasing = F)))

# Plot
p <- ggplot(boxplot_data, aes(x = IDENT, y = AUC_UPR_93_REACTOME, fill = IDENT)) +
  geom_boxplot() +
  labs(x = NULL, y = "AUC_UPR_93_REACTOME") +
  theme_minimal() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  NoLegend()

p

#ggsave(p, filename = paste0(output_figure_prefix, "Martirosyan_93_UPR_AUCell_ordered_boxplot_IDENT_Neurons.pdf"), width = 5, height = 5, bg = "white")
#ggsave(p, filename = paste0(output_figure_prefix, "Martirosyan_93_UPR_AUCell_ordered_boxplot_IDENT_Neurons.png"), width = 5, height = 5, dpi = 1000, bg = 'white')

IV.3.c Plot UMAP UPR

FeaturePlot(data.seurat, reduction = "UMAP", features = "AUC_UPR_93_REACTOME")

IV.4.a Plot Boxplot HSP90AA1

boxplot_data <- data.seurat@meta.data
boxplot_data$HSP90AA1 <- data.seurat@assays$RNA@data["HSP90AA1",]

# Filter out <5 cells per annotation
boxplot_data <- boxplot_data[boxplot_data$IDENT %in% names(which(table(boxplot_data$IDENT) >= 5)),]

# Calculating median values of "HSP90AA1" for each group
median_values <- tapply(boxplot_data$HSP90AA1, boxplot_data$IDENT, median)

# Sorting the metadata based on median values
boxplot_data$IDENT <- factor(boxplot_data$IDENT, levels = names(sort(median_values, decreasing = F)))

# Plot
p <- ggplot(boxplot_data, aes(x = IDENT, y = HSP90AA1, fill = IDENT)) +
  geom_boxplot() +
  labs(x = NULL, y = "HSP90AA1 Normalized expression") +
  theme_minimal() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  NoLegend()

p

#ggsave(p, filename = paste0(output_figure_prefix, "Martirosyan_HSP90AA1_ordered_boxplot_IDENT_Neurons.pdf"), width = 5, height = 5, bg = "white")
#ggsave(p, filename = paste0(output_figure_prefix, "Martirosyan_HSP90AA1_ordered_boxplot_IDENT_Neurons.png"), width = 5, height = 5, dpi = 1000, bg = 'white')

IV.4.b Plot UMAP HSP90AA1

FeaturePlot(data.seurat, reduction = "UMAP", features = "HSP90AA1")

IV.5.a Plot Boxplot HSPA8

boxplot_data <- data.seurat@meta.data
boxplot_data$HSPA8 <- data.seurat@assays$RNA@data["HSPA8",]

# Filter out <5 cells per annotation
boxplot_data <- boxplot_data[boxplot_data$IDENT %in% names(which(table(boxplot_data$IDENT) >= 5)),]

# Calculating median values of "HSPA8" for each group
median_values <- tapply(boxplot_data$HSPA8, boxplot_data$IDENT, median)

# Sorting the metadata based on median values
boxplot_data$IDENT <- factor(boxplot_data$IDENT, levels = names(sort(median_values, decreasing = F)))

# Plot
p <- ggplot(boxplot_data, aes(x = IDENT, y = HSPA8, fill = IDENT)) +
  geom_boxplot() +
  labs(x = NULL, y = "HSPA8 Normalized expression") +
  theme_minimal() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  NoLegend()

p

#ggsave(p, filename = paste0(output_figure_prefix, "Martirosyan_HSPA8_ordered_boxplot_IDENT_Neurons.pdf"), width = 5, height = 5, bg = "white")
#ggsave(p, filename = paste0(output_figure_prefix, "Martirosyan_HSPA8_ordered_boxplot_IDENT_Neurons.png"), width = 5, height = 5, dpi = 1000, bg = 'white')

IV.5.b Plot UMAP HSPA8

FeaturePlot(data.seurat, reduction = "UMAP", features = "HSPA8")

IV.6.a Plot Boxplot ATF4

boxplot_data <- data.seurat@meta.data
boxplot_data$ATF4 <- data.seurat@assays$RNA@data["ATF4",]

# Filter out <5 cells per annotation
boxplot_data <- boxplot_data[boxplot_data$IDENT %in% names(which(table(boxplot_data$IDENT) >= 5)),]

# Calculating median values of "ATF4" for each group
median_values <- tapply(boxplot_data$ATF4, boxplot_data$IDENT, median)

# Sorting the metadata based on median values
boxplot_data$IDENT <- factor(boxplot_data$IDENT, levels = names(sort(median_values, decreasing = F)))

# Plot
p <- ggplot(boxplot_data, aes(x = IDENT, y = ATF4, fill = IDENT)) +
  geom_boxplot() +
  labs(x = NULL, y = "ATF4 Normalized expression") +
  theme_minimal() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  NoLegend()

p

#ggsave(p, filename = paste0(output_figure_prefix, "Martirosyan_ATF4_ordered_boxplot_IDENT_Neurons.pdf"), width = 5, height = 5, bg = "white")
#ggsave(p, filename = paste0(output_figure_prefix, "Martirosyan_ATF4_ordered_boxplot_IDENT_Neurons.png"), width = 5, height = 5, dpi = 1000, bg = 'white')

IV.6.b Plot UMAP ATF4

FeaturePlot(data.seurat, reduction = "UMAP", features = "ATF4")