library(scater)
library(ggplot2)
library(BiocParallel)
library(blase)
library(patchwork)
RNGversion("3.5.0")
#> Warning in RNGkind("Mersenne-Twister", "Inversion", "Rounding"): non-uniform
#> 'Rounding' sampler used
SEED <- 7
set.seed(SEED)
N_CORES <- 2
bpparam <- MulticoreParam(N_CORES)
This article will show how BLASE can be used for annotating scRNA-seq data using existing bulk or microarray data. We make use of scRNA-seq Dogga et al. 2024 and microarray data from Painter et al. 2018. Code for generating the objects used here is available in the BLASE reproducibility repository.
First we’ll load in the data we’re using, pre-prepared from the BLASE reproducibility repository.
data(painter_microarray, package = "blase")
data(MCA_PF_SCE, package = "blase")
We can examine the true lifecycle stages, and also the calculated pseudotime trajectory (Slingshot).
plotUMAP(MCA_PF_SCE, colour_by = "STAGE_HR")
#| fig.alt: >
#| UMAP of Dogga et al. single-cell RNA-seq reference coloured by pseudotime,
#| starting with Rings, and ending with Schizonts.
plotUMAP(MCA_PF_SCE, colour_by = "slingPseudotime_1")
Now we’ll prepare BLASE for use.
First, we create the object, giving it the number of bins we want to use, and how to calculate them.
blase_data <- as.BlaseData(
MCA_PF_SCE,
pseudotime_slot = "slingPseudotime_1",
n_bins = 48,
split_by = "cells"
)
# Add these bins to the sc metadata
MCA_PF_SCE <- assign_pseudotime_bins(MCA_PF_SCE,
pseudotime_slot = "slingPseudotime_1",
n_bins = 48,
split_by = "cells"
)
Now we will select the genes we want to use, using BLASE’s gene peakedness metric.
gene_peakedness_info <- calculate_gene_peakedness(
MCA_PF_SCE,
window_pct = 5,
knots = 18,
BPPARAM = bpparam
)
genes_to_use <- gene_peakedness_spread_selection(
MCA_PF_SCE,
gene_peakedness_info,
genes_per_bin = 30,
n_gene_bins = 30
)
head(gene_peakedness_info)
#> gene peak_pseudotime mean_in_window mean_out_window ratio
#> 28 PF3D7-1401100 3.8311312 5.793085 1.528934 3.788969
#> 44 PF3D7-1401200 6.0203490 4.791942 1.921439 2.493934
#> 29 PF3D7-1401400 3.9679573 696.174509 39.689901 17.540344
#> 84 PF3D7-1401600 11.4933936 93.076303 7.631160 12.196875
#> 1 PF3D7-1401900 0.1368261 43.270750 3.321464 13.027613
#> 80 PF3D7-1402200 10.9460892 3.754015 1.554102 2.415553
#> window_start window_end deviance_explained
#> 28 3.4890659 4.1731965 0.1881312
#> 44 5.6782838 6.3624143 0.2035355
#> 29 3.6258920 4.3100226 0.2513795
#> 84 11.1513283 11.8354589 0.6891177
#> 1 -0.2052392 0.4788914 0.3846938
#> 80 10.6040239 11.2881544 0.1768597
By using the gene_peakedness_spread_selection
function, we can ensure
that genes with high ratios are selected from throughout the trajectory.
ggplot(gene_peakedness_info, aes(x = peak_pseudotime, y = ratio)) +
geom_point() +
ggtitle("All genes")
gene_peakedness_selected_genes <- gene_peakedness_info[
gene_peakedness_info$gene %in% genes_to_use,
]
ggplot(gene_peakedness_selected_genes, aes(x = peak_pseudotime, y = ratio)) +
geom_point() +
ggtitle("Selected genes")
Here, we add them to the BLASE object for mapping with.
genes(blase_data) <- genes_to_use
Now we can perform the actual mapping step, and review the results.
mapping_results <- map_all_best_bins(
blase_data,
painter_microarray,
BPPARAM = bpparam
)
plot_mapping_result_heatmap(mapping_results)
In the Painter et al. paper, the expected lifecycle stages are as follows:
Lifecycle Stage | Painter et al. HPI |
---|---|
Ring | 0-21 |
Trophozoite | 16-32 |
Schizont | 33-48 |
We can use this to transfer back to the scRNA-seq as follows:
plotUMAP(MCA_PF_SCE, colour_by = "pseudotime_bin")
annotated_sce <- annotate_sce(MCA_PF_SCE, mapping_results, include_stats = TRUE)
annotated_sce$BLASE_Annotation <- gsub(" hpi", "", annotated_sce$BLASE_Annotation)
annotated_sce$BLASE_Annotation <- as.numeric(annotated_sce$BLASE_Annotation)
annotation_plot <- plotUMAP(annotated_sce, colour_by = "BLASE_Annotation")
corr_plot <- plotUMAP(annotated_sce, colour_by = "BLASE_Annotation_Correlation") +
labs(color = "Correlation")
annotation_plot + corr_plot
annotated_sce$BLASE_Annotation_transfer <- ""
for (best_bulk in unique(annotated_sce$BLASE_Annotation)) {
mask <- colData(annotated_sce)[, "BLASE_Annotation"] == best_bulk
if (best_bulk %in% 0:15) {
annotated_sce[, mask]$BLASE_Annotation_transfer <- "Ring"
} else if (best_bulk %in% 16:21) {
annotated_sce[, mask]$BLASE_Annotation_transfer <- "Ring/Trophozoite"
} else if (best_bulk %in% 22:32) {
annotated_sce[, mask]$BLASE_Annotation_transfer <- "Trophozoite"
} else if (best_bulk %in% 33:48) {
annotated_sce[, mask]$BLASE_Annotation_transfer <- "Schizont"
}
}
plotUMAP(annotated_sce, colour_by = "BLASE_Annotation_transfer") +
plotUMAP(annotated_sce, colour_by = "STAGE_LR")
sessionInfo()
#> R version 4.5.1 Patched (2025-08-23 r88802)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.3 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.22-bioc/R/lib/libRblas.so
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0
#>
#> Random number generation:
#> RNG: Mersenne-Twister
#> Normal: Inversion
#> Sample: Rounding
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_GB LC_COLLATE=C
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: America/New_York
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats4 stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] patchwork_1.3.2 blase_0.99.3
#> [3] BiocParallel_1.43.4 scater_1.37.0
#> [5] ggplot2_4.0.0 scuttle_1.19.0
#> [7] SingleCellExperiment_1.31.1 SummarizedExperiment_1.39.2
#> [9] Biobase_2.69.1 GenomicRanges_1.61.5
#> [11] Seqinfo_0.99.2 IRanges_2.43.5
#> [13] S4Vectors_0.47.4 BiocGenerics_0.55.1
#> [15] generics_0.1.4 MatrixGenerics_1.21.0
#> [17] matrixStats_1.5.0 BiocStyle_2.37.1
#>
#> loaded via a namespace (and not attached):
#> [1] RColorBrewer_1.1-3 jsonlite_2.0.0 magrittr_2.0.4
#> [4] magick_2.9.0 spatstat.utils_3.2-0 ggbeeswarm_0.7.2
#> [7] farver_2.1.2 rmarkdown_2.30 vctrs_0.6.5
#> [10] ROCR_1.0-11 spatstat.explore_3.5-3 tinytex_0.57
#> [13] htmltools_0.5.8.1 S4Arrays_1.9.1 BiocNeighbors_2.3.1
#> [16] SparseArray_1.9.1 sass_0.4.10 sctransform_0.4.2
#> [19] parallelly_1.45.1 KernSmooth_2.23-26 bslib_0.9.0
#> [22] htmlwidgets_1.6.4 ica_1.0-3 plyr_1.8.9
#> [25] plotly_4.11.0 zoo_1.8-14 cachem_1.1.0
#> [28] igraph_2.1.4 mime_0.13 lifecycle_1.0.4
#> [31] pkgconfig_2.0.3 rsvd_1.0.5 Matrix_1.7-4
#> [34] R6_2.6.1 fastmap_1.2.0 fitdistrplus_1.2-4
#> [37] future_1.67.0 shiny_1.11.1 digest_0.6.37
#> [40] tensor_1.5.1 Seurat_5.3.0 RSpectra_0.16-2
#> [43] irlba_2.3.5.1 beachmat_2.25.5 labeling_0.4.3
#> [46] progressr_0.16.0 spatstat.sparse_3.1-0 mgcv_1.9-3
#> [49] polyclip_1.10-7 httr_1.4.7 abind_1.4-8
#> [52] compiler_4.5.1 withr_3.0.2 S7_0.2.0
#> [55] viridis_0.6.5 fastDummies_1.7.5 MASS_7.3-65
#> [58] DelayedArray_0.35.3 tools_4.5.1 vipor_0.4.7
#> [61] lmtest_0.9-40 beeswarm_0.4.0 httpuv_1.6.16
#> [64] future.apply_1.20.0 goftest_1.2-3 glue_1.8.0
#> [67] nlme_3.1-168 promises_1.3.3 grid_4.5.1
#> [70] Rtsne_0.17 cluster_2.1.8.1 reshape2_1.4.4
#> [73] spatstat.data_3.1-8 gtable_0.3.6 tidyr_1.3.1
#> [76] data.table_1.17.8 BiocSingular_1.25.0 ScaledMatrix_1.17.0
#> [79] sp_2.2-0 XVector_0.49.1 spatstat.geom_3.6-0
#> [82] RcppAnnoy_0.0.22 stringr_1.5.2 ggrepel_0.9.6
#> [85] RANN_2.6.2 pillar_1.11.1 spam_2.11-1
#> [88] RcppHNSW_0.6.0 later_1.4.4 splines_4.5.1
#> [91] dplyr_1.1.4 lattice_0.22-7 deldir_2.0-4
#> [94] survival_3.8-3 tidyselect_1.2.1 miniUI_0.1.2
#> [97] pbapply_1.7-4 knitr_1.50 gridExtra_2.3
#> [100] bookdown_0.45 scattermore_1.2 xfun_0.53
#> [103] stringi_1.8.7 boot_1.3-32 lazyeval_0.2.2
#> [106] yaml_2.3.10 evaluate_1.0.5 codetools_0.2-20
#> [109] tibble_3.3.0 BiocManager_1.30.26 cli_3.6.5
#> [112] uwot_0.2.3 xtable_1.8-4 reticulate_1.43.0
#> [115] jquerylib_0.1.4 dichromat_2.0-0.1 Rcpp_1.1.0
#> [118] spatstat.random_3.4-2 globals_0.18.0 png_0.1-8
#> [121] spatstat.univar_3.1-4 parallel_4.5.1 dotCall64_1.2
#> [124] sparseMatrixStats_1.21.0 listenv_0.9.1 viridisLite_0.4.2
#> [127] scales_1.4.0 ggridges_0.5.7 SeuratObject_5.2.0
#> [130] purrr_1.1.0 crayon_1.5.3 rlang_1.1.6
#> [133] cowplot_1.2.0