The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.
The vignette depends on Seurat package.
library(CIARA)
required <- c("Seurat")
if (!all(unlist(lapply(required, function(pkg) requireNamespace(pkg, quietly = TRUE)))))
knitr::opts_chunk$set(eval = FALSE)
In this vignette it is shown the analysis performed on single cell RNA seq human gastrula data from Tyser, R.C.V. et al., 2021.
We load the raw count matrix and umap coordinate defined in the original paper. Raw count matrix can be downloaded here
umap_elmir <- readRDS(system.file("extdata", "annot_umap.rds", package = "CIARA"))
coordinate_umap_human <- umap_elmir[, 2:3]
Obtain normalized count matrix and knn matrix ( euclidean distance on highly variable genes) using Seurat.
norm_human_data <- as.matrix(Seurat::GetAssayData(human_data_seurat, slot = "data", assay = "RNA"))
knn_human_data <- as.matrix(human_data_seurat@graphs$RNA_nn)
Cluster annotation provided in the original paper
original_cluster_human <- as.vector(umap_elmir$cluster_id)
names(original_cluster_human) <- names(umap_elmir$cell_name)
plot_umap(coordinate_umap_human, original_cluster_human)
#> [[1]]
CIARA (Cluster Independent Algorithm for the identification of RAre cell types) is a cluster independent approach that selects genes localized in a small number of neighboring cells from high dimensional PCA space. We don’t execute the CIARA algorithm and we directly load the result
ciara_genes <- row.names(result)[result[, 1]<1]
ciara_genes_top <- row.names(result)[order(as.numeric(result[, 1]))]
p <- list()
for(i in (ciara_genes_top)[1:5]) {
q <- plot_gene(norm_human_data_ciara, coordinate_umap_human, i, i)
p <- list(p, q)
}
p
#> [[1]]
#> [[1]][[1]]
#> [[1]][[1]][[1]]
#> [[1]][[1]][[1]][[1]]
#> [[1]][[1]][[1]][[1]][[1]]
#> list()
#>
#> [[1]][[1]][[1]][[1]][[2]]
#>
#>
#> [[1]][[1]][[1]][[2]]
#>
#>
#> [[1]][[1]][[2]]
#>
#>
#> [[1]][[2]]
#>
#>
#> [[2]]
It is also possible to explore wich genes are highly localized in which cells in an interactive way
norm_counts_small <- apply(norm_human_data_ciara, 1, function(x) {
y <- x/sum(x)
return(y)
})
gene_sum <- apply(norm_counts_small, 1, sum)
genes_name_text <- selection_localized_genes(norm_human_data_ciara, ciara_genes, min_number_cells = 4, max_number_genes = 4)
colnames(coordinate_umap_human) <- c("UMAP_1", "UMAP_2")
if ((requireNamespace("plotly", quietly = TRUE))) {
plot_interactive(coordinate_umap_human, gene_sum, genes_name_text, min_x = NULL, max_x = NULL, min_y = NULL, max_y = NULL)
}
We run louvain cluster analysis implemented in Seurat using as features the highly localized genes provided by CIARA
We can explore how much the number of clusters changes with the resolution. The original cluster 1 and 2 for resolution 0.01 are constantly present for higher values of resolution.
if ((requireNamespace("clustree", quietly = TRUE))) {
find_resolution(human_data_ciara, seq(0.01, 1, 0.1))
}
Cluster 2 (7 cells) is made up by primordial germ cells (PGCs). These cells expressed typical PGCs markers ad NANOS3, NANOG and DPPA5
` Merge the original cluster annotation with the one found using CIARA highly localized genes
final_cluster_human <- merge_cluster(original_cluster_human, ciara_cluster_human, 10)
#> Clusters with number of cells below 10 in the new partition: 1
Detect the clusters enriched in highly localized genes and then performs on these a sub-cluster analysis (using louvain algorithm implemented in Seurat)
result_test <- test_hvg(raw_counts_human_data, final_cluster_human, (ciara_genes), background, 100, 0.05)
raw_endoderm <- raw_counts_human_data[, as.vector(umap_elmir$cluster_id) == "Endoderm"]
raw_hemo <- raw_counts_human_data[, as.vector(umap_elmir$cluster_id) == "Hemogenic Endothelial Progenitors"]
raw_exe_meso <- raw_counts_human_data[, as.vector(umap_elmir$cluster_id) == "ExE Mesoderm"]
combined_endoderm <- cluster_analysis_sub(raw_endoderm, 0.2, 5, 30, "Endoderm")
combined_hemo <- cluster_analysis_sub(raw_hemo, 0.6, 5, 30, "Hemogenic Endothelial Progenitors")
combined_exe_meso <- cluster_analysis_sub(raw_exe_meso, 0.5, 5, 30, "ExE Mesoderm")
all_sub_cluster <- c(combined_endoderm$seurat_clusters, combined_hemo$seurat_clusters, combined_exe_meso$seurat_clusters)
final_cluster_human_version_sub <- merge_cluster(final_cluster_human, all_sub_cluster)
table(as.vector(final_cluster_human_version_sub))
#>
#> Advanced Mesoderm Axial Mesoderm
#> 164 23
#> Emergent Mesoderm Endoderm_0
#> 185 79
#> Endoderm_1 Endoderm_2
#> 45 11
#> Epiblast Erythroblasts
#> 133 32
#> ExE Mesoderm_0 ExE Mesoderm_1
#> 46 37
#> Hemogenic Endothelial Progenitors_0 Hemogenic Endothelial Progenitors_1
#> 33 27
#> Hemogenic Endothelial Progenitors_2 Hemogenic Endothelial Progenitors_3
#> 23 15
#> Hemogenic Endothelial Progenitors_4 Nascent Mesoderm
#> 13 98
#> Non-Neural Ectoderm PGC
#> 29 7
#> Primitive Streak
#> 195
Seurat::DefaultAssay(human_data_seurat) <- "RNA"
markers_human_final <- markers_cluster_seurat(human_data_seurat, final_cluster_human_version_sub, names(human_data_seurat$RNA_snn_res.0.1), 5)
markers_human_top_final <- markers_human_final[[1]]
markers_human_all_final <- markers_human_final[[3]]
white_black_markers <- white_black_markers(final_cluster_human_version_sub, "Hemogenic Endothelial Progenitors_4", norm_human_data, markers_human_all_final, 0)
sum(white_black_markers)
white_black_markers <- white_black_markers(final_cluster_human_version_sub, "Endoderm_2", norm_human_data, markers_human_all_final, 0)
sum(white_black_markers)
white_black_markers <- white_black_markers(final_cluster_human_version_sub, "ExE Mesoderm_0", norm_human_data, markers_human_all_final, 0)
sum(white_black_markers)
top_endo <- white_black_markers(final_cluster_human_version_sub, "Endoderm_2", norm_human_data, markers_human_all_final, 0)
top_endo <- names(top_endo)[top_endo]
mean_top_endo <- apply(norm_human_data[top_endo, final_cluster_human_version_sub == "Endoderm_2"], 1, mean)
mean_top_endo <- sort(mean_top_endo, decreasing = T)
top_endo <- names(mean_top_endo)
names(top_endo) <- rep("Endoderm_2", length(top_endo))
top_hemo <- white_black_markers(final_cluster_human_version_sub, "Hemogenic Endothelial Progenitors_4", norm_human_data, markers_human_all_final, 0)
top_hemo <- names(top_hemo)[top_hemo]
mean_top_hemo <- apply(norm_human_data[top_hemo, final_cluster_human_version_sub == "Hemogenic Endothelial Progenitors_4"], 1, mean)
mean_top_hemo <- sort(mean_top_hemo, decreasing = T)
top_hemo <- names(mean_top_hemo)
names(top_hemo) <- rep("Hemogenic Endothelial Progenitors_4", length(top_hemo))
top_meso <- white_black_markers(final_cluster_human_version_sub, "ExE Mesoderm_1", norm_human_data, markers_human_all_final, 0)
top_meso <- names(top_meso)[top_meso]
mean_top_meso <- apply(norm_human_data[top_meso, final_cluster_human_version_sub == "ExE Mesoderm_1"], 1, mean)
mean_top_meso <- sort(mean_top_meso, decreasing = T)
top_meso <- names(mean_top_meso)
names(top_meso) <- rep("ExE Mesoderm_1", length(top_meso))
load(system.file("extdata", "norm_human_data_plot.Rda", package = "CIARA"))
load(system.file("extdata", "top_meso.Rda", package = "CIARA"))
load(system.file("extdata", "top_endo.Rda", package = "CIARA"))
load(system.file("extdata", "top_hemo.Rda", package = "CIARA"))
toMatch <- c("Endoderm")
plot_balloon_marker(norm_human_data_plot[, grep(paste(toMatch, collapse="|"), final_cluster_human)], final_cluster_human_version_sub[grep(paste(toMatch, collapse="|"), final_cluster_human)], top_endo, 20, max_size=5, text_size=10)
toMatch <- c("Hemogenic Endothelial Progenitors")
plot_balloon_marker(norm_human_data_plot[, grep(paste(toMatch, collapse = "|"), final_cluster_human)], final_cluster_human_version_sub[grep(paste(toMatch, collapse = "|"), final_cluster_human)], top_hemo, 20, max_size = 5, text_size = 8)
toMatch <- c("ExE Mesoderm")
plot_balloon_marker(norm_human_data_plot[, grep(paste(toMatch, collapse = "|"), final_cluster_human)], final_cluster_human_version_sub[grep(paste(toMatch, collapse = "|"), final_cluster_human)], top_meso, length(top_meso), max_size = 5, text_size = 8)
Expression of some of the highly localized genes detected by CIARA that are markers of the three rare populations of cells.
utils::sessionInfo()
#> R version 4.0.2 (2020-06-22)
#> Platform: x86_64-apple-darwin17.0 (64-bit)
#> Running under: macOS Mojave 10.14.6
#>
#> Matrix products: default
#> BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
#>
#> locale:
#> [1] C/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] CIARA_0.1.0
#>
#> loaded via a namespace (and not attached):
#> [1] Seurat_4.0.5 Rtsne_0.15 colorspace_2.0-2
#> [4] deldir_1.0-6 ellipsis_0.3.2 ggridges_0.5.3
#> [7] spatstat.data_2.1-0 leiden_0.3.9 listenv_0.8.0
#> [10] farver_2.1.0 graphlayouts_0.8.0 ggrepel_0.9.1
#> [13] fansi_0.5.0 codetools_0.2-16 splines_4.0.2
#> [16] knitr_1.37 polyclip_1.10-0 jsonlite_1.7.2
#> [19] ica_1.0-2 cluster_2.1.0 png_0.1-7
#> [22] uwot_0.1.11 ggforce_0.3.3 shiny_1.7.1
#> [25] sctransform_0.3.2 spatstat.sparse_2.0-0 compiler_4.0.2
#> [28] httr_1.4.2 assertthat_0.2.1 SeuratObject_4.0.4
#> [31] Matrix_1.3-4 fastmap_1.1.0 lazyeval_0.2.2
#> [34] later_1.3.0 tweenr_1.0.2 htmltools_0.5.2
#> [37] tools_4.0.2 igraph_1.2.9 gtable_0.3.0
#> [40] glue_1.6.0 RANN_2.6.1 reshape2_1.4.4
#> [43] dplyr_1.0.7 Rcpp_1.0.7 scattermore_0.7
#> [46] jquerylib_0.1.4 vctrs_0.3.8 nlme_3.1-149
#> [49] crosstalk_1.2.0 ggraph_2.0.5 lmtest_0.9-39
#> [52] xfun_0.29 stringr_1.4.0 globals_0.14.0
#> [55] mime_0.12 miniUI_0.1.1.1 lifecycle_1.0.1
#> [58] irlba_2.3.5 goftest_1.2-3 future_1.23.0
#> [61] MASS_7.3-52 zoo_1.8-9 scales_1.1.1
#> [64] tidygraph_1.2.0 spatstat.core_2.3-2 promises_1.2.0.1
#> [67] spatstat.utils_2.2-0 parallel_4.0.2 RColorBrewer_1.1-2
#> [70] yaml_2.2.1 reticulate_1.22 pbapply_1.5-0
#> [73] gridExtra_2.3 ggplot2_3.3.5 sass_0.4.0
#> [76] rpart_4.1-15 stringi_1.7.6 highr_0.9
#> [79] rlang_0.4.12 pkgconfig_2.0.3 matrixStats_0.61.0
#> [82] evaluate_0.14 lattice_0.20-41 ROCR_1.0-11
#> [85] purrr_0.3.4 tensor_1.5 labeling_0.4.2
#> [88] patchwork_1.1.1 htmlwidgets_1.5.4 cowplot_1.1.1
#> [91] tidyselect_1.1.1 parallelly_1.29.0 RcppAnnoy_0.0.19
#> [94] plyr_1.8.6 magrittr_2.0.1 R6_2.5.1
#> [97] generics_0.1.1 DBI_1.1.2 mgcv_1.8-32
#> [100] pillar_1.6.4 fitdistrplus_1.1-6 survival_3.2-3
#> [103] abind_1.4-5 tibble_3.1.6 future.apply_1.8.1
#> [106] crayon_1.4.2 KernSmooth_2.23-17 utf8_1.2.2
#> [109] spatstat.geom_2.3-0 plotly_4.10.0 rmarkdown_2.11
#> [112] viridis_0.6.2 grid_4.0.2 data.table_1.14.2
#> [115] digest_0.6.29 xtable_1.8-4 tidyr_1.1.4
#> [118] httpuv_1.6.3 munsell_0.5.0 viridisLite_0.4.0
#> [121] bslib_0.3.1
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.