Compare the SCTransform model against Seurat, per gene

[shanikawm]

Closes the highest-risk gap from the fresh-install audit. Tutorials 1–4 were compared by eye — figures only, no numeric handoff — and SCTransform is the one where that has actually cost something: this model was once wrong in four separate ways at once (theta anti-correlated with R's at −0.89, residual variance ranking genes near randomly) and the tutorial still drew a perfectly plausible UMAP. A picture cannot fail.

What it compares

pbmc3k_sctransform_verify.R now dumps Seurat's own SCTModel.list[[1]]@feature.attributes, the ranked variable features, and the clips. The Python side writes the same table, and --report puts them side by side over all 12,572 modelled genes.

To make that possible, sctransform stores the whole fitted model on the SCT assay's meta_data under Seurat's column names — residual_mean, detection_rate, (Intercept) and log_umi join the existing residual_variance, theta and gmean. The regularized intercept and slope, which are the model, previously could not be inspected at all.

Results

quantity	agreement
detection_rate · gmean	max abs diff 5.6e-16 · 1.2e-12
(Intercept) · theta	Spearman 1.0000
non-overdispersed genes	3,848 on both — the same set (Jaccard 1.0000)
residual_variance	Spearman 0.9986, Pearson 0.9996
variable features	2,913/3,000; 98 of the top 100
clips, cell/gene counts, cluster counts	all identical

residual_mean is the one column that does not track by rank (Spearman 0.71, Pearson 0.99) — and the one nothing downstream reads, since Seurat records it but selects on residual_variance. Reported rather than hidden; the disagreement sits in genes whose residual mean is ~1e-3 or smaller.

Two things the comparison surfaced

The vignette's ±1 cluster gap is gone. It described shanuz resolving 13 clusters against Seurat's 12, blamed on the RNG and the differing clustering libraries. Both arms now agree exactly — 12 under SCTransform, 11 under LogNormalize. PR #55's graph work closed it, not anything in SCTransform.

The module docstring was overstating fidelity. It claimed 99.7% variable-feature agreement and theta at 0.96. The reproducible comparison measures 97.1% and Spearman 1.0000. The docstring now carries numbers the tutorial prints instead of hand-measured ones. (The vignette's table was already correct at 2,913/97.1% — only the docstring had drifted.)

Verification

• 803 passed, 25 skipped (up 4)
• Four new tests, all mutation-tested. One was decorative on the first attempt: the fixture's largest residual was 6.2 against a clip of 17.3, so clipped and unclipped agreed and the clip assertion proved nothing. It now uses a gene carrying a spike the depth model cannot absorb, and asserts the clip is actually reached before relying on it.
• Full chain re-run from deleted artifacts and reproduces exactly.
• ruff clean on shanuz/sctransform.py; the tutorial's 10 E402 are the sys.path idiom every tutorial uses and are identical on main.

🤖 Generated with Claude Code