sparse-count-pca

Read the documentation →

Single-cell transcriptomics data reaches the user as a large, very sparse matrix of integer counts. Many normalization methods assign nonzero values to observed zeros, making the transformed matrix dense and potentially increasing entry storage by an order of magnitude or more.

Most workflows do not need every normalized entry at once, because their primary consumer is PCA. sparse-count-pca represents supported count transforms exactly as a sparse matrix plus a small low-rank term. This factorization provides efficient matrix-vector and matrix-matrix products without materializing the full transformed matrix. The package uses those products to compute PCA of the column-centered transform through SciPy's LinearOperator interface and truncated SVD.

The implementation is pure Python and is tested against independent dense-matrix implementations and pinned external reference outputs. It covers PCA of GLM model residuals, size-factor-normalized log1p counts, shifted centered log-ratio (CLR) coordinates, classical correspondence analysis, and several more experimental options.

Installation

Python 3.11 or newer is required.

pip install sparse-count-pca

Quick start

The primary API writes Scanpy-compatible PCA outputs to an AnnData object:

import sparse_count_pca as scp

scp.residual_pca(
    adata,
    layer="counts",
    n_comps=50,
    model="poisson",
    residual="pearson",
)

Cell embeddings/scores are written to adata.obsm["X_pca"], component vectors to adata.varm["PCs"], and variance statistics and parameters to adata.uns["pca"]. The full matrix of normalized values is not stored, by design. For plotting those values and other uses, see normalized values and transform reuse.

Pass raw counts: finite, nonnegative, and integer-valued. Floating-point counts are accepted within 1e-8 of an integer, and check_values=False skips that check. Remove cells with no counts first, for example with sc.pp.filter_cells(adata, min_counts=1). See input validation.

Run a complete example with simulated data

Supported transforms

Compare transform assumptions, formulas, provenance, and API maturity

Already using Scanpy? See coming from Scanpy for what each call maps to, including the normalize_total + log1p + pca pipeline.

scaled_nb is a package-specific name for a negative-binomial model that scales overdispersion inversely with cell depth, as in the sSeq model from Yu, Huber, and Vitek (2013). It is not an implementation of SCTransform; its Pearson residuals agree with SCTransform's only under conditions that ordinary data does not meet. See the scaled-NB null model.

There is also an experimental scaled-NB extension to correspondence analysis; see the correspondence-analysis guide.

How it works

Supported dense transforms have an exact representation,

transformed matrix = sparse matrix + low-rank baseline,

which allows efficient factored matrix multiplication without materializing the actual matrix. Selecting PCA variables preserves this form, and column centering adds one rank-one term. This package uses scipy.sparse.linalg.LinearOperator to expose matrix-vector and matrix-matrix products to scipy.sparse.linalg.svds, using the ARPACK solver. See sparse plus low rank for more details.

Performance

Residual PCA never stores the dense residual matrix, so its memory grows with the number of stored counts rather than with cells × genes. Scanpy's normalized log1p is also sparse, so there the two packages perform similarly.

On 22,111 cells × 29,157 genes with 72 million stored counts, pooled from 10x's 5k Human PBMC 3' GEM-X donors 1, 2, 3, and 4, sparse-count-pca 1.0.0 and Scanpy 1.12.4 computed 30 components in float64 on one thread of a Linux server:

Figures here and below are medians of three runs. Memory is the process's peak resident memory, including the loaded counts. Residuals are clipped at Scanpy's ±√n_cells threshold in both packages.

To measure scaling, each package selected 5,000 highly variable genes by Pearson residual variance and computed PCA on them, in float32 with a 26 GiB memory limit. Inputs are random subsets of 10x's 1.3 million brain cells from E18 mice; quality filtering removed about 0.3% of each subset's cells.

API structure

• One-step AnnData API. Functions such as residual_pca, shifted_clr_pca, and correspondence_analysis write results directly into an AnnData object, as in the quick start above. Counts may come from .X or a layer. See AnnData workflows.

• One-step matrix API. The corresponding _matrix functions take dense or sparse matrices outside AnnData and return a result object:

  result = scp.residual_pca_matrix(counts, ...)
  
  result.scores       # observations by components
  result.components   # components by variables
  result.components.T # variables by components
  

• Two-step transform API. transform fits a normalization once, so normalized values can be plotted or inspected in bounded slices, per-gene variances can rank highly variable genes, or several PCA masks can share one fitted state. See normalized values and transform reuse.

Important distinctions

• A gene mask usually selects PCA columns, not normalization inputs. Slice the count matrix first if excluded genes should not affect cell totals, gene proportions, or CLR row means. Correspondence analysis is the exception: its mask defines the contingency table, so margins are recomputed after masking. See normalization, masking, and centering.
• Scaled-NB residuals are not a drop-in SCTransform v2 implementation. Do not expect default outputs to match. See the compatibility reference for the conditions under which equality can be tested.
• Residuals are clipped by default, as Seurat and Scanpy do, at Seurat's threshold ±sqrt(n_obs / 30) (clip="seurat"), so a few extreme values cannot dominate PCA. This symmetric clipping can make memory use depend on the data values, not just the matrix's size and sparsity, so the package caps that growth and raises RuntimeError beyond it. Pass clip_max_nnz_ratio=None to remove the cap. See clipping and precision for details and alternatives.
• Keep the default dtype="float64" unless memory requires otherwise; see calculation dtype for what float32 changes.

Related work

The cleartools projects provide sparse Rust-backed tooling for shifted CLR (PFlog) workflows: the Python scclr package, built on the runorm normalization crate. See the compatibility reference for the scope of the relationship.

10XGenomics/scan-rs is a Rust library used by Cell Ranger and contains similar matrix-representation and normalization machinery.

Citation

If you use sparse-count-pca in published work, please cite the version you used. Citation metadata is in CITATION.cff; the Cite this repository button on the GitHub page formats it as APA or BibTeX. Each GitHub release is archived on Zenodo with its own DOI. Please also cite the methods behind the transforms you use; the transform catalog lists each one's origin.

Development

uv sync --extra test --extra docs
uv run python -m pytest
uv run ruff check .
uv run ruff format --check .
uv run mkdocs build --strict

See CONTRIBUTING.md for project invariants and documentation expectations.