v4 PA_mod0.6x annotation - sens 3 with the soup

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_sens3")
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 1518 1275 1163 1064 1050 671 627 580 488 452 450 365 245 219 211
##    15  16  17  18  19  20  21  22 23 24 25 26
## 1 203 171 154 149 144 129 121 115 67 67 26 22

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 1518

# 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] "ST18"          "MOBP"          "TMEM144"       "TF"           
##  [5] "ENPP2"         "SLCO1A2"       "CNP"           "RP11-654K19.6"
##  [9] "CLMN"          "LINC01608"

Cluster 1 - astrocyte

count_per_cluster[,c(1,3)]

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

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] "LINC00499"     "RP11-134O15.3" "OBI1-AS1"      "PPP1R9A-AS1"  
##  [5] "RANBP3L"       "ATP1A2"        "NHSL1"         "GLIS3"        
##  [9] "PRODH"         "CABLES1"

Cluster 2 - neuron

count_per_cluster[,c(1,4)]

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

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] "AC114765.1"   "AC067956.1"   "CBLN2"        "CDH12"        "RP4-809F18.1"
##  [6] "RP11-191L9.4" "LY86-AS1"     "PDZRN4"       "THEMIS"       "CLSTN2"

Cluster 3 - neuron

count_per_cluster[,c(1,4)]

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

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] "CTC-552D5.1" "CBLN2"       "CDH9"        "GRM1"        "AC011288.2" 
##  [6] "NDST3"       "FSTL4"       "GNAL"        "CRACDL"      "EPHA6"

Cluster 4 - polydendrocyte

count_per_cluster[,c(1,5)]

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

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] "RP4-668E10.4" "VCAN"         "SMOC1"        "NXPH1"        "LHFPL3"      
##  [6] "MEGF11"       "GRIK1"        "TMEM132C"     "PTPRZ1"       "DISC1"

Cluster 5 - neuron

count_per_cluster[,c(1,6)]

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

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] "RP11-320L2.1" "CTC-340A15.2" "CTC-535M15.2" "VIM2P"        "RP11-197K6.1"
##  [6] "GRIN3A"       "TRABD2A"      "LINC01821"    "NRG1"         "TAFA1"

Cluster 6 - neuron - RP high

count_per_cluster[,c(1,7)]

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

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" "RP11-33A14.1"  "COL5A2"        "CTC-552D5.1"  
##  [5] "CUX2"          "RP11-78F17.1"  "LINC02822"     "RP11-197K6.1" 
##  [9] "CNGB1"         "LINC00326"

Cluster 7 - astrocyte

count_per_cluster[,c(1,8)]

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

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] "RP11-6L16.1"   "GLIS3"         "RP11-134O15.3" "RP11-627D16.1"
##  [5] "RFX4"          "ADGRV1"        "AQP4"          "GFAP"         
##  [9] "LINC00299"     "ZFP36L1"

Cluster 8 - neuron

count_per_cluster[,c(1,9)]

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

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] "CTD-2337L2.1" "HS3ST4"       "PCDH11X"      "FRMD3"        "KIAA1217"    
##  [6] "MCTP2"        "AC002331.2"   "RP1-133I17.1" "DPP4"         "RASEF"

Cluster 9 - neuron

count_per_cluster[,c(1,10)]

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

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] "DLX6-AS1"  "PWRN1"     "RGS12"     "SYNPR"     "CALB2"     "BTBD11"   
##  [7] "ADARB2"    "NR2F2-AS1" "LAMA3"     "PROX1-AS1"

Cluster 10 - neuron

count_per_cluster[,c(1,11)]

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

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

Cluster 11 - interneuron

count_per_cluster[,c(1,12)]

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

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"         "GAD2"          "SCN1A-AS1"     "LHX6"         
##  [5] "MEPE"          "RP11-656G20.1" "PTCHD4"        "ARX"          
##  [9] "SAMD5"         "ANK1"

Cluster 12 - interneuron

count_per_cluster[,c(1,13)]

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

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] "FGF13"   "PTCHD4"  "EYA4"    "GAD2"    "KIT"     "TACR1"   "SV2C"   
##  [8] "KIRREL1" "NXPH2"   "NRIP3"

Cluster 13 - neuron

count_per_cluster[,c(1,14)]

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

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] "RP11-20A20.2" "CTC-535M15.2" "CTC-340A15.2" "PLCH1"        "LINC01164"   
##  [6] "RORB"         "RMST"         "TRHDE"        "RP11-197K6.1" "CPNE4"

Cluster 14 - oligodendrocyte

count_per_cluster[,c(1,15)]

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

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] "LRP2"          "LRRC63"        "RP1-223B1.1"   "SLC5A11"      
##  [5] "RP11-314E23.2" "ITGA2"         "LINC01445"     "RP4-630C24.3" 
##  [9] "CARNS1"        "TMEM140"

Cluster 15 - neuron

count_per_cluster[,c(1,16)]

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

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] "COL5A2"        "SYT2"          "AJ006998.2"    "ATP6V1C2"     
##  [5] "SCN1A-AS1"     "RP11-153I24.5" "RP11-335E8.3"  "BCL11B"       
##  [9] "GRM8"          "COL24A1"

Cluster 16 - noise

count_per_cluster[,c(1,17)]

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

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] "HHIP"          "OPALIN"        "BOK"           "GPR37"        
##  [5] "PLPP2"         "MYRF"          "RP11-603J24.9" "MOG"          
##  [9] "MID1IP1"       "PRMT6"

Cluster 17 - neuron

count_per_cluster[,c(1,18)]

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

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] "CXCL14"        "NR2F2"         "DLX6-AS1"      "RELN"         
##  [5] "NR2F2-AS1"     "RP11-12B13.1"  "GAD2"          "RP11-406A20.1"
##  [9] "RGS12"         "GRIK1"

Cluster 18 - endothelial

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"        "SMYD1"         "RP11-47B24.1"  "NR4A2"        
##  [5] "PLA2G4A"       "SMIM32"        "RP11-632B21.2" "OLFML2B"      
##  [9] "STPG2-AS1"     "DCSTAMP"

Cluster 19 - microglia

count_per_cluster[,c(1,19)]

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

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] "ADAM28"   "TBXAS1"   "FYB1"     "LNCAROD"  "DOCK8"    "APBB1IP" 
##  [7] "INPP5D"   "ARHGAP15" "CD74"     "SLCO2B1"

Cluster 20 - neuron - RP high

count_per_cluster[,c(1,20)]

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

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] "RP11-897M7.4"  "RP11-79E3.2"   "RP11-404H1.1"  "CTD-2195H9.1" 
##  [5] "LINC02364"     "AC005150.1"    "RP11-582C12.1" "RORB"         
##  [9] "IL1RAPL2"      "CEMIP"

Cluster 21 - polydendrocyte

count_per_cluster[,c(1,21)]

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

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] "GPNMB"         "RP11-739N10.1" "LINC02552"     "LINC02315"    
##  [5] "RP4-668E10.4"  "COL9A1"        "LINC02058"     "CSPG4"        
##  [9] "LMO1"          "MIR3681HG"

Cluster 22 - neuron

count_per_cluster[,c(1,22)]

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

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] "HTR2C"         "NPSR1-AS1"     "IFNG-AS1"      "FER1L6-AS2"   
##  [5] "NXPH2"         "ENPP7P1"       "LINC02218"     "RP11-354I13.2"
##  [9] "CHRM2"         "COL12A1"

Cluster 23 - astrocyte / oligodendrocyte

count_per_cluster[,c(1,23)]

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

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] "ETNPPL"       "OBI1-AS1"     "LINC00499"    "RP11-268P4.5" "GLIS3"       
##  [6] "PPP1R9A-AS1"  "GJA1"         "MOBP"         "ADGRV1"       "ST18"

Cluster 24 - neuron

count_per_cluster[,c(1,24)]

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

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] "SCUBE3"        "NOG"           "LHX6"          "SP9"          
##  [5] "RP11-238K6.2"  "GPR149"        "RP11-519H15.1" "COL15A1"      
##  [9] "SCN1A-AS1"     "ANK1"

Cluster 25 - mural

count_per_cluster[,c(1,25)]

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

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] "HIGD1B" "HEYL"   "GPER1"  "MYL9"   "GJA4"   "TBX2"   "SIK1B"  "FOXD1" 
##  [9] "SIK1"   "PTGDR2"

Cluster 26 - endothelial

count_per_cluster[,c(1,26)]

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

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

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

cluster26$gene[1:10]

##  [1] "CLDN5"         "CFH"           "LRRC32"        "RP11-326C3.17"
##  [5] "CAVIN2"        "TGM2"          "GNG11"         "IFITM1"       
##  [9] "AC000082.4"    "SRARP"

Assign identities

Individual

dataObject.annotated <- RenameIdents(object = dataObject, 
                               "0" = "oligodendrocyte 1",
                               "1" = "astrocyte 1",
                               "2" = "neuron 1",
                               "3" = "neuron 2",
                               "4" = "polydendrocyte 1",
                               "5" = "neuron 3",
                               "6" = "neuron - RP 1",
                               "7" = "astrocyte 2",
                               "8" = "neuron 4",
                               "9" = "neuron 5",
                               "10" = "neuron 6",
                               "11" = "interneuron 1",
                               "12" = "interneuron 2",
                               "13" = "neuron 7",
                               "14" = "oligodendrocyte 2",
                               "15" = "neuron 8",
                               "16" = "noise",
                               "17" = "neuron 9",
                               "18" = "endothelial",
                               "19" = "microglia",
                               "20" = "neuron - RP 2",
                               "21" = "polydendrocyte 2",
                               "22" = "neuron 10",
                               "23" = "astrocyte 3",
                               "24" = "neuron 11",
                               "25" = "mural", 
                               "26" = "endothelial")
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 neuron 1 neuron 2 polydendrocyte 1
## 1 v4PA_mod0.6x              1518        1275     1163     1064             1050
##   neuron 3 neuron - RP 1 astrocyte 2 neuron 4 neuron 5 neuron 6 interneuron 1
## 1      671           627         580      488      452      450           365
##   interneuron 2 neuron 7 oligodendrocyte 2 neuron 8 noise neuron 9 endothelial
## 1           245      219               211      203   171      154         171
##   microglia neuron - RP 2 polydendrocyte 2 neuron 10 astrocyte 3 neuron 11
## 1       144           129              121       115          67        67
##   mural
## 1    26

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",
"neuron 1" = "neuron",
"neuron 2" = "neuron",
"polydendrocyte 1" = "polydendrocyte",
"neuron 3" = "neuron",
"neuron - RP 1" = "neuron",
"astrocyte 2" = "astrocyte",
"neuron 4" = "neuron",
"neuron 5" = "neuron",
"neuron 6" = "neuron",
"interneuron 1" = "neuron",
"interneuron 2" = "neuron",
"neuron 7" = "neuron",
"oligodendrocyte 2" = "oligodendrocyte",
"neuron 8" = "neuron",
"noise" = "noise",
"neuron 9" = "neuron",
"endothelial" = "endothelial",
"microglia" = "microglia",
"neuron - RP 2" = "neuron",
"polydendrocyte 2" = "polydendrocyte",
"neuron 10" = "neuron",
"astrocyte 3" = "astrocyte",
"neuron 11" = "neuron",
"mural" =   "mural", 
"endothelial" = "endothelial")
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 neuron polydendrocyte noise
## 1 v4PA_mod0.6x            1729      1922   6412           1171   171
##   endothelial microglia mural
## 1         171       144    26

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"))