Skip to content
Snippets Groups Projects
Commit 19c00ea7 authored by John Lafin's avatar John Lafin
Browse files

Complete annotation script

parent e72babe6
Branches
No related merge requests found
# Annotation and stress score filtering # Annotation and stress score filtering
# Author: John T Lafin # Author: John T Lafin
# Updated: 2024-10-15 # Updated: 2024-10-16
print("Script begin at: ") print("Script begin at: ")
Sys.time() Sys.time()
...@@ -70,12 +70,12 @@ dir.create(annotdir, showWarnings=FALSE) ...@@ -70,12 +70,12 @@ dir.create(annotdir, showWarnings=FALSE)
## Create list of labels ## Create list of labels
annot.list <- list() annot.list <- list()
annot.list$Endothelia <- c("Venous endothelia", annot.list$Endothelia <- c("venous endothelia",
"Venous endothelia", "venous endothelia",
"Capillary endothelia", "capillary endothelia",
"Capillary endothelia", "capillary endothelia",
"Arterial endothelia", "arterial endothelia",
"Lymphatic endothelia") "lymphatic endothelia")
annot.list$Epithelia <- c("transitional basal epithelia of urethra", annot.list$Epithelia <- c("transitional basal epithelia of urethra",
"basal epithelia of prostate", "basal epithelia of prostate",
"transitional basal epithelia of urethra", "transitional basal epithelia of urethra",
...@@ -92,7 +92,7 @@ annot.list$Epithelia <- c("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",
"squamous intermediate epithelia of urethra", "squamous intermediate epithelia of urethra",
"basal epithelia of prostate", "basal epithelia of prostate",
"transitional basal epithelia of urethra", "squamous basal epithelia of urethra",
"basal epithelia of prostate", "basal epithelia of prostate",
"proliferative cells", "proliferative cells",
"squamous luminal epithelia of urethra", "squamous luminal epithelia of urethra",
...@@ -100,15 +100,15 @@ annot.list$Epithelia <- c("transitional basal epithelia of urethra", ...@@ -100,15 +100,15 @@ annot.list$Epithelia <- c("transitional basal epithelia of urethra",
"transitional basal epithelia of urethra", "transitional basal epithelia of urethra",
"basal epithelia of prostate") "basal epithelia of prostate")
annot.list$Leukocyte <- c("T cell", annot.list$Leukocyte <- c("T cell",
"Macrophage", "macrophage",
"Mast cell", "mast cell",
"T cell", "T cell",
"T cell", "T cell",
"NK cell", "NK cell",
"Proliferating leukocyte", "proliferating leukocyte",
"B cell", "B cell",
"Dendritic cell", "dendritic cell",
"Plasma cell") "plasma cell")
annot.list$Stroma <- c("pericyte", annot.list$Stroma <- c("pericyte",
"major interstitial fibroblast", "major interstitial fibroblast",
"smooth muscle", "smooth muscle",
...@@ -150,7 +150,7 @@ VlnPlot(lin.list$Epithelia, "NE.score", group.by = "orig.ident") + ...@@ -150,7 +150,7 @@ VlnPlot(lin.list$Epithelia, "NE.score", group.by = "orig.ident") +
ggsave(paste0(annotdir, "NE_score_cutoff.png")) ggsave(paste0(annotdir, "NE_score_cutoff.png"))
lin.list$Epithelia@meta.data <- lin.list$Epithelia@meta.data %>% 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 ## Add ionocytes
VlnPlot(lin.list$Epithelia, "FOXI1", group.by = "orig.ident") + VlnPlot(lin.list$Epithelia, "FOXI1", group.by = "orig.ident") +
...@@ -159,10 +159,7 @@ 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")) ggsave(paste0(annotdir, "FOXI1_score_cutoff.png"))
ionocyte.barcodes <- WhichCells(lin.list$Epithelia, expression = FOXI1 > 0.3) ionocyte.barcodes <- WhichCells(lin.list$Epithelia, expression = FOXI1 > 0.3)
lin.list$Epithelia@meta.data[ionocyte.barcodes, "cellannotation"] <- "Ionocyte" lin.list$Epithelia@meta.data[ionocyte.barcodes, "cellannotation"] <- "ionocyte"
## Save updates to object
lin.list$Epithelia <- StashIdent(lin.list$Epithelia, "cellannotation")
## Split periepithelial fibroblasts by sex ## Split periepithelial fibroblasts by sex
fib_rename <- WhichCells(lin.list$Stroma, fib_rename <- WhichCells(lin.list$Stroma,
...@@ -170,7 +167,10 @@ fib_rename <- WhichCells(lin.list$Stroma, ...@@ -170,7 +167,10 @@ fib_rename <- WhichCells(lin.list$Stroma,
expression = Sex == "male") expression = Sex == "male")
levels(lin.list$Stroma$cellannotation) <- c(levels(lin.list$Stroma$cellannotation), "periepithelial fibroblast of prostate") 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" 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 ## Save labeled UMAP plots
for (lin in lins) { for (lin in lins) {
...@@ -192,9 +192,7 @@ Idents(lin.list$Epithelia) <- "RNA_snn_res.1.5" ...@@ -192,9 +192,7 @@ Idents(lin.list$Epithelia) <- "RNA_snn_res.1.5"
for (lin in names(lin.list)) { for (lin in names(lin.list)) {
VlnPlot(lin.list[[lin]], "Stress.score", pt.size=0) + VlnPlot(lin.list[[lin]], "Stress.score", pt.size=0) +
NoLegend() NoLegend()
ggsave(paste0(stressdir, lin, "_pre-filter.png"), ggsave(paste0(stressdir, lin, "_pre-filter.png"))
width=11,
height=8.5)
} }
## Define stressed clusters based on violin plots ## Define stressed clusters based on violin plots
...@@ -209,24 +207,59 @@ for (lin in names(lin.list)) { ...@@ -209,24 +207,59 @@ for (lin in names(lin.list)) {
lin.list[[lin]] <- subset(lin.list[[lin]], lin.list[[lin]] <- subset(lin.list[[lin]],
idents = stressed.cells[[lin]], idents = stressed.cells[[lin]],
invert = TRUE) 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 stressed cell removal in stroma
# Second round of stress filtering on stroma Idents(lin.list$Stroma) <- "RNA_snn_res.1"
# Revisit epithelia to relabel lin.list$Stroma <- subset(lin.list$Stroma,
# Make sure that ionocytes and NE cells are present in epi labels 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 seu <- merge(lin.list[[1]], lin.list[2:length(lin.list)])
for (lin in names(lin.list)) { plan(sequential)
VlnPlot(lin.list[[lin]], "Stress.score", pt.size=0) + seu <- preprocess(seu, default.res="cellannotation")
NoLegend()
ggsave(paste0(stressdir, lin, "_first_filter.png"), ## Save UMAP plot
width=11, DimPlot(seu, label=TRUE, repel=TRUE, label.size=3) +
height=8.5) 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"))
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment