Pseudobulk cell size rescaling example (original) (raw)
This vignette shows an example pseudobulk experiment testing cell size scale factors using a small example dataset of single-nucleus RNA-seq data (snRNA-seq) from human cortex (Darmanis et al. (2015)). Predictions are made using lute, and results plots are generated using ggplot2.
Experiment setup
Load example data
In this example, we source a real snRNA-seq dataset of human brain, including cortex and hippocampus published in darmanis_survey_2015. The full data along with other real single-cell RNA-seq datasets may be accessed from the scRNAseqpackage.
Load a stored subset of the example dataset with the following.
path <- system.file("extdata", "scRNAseq/darmanis_example.rda", package="lute")
data <- get(load(path))The loaded dataset is of type SingleCellExperiment, which is handled by thelute() function (see ?lute for details). Before calling the framework function, it needs to be processed to (1) define cell types and samples of interest (2) subset on cell type markers, and (3) define pseudobulks for each available sample.
For this experiment, we will consider two principal cell types for brain, neuron and glial cells (a.k.a. “K2”).
Define cell types of interest
First, identify nuclei labeled from only these types and remove the rest. Then define a new label "k2" using the valid remaining nuclei.
sampleIdVariable <- "experiment_sample_name"
oldTypes <- "cell.type"; newTypes <- "k2"
# remove non-k2 types
filterK2 <- data[[oldTypes]] %in%
c("neurons", "oligodendrocytes", "astrocytes", "OPC", "microglia")
data <- data[,filterK2]
# define new k2 variable
data[[newTypes]] <- ifelse(data[[oldTypes]]=="neurons", "neuron", "glial")
data[[newTypes]] <- factor(data[[newTypes]])Filter samples
Next, define the samples of interest for the experiment. We will select samples having at least 20 nuclei.
minNuclei <- 20
nucleiPerSample <- table(data[[sampleIdVariable]])
sampleIdVector <- unique(data[[sampleIdVariable]])
sampleIdVector <- sampleIdVector[nucleiPerSample >= minNuclei]
sampleIdVector # view## [1] "AB_S8" "AB_S11" "AB_S3" "AB_S4" "AB_S5" "AB_S7"Next, save samples having non-zero amounts of neuron and glial cells.
sampleIdVector <- unlist(lapply(sampleIdVector, function(sampleId){
numTypes <- length(
unique(
data[,data[[sampleIdVariable]]==sampleId][[newTypes]]))
if(numTypes==2){sampleId}
}))
sampleIdVector## [1] "AB_S11" "AB_S4" "AB_S5" "AB_S7"View the summaries by sample id, then save these as the true cell type proportions. These will be used later to assess the predictions.
proportionsList <- lapply(sampleIdVector, function(sampleId){
prop.table(table(data[,data$experiment_sample_name==sampleId]$k2))
})
dfProportions <- do.call(rbind, proportionsList)
rownames(dfProportions) <- sampleIdVector
colnames(dfProportions) <- paste0(colnames(dfProportions), ".true")
dfProportions <- as.data.frame(dfProportions)
knitr::kable(dfProportions) # viewMake pseudobulks with cell size scale factors
Define the cell size scale factors and use these to make the pseudobulks. For demonstration we set these to have large difference (i.e. neuron/glial > 3). While we set these manually, the cell scale factors could also be defined from library sizes or by referencing the cellScaleFactors package (link).
cellScalesVector <- c("glial" = 3, "neuron" = 10)Next make the pseudobulk datasets.
assayName <- "counts"
pseudobulkList <- lapply(sampleIdVector, function(sampleId){
dataIteration <- data[,data[[sampleIdVariable]]==sampleId]
ypb_from_sce(
singleCellExperiment = dataIteration, assayName = assayName,
cellTypeVariable = newTypes, cellScaleFactors = cellScalesVector)
})
dfPseudobulk <- do.call(cbind, pseudobulkList)
dfPseudobulk <- as.data.frame(dfPseudobulk)
colnames(dfPseudobulk) <- sampleIdVector
knitr::kable(head(dfPseudobulk))Get predictions
Predictions with scaling
Predict the neuron proportions using non-negative least squares (NNLS), the default deconvolution algorithm used by lute(). First, get the scaled proportions by setting the argument cellScaleFactors = cellScalesVector.
scaledResult <- lute(
singleCellExperiment = data,
bulkExpression = as.matrix(dfPseudobulk),
cellScaleFactors = cellScalesVector,
typemarkerAlgorithm = NULL,
cellTypeVariable = newTypes,
assayName = assayName)## Parsing deconvolution arguments...## Using NNLS...proportions.scaled <- scaledResult[[1]]@predictionsTable
knitr::kable(proportions.scaled) # viewPredictions without scaling
Next, get the unscaled result without setting s.
unscaledResult <- lute(
singleCellExperiment = data,
bulkExpression = as.matrix(dfPseudobulk),
typemarkerAlgorithm = NULL,
cellTypeVariable = newTypes,
assayName = assayName)## Parsing deconvolution arguments...## Using NNLS...proportionsUnscaled <- unscaledResult[[1]]@predictionsTable
knitr::kable(proportionsUnscaled) # viewNote proportions didn’t change for samples which had all glial or all neuron (AB_S8 and AB_S3).
Plot differences
Get the plot data tables
We will show the outcome of performing the cell scale factor adjustments using scatterplots and boxplots. Begin by appending the neuron proportion predictions from scaling treatments (scaled and unscaled) to the true proportions tabledfProportions.
dfProportions$neuron.unscaled <- proportionsUnscaled$neuron
dfProportions$neuron.scaled <- proportions.scaled$neuron
knitr::kable(dfProportions) # viewCalculate bias as the difference between true and predicted neuron proportions. Then calculate the error as the absolute of the bias thus defined.
# get bias
dfProportions$bias.neuron.unscaled <-
dfProportions$neuron.true-dfProportions$neuron.unscaled
dfProportions$bias.neuron.scaled <-
dfProportions$neuron.true-dfProportions$neuron.scaled
# get error
dfProportions$error.neuron.unscaled <-
abs(dfProportions$bias.neuron.unscaled)
dfProportions$error.neuron.scaled <-
abs(dfProportions$bias.neuron.scaled)Make the tall version of dfProportions in order to generate a plot with facets on the scale treatment (either “scaled” or “unscaled”).
dfPlotTall <- rbind(
data.frame(true = dfProportions$neuron.true,
predicted = dfProportions$neuron.scaled,
error = dfProportions$error.neuron.scaled,
sampleId = rownames(dfProportions),
type = rep("scaled", nrow(dfProportions))),
data.frame(true = dfProportions$neuron.true,
predicted = dfProportions$neuron.unscaled,
error = dfProportions$error.neuron.unscaled,
sampleId = rownames(dfProportions),
type = rep("unscaled", nrow(dfProportions)))
)
dfPlotTall <- as.data.frame(dfPlotTall)Make scatterplots of true versus predicted neuron proportions
Show sample results scatterplots of true (x-axis) by predicted (y-axis) neuron proportions. Also include a reference line (slope = 1, yintercept = 0) showing where agreement is absolute between proportions.
Also shows RMSE in plot titles.
dfPlotTallNew <- dfPlotTall
rmseScaled <-
rmse(
dfPlotTallNew[dfPlotTallNew$type=="scaled",]$true,
dfPlotTall[dfPlotTall$type=="scaled",]$predicted, "mean")
rmseUnscaled <-
rmse(
dfPlotTallNew[dfPlotTallNew$type=="unscaled",]$true,
dfPlotTallNew[dfPlotTallNew$type=="unscaled",]$predicted, "mean")
dfPlotTallNew$type <-
ifelse(grepl("un.*", dfPlotTallNew$type),
paste0(dfPlotTallNew$type,
" (RMSE = ", round(rmseScaled, 3), ")"),
paste0(dfPlotTallNew$type,
" (RMSE = ", round(rmseUnscaled, 3), ")"))
textSize <- 15
ggplot(dfPlotTallNew, aes(x = true, y = predicted)) +
geom_point(size = 4, alpha = 0.5) + geom_abline(slope = 1, intercept = 0) +
xlim(0, 1) + ylim(0, 1) + facet_wrap(~type) + theme_bw() +
xlab("True") + ylab("Predicted") +
theme(text = element_text(size = textSize),
axis.text.x = element_text(angle = 45, hjust = 1))Make boxplots with jittered points for neuron errors
Show jitters and boxplots by sample, depicting the neuron error (y-axis) by scale treatment (x-axis). The sample IDs are depicted by the point colors.
ggplot(dfPlotTall, aes(x = type, y = error, color = sampleId)) +
geom_jitter(alpha = 0.5, size = 4) + theme_bw() +
geom_boxplot(alpha = 0, color = "cyan") +
theme(text = element_text(size = textSize),
axis.text.x = element_text(angle = 45, hjust = 1)) +
xlab("Type") + ylab("Error (Neuron)")This process could be readily repeated for the remaining cell types, or just glial cells in this case.
Conclusions
This vignette showed how to conduct a basic pseudobulk experiment using cell size scale factors and an example snRNAseq dataset from human brainDarmanis et al. (2015). Some key details include sourcing and snRNA-seq data, defining a new cell type variable, setting the scale factors, making predictions, and performing comparative analyses of the prediction results. Further details about the importance of cell size scale factors are discussed inMaden et al. (2023), and examples of their utilizations may be found inMonaco et al. (2019), Racle and Gfeller (2020), and Sosina et al. (2021).
Session info
sessionInfo()## R version 4.5.0 RC (2025-04-04 r88126)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.2 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.21-bioc/R/lib/libRblas.so
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0
##
## 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] ggplot2_3.5.2 lute_1.4.0
## [3] SingleCellExperiment_1.30.0 SummarizedExperiment_1.38.0
## [5] Biobase_2.68.0 GenomicRanges_1.60.0
## [7] GenomeInfoDb_1.44.0 IRanges_2.42.0
## [9] S4Vectors_0.46.0 BiocGenerics_0.54.0
## [11] generics_0.1.3 MatrixGenerics_1.20.0
## [13] matrixStats_1.5.0 BiocStyle_2.36.0
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.2.1 dplyr_1.1.4 farver_2.1.2
## [4] fastmap_1.2.0 bluster_1.18.0 digest_0.6.37
## [7] rsvd_1.0.5 lifecycle_1.0.4 cluster_2.1.8.1
## [10] statmod_1.5.0 magrittr_2.0.3 compiler_4.5.0
## [13] rlang_1.1.6 sass_0.4.10 tools_4.5.0
## [16] igraph_2.1.4 yaml_2.3.10 knitr_1.50
## [19] S4Arrays_1.8.0 labeling_0.4.3 dqrng_0.4.1
## [22] DelayedArray_0.34.0 abind_1.4-8 BiocParallel_1.42.0
## [25] withr_3.0.2 grid_4.5.0 beachmat_2.24.0
## [28] colorspace_2.1-1 edgeR_4.6.0 scales_1.3.0
## [31] tinytex_0.57 cli_3.6.4 rmarkdown_2.29
## [34] crayon_1.5.3 metapod_1.16.0 httr_1.4.7
## [37] scuttle_1.18.0 cachem_1.1.0 nnls_1.6
## [40] parallel_4.5.0 BiocManager_1.30.25 XVector_0.48.0
## [43] vctrs_0.6.5 Matrix_1.7-3 jsonlite_2.0.0
## [46] bookdown_0.43 BiocSingular_1.24.0 BiocNeighbors_2.2.0
## [49] irlba_2.3.5.1 magick_2.8.6 locfit_1.5-9.12
## [52] limma_3.64.0 jquerylib_0.1.4 glue_1.8.0
## [55] codetools_0.2-20 gtable_0.3.6 UCSC.utils_1.4.0
## [58] ScaledMatrix_1.16.0 munsell_0.5.1 tibble_3.2.1
## [61] pillar_1.10.2 htmltools_0.5.8.1 GenomeInfoDbData_1.2.14
## [64] R6_2.6.1 evaluate_1.0.3 lattice_0.22-7
## [67] scran_1.36.0 bslib_0.9.0 Rcpp_1.0.14
## [70] SparseArray_1.8.0 xfun_0.52 pkgconfig_2.0.3Works cited
Darmanis, Spyros, Steven A. Sloan, Ye Zhang, Martin Enge, Christine Caneda, Lawrence M. Shuer, Melanie G. Hayden Gephart, Ben A. Barres, and Stephen R. Quake. 2015. “A Survey of Human Brain Transcriptome Diversity at the Single Cell Level.” Proceedings of the National Academy of Sciences of the United States of America 112 (23): 7285–90. https://doi.org/10.1073/pnas.1507125112.
Maden, Sean K., Sang Ho Kwon, Louise A. Huuki-Myers, Leonardo Collado-Torres, Stephanie C. Hicks, and Kristen R. Maynard. 2023. “Challenges and Opportunities to Computationally Deconvolve Heterogeneous Tissue with Varying Cell Sizes Using Single Cell RNA-Sequencing Datasets.” arXiv. https://doi.org/10.48550/arXiv.2305.06501.
Monaco, Gianni, Bernett Lee, Weili Xu, Seri Mustafah, You Yi Hwang, Christophe Carré, Nicolas Burdin, et al. 2019. “RNA-Seq Signatures Normalized by mRNA Abundance Allow Absolute Deconvolution of Human Immune Cell Types.” Cell Reports 26 (6): 1627–1640.e7. https://doi.org/10.1016/j.celrep.2019.01.041.
Racle, Julien, and David Gfeller. 2020. “EPIC: A Tool to Estimate the Proportions of Different Cell Types from Bulk Gene Expression Data.” Edited by Sebastian Boegel, Methods in Molecular Biology,, 233–48. https://doi.org/10.1007/978-1-0716-0327-7_17.
Sosina, Olukayode A., Matthew N. Tran, Kristen R. Maynard, Ran Tao, Margaret A. Taub, Keri Martinowich, Stephen A. Semick, et al. 2021. “Strategies for Cellular Deconvolution in Human Brain RNA Sequencing Data,” no. 10:750 (August). https://doi.org/10.12688/f1000research.50858.1.