library(blase)
library(SingleCellExperiment)
library(tradeSeq)
library(slingshot)
library(scran)
library(scater)
library(BiocParallel)
library(ami)
library(patchwork)
library(reshape2)
RNGversion("3.5.0")
#> Warning in RNGkind("Mersenne-Twister", "Inversion", "Rounding"): non-uniform
#> 'Rounding' sampler used
SEED <- 7
set.seed(SEED)
N_CORES <- 4
if (ami::using_ci()) {
N_CORES <- 2
}
BLASE is a tool for matching bulk RNA-seq samples to the pseudotime of a
trajectory in a scRNA-seq dataset. BLASE integrates with the Bioconductor
ecosystem by integrating closely with SingleCellExperiment objects. (BLASE
can be used with Seurat objects that have been converted using
as.SingleCellExperiment
.)
BLASE can be used for the following use-cases: * Identifying the progress of a bulk like RNA-seq sample through scRNA-seq reference trajectory. * Annotating scRNA-seq trajectories with reference bulk RNA-seq datasets * Correcting development rate differences in comparative RNA-seq samples.
We believe that in many trajectories that cells undergo, there will be a curve that shows the expression of some key genes. In the plot below, we can see that by looking at these three “genes” which have 2 peaks each over “pseudotime”, we can identify 6 unique states, which, as long as the genes correspond to some aspect of the trajectory, can be used to infer which stage of the trajectory a cell is currently in. This is of course a constructed example and unlikely to be exactly repeated in real life, but hopefully demonstrates the principle BLASE works on.
# Adapted from
# https://codepal.ai/code-generator/query/n4dEA6I9/plot-sine-wave-ggplot2-r
x <- seq(0, 2 * pi, length.out = 500)
gene1_expr <- 10 * (sin(x + 0.1) + 1)
gene2_expr <- 11 * (sin(x - pi) + 1)
gene3_expr <- 7 * (sin(x + pi / 2) + 1)
genes <- data.frame(gene1 = gene1_expr, gene2 = gene2_expr, gene3 = gene3_expr)
genes_melt <- melt(t(genes), varnames = c("gene", "x"))
genes_melt$x <- as.numeric(genes_melt$x)
ggplot(genes_melt, aes = aes()) +
geom_line(aes(x = x, y = value, color = gene)) +
geom_vline(xintercept = 0, linetype = "dashed") +
geom_vline(xintercept = 24, linetype = "dashed") +
geom_vline(xintercept = 181, linetype = "dashed") +
geom_vline(xintercept = 243, linetype = "dashed") +
geom_vline(xintercept = 315, linetype = "dashed") +
geom_vline(xintercept = 479, linetype = "dashed") +
geom_vline(xintercept = 500, linetype = "dashed") +
labs(
x = "pseudotime",
y = "Gene Expression",
title = "Variable Genes Over Pseudotime"
)
#> Warning in fortify(data, ...): Arguments in `...` must be used.
#> ✖ Problematic argument:
#> • aes = aes()
#> ℹ Did you misspell an argument name?
In this vignette, we will use BLASE to try to assign bulk RNA-seq samples to the pseudotime trajectory in a scRNA-seq data. We need to have the following to use Blase in this case:
fitGAM
and associationTest
To install BLASE from Bioconductor:
if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install("blase")
First, let’s load a Single Cell Experiment to use the tool with, based on the tradeSeq vignette. We have added pseudotime to the object, and used tradeSeq as shown in their vignette.
data(tradeSeq_BLASE_example_sce, package = "blase")
sce <- tradeSeq_BLASE_example_sce
plotUMAP(sce, text_by = "celltype", colour_by = "celltype") +
plotUMAP(sce, colour_by = "pseudotime") +
patchwork::plot_layout(ncol = 1, axis_title = "collect")
Now we’ll find the genes we want to use. We select the top 200 so that we can do some parameter tuning below.
# Use consecutive for genes that change over time
assoRes <- associationTest(
sce,
lineages = TRUE,
global = FALSE,
contrastType = "consecutive"
)
genelist <- blase::get_top_n_genes(assoRes, n_genes = 200, lineage = 1)
Please note that other methods of selecting genes are available. Please see other vignettes for details.
When using BLASE, it can be a good idea to tune the number of genes and bins used.
We can do this using the following commands, provided by BLASE:
res <- blase::find_best_params(
sce,
genelist,
split_by = "pseudotime_range",
pseudotime_slot = "pseudotime",
bins_count_range = c(5, 10),
gene_count_range = c(40, 80)
)
blase::plot_find_best_params_results(res)
It looks like 80 genes and 10 bins will give us good specificity, but let’s double check. In a non trivial dataset, this may take some repetition. We ignore the 5 bin result because this might not be enough resolution - this depends on your dataset. In general, more bins will reduce specificity because each bin will have a more similar cell composition. In general, aim for about as many bins as you have clusters.
Now we can check if this is a good fit by showing how other bins from the SC dataset map using these genes. Ideally, each bin is very specific, with a high “worst specificity”.
blase_data <- as.BlaseData(sce, pseudotime_slot = "pseudotime", n_bins = 10)
genes(blase_data) <- genelist[1:80]
evaluate_parameters(blase_data, make_plot = TRUE)
#> [1] 0.00140 0.04096 20.00000
We can also see how the genes change over pseudotime, by plotting expression of the top genes changing over each pseudotime bin.
evaluate_top_n_genes(blase_data)
<img src=“/tmp/RtmpYMhEX0/Rbuild3eacdfb8f9d8/blase/vignettes/assign-bulk-to-pseudotime_files/figure-html/evaluateTopGenes-1.png” alt=“Plot showing a the”best" genes selected for use with BLASE, according to how related their expression is to the pseudotime of the cell." width=“100%” />
If we’re happy, now we can try to map a bulk sample onto the single cell. We’ll do this by pseudobulking cell types in the SingleCellExperiment but in reality you should use a real bulk dataset. See some of our other articles for examples.
bulks_df <- DataFrame(row.names = rownames(counts(sce)))
for (type in unique(sce$celltype)) {
bulks_df <- cbind(
bulks_df,
rowSums(normcounts(subset(sce, , celltype == type)))
)
}
colnames(bulks_df) <- unique(sce$celltype)
multipotent_progenitors_result <- map_best_bin(
blase_data, "Multipotent progenitors", bulks_df
)
multipotent_progenitors_result
#> MappingResult for 'Multipotent progenitors': best_bin=1 correlation=0.984458509142053 top_2_distance=0.0799
#> Confident Result: TRUE (next max upper 0.936108081483306 )
#> with history for scores against 10 bins
#> Bootstrapped with 200 iterations
basophils_result <- map_best_bin(blase_data, "Basophils", bulks_df)
basophils_result
#> MappingResult for 'Basophils': best_bin=6 correlation=0.961556493202063 top_2_distance=0.0011
#> Confident Result: FALSE (next max upper 0.975955335569683 )
#> with history for scores against 10 bins
#> Bootstrapped with 200 iterations
neutrophils_result <- map_best_bin(blase_data, "Neutrophils", bulks_df)
neutrophils_result
#> MappingResult for 'Neutrophils': best_bin=10 correlation=0.989006094702297 top_2_distance=0.0764
#> Confident Result: TRUE (next max upper 0.949556183530129 )
#> with history for scores against 10 bins
#> Bootstrapped with 200 iterations
We can also generate this result with map_all_best_bins
, which can
also allow for parallelisation of the process for large numbers of
bulks:
map_all_best_bins(blase_data, bulks_df)
#> [[1]]
#> MappingResult for 'Multipotent progenitors': best_bin=1 correlation=0.984458509142053 top_2_distance=0.0799
#> Confident Result: TRUE (next max upper 0.938084125848705 )
#> with history for scores against 10 bins
#> Bootstrapped with 200 iterations
#>
#> [[2]]
#> MappingResult for 'Monocytes': best_bin=9 correlation=0.979817158931083 top_2_distance=0.0148
#> Confident Result: FALSE (next max upper 0.980156449740228 )
#> with history for scores against 10 bins
#> Bootstrapped with 200 iterations
#>
#> [[3]]
#> MappingResult for 'Neutrophils': best_bin=10 correlation=0.989006094702297 top_2_distance=0.0764
#> Confident Result: TRUE (next max upper 0.943140635445621 )
#> with history for scores against 10 bins
#> Bootstrapped with 200 iterations
#>
#> [[4]]
#> MappingResult for 'Basophils': best_bin=6 correlation=0.961556493202063 top_2_distance=0.0011
#> Confident Result: FALSE (next max upper 0.97610197117696 )
#> with history for scores against 10 bins
#> Bootstrapped with 200 iterations
#>
#> [[5]]
#> MappingResult for 'Megakaryocytes': best_bin=2 correlation=0.836872948898265 top_2_distance=0.056
#> Confident Result: FALSE (next max upper 0.90202722508178 )
#> with history for scores against 10 bins
#> Bootstrapped with 200 iterations
#>
#> [[6]]
#> MappingResult for 'GMP': best_bin=4 correlation=0.98187998124707 top_2_distance=0.0325
#> Confident Result: FALSE (next max upper 0.978960454097669 )
#> with history for scores against 10 bins
#> Bootstrapped with 200 iterations
To see how well we plotted the differences between these cell types, we can plot a heatmap of all these results. The neutrophils and Erythrocytes are confidently mapped, but the GMP population doesn’t map as well.
plot_mapping_result_heatmap(list(
multipotent_progenitors_result,
basophils_result,
neutrophils_result
), annotate_correlation = TRUE, text_background = TRUE)
To look into the Basophils mapping, we can plot the correlation over each bin. Because the the lower boundary of the best selection (bin 3) is higher than the higher bounds of the next-best bin (bin 4), BLASE makes a confident call.
plot_mapping_result_corr(basophils_result)
Now we can do some more detailed plotting about where the bins mapped onto the Single Cell data, and the cell type proportions we expect , but we need to add the pseudotime bins to the metadata of the SCE. We can see that the true proportions of cell types in the mapped bin do indeed map to the cell type we expect.
binned_sce <- assign_pseudotime_bins(
sce,
split_by = "pseudotime_range",
n_bins = 10,
pseudotime_slot = "pseudotime"
)
plot_mapping_result(
binned_sce,
multipotent_progenitors_result,
group_by_slot = "celltype"
)
To view the bin population chart in full, use:
plot_bin_population(binned_sce, 1, group_by_slot = "celltype")
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] reshape2_1.4.4 ami_0.2.1
#> [3] scran_1.37.0 slingshot_2.17.0
#> [5] TrajectoryUtils_1.17.0 princurve_2.1.6
#> [7] tradeSeq_1.23.0 ggVennDiagram_1.5.4
#> [9] limma_3.65.5 patchwork_1.3.2
#> [11] blase_0.99.3 BiocParallel_1.43.4
#> [13] scater_1.37.0 ggplot2_4.0.0
#> [15] scuttle_1.19.0 SingleCellExperiment_1.31.1
#> [17] SummarizedExperiment_1.39.2 Biobase_2.69.1
#> [19] GenomicRanges_1.61.5 Seqinfo_0.99.2
#> [21] IRanges_2.43.5 S4Vectors_0.47.4
#> [23] BiocGenerics_0.55.1 generics_0.1.4
#> [25] MatrixGenerics_1.21.0 matrixStats_1.5.0
#> [27] BiocStyle_2.37.1
#>
#> loaded via a namespace (and not attached):
#> [1] RcppAnnoy_0.0.22 splines_4.5.1 later_1.4.4
#> [4] tibble_3.3.0 polyclip_1.10-7 fastDummies_1.7.5
#> [7] lifecycle_1.0.4 edgeR_4.7.5 globals_0.18.0
#> [10] lattice_0.22-7 MASS_7.3-65 magrittr_2.0.4
#> [13] plotly_4.11.0 sass_0.4.10 rmarkdown_2.30
#> [16] jquerylib_0.1.4 yaml_2.3.10 metapod_1.17.0
#> [19] httpuv_1.6.16 Seurat_5.3.0 sctransform_0.4.2
#> [22] spam_2.11-1 sp_2.2-0 spatstat.sparse_3.1-0
#> [25] reticulate_1.43.0 cowplot_1.2.0 pbapply_1.7-4
#> [28] RColorBrewer_1.1-3 abind_1.4-8 Rtsne_0.17
#> [31] purrr_1.1.0 ggrepel_0.9.6 irlba_2.3.5.1
#> [34] listenv_0.9.1 spatstat.utils_3.2-0 goftest_1.2-3
#> [37] RSpectra_0.16-2 dqrng_0.4.1 spatstat.random_3.4-2
#> [40] fitdistrplus_1.2-4 parallelly_1.45.1 codetools_0.2-20
#> [43] DelayedArray_0.35.3 tidyselect_1.2.1 farver_2.1.2
#> [46] ScaledMatrix_1.17.0 viridis_0.6.5 spatstat.explore_3.5-3
#> [49] jsonlite_2.0.0 BiocNeighbors_2.3.1 progressr_0.16.0
#> [52] ggridges_0.5.7 survival_3.8-3 tools_4.5.1
#> [55] ica_1.0-3 Rcpp_1.1.0 glue_1.8.0
#> [58] gridExtra_2.3 SparseArray_1.9.1 xfun_0.53
#> [61] mgcv_1.9-3 dplyr_1.1.4 withr_3.0.2
#> [64] BiocManager_1.30.26 fastmap_1.2.0 bluster_1.19.0
#> [67] boot_1.3-32 digest_0.6.37 rsvd_1.0.5
#> [70] R6_2.6.1 mime_0.13 scattermore_1.2
#> [73] tensor_1.5.1 dichromat_2.0-0.1 spatstat.data_3.1-8
#> [76] tidyr_1.3.1 data.table_1.17.8 httr_1.4.7
#> [79] htmlwidgets_1.6.4 S4Arrays_1.9.1 uwot_0.2.3
#> [82] pkgconfig_2.0.3 gtable_0.3.6 lmtest_0.9-40
#> [85] S7_0.2.0 XVector_0.49.1 htmltools_0.5.8.1
#> [88] dotCall64_1.2 bookdown_0.45 SeuratObject_5.2.0
#> [91] scales_1.4.0 png_0.1-8 spatstat.univar_3.1-4
#> [94] knitr_1.50 nlme_3.1-168 cachem_1.1.0
#> [97] zoo_1.8-14 stringr_1.5.2 KernSmooth_2.23-26
#> [100] parallel_4.5.1 miniUI_0.1.2 vipor_0.4.7
#> [103] pillar_1.11.1 grid_4.5.1 vctrs_0.6.5
#> [106] RANN_2.6.2 promises_1.3.3 BiocSingular_1.25.0
#> [109] beachmat_2.25.5 xtable_1.8-4 cluster_2.1.8.1
#> [112] beeswarm_0.4.0 evaluate_1.0.5 tinytex_0.57
#> [115] magick_2.9.0 locfit_1.5-9.12 cli_3.6.5
#> [118] compiler_4.5.1 rlang_1.1.6 crayon_1.5.3
#> [121] future.apply_1.20.0 labeling_0.4.3 plyr_1.8.9
#> [124] ggbeeswarm_0.7.2 stringi_1.8.7 viridisLite_0.4.2
#> [127] deldir_2.0-4 lazyeval_0.2.2 spatstat.geom_3.6-0
#> [130] Matrix_1.7-4 RcppHNSW_0.6.0 sparseMatrixStats_1.21.0
#> [133] future_1.67.0 statmod_1.5.0 shiny_1.11.1
#> [136] ROCR_1.0-11 igraph_2.1.4 bslib_0.9.0