sctransform

Pure-python implementation of R sctransform — variance-stabilizing transformation for single-cell UMI count data. No R, no rpy2 required.

Implements both model flavors of R sctransform:

• v2 model (vst.flavor="v2" in R): negative-binomial GLM with library-size offset (slope fixed at log(10)), per-gene intercept + overdispersion estimation, Poisson-gene exclusion.
• v1 model (original sctransform): per-gene Poisson GLM with free slope on log10_umi, theta via a port of MASS::theta.ml.

Both use kernel-smoothed regularization of parameters — all validated end-to-end against the R reference implementation (with glmGamPoi).

Why

Existing python ports either wrap R via rpy2 (defeats the purpose) or delegate the v2 model to R's glmGamPoi. This package re-implements everything natively and numerically matches R:

See benchmark/ for the full executed notebook and reproduction scripts.

Installation

pip install sctransform
# or with the scanpy interface:
pip install sctransform[scanpy]

Quickstart

import scanpy as sc
from sctransform import SCTransform

# adata with raw counts in adata.layers["counts"]
residuals = SCTransform(adata, layer="counts", vst_flavor="v2", var_features_n=3000)

# downstream, Seurat-style: PCA on Pearson residuals of variable features
adata.obsm["X_pca"] = sc.pp.pca(residuals.values, n_comps=30)
sc.pp.neighbors(adata, use_rep="X_pca")
sc.tl.umap(adata)

Low-level API (R-like, genes x cells matrix in, full vst output out):

from sctransform import vst
vst_out = vst(umi_matrix, gene_names=genes, vst_flavor="v2")
vst_out["y"]                    # residuals (genes x cells DataFrame)
vst_out["model_parameters_fit"] # regularized per-gene theta / intercept / slope
vst_out["gene_attr"]            # residual_variance etc. (for HVG ranking)

Scope

• v2 model (glmGamPoi_offset): fixed slope, intercept + theta learned, exclude_poisson
• v1 model (poisson/theta_ml, free slope) — exact port of qpois_reg IRLS + MASS::theta.ml
• Pearson & deviance residuals, min_variance ("umi_median"), residual clipping
• Cell/gene subsampling (density-based gene sampling, as in R)
• batch_var support — roadmap
• corrected counts (return_corrected_umi) — roadmap

Relationship to Seurat

SCTransform(adata, vst_flavor="v2") replaces NormalizeData + FindVariableFeatures + ScaleData: use the returned residuals as scale.data for PCA (as shown above). HVG flags and regularized model parameters are stored in adata.var / adata.uns["sct"].

Tutorials

• tutorials/pbmc3k_sctransform_vs_lognormalize.ipynb - end-to-end pbmc3k analysis, side-by-side vs the traditional log-normalize workflow (variance-stabilization QC, UMAP, cluster agreement, markers, depth effects).

License & attribution

GPL-3.0. Contains code ported from R (GPL-2+, R Core Team) and algorithms re-implemented from R sctransform (GPL-3, Christoph Hafemeister / Rahul Satija lab).

If you use this package, please cite the original publications:

• Hafemeister & Satija, Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression, Genome Biology 2019
• Choudhary & Satija, Comparison and evaluation of statistical error models for scRNA-seq, Genome Biology 2022
• Ahlmann-Eltze & Huber, glmGamPoi: Fitting Gamma-Poisson Generalized Linear Models on Single Cell Count Data, Bioinformatics 2021