v3 annotations
Kimberly Olney, PhD
2023-12-26
Pipseq chemistry kit version 3
Sample NPID control female NA04-324
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("v3")
sample_color.panel <- c("#4682B4")
color.panel <- dittoColors()Read in data
See script 01 for pre-processing and QC filtering. script 02 runs FindAllMarkers. Output is Wilcox and ROC.
dataObject <- readRDS(file = paste0("../../rObjects/", sampleID, ".filtered.rds"))
markers.strict <- readRDS(file = paste0("../../results/dimensionality_reduction/markers/",
sampleID, "_FindAllMarkers_strict_logFC1_FDR0.05.rds"))
ROC.markers <- readRDS(file = paste0("../../results/dimensionality_reduction/markers/",
sampleID, "_ROC_markers.rds"))Unannotated
UMAP
ditto_umap <- dittoDimPlot(object = dataObject,
var = "seurat_clusters",
reduction.use = "umap",
do.label = TRUE,
labels.highlight = TRUE)
ditto_umap
Nuclei count per cluster
count_per_cluster <- FetchData(dataObject,
vars = c("ident", "orig.ident")) %>%
dplyr::count(ident, orig.ident) %>%
tidyr::spread(ident, n)
count_per_cluster## orig.ident 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
## 1 v3 1684 1631 1334 1290 1234 1001 953 947 766 715 414 363 228 215 161
## 15
## 1 100count_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 + 100 # 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
Finding differentially expressed features (cluster biomarkers)
Seurat can help you find markers that define clusters via differential expression (DE). By default, it identifies positive and negative markers of a single cluster (specified in ident.1), compared to all other cells. FindAllMarkers() automates this process for all clusters, but you can also test groups of clusters vs. each other, or against all cells.
The default is wilcox which identifies differentially expressed genes between two groups of cells using a Wilcoxon Rank Sum test. ## FindAllMarkers default wilxcon
# find markers for every cluster compared to all remaining cells,
# report only the positive ones
all.markers <- FindAllMarkers(dataObject, only.pos = TRUE)
# rearrange to order by cluster & filter to only include log2FC > 1 & FDR < 0.05
markers.strict <- all.markers %>%
group_by(cluster) %>%
dplyr::filter(avg_log2FC > 1 & p_val_adj < 0.05)
# compare
table(all.markers$cluster)
table(markers.strict$cluster)Get markers for each cluster
# subset
cluster0 <- markers.strict[markers.strict$cluster == 0,]
cluster1 <- markers.strict[markers.strict$cluster == 1,]
cluster2 <- markers.strict[markers.strict$cluster == 2,]
cluster3 <- markers.strict[markers.strict$cluster == 3,]
cluster4 <- markers.strict[markers.strict$cluster == 4,]
cluster5 <- markers.strict[markers.strict$cluster == 5,]
cluster6 <- markers.strict[markers.strict$cluster == 6,]
cluster7 <- markers.strict[markers.strict$cluster == 7,]
cluster8 <- markers.strict[markers.strict$cluster == 8,]
cluster9 <- markers.strict[markers.strict$cluster == 9,]
cluster10 <- markers.strict[markers.strict$cluster == 10,]
cluster11 <- markers.strict[markers.strict$cluster == 11,]
cluster12 <- markers.strict[markers.strict$cluster == 12,]
cluster13 <- markers.strict[markers.strict$cluster == 13,]
cluster14 <- markers.strict[markers.strict$cluster == 14,]
cluster15 <- markers.strict[markers.strict$cluster == 15,]
cluster16 <- markers.strict[markers.strict$cluster == 16,]ROC
Seurat has several tests for differential expression which can be set with the test.use parameter (see https://satijalab.org/seurat/articles/de_vignette for details). For example, the ROC test returns the ‘classification power’ for any individual marker (ranging from 0 - random, to 1 - perfect).
Best if run in a script as the compute takes ~3 hours to run.
ROC.markers <-
FindAllMarkers(
dataObject,
logfc.threshold = 0.25,
test.use = "roc",
only.pos = TRUE
)Get markers for each cluster
# subset
ROC.cluster0 <- ROC.markers[ROC.markers$cluster == 0,]
ROC.cluster1 <- ROC.markers[ROC.markers$cluster == 1,]
ROC.cluster2 <- ROC.markers[ROC.markers$cluster == 2,]
ROC.cluster3 <- ROC.markers[ROC.markers$cluster == 3,]
ROC.cluster4 <- ROC.markers[ROC.markers$cluster == 4,]
ROC.cluster5 <- ROC.markers[ROC.markers$cluster == 5,]
ROC.cluster6 <- ROC.markers[ROC.markers$cluster == 6,]
ROC.cluster7 <- ROC.markers[ROC.markers$cluster == 7,]
ROC.cluster8 <- ROC.markers[ROC.markers$cluster == 8,]
ROC.cluster9 <- ROC.markers[ROC.markers$cluster == 9,]
ROC.cluster10 <- ROC.markers[ROC.markers$cluster == 10,]
ROC.cluster11 <- ROC.markers[ROC.markers$cluster == 11,]
ROC.cluster12 <- ROC.markers[ROC.markers$cluster == 12,]
ROC.cluster13 <- ROC.markers[ROC.markers$cluster == 13,]
ROC.cluster14 <- ROC.markers[ROC.markers$cluster == 14,]
ROC.cluster15 <- ROC.markers[ROC.markers$cluster == 15,]Top variable genes
# print top variable genes
top100 <- dataObject@assays$SCT@var.features[1:100]
top100## [1] "CHI3L1" "DPP10" "GFAP" "AQP1"
## [5] "ROBO2" "TNR" "AQP4" "CXCL14"
## [9] "KCNIP4" "H2AC18" "CEMIP" "MT2A"
## [13] "SLC1A3" "ADAM28" "H2AC19" "LHFPL3"
## [17] "DDIT3" "SYT1" "DTNA" "SLC11A1"
## [21] "RP11-386I14.4" "CH507-528H12.1" "CLDN5" "GJA1"
## [25] "FAM189A2" "SLC1A2" "ADAMTS9" "AHCYL1"
## [29] "PTPRZ1" "MT1E" "CSMD1" "TNC"
## [33] "OPCML" "EGR1" "DSCAM" "F3"
## [37] "LRRTM4" "SLC14A1" "HIF1A-AS3" "CCK"
## [41] "PLP1" "GRIK1" "GADD45B" "ATP1A2"
## [45] "VCAN" "ZSCAN31" "NRGN" "ZFP36L1"
## [49] "VMP1" "CD44" "NRG3" "ADCY2"
## [53] "GLIS3" "HES1" "FLT1" "RGS1"
## [57] "CLU" "RP11-6L16.1" "MT1G" "ARHGAP24"
## [61] "ARC" "DOCK8" "MMP16" "LRMDA"
## [65] "USP53" "SORBS1" "FOS" "LAMA2"
## [69] "RFX4" "MT1M" "GLUL" "BEST3"
## [73] "FGF14" "SPARCL1" "CCN1" "RP11-358F13.1"
## [77] "ANGPTL4" "CNTNAP2" "PVT1" "LINC01010"
## [81] "GPC5" "NRXN1" "CNTN5" "VEGFA"
## [85] "RBFOX1" "INHBA" "CH507-513H4.1" "BCYRN1"
## [89] "ADM" "OLR1" "KIAA1217" "SYNPR"
## [93] "SNTG1" "TLR2" "RP11-909M7.3" "RALYL"
## [97] "GPX3" "CST3" "IGFBP7" "WWC1"write.table(top100, paste0("../../results/variance/", sampleID, "_top100_variable_features.txt"),
sep = "\t", row.names = FALSE, quote = FALSE)PCs
# Printing out the most variable genes driving PCs
print(x = dataObject[["pca"]],
dims = 1:5,
nfeatures = 10)## PC_ 1
## Positive: GFAP, DPP10, AQP4, AQP1, DTNA, SLC1A3, SLC1A2, AHCYL1, FAM189A2, CHI3L1
## Negative: PLP1, CNP, TF, MBP, PIP4K2A, IL1RAPL1, ST18, CNDP1, SIK3, CLDND1
## PC_ 2
## Positive: ROBO2, SYT1, KCNIP4, CSMD1, OPCML, LRRTM4, RP11-909M7.3, RBFOX1, NRXN1, FGF14
## Negative: GFAP, AQP4, AQP1, SLC1A3, CHI3L1, DTNA, AHCYL1, GJA1, F3, FAM189A2
## PC_ 3
## Positive: ADAM28, SLC11A1, ARHGAP24, DOCK8, LRMDA, H2AC18, H2AC19, TLR2, APBB1IP, DENND3
## Negative: PLP1, CNP, SYT1, NRGN, CLU, TF, AQP4, AQP1, FTH1, SNAP25
## PC_ 4
## Positive: TNR, LHFPL3, CH507-528H12.1, PTPRZ1, CH507-513H4.1, DSCAM, VCAN, IL1RAPL1, SIK3, MMP16
## Negative: ADAM28, H2AC18, SLC11A1, H2AC19, DOCK8, DENND3, PLP1, SPP1, LRMDA, NRGN
## PC_ 5
## Positive: ROBO2, SYT1, DPP10, CH507-528H12.1, CH507-513H4.1, NRGN, MT-RNR2, BCYRN1, RBFOX1, FRMPD4
## Negative: TNR, LHFPL3, PTPRZ1, DSCAM, VCAN, MMP16, BEST3, MEGF11, CA10, OPCMLPlot markers
Heatmap
# get the top 10 markers for each cluster
top10_Markers <- markers.strict %>%
group_by(cluster) %>%
dplyr::filter(avg_log2FC > 2) %>%
slice_head(n = 5) %>%
ungroup()
# DoHeatmap function for plotting
heat_markers <-
DoHeatmap(dataObject, features = top10_Markers$gene) + NoLegend()
# downsample to 100 nuclei per cluster.
heat_markers_downsample <-
DoHeatmap(subset(dataObject, downsample = 100),
features = top10_Markers$gene,
size = 3)
top10_ROC <- ROC.markers %>%
group_by(cluster) %>%
dplyr::filter(avg_log2FC > 2 & myAUC > 0.8) %>%
slice_head(n = 5) %>%
ungroup()
heat_ROC <-
DoHeatmap(dataObject, features = top10_ROC$gene) + NoLegend()
heat_ROC_downsample <-
DoHeatmap(subset(dataObject, downsample = 100),
features = top10_ROC$gene,
size = 3)
heat_markers
heat_markers_downsample
heat_ROC
heat_ROC_downsample
Ridge plot
# MALAT1 is highly expressed in brain nuclei
RidgePlot(dataObject, features = c("XIST", "MALAT1"), ncol = 2)
RidgePlot(dataObject, features = c("GFAP"), ncol = 1)
RidgePlot(dataObject, features = c("VCAN"), ncol = 1)
Feature plot
# UMAP showing the expression of select features
umap_feature <-
FeaturePlot(dataObject,
features = c("TYROBP", "MOG", "AQP4", "ENO2", "XIST", "MALAT1"))
umap_feature
Dot plot markers
markers.to.plot <-
c(
"SLC1A2",
"GJA1",
"CLDN5",
"TGM2",
"STAB1",
"CD74",
"VCAN",
"TNR",
"S100B",
"ABCA2",
"NRG1",
"MEG3",
"ERBB4",
"GAD1"
)
dot <- DotPlot(dataObject, features = markers.to.plot,
dot.scale = 8) + RotatedAxis()
dot
Violins - Canonical cell-type markers
Markers obtained from dropviz
# Astrocytes
VlnPlot(dataObject,
features = c("AQP4", "GJA1", "F3"))
# Endothelial
VlnPlot(dataObject,
features = c("CLDN5", "ADGRF5", "FLT1"))
# Fibroblast-Like_Dcn
VlnPlot(dataObject,
features = c("DCN", "SLC6A13", "IGF2"))
# Oligodendrocyte
VlnPlot(dataObject,
features = c("MOG", "OLIG1", "MBP"))
# Polydendrocyte
VlnPlot(dataObject,
features = c("VCAN", "TNR", "C1QL1"))
# Smooth muscle
VlnPlot(dataObject,
features = c("ACTA2", "RGS5"))
# Microglia / Marchophage
VlnPlot(dataObject,
features = c("CTSS", "C1QC", "C1QB"))
# Microglia / Marchophage
VlnPlot(dataObject,
features = c("TYROBP","P2RY12", "AIF1"))
# Interneuron
VlnPlot(dataObject,
features = c("VIP", "CRH", "KIT"))
# Interneuron - GAD1 & GAD2
VlnPlot(dataObject,
features = c("GAD1", "GAD2"))
Cluster Annotation
Cluster 0 - neuron
CTNNA3 - mural
POLR2F - neurogenesis/ polydendrocyte
IL1RAPL1 - neuron
RNF220 - neuron
CCDC88A - astrocyte
# Number of cells per condition
count_per_cluster[,c(1,2)]## orig.ident 0
## 1 v3 1684# 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")
VlnPlot(dataObject,
features = ROC.cluster0$gene[1:10],
cols = color.panel,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
cluster0$gene[1:10]## [1] "ENSG10010138868.1" "CH507-513H4.1" "CTNNA3"
## [4] "CCDC88A" "IL1RAPL1" "CH507-528H12.1"
## [7] "MT-RNR2" "POLR2F" "SLCO5A1"
## [10] "RNF220"ROC.cluster0$gene[1:10]## [1] "CH507-528H12.1" "CH507-513H4.1" "ENSG10010138868.1"
## [4] "MT-RNR2" "CTNNA3" "IL1RAPL1"
## [7] "CCDC88A" "RNF220" "POLR2F"
## [10] "PALM2AKAP2"Cluster 1 - oligodendrocyte
VWA1 - endothelial
PLP1 - oligodendrocyte
MAG - polydendrocyte
CLDN1 - polydendrocyte / oligodendrocyte
APOD - fibroblast-like / oligodendrocyte
CNP - polydendrocyte / oligodendrocyte
# Number of cells per condition
count_per_cluster[,c(1,3)]## orig.ident 1
## 1 v3 1631DimPlot(object = subset(dataObject, seurat_clusters == "1"),
reduction = "umap",
label = TRUE,
label.box = TRUE,
label.size = 3,
repel = TRUE,
cols = "#56B4E9")
VlnPlot(dataObject,
features = cluster1$gene[1:10],
cols = color.panel,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
VlnPlot(dataObject,
features = ROC.cluster1$gene[1:10],
cols = color.panel,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
cluster1$gene[1:10]## [1] "VWA1" "LINC00844" "APOD" "TSC22D4" "RNASE1" "MAG"
## [7] "CBR1" "FTH1P10" "SELENOP" "CD9"ROC.cluster1$gene[1:10]## [1] "CNP" "PLP1" "FTH1" "APOD" "CLDND1" "SUN2" "MAG" "APLP1"
## [9] "GPRC5B" "GPM6B"Cluster 2 - oligodendrocyte, non-specific
AATK - polydendrocyte
DYSF - mural/endothelial
# Number of cells per condition
count_per_cluster[,c(1,4)]## orig.ident 2
## 1 v3 1334DimPlot(object = subset(dataObject, seurat_clusters == "2"),
reduction = "umap",
label = TRUE,
label.box = TRUE,
label.size = 3,
repel = TRUE,
cols = "#009E73")
VlnPlot(dataObject,
features = cluster2$gene[1:10],
cols = color.panel,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
VlnPlot(dataObject,
features = ROC.cluster2$gene[1:10],
cols = color.panel,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
cluster2$gene[1:10]## [1] "AATK" "MT-ND2" "MT-CO2" "DYSF"
## [5] "LINC00685" "NMNAT2" "MTCO2P12" "RP11-506B6.7"
## [9] "CASC6" "RP11-587D21.4"ROC.cluster2$gene[1:10]## [1] "AATK" "MT-ND2" "MT-CO2" NA NA NA NA NA
## [9] NA NACluster 3 - oligodendrocyte / polydendrocyte - non-specific
SLC5A11 - oligodendrocyte / polydendrocyte / astrocyte
SAMD12 - neuron endothelial
# Number of cells per condition
count_per_cluster[,c(1,4)]## orig.ident 2
## 1 v3 1334DimPlot(object = subset(dataObject, seurat_clusters == "3"),
reduction = "umap",
label = TRUE,
label.box = TRUE,
label.size = 3,
repel = TRUE,
cols = "#F0E442")
VlnPlot(dataObject,
features = cluster3$gene[1:10],
cols = color.panel,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
VlnPlot(dataObject,
features = ROC.cluster3$gene[1:10],
cols = color.panel,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
cluster3$gene[1:10]## [1] "SLC5A11" "SAMD12" "LDB3" "LINC00609" "PDIA2"
## [6] "RP1-223B1.1" "SEC14L5" "PLCH2" "NMNAT2" "TMEM63A"ROC.cluster3$gene[1:10]## [1] "SLC5A11" "AATK" NA NA NA NA NA
## [8] NA NA NACluster 4 - fibroblast-like
GADD45B - fibroblast-like
FGFR1 - fibroblast-like
EGR1 - neuron
GLUL - astrocyte
AKAP6 - neuron
# Number of cells per condition
count_per_cluster[,c(1,5)]## orig.ident 3
## 1 v3 1290DimPlot(object = subset(dataObject, seurat_clusters == "4"),
reduction = "umap",
label = TRUE,
label.box = TRUE,
label.size = 3,
repel = TRUE,
cols = "#0072B2")
VlnPlot(dataObject,
features = cluster4$gene[1:10],
cols = color.panel,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
VlnPlot(dataObject,
features = ROC.cluster4$gene[1:10],
cols = color.panel,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
cluster4$gene[1:10]## [1] "GADD45B" "FGFR1" "EGR1" "DDIT3" "GLUL" "C11orf96"
## [7] "ARC" "NEAT1" "FOS" "PRUNE2"ROC.cluster4$gene[1:10]## [1] "GLUL" "GADD45B" "NEAT1" "PRUNE2" "FGFR1" "AKAP6" "CNDP1"
## [8] "ELL2" "MAT2A" "EGR1"Cluster 5 - oligodendrocyte
S100B - oligodendrocyte
STMN1 - oligodendrocyte
STMN4 - oligodendrocyte
TSC22D4 - BergmannGlia / astrocyte
SPP1 - fibroblast-like
# Number of cells per condition
count_per_cluster[,c(1,6)]## orig.ident 4
## 1 v3 1234DimPlot(object = subset(dataObject, seurat_clusters == "5"),
reduction = "umap",
label = TRUE,
label.box = TRUE,
label.size = 3,
repel = TRUE,
cols = "#D55E00")
VlnPlot(dataObject,
features = cluster5$gene[1:10],
cols = color.panel,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
cluster5$gene[1:10]## [1] "S100B" "STMN1" "STMN4" "SYS1-DBNDD2" "TSC22D4"
## [6] "PPP1R14A" "DBNDD2" "FTH1P10" "MARCKSL1" "SPP1"VlnPlot(dataObject,
features = ROC.cluster5$gene[1:10],
cols = color.panel,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
cluster5$gene[1:10]## [1] "S100B" "STMN1" "STMN4" "SYS1-DBNDD2" "TSC22D4"
## [6] "PPP1R14A" "DBNDD2" "FTH1P10" "MARCKSL1" "SPP1"ROC.cluster5$gene[1:10]## [1] "TF" "SPP1" "S100B" "TUBA1A" "MARCKSL1" "PLP1"
## [7] "FTL" "CNP" "FTH1" "CRYAB"Cluster 6 - astrocyte
SLC1A2 - astrocyte  AQP4 - astrocyte
SLC1A3 - astrocyte
GFAP - astrocyte
DTNA - astrocyte
# Number of cells per condition
count_per_cluster[,c(1,7)]## orig.ident 5
## 1 v3 1001DimPlot(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")
VlnPlot(dataObject,
features = ROC.cluster6$gene[1:10],
cols = color.panel,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
cluster6$gene[1:10]## [1] "AQP4" "FAM189A2" "ADCY2" "SORBS1" "GFAP" "RFX4"
## [7] "F3" "MGST1" "NRG3" "SLC1A2"ROC.cluster6$gene[1:10]## [1] "DTNA" "SORBS1" "GFAP" "AQP4" "AHCYL1" "SLC1A3"
## [7] "ADCY2" "CLU" "FAM189A2" "SPARCL1"Cluster 7 - oligodendrocyte/ neuron
PCDH9 - oligodendrocyte
IL1RAPL1 - neuron
CTNNA3 - mural
ST18 - granular neuron
UNC5C - interneurons / neuron
# Number of cells per condition
count_per_cluster[,c(1,8)]## orig.ident 6
## 1 v3 953DimPlot(object = subset(dataObject, seurat_clusters == "7"),
reduction = "umap",
label = TRUE,
label.box = TRUE,
label.size = 3,
repel = TRUE,
cols = "#666666")
VlnPlot(dataObject,
features = cluster7$gene[1:10],
cols = color.panel,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
VlnPlot(dataObject,
features = ROC.cluster7$gene[1:10],
cols = color.panel,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
cluster7$gene[1:10]## [1] "PCDH9" "IL1RAPL1" "CTNNA3"
## [4] "ST18" "ENSG10010138868.1" "UNC5C"
## [7] "MAN2A1" "PLCL1" "PDE1C"
## [10] "PEX5L"ROC.cluster7$gene[1:10]## [1] "MALAT1" "PCDH9" "IL1RAPL1"
## [4] "CTNNA3" "ST18" "CH507-528H12.1"
## [7] "CH507-513H4.1" "DLG2" "UNC5C"
## [10] "ENSG10010138868.1"Cluster 8 - oligodendrocyte
APLP1 - oligodendrocyte
ABCA2 - oligodendrocyte
SUN2 - microglia
TMEM151A - oligodendrocyte
# Number of cells per condition
count_per_cluster[,c(1,9)]## orig.ident 7
## 1 v3 947DimPlot(object = subset(dataObject, seurat_clusters == "8"),
reduction = "umap",
label = TRUE,
label.box = TRUE,
label.size = 3,
repel = TRUE,
cols = "#AD7700")
VlnPlot(dataObject,
features = cluster8$gene[1:10],
cols = color.panel,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
VlnPlot(dataObject,
features = ROC.cluster8$gene[1:10],
cols = color.panel,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
cluster8$gene[1:10]## [1] "ABCA2" "APLP1" "SUN2" "TUBB4A"
## [5] "CARNS1" "RP11-339B21.12" "NKX6-2" "HAPLN2"
## [9] "LGI3" "TMEM151A"ROC.cluster8$gene[1:10]## [1] "PLP1" "ABCA2" "CNP" "APLP1"
## [5] "CLDND1" "TF" "RP11-229P13.28" "GPM6B"
## [9] "SUN2" "MBP"Cluster 9 - neuron
PRUNE2 - neuron
SIK3 - neuron
AKAP6 - purkinje neuron
SH3RF1 - neuron / polydendrocyte
# Number of cells per condition
count_per_cluster[,c(1,10)]## orig.ident 8
## 1 v3 766DimPlot(object = subset(dataObject, seurat_clusters == "9"),
reduction = "umap",
label = TRUE,
label.box = TRUE,
label.size = 3,
repel = TRUE,
cols = "#1C91D4")
VlnPlot(dataObject,
features = cluster9$gene[1:10],
cols = color.panel,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
VlnPlot(dataObject,
features = ROC.cluster9$gene[1:10],
cols = color.panel,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
cluster9$gene[1:10]## [1] "PRUNE2" "SIK3" "ENSG10010138868.1"
## [4] "AKAP6" "SH3RF1" "SH3PXD2B"
## [7] "JMJD1C" "TMTC2" "GTF2H2B"
## [10] "DLG1"ROC.cluster9$gene[1:10]## [1] "SIK3" "PRUNE2" "MALAT1"
## [4] "ENSG10010138868.1" "AKAP6" "TMTC2"
## [7] "JMJD1C" "DLG1" "CH507-528H12.1"
## [10] "CH507-513H4.1"Cluster 10 - neuron
SNAP25 - neuron
SYT1 - neuron
VSNL1 - neuron
CHN1 - neuron
# Number of cells per condition
count_per_cluster[,c(1,11)]## orig.ident 9
## 1 v3 715DimPlot(object = subset(dataObject, seurat_clusters == "10"),
reduction = "umap",
label = TRUE,
label.box = TRUE,
label.size = 3,
repel = TRUE,
cols = "#007756")
VlnPlot(dataObject,
features = cluster10$gene[1:10],
cols = color.panel,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
VlnPlot(dataObject,
features = ROC.cluster10$gene[1:10],
cols = color.panel,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
cluster10$gene[1:10]## [1] "SNAP25" "SYT1" "VSNL1" "CHN1" "STMN2" "NRGN" "SLC17A7"
## [8] "OLFM1" "BASP1" "NELL2"ROC.cluster10$gene[1:10]## [1] "SYT1" "SNAP25" "CHN1" "VSNL1" "LMO4" "NRGN" "BCYRN1"
## [8] "STMN2" "RTN1" "SLC17A7"Cluster 11 - polydendrocyte
TNR - polydendrocyte
VCAN - polydendrocyte
PTPRZ1 - polydendrocyte
LHFPL3 - polydendrocyte
DSCAM - polydendrocyte
# Number of cells per condition
count_per_cluster[,c(1,12)]## orig.ident 10
## 1 v3 414DimPlot(object = subset(dataObject, seurat_clusters == "11"),
reduction = "umap",
label = TRUE,
label.box = TRUE,
label.size = 3,
repel = TRUE,
cols = "#D5C711")
VlnPlot(dataObject,
features = cluster11$gene[1:10],
cols = color.panel,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
VlnPlot(dataObject,
features = ROC.cluster11$gene[1:10],
cols = color.panel,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
cluster11$gene[1:10]## [1] "TNR" "VCAN" "PTPRZ1" "LHFPL3" "DSCAM" "OPCML" "NRCAM" "MEGF11"
## [9] "MMP16" "PCDH15"ROC.cluster11$gene[1:10]## [1] "TNR" "PTPRZ1" "LHFPL3" "VCAN" "DSCAM" "LRRC4C" "OPCML" "NRCAM"
## [9] "CSMD1" "MEGF11"Cluster 12 - neuron
GRIP1 - interneuron
GRIK2 - interneuron
FGF14 - interneuron
PAK3 - neuron
# Number of cells per condition
count_per_cluster[,c(1,13)]## orig.ident 11
## 1 v3 363DimPlot(object = subset(dataObject, seurat_clusters == "12"),
reduction = "umap",
label = TRUE,
label.box = TRUE,
label.size = 3,
repel = TRUE,
cols = "#005685")
VlnPlot(dataObject,
features = cluster12$gene[1:10],
cols = color.panel,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
VlnPlot(dataObject,
features = ROC.cluster12$gene[1:10],
cols = color.panel,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
cluster12$gene[1:10]## [1] "RP11-909M7.3" "GRIP1" "GRIK2" "FGF14" "PAK3"
## [6] "GRIA1" "CELF4" "MYT1L" "CACNA1B" "GRIN2B"ROC.cluster12$gene[1:10]## [1] "SNHG14" "CNTNAP2" "DLGAP1" "GRIK2" "RP11-909M7.3"
## [6] "SYT1" "PLCB1" "GRIP1" "FGF14" "NRXN3"Cluster 13 - microglia
SLC11A1 - microglia / macrophage
ADAM28 - interneuron
ARHGAP24 - polydendroctes
PLXDC2 - microglia/ macrophage
DENND3 - enothelial
# Number of cells per condition
count_per_cluster[,c(1,14)]## orig.ident 12
## 1 v3 228DimPlot(object = subset(dataObject, seurat_clusters == "13"),
reduction = "umap",
label = TRUE,
label.box = TRUE,
label.size = 3,
repel = TRUE,
cols = "#A04700")
VlnPlot(dataObject,
features = cluster13$gene[1:10],
cols = color.panel,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
VlnPlot(dataObject,
features = ROC.cluster13$gene[1:10],
cols = color.panel,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
cluster13$gene[1:10]## [1] "SLC11A1" "ADAM28" "LRMDA" "ARHGAP24" "DENND3" "DOCK8"
## [7] "RNASET2" "CHST11" "INPP5D" "APBB1IP"ROC.cluster13$gene[1:10]## [1] "SLC11A1" "ADAM28" "PLXDC2" "LRMDA" "ARHGAP24" "SRGAP2"
## [7] "DENND3" "ARHGAP26" "CELF2" "SRGAP2B"Cluster 14 - neuron
HS3ST4 - neuron
CELF4 - neuron
RALYL - interneuron
SYT1 - neuron  SNHG14 - neuron
# Number of cells per condition
count_per_cluster[,c(1,15)]## orig.ident 13
## 1 v3 215DimPlot(object = subset(dataObject, seurat_clusters == "14"),
reduction = "umap",
label = TRUE,
label.box = TRUE,
label.size = 3,
repel = TRUE,
cols = "#B14380")
VlnPlot(dataObject,
features = cluster14$gene[1:10],
cols = color.panel,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
VlnPlot(dataObject,
features = ROC.cluster14$gene[1:10],
cols = color.panel,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
cluster14$gene[1:10]## [1] "HS3ST4" "CELF4" "KIAA1217" "RALYL" "SCN2A" "ARPP21"
## [7] "GABRB2" "MYT1L" "SYT7" "MIAT"ROC.cluster14$gene[1:10]## [1] "SYT1" "SNHG14" "CSMD1" "RBFOX1" "RP11-909M7.3"
## [6] "MEG3" "MDGA2" "HS3ST4" "KIAA1217" "NRXN1"Cluster 15 - fibroblast-like
NID1 - fibroblast-like
LEF1 - endothelial
IFITM1 - mural / fibroblast-like
CFH - fibroblast-like
DCN - fibroblast-like
# Number of cells per condition
count_per_cluster[,c(1,16)]## orig.ident 14
## 1 v3 161DimPlot(object = subset(dataObject, seurat_clusters == "15"),
reduction = "umap",
label = TRUE,
label.box = TRUE,
label.size = 3,
repel = TRUE,
cols = "#4D4D4D")
VlnPlot(dataObject,
features = cluster15$gene[1:10],
cols = color.panel,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
VlnPlot(dataObject,
features = ROC.cluster15$gene[1:10],
cols = color.panel,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
cluster15$gene[1:10]## [1] "NID1" "LEF1" "IFITM1" "CFH" "DCN" "COL1A2" "ADH1B" "C7"
## [9] "CLDN5" "PDLIM1"ROC.cluster15$gene[1:10]## [1] "RBMS3" "PTPRG" "USP53" "RP11-9J18.1" "RBPMS"
## [6] "ARHGAP29" "CRIM1" "NID1" "EPS8" "SVIL"Assign identities
# Rename all identities
dataObject.annotated <- RenameIdents(object = dataObject,
"0" = "Neurons",
"1" = "Oligodendrocytes",
"2" = "Oligodendrocytes",
"3" = "Oligodendrocytes / Polydendrocyte",
"4" = "Fibroblast-like",
"5" = "Oligodendrocytes",
"6" = "Astrocytes",
"7" = "Oligodendrocytes",
"8" = "Oligodendrocytes",
"9" = "Neurons",
"10" = "Neurons",
"11" = "Polydendrocytes",
"12" = "Neurons",
"13" = "Microglia",
"14" = "Neurons",
"15" = "Fibroblast-like")
dataObject.annotated$annotated_clusters <- factor(Idents(dataObject.annotated),
levels = c("Astrocytes",
"Fibroblast-like",
"Microglia",
"Oligodendrocytes",
"Polydendrocytes",
"Oligodendrocytes / Polydendrocyte",
"Neurons" ))
Idents(dataObject.annotated) <- "annotated_clusters"Unannotated
UMAP
ditto_umap <- dittoDimPlot(object = dataObject.annotated,
var = "annotated_clusters",
reduction.use = "umap",
do.label = TRUE,
labels.highlight = TRUE)
ditto_umap
Dot plot markers
markers.to.plot <-
c(
"SLC1A2",
"GJA1",
"CLDN5",
"TGM2",
"STAB1",
"CD74",
"VCAN",
"TNR",
"S100B",
"ABCA2",
"NRG1",
"MEG3",
"ERBB4",
"GAD1"
)
dot <- DotPlot(dataObject.annotated, features = markers.to.plot,
dot.scale = 8) + RotatedAxis()
dot
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 Astrocytes Fibroblast-like Microglia Oligodendrocytes
## 1 v3 953 1334 215 5679
## Polydendrocytes Oligodendrocytes / Polydendrocyte Neurons
## 1 363 1290 3202count_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
## 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() %>%
ggplot(aes(x=orig.ident,y=percent, fill=annotated_clusters)) +
geom_col() +
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))
relative_abundance
Save RDS
saveRDS(dataObject.annotated, paste0("../../rObjects/", sampleID, "_annotated.rds"))