diff --git a/scripts/04_annotation.R b/scripts/04_annotation.R
index 5c1d175c8f0a5c2609bbca99383da55b92ae9228..1e09dc673415291b7f671f82b65f023434640126 100644
--- a/scripts/04_annotation.R
+++ b/scripts/04_annotation.R
@@ -1,6 +1,6 @@
# Annotation and stress score filtering
# Author: John T Lafin
-# Updated: 2024-10-15
+# Updated: 2024-10-16
print("Script begin at: ")
Sys.time()
@@ -70,12 +70,12 @@ dir.create(annotdir, showWarnings=FALSE)
## Create list of labels
annot.list <- list()
-annot.list$Endothelia <- c("Venous endothelia",
- "Venous endothelia",
- "Capillary endothelia",
- "Capillary endothelia",
- "Arterial endothelia",
- "Lymphatic endothelia")
+annot.list$Endothelia <- c("venous endothelia",
+ "venous endothelia",
+ "capillary endothelia",
+ "capillary endothelia",
+ "arterial endothelia",
+ "lymphatic endothelia")
annot.list$Epithelia <- c("transitional basal epithelia of urethra",
"basal epithelia of prostate",
"transitional basal epithelia of urethra",
@@ -92,7 +92,7 @@ annot.list$Epithelia <- c("transitional basal epithelia of urethra",
"squamous intermediate epithelia of urethra",
"squamous intermediate epithelia of urethra",
"basal epithelia of prostate",
- "transitional basal epithelia of urethra",
+ "squamous basal epithelia of urethra",
"basal epithelia of prostate",
"proliferative cells",
"squamous luminal epithelia of urethra",
@@ -100,15 +100,15 @@ annot.list$Epithelia <- c("transitional basal epithelia of urethra",
"transitional basal epithelia of urethra",
"basal epithelia of prostate")
annot.list$Leukocyte <- c("T cell",
- "Macrophage",
- "Mast cell",
+ "macrophage",
+ "mast cell",
"T cell",
"T cell",
"NK cell",
- "Proliferating leukocyte",
+ "proliferating leukocyte",
"B cell",
- "Dendritic cell",
- "Plasma cell")
+ "dendritic cell",
+ "plasma cell")
annot.list$Stroma <- c("pericyte",
"major interstitial fibroblast",
"smooth muscle",
@@ -150,7 +150,7 @@ VlnPlot(lin.list$Epithelia, "NE.score", group.by = "orig.ident") +
ggsave(paste0(annotdir, "NE_score_cutoff.png"))
lin.list$Epithelia@meta.data <- lin.list$Epithelia@meta.data %>%
- mutate(cellannotation = if_else(NE.score > 0.3, "Neuroendocrine cell", cellannotation))
+ mutate(cellannotation = if_else(NE.score > 0.3, "neuroendocrine cell", cellannotation))
## Add ionocytes
VlnPlot(lin.list$Epithelia, "FOXI1", group.by = "orig.ident") +
@@ -159,10 +159,7 @@ VlnPlot(lin.list$Epithelia, "FOXI1", group.by = "orig.ident") +
ggsave(paste0(annotdir, "FOXI1_score_cutoff.png"))
ionocyte.barcodes <- WhichCells(lin.list$Epithelia, expression = FOXI1 > 0.3)
-lin.list$Epithelia@meta.data[ionocyte.barcodes, "cellannotation"] <- "Ionocyte"
-
-## Save updates to object
-lin.list$Epithelia <- StashIdent(lin.list$Epithelia, "cellannotation")
+lin.list$Epithelia@meta.data[ionocyte.barcodes, "cellannotation"] <- "ionocyte"
## Split periepithelial fibroblasts by sex
fib_rename <- WhichCells(lin.list$Stroma,
@@ -170,7 +167,10 @@ fib_rename <- WhichCells(lin.list$Stroma,
expression = Sex == "male")
levels(lin.list$Stroma$cellannotation) <- c(levels(lin.list$Stroma$cellannotation), "periepithelial fibroblast of prostate")
lin.list$Stroma@meta.data[fib_rename, "cellannotation"] <- "periepithelial fibroblast of prostate"
-Idents(stroma) <- "cellannotation"
+
+## Reset idents
+Idents(lin.list$Epithelia) <- "cellannotation"
+Idents(lin.list$Stroma) <- "cellannotation"
## Save labeled UMAP plots
for (lin in lins) {
@@ -192,9 +192,7 @@ Idents(lin.list$Epithelia) <- "RNA_snn_res.1.5"
for (lin in names(lin.list)) {
VlnPlot(lin.list[[lin]], "Stress.score", pt.size=0) +
NoLegend()
- ggsave(paste0(stressdir, lin, "_pre-filter.png"),
- width=11,
- height=8.5)
+ ggsave(paste0(stressdir, lin, "_pre-filter.png"))
}
## Define stressed clusters based on violin plots
@@ -209,24 +207,59 @@ for (lin in names(lin.list)) {
lin.list[[lin]] <- subset(lin.list[[lin]],
idents = stressed.cells[[lin]],
invert = TRUE)
+ lin.list[[lin]] <- preprocess(lin.list[[lin]], default.res="cellannotation")
}
-lin.list <- lapply(lin.list, preprocess)
+## Save violin plots after first round of filtering
+for (lin in names(lin.list)) {
+ VlnPlot(lin.list[[lin]], "Stress.score", pt.size=0, group.by="RNA_snn_res.1") +
+ NoLegend()
+ ggsave(paste0(stressdir, lin, "_first_filter.png"))
+}
-### TO-DO
-# Second round of stress filtering on stroma
-# Revisit epithelia to relabel
-# Make sure that ionocytes and NE cells are present in epi labels
+## Second round of stressed cell removal in stroma
+Idents(lin.list$Stroma) <- "RNA_snn_res.1"
+lin.list$Stroma <- subset(lin.list$Stroma,
+ idents = 8,
+ invert = TRUE) %>%
+ preprocess(default.res="cellannotation")
+## Save violin plot
+VlnPlot(lin.list$Stroma, "Stress.score", pt.size=0, group.by="RNA_snn_res.1") +
+ NoLegend()
+ggsave(paste0(stressdir, lin, "_second_filter.png"))
+## Reannotate endothelia
+lin.list$Endothelia@meta.data <- lin.list$Endothelia@meta.data %>%
+ mutate(cellannotation = if_else(RNA_snn_res.1 == 7, "venous endothelia", cellannotation))
+Idents(lin.list$Endothelia) <- "cellannotation"
+## Save labeled UMAP plots
+for (lin in lins) {
+ DimPlot(lin.list[[lin]])
+ ggsave(paste0(annotdir, lin,"_umap-2.png"))
+}
+# Merge -------------------------------------------------------------------
-# Produce post-filter stress plots
-for (lin in names(lin.list)) {
- VlnPlot(lin.list[[lin]], "Stress.score", pt.size=0) +
- NoLegend()
- ggsave(paste0(stressdir, lin, "_first_filter.png"),
- width=11,
- height=8.5)
+seu <- merge(lin.list[[1]], lin.list[2:length(lin.list)])
+plan(sequential)
+seu <- preprocess(seu, default.res="cellannotation")
+
+## Save UMAP plot
+DimPlot(seu, label=TRUE, repel=TRUE, label.size=3) +
+ NoLegend()
+ggsave(paste0(annotdir, "global_umap.png"))
+
+# Save outputs ------------------------------------------------------------
+
+objdir <- paste0(outdir, "objects/")
+dir.create(objdir, showWarnings=FALSE)
+
+saveRDS(seu, paste0(objdir, "global.rds"))
+for (lin in names(lin.list)){
+ saveRDS(lin.list[[lin]], paste0(objdir, lin, ".rds"))
}
+
+## Save session info
+capture.output(sessionInfo(), file = paste0(outdir, "session_info.txt"))