Single nucleus

Annotation of data from “Human and mouse single-nucleus transcriptomics reveal TREM2-dependent and TREM2-independent cellular responses in Alzheimer’s disease” Zhou et al. 2020
https://www.nature.com/articles/s41591-019-0695-9#data-availability

Set working directory

Load libraries

save function

Set variables and thresholds

pathToRef <- "/research/labs/neurology/fryer/projects/references/mouse/refdata-gex-mm10-2020-A"
pathToTestType <- "7months/"
treatment <- "mouse"
tissue <- "Brain"
seqType <- "nucleus"
tool <- "cellranger"

# replicate groups
replicate <- "all"
replicate1 <- "replicate1"
replicate2 <- "replicate2"
replicate3 <- "replicate3"

# genotype
WT <- "WT"
Trem2_KO <- "Trem2_KO"
WT_5XFAD <- "WT_5XFAD"
Trem2_KO_5XFAD <- "Trem2_KO_5XFAD"
sample_order <- c(
      "WT_1",
      "WT_2",
      "WT_3",
      "Trem2_KO_1",
      "Trem2_KO_2",
      "Trem2_KO_3",
      "WT_5XFAD_1",
      "WT_5XFAD_2",
      "WT_5XFAD_3",
      "Trem2_KO_5XFAD_1",
      "Trem2_KO_5XFAD_2",
      "Trem2_KO_5XFAD_3"
)

Read data

# read object
data.unannotated <- readRDS(paste0("../../rObjects/", pathToTestType, 
                                  tolower(tissue), "_unannotated.rds"))
# Run TSNE to have for shiny app 
data.unannotated <- RunTSNE(data.unannotated)
data.unannotated$seurat_clusters <- data.unannotated$SCT_snn_res.0.1
Idents(data.unannotated) <- "seurat_clusters"
DefaultAssay(data.unannotated) <- "RNA"
data.unannotated <- NormalizeData(data.unannotated)

data.unannotated
## An object of class Seurat 
## 30037 features across 80588 samples within 2 assays 
## Active assay: RNA (15019 features, 0 variable features)
##  1 other assay present: SCT
##  4 dimensional reductions calculated: pca, harmony, umap, tsne

colors

# get colors
# dittoColors()
cluster_colors <- c("#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7", "#666666", "#AD7700", "#1C91D4", "#007756", "#D5C711", "#005685")
dittoDimPlot(object = data.unannotated,
             var = "seurat_clusters",
             reduction.use = "umap",
             do.label = TRUE,
             labels.highlight = TRUE)

dittoDimPlot(object = data.unannotated,
             var = "seurat_clusters",
             reduction.use = "tsne",
             do.label = TRUE,
             labels.highlight = TRUE)

Cluster tree

  • Cluster trees are helpful in deciding what clusters to merge.
Idents(data.unannotated) <- "seurat_clusters"
data.unannotated <- BuildClusterTree(object = data.unannotated,
                                     dims = 1:15, # min.pc in processing script
                                     reorder = FALSE,
                                     reorder.numeric = FALSE)
tree <- data.unannotated@tools$BuildClusterTree
tree$tip.label <- paste0(tree$tip.label)

ggtree::ggtree(tree, aes(x, y)) +
  scale_y_reverse() +
  ggtree::geom_tree() +
  ggtree::theme_tree() +
  ggtree::geom_tiplab(offset = 1) +
  ggtree::geom_tippoint(shape = 16, size = 5) +
  ggtree::geom_tippoint(color = cluster_colors[1:length(tree$tip.label)], shape = 16, size = 5) +
  coord_cartesian(clip = 'off') +
  theme(plot.margin = unit(c(0,2.5,0,0), 'cm'))

Potential Markers

Highly variable PCs

# Printing out the most variable genes driving PCs
print(x = data.unannotated[["harmony"]], 
      dims = 1:10, 
      nfeatures = 10)
## harmony_ 1 
## Positive:  Plp1, Mbp, Neat1, Ptgds, Mag, Cldn11, Mog, Mal, Apod, Gatm 
## Negative:  Celf2, Syt1, Gria2, Arpp21, Pde10a, Kcnip4, Kalrn, Phactr1, Csmd1, Ahi1 
## harmony_ 2 
## Positive:  Atp1a2, Slc1a3, Slc1a2, Cst3, Apoe, Plpp3, Gpr37l1, Bcan, Ptprz1, Ntsr2 
## Negative:  Plp1, Mbp, Ptprd, Mag, Cldn11, Mog, Mal, Enpp2, Pex5l, Ptgds 
## harmony_ 3 
## Positive:  Slc1a2, Atp1a2, Slc1a3, Plpp3, Gpr37l1, Bcan, Ntsr2, Ptprz1, Sparcl1, Clu 
## Negative:  Hexb, Csf1r, C1qc, P2ry12, C1qa, Laptm5, C1qb, Unc93b1, Arhgap45, Cx3cr1 
## harmony_ 4 
## Positive:  Ptprd, Tcf4, Mef2c, Dlgap1, Miat, Gm28928, Car10, Opcml, Syt1, Slc17a7 
## Negative:  Pde10a, Rgs9, Penk, Phactr1, Gpr88, Adcy5, Rarb, Gprin3, Adora2a, Meis2 
## harmony_ 5 
## Positive:  Pdgfra, Ptprz1, Lhfpl3, Vcan, Cspg4, Adarb2, Nxph1, Serpine2, Tnr, C1ql1 
## Negative:  Slc1a2, Ptprd, Ntsr2, Glul, Kalrn, Slc1a3, R3hdm1, S1pr1, Arpp21, Gja1 
## harmony_ 6 
## Positive:  Pdgfra, Lhfpl3, Cspg4, Vcan, Ptprz1, Ptprd, Tnr, Matn4, Arpp21, Gpr17 
## Negative:  Erbb4, Dlx6os1, Adarb2, Gad1, Slc6a1, Atp1b1, Kcnip1, Slc32a1, Galntl6, Nxph1 
## harmony_ 7 
## Positive:  Ptgds, Slc6a20a, mt-Co1, Vtn, mt-Co3, mt-Cytb, mt-Atp8, mt-Co2, mt-Nd4l, Aebp1 
## Negative:  Ptprz1, Bcan, Luzp2, Gpr37l1, Erbb4, Slc6a1, Adarb2, Nxph1, Ntsr2, Gm3764 
## harmony_ 8 
## Positive:  A830036E02Rik, Mef2c, Lingo2, Fam19a1, BC030499, Kcnip4, Ptprd, Rorb, 2900055J20Rik, Cntn5 
## Negative:  Thsd7b, Zfpm2, Sdk2, Gm28928, Foxp2, mt-Co3, mt-Co1, Slc1a2, mt-Atp8, mt-Cytb 
## harmony_ 9 
## Positive:  Sparcl1, mt-Co1, Syt1, mt-Co3, mt-Atp8, Chgb, mt-Co2, Pde10a, mt-Cytb, Atp2b2 
## Negative:  Nnat, C130073E24Rik, Hap1, Cpne7, AC124490.1, Camk2d, Baiap3, Ly6h, Celf4, Nrp2 
## harmony_ 10 
## Positive:  Cemip, Nxph1, Coro6, Moxd1, Etl4, Kcnc1, Sox6, Me3, Kcnc2, Cntnap4 
## Negative:  Adarb2, Vip, Cnr1, mt-Co1, mt-Co3, mt-Atp8, mt-Co2, mt-Cytb, mt-Nd4l, Kit

Top 100 variable genes

# print top variable genes
DefaultAssay(data.unannotated) <- "SCT"
VariableFeatures(data.unannotated)[1:100]
##   [1] "Plp1"          "Ptgds"         "Atp1a2"        "Mbp"          
##   [5] "Vip"           "Adarb2"        "Penk"          "Slc1a3"       
##   [9] "Sst"           "Apoe"          "Slc1a2"        "Neat1"        
##  [13] "Vtn"           "Hexb"          "Cst3"          "Ptprz1"       
##  [17] "Npy"           "Slc6a20a"      "Slc7a11"       "Apod"         
##  [21] "Pdgfra"        "Nnat"          "Cldn11"        "Mal"          
##  [25] "Mag"           "Plpp3"         "Pde10a"        "Gm42418"      
##  [29] "Slco1a4"       "Csf1r"         "Mog"           "P2ry12"       
##  [33] "C1qc"          "Tac2"          "Ntsr2"         "Gpr37l1"      
##  [37] "Cnr1"          "Bcan"          "Gjc3"          "Ndnf"         
##  [41] "Flt1"          "Aspa"          "Gm26917"       "Clu"          
##  [45] "C1qa"          "Slc47a1"       "Ctgf"          "Pltp"         
##  [49] "Gatm"          "Bcas1"         "9630013A20Rik" "Atp13a5"      
##  [53] "Car2"          "Enpp2"         "Gpr17"         "C1qb"         
##  [57] "Cx3cr1"        "Nr4a2"         "Cyr61"         "Sparcl1"      
##  [61] "Crh"           "Laptm5"        "Cldn5"         "Glul"         
##  [65] "Dcn"           "Mfge8"         "Fn1"           "Ranbp3l"      
##  [69] "Ly86"          "Vcan"          "Gja1"          "Siglech"      
##  [73] "Slc6a13"       "Mrc1"          "F13a1"         "Fcrls"        
##  [77] "Ctss"          "Arhgap45"      "Selplg"        "Dlx6os1"      
##  [81] "Luzp2"         "Mobp"          "Ugt8a"         "Slc7a10"      
##  [85] "Stab1"         "Ly6c1"         "Rgs5"          "C030029H02Rik"
##  [89] "Unc93b1"       "Hapln2"        "Hpgds"         "Cspg4"        
##  [93] "7630403G23Rik" "Cyth4"         "Adgrl4"        "Slc13a4"      
##  [97] "Ptprd"         "S1pr1"         "Inpp5d"        "Ltc4s"
DefaultAssay(data.unannotated) <- "RNA"

Astrocytes

Idents(data.unannotated) <- "seurat_clusters"
VlnPlot(data.unannotated,
        features = "Gfap",
        cols = cluster_colors,
        split.by = "seurat_clusters",
        group.by = "seurat_clusters")
## The default behaviour of split.by has changed.
## Separate violin plots are now plotted side-by-side.
## To restore the old behaviour of a single split violin,
## set split.plot = TRUE.
##       
## This message will be shown once per session.

VlnPlot(data.unannotated,
        features = "Aldh1l1",
        cols = cluster_colors,
        split.by = "seurat_clusters",
        group.by = "seurat_clusters")

VlnPlot(data.unannotated,
        features = "Slc1a3",
        cols = cluster_colors,
        split.by = "seurat_clusters",
        group.by = "seurat_clusters")

VlnPlot(data.unannotated,
        features = "Clu",
        cols = cluster_colors,
        split.by = "seurat_clusters",
        group.by = "seurat_clusters")

Endothelial cells

VlnPlot(data.unannotated,
        features = "Flt1",
        cols = cluster_colors)

VlnPlot(data.unannotated,
        features = "Ly6c1",
        cols = cluster_colors)

VlnPlot(data.unannotated,
        features = "Ocln",
        cols = cluster_colors)

Microglia / Macrophages

VlnPlot(data.unannotated,
        features = "P2ry12",
        cols = cluster_colors)

VlnPlot(data.unannotated,
        features = "Tmem119",
        cols = cluster_colors)

VlnPlot(data.unannotated,
        features = "Itgam", # aka Cd11b
        cols = cluster_colors)

VlnPlot(data.unannotated,
        features = "Ptprc", # aka Cd45
        cols = cluster_colors)

VlnPlot(data.unannotated,
        features = "Mrc1",
        cols = cluster_colors)

VlnPlot(data.unannotated,
        features = "Ms4a7",
        cols = cluster_colors)

VlnPlot(data.unannotated,
        features = "Pf4",
        cols = cluster_colors)

VlnPlot(data.unannotated,
        features = "Cd14",
        cols = cluster_colors)

VlnPlot(data.unannotated,
        features = "Cd68",
        cols = cluster_colors)

VlnPlot(data.unannotated,
        features = "Ccr5",
        cols = cluster_colors)

Neurons

VlnPlot(data.unannotated,
        features = "Gad1",
        cols = cluster_colors)

VlnPlot(data.unannotated,
        features = "Gad2",
        cols = cluster_colors)

VlnPlot(data.unannotated,
        features = "Rbfox3",
        cols = cluster_colors)

VlnPlot(data.unannotated,
        features = "Nrgn",
        cols = cluster_colors)

VlnPlot(data.unannotated,
        features = "Slc17a7",
        cols = cluster_colors)

Oligodendroyctes

VlnPlot(data.unannotated,
        features = "Mbp",
        cols = cluster_colors)

VlnPlot(data.unannotated,
        features = "Mog",
        cols = cluster_colors)

VlnPlot(data.unannotated,
        features = "Plp1",
        cols = cluster_colors)

VlnPlot(data.unannotated,
        features = "Gpr17",
        cols = cluster_colors)

VlnPlot(data.unannotated,
        features = "Olig1",
        cols = cluster_colors)

VlnPlot(data.unannotated,
        features = "Olig2",
        cols = cluster_colors)

Fibroblasts

VlnPlot(data.unannotated,
        features = "Col1a1",
        cols = cluster_colors)

VlnPlot(data.unannotated,
        features = "Col1a2",
        cols = cluster_colors)

Choroid contamination

VlnPlot(data.unannotated,
        features = "Ttr",
        cols = cluster_colors)

Detect markers

# Find markers for each cluster
# Compares single cluster vs all other clusters
# Default is logfc.threshold = 0.25, min.pct = 0.5
Idents(data.unannotated) <- data.unannotated$SCT_snn_res.0.1
all.markers <- FindAllMarkers(object = data.unannotated,
                              assay = "RNA",
                              test.use = "MAST",
                              verbose = TRUE,
                              return.thresh = 0.05,
                              densify = TRUE)

# add column
all.markers$delta_pct <- abs(all.markers$pct.1 - all.markers$pct.2)

# rename columns and rows
rownames(all.markers) <- 1:nrow(all.markers)
all.markers <- all.markers[,c(6,7,1,5,2:4,8)]
colnames(all.markers)[c(6,7)] <- c("pct_1","pct_2")

# save
saveRDS(all.markers, paste0("../../rObjects/", pathToTestType, "data_all_markers.rds"))
# more stringent filtering
all.markers <- all.markers[all.markers$p_val_adj < 0.01,]

# compare 
table(all.markers$cluster)
## 
##    0    1    2    3    4    5    6    7    8    9   10   11   12 
## 2027 1617 3587 1676 3414 2012 4043 2603 1745 3350  573 2527  910
# subset
cluster0 <- all.markers[all.markers$cluster == 0,]
cluster1 <- all.markers[all.markers$cluster == 1,]
cluster2 <- all.markers[all.markers$cluster == 2,]
cluster3 <- all.markers[all.markers$cluster == 3,]
cluster4 <- all.markers[all.markers$cluster == 4,]
cluster5 <- all.markers[all.markers$cluster == 5,]
cluster6 <- all.markers[all.markers$cluster == 6,]
cluster7 <- all.markers[all.markers$cluster == 7,]
cluster8 <- all.markers[all.markers$cluster == 8,]
cluster9 <- all.markers[all.markers$cluster == 9,]
cluster10 <- all.markers[all.markers$cluster == 10,]
cluster11 <- all.markers[all.markers$cluster == 11,]
cluster12 <- all.markers[all.markers$cluster == 12,]

# save
write.table(all.markers,
            paste0("../../results/", pathToTestType,
                   "markers/putative_cluster_markers.tsv"), 
            sep = "\t", quote = FALSE, row.names = FALSE)

Cluster Annotations

Cluster 0 & 1: Pan-neuronal

# Top 40 markers based on p_val_adj and then avg_log2FC
VlnPlot(data.unannotated,
        features = cluster0$gene[1:20],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "SCT_snn_res.0.1")

VlnPlot(data.unannotated,
        features = cluster0$gene[21:40],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "SCT_snn_res.0.1")

VlnPlot(data.unannotated,
        features = cluster1$gene[1:20],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "SCT_snn_res.0.1")

VlnPlot(data.unannotated,
        features = cluster1$gene[21:40],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "SCT_snn_res.0.1")

VlnPlot(data.unannotated,
        features = "Grin1",
        cols = cluster_colors,
        split.by = "SCT_snn_res.0.1")

VlnPlot(data.unannotated,
        features = "Rbfox3",
        cols = cluster_colors,
        split.by = "SCT_snn_res.0.1")

u1 <- DimPlot(data.unannotated,
              cols = cluster_colors,
              label = TRUE) + theme(legend.position = "none")
u2 <- DimPlot(data.unannotated,
              dims = c(2,3),
              label = TRUE,
              cols = cluster_colors) + theme(legend.position = "none")
f1 <- FeaturePlot(data.unannotated, 
                  features = "Grin1",
                  pt.size = 0.4) + 
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
f2 <- FeaturePlot(data.unannotated, 
                  features = "Rbfox3",
                  pt.size = 0.4) + 
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
# tsne 
t1 <- FeaturePlot(data.unannotated, 
                  features = "Grin1",
                  reduction = "tsne",
                  pt.size = 0.4) +
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
t2 <- FeaturePlot(data.unannotated, 
                  features = "Rbfox3",
                  reduction = "tsne",
                  pt.size = 0.4) +
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
plots <- list(u1,u2,f1,f2)
layout <- rbind(c(1,2),c(3,4))
grid <- grid.arrange(grobs = plots, layout_matrix = layout)

## Cluster 3 & 10: Excitatory neuron

# Top 40 markers based on p_val_adj and then avg_log2FC
VlnPlot(data.unannotated,
        features = cluster3$gene[1:20],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "SCT_snn_res.0.1")

VlnPlot(data.unannotated,
        features = cluster3$gene[21:40],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "SCT_snn_res.0.1")

VlnPlot(data.unannotated,
        features = cluster10$gene[1:20],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "SCT_snn_res.0.1")

VlnPlot(data.unannotated,
        features = cluster10$gene[21:40],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "SCT_snn_res.0.1")

VlnPlot(data.unannotated,
        features = "Slc17a7",
        cols = cluster_colors,
        split.by = "SCT_snn_res.0.1")

VlnPlot(data.unannotated,
        features = "Satb2",
        cols = cluster_colors,
        split.by = "SCT_snn_res.0.1")

VlnPlot(data.unannotated,
        features = "Col5a1",
        cols = cluster_colors,
        split.by = "SCT_snn_res.0.1")

VlnPlot(data.unannotated,
        features = "Sdk2",
        cols = cluster_colors,
        split.by = "SCT_snn_res.0.1")

u1 <- DimPlot(data.unannotated,
              cols = cluster_colors,
              label = TRUE) + theme(legend.position = "none")
u2 <- DimPlot(data.unannotated,
              dims = c(2,3),
              label = TRUE,
              cols = cluster_colors) + theme(legend.position = "none")
f1 <- FeaturePlot(data.unannotated, 
                  features = "Slc17a7",
                  pt.size = 0.4) + 
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
f2 <- FeaturePlot(data.unannotated, 
                  features = "Satb2",
                  pt.size = 0.4) + 
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
f3 <- FeaturePlot(data.unannotated, 
                  features = "Col5a1",
                  pt.size = 0.4) + 
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
f4 <- FeaturePlot(data.unannotated, 
                  features = "Sdk2",
                  pt.size = 0.4) + 
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
# tsne 
t1 <- FeaturePlot(data.unannotated, 
                  features = "Slc17a7",
                  reduction = "tsne",
                  pt.size = 0.4) +
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
t2 <- FeaturePlot(data.unannotated, 
                  features = "Satb2",
                  reduction = "tsne",
                  pt.size = 0.4) +
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
t3 <- FeaturePlot(data.unannotated, 
                  features = "Col5a1",
                  reduction = "tsne",
                  pt.size = 0.4) +
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
t4 <- FeaturePlot(data.unannotated, 
                  features = "Sdk2",
                  reduction = "tsne",
                  pt.size = 0.4) +
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
plots <- list(u1,u2,f1,f2)
layout <- rbind(c(1,2),c(3,4))
grid <- grid.arrange(grobs = plots, layout_matrix = layout)

plots <- list(f1,f2,t1,t2)
layout <- rbind(c(1,2),c(3,4))
grid <- grid.arrange(grobs = plots, layout_matrix = layout)

plots <- list(f3,f4,t3,t4)
layout <- rbind(c(1,2),c(3,4))
grid <- grid.arrange(grobs = plots, layout_matrix = layout)

## Cluster 5,8,12: Inhibitory neuron

# Top 20 markers based on p_val_adj and then avg_log2FC
VlnPlot(data.unannotated,
        features = cluster5$gene[1:20],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "SCT_snn_res.0.1")

VlnPlot(data.unannotated,
        features = cluster8$gene[1:20],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "SCT_snn_res.0.1")

VlnPlot(data.unannotated,
        features = cluster12$gene[1:20],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "SCT_snn_res.0.1")

VlnPlot(data.unannotated,
        features = "Gad1",
        cols = cluster_colors,
        split.by = "SCT_snn_res.0.1")

VlnPlot(data.unannotated,
        features = "Gad2",
        cols = cluster_colors,
        split.by = "SCT_snn_res.0.1")

VlnPlot(data.unannotated,
        features = "Tac1",
        cols = cluster_colors,
        split.by = "SCT_snn_res.0.1")

VlnPlot(data.unannotated,
        features = "Penk",
        cols = cluster_colors,
        split.by = "SCT_snn_res.0.1")

VlnPlot(data.unannotated,
        features = "Sst",
        cols = cluster_colors,
        split.by = "SCT_snn_res.0.1")

VlnPlot(data.unannotated,
        features = "Npy",
        cols = cluster_colors,
        split.by = "SCT_snn_res.0.1")

u1 <- DimPlot(data.unannotated,
              cols = cluster_colors,
              label = TRUE) + theme(legend.position = "none")
u2 <- DimPlot(data.unannotated,
              dims = c(2,3),
              label = TRUE,
              cols = cluster_colors) + theme(legend.position = "none")
f1 <- FeaturePlot(data.unannotated, 
                  features = "Gad1",
                  pt.size = 0.4) + 
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
f2 <- FeaturePlot(data.unannotated, 
                  features = "Gad2",
                  pt.size = 0.4) + 
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
f3 <- FeaturePlot(data.unannotated, 
                  features = "Tac1",
                  pt.size = 0.4) + 
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
f4 <- FeaturePlot(data.unannotated, 
                  features = "Penk",
                  pt.size = 0.4) + 
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
f5 <- FeaturePlot(data.unannotated, 
                  features = "Sst",
                  pt.size = 0.4) + 
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
f6 <- FeaturePlot(data.unannotated, 
                  features = "Npy",
                  pt.size = 0.4) + 
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
# tsne 
t1 <- FeaturePlot(data.unannotated, 
                  features = "Gad1",
                  reduction = "tsne",
                  pt.size = 0.4) +
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
t2 <- FeaturePlot(data.unannotated, 
                  features = "Gad2",
                  reduction = "tsne",
                  pt.size = 0.4) +
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
t3 <- FeaturePlot(data.unannotated, 
                  features = "Tac1",
                  reduction = "tsne",
                  pt.size = 0.4) +
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
t4 <- FeaturePlot(data.unannotated, 
                  features = "Penk",
                  reduction = "tsne",
                  pt.size = 0.4) +
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
t5 <- FeaturePlot(data.unannotated, 
                  features = "Sst",
                  reduction = "tsne",
                  pt.size = 0.4) +
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
t6 <- FeaturePlot(data.unannotated, 
                  features = "Npy",
                  reduction = "tsne",
                  pt.size = 0.4) +
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
plots <- list(u1,u2,f1,f2)
layout <- rbind(c(1,2),c(3,4))
grid <- grid.arrange(grobs = plots, layout_matrix = layout)

plots <- list(f1,f2,t1,t2)
layout <- rbind(c(1,2),c(3,4))
grid <- grid.arrange(grobs = plots, layout_matrix = layout)

plots <- list(f3,f4,t3,t4)
layout <- rbind(c(1,2),c(3,4))
grid <- grid.arrange(grobs = plots, layout_matrix = layout)

plots <- list(f5,f6,t5,t6)
layout <- rbind(c(1,2),c(3,4))
grid <- grid.arrange(grobs = plots, layout_matrix = layout)

## Cluster 6, 11, & 12: Microglia

# Top 40 markers based on p_val_adj and then avg_log2FC
VlnPlot(data.unannotated,
        features = cluster6$gene[1:20],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "SCT_snn_res.0.1")

VlnPlot(data.unannotated,
        features = cluster6$gene[21:40],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "SCT_snn_res.0.1")

VlnPlot(data.unannotated,
        features = "C1qa",
        cols = cluster_colors,
        split.by = "SCT_snn_res.0.1")

VlnPlot(data.unannotated,
        features = "Csf1r",
        cols = cluster_colors,
        split.by = "SCT_snn_res.0.1")

u1 <- DimPlot(data.unannotated,
              cols = cluster_colors,
              label = TRUE) + theme(legend.position = "none")
u2 <- DimPlot(data.unannotated,
              dims = c(2,3),
              label = TRUE,
              cols = cluster_colors) + theme(legend.position = "none")
f1 <- FeaturePlot(data.unannotated, 
                  features = "C1qa",
                  pt.size = 0.4) + 
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
f2 <- FeaturePlot(data.unannotated, 
                  features = "Csf1r",
                  pt.size = 0.4) + 
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
# tsne 
t1 <- FeaturePlot(data.unannotated, 
                  features = "C1qa",
                  reduction = "tsne",
                  pt.size = 0.4) +
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
t2 <- FeaturePlot(data.unannotated, 
                  features = "Csf1r",
                  reduction = "tsne",
                  pt.size = 0.4) +
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
plots <- list(u1,u2,f1,f2)
layout <- rbind(c(1,2),c(3,4))
grid <- grid.arrange(grobs = plots, layout_matrix = layout)

plots <- list(f1,f2,t1,t2)
layout <- rbind(c(1,2),c(3,4))
grid <- grid.arrange(grobs = plots, layout_matrix = layout)

## Cluster 4: Astrocyte

# Top 40 markers based on p_val_adj and then avg_log2FC
VlnPlot(data.unannotated,
        features = cluster4$gene[1:20],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "SCT_snn_res.0.1")

VlnPlot(data.unannotated,
        features = cluster4$gene[21:40],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "SCT_snn_res.0.1")

VlnPlot(data.unannotated,
        features = "Aqp4",
        cols = cluster_colors,
        split.by = "SCT_snn_res.0.1")

VlnPlot(data.unannotated,
        features = "Gja1",
        cols = cluster_colors,
        split.by = "SCT_snn_res.0.1")

u1 <- DimPlot(data.unannotated,
              cols = cluster_colors,
              label = TRUE) + theme(legend.position = "none")
u2 <- DimPlot(data.unannotated,
              dims = c(2,3),
              label = TRUE,
              cols = cluster_colors) + theme(legend.position = "none")
f1 <- FeaturePlot(data.unannotated, 
                  features = "Aqp4",
                  pt.size = 0.4) + 
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
f2 <- FeaturePlot(data.unannotated, 
                  features = "Gja1",
                  pt.size = 0.4) + 
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
# tsne 
t1 <- FeaturePlot(data.unannotated, 
                  features = "Aqp4",
                  reduction = "tsne",
                  pt.size = 0.4) +
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
t2 <- FeaturePlot(data.unannotated, 
                  features = "Gja1",
                  reduction = "tsne",
                  pt.size = 0.4) +
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
plots <- list(u1,u2,f1,f2)
layout <- rbind(c(1,2),c(3,4))
grid <- grid.arrange(grobs = plots, layout_matrix = layout)

plots <- list(f1,f2,t1,t2)
layout <- rbind(c(1,2),c(3,4))
grid <- grid.arrange(grobs = plots, layout_matrix = layout)

Cluster 2: Oligodendrocyte

# Top 40 markers based on p_val_adj and then avg_log2FC
VlnPlot(data.unannotated,
        features = cluster2$gene[1:20],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "SCT_snn_res.0.1")

VlnPlot(data.unannotated,
        features = cluster2$gene[21:40],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "SCT_snn_res.0.1")

VlnPlot(data.unannotated,
        features = "Mog",
        cols = cluster_colors,
        split.by = "SCT_snn_res.0.1")

VlnPlot(data.unannotated,
        features = "Cldn11",
        cols = cluster_colors,
        split.by = "SCT_snn_res.0.1")

u1 <- DimPlot(data.unannotated,
              cols = cluster_colors,
              label = TRUE) + theme(legend.position = "none")
u2 <- DimPlot(data.unannotated,
              dims = c(2,3),
              label = TRUE,
              cols = cluster_colors) + theme(legend.position = "none")
f1 <- FeaturePlot(data.unannotated, 
                  features = "Mog",
                  pt.size = 0.4) + 
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
f2 <- FeaturePlot(data.unannotated, 
                  features = "Cldn11",
                  pt.size = 0.4) + 
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
# tsne 
t1 <- FeaturePlot(data.unannotated, 
                  features = "Mog",
                  reduction = "tsne",
                  pt.size = 0.4) +
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
t2 <- FeaturePlot(data.unannotated, 
                  features = "Cldn11",
                  reduction = "tsne",
                  pt.size = 0.4) +
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
plots <- list(u1,u2,f1,f2)
layout <- rbind(c(1,2),c(3,4))
grid <- grid.arrange(grobs = plots, layout_matrix = layout)

plots <- list(f1,f2,t1,t2)
layout <- rbind(c(1,2),c(3,4))
grid <- grid.arrange(grobs = plots, layout_matrix = layout)

## Cluster 7: OPC

# Top 40 markers based on p_val_adj and then avg_log2FC
VlnPlot(data.unannotated,
        features = cluster7$gene[1:20],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "SCT_snn_res.0.1")

VlnPlot(data.unannotated,
        features = cluster7$gene[21:40],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "SCT_snn_res.0.1")

VlnPlot(data.unannotated,
        features = "Pdgfra",
        cols = cluster_colors,
        split.by = "SCT_snn_res.0.1")

VlnPlot(data.unannotated,
        features = "Vcan",
        cols = cluster_colors,
        split.by = "SCT_snn_res.0.1")

u1 <- DimPlot(data.unannotated,
              cols = cluster_colors,
              label = TRUE) + theme(legend.position = "none")
u2 <- DimPlot(data.unannotated,
              dims = c(2,3),
              label = TRUE,
              cols = cluster_colors) + theme(legend.position = "none")
f1 <- FeaturePlot(data.unannotated, 
                  features = "Pdgfra",
                  pt.size = 0.4) + 
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
f2 <- FeaturePlot(data.unannotated, 
                  features = "Vcan",
                  pt.size = 0.4) + 
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
# tsne 
t1 <- FeaturePlot(data.unannotated, 
                  features = "Pdgfra",
                  reduction = "tsne",
                  pt.size = 0.4) +
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
t2 <- FeaturePlot(data.unannotated, 
                  features = "Vcan",
                  reduction = "tsne",
                  pt.size = 0.4) +
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
plots <- list(u1,u2,f1,f2)
layout <- rbind(c(1,2),c(3,4))
grid <- grid.arrange(grobs = plots, layout_matrix = layout)

plots <- list(f1,f2,t1,t2)
layout <- rbind(c(1,2),c(3,4))
grid <- grid.arrange(grobs = plots, layout_matrix = layout)

Cluster 9: Endothelial cell

# Top 40 markers based on p_val_adj and then avg_log2FC
VlnPlot(data.unannotated,
        features = cluster9$gene[1:20],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "SCT_snn_res.0.1")

VlnPlot(data.unannotated,
        features = cluster9$gene[21:40],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "SCT_snn_res.0.1")

VlnPlot(data.unannotated,
        features = "Vtn",
        cols = cluster_colors,
        split.by = "SCT_snn_res.0.1")

VlnPlot(data.unannotated,
        features = "Itm2a",
        cols = cluster_colors,
        split.by = "SCT_snn_res.0.1")

u1 <- DimPlot(data.unannotated,
              cols = cluster_colors,
              label = TRUE) + theme(legend.position = "none")
u2 <- DimPlot(data.unannotated,
              dims = c(2,3),
              label = TRUE,
              cols = cluster_colors) + theme(legend.position = "none")
f1 <- FeaturePlot(data.unannotated, 
                  features = "Vtn",
                  pt.size = 0.4) + 
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
f2 <- FeaturePlot(data.unannotated, 
                  features = "Itm2a",
                  pt.size = 0.4) + 
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
# tsne 
t1 <- FeaturePlot(data.unannotated, 
                  features = "Vtn",
                  reduction = "tsne",
                  pt.size = 0.4) +
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
t2 <- FeaturePlot(data.unannotated, 
                  features = "Itm2a",
                  reduction = "tsne",
                  pt.size = 0.4) +
  scale_colour_gradientn(colours = c("grey","red","maroon"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
plots <- list(u1,u2,f1,f2)
layout <- rbind(c(1,2),c(3,4))
grid <- grid.arrange(grobs = plots, layout_matrix = layout)

plots <- list(f1,f2,t1,t2)
layout <- rbind(c(1,2),c(3,4))
grid <- grid.arrange(grobs = plots, layout_matrix = layout)

# Assign indiviudal cluster names

# Rename all identities
data.annotated <- RenameIdents(object = data.unannotated, 
                               "0" = "neuron",
                               "1" = "neuron",
                               "2" = "oligodendrocyte",
                               "3" = "neuron",
                               "4" = "astrocyte",
                               "5" = "neuron",
                               "6" = "microglia",
                               "7" = "OPC",
                               "8" = "neuron",
                               "9" = "endotehelial",
                               "10" = "neuron",
                               "11" = "microglia",
                               "12" = "microglia"
                               )
data.annotated$seurat_clusters <- Idents(data.annotated)
data.annotated$individual_clusters <- Idents(data.annotated)
# umap
umap <- DimPlot(object = data.annotated, 
                group.by = "individual_clusters",
                reduction = "umap",
                cols = cluster_colors)
umap

tann <- DimPlot(object = data.annotated, 
                group.by = "individual_clusters",
                reduction = "tsne",
                cols = cluster_colors)
tann

## Percent cells per cluster

# replicate
breplicate <- data.annotated@meta.data %>%
  group_by(individual_clusters, replicate) %>%
  dplyr::count() %>%
  group_by(individual_clusters) %>%
  dplyr::mutate(percent = 100*n/sum(n)) %>%
  ungroup() %>%
  ggplot(aes(x=individual_clusters,y=percent, fill=replicate)) +
  geom_col() +
  ggtitle("Percentage of replicate per cluster") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
breplicate

# genotype
bgenotype <- data.annotated@meta.data %>%
  group_by(individual_clusters, genotype) %>%
  dplyr::count() %>%
  group_by(individual_clusters) %>%
  dplyr::mutate(percent = 100*n/sum(n)) %>%
  ungroup() %>%
  ggplot(aes(x=individual_clusters,y=percent, fill=genotype)) +
  geom_col() +
  ggtitle("Percentage of genotype per cluster") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
bgenotype

# sample
bsample <- data.annotated@meta.data %>%
  group_by(individual_clusters, sample) %>%
  dplyr::count() %>%
  group_by(individual_clusters) %>%
  dplyr::mutate(percent = 100*n/sum(n)) %>%
  ungroup() %>%
  ggplot(aes(x=individual_clusters,y=percent, fill=sample)) +
  geom_col() +
  ggtitle("Percentage of sample per cluster") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
bsample

Cluster tree

Idents(data.annotated) <- data.annotated$individual_clusters
data.annotated <- BuildClusterTree(object = data.annotated,
                                     dims = 1:15,
                                     reorder = FALSE,
                                     reorder.numeric = FALSE)

tree <- data.annotated@tools$BuildClusterTree
tree$tip.label <- paste0(tree$tip.label)

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

## Genes of interest

goi1 <- c("Grin1","Mog","Aqp4","C1qa","Pdgfra","Vtn")

v1 <- VlnPlot(data.annotated, 
        features = goi1,
        split.by = 'individual_clusters',
        cols = cluster_colors,
        flip = TRUE,
        stack = TRUE) +
    theme(legend.position = "none", axis.text.x = element_text(angle = 45, hjust = 1))
v1
# load libraries
library(ShinyCell)

Idents(data.annotated) <- "individual_clusters"
sc.config <- createConfig(data.annotated)
makeShinyApp(data.annotated, sc.config, gene.mapping = TRUE, shiny.title = "Colonna 7month")