R vs Python

How SpaNorm (R) and spanormpy (Python) correspond, and where they intentionally differ

spanormpy is a faithful port of SpaNorm 1.7.15. The two packages implement the same algorithms, and a large test suite compares spanormpy with reference values computed by the R package. This page lists every exported R function with its Python counterpart, and the places where spanormpy deliberately or unavoidably differs.

Names are snake_case in Python with R’s defaults (sample.p is sample_p), and R’s BPPARAM is n_jobs. The main differences in how the packages are used:

R (SpaNorm) Python (spanormpy)
Container SpatialExperiment (also Seurat) AnnData (also SpatialData)
Normalised values logcounts assay layers["logcounts"]
Model fit metadata(spe)$SpaNorm var, varm, obsm, uns["spanorm"]
Updates the object returns a modified copy in place (copy=True to return a copy)
Plots ggplot2 matplotlib
Parallelism BPPARAM (BiocParallel) n_jobs

R to Python mapping

Every export of R SpaNorm 1.7.16 (commit f233b7b, NAMESPACE; the export list is unchanged from 1.7.15) and its spanormpy counterpart. “internal” rows are building blocks that R exports for the GPU and polishing engines; they live in private modules of spanormpy (no stability guarantee), or are plain torch operations where R needed a wrapper around torch for CPU/GPU dispatch.

Argument names are snake_case with R’s defaults (sample.p is sample_p); R’s BPPARAM is n_jobs. Intentional deviations are listed in r-discrepancies and summarised in the README (“Differences from R”).

R export Python Notes
SpaNorm spanormpy.spanorm in place on AnnData (copy=True to return a copy); output in layers["logcounts"]
SpaNormPCA spanormpy.spanorm_pca result in obsm["X_spanorm_pca"]
SpaNormSVG spanormpy.spanorm_svg statistics in var["svg_F"], var["svg_p"], var["svg_fdr"]
add_vec_mat_gpu internal: torch.add (broadcasting) R wrapper not needed; torch tensors are used directly
calculateMu spanormpy.calculate_mu
checkGPU internal: spanormpy._backend.check_gpu
diag_mat internal: torch.diag_embed R wrapper not needed
dnbinom_gpu internal: spanormpy._nbdist.dnbinom
fastSizeFactors spanormpy.fast_size_factors writes obs["spanorm_size_factors"]
filterGenes spanormpy.filter_genes deviation: the default inplace=True subsets the genes of adata and returns None; inplace=False returns R’s logical vector (see r-discrepancies.md)
fitNB spanormpy.fit_nb
getBackendDevice internal: spanormpy._backend.get_backend_device
getBackendDtype internal: spanormpy._backend.get_backend_dtype
getGPUMemoryBudget internal: spanormpy._backend.get_gpu_memory_budget
gpuDtypeBytes internal: spanormpy._backend.gpu_dtype_bytes
hasBadValues internal: spanormpy._backend.has_bad_values
invert_mat internal: torch.linalg.inv R wrapper not needed
invert_mat_batched internal: torch.linalg.inv (batched) R wrapper not needed
isPolished spanormpy.is_polished takes a stored SpaNormFit
is_torch_tensor internal: torch.is_tensor R wrapper not needed
matmul_gpu internal: torch.matmul R wrapper not needed
mult_vec_mat_gpu internal: torch.mul (broadcasting) R wrapper not needed
nbAbsorbGramBatch internal: spanormpy._solver.nb_absorb_gram_batch
nbBlockDesign internal: spanormpy._block_design.nb_block_design S3 methods (as.matrix, dim, dimnames, print) are methods of NBBlockDesign
nbDevianceMoments internal: spanormpy._quasilik.nb_deviance_moments
nbGramBatch internal: spanormpy._solver.nb_gram_batch
nbMuFloor spanormpy.nb_mu_floor
nbNewtonSolver internal: spanormpy._solver.nb_newton_solver
nbNewtonSolverBatch internal: spanormpy._solver.nb_newton_solver_batch
nbProfilePsi spanormpy.nb_profile_psi
nbUnitDeviance internal: spanormpy._quasilik.nb_unit_deviance
plotCovariate spanormpy.pl.plot_covariate matplotlib; returns Axes
plotSpatial spanormpy.pl.plot_spatial matplotlib; returns Axes
polishNB spanormpy.polish_nb
polishSpaNorm spanormpy.polish_spanorm
qlDispersion internal: spanormpy._quasilik.ql_dispersion
qlMomentTable internal: spanormpy._quasilik.ql_moment_table
resetGPUCache internal: spanormpy._backend.reset_gpu_cache
rowSums_gpu internal: torch.sum (dim=1) R wrapper not needed
setGPUMemoryBudget internal: spanormpy._backend.set_gpu_memory_budget also spanormpy.settings for the polish memory budget
tcrossprod_gpu internal: torch.matmul (A @ B.T) R wrapper not needed
toGPUMatrix internal: spanormpy._backend.as_tensor
toGPUVector internal: spanormpy._backend.as_tensor
toRMatrix internal: torch.Tensor.numpy R wrapper not needed
topSVGs spanormpy.top_svgs returns a DataFrame
tpsBasis internal: spanormpy._tps.tps_basis
SpaNormFit internal: spanormpy._fit.SpaNormFit exported class; stored in AnnData var/varm/obsm/uns["spanorm"], the $ method is attribute access

Discrepancies from R SpaNorm

geneBlockCount integer overflow

geneBlockCount() (R/gpuFunctions.R) does as.integer(ceiling(needed.bytes / budget.bytes)) before clamping to n_genes. When the ratio exceeds 2^31 (for example a 1-byte budget on a large problem) as.integer yields NA with a warning and the function returns NA. spanormpy._backend.gene_block_count clamps in floating point and returns n_genes in that case. The two agree everywhere R returns a value.

pnbinom / qnbinom on accelerators

On CPU tensors pnbinom/qnbinom are exact (scipy incomplete beta; for size >= 1000 and q <= 5000 the CDF is summed from the pmf in extended precision because scipy’s betainc loses ~5e-11 relative accuracy at size 1e6), and qnbinom uses R’s 64*eps fuzz. On accelerators they use R’s GPU algorithm (cumulative pmf sum), whose CDF is accurate to ~1e-11 relative in far tails, so qnbinom there uses a 1e-10 relative fuzz instead. A p lying exactly on a CDF step can therefore differ from R by one count on accelerators, never on CPU. If the search cap (1e6) is reached, unfinished elements are NaN with a warning.

Nested absorption with a group holding no cells (.newtonSolver)

In the diagonal nested path .newtonSolver() sums weights and X*w per group with rowsum() on a factor, which drops a group with no cells. cvec and B then come back short, pen_z recycles against them and the Schur complement is misaligned (or errors). The batched .absorbBatch() already keeps the slot through .segmentSum(); R’s own test notes the case is unreachable from the package’s designs. spanormpy._solver._newton_solver uses the same zero-filled segment sum for both, so an empty group keeps its (penalty-only) slot. Results agree with R wherever R returns a valid answer.

Singularity test in the dense solvers

R’s solve() signals a singular system through LAPACK dgesv (exact zero pivot) or a reciprocal condition number below .Machine$double.eps (estimated by dgecon). _solver._lu_checked applies the same two tests with the exact 1-norm condition number, so a matrix whose estimated and exact rcond straddle machine epsilon could in principle be classified differently.

Batched Cholesky on a failed slice

R’s torch branch of .newtonSolverBatch()$xcov calls torch_cholesky_inverse() on the whole stack, which raises on a singular slice; its base-R branch returns NA for that slice. Python follows the base-R behaviour: failed slices are NaN and their neighbours are unaffected.

edgeR compute_unit_nb_deviance integer division (matched, not fixed)

edgeR’s C routine evaluates the small-dispersion (phi < 1e-4) Poisson-like branch with 2/3*resid, which is C integer division and equals 0. spanormpy._edger._unit_nb_deviance reproduces this exactly (0.0 * resid), so deviances and the residual variances fed to squeezeVar agree with R. It only matters for trended dispersions below 1e-4.

edgeR one-group Newton fit: unassigned coefficient on non-convergence (matched)

In glm_one_group_vec the output *beta is only assigned when the Newton iteration converges within maxit (50); fit_one_group_mat / average_log_cpm then copy the uninitialised local ocoef into the result. R observes ~4.6e-310 (effectively 0) for the first such row and the previous row’s value afterwards (same stack slot). spanormpy._edger._glm_one_group reproduces that (0.0 for a leading row, otherwise the previous Newton-fitted row’s coefficient). This is undefined behaviour in edgeR, matched empirically on the R-vs-Python fuzz (all non-converging scenarios reproduce); it affects very sparse, strongly overdispersed genes at large grid dispersions.

cleanSizeFactors: nls stopping point

scuttle::cleanSizeFactors (the positive = TRUE repair of non-positive pooled estimates) fits a nonlinear trend with nls, which stops at a 1e-5 relative-offset criterion. The Python Gauss-Newton uses LAPACK QR where R uses LINPACK dqrdc2, so the replaced entries agree to ~1e-9 relative rather than 1e-10; all other size factors agree to 1e-10.

fitNB / fitSpaNormNB: tol reaches only the outer loop

fitSpaNormNB(..., tol = ) has tol as a formal argument, so it is not part of the ... forwarded to fitNBGivenPsi; the inner IRLS loop therefore always uses its default tol = 1e-4 and fitNB(tol = ) only changes the outer dispersion loop. Matched in fit_spanorm_nb (pinned by test_tol_reaches_only_outer_loop).

fitNBGivenPsi / fitSpaNormNB: NaN in step-halving and convergence comparisons

R evaluates if (is.na(degener) | degener) (NaN log-likelihood sums are treated as degenerate, matched), but step[check.gene] <- ... with an NA in the per-cell logical index, and ifelse((l1 - l2) / abs(l2) < tol, TRUE, FALSE) with NaN inputs, error or give NA in R. Python treats NaN comparisons as False: such cells keep their step size, and a NaN convergence ratio counts as not converged.

normalise*: a single NaN dispersion (fixed in SpaNorm 1.7.16; spanormpy matches)

Up to SpaNorm 1.7.15 a NaN in psi made psi.max = exp(median(log psi) + 4 * mad(log psi)) NaN and so every gene’s logPAC, Pearson and medbio output NaN; a psi of 0 (log = -Inf) dragged the median to -Inf. SpaNorm 1.7.16 computes the cap on the finite log-dispersions only (Inf, i.e. no winsorisation, if none are finite), so only the NaN gene is NaN and psi = 0 genes are excluded from the median/MAD. compute_psi_max (used by winsorise_psi and global_psi_max) is the same rule; the source-of-truth R version for this behaviour is now 1.7.16 (pinned by test_degenerate_psi_matches_r, with references regenerated from 1.7.16).

Batch factor levels: locale collation

R’s factor() sorts character batch labels with the session locale (en_US.UTF-8: a < A < b < B, punctuation < digits < letters); pandas and numpy sort by codepoint (A < B < a < b). spanormpy._io.check_batch sorts strings with an approximation of the glibc en_US.UTF-8 collation (punctuation < digits < letters, case-insensitive first, lower case before upper case on ties). Exotic characters (accents, non-Latin scripts) may still order differently from R. Which level is the reference only changes the parametrisation of the batch columns, not the fitted values or the normalised data. For full control pass a pandas.Categorical, whose category order is used as given. Numeric vectors are one covariate column and booleans are factors, as in model.matrix. A single-level batch raises R’s “contrasts can be applied only to factors with 2 or more levels”; unused categorical levels give all-zero columns, as in R.

Adjustment backend

R’s SpaNorm() calls the normalise* functions without a backend (they use the accelerator whenever one is available). spanorm(backend=...) passes the user’s backend to the adjustment too, so backend="cpu" really keeps everything on the CPU.

Output sparsity

R always stores logcounts as a sparse matrix. spanorm writes a dense float64 layer for dense input and a CSR layer for sparse input; values are identical.

Default pipeline (own size factors and sampling)

Python’s quick_cluster and RNG differ from R’s, so the default pipeline is checked for closeness only (see tests/test_spanorm.py::test_default_pipeline_close).

Rerun after a gene or spot subset

R keeps metadata(spe)$SpaNorm when an SPE is subset, and .spaNormCore refits because fit$ngenes != nrow(emat) (genes) or fit$ncells != ncol(emat) (spots). spanorm on AnnData differs for genes only: per-gene fit data live in var/varm, so a gene subset carries its slice of the fit and is reused (with the whole fit’s psi_max, as R’s subsetFitGenes), which equals R’s adjustment of the full fit on those genes. R’s refit behaviour is reachable with overwrite=True (pinned against R in tests/test_api_robustness.py). A spot subset refits, as R does (detected through the ncells stored in uns["spanorm"]["spanorm"]).

Degenerate inputs

An all-zero gene and a zero-total spot with supplied size factors are handled by R (finite logcounts, size factor floored at 1e-8) and identically here. Without supplied size factors R’s calculateSumFactors fails on a zero-total spot (“size factors should be positive”); so does spanorm. A batch factor with an unused level gives an all-zero unregularised column and a singular system in both (R: dgesv ... exactly singular; here torch’s LinAlgError, a different exception type).

qnbinom in the far upper tail (CPU)

R >= 4.2 fuzzes qnbinom much less than the older 64*eps when the tail probability is tiny: at p = 1 - 1e-10 the 64*eps fuzz (1.4e-14) is ~1e-4 of the tail and moves the quantile by 2-3 counts. qnbinom on CPU now searches the survival function with sf(k) <= (1-p) + 8*eps*p when 1 - p < 1e-6, and keeps the 64*eps CDF fuzz elsewhere (the CDF R feeds back through p carries errors up to ~1e-14 that only the larger fuzz absorbs). This is needed for nb_deviance_moments to reproduce R’s summation window to 1e-10. At 1 - p ~ 1e-13 R’s own CDF is quantised at 1.1e-16 and its quantile can differ from the exact one by a few counts; not reproducible.

filter_genes returns None by default

R’s filterGenes() returns a logical vector and leaves the object alone. filter_genes follows the scanpy convention: the default inplace=True removes the genes failing the filter from adata and returns None. filter_genes(adata, inplace=False) leaves adata unchanged and returns R’s boolean mask (genes with rowMeans(counts > 0) >= prop).

polish_spanorm after a gene or spot subset

polishSpaNorm calls getSpaNormFit, which in R errors when fit$ngenes/fit$ncells differ from the object. polish_spanorm follows R for spot subsets (detected through the ncells stored in uns["spanorm"][name], _fit.fit_is_stale_for_spots): it raises “… model dimensions do not match the data. Please rerun ‘spanorm’.” (or 'spanorm_svg' for the null fit), also for a gene plus spot subset, and writes nothing. A gene subset reuses the stored fits (see “Rerun after a gene or spot subset”): the genes are polished independently, so the result equals the full run’s genes. The polished subset fit keeps the whole fit’s psi_max (threshold from the unpolished dispersions), because the other genes’ polished dispersions are not available. A spot reorder keeps the count and reuses the fit. get_spanorm_fit(validate=True) raises the same error.

polish_spanorm bookkeeping

adata.uns["spanorm"]["adj_method"] is updated to the method used to rewrite logcounts. SVG columns removed on polishing are svg_F, svg_p, svg_fdr, svg_F_raw (R: svg.F, svg.p, svg.fdr). The polish settings key SpaNorm (R package version) is spanorm (spanormpy version). Per-gene held_out is stored as a categorical var column (h5ad cannot hold missing strings in an object column). R’s Seurat/DelayedArray branches and the pre-polish-slot RDS fixture have no Python counterpart.

spanorm_svg bookkeeping and subsets

SpaNormSVG writes svg.F, svg.p, svg.fdr and returns the unclamped statistic only as an attribute of svgTest()’s table, which rowData<- drops; spanorm_svg keeps it as var["svg_F_raw"] (svg_test(...).attrs["F_raw"] for the table). The “SVG results exist” warning fires when all three of svg_F, svg_p, svg_fdr are present, and all four columns are replaced. Like polish_spanorm, a spot subset of an object carrying a fit is an error (“… model dimensions do not match …”), and a gene subset reuses the stored full and null fits (the null of a gene subset keeps the whole fit’s dispersion-winsorisation threshold), so the subset’s statistics equal the full run’s rows. A null that has to be fitted on a gene subset winsorises with the subset’s own dispersions. spanorm_svg gained a random_state argument (seed of the null’s dispersion subsampling, R: the global RNG). The paired null reuses polish_spanorm_fit, with the full fit’s recorded psi_method, ls, cells, maxit, tol and psi_range, on the CPU backend (R does not pass backend to the polish either). top_svgs follows R for NA FDRs and for the counts check; see “top_svgs and NaN FDRs” below.

top_svgs and NaN FDRs

R’s topSVGs subsets with results[results$svg.fdr <= fdr, ]; a NA FDR makes that an all-NA row (named NA, NA.1, …), which order() puts last and min(n, nrow) counts. This is an indexing artefact, but top_svgs reproduces the row count and order: genes with a NaN FDR are returned after the sorted finite rows, in gene order, counted toward n, with all statistics NaN. Only the names differ: they keep their real gene names instead of R’s NA junk names. top_svgs also takes layer and runs the counts check R’s checkSPE does (it supersedes the earlier note that counts are not checked).

SpaNormPCA with NaN-FDR genes

topSVGs returns all-NA rows named NA, NA.1, … for genes whose FDR is NaN (see above), and R’s getResiduals(...)[svgs, ] then fails with “subscript out of bounds”. spanorm_pca raises a ValueError in the same situation (use a smaller nsvgs or set the FDRs). percentVar/variance cover only the ncomponents computed components (as in R), so percent_var sums to 100 over them, not over all components.

spanorm_pca(n_jobs=...) stands in for R’s BPPARAM but has no effect (the decomposition uses BLAS threads). random_state seeds the optional randomised SVD (R’s BSPARAM algorithms have their own RNG use).

Plotting

plot_spatial(crop=True) only acts together with img=True (R: crop & img); R drops image pixels outside the spot extent, here the axis limits are set to the extent padded by the spot radius so edge spots are not clipped. With circles=True the default radius is R’s fixed 100 (times the image scale factor). show=None means False, and the Axes is always returned.