Sccomp is a generalised method for differential composition and variability analyses.
Bioconductor
if (!requireNamespace("BiocManager")) install.packages("BiocManager")
BiocManager::install("sccomp")
Github
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.
Of the output table, the estimate columns start with the prefix
c_
indicate composition
, or with
v_
indicate variability
(when
formula_variability is set).
sccomp_result =
counts_obj |>
sccomp_estimate(
formula_composition = ~ type,
.sample = sample,
.cell_group = cell_group,
.count = count,
bimodal_mean_variability_association = TRUE,
cores = 1
) |>
sccomp_remove_outliers(cores = 1) |> # Optional
sccomp_test()
## Chain 1: ------------------------------------------------------------
## Chain 1: EXPERIMENTAL ALGORITHM:
## Chain 1: This procedure has not been thoroughly tested and may be unstable
## Chain 1: or buggy. The interface is subject to change.
## Chain 1: ------------------------------------------------------------
## Chain 1:
## Chain 1:
## Chain 1:
## Chain 1: Gradient evaluation took 0.000453 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 4.53 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Begin eta adaptation.
## Chain 1: Iteration: 1 / 250 [ 0%] (Adaptation)
## Chain 1: Iteration: 50 / 250 [ 20%] (Adaptation)
## Chain 1: Iteration: 100 / 250 [ 40%] (Adaptation)
## Chain 1: Iteration: 150 / 250 [ 60%] (Adaptation)
## Chain 1: Iteration: 200 / 250 [ 80%] (Adaptation)
## Chain 1: Success! Found best value [eta = 1] earlier than expected.
## Chain 1:
## Chain 1: Begin stochastic gradient ascent.
## Chain 1: iter ELBO delta_ELBO_mean delta_ELBO_med notes
## Chain 1: 100 -4138.208 1.000 1.000
## Chain 1: 200 -3725.386 0.555 1.000
## Chain 1: 300 -3712.067 0.371 0.111
## Chain 1: 400 -3709.968 0.279 0.111
## Chain 1: 500 -3708.724 0.223 0.004 MEDIAN ELBO CONVERGED
## Chain 1:
## Chain 1: Drawing a sample of size 1000 from the approximate posterior...
## Chain 1: COMPLETED.
## Chain 1: ------------------------------------------------------------
## Chain 1: EXPERIMENTAL ALGORITHM:
## Chain 1: This procedure has not been thoroughly tested and may be unstable
## Chain 1: or buggy. The interface is subject to change.
## Chain 1: ------------------------------------------------------------
## Chain 1:
## Chain 1:
## Chain 1:
## Chain 1: Gradient evaluation took 0.000304 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 3.04 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Begin eta adaptation.
## Chain 1: Iteration: 1 / 250 [ 0%] (Adaptation)
## Chain 1: Iteration: 50 / 250 [ 20%] (Adaptation)
## Chain 1: Iteration: 100 / 250 [ 40%] (Adaptation)
## Chain 1: Iteration: 150 / 250 [ 60%] (Adaptation)
## Chain 1: Iteration: 200 / 250 [ 80%] (Adaptation)
## Chain 1: Success! Found best value [eta = 1] earlier than expected.
## Chain 1:
## Chain 1: Begin stochastic gradient ascent.
## Chain 1: iter ELBO delta_ELBO_mean delta_ELBO_med notes
## Chain 1: 100 -3841.936 1.000 1.000
## Chain 1: 200 -3370.631 0.570 1.000
## Chain 1: 300 -3347.241 0.382 0.140
## Chain 1: 400 -3336.206 0.288 0.140
## Chain 1: 500 -3329.117 0.230 0.007 MEDIAN ELBO CONVERGED
## Chain 1:
## Chain 1: Drawing a sample of size 1000 from the approximate posterior...
## Chain 1: COMPLETED.
## Chain 1: ------------------------------------------------------------
## Chain 1: EXPERIMENTAL ALGORITHM:
## Chain 1: This procedure has not been thoroughly tested and may be unstable
## Chain 1: or buggy. The interface is subject to change.
## Chain 1: ------------------------------------------------------------
## Chain 1:
## Chain 1:
## Chain 1:
## Chain 1: Gradient evaluation took 0.000297 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.97 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Begin eta adaptation.
## Chain 1: Iteration: 1 / 250 [ 0%] (Adaptation)
## Chain 1: Iteration: 50 / 250 [ 20%] (Adaptation)
## Chain 1: Iteration: 100 / 250 [ 40%] (Adaptation)
## Chain 1: Iteration: 150 / 250 [ 60%] (Adaptation)
## Chain 1: Iteration: 200 / 250 [ 80%] (Adaptation)
## Chain 1: Success! Found best value [eta = 1] earlier than expected.
## Chain 1:
## Chain 1: Begin stochastic gradient ascent.
## Chain 1: iter ELBO delta_ELBO_mean delta_ELBO_med notes
## Chain 1: 100 -4216.054 1.000 1.000
## Chain 1: 200 -3659.747 0.576 1.000
## Chain 1: 300 -3643.644 0.385 0.152
## Chain 1: 400 -3635.401 0.290 0.152
## Chain 1: 500 -3627.537 0.232 0.004 MEDIAN ELBO CONVERGED
## Chain 1:
## Chain 1: Drawing a sample of size 1000 from the approximate posterior...
## Chain 1: COMPLETED.
## Joining with `by = join_by(cell_group, sample)`
## Joining with `by = join_by(cell_group, type)`
A plot of group proportion, faceted by groups. The blue boxplots
represent the posterior predictive check. If the model is likely to be
descriptively adequate to the data, the blue box plot should roughly
overlay with the black box plot, which represents the observed data. The
outliers are coloured in red. A box plot will be returned for every
(discrete) covariate present in formula_composition
. The
colour coding represents the significant associations for composition
and/or variability.
## [[1]]
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 bigger than the minimal effect according to the 95% credible interval. Facets represent the covariates in the model.
We can plot the relationship between abundance and variability. As we can see below, they are positively correlated, you also appreciate that this relationship is by model for single cell RNA sequencing data.
sccomp
models, these 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 for Bayesian
inference.
seurat_obj |>
sccomp_estimate(
formula_composition = ~ 0 + type,
.sample = sample,
.cell_group = cell_group,
bimodal_mean_variability_association = TRUE,
cores = 1
) |>
sccomp_test( contrasts = c("typecancer - typehealthy", "typehealthy - typecancer"))
## Chain 1: ------------------------------------------------------------
## Chain 1: EXPERIMENTAL ALGORITHM:
## Chain 1: This procedure has not been thoroughly tested and may be unstable
## Chain 1: or buggy. The interface is subject to change.
## Chain 1: ------------------------------------------------------------
## Chain 1:
## Chain 1:
## Chain 1:
## Chain 1: Gradient evaluation took 0.000271 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.71 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Begin eta adaptation.
## Chain 1: Iteration: 1 / 250 [ 0%] (Adaptation)
## Chain 1: Iteration: 50 / 250 [ 20%] (Adaptation)
## Chain 1: Iteration: 100 / 250 [ 40%] (Adaptation)
## Chain 1: Iteration: 150 / 250 [ 60%] (Adaptation)
## Chain 1: Iteration: 200 / 250 [ 80%] (Adaptation)
## Chain 1: Success! Found best value [eta = 1] earlier than expected.
## Chain 1:
## Chain 1: Begin stochastic gradient ascent.
## Chain 1: iter ELBO delta_ELBO_mean delta_ELBO_med notes
## Chain 1: 100 -3198.862 1.000 1.000
## Chain 1: 200 -3158.895 0.506 1.000
## Chain 1: 300 -3158.482 0.338 0.013
## Chain 1: 400 -3156.540 0.253 0.013
## Chain 1: 500 -3159.343 0.203 0.001 MEDIAN ELBO CONVERGED
## Chain 1:
## Chain 1: Drawing a sample of size 1000 from the approximate posterior...
## Chain 1: COMPLETED.
## # A tibble: 60 × 18
## cell_group parameter factor c_lower c_effect c_upper c_pH0 c_FDR c_n_eff
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 B immature typecanc… <NA> -1.80 -1.36 -0.901 0 0 NA
## 2 B immature typeheal… <NA> 0.901 1.36 1.80 0 0 NA
## 3 B mem typecanc… <NA> -2.48 -1.87 -1.22 0 0 NA
## 4 B mem typeheal… <NA> 1.22 1.87 2.48 0 0 NA
## 5 CD4 cm S100… typecanc… <NA> -1.41 -1.07 -0.724 0 0 NA
## 6 CD4 cm S100… typeheal… <NA> 0.724 1.07 1.41 0 0 NA
## 7 CD4 cm high… typecanc… <NA> 0.972 1.99 2.96 0 0 NA
## 8 CD4 cm high… typeheal… <NA> -2.96 -1.99 -0.972 0 0 NA
## 9 CD4 cm ribo… typecanc… <NA> 0.150 0.824 1.52 0.0330 0.00671 NA
## 10 CD4 cm ribo… typeheal… <NA> -1.52 -0.824 -0.150 0.0330 0.00671 NA
## # ℹ 50 more rows
## # ℹ 9 more variables: c_R_k_hat <dbl>, v_lower <dbl>, v_effect <dbl>,
## # v_upper <dbl>, v_pH0 <dbl>, v_FDR <dbl>, v_n_eff <dbl>, v_R_k_hat <dbl>,
## # count_data <list>
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,
bimodal_mean_variability_association = TRUE,
cores = 1,
enable_loo = TRUE
)
# Fit second model
model_without_association =
seurat_obj |>
sccomp_estimate(
formula_composition = ~ 1,
.sample = sample,
.cell_group = cell_group,
bimodal_mean_variability_association = TRUE,
cores = 1 ,
enable_loo = TRUE
)
# Compare models
loo_compare(
model_with_factor_association |> attr("fit") |> loo(),
model_without_association |> attr("fit") |> loo()
)
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,
bimodal_mean_variability_association = TRUE,
cores = 1
)
## Chain 1: ------------------------------------------------------------
## Chain 1: EXPERIMENTAL ALGORITHM:
## Chain 1: This procedure has not been thoroughly tested and may be unstable
## Chain 1: or buggy. The interface is subject to change.
## Chain 1: ------------------------------------------------------------
## Chain 1:
## Chain 1:
## Chain 1:
## Chain 1: Gradient evaluation took 0.000274 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.74 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Begin eta adaptation.
## Chain 1: Iteration: 1 / 250 [ 0%] (Adaptation)
## Chain 1: Iteration: 50 / 250 [ 20%] (Adaptation)
## Chain 1: Iteration: 100 / 250 [ 40%] (Adaptation)
## Chain 1: Iteration: 150 / 250 [ 60%] (Adaptation)
## Chain 1: Iteration: 200 / 250 [ 80%] (Adaptation)
## Chain 1: Success! Found best value [eta = 1] earlier than expected.
## Chain 1:
## Chain 1: Begin stochastic gradient ascent.
## Chain 1: iter ELBO delta_ELBO_mean delta_ELBO_med notes
## Chain 1: 100 -3232.395 1.000 1.000
## Chain 1: 200 -3182.012 0.508 1.000
## Chain 1: 300 -3170.789 0.340 0.016
## Chain 1: 400 -3173.175 0.255 0.016
## Chain 1: 500 -3168.396 0.204 0.004 MEDIAN ELBO CONVERGED
## Chain 1:
## Chain 1: Drawing a sample of size 1000 from the approximate posterior...
## Chain 1: COMPLETED.
## # A tibble: 60 × 14
## cell_group parameter factor c_lower c_effect c_upper c_n_eff c_R_k_hat
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 B immature (Interce… <NA> 0.565 0.925 1.25 NaN 3.02
## 2 B immature typeheal… type 0.897 1.41 1.88 NaN 3.02
## 3 B mem (Interce… <NA> -1.24 -0.859 -0.443 NaN 2.93
## 4 B mem typeheal… type 1.49 2.03 2.58 NaN 2.99
## 5 CD4 cm S100A4 (Interce… <NA> 1.42 1.69 2.01 NaN 2.99
## 6 CD4 cm S100A4 typeheal… type 0.492 1.02 1.48 NaN 3.03
## 7 CD4 cm high cyto… (Interce… <NA> -1.13 -0.572 -0.0278 NaN 2.98
## 8 CD4 cm high cyto… typeheal… type -2.34 -1.57 -0.862 NaN 3.02
## 9 CD4 cm ribosome (Interce… <NA> -0.0121 0.326 0.659 NaN 3.04
## 10 CD4 cm ribosome typeheal… type -1.41 -0.917 -0.473 NaN 2.94
## # ℹ 50 more rows
## # ℹ 6 more variables: v_lower <dbl>, v_effect <dbl>, v_upper <dbl>,
## # v_n_eff <dbl>, v_R_k_hat <dbl>, count_data <list>
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).
We recommend setting
bimodal_mean_variability_association = FALSE
(Default).
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.
Plot 1D significance plot
## Joining with `by = join_by(cell_group, sample)`
## Joining with `by = join_by(cell_group, type)`
Plot 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. The for each category of interest, the composition and
variability effects should be generally uncorrelated.
## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04 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
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 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: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] rstan_2.32.6 StanHeaders_2.32.10 tidyr_1.3.1
## [4] forcats_1.0.0 ggplot2_3.5.1 sccomp_1.9.0
## [7] dplyr_1.1.4 prettydoc_0.4.1
##
## loaded via a namespace (and not attached):
## [1] gridExtra_2.3 inline_0.3.19
## [3] rlang_1.1.4 magrittr_2.0.3
## [5] matrixStats_1.3.0 compiler_4.4.1
## [7] loo_2.8.0 vctrs_0.6.5
## [9] stringr_1.5.1 pkgconfig_2.0.3
## [11] crayon_1.5.3 fastmap_1.2.0
## [13] XVector_0.45.0 labeling_0.4.3
## [15] utf8_1.2.4 rmarkdown_2.28
## [17] tzdb_0.4.0 UCSC.utils_1.1.0
## [19] purrr_1.0.2 xfun_0.47
## [21] zlibbioc_1.51.1 cachem_1.1.0
## [23] GenomeInfoDb_1.41.1 jsonlite_1.8.8
## [25] highr_0.11 DelayedArray_0.31.11
## [27] parallel_4.4.1 R6_2.5.1
## [29] RColorBrewer_1.1-3 bslib_0.8.0
## [31] stringi_1.8.4 parallelly_1.38.0
## [33] boot_1.3-31 GenomicRanges_1.57.1
## [35] jquerylib_0.1.4 Rcpp_1.0.13
## [37] SummarizedExperiment_1.35.1 knitr_1.48
## [39] future.apply_1.11.2 readr_2.1.5
## [41] IRanges_2.39.2 Matrix_1.7-0
## [43] tidyselect_1.2.1 abind_1.4-5
## [45] yaml_2.3.10 codetools_0.2-20
## [47] curl_5.2.2 listenv_0.9.1
## [49] pkgbuild_1.4.4 lattice_0.22-6
## [51] tibble_3.2.1 Biobase_2.65.1
## [53] withr_3.0.1 evaluate_0.24.0
## [55] future_1.34.0 RcppParallel_5.1.9
## [57] pillar_1.9.0 MatrixGenerics_1.17.0
## [59] stats4_4.4.1 generics_0.1.3
## [61] sp_2.1-4 S4Vectors_0.43.2
## [63] hms_1.1.3 rstantools_2.4.0
## [65] munsell_0.5.1 scales_1.3.0
## [67] globals_0.16.3 glue_1.7.0
## [69] maketools_1.3.0 tools_4.4.1
## [71] sys_3.4.2 buildtools_1.0.0
## [73] dotCall64_1.1-1 grid_4.4.1
## [75] QuickJSR_1.3.1 colorspace_2.1-1
## [77] SingleCellExperiment_1.27.2 GenomeInfoDbData_1.2.12
## [79] patchwork_1.2.0 cli_3.6.3
## [81] spam_2.10-0 fansi_1.0.6
## [83] S4Arrays_1.5.7 V8_5.0.0
## [85] gtable_0.3.5 sass_0.4.9
## [87] digest_0.6.37 progressr_0.14.0
## [89] BiocGenerics_0.51.0 SparseArray_1.5.31
## [91] ggrepel_0.9.5 SeuratObject_5.0.2
## [93] farver_2.1.2 htmltools_0.5.8.1
## [95] lifecycle_1.0.4 httr_1.4.7