Overview of the sccomp package
Stefano Mangiola
2024-10-15
Source:vignettes/introduction.Rmd
introduction.RmdAbstract
Sccomp is a generalised method for differential composition and variability analyses. It tests differences in cell type proportion, for count data. It is able of flexible continuous/discrete models, removal of unwanted variation, removal of outliers, multilevel random/effect modelling.
sccomp
tests differences in cell type proportions from single-cell data. It is
robust against outliers, it models continuous and discrete factors, and
capable of random-effect/intercept modelling.
Please cite PNAS - sccomp: Robust differential composition and variability analysis for single-cell data
Characteristics
- Complex linear models with continuous and categorical covariates
- Multilevel modelling, with population fixed and random effects/intercept
- Modelling data from counts
- Testing differences in cell-type proportionality
- Testing differences in cell-type specific variability
- Cell-type information share for variability adaptive shrinkage
- Testing differential variability
- Probabilistic outlier identification
- Cross-dataset learning (hyperpriors).
Installation
sccomp is based on cmdstanr which provides
the latest version of cmdstan the Bayesian modelling tool.
cmdstanr is not on CRAN, so we need to have 3 simple step
process (that will be prompted to the user is forgot).
- R installation of
sccomp - R installation of
cmdstanr -
cmdstanrcall tocmdstaninstallation
Bioconductor
if (!requireNamespace("BiocManager")) install.packages("BiocManager")
# Step 1
BiocManager::install("sccomp")
# Step 2
install.packages("cmdstanr", repos = c("https://stan-dev.r-universe.dev/", getOption("repos")))
# Step 3
cmdstanr::check_cmdstan_toolchain(fix = TRUE) # Just checking system setting
cmdstanr::install_cmdstan()Github
# Step 1
devtools::install_github("MangiolaLaboratory/sccomp")
# Step 2
install.packages("cmdstanr", repos = c("https://stan-dev.r-universe.dev/", getOption("repos")))
# Step 3
cmdstanr::check_cmdstan_toolchain(fix = TRUE) # Just checking system setting
cmdstanr::install_cmdstan()| Function | Description |
|---|---|
sccomp_estimate |
Fit the model onto the data, and estimate the coefficients |
sccomp_remove_outliers |
Identify outliers probabilistically based on the model fit, and exclude them from the estimation |
sccomp_test |
Calculate the probability that the coefficients are outside the H0 interval (i.e. test_composition_above_logit_fold_change) |
sccomp_replicate |
Simulate data from the model, or part of the model |
sccomp_predict |
Predicts proportions, based on the model, or part of the model |
sccomp_remove_unwanted_variation |
Removes the variability for unwanted factors |
plot |
Plots summary plots to asses significance |
Analysis
library(dplyr)
library(sccomp)
library(ggplot2)
library(forcats)
library(tidyr)
data("seurat_obj")
data("sce_obj")
data("counts_obj")sccomp can model changes in composition and variability.
By default, the formula for variability is either ~1, which
assumes that the cell-group variability is independent of any covariate
or ~ factor_of_interest, which assumes that the model is
dependent on the factor of interest only. The variability model must be
a subset of the model for composition.
Binary factor
Of the output table, the estimate columns start with the prefix
c_ indicate composition, or with
v_ indicate variability (when
formula_variability is set).
From Seurat, SingleCellExperiment, metadata objects
sccomp_result =
sce_obj |>
sccomp_estimate(
formula_composition = ~ type,
.sample = sample,
.cell_group = cell_group,
cores = 1
) |>
sccomp_remove_outliers(cores = 1) |> # Optional
sccomp_test()From counts
sccomp_result =
counts_obj |>
sccomp_estimate(
formula_composition = ~ type,
.sample = sample,
.cell_group = cell_group,
.count = count,
cores = 1, verbose = FALSE
) |>
sccomp_remove_outliers(cores = 1, verbose = FALSE) |> # Optional
sccomp_test()Here you see the results of the fit, the effects of the factor on composition and variability. You also can see the uncertainty around those effects.
The output is a tibble containing the Following columns
-
cell_group- The cell groups being tested. -
parameter- The parameter being estimated from the design matrix described by the inputformula_compositionandformula_variability. -
factor- The covariate factor in the formula, if applicable (e.g., not present for Intercept or contrasts). -
c_lower- Lower (2.5%) quantile of the posterior distribution for a composition (c) parameter. -
c_effect- Mean of the posterior distribution for a composition (c) parameter. -
c_upper- Upper (97.5%) quantile of the posterior distribution for a composition (c) parameter. -
c_pH0- Probability of the null hypothesis (no difference) for a composition (c). This is not a p-value. -
c_FDR- False-discovery rate of the null hypothesis for a composition (c). -
v_lower- Lower (2.5%) quantile of the posterior distribution for a variability (v) parameter. -
v_effect- Mean of the posterior distribution for a variability (v) parameter. -
v_upper- Upper (97.5%) quantile of the posterior distribution for a variability (v) parameter. -
v_pH0- Probability of the null hypothesis for a variability (v). -
v_FDR- False-discovery rate of the null hypothesis for a variability (v). -
count_data- Nested input count data.
sccomp_resultSummary plots
A plot of group proportions, faceted by groups. The blue boxplots represent the posterior predictive check. If the model is descriptively adequate for the data, the blue boxplots should roughly overlay the black boxplots, which represent the observed data. The outliers are coloured in red. A boxplot will be returned for every (discrete) covariate present in formula_composition. The colour coding represents the significant associations for composition and/or variability.
sccomp_result |>
sccomp_boxplot(factor = "type")A plot of estimates of differential composition (c_) on the x-axis and differential variability (v_) on the y-axis. The error bars represent 95% credible intervals. The dashed lines represent the minimal effect that the hypothesis test is based on. An effect is labelled as significant if it exceeds the minimal effect according to the 95% credible interval. Facets represent the covariates in the model.
sccomp_result |>
plot_1D_intervals()We can plot the relationship between abundance and variability. As we can see below, they are positively correlated. sccomp models this relationship to obtain a shrinkage effect on the estimates of both the abundance and the variability. This shrinkage is adaptive as it is modelled jointly, thanks to Bayesian inference.
sccomp_result |>
plot_2D_intervals()You can produce the series of plots calling the plot
method.
sccomp_result |> plot() Model proportions directly (e.g. from deconvolution)
Note: If counts are available, we strongly discourage the use of proportions, as an important source of uncertainty (i.e., for rare groups/cell types) is not modeled.
The use of proportions is better suited for modelling deconvolution results (e.g., of bulk RNA data), in which case counts are not available.
Proportions should be greater than 0. Assuming that zeros derive from a precision threshold (e.g., deconvolution), zeros are converted to the smallest non-zero value.
Continuous factor
sccomp is able to fit erbitrary complex models. In this
example we have a continuous and binary covariate.
res =
seurat_obj |>
sccomp_estimate(
formula_composition = ~ type + continuous_covariate,
.sample = sample, .cell_group = cell_group,
cores = 1, verbose=FALSE
)
resRandom Effect Modeling
sccomp supports multilevel modeling by allowing the
inclusion of random effects in the compositional and variability
formulas. This is particularly useful when your data has hierarchical or
grouped structures, such as measurements nested within subjects,
batches, or experimental units. By incorporating random effects, sccomp
can account for variability at different levels of your data, improving
model fit and inference accuracy.
Random Intercept Model
In this example, we demonstrate how to fit a random intercept model using sccomp. We’ll model the cell-type proportions with both fixed effects (e.g., treatment) and random effects (e.g., subject-specific variability).
Here is the input data
seurat_obj[[]] |> as_tibble()## Loading required package: SeuratObject## Loading required package: sp## 'SeuratObject' was built under R 4.4.0 but the current version is
## 4.4.1; it is recomended that you reinstall 'SeuratObject' as the ABI
## for R may have changed##
## Attaching package: 'SeuratObject'## The following objects are masked from 'package:base':
##
## intersect, t## # A tibble: 106,297 × 9
## cell_group nCount_RNA nFeature_RNA group__ group__wrong sample type group2__
## <chr> <dbl> <int> <chr> <chr> <chr> <chr> <chr>
## 1 CD4 naive 0 0 GROUP1 1 SI-GA… canc… GROUP21
## 2 Mono clas… 0 0 GROUP1 1 SI-GA… canc… GROUP21
## 3 CD4 cm S1… 0 0 GROUP1 1 SI-GA… canc… GROUP21
## 4 B immature 0 0 GROUP1 1 SI-GA… canc… GROUP21
## 5 CD8 naive 0 0 GROUP1 1 SI-GA… canc… GROUP21
## 6 CD4 naive 0 0 GROUP1 1 SI-GA… canc… GROUP21
## 7 Mono clas… 0 0 GROUP1 1 SI-GA… canc… GROUP21
## 8 CD4 cm S1… 0 0 GROUP1 1 SI-GA… canc… GROUP21
## 9 CD4 cm hi… 0 0 GROUP1 1 SI-GA… canc… GROUP21
## 10 B immature 0 0 GROUP1 1 SI-GA… canc… GROUP21
## # ℹ 106,287 more rows
## # ℹ 1 more variable: continuous_covariate <dbl>
res =
seurat_obj |>
sccomp_estimate(
formula_composition = ~ type + (1 | group__),
.sample = sample,
.cell_group = cell_group,
bimodal_mean_variability_association = TRUE,
cores = 1, verbose = FALSE
)
resRandom Effect Model (random slopes)
sccomp can model random slopes. We providean example
below.
res =
seurat_obj |>
sccomp_estimate(
formula_composition = ~ type + (type | group__),
.sample = sample,
.cell_group = cell_group,
bimodal_mean_variability_association = TRUE,
cores = 1, verbose = FALSE
)
resNested Random Effects
If you have a more complex hierarchy, such as measurements nested
within subjects and subjects nested within batches, you can include
multiple grouping variables. Here group2__ is nested within
group__.
res =
seurat_obj |>
sccomp_estimate(
formula_composition = ~ type + (type | group__) + (1 | group2__),
.sample = sample,
.cell_group = cell_group,
bimodal_mean_variability_association = TRUE,
cores = 1, verbose = FALSE
)
resAn aid to result interpretation and communication
The estimated effects are expressed in the unconstrained space of the parameters, similar to differential expression analysis that expresses changes in terms of log fold change. However, for differences in proportion, logit fold change must be used, which is harder to interpret and understand.
Therefore, we provide a more intuitive proportional fold change that can be more easily understood. However, these cannot be used to infer significance (use sccomp_test() instead), and a lot of care must be taken given the nonlinearity of these measures (a 1-fold increase from 0.0001 to 0.0002 carries a different weight than a 1-fold increase from 0.4 to 0.8).
From your estimates, you can specify which effects you are interested in (this can be a subset of the full model if you wish to exclude unwanted effects), and the two points you would like to compare.
In the case of a categorical variable, the starting and ending points are categories.
sccomp_result |>
sccomp_proportional_fold_change(
formula_composition = ~ type,
from = "healthy",
to = "cancer"
) |>
select(cell_group, statement)Contrasts
seurat_obj |>
sccomp_estimate(
formula_composition = ~ 0 + type,
.sample = sample,
.cell_group = cell_group,
cores = 1, verbose = FALSE
) |>
sccomp_test( contrasts = c("typecancer - typehealthy", "typehealthy - typecancer"))Categorical factor (e.g. Bayesian ANOVA)
This is achieved through model comparison with loo. In
the following example, the model with association with factors better
fits the data compared to the baseline model with no factor association.
For comparisons check_outliers must be set to FALSE as the
leave-one-out must work with the same amount of data, while outlier
elimination does not guarantee it.
If elpd_diff is away from zero of > 5
se_diff difference of 5, we are confident that a model is
better than the other reference.
In this case, -79.9 / 11.5 = -6.9, therefore we can conclude that model
one, the one with factor association, is better than model two.
library(loo)
# Fit first model
model_with_factor_association =
seurat_obj |>
sccomp_estimate(
formula_composition = ~ type,
.sample = sample,
.cell_group = cell_group,
inference_method = "hmc",
enable_loo = TRUE
)
# Fit second model
model_without_association =
seurat_obj |>
sccomp_estimate(
formula_composition = ~ 1,
.sample = sample,
.cell_group = cell_group,
inference_method = "hmc",
enable_loo = TRUE
)
# Compare models
loo_compare(
attr(model_with_factor_association, "fit")$loo(),
attr(model_without_association, "fit")$loo()
)Differential variability, binary factor
We can model the cell-group variability also dependent on the type, and so test differences in variability
res =
seurat_obj |>
sccomp_estimate(
formula_composition = ~ type,
formula_variability = ~ type,
.sample = sample,
.cell_group = cell_group,
cores = 1, verbose = FALSE
)
resPlot 1D significance plot
plots = res |> sccomp_test() |> plot()
plots$credible_intervals_1DPlot 2D significance plot Data points are cell groups. Error bars are the 95% credible interval. The dashed lines represent the default threshold fold change for which the probabilities (c_pH0, v_pH0) are calculated. pH0 of 0 represent the rejection of the null hypothesis that no effect is observed.
This plot is provided only if differential variability has been
tested. The differential variability estimates are reliable only if the
linear association between mean and variability for
(intercept) (left-hand side facet) is satisfied. A
scatterplot (besides the Intercept) is provided for each category of
interest. For each category of interest, the composition and variability
effects should be generally uncorrelated.
plots$credible_intervals_2DSuggested settings
For single-cell RNA sequencing
We recommend setting
bimodal_mean_variability_association = TRUE. The
bimodality of the mean-variability association can be confirmed from the
plots$credible_intervals_2D (see below).
For CyTOF and microbiome data
We recommend setting
bimodal_mean_variability_association = FALSE
(Default).
Visualisation of the MCMC chains from the posterior distribution
It is possible to directly evaluate the posterior distribution. In this example, we plot the Monte Carlo chain for the slope parameter of the first cell type. We can see that it has converged and is negative with probability 1.
library(cmdstanr)
library(posterior)
library(bayesplot)
# Assuming res contains the fit object from cmdstanr
fit <- res |> attr("fit")
# Extract draws for 'beta[2,1]'
draws <- as_draws_array(fit$draws("beta[2,1]"))
# Create a traceplot for 'beta[2,1]'
mcmc_trace(draws, pars = "beta[2,1]")## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.5 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=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [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: UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] SeuratObject_5.0.2 sp_2.1-4 tidyr_1.3.1 forcats_1.0.0
## [5] ggplot2_3.5.1 sccomp_1.9.2 dplyr_1.1.4
##
## loaded via a namespace (and not attached):
## [1] rlang_1.1.4 magrittr_2.0.3
## [3] matrixStats_1.4.1 compiler_4.4.1
## [5] systemfonts_1.1.0 callr_3.7.6
## [7] vctrs_0.6.5 stringr_1.5.1
## [9] pkgconfig_2.0.3 crayon_1.5.3
## [11] fastmap_1.2.0 backports_1.5.0
## [13] XVector_0.45.0 utf8_1.2.4
## [15] cmdstanr_0.8.0 rmarkdown_2.28
## [17] tzdb_0.4.0 UCSC.utils_1.1.0
## [19] ps_1.8.0 ragg_1.3.3
## [21] purrr_1.0.2 xfun_0.48
## [23] zlibbioc_1.51.1 cachem_1.1.0
## [25] GenomeInfoDb_1.41.2 jsonlite_1.8.9
## [27] DelayedArray_0.31.14 parallel_4.4.1
## [29] R6_2.5.1 bslib_0.8.0
## [31] stringi_1.8.4 parallelly_1.38.0
## [33] GenomicRanges_1.57.2 jquerylib_0.1.4
## [35] Rcpp_1.0.13 SummarizedExperiment_1.35.4
## [37] knitr_1.48 future.apply_1.11.2
## [39] readr_2.1.5 IRanges_2.39.2
## [41] Matrix_1.7-0 tidyselect_1.2.1
## [43] abind_1.4-8 yaml_2.3.10
## [45] codetools_0.2-20 processx_3.8.4
## [47] listenv_0.9.1 lattice_0.22-6
## [49] tibble_3.2.1 Biobase_2.65.1
## [51] withr_3.0.1 posterior_1.6.0
## [53] evaluate_1.0.1 future_1.34.0
## [55] desc_1.4.3 pillar_1.9.0
## [57] MatrixGenerics_1.17.0 tensorA_0.36.2.1
## [59] checkmate_2.3.2 stats4_4.4.1
## [61] distributional_0.5.0 generics_0.1.3
## [63] S4Vectors_0.43.2 hms_1.1.3
## [65] munsell_0.5.1 scales_1.3.0
## [67] globals_0.16.3 glue_1.8.0
## [69] tools_4.4.1 dotCall64_1.2
## [71] fs_1.6.4 grid_4.4.1
## [73] instantiate_0.2.3.9000 colorspace_2.1-1
## [75] SingleCellExperiment_1.27.2 GenomeInfoDbData_1.2.13
## [77] patchwork_1.3.0 cli_3.6.3
## [79] textshaping_0.4.0 spam_2.11-0
## [81] fansi_1.0.6 S4Arrays_1.5.10
## [83] gtable_0.3.5 sass_0.4.9
## [85] digest_0.6.37 progressr_0.14.0
## [87] BiocGenerics_0.51.3 SparseArray_1.5.44
## [89] ggrepel_0.9.6 htmlwidgets_1.6.4
## [91] farver_2.1.2 htmltools_0.5.8.1
## [93] pkgdown_2.1.1 lifecycle_1.0.4
## [95] prettydoc_0.4.1 httr_1.4.7