Annotation cluster res 0.2
NadineBestard
30/03/2021
library(here) #reproducible paths
library(scater) #feature plots
library(patchwork) # agregate plots
library(scran) # for findmarkers
library(pals) # for palettes with large n #kelly()22, #polychrome()#36, cols
library(tidyverse) # edit project <- "fire-mice"cols25 <- unname(cols25())
# remove the black and white from the pallete, and the similars to cols25
# assessed with pal.bands and pal.cube
kelly_col <- unname(kelly()[-c(1,2,5,6)])
# remove the colours that are similar to the
cols25 <-
cols <- c(kelly_col, cols25)if (!file.exists(here("processed", project, "sce_anno_02_res02.RDS"))) {
sce <- readRDS(here("processed", project, "sce_clusters_02.RDS"))
}else{
sce <- readRDS(here("processed", project, "sce_anno_02_res02.RDS"))
}Annotation
First annotation with known markers
plotTSNE(sce, colour_by = "originalexp_snn_res.0.2", text_by = "originalexp_snn_res.0.2", force = 0) + scale_color_manual(values = cols)## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
Click to expand the Neurons marker plots
#Neurons:
list_plots <- lapply(c("Snap25", "Stmn2", "Rbfox3", "Gabrb2"),
function(x)plotTSNE(sce, colour_by = x ))
wrap_plots(list_plots) + plot_annotation(title = "Neurons")
plotExpression(sce, features=c("Snap25", "Stmn2", "Rbfox3", "Gabrb2") ,
x="originalexp_snn_res.0.2", colour_by = "originalexp_snn_res.0.2", ncol=1) + scale_color_manual(values = cols)## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
#Inhibitory neurons
list_plots <- lapply(c("Gad1", "Gad2", "Slc32a1", "Pvalb"),
function(x)plotTSNE(sce, colour_by = x ))
wrap_plots(list_plots) + plot_annotation(title = "Inhibitory Neurons")
plotExpression(sce, features=c("Gad1", "Gad2", "Slc32a1", "Pvalb"),
x="originalexp_snn_res.0.2", colour_by = "originalexp_snn_res.0.2", ncol=1) + scale_color_manual(values = cols)## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
# Excitatory neurons
list_plots <- lapply(c("Satb2", "Slc12a6", "Slc17a7"),
function(x)plotTSNE(sce, colour_by = x ))
wrap_plots(list_plots) + plot_annotation(title = "Exitatory Neurons")
plotExpression(sce, features=c("Satb2", "Slc12a6", "Slc17a7"),
x="originalexp_snn_res.0.2", colour_by = "originalexp_snn_res.0.2", ncol=1) + scale_color_manual(values = cols)## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
# RBFOX3 (and other Granulate cell markers)
list_plots <- lapply(c("Cdh15", "Calb2", "Rbfox3", "Reln"),
function(x)plotTSNE(sce, colour_by = x ))
wrap_plots(list_plots) + plot_annotation(title = "RBFO3+ Neurons")
plotExpression(sce, features=c("Cdh15", "Calb2", "Rbfox3", "Reln"),
x="originalexp_snn_res.0.2", colour_by = "originalexp_snn_res.0.2", ncol=1) + scale_color_manual(values = cols)## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
Click to expand the Stromal marker plots
#Stromal
list_plots <- lapply(c( "Lamb1" ,
"Hspg2",
"Col4a1",
"Fn1",
"Lama2"),
function(x) plotTSNE(sce, colour_by = x ))
wrap_plots(list_plots) + plot_annotation(title = "Stromal")
plotExpression(sce, features=c( "Lamb1" ,
"Hspg2",
"Col4a1",
"Fn1",
"Lama2"),
x="originalexp_snn_res.0.2", colour_by = "originalexp_snn_res.0.2", ncol=1) + scale_color_manual(values = cols)## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
# Endothelial cells and pericytes
list_plots <- lapply(c( "Cldn5",
"Icam2",
"Pdgfrb",
"Notch3",
"Vwf",
"Flt1",
"Mecom"),
function(x)plotTSNE(sce, colour_by = x ))
wrap_plots(list_plots) + plot_annotation(title = " Endothelial cells and pericytes")
plotExpression(sce, features=c( "Cldn5",
"Icam2",
"Pdgfrb",
"Notch3",
"Vwf",
"Flt1",
"Mecom"),
x="originalexp_snn_res.0.2", colour_by = "originalexp_snn_res.0.2", ncol=1) + scale_color_manual(values = cols)## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
Click to expand the Epithelial cells
list_plots <- lapply(c("Ttr", "Kcnj13", "Krt18"),
function(x)plotTSNE(sce, colour_by = x ))
wrap_plots(list_plots) + plot_annotation(title = "Choroid plexus epithelial cells")
plotExpression(sce, features=c("Ttr", "Kcnj13", "Krt18"),
x="originalexp_snn_res.0.2", colour_by = "originalexp_snn_res.0.2", ncol=1) + scale_color_manual(values = cols)## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
Click to expand the Astrocytes marker plots
# Astrocytes
list_plots <- lapply(c("Gja1",
"Aqp4",
"Glul",
"Sox9",
"Ndrg2",
"Gfap",
"Aldh1a1",
"Aldh1l1",
"Vim",
"Apoe",
"Fgfr3"),
function(x)plotTSNE(sce, colour_by = x ))
wrap_plots(list_plots) + plot_annotation(title = "Astrocyte")
plotExpression(sce, features=c("Gja1",
"Aqp4",
"Glul",
"Sox9",
"Ndrg2",
"Gfap",
"Aldh1a1",
"Aldh1l1",
"Vim",
"Apoe",
"Fgfr3"),
x="originalexp_snn_res.0.2", colour_by = "originalexp_snn_res.0.2", ncol=1) + scale_color_manual(values = cols)## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
#Astrocyte markers as described in Zeisel et al. 2018 in the mouse for telencephalon and non-telencephalon astrocytes
list_plots <- lapply(c( "Agt",
"Mfge8",
"Slc6a11",
"Slc6a9",
"Gdf10",
"Islr",
"Gfap",
"Aqp4") ,
function(x)plotTSNE(sce, colour_by = x ))
wrap_plots(list_plots) + plot_annotation(title = "Astrocyte mouse telecephalon")
Click to expand the Immune cells plots
#Microglia and macrophages
list_plots <- lapply(c( "Cd74",
"Spi1",
"Mrc1",
"Tmem119",
"Cx3cr1",
"Aif1",
"P2ry12",
"C1qc",
"C1qa"),
function(x)plotTSNE(sce, colour_by = x ))
wrap_plots(list_plots) + plot_annotation(title = "Microglia and macrophages")
plotExpression(sce, features=c( "Cd74",
"Spi1",
"Mrc1",
"Tmem119",
"Cx3cr1",
"Aif1",
"P2ry12",
"C1qc",
"C1qa"),
x="originalexp_snn_res.0.2", colour_by = "originalexp_snn_res.0.2", ncol=1) + scale_color_manual(values = cols)## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
# Border associated mcrophages
list_plots <- lapply(c( "Mrc1", "Ms4a7", "Apoe"),
function(x)plotTSNE(sce, colour_by = x ))
wrap_plots(list_plots) + plot_annotation(title = "Border associated macrophages")
plotExpression(sce, features=c( "Mrc1", "Ms4a7", "Apoe"),
x="originalexp_snn_res.0.2", colour_by = "originalexp_snn_res.0.2", ncol=1) + scale_color_manual(values = cols)## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
# Immune cells
list_plots <- lapply(c(
# Tcells
"Cd3e",
# Bcells
"Cd19",
# Natural killer
"Klrb1c", "Cd209a",
# all immune cells (CD45)
"Ptprc"),
function(x)plotTSNE(sce,
colour_by = x ,
point_alpha=0.3,
point_size = 0.5))
wrap_plots(list_plots) + plot_annotation(title = "Immune cells")
plotExpression(sce, features=c(
# Tcells
"Cd3e",
# Bcells
"Cd19",
# Natural killer
"Klrb1c", "Cd209a",
# all immune cells (CD45)
"Ptprc"),
x="originalexp_snn_res.0.2", colour_by = "originalexp_snn_res.0.2", ncol=1) + scale_color_manual(values = cols)## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
#monocybes/neutrophils
list_plots <- lapply(c(
# monocytes
"S100a9",
#neutrophils
"Ly6g", "Camp"),
function(x)plotTSNE(sce,
colour_by = x ,
point_alpha=0.3,
point_size = 0.5))
wrap_plots(list_plots) + plot_annotation(title = "Monocytes and Neutorphils cells")
plotExpression(sce, features=c(
# monocytes
"S100a9",
#neutrophils
"Ly6g", "Camp"),
x="originalexp_snn_res.0.2", colour_by = "originalexp_snn_res.0.2", ncol=1) + scale_color_manual(values = cols)## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
Click to expand the Oligodendroglia marker plots
#OPCs
list_plots <- lapply(c("Pdgfra",
"Cspg4",
"Gpr17",
"Ptprz1",
"Olig1",
"Olig2",
"Pcdh15",
"Ptgds",
"Bcan"),
function(x)plotTSNE(sce, colour_by = x ))
wrap_plots(list_plots) + plot_annotation(title = "OPCS")
plotExpression(sce, features=c("Pdgfra",
"Cspg4",
"Gpr17",
"Ptprz1",
"Olig1",
"Olig2",
"Pcdh15",
"Ptgds",
"Bcan"),
x="originalexp_snn_res.0.2", colour_by = "originalexp_snn_res.0.2", ncol=1) + scale_color_manual(values = cols)## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
#Oligodendrocytes
list_plots <- lapply(c("Plp1",
"Cnp",
"Mag",
"Mog",
"Mobp",
"Mbp",
"Sox10" ),
function(x)plotTSNE(sce, colour_by = x ))
wrap_plots(list_plots) + plot_annotation(title = "Oligodendrocytes")
plotExpression(sce, features=c("Plp1",
"Cnp",
"Mag",
"Mog",
"Mobp",
"Mbp",
"Sox10" ),
x="originalexp_snn_res.0.2", colour_by = "originalexp_snn_res.0.2", ncol=1) + scale_color_manual(values = cols)## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
Click to expand the ependymal marker plots
list_plots <- lapply(c("Vit", "Sox9", "Dynlrb2", "Ccdc153", "Rsph1", "Tm4sf1", "Pcp4l1", "Pcp4", "Hspa2", "Cd24a", "Mt2", "Chchd10"),
function(x)plotTSNE(sce, colour_by = x ))
wrap_plots(list_plots) + plot_annotation(title = "Ependymal")
Top markers and differential expression
In order to help with the identification of the different clusters we pull out the most expressed genes for each cluster.
if(!file.exists(here("outs", project, "top_genes_res0.2.tsv"))){
# Create object where instead of cells in cols we have the clusters with counts aggregated
summed_res0.2 <- aggregateAcrossCells(sce, id=colData(sce)[,c("originalexp_snn_res.0.2")])
# # Delete ribosomal and mitochondrial genes
# is_ribo <- grepl("^Rp[sl]", rownames(summed_res0.2))
# is_mito <- grepl("^mt-", rownames(summed_res0.2))
# keep <- !(is_ribo | is_mito)
# summed_res0.2 <- summed_res0.2[keep,]
# for each column (aggregated counts) of the object extract the top genes
top_genes_res0.2 <- apply(counts(summed_res0.2),
MARGIN = 2,
function(x){
# get the gene names from the top 50 expressed genes
names(sort(x, decreasing = T)[1:100])
}
)
# save result
write.table(top_genes_res0.2, file = here("outs", project, "top_genes_res0.2.tsv"), sep = "\t", quote = FALSE, col.names = TRUE, row.names = FALSE)
}We compute here a pairwise differential expression between all the clusters. The results are saved in a list of dataframes, “markers_res0.2/markers_res0.2.RDS” one for each cluster. A .csv file (that can be opened with a spreadsheet program) for each cluster is available in “markers_res0.2”.
The default philosophy of findMarkers() is to identify a combination of marker genes that - together - uniquely define one cluster against the rest. To this end, we collect the top DE genes from each pairwise comparison involving a particular cluster to assemble a set of candidate markers for that cluster. Of particular interest is the Top field. The set of genes with Top \(≤X\) is the union of the top \(X\) genes (ranked by p-value) from each pairwise comparison involving the cluster of interest. For example, the set of all genes with Top values of 1 contains the gene with the lowest p-value from each comparison. Similarly, the set of genes with Top values less than or equal to 10 contains the top 10 genes from each comparison.
# compute markers
if(!file.exists(here("outs", project, "markers_res0.2", "markers_res0.2.RDS"))){
dir.create(here("outs", project, "markers_res0.2"))
markers <- findMarkers(sce, groups = sce$originalexp_snn_res.0.2, direction="up", lfc=1)
saveRDS(markers, here("outs", project,"markers_res0.2", "markers_res0.2.RDS"))
# save the top 100 genes for each one of the df
lapply(names(markers), function(x){
top_markers <- head(markers[[x]], 100)
write.csv(top_markers, here("outs", project, "markers_res0.2", paste0("top_markers_", x, "_res0.2.csv")), quote = FALSE)
}
)
}Rename the Clusters with assigned Cell Types
With the help of the known markers plotted above, the DE between clusters and the top expressed genes the celltype identity of every cluster has been identified. # in progress
plotTSNE(sce, colour_by = "originalexp_snn_res.0.2", text_by = "originalexp_snn_res.0.2", force = 0) + scale_color_manual(values = cols)## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
if (!file.exists(here("processed", project, "sce_anno_02_res02.RDS"))) {
annotation <- read.delim(here("data", "young_data_cluster_names_anno_02_res02.tsv"))
#celltype <- annotation$Celltype #TODO , add celltype if we keep it
clusters_named <- annotation$Cluster_name
#sce$celltype <- sce$originalexp_snn_res.0.2
sce$clusters_named <- sce$originalexp_snn_res.0.2
# add the celltypes as described above in the correct order to replace the levels.
#levels(sce$celltype) <- celltype
# and the clusters named to keep the resolution
levels(sce$clusters_named) <- clusters_named
# relevel to get the appropriate order in the plots
sce$genotype <- relevel(factor(sce$genotype), "WT")
# TODO put in the perfect order, for now just alphabetical
sce$clusters_named <-
factor(sce$clusters_named,
levels = clusters_named[order(clusters_named)])
#
# sce$celltype <-
# factor(sce$celltype,
# levels = c())
saveRDS(sce, here("processed", project, "sce_anno_02_res02.RDS"))
}
plotTSNE(sce, colour_by = "clusters_named", text_by = "clusters_named", text_size = 3, force = 1) + scale_color_manual(values = cols)## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
Marker plots
Exploring the shiny app and the DE between clusters specific markers for each one of the subtypes is identified.
annotation <- read.delim(here("data", "young_data_cluster_names_anno_02_res02.tsv"))
annotation <- annotation %>%
arrange(Cluster_name) %>%
unite(col = "markers", c("Celltype_marker","Subcluster_marker"), sep = ",", remove = FALSE, na.rm = TRUE)
markers <- annotation$markers %>% strsplit(",") %>% unlist() %>% unique()
plotExpression(sce, features=markers,
x="clusters_named", colour_by = "clusters_named", ncol=1, scales = "free_y") +
scale_color_manual(values = cols) + theme(axis.text.x = element_text(angle = 45, hjust=1, vjust = 1)) ## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
Session Info
Click to expand
sessionInfo()## R version 4.1.1 (2021-08-10)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19043)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United Kingdom.1252
## [2] LC_CTYPE=English_United Kingdom.1252
## [3] LC_MONETARY=English_United Kingdom.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United Kingdom.1252
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] forcats_0.5.1 stringr_1.4.0
## [3] dplyr_1.0.8 purrr_0.3.4
## [5] readr_2.1.2 tidyr_1.2.0
## [7] tibble_3.1.6 tidyverse_1.3.1
## [9] pals_1.7 scran_1.22.1
## [11] patchwork_1.1.1 scater_1.23.5
## [13] ggplot2_3.3.5 scuttle_1.4.0
## [15] SingleCellExperiment_1.16.0 SummarizedExperiment_1.24.0
## [17] Biobase_2.54.0 GenomicRanges_1.46.1
## [19] GenomeInfoDb_1.30.1 IRanges_2.28.0
## [21] S4Vectors_0.32.3 BiocGenerics_0.40.0
## [23] MatrixGenerics_1.6.0 matrixStats_0.61.0
## [25] here_1.0.1
##
## loaded via a namespace (and not attached):
## [1] ggbeeswarm_0.6.0 colorspace_2.0-2
## [3] ellipsis_0.3.2 rprojroot_2.0.2
## [5] bluster_1.4.0 XVector_0.34.0
## [7] fs_1.5.2 BiocNeighbors_1.12.0
## [9] dichromat_2.0-0 rstudioapi_0.13
## [11] farver_2.1.0 ggrepel_0.9.1
## [13] fansi_1.0.2 lubridate_1.8.0
## [15] xml2_1.3.3 sparseMatrixStats_1.6.0
## [17] knitr_1.37 jsonlite_1.7.3
## [19] broom_0.7.12 cluster_2.1.2
## [21] dbplyr_2.1.1 mapproj_1.2.8
## [23] compiler_4.1.1 httr_1.4.2
## [25] dqrng_0.3.0 backports_1.4.1
## [27] assertthat_0.2.1 Matrix_1.4-0
## [29] fastmap_1.1.0 limma_3.50.0
## [31] cli_3.2.0 BiocSingular_1.10.0
## [33] htmltools_0.5.2 tools_4.1.1
## [35] rsvd_1.0.5 igraph_1.2.11
## [37] gtable_0.3.0 glue_1.6.1
## [39] GenomeInfoDbData_1.2.7 maps_3.4.0
## [41] Rcpp_1.0.8 cellranger_1.1.0
## [43] jquerylib_0.1.4 vctrs_0.3.8
## [45] DelayedMatrixStats_1.16.0 xfun_0.29
## [47] rvest_1.0.2 beachmat_2.10.0
## [49] lifecycle_1.0.1 irlba_2.3.5
## [51] statmod_1.4.36 edgeR_3.36.0
## [53] zlibbioc_1.40.0 scales_1.1.1
## [55] hms_1.1.1 parallel_4.1.1
## [57] yaml_2.2.2 gridExtra_2.3
## [59] sass_0.4.0 stringi_1.7.6
## [61] highr_0.9 ScaledMatrix_1.2.0
## [63] BiocParallel_1.28.3 rlang_1.0.1
## [65] pkgconfig_2.0.3 bitops_1.0-7
## [67] evaluate_0.14 lattice_0.20-45
## [69] labeling_0.4.2 cowplot_1.1.1
## [71] tidyselect_1.1.1 magrittr_2.0.2
## [73] R6_2.5.1 generics_0.1.2
## [75] metapod_1.2.0 DelayedArray_0.20.0
## [77] DBI_1.1.2 pillar_1.7.0
## [79] haven_2.4.3 withr_2.4.3
## [81] RCurl_1.98-1.6 modelr_0.1.8
## [83] crayon_1.5.0 utf8_1.2.2
## [85] tzdb_0.2.0 rmarkdown_2.11
## [87] viridis_0.6.2 locfit_1.5-9.4
## [89] grid_4.1.1 readxl_1.3.1
## [91] reprex_2.0.1 digest_0.6.29
## [93] munsell_0.5.0 beeswarm_0.4.0
## [95] viridisLite_0.4.0 vipor_0.4.5
## [97] bslib_0.3.1