Lifecycle:maturing R build status

Cellular omics such as single-cell genomics, proteomics, and microbiomics allow the characterization of tissue and microbial community composition, which can be compared between conditions to identify biological drivers. This strategy has been critical to unveiling markers of disease progression in conditions such as cancer and pathogen infections.

For cellular omic data, no method for differential variability analysis exists, and methods for differential composition analysis only take a few fundamental data properties into account. Here we introduce sccomp, a generalised method for differential composition and variability analyses capable of jointly modelling data count distribution, compositionality, group-specific variability, and proportion mean-variability association, while being robust to outliers.

sccomp is an extensive analysis framework that allows realistic data simulation and cross-study knowledge transfer. We demonstrate that mean-variability association is ubiquitous across technologies, highlighting the inadequacy of the very popular Dirichlet-multinomial modeling and providing essential principles for differential variability analysis.

0.0.1 Comparison with other methods

Method Year Model I II III IV V VI
sccomp 2023 Sum-constrained Beta-binomial
scCODA 2021 Dirichlet-multinomial
quasi-binom. 2021 Quasi-binomial
rlm 2021 Robust-log-linear
propeller 2021 Logit-linear + limma
ANCOM-BC 2020 Log-linear
corncob 2020 Beta-binomial
scDC 2019 Log-linear
dmbvs 2017 Dirichlet-multinomial
MixMC 2016 Zero-inflated Log-linear
ALDEx2 2014 Dirichlet-multinomial

0.0.2 Cite

Mangiola, Stefano, Alexandra J. Roth-Schulze, Marie Trussart, Enrique Zozaya-Valdés, Mengyao Ma, Zijie Gao, Alan F. Rubin, Terence P. Speed, Heejung Shim, and Anthony T. Papenfuss. 2023. “Sccomp: Robust Differential Composition and Variability Analysis for Single-Cell Data.” Proceedings of the National Academy of Sciences of the United States of America 120 (33): e2203828120. https://doi.org/10.1073/pnas.2203828120 PNAS - sccomp: Robust differential composition and variability analysis for single-cell data

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

  1. R installation of sccomp
  2. R installation of cmdstanr
  3. cmdstanr call to cmdstan installation

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

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

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 = 
  sce_obj |>
  sccomp_estimate( 
    formula_composition = ~ type, 
    .sample =  sample, 
    .cell_group = cell_group, 
    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, 
    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 input formula_composition and formula_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_result

2.2 Summary 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() 

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

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

res

2.5 Random 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.

2.5.1 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 namespace: SeuratObject