Evaluation#

Abstract#

pyglmGamPoi is a pure Python/C++ reimplementation of sctransform::fit_glmGamPoi_offset, the negative-binomial offset model at the core of Seurat SCTransform v2. We evaluate numerical parity against R on both synthetic count matrices (100–2000 genes × 200–2000 cells) and three real 10x Genomics samples from GSE288946 (~2000 × 2000 SCT step-1 subsamples each).

Across all benchmarks, finite-gene theta estimates agree with R at machine-precision levels (median |Δ| ≈ 2×10-7), intercepts match at ~10:sup:`-9` (100% within 0.01 with shrinkage enabled), and Poisson-boundary theta = Inf assignments match R exactly (463/463 pooled). With OpenMP (8 threads) and a native CSR sparse path, Python completes the real-data step-1 workload in 7.1 s vs R 46.5 s (~**6.4×** faster), while eliminating the R runtime dependency for trackcell SCT v2 pipelines.

Motivation#

Single-cell RNA-seq normalization via SCTransform v2 fits, for every gene, a negative-binomial generalized linear model with a sequencing-depth offset:

\[y_{gc} \sim \mathrm{NB}(\mu_{gc}, \theta_g), \quad \log \mu_{gc} = \beta_g + \log(10^{\mathrm{log10\_umi}_c})\]

Seurat implements this through the R package glmGamPoi (Hafemeister & Satija, 2019). Python pipelines such as trackcell historically bridged to R via Rscript subprocesses, adding environment complexity (conda st, CRAN/Bioconductor versions, TRACKCELL_RSCRIPT overrides) and blocking fully native workflows.

pyglmGamPoi reimplements the glmGamPoi algorithms in C++17 (Newton–Raphson intercept fit, Cox–Reid overdispersion MLE, optional cross-gene overdispersion_shrinkage) and exposes a drop-in API for trackcell:

from pyglmGamPoi import fit_offset_model, is_available

if is_available():
    model_pars = fit_offset_model(umi, regressor_data, gene_index)

This page documents the replacement effect: what R call is replaced, how we measure parity, and what accuracy/speed users should expect.

What is replaced (and what is not)#

Component

pyglmGamPoi

Notes

fit_glmGamPoi_offset

Replaced

Per-gene NB offset fit; primary scope

Cox–Reid overdispersion MLE

Replaced

allow_inf_theta → theta = Inf at Poisson boundary

Cross-gene overdispersion_shrinkage

Replaced

loc_median_fit trend + intercept refit (R glm_gp default)

vst() subsampling (cells/genes step1)

trackcell

RNG differs between R sample() and NumPy

_reg_model_pars ksmooth / outlier logic

trackcell

Kernel smoothing after initial fit

Pearson residuals, HVG, integration

trackcell

Downstream of model_pars

Replacing only the NB fit isolates a well-defined numerical target and avoids conflating RNG subsampling differences with algorithmic ones.

Benchmark design#

Two complementary suites:

Synthetic scaling (scripts/generate_benchmark_report.py)

Poisson-distributed UMI matrices at five sizes, same API as production. Validates that parity holds from smoke tests (100 genes) up to SCT step-1 scale (2000 × 2000).

Real data — Tier B (scripts/run_tcl_benchmark.py)

Three GSE288946 samples (GSM8779707–709), filtered 10x counts aligned to R cells_step1 / genes_step1 exports from Seurat vst(seed = 1448145). This Tier-B protocol (see trackcell tcl_test README) removes subsampling RNG as a confounder so differences reflect the fitter only. Count matrices are passed as CSR sparse (typical 10x input) without densifying.

Reference implementation

Fresh sctransform::fit_glmGamPoi_offset(umi, model_str = "y ~ log_umi", ...) in conda env st (Seurat 4.4 + glmGamPoi). Cached model_pars_step1.csv from vst() is not used as ground truth — R itself differs from fit_glmGamPoi_offset on the intercept scale (~0.5).

Environment

Benchmarks run with OMP_NUM_THREADS=8 unless noted. Regenerate with python scripts/run_all_benchmarks.py.

Metrics

Per gene: theta, Intercept, log_umi (= log(10)). Reported: median absolute difference, Pearson r, fraction within 0.01, Inf-theta match count, wall time.

Results#

Synthetic scaling#

Synthetic scaling benchmark (Python vs R fit_glmGamPoi_offset)#

Genes

Cells

Intercept med |diff|

Theta med |diff|

Intercept r

Theta r

Inf θ match

Speedup

100

200

0.0000

4.11e-07

1.000000

1.000000

0/0

164.39x

200

500

0.0000

4.74e-07

1.000000

1.000000

0/0

76.89x

500

1000

0.0000

5.77e-07

1.000000

1.000000

0/0

36.03x

1000

1500

0.0000

6.63e-07

1.000000

1.000000

0/0

17.33x

2000

2000

0.0000

6.32e-07

1.000000

1.000000

0/0

8.79x

At SCT step-1 scale (2000×2000 synthetic), Python is ~8.8× faster than R with identical theta / Intercept (med |Δ| ≈ 10-11 / 10-7).

Real data (GSE288946)#

Real-data parity (GSE288946 SCT step1, Tier-B aligned subsets)#

Sample

Genes

Cells

Intercept med |diff|

Theta med |diff|

Intercept r

Theta r

Inf θ match

Speedup

GSM8779707

1995

2000

1.55e-09

2.44e-07

1.000000

0.028824

217/217

6.24x

GSM8779708

1997

2000

2.52e-09

2.02e-07

1.000000

-0.000408

188/188

6.95x

GSM8779709

1996

2000

6.19e-09

1.68e-07

1.000000

-0.000520

58/58

6.38x
Pooled across samples (vs R live)#

Metric

Value

Samples / total genes

3 / 5,988

Intercept med |diff| (median across samples)

2.52e-09

Genes with |Δ Intercept| < 0.01 (median %)

100.0%

Theta med |diff| (median across samples)

2.02e-07

Inf θ match (pooled)

463/463

Wall time Python / R (total)

7.1s / 46.5s

Speedup (median)

6.38x

Note

Tier-B alignment: Python and R use identical cells_step1 / genes_step1 subsamples exported from Seurat vst() (seed 1448145). Comparison is against a fresh sctransform::fit_glmGamPoi_offset call (conda env st), not cached model_pars_step1.csv — the latter differs from fit_glmGamPoi_offset by ~0.5 on the intercept scale.

Performance summary#

Speed vs R (fit_glmGamPoi_offset, 8 OpenMP threads)#

Workload

Python

R (conda st)

Speedup

Synthetic 2000×2000

2.0 s

17.5 s

8.8×

Real GSM8779707 (1995×2000)

2.5 s

15.5 s

6.2×

Real 3-sample pooled

7.1 s

46.5 s

6.4×

Implementation milestones (same Tier-B data, GSM8779707):

Configuration

Wall time

vs R

Initial single-thread C++ (dense)

~15 s

~1.0×

  • OpenMP (8 threads)

~2.2 s

~7×

  • CSR sparse path (no toarray())

~2.2 s

~7× (parity unchanged)

The CSR path avoids materializing ~90% zeros in 10x matrices; numerical results match dense after sort_indices() (required when gene/cell slicing leaves unsorted CSR rows).

Representative gene-level agreement (GSM8779707, vs R live)#

Gene

Python θ

R θ

|Δ Intercept|

RAB7A

4.859

4.859

~10-9

RPS7

7.559

7.559

~10-9

TTK (multimodal OD)

0.00774

0.00774

~10-8

AL392046.1 (Poisson)

Inf

Inf

exact

Downstream effect (trackcell integration)#

When Python uses R-exported model_pars_fit (ksmooth applied in R) and computes Pearson residuals, residual correlation vs R reaches r ≈ 1.0 (median |diff| ≈ 10-13 on GSM8779707). This confirms that once model parameters match, downstream SCT quantities are unchanged within floating-point noise.

Discussion#

Accuracy. Finite theta values are effectively identical to R (median |Δ| ~ 2×10-7 on real data). With default cross-gene shrinkage, intercepts match R at machine precision (median |Δ| ~ 2.5×10-9; 100% of genes within 0.01). Shrinkage gene_means now follow R glm_gp first-pass rowMeans(Mu) (rough-dispersion beta), not the post-OD intercept. theta = Inf genes (Poisson boundary under allow_inf_theta = TRUE) match R exactly (217/217, 188/188, 58/58 across the three samples). Pearson r for theta on real data is not informative when many genes share Inf status; use median |diff| instead.

Speed. OpenMP parallelizes the per-gene C++ loop (schedule(dynamic, 16)). The native CSR path reads 10x count matrices in place (implicit zeros still contribute to the likelihood). On GSE288946 step1 (2000×2000, 8 threads), Python is ~6–7× faster than R while preserving identical theta / Intercept parity. Set OMP_NUM_THREADS before importing pyglmGamPoi.

Why not compare to ``vst()`` exports directly? Seurat’s vst() pipeline applies additional steps before/after fit_glmGamPoi_offset; cached model_pars_step1.csv shows intercept offset ~0.5 vs a direct fit_glmGamPoi_offset call even in R. Our comparison target is the function trackcell actually replaces.

Integration path. In trackcell fit_offset_model, priority becomes: (1) pyglmGamPoi native, (2) R fit_glmGamPoi_offset, (3) statsmodels nb_offset fallback. Users install with pip install pyglmGamPoi — no conda env st required for the NB fit step.

Conclusion#

pyglmGamPoi provides a drop-in, R-free replacement for sctransform::fit_glmGamPoi_offset with demonstrated numerical parity on synthetic and real scRNA-seq data, and ~6–7× wall-clock speedup on production-scale SCT step-1 workloads (8 threads). For trackcell SCT v2 users, this removes the largest R dependency in the normalization path while preserving downstream compatibility with Seurat-derived workflows.

Reproducing#

pip install -e ".[test]"
export OMP_NUM_THREADS=8
python scripts/run_all_benchmarks.py
pytest tests/test_theta_od_parity.py -q
pytest tests/test_parity_r.py -m parity_r   # optional, needs R

Environment variable PYGLMGAMPOI_RSCRIPT overrides the Rscript path for parity tests (default: conda run -n st Rscript).