Introducing mRNA bias into simulations with scaling factors
Alexander Dietrich
11/29/2021
Source:vignettes/simulator_scaling_factors.Rmd
simulator_scaling_factors.Rmd
Using scaling factors
This package allows the user to introduce an mRNA bias into
pseudo-bulk RNA-seq samples. Different cell-types contain different
amounts of mRNA (Dendritic cells for examples contain much more than
Neutrophils); this bias can be added into simulations artificially in
different ways.
The scaling factors will always be applied on the single-cell dataset
first, altering its expression profiles accordingly, and then the
pseudo-bulk samples are generated by summing up the count data from the
sampled cells.
# Example data
counts <- Matrix::Matrix(matrix(stats::rpois(3e5, 5), ncol = 300), sparse = TRUE)
tpm <- Matrix::Matrix(matrix(stats::rpois(3e5, 5), ncol = 300), sparse = TRUE)
tpm <- Matrix::t(1e6 * Matrix::t(tpm) / Matrix::colSums(tpm))
colnames(counts) <- paste0("cell_", rep(1:300))
colnames(tpm) <- paste0("cell_", rep(1:300))
rownames(counts) <- paste0("gene_", rep(1:1000))
rownames(tpm) <- paste0("gene_", rep(1:1000))
annotation <- data.frame(
"ID" = paste0("cell_", rep(1:300)),
"cell_type" = c(rep("T cells CD4", 150), rep("T cells CD8", 150)),
"spikes" = stats::runif(300),
"add_1" = stats::runif(300),
"add_2" = stats::runif(300)
)
ds <- SimBu::dataset(
annotation = annotation,
count_matrix = counts,
name = "test_dataset"
)
#> Filtering genes...
#> Created dataset.
Pre-defined scaling factors
Some studies have proposed scaling factors for immune cells, such as EPIC (Racle et al. 2017) or quanTIseq (Finotello et al. 2019), deconvolution tools which correct for the mRNA bias internally using these values:
epic <- data.frame(
type = c(
"B cells",
"Macrophages",
"Monocytes",
"Neutrophils",
"NK cells",
"T cells",
"T cells CD4",
"T cells CD8",
"T helper cells",
"T regulatory cells",
"otherCells",
"default"
),
mRNA = c(
0.4016,
1.4196,
1.4196,
0.1300,
0.4396,
0.3952,
0.3952,
0.3952,
0.3952,
0.3952,
0.4000,
0.4000
)
)
epic
#> type mRNA
#> 1 B cells 0.4016
#> 2 Macrophages 1.4196
#> 3 Monocytes 1.4196
#> 4 Neutrophils 0.1300
#> 5 NK cells 0.4396
#> 6 T cells 0.3952
#> 7 T cells CD4 0.3952
#> 8 T cells CD8 0.3952
#> 9 T helper cells 0.3952
#> 10 T regulatory cells 0.3952
#> 11 otherCells 0.4000
#> 12 default 0.4000
quantiseq <- data.frame(
type = c(
"B cells",
"Macrophages",
"MacrophagesM2",
"Monocytes",
"Neutrophils",
"NK cells",
"T cells CD4",
"T cells CD8",
"T regulatory cells",
"Dendritic cells",
"T cells"
),
mRNA = c(
65.66148,
138.11520,
119.35447,
130.65455,
27.73634,
117.71584,
63.87200,
70.25659,
72.55110,
140.76091,
68.89323
)
)
quantiseq
#> type mRNA
#> 1 B cells 65.66148
#> 2 Macrophages 138.11520
#> 3 MacrophagesM2 119.35447
#> 4 Monocytes 130.65455
#> 5 Neutrophils 27.73634
#> 6 NK cells 117.71584
#> 7 T cells CD4 63.87200
#> 8 T cells CD8 70.25659
#> 9 T regulatory cells 72.55110
#> 10 Dendritic cells 140.76091
#> 11 T cells 68.89323
When you want to apply one of these scaling factors into your
simulation (therefore in-/decreasing the expression signals for the
cell-types), we can use the scaling_factor
parameter. Note,
that these pre-defined scaling factors only offer values for a certain
number of cell types, and your annotation in the provided dataset has to
match these names 1:1. All cell types from your dataset which are not
present in this scaling factor remain unscaled and a warning message
will appear.
sim_epic <- SimBu::simulate_bulk(
data = ds,
scenario = "random",
scaling_factor = "epic",
nsamples = 10,
ncells = 100,
BPPARAM = BiocParallel::MulticoreParam(workers = 4),
run_parallel = TRUE
)
#> Using parallel generation of simulations.
#> Using EPIC scaling factors.
#> Finished simulation.
We can also try out some custom scaling factors, for example
increasing the expression levels for a single cell-type (T cells CD8) by
10-fold compared to the rest. All cell-types which are not mentioned in
the named list given to custom_scaling_vector
will be
transformed with a scaling factor of 1, meaning nothing changes for
them.
sim_extreme_b <- SimBu::simulate_bulk(
data = ds,
scenario = "random",
scaling_factor = "custom",
custom_scaling_vector = c("T cells CD8" = 10),
nsamples = 10,
ncells = 100,
BPPARAM = BiocParallel::MulticoreParam(workers = 4),
run_parallel = TRUE
)
#> Using parallel generation of simulations.
#> Using custom scaling factors.
#> Warning in merge_scaling_factor(data = data, scaling_factor_values =
#> custom_scaling_vector, : For some cell type(s) in the dataset, no scaling
#> factor is available when using custom: T cells CD4. This cell type will not be
#> re-scaled.
#> Finished simulation.
Important: Watch out that the cell-type annotation names in your dataset are the same as in the scaling factor! Otherwise the scaling factor will not be applied or even worse, applied to a different cell-type.
Dataset specific scaling factors
You can also choose to calculate scaling factors, which are depending on your provided single-cell dataset. Compared to the previous section, this will give a unique value for each cell rather than a cell-type, making it possibly more sensitive.
Reads and genes
Two straight forward approaches would be the number of reads or number of expressed genes/features. As these values are easily obtainable from the provided count data, SimBu already calculates them during dataset generation.
sim_reads <- SimBu::simulate_bulk(
data = ds,
scenario = "random",
scaling_factor = "read_number", # use number of reads as scaling factor
nsamples = 10,
ncells = 100,
BPPARAM = BiocParallel::MulticoreParam(workers = 4),
run_parallel = TRUE
)
#> Using parallel generation of simulations.
#> Finished simulation.
sim_genes <- SimBu::simulate_bulk(
data = ds,
scenario = "random",
scaling_factor = "expressed_genes", # use number of expressed genes column as scaling factor
nsamples = 10,
ncells = 100,
BPPARAM = BiocParallel::MulticoreParam(workers = 4),
run_parallel = TRUE
)
#> Using parallel generation of simulations.
#> Finished simulation.
These options would also allow you to use other numerical
measurements you have for the single cells as scaling factors, such as
weight or size for example. Lets pretend, add_1
and
add_2
are such measurements. With the
additional_cols
parameter, they can be added to the SimBu
dataset and we can use them as scaling factor as well:
utils::head(annotation)
#> ID cell_type spikes add_1 add_2
#> 1 cell_1 T cells CD4 0.79817472 0.5437374 0.3959318
#> 2 cell_2 T cells CD4 0.58401273 0.2811142 0.7472875
#> 3 cell_3 T cells CD4 0.08633433 0.1951225 0.1903732
#> 4 cell_4 T cells CD4 0.94282994 0.7435997 0.3436296
#> 5 cell_5 T cells CD4 0.60404503 0.5043812 0.3328758
#> 6 cell_6 T cells CD4 0.02199540 0.5871896 0.7028115
ds <- SimBu::dataset(
annotation = annotation,
count_matrix = counts,
name = "test_dataset",
additional_cols = c("add_1", "add_2")
) # add columns to dataset
#> Filtering genes...
#> Created dataset.
sim_genes <- SimBu::simulate_bulk(
data = ds,
scenario = "random",
scaling_factor = "add_1", # use add_1 column as scaling factor
nsamples = 10,
ncells = 100,
BPPARAM = BiocParallel::MulticoreParam(workers = 4),
run_parallel = TRUE
)
#> Using parallel generation of simulations.
#> Scaling by custom column in annotation table; if no scaling is wished instead, use 'NONE'.
#> Finished simulation.
Spike-ins
One other numerical measurement can be spike-ins. Usually the number
of reads mapped to spike-in molecules per cell is given in the cell
annotations. If this is the case, they can be stored in the dataset
annotation using the spike_in_col
parameter, where you
indicate the name of the column from the annotation
dataframe in which the spike-in information is stored. To calculate a
scaling factor from this, the number of reads are also necessary, so we
will add this information as well (as above using the
read_number_col
parameter).
The scaling factor with spike-ins is calculated as the “% of reads
NOT mapped to spike-in reads”, or:
(n_reads - n_spike_in)/n_reads
for each cell. We apply it
like this:
ds <- SimBu::dataset(
annotation = annotation,
count_matrix = counts,
name = "test_dataset",
spike_in_col = "spikes"
) # give the name in the annotation file, that contains spike-in information
sim_spike <- SimBu::simulate_bulk(
data = ds,
scenario = "random",
scaling_factor = "spike_in", # use spike-in scaling factor
nsamples = 10,
ncells = 100,
BPPARAM = BiocParallel::MulticoreParam(workers = 4),
run_parallel = TRUE
)
Census - estimate mRNA counts per cell
Census
is an approach which tries to convert TPM counts into relative
transcript counts. This basically means, you get the mRNA counts per
cell, which can differ between cell-types.
(Qiu et al.
2017) state in their paper, that it should only be applied to
TPM/FPKM normalized data, but I tried it out with raw expression counts
as well, which worked as well.
Census calculates a vector with a scaling value for each cell in a
sample. You can switch this feature on, by setting the
scaling_factor
parameter to census
.
ds <- SimBu::dataset(
annotation = annotation,
count_matrix = counts,
tpm_matrix = tpm,
name = "test_dataset"
)
#> Filtering genes...
#> Created dataset.
sim_census <- SimBu::simulate_bulk(
data = ds,
scenario = "random",
scaling_factor = "census", # use census scaling factor
nsamples = 10,
ncells = 100,
BPPARAM = BiocParallel::MulticoreParam(workers = 4),
run_parallel = TRUE
)
#> Using parallel generation of simulations.
#> Finished simulation.
In our analysis we found, that Census is basically a complicated way of estimating the number of expressed genes per cell. It will remain to the user to decide if he/she wants to use census or simply the number of expressed genes (as shown above) as scaling factor.
utils::sessionInfo()
#> R version 4.3.3 (2024-02-29)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 22.04.4 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.20.so; LAPACK version 3.10.0
#>
#> 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] stats graphics grDevices utils datasets methods base
#>
#> loaded via a namespace (and not attached):
#> [1] sass_0.4.9 SparseArray_1.2.4
#> [3] bitops_1.0-7 lattice_0.22-5
#> [5] digest_0.6.35 magrittr_2.0.3
#> [7] evaluate_0.23 sparseMatrixStats_1.14.0
#> [9] grid_4.3.3 fastmap_1.1.1
#> [11] jsonlite_1.8.8 Matrix_1.6-5
#> [13] GenomeInfoDb_1.38.8 proxyC_0.3.4
#> [15] SimBu_1.5.4 purrr_1.0.2
#> [17] codetools_0.2-19 textshaping_0.3.7
#> [19] jquerylib_0.1.4 abind_1.4-5
#> [21] cli_3.6.2 rlang_1.1.3
#> [23] crayon_1.5.2 XVector_0.42.0
#> [25] Biobase_2.62.0 cachem_1.0.8
#> [27] DelayedArray_0.28.0 yaml_2.3.8
#> [29] S4Arrays_1.2.1 parallel_4.3.3
#> [31] tools_4.3.3 BiocParallel_1.36.0
#> [33] memoise_2.0.1 GenomeInfoDbData_1.2.11
#> [35] SummarizedExperiment_1.32.0 BiocGenerics_0.48.1
#> [37] vctrs_0.6.5 R6_2.5.1
#> [39] matrixStats_1.2.0 stats4_4.3.3
#> [41] lifecycle_1.0.4 zlibbioc_1.48.2
#> [43] S4Vectors_0.40.2 fs_1.6.3
#> [45] IRanges_2.36.0 ragg_1.3.0
#> [47] desc_1.4.3 pkgdown_2.0.7
#> [49] RcppParallel_5.1.7 bslib_0.6.2
#> [51] Rcpp_1.0.12 systemfonts_1.0.6
#> [53] xfun_0.43 GenomicRanges_1.54.1
#> [55] MatrixGenerics_1.14.0 knitr_1.45
#> [57] htmltools_0.5.8 rmarkdown_2.26
#> [59] compiler_4.3.3 RCurl_1.98-1.14