v4 Fluent annotation

Pipseq chemistry kit comparison Sample NPID control female NA05-055

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_Fluent")
sample_color.panel <- c("pink")
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
## 1   v4Fluent 2671 2327 2001 1780 1746 1472 1321 1125 1025 863 731 695 694 506
##    14  15  16  17  18  19  20 21
## 1 497 483 399 365 305 226 157 66

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   v4Fluent 2671

# 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"          "CNP"           "PLP1"          "CNDP1"        
##  [5] "CLMN"          "RNF220"        "TMEM144"       "RP11-654K19.6"
##  [9] "KCNH8"         "SUN2"

Cluster 1 - oligodendrocyte

count_per_cluster[,c(1,3)]

##   orig.ident    1
## 1   v4Fluent 2327

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] "ST18"          "RNF220"        "CLMN"          "CTNNA3"       
##  [5] "KCNH8"         "RP11-654K19.6" "FRMD4B"        "TMEM144"      
##  [9] "SLC7A14-AS1"   "MOBP"

Cluster 2 - oligodendrocyte

count_per_cluster[,c(1,4)]

##   orig.ident    2
## 1   v4Fluent 2001

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

Cluster 3 - astrocyte

count_per_cluster[,c(1,4)]

##   orig.ident    2
## 1   v4Fluent 2001

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] "RP11-6L16.1" "RFX4"        "CD44"        "FAM189A2"    "GFAP"       
##  [6] "GLIS3"       "AQP4"        "LINC02649"   "ZFP36L1"     "ITPR2"

Cluster 4 - oligodendrocyte

count_per_cluster[,c(1,5)]

##   orig.ident    3
## 1   v4Fluent 1780

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] "SLC5A11"     "CNP"         "TF"          "SPP1"        "ENPP2"      
##  [6] "MOBP"        "RP1-223B1.1" "CARNS1"      "CLDND1"      "SYNJ2"

Cluster 5 - astrocyte

count_per_cluster[,c(1,6)]

##   orig.ident    4
## 1   v4Fluent 1746

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] "GNA14"    "IRAG1"    "CABLES1"  "OBI1-AS1" "BMPR1B"   "SLC14A1" 
##  [7] "ADGRV1"   "SLC1A2"   "RFX4"     "FAM189A2"

Cluster 6 - noise

count_per_cluster[,c(1,7)]

##   orig.ident    5
## 1   v4Fluent 1472

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] "GADD45B"        "CNDP1"          "PRUNE2"         "GUSBP2"        
##  [5] "GLUL"           "FGFR1"          "PLA2G4C"        "ELL2"          
##  [9] "IQGAP1"         "RP11-1023L17.1"

Cluster 7 - neuron

count_per_cluster[,c(1,8)]

##   orig.ident    6
## 1   v4Fluent 1321

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] "CDH9"       "CBLN2"      "MIR137HG"   "NDST3"      "CRACDL"    
##  [6] "LDB2"       "EPHA4"      "GRM1"       "AC011288.2" "NELL2"

Cluster 8 - polydendrocyte

count_per_cluster[,c(1,9)]

##   orig.ident    7
## 1   v4Fluent 1125

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] "MEGF11"       "LHFPL3"       "VCAN"         "TNR"          "LUZP2"       
##  [6] "NXPH1"        "PTPRZ1"       "BRINP3"       "RP4-668E10.4" "PCDH15"

Cluster 9 - neuron

count_per_cluster[,c(1,10)]

##   orig.ident    8
## 1   v4Fluent 1025

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] "MIR137HG"      "RP11-213B3.1"  "ST6GALNAC5"    "RP11-673E1.1" 
##  [5] "CBLN2"         "LDB2"          "AC011288.2"    "RP11-147G16.1"
##  [9] "EPHA4"         "CRACDL"

Cluster 10 - microglia

count_per_cluster[,c(1,11)]

##   orig.ident   9
## 1   v4Fluent 863

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] "ADAM28"   "SLC11A1"  "DOCK8"    "LRMDA"    "ARHGAP24" "PTPRC"   
##  [7] "APBB1IP"  "ARHGAP15" "TBXAS1"   "FYB1"

Cluster 11 - neuron

count_per_cluster[,c(1,12)]

##   orig.ident  10
## 1   v4Fluent 731

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] "SYNPR"   "RGS12"   "GALNTL6" "KCNT2"   "ADARB2"  "DSCAM"   "CHRM3"  
##  [8] "MYT1L"   "GRIP1"   "ASIC2"

Cluster 12 - neuron

count_per_cluster[,c(1,13)]

##   orig.ident  11
## 1   v4Fluent 695

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] "NRG1"     "CLSTN2"   "RALYL"    "NELL2"    "MIR137HG" "KCNH7"   
##  [7] "LDB2"     "HECW1"    "ARPP21"   "LRFN5"

Cluster 13 - neuron

count_per_cluster[,c(1,14)]

##   orig.ident  12
## 1   v4Fluent 694

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] "KCNIP4"       "RP11-909M7.3" "MEG3"         "RALYL"        "BCYRN1"      
##  [6] "AC011288.2"   "NELL2"        "CHN1"         "RP11-213B3.1" "NDST3"

Cluster 14 - neuron

count_per_cluster[,c(1,15)]

##   orig.ident  13
## 1   v4Fluent 506

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] "NXPH1"         "PTCHD4"        "ZNF385D"       "ZNF804A"      
##  [5] "SCN1A-AS1"     "GRIP1"         "GABRG3"        "RP11-123O10.3"
##  [9] "KCNC2"         "MYO16"

Cluster 15 - neuron

count_per_cluster[,c(1,16)]

##   orig.ident  14
## 1   v4Fluent 497

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] "CBLN2"        "RP11-909M7.3" "RALYL"        "AC011288.2"   "MEG3"        
##  [6] "KCNIP4"       "KCNQ5"        "TESPA1"       "STXBP5L"      "CDH18"

Cluster 16 - neuron

count_per_cluster[,c(1,17)]

##   orig.ident  15
## 1   v4Fluent 483

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] "CXCL14"    "FSTL5"     "GRIK1"     "RGS12"     "GRIP1"     "NR2F2-AS1"
##  [7] "ZNF385D"   "ADARB2"    "NR2F2"     "C8orf34"

Cluster 17 - neuron

count_per_cluster[,c(1,18)]

##   orig.ident  16
## 1   v4Fluent 399

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] "GRIK1"     "LINC01322" "NXPH1"     "KIAA1217"  "SYNPR"     "KIF26B"   
##  [7] "GRIK3"     "ARX"       "MAFB"      "KLF5"

Cluster 18 - neuron

count_per_cluster[,c(1,1)]

##   orig.ident orig.ident.1
## 1   v4Fluent     v4Fluent

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] "FGF13"         "RP11-123O10.3" "PTCHD4"        "GRIK1"        
##  [5] "MYO16"         "NXPH1"         "EYA4"          "FSTL5"        
##  [9] "GRIP1"         "SGCZ"

Cluster 19 - noise

count_per_cluster[,c(1,19)]

##   orig.ident  17
## 1   v4Fluent 365

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] "SVIL"          "ARHGAP29"      "RP11-649A16.1" "CFH"          
##  [5] "EBF1"          "RP11-326C3.17" "LEF1"          "C7"           
##  [9] "GGT5"          "COL1A2"

Cluster 20 - neuron

count_per_cluster[,c(1,20)]

##   orig.ident  18
## 1   v4Fluent 305

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] "HS3ST4"       "CTD-2337L2.1" "RAB3B"        "EGFEM1P"      "CDH6"        
##  [6] "SULF1"        "DIRAS2"       "ZNF385B"      "FEZF2"        "SEMA3E"

Cluster 21 - mural

count_per_cluster[,c(1,21)]

##   orig.ident  19
## 1   v4Fluent 226

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] "MYL9"   "ACTA2"  "TAGLN"  "TPM2"   "DES"    "TGM2"   "COL8A1" "ACTG2" 
##  [9] "CARMN"  "ITGA5"

Assign identities

Individual

dataObject.annotated <- RenameIdents(object = dataObject, 
                               "0" = "oligodendrocyte 1",
                               "1" = "oligodendrocyte 2",
                               "2" = "oligodendrocyte 3",
                               "3" = "astrocyte 1",
                               "4" = "oligodendrocyte 4",
                               "5" = "astrocyte 2",
                               "6" = "noise 1",
                               "7" = "neuron 1",
                               "8" = "polydendrocyte",
                               "9" = "neuron 2",
                               "10" = "microglia",
                               "11" = "neuron 3",
                               "12" = "neuron 4",
                               "13" = "neuron 5",
                               "14" = "neuron 6",
                               "15" = "neuron 7",
                               "16" = "neuron 8",
                               "17" = "neuron 9",
                               "18" = "neuron 10",
                               "19" = "noise 2",
                               "20" = "neuron 11",
                               "21" = "mural")
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 oligodendrocyte 2 oligodendrocyte 3 astrocyte 1
## 1   v4Fluent              2671              2327              2001        1780
##   oligodendrocyte 4 astrocyte 2 noise 1 neuron 1 polydendrocyte neuron 2
## 1              1746        1472    1321     1125           1025      863
##   microglia neuron 3 neuron 4 neuron 5 neuron 6 neuron 7 neuron 8 neuron 9
## 1       731      695      694      506      497      483      399      365
##   neuron 10 noise 2 neuron 11 mural
## 1       305     226       157    66

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, 
                               "0" = "oligodendrocyte",
                               "1" = "oligodendrocyte",
                               "2" = "oligodendrocyte",
                               "3" = "astrocyte",
                               "4" = "oligodendrocyte",
                               "5" = "astrocyte",
                               "6" = "noise",
                               "7" = "neuron",
                               "8" = "polydendrocyte",
                               "9" = "neuron",
                               "10" = "microglia",
                               "11" = "neuron",
                               "12" = "neuron",
                               "13" = "neuron",
                               "14" = "neuron",
                               "15" = "neuron",
                               "16" = "neuron",
                               "17" = "neuron",
                               "18" = "neuron",
                               "19" = "noise",
                               "20" = "neuron",
                               "21" = "mural")
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 noise neuron polydendrocyte microglia
## 1   v4Fluent            8745      3252  1547   6089           1025       731
##   mural
## 1    66

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