path.to.sample.with.cell.cycle.regressed <- "path_to_cell_cycle_regressed_object"
sample <- readRDS(path.to.sample.with.cell.cycle.regressed)
path.to.supervised.annotation.set <- "path_to_supervised_annotation_set"
markers <- readRDS(path.to.supervised.annotation.set)
colors.use.reduced <- c("Unannotated" = "#C0C0C0",
"TME" = "#70798c",
"IPC-like" = "#be920e",
"CP-like" = "#be660e",
"Cilia-like" = "#be0e0e",
"Mesenchymal-like" = "#0ebe66",
"RG-like" = "#0497c8",
"NPC-like" = "#0466c8",
"OPC-like" = "#0435c8",
"Immune-like" = "#920ebe",
"Hypoxic" = "#92be0e")
# # Get S phase genes.
s.genes <- Seurat::cc.genes.updated.2019$s.genes
# Get G2-M phase genes.
g2m.genes <- Seurat::cc.genes.updated.2019$g2m.genes
regress_out_vars <- c("nCount_RNA", "nFeature_RNA", "percent.mt")
normalization_batch <- "orig.ident"
integration_batch <- "technology"
# Change the assay to Seurat V5.
sample[["RNA"]] <- as(object = sample[["RNA"]], Class = "Assay5")
regress_out_vars <- c("nCount_RNA", "nFeature_RNA", "percent.mt", "S.Score", "G2M.Score")
normalization_batch <- "orig.ident"
integration_batch <- "technology"
sample[["RNA"]] <- split(x = sample[["RNA"]],
f = sample@meta.data[, normalization_batch])
# Perform normalization.
sample <- Seurat::NormalizeData(sample)
sample <- Seurat::FindVariableFeatures(sample)
sample <- Seurat::ScaleData(sample, vars.to.regress = regress_out_vars)
sample <- Seurat::RunPCA(sample)
sample <- SeuratObject::JoinLayers(sample)
sample@assays[["RNA"]]@layers <- sample@assays[["RNA"]]@layers[c("counts", "data", "scale.data")]
sample <- Seurat::CellCycleScoring(sample, s.features = s.genes, g2m.features = g2m.genes, set.ident = FALSE)
sample[["RNA"]] <- split(x = sample[["RNA"]],
f = sample@meta.data[, normalization_batch])
# Perform normalization.
sample <- Seurat::NormalizeData(sample)
sample <- Seurat::FindVariableFeatures(sample)
sample <- Seurat::ScaleData(sample, vars.to.regress = c(regress_out_vars, "S.Score", "G2M.Score"))
sample <- Seurat::RunPCA(sample)
sample <- SeuratObject::JoinLayers(sample)
sample@assays[["RNA"]]@layers <- sample@assays[["RNA"]]@layers[c("counts", "data", "scale.data")]
sample <- harmony::RunHarmony(sample,
assay = "RNA",
group.by.vars = c(normalization_batch, integration_batch),
theta = c(1, 2))
sample <- Seurat::FindNeighbors(sample, reduction = "harmony", dims = 1:30)
sample <- Seurat::FindClusters(sample, cluster.name = "harmony_clusters")
sample <- Seurat::RunUMAP(sample, reduction = "harmony", dims = 1:30, reduction.name = "umap.harmony", reduction.key = "UMAPHARMONY_")
sample <- Seurat::CellCycleScoring(sample, s.features = s.genes, g2m.features = g2m.genes, set.ident = FALSE)
sample$Annotation <- as.character(sample$Annotation)
sample$Annotation[sample$Annotation == "Tumor"] <- "Unannotated"
sample$Annotation <- factor(sample$Annotation, levels = c("IPC-like", "CP-like", "Cilia-like", "OPC-like", "NPC-like", "RG-like", "Mesenchymal-like", "Hypoxic", "Immune-like", "TME", "Unannotated"))
out <- SCpubr::do_EnrichmentHeatmap(sample, input_gene_list = markers, flavor = "UCell", return_object = TRUE)
sample.use <- out$Object
sample.use <- Seurat::ScaleData(sample.use)
sample.use <- Seurat::RunPCA(sample.use, features = rownames(sample.use))
sample.use <- Seurat::RunUMAP(sample.use, reduction = "pca", dims = 1:30)
sample.use <- Seurat::FindNeighbors(sample.use, reduction = "pca", dims = 1:30)
sample.use <- Seurat::FindClusters(sample.use, reduction = "pca", dims = 1:30)
p1 <- SCpubr::do_DimPlot(sample,
group.by = "Annotation",
reduction = "umap.harmony",
idents.keep = "IPC-like",
colors.use = colors.use.reduced,
label = TRUE,
repel = TRUE,
legend.position = "none",
raster = TRUE,
raster.dpi = 2048,
na.value = "grey90",
pt.size = 8,
font.size = 16,
plot.title = "Current Annotation")
p2 <- SCpubr::do_DimPlot(sample.use,
reduction = "umap",
label = TRUE,
repel = TRUE,
legend.position = "none",
raster = TRUE,
raster.dpi = 2048,
na.value = "grey90",
pt.size = 8,
font.size = 16)
p3 <- SCpubr::do_FeaturePlot(sample.use,
reduction = "umap",
features = c("Neuronal.IPC", "Cycle"),
ncol = 1,
raster = TRUE,
raster.dpi = 2048,
na.value = "grey90",
pt.size = 8,
font.size = 16)
sample.use$New_Annotation <- as.character(sample.use$seurat_clusters)
sample.use$New_Annotation[sample.use$New_Annotation == "14"] <- "IPC-like"
sample$New_Annotation <- sample.use$New_Annotation
p4 <- SCpubr::do_DimPlot(sample,
group.by = "New_Annotation",
reduction = "umap.harmony",
idents.keep = "IPC-like",
colors.use = colors.use,
label = TRUE,
repel = TRUE,
legend.position = "none",
raster = TRUE,
raster.dpi = 2048,
na.value = "grey90",
pt.size = 8,
font.size = 16)
data.use <- list("Old Annotation" = names(sample$Annotation[sample$Annotation == "IPC-like"]),
"New Annotation" = names(sample$New_Annotation[sample$New_Annotation == "IPC-like"]))
p5 <- SCpubr::do_BarPlot(sample = sample,
split.by = "New_Annotation",
group.by = "Annotation",
position = "fill",
flip = TRUE,
colors.use = colors.use.reduced,
legend.ncol = 3,
xlab = "",
order = TRUE,
order.by = "IPC-like",
font.size = 16)
layout <- "ABC
DBE"
p <- patchwork::wrap_plots(A = p1,
B = p3,
C = p4,
D = p2,
E = p5,
design = layout)