Downstream analysis
In the following, we visualize the CAESAR annotation results.
# set up colors
cols_manual <- setNames(
c(
"#4374A5", "#FCDDDE", "#2AB67F", "#F08A21", "#737373"
),
c(
"GCL", "MCL", "ONL", "GL", "Unknown"
)
)
celltypes_manual <- c("GCL", "MCL", "ONL", "GL", "Unknown")
cols <- setNames(
c(
"#4374A5", "#FCDDDE", "#2AB673", "#F08A21", "#E04D50", "#737373"
),
c(
"GC", "M/TC", "OSNs", "PGC", "EPL-IN", "unassigned"
)
)
celltypes <- c("GC", "M/TC", "OSNs", "PGC", "EPL-IN", "unassigned")
colnames(pos) <- paste0("pos", 1:2)
MOB_ST@reductions[["pos"]] <- CreateDimReducObject(
embeddings = as.matrix(pos),
key = paste0("pos", "_"), assay = "RNA"
)First, we visualize the manual annotation.
Idents(MOB_ST) <- factor(MOB_ST$manual_annotation, levels = celltypes_manual)
DimPlot(MOB_ST, reduction = "pos", cols = cols_manual, pt.size = 8)
Then, we visualize the CAESAR annotation without account for ‘unassigned’.
Idents(MOB_ST) <- factor(MOB_ST$CAESAR, levels = celltypes)
DimPlot(MOB_ST, reduction = "pos", cols = cols, pt.size = 8)
And visualize the CAESAR annotation account for ‘unassigned’.
Idents(MOB_ST) <- factor(MOB_ST$CAESARunasg, levels = celltypes)
DimPlot(MOB_ST, reduction = "pos", cols = cols, pt.size = 8)
The confidence level of the CAESAR annotation can be visualized by
FeaturePlot(
MOB_ST,
reduction = "pos", features = "CAESARconf", pt.size = 8,
cols = c("blue", "lightgrey"), min.cutoff = 0.0, max.cutoff = 1.0
)
CAESAR provides the cell mixing proportion for each cell type, which can be visualized by
caesar_prob <- colnames(MOB_ST@meta.data)[15:19]
print(caesar_prob)
#> [1] "prob_EPL.IN" "prob_GC" "prob_M.TC" "prob_OSNs" "prob_PGC"
plots <- lapply(caesar_prob, function(feature) {
FeaturePlot(MOB_ST, features = feature, reduction = "pos", pt.size = 3.5) +
scale_color_gradientn(
colors = c("#f6eff7", "#feebe2", "#f768a1", "#7a0177", "#6e016b"),
values = scales::rescale(c(0.0, 0.125, 0.25, 0.375, 0.50)),
limits = c(0.0, 0.50)
) + labs(title = feature)
})
#> Scale for colour is already present.
#> Adding another scale for colour, which will replace the existing scale.
#> Scale for colour is already present.
#> Adding another scale for colour, which will replace the existing scale.
#> Scale for colour is already present.
#> Adding another scale for colour, which will replace the existing scale.
#> Scale for colour is already present.
#> Adding another scale for colour, which will replace the existing scale.
#> Scale for colour is already present.
#> Adding another scale for colour, which will replace the existing scale.
cowplot::plot_grid(plotlist = plots, ncol = 2)
The annotation accuracy is calculated by
acc_st <- function(manual_annotation, pred) {
manual_annotation <- as.character(manual_annotation)
pred <- as.character(pred)
manual_annotation[manual_annotation == "GCL"] <- "GC"
manual_annotation[manual_annotation == "MCL"] <- "M/TC"
manual_annotation[manual_annotation == "ONL"] <- "OSNs"
manual_annotation[manual_annotation == "GL"] <- "PGC"
return(mean(manual_annotation == pred))
}
print(paste0(
"The ACC of CAESAR annotation is ",
acc_st(MOB_ST$manual_annotation, MOB_ST$CAESARunasg)
))
#> [1] "The ACC of CAESAR annotation is 0.717857142857143"Next, we detect and visualize the signature genes for each cell type.
Idents(MOB_ST) <- factor(MOB_ST$CAESARunasg, celltypes)
sg_List <- find.sig.genes(MOB_ST)
str(sg_List)
#> List of 4
#> $ GC :'data.frame': 13895 obs. of 5 variables:
#> ..$ distance : num [1:13895] 5.49 5.49 5.52 5.53 5.56 ...
#> ..$ expr.prop : num [1:13895] 0.344 0.422 0.322 0.322 0.256 ...
#> ..$ expr.prop.others: num [1:13895] 0.0842 0.0947 0.0579 0.0421 0.0421 ...
#> ..$ label : chr [1:13895] "GC" "GC" "GC" "GC" ...
#> ..$ gene : chr [1:13895] "Cdkn1a" "Rasgrp2" "Kcnh5" "Kcnh1" ...
#> $ M/TC:'data.frame': 13895 obs. of 5 variables:
#> ..$ distance : num [1:13895] 6.69 6.74 6.74 6.75 6.75 ...
#> ..$ expr.prop : num [1:13895] 0.3494 0.0964 0.6627 0.7711 0.3012 ...
#> ..$ expr.prop.others: num [1:13895] 0.1066 0.0152 0.2792 0.3147 0.1015 ...
#> ..$ label : chr [1:13895] "M/TC" "M/TC" "M/TC" "M/TC" ...
#> ..$ gene : chr [1:13895] "Raver2" "Slc16a10" "Zfp804a" "Nefl" ...
#> $ OSNs:'data.frame': 13895 obs. of 5 variables:
#> ..$ distance : num [1:13895] 9.69 9.72 9.73 9.73 9.73 ...
#> ..$ expr.prop : num [1:13895] 0.836 0.448 0.612 0.701 0.493 ...
#> ..$ expr.prop.others: num [1:13895] 0.2582 0.0892 0.1643 0.1315 0.1174 ...
#> ..$ label : chr [1:13895] "OSNs" "OSNs" "OSNs" "OSNs" ...
#> ..$ gene : chr [1:13895] "Hmgcs2" "Alx4" "Afap1l2" "Igfbp4" ...
#> $ PGC :'data.frame': 13895 obs. of 5 variables:
#> ..$ distance : num [1:13895] 7.71 7.74 7.74 7.77 7.77 ...
#> ..$ expr.prop : num [1:13895] 0.55 0.825 0.6 0.35 0.5 0.075 0.375 0.2 0.1 0.85 ...
#> ..$ expr.prop.others: num [1:13895] 0.183 0.412 0.267 0.188 0.192 ...
#> ..$ label : chr [1:13895] "PGC" "PGC" "PGC" "PGC" ...
#> ..$ gene : chr [1:13895] "Cnr1" "AI593442" "Shisa6" "Kcnj12" ...We visualize the top three signature genes by a dot plot.
# obtain the top three signature genes
celltypes_plot <- setdiff(names(sg_List), "unassigned")
top3sgs <- Intsg(list(sg_List), 3)[celltypes_plot]
print(top3sgs)
#> $GC
#> [1] "Cdkn1a" "Rasgrp2" "Kcnh5"
#>
#> $`M/TC`
#> [1] "Raver2" "Zfp804a" "Nefl"
#>
#> $OSNs
#> [1] "Hmgcs2" "Alx4" "Afap1l2"
#>
#> $PGC
#> [1] "Cnr1" "AI593442" "Shisa6"
sg_features <- unname(unlist(top3sgs))
DotPlot(
MOB_ST,
idents = celltypes_plot, col.min = -1, col.max = 2, dot.scale = 7,
features = sg_features, scale.min = 0, scale.max = 30
) + theme(axis.text.x = element_text(face = "italic", angle = 45, vjust = 1, hjust = 1))
#> Warning: Scaling data with a low number of groups may produce misleading
#> results
Next, we calculate the UMAP projections of co-embeddings of cells and the selected signature genes.
# calculate coumap
MOB_ST <- CoUMAP(
MOB_ST, reduction = "caesar", reduction.name = "caesarUMAP",
gene.set = sg_features
)
df_gene_label <- data.frame(
gene = unlist(top3sgs),
label = rep(names(top3sgs), each = 3)
)
CoUMAP.plot(
MOB_ST, reduction = "caesarUMAP", gene_txtdata = df_gene_label,
cols = c("gene" = "#000000", cols)
)
Session Info
sessionInfo()
#> R version 4.3.1 (2023-06-16 ucrt)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 11 x64 (build 22631)
#>
#> Matrix products: default
#>
#>
#> locale:
#> [1] LC_COLLATE=C
#> [2] LC_CTYPE=Chinese (Simplified)_China.utf8
#> [3] LC_MONETARY=Chinese (Simplified)_China.utf8
#> [4] LC_NUMERIC=C
#> [5] LC_TIME=Chinese (Simplified)_China.utf8
#>
#> time zone: Asia/Shanghai
#> tzcode source: internal
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] purrr_1.0.2 ggplot2_3.5.1 ProFAST_1.4 gtools_3.9.5
#> [5] Seurat_4.4.0 CAESAR.Suite_0.1.0 SeuratObject_4.1.4 sp_2.1-4
#>
#> loaded via a namespace (and not attached):
#> [1] matrixStats_1.3.0 spatstat.sparse_3.1-0
#> [3] bitops_1.0-8 httr_1.4.7
#> [5] RColorBrewer_1.1-3 tools_4.3.1
#> [7] sctransform_0.4.1 backports_1.5.0
#> [9] utf8_1.2.4 R6_2.5.1
#> [11] lazyeval_0.2.2 uwot_0.2.2
#> [13] withr_3.0.1 prettyunits_1.2.0
#> [15] gridExtra_2.3 progressr_0.14.0
#> [17] cli_3.6.3 Biobase_2.62.0
#> [19] spatstat.explore_3.3-2 labeling_0.4.3
#> [21] sass_0.4.9 mvtnorm_1.2-6
#> [23] spatstat.data_3.1-2 proxy_0.4-27
#> [25] ggridges_0.5.6 pbapply_1.7-2
#> [27] harmony_1.2.1 scater_1.30.1
#> [29] parallelly_1.38.0 readxl_1.4.3
#> [31] rstudioapi_0.16.0 RSQLite_2.3.7
#> [33] FNN_1.1.4 generics_0.1.3
#> [35] ica_1.0-3 spatstat.random_3.3-1
#> [37] car_3.1-2 dplyr_1.1.4
#> [39] Matrix_1.6-1.1 ggbeeswarm_0.7.2
#> [41] fansi_1.0.6 DescTools_0.99.56
#> [43] S4Vectors_0.40.2 abind_1.4-5
#> [45] lifecycle_1.0.4 yaml_2.3.10
#> [47] CompQuadForm_1.4.3 carData_3.0-5
#> [49] SummarizedExperiment_1.32.0 SparseArray_1.2.4
#> [51] BiocFileCache_2.10.2 Rtsne_0.17
#> [53] grid_4.3.1 blob_1.2.4
#> [55] promises_1.3.0 crayon_1.5.3
#> [57] GiRaF_1.0.1 miniUI_0.1.1.1
#> [59] lattice_0.21-8 beachmat_2.18.1
#> [61] cowplot_1.1.3 KEGGREST_1.42.0
#> [63] pillar_1.9.0 knitr_1.48
#> [65] GenomicRanges_1.54.1 boot_1.3-28.1
#> [67] gld_2.6.6 future.apply_1.11.2
#> [69] codetools_0.2-19 leiden_0.4.3.1
#> [71] glue_1.7.0 spatstat.univar_3.0-0
#> [73] data.table_1.16.0 vctrs_0.6.5
#> [75] png_0.1-8 spam_2.10-0
#> [77] org.Mm.eg.db_3.18.0 cellranger_1.1.0
#> [79] gtable_0.3.5 cachem_1.1.0
#> [81] xfun_0.47 S4Arrays_1.2.1
#> [83] mime_0.12 survival_3.5-5
#> [85] SingleCellExperiment_1.24.0 fitdistrplus_1.2-1
#> [87] ROCR_1.0-11 nlme_3.1-162
#> [89] bit64_4.0.5 progress_1.2.3
#> [91] filelock_1.0.3 RcppAnnoy_0.0.22
#> [93] GenomeInfoDb_1.38.8 bslib_0.8.0
#> [95] irlba_2.3.5.1 vipor_0.4.7
#> [97] KernSmooth_2.23-21 colorspace_2.1-1
#> [99] BiocGenerics_0.48.1 DBI_1.2.3
#> [101] ade4_1.7-22 Exact_3.3
#> [103] tidyselect_1.2.1 DR.SC_3.4
#> [105] bit_4.0.5 compiler_4.3.1
#> [107] curl_5.2.2 BiocNeighbors_1.20.2
#> [109] expm_1.0-0 xml2_1.3.6
#> [111] DelayedArray_0.28.0 plotly_4.10.4
#> [113] scales_1.3.0 lmtest_0.9-40
#> [115] rappdirs_0.3.3 stringr_1.5.1
#> [117] digest_0.6.37 goftest_1.2-3
#> [119] spatstat.utils_3.1-0 rmarkdown_2.28
#> [121] XVector_0.42.0 htmltools_0.5.8.1
#> [123] pkgconfig_2.0.3 sparseMatrixStats_1.14.0
#> [125] MatrixGenerics_1.14.0 highr_0.11
#> [127] dbplyr_2.5.0 fastmap_1.2.0
#> [129] rlang_1.1.4 htmlwidgets_1.6.4
#> [131] ggthemes_5.1.0 shiny_1.9.1
#> [133] DelayedMatrixStats_1.24.0 farver_2.1.2
#> [135] jquerylib_0.1.4 zoo_1.8-12
#> [137] jsonlite_1.8.8 BiocParallel_1.36.0
#> [139] mclust_6.1.1 BiocSingular_1.18.0
#> [141] RCurl_1.98-1.16 magrittr_2.0.3
#> [143] scuttle_1.12.0 GenomeInfoDbData_1.2.11
#> [145] dotCall64_1.1-1 patchwork_1.2.0
#> [147] munsell_0.5.1 Rcpp_1.0.13
#> [149] viridis_0.6.5 reticulate_1.38.0
#> [151] furrr_0.3.1 stringi_1.8.4
#> [153] rootSolve_1.8.2.4 zlibbioc_1.48.2
#> [155] MASS_7.3-60 org.Hs.eg.db_3.18.0
#> [157] plyr_1.8.9 parallel_4.3.1
#> [159] PRECAST_1.6.5 listenv_0.9.1
#> [161] ggrepel_0.9.5 lmom_3.0
#> [163] deldir_2.0-4 Biostrings_2.70.3
#> [165] splines_4.3.1 tensor_1.5
#> [167] hms_1.1.3 igraph_2.0.3
#> [169] ggpubr_0.6.0 spatstat.geom_3.3-2
#> [171] ggsignif_0.6.4 reshape2_1.4.4
#> [173] biomaRt_2.58.2 stats4_4.3.1
#> [175] ScaledMatrix_1.10.0 XML_3.99-0.17
#> [177] evaluate_0.24.0 httpuv_1.6.15
#> [179] RANN_2.6.2 tidyr_1.3.1
#> [181] polyclip_1.10-7 future_1.34.0
#> [183] scattermore_1.2 rsvd_1.0.5
#> [185] broom_1.0.6 xtable_1.8-4
#> [187] e1071_1.7-14 rstatix_0.7.2
#> [189] later_1.3.2 viridisLite_0.4.2
#> [191] class_7.3-22 tibble_3.2.1
#> [193] memoise_2.0.1 beeswarm_0.4.0
#> [195] AnnotationDbi_1.64.1 IRanges_2.36.0
#> [197] cluster_2.1.4 globals_0.16.3