Lifecycle:maturing R build status

Sccomp is a generalised method for differential composition and variability analyses.

0.1 Characteristics

1 Installation

Bioconductor

if (!requireNamespace("BiocManager")) install.packages("BiocManager")
BiocManager::install("sccomp")

Github

devtools::install_github("stemangiola/sccomp")

2 Analysis

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.

2.1 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).

2.1.1 From Seurat, SingleCellExperiment, metadata objects

sccomp_result = 
  single_cell_object |>
  sccomp_estimate( 
    formula_composition = ~ type, 
    .sample =  sample, 
    .cell_group = cell_group, 
    bimodal_mean_variability_association = TRUE,
    cores = 1 
  ) |> 
  sccomp_remove_outliers(cores = 1) |> # Optional
  sccomp_test()

2.1.2 From counts

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.000557 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 5.57 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 
## Further attempt with Variational Bayes: Error in eval(expr, envir, enclos): stan::variational::normal_meanfield::calc_grad: The number of dropped evaluations has reached its maximum amount (10). Your model may be either severely ill-conditioned or misspecified.
## 
## 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.000544 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 5.44 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: Iteration: 250 / 250 [100%]  (Adaptation)
## Chain 1: Success! Found best value [eta = 0.1].
## Chain 1: 
## Chain 1: Begin stochastic gradient ascent.
## Chain 1:   iter             ELBO   delta_ELBO_mean   delta_ELBO_med   notes 
## Chain 1:    100        -4848.007             1.000            1.000
## Chain 1:    200        -4330.974             0.560            1.000
## Chain 1:    300        -4120.439             0.390            0.119
## Chain 1:    400        -4021.652             0.299            0.119
## Chain 1:    500        -3946.134             0.243            0.051
## Chain 1:    600        -3887.462             0.205            0.051
## Chain 1:    700        -3844.630             0.177            0.025
## Chain 1:    800        -3813.049             0.156            0.025
## Chain 1:    900        -3787.362             0.139            0.019
## Chain 1:   1000        -3770.283             0.126            0.019
## Chain 1:   1100        -3751.395             0.027            0.015
## Chain 1:   1200        -3738.448             0.015            0.011
## Chain 1:   1300        -3729.826             0.010            0.008   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.0006 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 6 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        -3735.308             1.000            1.000
## Chain 1:    200        -3536.380             0.528            1.000
## Chain 1:    300        -3512.206             0.354            0.056
## Chain 1:    400        -3504.972             0.266            0.056
## Chain 1:    500        -3505.723             0.213            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.000591 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 5.91 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        -9408.679             1.000            1.000
## Chain 1:    200        -4492.127             1.047            1.094
## Chain 1:    300        -3733.221             0.766            1.000
## Chain 1:    400        -3690.838             0.577            1.000
## Chain 1:    500        -3683.336             0.462            0.203
## Chain 1:    600        -3681.033             0.385            0.203
## Chain 1:    700        -3680.346             0.330            0.011
## Chain 1:    800        -3676.185             0.289            0.011
## Chain 1:    900        -3672.408             0.257            0.002   MEDIAN ELBO CONVERGED
## Chain 1: 
## Chain 1: Drawing a sample of size 1000 from the approximate posterior... 
## Chain 1: COMPLETED.

2.2 Summary plots

plots = sccomp_result |> plot() 
## 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.

plots$boxplot
## [[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.

plots$credible_intervals_1D

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.

plots$credible_intervals_2D

2.3 Contrasts

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.000476 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 4.76 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        -3293.312             1.000            1.000
## Chain 1:    200        -3162.575             0.521            1.000
## Chain 1:    300        -3162.088             0.347            0.041
## Chain 1:    400        -3159.455             0.261            0.041
## Chain 1:    500        -3160.147             0.209            0.001   MEDIAN ELBO CONVERGED
## Chain 1: 
## Chain 1: Drawing a sample of size 1000 from the approximate posterior... 
## Chain 1: COMPLETED.

2.4 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, 
    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()
)

2.5 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,
    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.000467 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 4.67 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        -3269.338             1.000            1.000
## Chain 1:    200        -3179.795             0.514            1.000
## Chain 1:    300        -3174.425             0.343            0.028
## Chain 1:    400        -3175.979             0.258            0.028
## Chain 1:    500        -3165.542             0.207            0.003   MEDIAN ELBO CONVERGED
## Chain 1: 
## Chain 1: Drawing a sample of size 1000 from the approximate posterior... 
## Chain 1: COMPLETED.
res

3 Suggested settings

3.1 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).

3.2 For CyTOF and microbiome data

We recommend setting bimodal_mean_variability_association = FALSE (Default).

3.3 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.

res %>% attr("fit") %>% rstan::traceplot("beta[2,1]")

Plot 1D significance plot

plots = res |> sccomp_test() |> plot()
## Joining with `by = join_by(cell_group, sample)`
## Joining with `by = join_by(cell_group, type)`
plots$credible_intervals_1D

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.

plots$credible_intervals_2D

sessionInfo()
## R version 4.4.0 RC (2024-04-16 r86468)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.4 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.20-bioc/R/lib/libRblas.so 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.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] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] rstan_2.32.6       StanHeaders_2.32.7 tidyr_1.3.1        forcats_1.0.0     
## [5] ggplot2_3.5.1      sccomp_1.9.0       dplyr_1.1.4       
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.2.1            farver_2.1.1               
##  [3] loo_2.7.0                   fastmap_1.1.1              
##  [5] SingleCellExperiment_1.27.0 digest_0.6.35              
##  [7] dotCall64_1.1-1             lifecycle_1.0.4            
##  [9] SeuratObject_5.0.1          magrittr_2.0.3             
## [11] compiler_4.4.0              rlang_1.1.3                
## [13] sass_0.4.9                  tools_4.4.0                
## [15] utf8_1.2.4                  yaml_2.3.8                 
## [17] knitr_1.46                  labeling_0.4.3             
## [19] S4Arrays_1.5.0              curl_5.2.1                 
## [21] sp_2.1-4                    pkgbuild_1.4.4             
## [23] DelayedArray_0.31.0         RColorBrewer_1.1-3         
## [25] abind_1.4-5                 withr_3.0.0                
## [27] purrr_1.0.2                 BiocGenerics_0.51.0        
## [29] grid_4.4.0                  stats4_4.4.0               
## [31] fansi_1.0.6                 colorspace_2.1-0           
## [33] future_1.33.2               inline_0.3.19              
## [35] progressr_0.14.0            globals_0.16.3             
## [37] scales_1.3.0                SummarizedExperiment_1.35.0
## [39] cli_3.6.2                   rmarkdown_2.26             
## [41] crayon_1.5.2                generics_0.1.3             
## [43] RcppParallel_5.1.7          future.apply_1.11.2        
## [45] httr_1.4.7                  tzdb_0.4.0                 
## [47] cachem_1.0.8                stringr_1.5.1              
## [49] zlibbioc_1.51.0             parallel_4.4.0             
## [51] XVector_0.45.0              matrixStats_1.3.0          
## [53] vctrs_0.6.5                 V8_4.4.2                   
## [55] boot_1.3-30                 Matrix_1.7-0               
## [57] jsonlite_1.8.8              IRanges_2.39.0             
## [59] hms_1.1.3                   patchwork_1.2.0            
## [61] S4Vectors_0.43.0            ggrepel_0.9.5              
## [63] listenv_0.9.1               jquerylib_0.1.4            
## [65] glue_1.7.0                  parallelly_1.37.1          
## [67] spam_2.10-0                 codetools_0.2-20           
## [69] stringi_1.8.3               gtable_0.3.5               
## [71] QuickJSR_1.1.3              GenomeInfoDb_1.41.0        
## [73] GenomicRanges_1.57.0        UCSC.utils_1.1.0           
## [75] munsell_0.5.1               tibble_3.2.1               
## [77] pillar_1.9.0                htmltools_0.5.8.1          
## [79] GenomeInfoDbData_1.2.12     R6_2.5.1                   
## [81] evaluate_0.23               lattice_0.22-6             
## [83] Biobase_2.65.0              highr_0.10                 
## [85] readr_2.1.5                 bslib_0.7.0                
## [87] rstantools_2.4.0            Rcpp_1.0.12                
## [89] gridExtra_2.3               SparseArray_1.5.0          
## [91] xfun_0.43                   MatrixGenerics_1.17.0      
## [93] prettydoc_0.4.1             pkgconfig_2.0.3