BLASE for excluding developmental genes from bulk RNA-seq
BLASE-for-excluding-developmental-genes-from-bulk-RNA-seq.Rmd
library(scater)
library(ggplot2)
library(BiocParallel)
library(blase)
library(limma)
library(ggVennDiagram)
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
}
bpparam <- MulticoreParam(N_CORES)This article will show how BLASE can be used for excluding the effect of developmental differences when analysing bulk RNA-seq data. We make use of scRNA-seq Dogga et al. 2024 and bulk RNA-seq data from Zhang et al. 2021. Code for generating the objects used here is available in the BLASE reproducibility repository.
Load Data
First we’ll load in the data we’re using, pre-prepared from the BLASE reproducibility repository.
zhang_bulk <- readRDS("Data/zhang_2021_heat_shock_bulk.rds")
dogga_sc <- readRDS("Data/MCA_PF_SCE.rds")We can examine the true lifecycle stages, and also the calculated pseudotime trajectory (Slingshot).
plotUMAP(dogga_sc, 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(dogga_sc, colour_by = "slingPseudotime_1")
## Prepare BLASE
Now we’ll prepare BLASE for use.
Create BLASE data object
First, we create the object, giving it the number of bins we want to use, and how to calculate them.
blase_data <- as.BlaseData(
dogga_sc,
pseudotime_slot = "slingPseudotime_1",
n_bins = 8,
split_by = "pseudotime_range"
)
# Add these bins to the sc metadata
dogga_sc <- blase::assign_pseudotime_bins(dogga_sc,
pseudotime_slot = "slingPseudotime_1",
n_bins = 8,
split_by = "pseudotime_range"
)Select Genes
Now we will select the genes we want to use, using BLASE’s gene peakedness metric.
gene_peakedness_info <- blase::calculate_gene_peakedness(
dogga_sc,
window_pct = 5,
knots = 18,
BPPARAM = bpparam
)
genes_to_use <- blase::gene_peakedness_spread_selection(
dogga_sc,
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
#> 72 PF3D7-1401100 9.811180 5.946666 1.576877 3.771165
#> 57 PF3D7-1401200 7.767184 4.928158 1.981323 2.487307
#> 59 PF3D7-1401300 8.039717 2.207782 1.285688 1.717199
#> 721 PF3D7-1401400 9.811180 942.073177 49.944485 18.862406
#> 15 PF3D7-1401600 2.043996 96.706165 8.207937 11.782031
#> 19 PF3D7-1402200 2.589061 3.516862 1.549402 2.269819
#> window_start window_end deviance_explained
#> 72 9.470514 10.151846 0.16478147
#> 57 7.426518 8.107850 0.19846899
#> 59 7.699051 8.380383 0.04479178
#> 721 9.470514 10.151846 0.27186813
#> 15 1.703330 2.384662 0.64855890
#> 19 2.248395 2.929727 0.17093181Here, we add them to the BLASE object for mapping with.
genes(blase_data) <- genes_to_useCalculate Mappings
Now we can perform the actual mapping step, and review the results.
mapping_results <- blase::map_all_best_bins(
blase_data,
zhang_bulk,
BPPARAM = bpparam
)
blase::plot_mapping_result_heatmap(mapping_results)
Calculate DE genes
We calculate DE genes using Limma (ref) as we have access to only normalised counts.
metadata <- data.frame(row.names = seq_len(12))
metadata$strain <- c(
rep("NF54", 4),
rep("PB4", 4),
rep("PB31", 4)
)
metadata$growth_conditions <- rep(c("Normal", "Normal", "HS", "HS"), 3)
metadata$group <- paste0(metadata$strain, "_", metadata$growth_conditions)
rownames(metadata) <- c(
"NF54_37_1",
"NF54_37_2",
"NF54_41_1",
"NF54_41_2",
"PB4_37_1",
"PB4_37_2",
"PB4_41_1",
"PB4_41_2",
"PB31_37_1",
"PB31_37_2",
"PB31_41_1",
"PB31_41_2"
)
design <- model.matrix(~ 0 + group, metadata)
colnames(design) <- gsub("group", "", colnames(design))
contr.matrix <- makeContrasts(
WT_normal_v_HS = NF54_Normal - NF54_HS,
PB31_normal_v_HS = PB31_Normal - PB31_HS,
normal_NF54_v_PB31 = NF54_Normal - PB31_Normal,
levels = colnames(design)
)Phenotype + Development (Bulk DE)
vfit <- lmFit(zhang_bulk, design)
vfit <- contrasts.fit(vfit, contrasts = contr.matrix)
efit <- eBayes(vfit)
summary(decideTests(efit))
#> WT_normal_v_HS PB31_normal_v_HS normal_NF54_v_PB31
#> Down 289 328 229
#> NotSig 2188 1884 2231
#> Up 80 345 97
PB31_normal_v_HS_BULK_DE <- topTable(efit, n = Inf, coef = 2)
PB31_normal_v_HS_BULK_DE <- PB31_normal_v_HS_BULK_DE[PB31_normal_v_HS_BULK_DE$adj.P.Val < 0.05 & PB31_normal_v_HS_BULK_DE$logFC > 0, ]
PB31_normal_v_HS_BULK_DE <- PB31_normal_v_HS_BULK_DE[order(-PB31_normal_v_HS_BULK_DE$logFC), ]Development (SC DE)
First, let’s pseudobulk based on the pseudotime bins
pseudobulked_dogga_sc <- data.frame(row.names = rownames(dogga_sc))
for (r in mapping_results) {
pseudobulked_dogga_sc[bulk_name(r)] <- rowSums(counts(dogga_sc[, dogga_sc$pseudotime_bin == best_bin(r)]))
}
pseudobulked_dogga_sc <- as.matrix(pseudobulked_dogga_sc)And now calculate DE genes for these
## Normalise, not needed for true bulks
par(mfrow = c(1, 2))
v <- voom(pseudobulked_dogga_sc, design, plot = FALSE)
vfit_sc <- lmFit(v, design)
vfit_sc <- contrasts.fit(vfit_sc, contrasts = contr.matrix)
efit_sc <- eBayes(vfit_sc)
summary(decideTests(efit_sc))
#> WT_normal_v_HS PB31_normal_v_HS normal_NF54_v_PB31
#> Down 0 0 0
#> NotSig 1768 1768 1768
#> Up 0 0 0Remove Developmental Genes
How many of these genes overlap? We can look at the intersection using a Venn Diagram:
ggVennDiagram(list(
Phenotype = rownames(PB31_normal_v_HS_BULK_DE),
Development = rownames(PB31_normal_v_HS_SC_DE)
)) + ggtitle("PB31 Normal vs Heat Shock")
Great, there are genes which we can correct. We can now remove these from the original DE genes list.
NB: This is a primitive approach, and there are likely many other and better ways to calculate which genes should be kept or removed.
PB31_normal_v_HS_corrected_DE <- rownames(PB31_normal_v_HS_BULK_DE[!(rownames(PB31_normal_v_HS_BULK_DE) %in% rownames(PB31_normal_v_HS_SC_DE)), ])
head(PB31_normal_v_HS_corrected_DE)
#> [1] "PF3D7-1003600" "PF3D7-1310700" "PF3D7-1035200" "PF3D7-0714000"
#> [5] "PF3D7-0725400" "PF3D7-1035500"Further downstream analysis can now be performed, as desired by the researcher, for example Gene Ontology term enrichment.
Session Info
sessionInfo()
#> R version 4.5.1 (2025-06-13)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.3 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
#>
#> Random number generation:
#> RNG: Mersenne-Twister
#> Normal: Inversion
#> Sample: Rounding
#>
#> locale:
#> [1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
#> [4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
#> [7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
#> [10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats4 stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] ggVennDiagram_1.5.4 limma_3.64.3
#> [3] blase_0.99.0 BiocParallel_1.42.1
#> [5] scater_1.36.0 ggplot2_3.5.2
#> [7] scuttle_1.18.0 SingleCellExperiment_1.30.1
#> [9] SummarizedExperiment_1.38.1 Biobase_2.68.0
#> [11] GenomicRanges_1.60.0 GenomeInfoDb_1.44.2
#> [13] IRanges_2.42.0 S4Vectors_0.46.0
#> [15] BiocGenerics_0.54.0 generics_0.1.4
#> [17] MatrixGenerics_1.20.0 matrixStats_1.5.0
#>
#> loaded via a namespace (and not attached):
#> [1] tidyselect_1.2.1 ami_0.2.1 viridisLite_0.4.2
#> [4] dplyr_1.1.4 vipor_0.4.7 farver_2.1.2
#> [7] viridis_0.6.5 fastmap_1.2.0 digest_0.6.37
#> [10] rsvd_1.0.5 lifecycle_1.0.4 statmod_1.5.0
#> [13] magrittr_2.0.3 compiler_4.5.1 rlang_1.1.6
#> [16] sass_0.4.10 tools_4.5.1 yaml_2.3.10
#> [19] knitr_1.50 labeling_0.4.3 S4Arrays_1.8.1
#> [22] htmlwidgets_1.6.4 DelayedArray_0.34.1 RColorBrewer_1.1-3
#> [25] abind_1.4-8 withr_3.0.2 desc_1.4.3
#> [28] grid_4.5.1 beachmat_2.24.0 scales_1.4.0
#> [31] cli_3.6.5 rmarkdown_2.29 crayon_1.5.3
#> [34] ragg_1.5.0 httr_1.4.7 ggbeeswarm_0.7.2
#> [37] cachem_1.1.0 splines_4.5.1 parallel_4.5.1
#> [40] XVector_0.48.0 vctrs_0.6.5 boot_1.3-31
#> [43] Matrix_1.7-3 jsonlite_2.0.0 BiocSingular_1.24.0
#> [46] patchwork_1.3.2 BiocNeighbors_2.2.0 ggrepel_0.9.6
#> [49] irlba_2.3.5.1 beeswarm_0.4.0 systemfonts_1.2.3
#> [52] jquerylib_0.1.4 glue_1.8.0 pkgdown_2.1.3
#> [55] codetools_0.2-20 cowplot_1.2.0 gtable_0.3.6
#> [58] UCSC.utils_1.4.0 ScaledMatrix_1.16.0 tibble_3.3.0
#> [61] pillar_1.11.0 htmltools_0.5.8.1 GenomeInfoDbData_1.2.14
#> [64] R6_2.6.1 sparseMatrixStats_1.20.0 textshaping_1.0.3
#> [67] evaluate_1.0.5 lattice_0.22-7 bslib_0.9.0
#> [70] Rcpp_1.1.0 gridExtra_2.3 SparseArray_1.8.1
#> [73] nlme_3.1-168 mgcv_1.9-3 xfun_0.53
#> [76] fs_1.6.6 pkgconfig_2.0.3