Contents

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
}

0.1 Introduction

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.

0.1.1 Principle

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?

Plot showing gene expression varying over time for three genes, with lines marking different areas which can be uniquely identified by the rank-order of the expression of the genes.

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:

  • A Single Cell Experiment with pseudotime.
  • A list of genes that explain the trajectory well, in this case generated by tradeSeq’s fitGAM and associationTest

0.2 Installing BLASE

To install BLASE from Bioconductor:

if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager")

BiocManager::install("blase")

0.3 Setting up the Single Cell Experiment

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")

Plots showing a UMAP of the single-cell data, coloured by celltype and pseudotime. The UMAP is shaped like a horizontal line, showing the differentiation of a set of immune cell types.

0.4 Finding the most descriptive genes with tradeSeq

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.

0.5 Parameter Tuning for BLASE

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)

Figures showing the minimum and mean convexity, as well as % confident mappings different numbers of genes and bins calculated.

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.

0.5.1 Inspect Bin Choice

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)

These plots show the correlations of each bin to itself and all other bins. This shows how distinct each bin is compared to others, based on the number of bins and genes selected.

#> [1]  0.00140  0.04096 20.00000

0.5.2 Inspect Genes Choice

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%” />

0.6 Mapping Bulk Samples to SC with BLASE

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

0.6.1 Plotting Heatmap of Mappings

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)

A heatmap showing the correlations BLASE found for a subset of cell types in this data. Confident mappings are indicated with asterisks over the mapped tile. This information is also available as text from the results of the previous code block.

0.6.2 Plotting Detailed Correlation Maps

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)

A plot of the correlation of each bin with the basophils pseudobulk sample, also showing upper and lower bounds as calculated by BLASE. The best mapping is for bin 3, and bin 4 is the next most similar.

0.6.3 Plotting Summary Plots of Mappings

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"
)

An overview of mapping for Multipotent progenitors, including UMAPs of all cells and just those in the mapped pseudotime bin coloured by pseudotime bin, and cell-type. Also includes a line plot of correlation of the mapping over every bin, and a bar plot of the different cell-type populations in that bin, showing that the bulk was correctly mapped to a bin consisting mostly of Multipotent progenitor cells.

To view the bin population chart in full, use:

plot_bin_population(binned_sce, 1, group_by_slot = "celltype")

A full size bar plot of the cell-type populations of bin 1.

0.7 Session Info

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