v4 PA_mod0.6x annotation - SoupX

libraries & paths

library(dplyr)
library(Seurat)
library(patchwork)
library(knitr)
library(dittoSeq)
library(reshape2)
Sys.setenv(RSTUDIO_PANDOC="/usr/local/biotools/pandoc/3.1.2/bin")

sampleID <- c("v4_PA_mod0.6x_SoupX")
sample_color.panel <- c("#009E73")
color.panel <- dittoColors()

Read in object

# read object
dataObject <- readRDS(file = paste0("../../../rObjects/", sampleID, ".filtered.rds"))
markers.strict <- readRDS(file = paste0("../../../rObjects/", 
                          sampleID, "_FindAllMarkers_strict_logFC1_FDR0.05.rds"))

Unannotated

UMAP

ditto_umap <- dittoDimPlot(object = dataObject,
             var = "seurat_clusters",
             reduction.use = "umap",
             do.label = TRUE,
             labels.highlight = TRUE)
ditto_umap

Nuclei count per cluster

count_per_cluster <- FetchData(dataObject,
                               vars = c("ident", "orig.ident")) %>%
  dplyr::count(ident, orig.ident) %>%
  tidyr::spread(ident, n)
count_per_cluster

##     orig.ident    0    1    2    3   4   5   6   7   8   9  10  11  12  13  14
## 1 v4PA_mod0.6x 1843 1370 1203 1112 980 913 599 535 517 509 496 485 436 374 300
##    15  16  17  18  19  20  21 22 23 24 25
## 1 257 227 172 159 153 146 108 82 55 37 29

count_melt <- reshape2::melt(count_per_cluster)
colnames(count_melt) <- c("ident", "cluster", "number of nuclei")
count_max <- count_melt[which.max(count_melt$`number of nuclei`), ]
count_max_value <- count_max$`number of nuclei`
cellmax <- count_max_value + 200 # so that the figure doesn't cut off the text
count_bar <- ggplot(count_melt, aes(x = factor(cluster), y = `number of nuclei`, fill = `ident`)) +
  geom_bar(
    stat = "identity",
    colour = "black",
    width = 1,
    position = position_dodge(width = 0.8)
  ) +
  geom_text(
    aes(label = `number of nuclei`),
    position = position_dodge(width = 0.9),
    vjust = -0.25,
    angle = 45,
    hjust = -.01
  ) +
  theme_classic() + scale_fill_manual(values = sample_color.panel) +
  ggtitle("Number of nuclei per cluster") +  xlab("cluster") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_y_continuous(limits = c(0, cellmax))
count_bar

Cluster tree

dataObject <- BuildClusterTree(
  object = dataObject,
  dims = 1:30,
  reorder = FALSE,
  reorder.numeric = FALSE
)

tree <- dataObject@tools$BuildClusterTree
tree$tip.label <- paste0("Cluster ", tree$tip.label)
nClusters <- length(tree$tip.label)

tree_graph <- ggtree::ggtree(tree, aes(x, y)) +
  scale_y_reverse() +
  ggtree::geom_tree() +
  ggtree::theme_tree() +
  ggtree::geom_tiplab(offset = 1) +
  ggtree::geom_tippoint(color = color.panel[1:nClusters],
                        shape = 16,
                        size = 5) +
  coord_cartesian(clip = 'off') +
  theme(plot.margin = unit(c(0, 2.5, 0, 0), 'cm'))
tree_graph

Violins - Canonical cell-type markers

Markers obtained from dropviz

# Astrocytes
VlnPlot(dataObject,
        features = c("AQP4", "GJA1", "CLU", "GFAP"))

# Endothelial
VlnPlot(dataObject,
        features = c("CLDN5", "ADGRF5", "FLT1"))

# Fibroblast-Like_Dcn
VlnPlot(dataObject,
        features = c("COL1A1", "COL1A2", "DCN"))

# Microglia / Marchophage 
VlnPlot(dataObject,
        features = c("C1QA", "C1QC", "C1QB"))

# Microglia / Marchophage 
VlnPlot(dataObject,
        features = c( "HEXB", "TMEM119", "ITGAM", "TYROBP","P2RY12", "AIF1"))

# Neurons
VlnPlot(dataObject,
        features = c("RBFOX3", "SNAP25","SYT1", "GAD1", "GAD2"))

# Oligodendrocyte
VlnPlot(dataObject,
        features = c("PLP1","MBP", "MOG"))

# OPC
VlnPlot(dataObject,
        features = c("OLIG1","PDGFRA", "VCAN", "TNR"))

# Smooth muscle 
VlnPlot(dataObject,
        features = c("ACTA2","RGS5", "VTN", "MYL5"))

Percent cell type - Canonical cell-type markers

# astrocyte
dataObject[["percent.astrocyte"]] <- PercentageFeatureSet(dataObject, 
                                     features = c("CLU", "GFAP", "AQP4", "GJA1"))
# endothelial
dataObject[["percent.endothelial"]] <- PercentageFeatureSet(dataObject, 
                                       features = c("CLDN5", "ADGRF5", "FLT1"))
# fibroblast 
dataObject[["percent.fibroblast"]] <- PercentageFeatureSet(dataObject, 
                                      features = c("COL1A1", "COL1A2", "DCN"))
# microglia
dataObject[["percent.microglia"]] <- PercentageFeatureSet(dataObject, 
                                     features = c("HEXB", "C1QA", "C1QC", "C1QB", 
                                     "TMEM119", "ITGAM", "TYROBP","P2RY12", "AIF1")) 
# neuron 
dataObject[["percent.neurons"]] <- PercentageFeatureSet(dataObject, 
                                   features = c("RBFOX3", "SNAP25","SYT1", "GAD1", "GAD2"))
# oligodendrocyte
dataObject[["percent.oligodendrocyte"]] <- PercentageFeatureSet(dataObject, 
                                           features = c("PLP1","MBP", "MOG"))
# oligodendrocyte precursor cells
dataObject[["percent.opc"]] <- PercentageFeatureSet(dataObject, 
                               features = c("OLIG1","PDGFRA", "VCAN", "TNR"))
# smooth muscle 
dataObject[["percent.smooth"]] <- PercentageFeatureSet(dataObject, 
                                  features = c("ACTA2","RGS5", "VTN", "MYL5"))

Feature plot percent cell type - Canonical cell-type markers

# astrocyte
FeaturePlot(dataObject, features = "percent.astrocyte", label = TRUE)  

# endothelial
FeaturePlot(dataObject, features = "percent.endothelial", label = TRUE) 

# fibroblast
FeaturePlot(dataObject, features = "percent.fibroblast", label = TRUE) 

# microglia
FeaturePlot(dataObject, features = "percent.microglia", label = TRUE)  

# neuron 
FeaturePlot(dataObject, features = "percent.neurons", label = TRUE) 

# oligodendrocyte
FeaturePlot(dataObject, features = "percent.oligodendrocyte", label = TRUE) 

# oligodendrocyte precursor cells
FeaturePlot(dataObject, features = "percent.opc", label = TRUE)  

# smooth
FeaturePlot(dataObject, features = "percent.smooth", label = TRUE) 

Markers per cluster

Get markers for each cluster

# unique clusters variable
unique_clusters <- unique(markers.strict$cluster)

# empty list to store individual cluster data frames
cluster_list <- list()

# loop through each cluster and create a data frame
for (i in unique_clusters) {
  cluster_name <- paste0("cluster", i)
  cluster_data <- markers.strict[markers.strict$cluster == i, ]
  assign(cluster_name, cluster_data)
  cluster_list[[cluster_name]] <- cluster_data
}

Feature plot

# UMAP showing the expression of select features
umap_feature <-
  FeaturePlot(dataObject,
              features = c("TYROBP", "MOG", "AQP4", "RBFOX3"))
umap_feature

Cluster Annotation

Cluster 0 - oligodendrocyte

# Number of cells per condition
count_per_cluster[,c(1,2)]

##     orig.ident    0
## 1 v4PA_mod0.6x 1843

# UMAP with only cluster 0
DimPlot(object = subset(dataObject, seurat_clusters == "0"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = color.panel[1])

VlnPlot(dataObject,
        features = cluster0$gene[1:10],
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

cluster0$gene[1:10]

##  [1] "PLP1"    "ST18"    "CTNNA3"  "TMEM144" "MOBP"    "MBP"     "TF"     
##  [8] "RNF220"  "CNP"     "ENPP2"

Cluster 1 - astrocyte

count_per_cluster[,c(1,3)]

##     orig.ident    1
## 1 v4PA_mod0.6x 1370

DimPlot(object = subset(dataObject, seurat_clusters == "1"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = color.panel[2])

VlnPlot(dataObject,
        features = cluster1$gene[1:10],
        cols = color.panel,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

cluster1$gene[1:10]

##  [1] "OBI1-AS1"      "ADGRV1"        "LINC00499"     "RP11-134O15.3"
##  [5] "GLIS3"         "PPP1R9A-AS1"   "NHSL1"         "ATP1A2"       
##  [9] "RANBP3L"       "CABLES1"

Cluster 2 - polydendrocyte

count_per_cluster[,c(1,4)]

##     orig.ident    2
## 1 v4PA_mod0.6x 1203

DimPlot(object = subset(dataObject, seurat_clusters == "2"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = color.panel[3])

VlnPlot(dataObject,
        features = cluster2$gene[1:10],
        cols = color.panel,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

cluster2$gene[1:10]

##  [1] "VCAN"         "RP4-668E10.4" "LHFPL3"       "SMOC1"        "MEGF11"      
##  [6] "NXPH1"        "GRIK1"        "PTPRZ1"       "TMEM132C"     "SOX6"

Cluster 3 - neuron

count_per_cluster[,c(1,4)]

##     orig.ident    2
## 1 v4PA_mod0.6x 1203

DimPlot(object = subset(dataObject, seurat_clusters == "3"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols =  color.panel[4])

VlnPlot(dataObject,
        features = cluster3$gene[1:10],
        cols = color.panel,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

cluster3$gene[1:10]

##  [1] "AC114765.1"   "CDH12"        "AC067956.1"   "CBLN2"        "RP4-809F18.1"
##  [6] "RP11-191L9.4" "CLSTN2"       "PDZRN4"       "LY86-AS1"     "RALYL"

Cluster 4 - neuron

count_per_cluster[,c(1,5)]

##     orig.ident    3
## 1 v4PA_mod0.6x 1112

DimPlot(object = subset(dataObject, seurat_clusters == "4"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = color.panel[5])

VlnPlot(dataObject,
        features = cluster4$gene[1:10],
        cols = color.panel,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

cluster4$gene[1:10]

##  [1] "CTC-340A15.2" "CTC-535M15.2" "CPNE4"        "CLSTN2"       "RP11-320L2.1"
##  [6] "RORB"         "NRG1"         "HTR1F"        "IL1RAPL2"     "RP11-197K6.1"

Cluster 5 - neuron

count_per_cluster[,c(1,6)]

##     orig.ident   4
## 1 v4PA_mod0.6x 980

DimPlot(object = subset(dataObject, seurat_clusters == "5"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = color.panel[6])

VlnPlot(dataObject,
        features = cluster5$gene[1:10],
        cols = color.panel,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

cluster5$gene[1:10]

##  [1] "CTC-552D5.1" "CBLN2"       "AC011288.2"  "CDH9"        "GRM1"       
##  [6] "NDST3"       "EPHA6"       "CNTN5"       "TRPC6"       "IQCJ-SCHIP1"

Cluster 6 - neuron - RP high

count_per_cluster[,c(1,7)]

##     orig.ident   5
## 1 v4PA_mod0.6x 913

DimPlot(object = subset(dataObject, seurat_clusters == "6"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = "#CC79A7")

VlnPlot(dataObject,
        features = cluster6$gene[1:10],
        cols = color.panel,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

cluster6$gene[1:10]

##  [1] "RP11-170M17.1" "CUX2"          "RP11-213B3.1"  "COL5A2"       
##  [5] "RP11-78F17.1"  "RP11-33A14.1"  "LINC02822"     "LINC00326"    
##  [9] "AJ006998.2"    "CNGB1"

Cluster 7 - neuron - noise

count_per_cluster[,c(1,8)]

##     orig.ident   6
## 1 v4PA_mod0.6x 599

DimPlot(object = subset(dataObject, seurat_clusters == "7"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = color.panel[7])

VlnPlot(dataObject,
        features = cluster7$gene[1:10],
        cols = color.panel,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

cluster7$gene[1:10]

##  [1] "DLX6-AS1"  "RGS12"     "ADARB2"    "PWRN1"     "SYNPR"     "BTBD11"   
##  [7] "CALB2"     "VWC2L"     "PROX1"     "NR2F2-AS1"

Cluster 8 - neuron

count_per_cluster[,c(1,9)]

##     orig.ident   7
## 1 v4PA_mod0.6x 535

DimPlot(object = subset(dataObject, seurat_clusters == "8"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = color.panel[8])

VlnPlot(dataObject,
        features = cluster8$gene[1:10],
        cols = color.panel,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

cluster8$gene[1:10]

##  [1] "HS3ST4"        "CTD-2337L2.1"  "PCDH11X"       "EGFEM1P"      
##  [5] "FRMD3"         "MSC-AS1"       "KIAA1217"      "ZBTB7C"       
##  [9] "RP11-563D10.1" "LINC02232"

Cluster 9 - neuron

count_per_cluster[,c(1,10)]

##     orig.ident   8
## 1 v4PA_mod0.6x 517

DimPlot(object = subset(dataObject, seurat_clusters == "9"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = color.panel[9])

VlnPlot(dataObject,
        features = cluster9$gene[1:10],
        cols = color.panel,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

cluster9$gene[1:10]

##  [1] "HS3ST4"        "CTD-2337L2.1"  "EGFEM1P"       "SEMA3E"       
##  [5] "FOXP2"         "RP11-640F22.1" "ADAMTSL1"      "RP11-17M24.3" 
##  [9] "CTC-304I17.6"  "SULF1"

Cluster 10 - astrocyte

count_per_cluster[,c(1,11)]

##     orig.ident   9
## 1 v4PA_mod0.6x 509

DimPlot(object = subset(dataObject, seurat_clusters == "10"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = color.panel[10])

VlnPlot(dataObject,
        features = cluster10$gene[1:10],
        cols = color.panel,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

cluster10$gene[1:10]

##  [1] "RP11-6L16.1"   "ADGRV1"        "GLIS3"         "RP11-134O15.3"
##  [5] "GFAP"          "RP11-627D16.1" "RFX4"          "AQP4"         
##  [9] "LINC00299"     "PARD3B"

Cluster 11 - interneuron - noise

count_per_cluster[,c(1,12)]

##     orig.ident  10
## 1 v4PA_mod0.6x 496

DimPlot(object = subset(dataObject, seurat_clusters == "11"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = color.panel[11])

VlnPlot(dataObject,
        features = cluster11$gene[1:10],
        cols = color.panel,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

cluster11$gene[1:10]

##  [1] "NXPH1"         "ZNF385D"       "RP11-123O10.3" "GRIK1"        
##  [5] "PTCHD4"        "SAMD5"         "SLIT2"         "GAD2"         
##  [9] "BTBD11"        "PAM"

Cluster 12 - interneuron

count_per_cluster[,c(1,13)]

##     orig.ident  11
## 1 v4PA_mod0.6x 485

DimPlot(object = subset(dataObject, seurat_clusters == "12"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = color.panel[12])

VlnPlot(dataObject,
        features = cluster12$gene[1:10],
        cols = color.panel,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

cluster12$gene[1:10]

##  [1] "CH507-528H12.1" "CTC-552D5.1"    "CH507-513H4.1"  "CUX2"          
##  [5] "CBLN2"          "AC011288.2"     "EPHA6"          "TRPC6"         
##  [9] "GRM1"           "KCNIP4"

Cluster 13 - neuron

count_per_cluster[,c(1,14)]

##     orig.ident  12
## 1 v4PA_mod0.6x 436

DimPlot(object = subset(dataObject, seurat_clusters == "13"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = color.panel[13])

VlnPlot(dataObject,
        features = cluster13$gene[1:10],
        cols = color.panel,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

cluster13$gene[1:10]

##  [1] "CH507-513H4.1"  "CH507-528H12.1" "CLSTN2"         "CTC-535M15.2"  
##  [5] "CTC-340A15.2"   "FRMPD4"         "HS3ST4"         "RYR2"          
##  [9] "SLIT3"          "MSC-AS1"

Cluster 14 - interneuron

count_per_cluster[,c(1,15)]

##     orig.ident  13
## 1 v4PA_mod0.6x 374

DimPlot(object = subset(dataObject, seurat_clusters == "14"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = color.panel[14])

VlnPlot(dataObject,
        features = cluster14$gene[1:10],
        cols = color.panel,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

cluster14$gene[1:10]

##  [1] "FGF13"         "NXPH1"         "GRIK1"         "PTCHD4"       
##  [5] "ALK"           "EYA4"          "RP11-123O10.3" "GAD2"         
##  [9] "UBASH3B"       "KIT"

Cluster 15 - neuron

count_per_cluster[,c(1,16)]

##     orig.ident  14
## 1 v4PA_mod0.6x 300

DimPlot(object = subset(dataObject, seurat_clusters == "15"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = color.panel[15])

VlnPlot(dataObject,
        features = cluster15$gene[1:10],
        cols = color.panel,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

cluster15$gene[1:10]

##  [1] "CH507-513H4.1"  "CH507-528H12.1" "THEMIS"         "RP11-191L9.4"  
##  [5] "AC067956.1"     "AC114765.1"     "CDH12"          "LINC02254"     
##  [9] "CBLN2"          "CNTN4"

Cluster 16 - neuron

count_per_cluster[,c(1,17)]

##     orig.ident  15
## 1 v4PA_mod0.6x 257

DimPlot(object = subset(dataObject, seurat_clusters == "16"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = color.panel[16])

VlnPlot(dataObject,
        features = cluster16$gene[1:10],
        cols = color.panel,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

cluster16$gene[1:10]

##  [1] "COL5A2"        "AJ006998.2"    "SYT2"          "NEFM"         
##  [5] "NEFH"          "SCN1A-AS1"     "ATP6V1C2"      "RP11-153I24.5"
##  [9] "RP11-335E8.3"  "GRM8"

Cluster 17 - astrocyte & oligodenrocyte?

count_per_cluster[,c(1,18)]

##     orig.ident  16
## 1 v4PA_mod0.6x 227

DimPlot(object = subset(dataObject, seurat_clusters == "17"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = color.panel[17])

VlnPlot(dataObject,
        features = cluster17$gene[1:10],
        cols = color.panel,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

cluster17$gene[1:10]

##  [1] "LINC00499"     "ADGRV1"        "PPP1R9A-AS1"   "OBI1-AS1"     
##  [5] "GLIS3"         "PRODH"         "GPC5"          "SLC1A2"       
##  [9] "AC084149.2"    "RP11-134O15.3"

Cluster 18 - neuron

count_per_cluster[,c(1,1)]

##     orig.ident orig.ident.1
## 1 v4PA_mod0.6x v4PA_mod0.6x

DimPlot(object = subset(dataObject, seurat_clusters == "18"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = color.panel[19])

VlnPlot(dataObject,
        features = cluster18$gene[1:10],
        cols = color.panel,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

cluster18$gene[1:10]

##  [1] "ATP10A"        "RP11-47B24.1"  "RP11-632B21.2" "SMYD1"        
##  [5] "NR4A2"         "STPG2-AS1"     "PLA2G4A"       "LINC00683"    
##  [9] "GALNT14"       "RP11-739G5.1"

Cluster 19 - neuron

count_per_cluster[,c(1,19)]

##     orig.ident  17
## 1 v4PA_mod0.6x 172

DimPlot(object = subset(dataObject, seurat_clusters == "19"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = color.panel[20])

VlnPlot(dataObject,
        features = cluster19$gene[1:10],
        cols = color.panel,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

cluster19$gene[1:10]

##  [1] "CXCL14"        "RELN"          "NR2F2-AS1"     "NR2F2"        
##  [5] "DLX6-AS1"      "RP11-406A20.1" "RP11-12B13.1"  "GAD2"         
##  [9] "CHRNA7"        "GRIK1"

Cluster 20 - neuron

count_per_cluster[,c(1,20)]

##     orig.ident  18
## 1 v4PA_mod0.6x 159

DimPlot(object = subset(dataObject, seurat_clusters == "20"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = color.panel[21])

VlnPlot(dataObject,
        features = cluster20$gene[1:10],
        cols = color.panel,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

cluster20$gene[1:10]

##  [1] "HTR2C"         "NPSR1-AS1"     "IFNG-AS1"      "FER1L6-AS2"   
##  [5] "NXPH2"         "ENPP7P1"       "LINC02218"     "CD200R1L"     
##  [9] "RP11-354I13.2" "CHRM2"

Cluster 21 - microglia

count_per_cluster[,c(1,21)]

##     orig.ident  19
## 1 v4PA_mod0.6x 153

DimPlot(object = subset(dataObject, seurat_clusters == "21"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = color.panel[22])

VlnPlot(dataObject,
        features = cluster21$gene[1:10],
        cols = color.panel,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

cluster21$gene[1:10]

##  [1] "DOCK8"    "FYB1"     "LNCAROD"  "TBXAS1"   "ADAM28"   "INPP5D"  
##  [7] "ARHGAP15" "APBB1IP"  "SLCO2B1"  "CD74"

Cluster 22 - neuron - noise

count_per_cluster[,c(1,22)]

##     orig.ident  20
## 1 v4PA_mod0.6x 146

DimPlot(object = subset(dataObject, seurat_clusters == "22"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = color.panel[23])

VlnPlot(dataObject,
        features = cluster22$gene[1:10],
        cols = color.panel,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

cluster22$gene[1:10]

##  [1] "SCUBE3"        "NOG"           "RP11-238K6.2"  "LHX6"         
##  [5] "LINC00943"     "NPNT"          "SCN1A-AS1"     "RP11-519H15.1"
##  [9] "RP11-407A16.3" "ANK1"

Cluster 23 - endothelial

count_per_cluster[,c(1,23)]

##     orig.ident  21
## 1 v4PA_mod0.6x 108

DimPlot(object = subset(dataObject, seurat_clusters == "23"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = color.panel[24])

VlnPlot(dataObject,
        features = cluster23$gene[1:10],
        cols = color.panel,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

cluster23$gene[1:10]

##  [1] "EBF1"          "IFITM3"        "RP11-326C3.17" "A2M"          
##  [5] "FLT1"          "CFH"           "SLC6A12"       "IFITM1"       
##  [9] "NOTCH3"        "LEF1"

Cluster 24 - polydendrocyte

count_per_cluster[,c(1,24)]

##     orig.ident 22
## 1 v4PA_mod0.6x 82

DimPlot(object = subset(dataObject, seurat_clusters == "24"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = color.panel[25])

VlnPlot(dataObject,
        features = cluster24$gene[1:10],
        cols = color.panel,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

cluster24$gene[1:10]

##  [1] "RP4-668E10.4"  "LHFPL3"        "VCAN"          "RP1-45C12.1"  
##  [5] "SMOC1"         "RP11-111G13.1" "COL9A1"        "MOBP"         
##  [9] "PLP1"          "CTNNA3"

Cluster 25 - microglia

count_per_cluster[,c(1,25)]

##     orig.ident 23
## 1 v4PA_mod0.6x 55

DimPlot(object = subset(dataObject, seurat_clusters == "25"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = color.panel[26])

VlnPlot(dataObject,
        features = cluster25$gene[1:10],
        cols = color.panel,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

cluster25$gene[1:10]

##  [1] "ADAM28"  "LNCAROD" "SYK"     "FYB1"    "CD74"    "CSF1R"   "TBXAS1" 
##  [8] "DOCK8"   "PTPRC"   "APBB1IP"

Assign identities

Individual

dataObject.annotated <- RenameIdents(object = dataObject, 
                               "0" = "oligodendrocyte 1",
                               "1" = "astrocyte 1",
                               "2" = "polydendrocyte 1",
                               "3" = "neuron 1",
                               "4" = "neuron 2",
                               "5" = "neuron 3",
                               "6" = "neuron 4",
                               "7" = "neuron/noise 1",
                               "8" = "neuron 6",
                               "9" = "neuron 7",
                               "10" = "astrocyte 2",
                               "11" = "interneuron/noise",
                               "12" = "interneuron 2",
                               "13" = "neuron 8",
                               "14" = "interneuron 3",
                               "15" = "neuron 9",
                               "16" = "neuron 10",
                               "17" = "astrocyte/noise",
                               "18" = "neuron 11",
                               "19" = "neuron 12",
                               "20" = "neuron 13",
                               "21" = "microglia 1",
                               "22" = "neuron/noise 2",
                               "23" = "endothelial",
                               "24" = "polydendrocyte 2",
                               "25" = "microglia 2")
dataObject.annotated$individual_clusters <- factor(Idents(dataObject.annotated))

UMAP_ind <- dittoDimPlot(object = dataObject.annotated,
             var = "individual_clusters",
             reduction.use = "umap",
             do.label = TRUE,
             labels.highlight = TRUE)
UMAP_ind

## png 
##   2

Nuclei count per cluster

count_per_cluster <- FetchData(dataObject.annotated,
                               vars = c("ident", "orig.ident")) %>%
  dplyr::count(ident, orig.ident) %>%
  tidyr::spread(ident, n)
count_per_cluster

##     orig.ident oligodendrocyte 1 astrocyte 1 polydendrocyte 1 neuron 1 neuron 2
## 1 v4PA_mod0.6x              1843        1370             1203     1112      980
##   neuron 3 neuron 4 neuron/noise 1 neuron 6 neuron 7 astrocyte 2
## 1      913      599            535      517      509         496
##   interneuron/noise interneuron 2 neuron 8 interneuron 3 neuron 9 neuron 10
## 1               485           436      374           300      257       227
##   astrocyte/noise neuron 11 neuron 12 neuron 13 microglia 1 neuron/noise 2
## 1             172       159       153       146         108             82
##   endothelial polydendrocyte 2 microglia 2
## 1          55               37          29

count_melt <- reshape2::melt(count_per_cluster)
colnames(count_melt) <- c("ident", "cluster", "number of nuclei")
count_max <- count_melt[which.max(count_melt$`number of nuclei`), ]
count_max_value <- count_max$`number of nuclei`
cellmax <- count_max_value + 500 # so that the figure doesn't cut off the text
count_bar <- ggplot(count_melt, aes(x = factor(cluster), y = `number of nuclei`, fill = `cluster`)) +
  geom_bar(
    stat = "identity",
    colour = "black",
    width = 1,
    position = position_dodge(width = 0.8)
  ) +
  geom_text(
    aes(label = `number of nuclei`),
    position = position_dodge(width = 0.9),
    vjust = -0.25,
    angle = 45,
    hjust = -.01
  ) +
  theme_classic() + scale_fill_manual(values = color.panel) +
  ggtitle("Number of nuclei per cluster") +  xlab("cluster") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_y_continuous(limits = c(0, cellmax))
count_bar

## png 
##   2

DotPlot

markers.to.plot <-
  c(
"CLU", 
"GFAP", 
"AQP4", 
"GJA1",
"CLDN5", 
"ADGRF5", 
"FLT1",
"COL1A1", 
"COL1A2", 
"DCN",
"HEXB", 
"C1QA", 
"C1QC", 
"C1QB", 
"TMEM119", 
"ITGAM", 
"TYROBP",
"P2RY12", 
"AIF1",
"RBFOX3", 
"SNAP25",
"SYT1", 
"GAD1", 
"GAD2",
"PLP1",
"MBP", 
"MOG",
"OLIG1",
"PDGFRA", 
"VCAN", 
"TNR",
"ACTA2",
"RGS5", 
"VTN", 
"MYL5"
  )

dot_ind <- DotPlot(dataObject, 
                   features = markers.to.plot, 
                   cluster.idents = TRUE,
                   dot.scale = 8) + RotatedAxis()
dot_ind

## png 
##   2

Merged

dataObject.annotated <- RenameIdents(object = dataObject.annotated, 
"oligodendrocyte 1" = "oligodendrocyte",
"astrocyte 1" = "astrocyte",
"polydendrocyte 1" = "polydendrocyte",
"neuron 1" = "neuron",
"neuron 2" = "neuron",
"neuron 3" = "neuron",
"neuron 4" = "neuron",
"neuron/noise 1" = "neuron/noise",
"neuron 6" = "neuron",
"neuron 7" = "neuron",
"astrocyte 2" = "astrocyte",
"interneuron/noise" = "interneuron/noise",
"interneuron 2" = "interneuron",
"neuron 8" = "neuron",
"interneuron 3" = "interneuron",
"neuron 9" = "neuron",
"neuron 10" = "neuron",
"astrocyte/noise" = "astrocyte/noise",
"neuron 11" = "neuron",
"neuron 12" = "neuron",
"neuron 13" = "neuron",
"microglia 1" = "microglia",
"neuron/noise 2" = "neuron/noise",
"endothelial" = "endothelial",
"polydendrocyte 2" = "polydendrocyte",
"microglia 2" = "microglia")
dataObject.annotated$annotated_clusters <- factor(Idents(dataObject.annotated))
Idents(dataObject.annotated) <- "annotated_clusters"

UMAP_merge <- dittoDimPlot(object = dataObject.annotated,
             var = "annotated_clusters",
             reduction.use = "umap",
             do.label = TRUE,
             labels.highlight = TRUE)
UMAP_merge

## png 
##   2

Nuclei count per cluster

count_per_cluster <- FetchData(dataObject.annotated,
                               vars = c("ident", "orig.ident")) %>%
  dplyr::count(ident, orig.ident) %>%
  tidyr::spread(ident, n)
count_per_cluster

##     orig.ident oligodendrocyte astrocyte polydendrocyte neuron neuron/noise
## 1 v4PA_mod0.6x            1843      1866           1240   5946          617
##   interneuron/noise interneuron astrocyte/noise microglia endothelial
## 1               485         736             172       137          55

count_melt <- reshape2::melt(count_per_cluster)
colnames(count_melt) <- c("ident", "cluster", "number of nuclei")
count_max <- count_melt[which.max(count_melt$`number of nuclei`), ]
count_max_value <- count_max$`number of nuclei`
cellmax <- count_max_value + 500 # so that the figure doesn't cut off the text
count_bar <- ggplot(count_melt, aes(x = factor(cluster), y = `number of nuclei`, fill = `cluster`)) +
  geom_bar(
    stat = "identity",
    colour = "black",
    width = 1,
    position = position_dodge(width = 0.8)
  ) +
  geom_text(
    aes(label = `number of nuclei`),
    position = position_dodge(width = 0.9),
    vjust = -0.25,
    angle = 45,
    hjust = -.01
  ) +
  theme_classic() + scale_fill_manual(values = color.panel) +
  ggtitle("Number of nuclei per cluster") +  xlab("cluster") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_y_continuous(limits = c(0, cellmax))
count_bar

## png 
##   2

Relative abundance

relative_abundance <- dataObject.annotated@meta.data %>%
  group_by(annotated_clusters, orig.ident) %>%
  dplyr::count() %>%
  group_by(orig.ident) %>%
  dplyr::mutate(percent = 100 * n / sum(n)) %>%
  ungroup()


rel_abun <- ggplot(relative_abundance, aes(x = orig.ident, y = percent, fill = annotated_clusters)) +
  geom_col() +
  geom_text(aes(label = paste0(round(percent), "%")), 
            position = position_stack(vjust = 0.5), size = 3, color = "white") +
  scale_fill_manual(values = color.panel) +
  ggtitle("Percentage of cell type") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))
rel_abun

DotPlot

markers.to.plot <-
  c(
"CLU", 
"GFAP", 
"AQP4", 
"GJA1",
"CLDN5", 
"ADGRF5", 
"FLT1",
"COL1A1", 
"COL1A2", 
"DCN",
"HEXB", 
"C1QA", 
"C1QC", 
"C1QB", 
"TMEM119", 
"ITGAM", 
"TYROBP",
"P2RY12", 
"AIF1",
"RBFOX3", 
"SNAP25",
"SYT1", 
"GAD1", 
"GAD2",
"PLP1",
"MBP", 
"MOG",
"OLIG1",
"PDGFRA", 
"VCAN", 
"TNR",
"ACTA2",
"RGS5", 
"VTN", 
"MYL5"
  )
dot_merge <- DotPlot(dataObject.annotated, 
               features = markers.to.plot, 
               cluster.idents = TRUE,
               dot.scale = 8) + RotatedAxis()
dot_merge

## png 
##   2

Save RDS

saveRDS(dataObject.annotated, paste0("../../../rObjects/", sampleID, ".annotated.rds"))