Usage
polishSpaNorm(
spe,
adj.method = c("auto", "logpac", "pearson", "medbio", "meanbio"),
scale.factor = 1,
psi.method = c("fixed", "profile"),
ls = c("fixed", "joint"),
cells = c("all", "fit"),
null = TRUE,
maxit = 50L,
tol = 1e-08,
engine = c("batch", "gene"),
batch.size = NULL,
backend = c("cpu", "auto", "gpu"),
BPPARAM = BiocParallel::SerialParam(),
overwrite = FALSE,
verbose = TRUE,
assay = NULL,
...
)
## S4 method for signature 'SpatialExperiment'
polishSpaNorm(
spe,
adj.method = c("auto", "logpac", "pearson", "medbio", "meanbio"),
scale.factor = 1,
psi.method = c("fixed", "profile"),
ls = c("fixed", "joint"),
cells = c("all", "fit"),
null = TRUE,
maxit = 50L,
tol = 1e-08,
engine = c("batch", "gene"),
batch.size = NULL,
backend = c("cpu", "auto", "gpu"),
BPPARAM = BiocParallel::SerialParam(),
overwrite = FALSE,
verbose = TRUE,
assay = NULL,
...
)
## S4 method for signature 'Seurat'
polishSpaNorm(
spe,
adj.method = c("auto", "logpac", "pearson", "medbio", "meanbio"),
scale.factor = 1,
psi.method = c("fixed", "profile"),
ls = c("fixed", "joint"),
cells = c("all", "fit"),
null = TRUE,
maxit = 50L,
tol = 1e-08,
engine = c("batch", "gene"),
batch.size = NULL,
backend = c("cpu", "auto", "gpu"),
BPPARAM = BiocParallel::SerialParam(),
overwrite = FALSE,
verbose = TRUE,
assay = NULL,
...
)
Details
SpaNorm’s mean for gene \(g\) in cell/spot \(i\) is
{gi} = {}g + a_1 w{i1} + {j } W_{ij} _{gj},
where \(w_{i1}\) is the log library size and its coefficient \(a_1\) is shared by every gene. The polish maximises each gene’s penalised negative binomial log-likelihood by damped Newton (see polishNB()), with \(a_1\) held at the fit’s value (ls = “fixed”), \(\bar{\mu}_g\) a per-gene unpenalised intercept, and the other coefficients ridge-penalised as SpaNorm() penalises them: by lambda.a for the biology and library-size splines, not at all for batch, times the number of cells/spots.
The polished coefficients are the optimum of the unwinsorised penalised likelihood. The normalisation (and SpaNormSVG()) still applies SpaNorm’s usual winsorisation, exactly as for an unpolished fit: it forms the mean through calculateMu(), which by default winsorises the upper tail of each gene’s log mean, and caps the dispersions at exp(median(log psi) + 4 MAD(log psi)).
A gene with no counts in the polished cells/spots has no finite optimum (its intercept runs to minus infinity), and nor has a gene with no counts in every polished cell/spot of one batch level (that level’s unpenalised coefficient runs to minus infinity), so neither is polished: it keeps its input coefficients and dispersion, and its diagnostics row reads polished = FALSE with iterations = 0 and held_out naming the reason, “all-zero” or “zero-batch-level” (NA for a gene that was polished). Genes held out for a batch level are counted in a warning. The levels are the groups of polished cells/spots that share one row of the model’s unpenalised columns, the intercept and the batch indicators: for a batch factor, its levels; for several batch variables, each combination of levels present, which can hold out a gene whose optimum under an additive batch design is finite. Groups are formed only when every unpenalised column other than the intercept is a 0/1 indicator; when one is not (a numeric batch covariate, or a spline left unpenalised by lambda.a = 0), only the all-zero rule applies. A gene the Newton engine cannot converge from either start also keeps its input fit (see polishNB()); a warning counts these genes, and how many had a singular information matrix. A batch matrix given with every level’s indicator makes the batch columns collinear with the intercept, and so every gene’s information singular: this is detected once, up front, by a QR of the unpenalised columns, and no gene is passed to the engine, so the outcome does not depend on rounding. The fit is marked polished either way.
With ls = “joint” the shared coefficient \(a_1\) is moved to the joint optimum of the total penalised log-likelihood over the polished genes: after the per-gene polish at the fit’s \(a_1\), Newton steps on the profile log-likelihood of \(a_1\) (each gene’s own coefficients profiled out; the information is the Fisher information of the polish) alternate with a warm re-polish of every gene at the candidate value, and a step is halved until the total does not fall. It stops when the standardised score \(|U|/\sqrt{I}\) (a1’s distance from its optimum in units of its own profiled SE) is below 1e-3, or after 10 steps; maxit and tol are the per-gene Newton’s, in the cold pass and every re-polish alike. With psi.method = “profile” each re-polish also re-profiles the dispersion. A gene whose information is singular at a step is left out of that step’s score, information and total. If the profiled information is not positive (\(a_1\) not identified: the log library size is collinear with the model’s other terms), the fit is the ls = “fixed” polish, with a warning. The joint estimate weights genes by their information, so bright genes dominate it, where SpaNorm()’s \(a_1\) is an unweighted mean over genes. Genes that are not polished (held out, or not convergeable) do not inform \(a_1\); they keep their input \(\bar{\mu}_g\), other coefficients and dispersion exactly, but take the new shared \(a_1\), since it is one value for every gene, so their fitted mean shifts by (a_1^{new} - a_1^{old}) w_{i1}. With no polished gene at all, \(a_1\) is left at the fit’s value with a warning.
The polished fit replaces ‘SpaNorm’ in the object’s metadata (@misc for Seurat) and the input fit is kept as ‘SpaNormUnpolished’. With null = TRUE a stored ‘SpaNormNull’ is polished with the same settings and the input kept as ‘SpaNormNullUnpolished’. SVG results in rowData were computed from the unpolished fit, so they are removed with a warning; rerun SpaNormSVG(). The normalised assay (‘logcounts’, or the ‘data’ layer for Seurat) is rewritten from the polished fit.
The polished fit’s polish slot (see isPolished()) holds settings (psi.method, ls, cells, maxit, tol, the penalty vector pen, the library-size coefficient before (a1.input) and after (a1), with psi.method = “profile” the dispersion’s search interval psi.range (as passed through …, else polishNB()’s default), and the SpaNorm version) and genes, one row per gene: polishNB()’s diagnostics plus loglik, the penalised log-likelihood over the polished cells at the returned fit (NA for a gene that was not polished), and held_out (see above). The fit’s loglik slot is left as the shared fit’s iteration trace. With ls = “joint”, settings also holds ls.iterations (accepted steps on \(a_1\)), ls.maxit and ls.tol (the cap on those steps and the standardised-score stop, 10 and 1e-3, the latter in units of a1’s profiled SE), ls.score (the final pooled score \(U\)), ls.se (\(1/\sqrt{I}\), the profiled standard error of \(a_1\)), ls.singular (genes left out of \(U\) and \(I\) at some step because their information was singular) and ls.converged, and the per-gene iterations, capped and singular include the warm re-polishes.
The negative binomial likelihood is defined on counts only, so a non-integer counts assay (for example a back-transform such as 2^logcounts - 1) is refused.