residual_transform: Residual-based Variance Stabilizing Transformation

View source: R/residual_transform.R

residual_transformR Documentation

Residual-based Variance Stabilizing Transformation


Fit an intercept Gamma-Poisson model that corrects for sequencing depth and return the residuals as variance stabilized results for further downstream application, for which no proper count-based method exist or is performant enough (e.g., clustering, dimensionality reduction).


  residual_type = c("randomized_quantile", "pearson", "analytic_pearson"),
  clipping = FALSE,
  overdispersion = 0.05,
  size_factors = TRUE,
  offset_model = TRUE,
  overdispersion_shrinkage = TRUE,
  ridge_penalty = 2,
  on_disk = NULL,
  return_fit = FALSE,
  verbose = FALSE,



any matrix-like object (e.g. matrix, dgCMatrix, DelayedArray, HDF5Matrix) with one column per sample and row per gene. It can also be an object of type glmGamPoi, in which case it is directly used to calculate the variance-stabilized values.


a string that specifies what kind of residual is returned as variance stabilized-value.


The discrete nature of count distribution stops simple transformations from obtaining a truly standard normal residuals. The trick of of quantile randomized residuals is to match the cumulative density function of the Gamma-Poisson and the Normal distribution. Due to the discrete nature of Gamma-Poisson distribution, a count does not correspond to a single quantile of the Normal distribution, but to a range of possible value. This is resolved by randomly choosing one of the mapping values from the Normal distribution as the residual. This ensures perfectly normal distributed residuals, for the cost of introducing randomness. More details are available in the documentation of statmod::qresiduals() and the corresponding publication by Dunn and Smyth (1996).


The Pearson residuals are defined as res = (y - m) / sqrt(m + m^2 * theta).


Similar to the method above, however, instead of estimating m using a GLM model fit, m is approximated by m_ij = (\sum_j y_{ij}) (\sum_i y_{ij}) / (\sum_{i,j} y_{ij}). For all details, see Lause et al. (2021). Note that overdispersion_shrinkage and ridge_penalty are ignored when fitting analytic Pearson residuals.

The two above options are the most common choices, however you can use any residual_type supported by glmGamPoi::residuals.glmGamPoi(). Default: "randomized_quantile"


a single boolean or numeric value specifying that all residuals are in the range ⁠[-clipping, +clipping]⁠. If clipping = TRUE, we use the default of clipping = sqrt(ncol(data)) which is the default behavior for sctransform. Default: FALSE, which means no clipping is applied.


the simplest count model is the Poisson model. However, the Poisson model assumes that variance = mean. For many applications this is too rigid and the Gamma-Poisson allows a more flexible mean-variance relation (variance = mean + mean^2 * overdispersion).
overdispersion can either be

  • a single boolean that indicates if an overdispersion is estimated for each gene.

  • a numeric vector of length nrow(data) fixing the overdispersion to those values.

  • the string "global" to indicate that one dispersion is fit across all genes.

Note that overdispersion = 0 and overdispersion = FALSE are equivalent and both reduce the Gamma-Poisson to the classical Poisson model. Default: 0.05 which is roughly the overdispersion observed on ostensibly homogeneous cell lines.


boolean to decide if \beta_1 in y = \beta_0 + \beta_1 log(sf), is set to 1 (i.e., treating the log of the size factors as an offset ) or is estimated per gene. From a theoretical point, it should be fine to treat \beta_1 as an offset, because a cell that is twice as big, should have twice as many counts per gene (without any gene-specific effects). However, sctransform suggested that it would be advantageous to nonetheless estimate \beta_0 as it may counter data artifacts. On the other side, Lause et al. (2020) demonstrated that the estimating \beta_0 and \beta_1 together can be difficult. If you still want to fit sctransform's model, you can set the ridge_penalty argument to a non-zero value, which shrinks \beta_1 towards 1 and resolves the degeneracy.
Default: TRUE.

overdispersion_shrinkage, size_factors

arguments that are passed to the underlying call to glmGamPoi::glm_gp(). Default for each: TRUE.


another argument that is passed to glmGamPoi::glm_gp(). It is ignored if offset_model = TRUE. Default: 2.


a boolean that indicates if the dataset is loaded into memory or if it is kept on disk to reduce the memory usage. Processing in memory can be significantly faster than on disk. Default: NULL which means that the data is only processed in memory if data is an in-memory data structure.


boolean to decide if the matrix of residuals is returned directly (return_fit = FALSE) or if in addition the glmGamPoi-fit is returned (return_fit = TRUE) . Default: FALSE.


boolean that decides if information about the individual steps are printed. Default: FALSE


additional parameters passed to glmGamPoi::glm_gp().


Internally, this method uses the glmGamPoi package. The function goes through the following steps

  1. fit model using glmGamPoi::glm_gp()

  2. plug in the trended overdispersion estimates

  3. call glmGamPoi::residuals.glmGamPoi() to calculate the residuals.


a matrix (or a vector if the input is a vector) with the transformed values. If return_fit = TRUE, a list is returned with two elements: fit and Residuals.


Ahlmann-Eltze, Constantin, and Wolfgang Huber. "glmGamPoi: Fitting Gamma-Poisson Generalized Linear Models on Single Cell Count Data." Bioinformatics (2020)

Dunn, Peter K., and Gordon K. Smyth. "Randomized quantile residuals." Journal of Computational and Graphical Statistics 5.3 (1996): 236-244.

Hafemeister, Christoph, and Rahul Satija. "Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression." Genome biology 20.1 (2019): 1-15.

Hafemeister, Christoph, and Rahul Satija. "Analyzing scRNA-seq data with the sctransform and offset models." (2020)

Lause, Jan, Philipp Berens, and Dmitry Kobak. "Analytic Pearson residuals for normalization of single-cell RNA-seq UMI data." Genome Biology (2021).

See Also

glmGamPoi::glm_gp(), glmGamPoi::residuals.glmGamPoi(), sctransform::vst(), statmod::qresiduals()


  # Load a single cell dataset
  sce <- TENxPBMCData::TENxPBMCData("pbmc4k")
  # Reduce size for this example
  sce_red <- sce[sample(which(rowSums2(counts(sce)) > 0), 1000),
                 sample(ncol(sce), 200)]
  counts(sce_red) <- as.matrix(counts(sce_red))

  # Residual Based Variance Stabilizing Transformation
  rq <- residual_transform(sce_red, residual_type = "randomized_quantile",
                           verbose = TRUE)
  pearson <- residual_transform(sce_red, residual_type = "pearson", verbose = TRUE)

  # Plot first two principal components
  pearson_pca <- prcomp(t(pearson), rank. = 2)
  rq_pca <- prcomp(t(rq), rank. = 2)
  plot(rq_pca$x, asp = 1)
  points(pearson_pca$x, col = "red")

const-ae/transformGamPoi documentation built on April 14, 2023, 11:33 p.m.