residual_transform: Residual-based Variance Stabilizing Transformation
In const-ae/transformGamPoi: Variance Stabilizing Transformation for Gamma-Poisson Models

View source: R/residual_transform.R

residual_transform

R Documentation

Residual-based Variance Stabilizing Transformation

Description

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

Usage

residual_transform(
  data,
  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,
  ...
)

Arguments

`data`	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.
`residual_type`	a string that specifies what kind of residual is returned as variance stabilized-value. `"randomized_quantile"` 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). `"pearson"` The Pearson residuals are defined as `res = (y - m) / sqrt(m + m^2 * theta)`. `"analytic_pearson"` 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"`
`clipping`	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.
`overdispersion`	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.
`offset_model`	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`.
`ridge_penalty`	another argument that is passed to `glmGamPoi::glm_gp()`. It is ignored if `offset_model = TRUE`. Default: `2`.
`on_disk`	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.
`return_fit`	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`.
`verbose`	boolean that decides if information about the individual steps are printed. Default: `FALSE`
`...`	additional parameters passed to `glmGamPoi::glm_gp()`.

Details

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

fit model using glmGamPoi::glm_gp()
plug in the trended overdispersion estimates
call glmGamPoi::residuals.glmGamPoi() to calculate the residuals.

Value

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.

References

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

Examples

  # Load a single cell dataset
  sce <- TENxPBMCData::TENxPBMCData("pbmc4k")
  # Reduce size for this example
  set.seed(1)
  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.