Read in processed filtered data. Remove doublets and reprocess. Check resolutions and save object.

Libraris, paths, colors

Read in object

# read object
dataObject <- readRDS(file = paste0("../rObjects/", treatment, "_unannotated_doublets_removed.rds"))
DefaultAssay(dataObject) <- "RNA"
dataObject <- NormalizeData(dataObject)
dataObject <- FindVariableFeatures(dataObject)
dataObject <- ScaleData(dataObject)
dataObject <- JoinLayers(dataObject)

# inspect
dataObject
## An object of class Seurat 
## 17728 features across 2498 samples within 2 assays 
## Active assay: RNA (8864 features, 2000 variable features)
##  3 layers present: counts, data, scale.data
##  1 other assay present: SCT
##  2 dimensional reductions calculated: pca, umap

Unannotated

UMAP

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

### 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

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
## 1 SeuratProject 578 521 373 348 286 187 145 60
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_colors) +
  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

Violins

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

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

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

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

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

Markers per cluster

markers <- SeuratWrappers::RunPrestoAll(
  object = dataObject,
  assay = "RNA",
  slot = "counts",
  only.pos = FALSE
)
write.table(markers, 
            paste0("../results/markers/", treatment, "_markers.tsv"),
            quote = FALSE,
            row.names = FALSE)
saveRDS(markers, paste0("../rObjects/", treatment, "_wilcox_markers.rds"))

# rearrange to order by cluster & filter to only include log2FC > 1 & FDR < 0.05
 all.markers.strict <- markers %>%
   group_by(cluster) %>%
   dplyr::filter(avg_log2FC > 1 & p_val_adj < 0.05)
saveRDS(all.markers.strict, paste0("../rObjects/", treatment,"_wilcox_markers_log2FC1_q0.01.rds"))

Get markers for each cluster

# unique clusters variable
unique_clusters <- unique(all.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 <- all.markers.strict[all.markers.strict$cluster == i, ]
  assign(cluster_name, cluster_data)
  cluster_list[[cluster_name]] <- cluster_data
}

Cluster Annotation

Cluster 0 - neuron

# Number of cells per condition
count_per_cluster[,c(1,2)]
##      orig.ident   0
## 1 SeuratProject 578
# 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] "VWA5B2"             "KCNJ6"              "ENSSSCG00000061658"
##  [4] "MYBPC1"             "ITGA2B"             "CPLX2"             
##  [7] NA                   NA                   NA                  
## [10] NA

Cluster 1 - neuron

count_per_cluster[,c(1,3)]
##      orig.ident   1
## 1 SeuratProject 521
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] "ENSSSCG00000048856" "CACNB2"             "GALNT17"           
##  [4] "SYT17"              "CA10"               "NXN"               
##  [7] "ENSSSCG00000058255" "TAFA1"              "CUX2"              
## [10] "RFX3"

Cluster 2 - noise/neuron/polydendrocyte?

count_per_cluster[,c(1,4)]
##      orig.ident   2
## 1 SeuratProject 373
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] "ENSSSCG00000055221" "ENSSSCG00000061137" "C10orf90"          
##  [4] "QKI"                "POLR2F"             "NCKAP5"            
##  [7] "ZNF536"             "RNF220"             "ENSSSCG00000058036"
## [10] "MBP"

Cluster 3 - neuron

count_per_cluster[,c(1,4)]
##      orig.ident   2
## 1 SeuratProject 373
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] "GPC6"   "HS6ST3" "NCAM2"  "PDZRN3" "ZMAT4"  "KLHL29" "RASSF3" "SNAP25"
##  [9] "ADGRV1" "OSBP2"

Cluster 4 - neuron

count_per_cluster[,c(1,5)]
##      orig.ident   3
## 1 SeuratProject 348
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] "KIAA1217" "HS3ST4"   "DPP10"    "GRIK3"    "MCTP1"    "ZFPM2"   
##  [7] "THSD7B"   "CLSTN2"   "FOXP2"    "SULF1"

Cluster 5 - neuron /interneuron

count_per_cluster[,c(1,6)]
##      orig.ident   4
## 1 SeuratProject 286
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] "GNG4"     "CADPS2"   "PCDH9"    "ZNF804B"  "MGAT5"    "ADAMTSL3"
##  [7] "CACNB4"   "KCTD16"   "GNB4"     "SGCZ"

Cluster 6 - interneuron

count_per_cluster[,c(1,7)]
##      orig.ident   5
## 1 SeuratProject 187
DimPlot(object = subset(dataObject, seurat_clusters == "6"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = color.panel[7])

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

cluster6$gene[1:10]
##  [1] "ERBB4"   "NXPH1"   "GALNTL6" "KCNIP1"  "PTPRM"   "UBASH3B" "ANK1"   
##  [8] "RBMS3"   "NPAS3"   "KCNC2"

Cluster 7 - neuron

count_per_cluster[,c(1,8)]
##      orig.ident   6
## 1 SeuratProject 145
DimPlot(object = subset(dataObject, seurat_clusters == "7"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = color.panel[8])

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

cluster6$gene[1:10]
##  [1] "ERBB4"   "NXPH1"   "GALNTL6" "KCNIP1"  "PTPRM"   "UBASH3B" "ANK1"   
##  [8] "RBMS3"   "NPAS3"   "KCNC2"

Assign identities

Individual

dataObject.annotated <- RenameIdents(object = dataObject, 
                               "0" = "neuron",
                               "1" = "neuron",
                               "2" = "noise",
                               "3" = "neuron",
                               "4" = "neuron",
                               "5" = "neuron",
                               "6" = "interneuron",
                               "7" = "neuron")
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

Remove noise

# filter
dataObject.no_noise <- subset(dataObject.annotated, subset = 
    (individual_clusters != "noise") 
)

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

Reprocess post-noise removal

# save singlets as the new dataObject
dataObject <- dataObject.no_noise

# transform
dataObject <- SCTransform(dataObject, verbose = FALSE)

# run PCA on the merged object
dataObject <- RunPCA(object = dataObject)
Idents(dataObject) <- "sample"

# Determine the K-nearest neighbor graph
dataObject <- FindNeighbors(object = dataObject, 
                                 assay = "SCT", 
                                 reduction = "pca",
                                 dims = 1:15)
# Run UMAP
dataObject <- RunUMAP(dataObject,
                           dims = 1:15,
                           reduction = "pca",
                           n.components = 3) 

# Determine the clusters for various resolutions
dataObject <- FindClusters(object = dataObject,
                                 algorithm = 1, # 1= Louvain
                                 resolution = seq(0.1,1,by=0.1))
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 2125
## Number of edges: 84926
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9463
## Number of communities: 4
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 2125
## Number of edges: 84926
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9022
## Number of communities: 4
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 2125
## Number of edges: 84926
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8678
## Number of communities: 5
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 2125
## Number of edges: 84926
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8438
## Number of communities: 6
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 2125
## Number of edges: 84926
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8201
## Number of communities: 6
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 2125
## Number of edges: 84926
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8014
## Number of communities: 7
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 2125
## Number of edges: 84926
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7836
## Number of communities: 7
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 2125
## Number of edges: 84926
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7658
## Number of communities: 7
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 2125
## Number of edges: 84926
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7480
## Number of communities: 7
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 2125
## Number of edges: 84926
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7303
## Number of communities: 7
## Elapsed time: 0 seconds

Explore resolutions

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

# 0.4
umap0.2 <- DimPlot(dataObject,
        group.by = "SCT_snn_res.0.2",
        label = TRUE)
umap0.2

# 0.4
umap0.4 <- DimPlot(dataObject,
        group.by = "SCT_snn_res.0.4",
        label = TRUE)
umap0.4

Set resolution

dataObject$seurat_clusters <- dataObject$SCT_snn_res.0.2
Idents(dataObject) <- dataObject$seurat_clusters
DefaultAssay(dataObject) <- "RNA"
dataObject <- NormalizeData(dataObject)
dataObject <- FindVariableFeatures(dataObject)
dataObject <- ScaleData(dataObject)
dataObject <- JoinLayers(dataObject)

Violins

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

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

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

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

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
## 1 SeuratProject 1509 286 187 143
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

Markers per cluster

markers <- SeuratWrappers::RunPrestoAll(
  object = dataObject,
  assay = "RNA",
  slot = "counts",
  only.pos = FALSE
)
write.table(markers, 
            paste0("../results/markers/", treatment, "_markers_noise_clusters_removed.tsv"),
            quote = FALSE,
            row.names = FALSE)
saveRDS(markers, paste0("../rObjects/", treatment, "_wilcox_markers_noise_clusters_removed.rds"))

# rearrange to order by cluster & filter to only include log2FC > 1 & FDR < 0.05
 all.markers.strict <- markers %>%
   group_by(cluster) %>%
   dplyr::filter(avg_log2FC > 1 & p_val_adj < 0.05)
saveRDS(all.markers.strict, paste0("../rObjects/", treatment,"_wilcox_markers_noise_clusters_removed_log2FC1_q0.01.rds"))

Get markers for each cluster

# unique clusters variable
unique_clusters <- unique(all.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 <- all.markers.strict[all.markers.strict$cluster == i, ]
  assign(cluster_name, cluster_data)
  cluster_list[[cluster_name]] <- cluster_data
}

Cluster Annotation

Cluster 0 - neuron

# Number of cells per condition
count_per_cluster[,c(1,2)]
##      orig.ident    0
## 1 SeuratProject 1509
# 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] "HS6ST3"  "TAFA1"   "GRM7"    "SYT17"   "GPC6"    "CA10"    "SLIT3"  
##  [8] "GALNT17" "TTLL11"  "KCNB2"

Cluster 1 - neuron

count_per_cluster[,c(1,3)]
##      orig.ident   1
## 1 SeuratProject 286
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] "KIAA1217" "HS3ST4"   "GRIK3"    "DPP10"    "MCTP1"    "ZFPM2"   
##  [7] "THSD7B"   "CLSTN2"   "FOXP2"    "SULF1"

Cluster 2 - neuron

count_per_cluster[,c(1,4)]
##      orig.ident   2
## 1 SeuratProject 187
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] "GNG4"     "PCDH9"    "ZNF804B"  "CADPS2"   "MGAT5"    "ADAMTSL3"
##  [7] "GNB4"     "KCTD16"   "CACNB4"   "SH3GL3"

Cluster 3 - interneuron

count_per_cluster[,c(1,4)]
##      orig.ident   2
## 1 SeuratProject 187
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] "ERBB4"    "NXPH1"    "GALNTL6"  "KCNIP1"   "PTPRM"    "ANK1"    
##  [7] "NPAS3"    "PTCHD4"   "ZNF536"   "TMEM132C"

Assign identities

Individual

dataObject.annotated <- RenameIdents(object = dataObject, 
                               "0" = "neuron",
                               "1" = "neuron",
                               "2" = "neuron",
                               "3" = "interneuron")
dataObject.annotated$annotated_clusters <- factor(Idents(dataObject.annotated))

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

## png 
##   2

Save

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