Parameters, performance and reproducibility

This article explains how to choose the settings of compareSpatial(), what determines its run time, and what makes its results reproducible. It assumes the concepts of How STcompare works. The examples use the simulated kidneys that ship with the package.

library(STcompare)
library(SpatialExperiment)
data(speKidney)
rast <- SEraster::rasterizeGeneExpression(speKidney, assay_name = "counts",
                                          resolution = 0.2, fun = "mean",
                                          square = FALSE)

The settings at a glance

Setting	Default	What it controls	When to change it
Rasterization resolution	set in SEraster	The pixel size, hence the number of shared pixels $N$, the run time and the spatial detail	Larger than the misalignment of the sections, smaller than the structures of interest
assay	the first assay	The values compared	Use a normalized, linear-scale assay (for example CPM) for the similarity
genes	the genes in both samples	Which genes are compared	To compare whole transcriptomes, choose the genes with a spatial pattern in both samples
surrogate	"remap"	How the surrogates take their values: the gene’s own values in the surrogate’s rank order ("remap"), or the nearly normal values the smoothing and the noise leave ("gaussian", the legacy functions)	"gaussian" to compare with the legacy functions and the published analyses
minDetected	0 with "remap"; $\sqrt{N}$ of the $N$ shared pixels with "gaussian"	Which genes are tested for correlation: those detected in at least this share of the pixels of both samples (0: every gene that is not constant)	With "gaussian", raise it (for example to 0.1) when testing many sparsely detected genes
delta	0.01, 0.05, 0.1, …, 0.9	The candidate smoothing bandwidths of the surrogates	When many genes have deltaGridEdge = TRUE
nPermutations	10000	The largest number of permutations per gene, hence the smallest p-value, $1 / (n + 1)$	Smaller for a quick look when many genes are significant; larger when few genes are significant among many
exceedances	10	When a gene stops early, hence the precision of large p-values	Rarely; Inf gives a fixed number of permutations
maxDistPrctile	0.25	The range of distances over which the variograms are matched	Rarely
foldChange	1	The similarity band, |log2(y / x)| <= foldChange (1 is two-fold)	log2(1.5) for 1.5-fold
minQuantile, minPixels	0.05, 0.1	Which pixels count for the similarity, and the least share of pixels kept	For genes expressed in a small part of the tissue
adjustMethod	"BH"	The adjustment of pValuePermute across the genes of the call, which gives pValuePermuteAdj	Any method of p.adjust()
seed	0	The random permutations	To check that a conclusion does not depend on the random draws
nThreads	getOption("STcompare.nThreads", 1)	The number of threads	Set it to the number of cores available

Preparing the input

Joint rasterization. Both samples must be rasterized in one call of SEraster::rasterizeGeneExpression(), after alignment, so that they share their pixels. compareSpatial() uses the pixels whose names are in both samples and stops if their coordinates differ.

Resolution. The resolution sets the number of shared pixels $N$. Larger pixels smooth over small misalignments and noise, but merge structures; smaller pixels resolve more detail, cost more run time per permutation and, when a pixel holds only one cell or spot, carry the noise of single measurements. As a guide, the case studies use pixels of about 10 Visium spots (kidney, 311 shared pixels) and of one Visium spot (brain, 2170 shared pixels).

Aggregation and normalization. Either aggregate the counts with fun = "sum" and normalize the pixels, for example to counts per million (CPM), as in the kidney case study, or normalize the cells or spots first and aggregate them with fun = "mean", as in the brain case study. The correlation can also be computed on log-transformed values. The similarity compares fold changes and needs a linear scale.

Genes. A gene that is constant in a sample (for example never detected) cannot be tested for correlation and gets status = "skipped". Every other gene is tested: the surrogates take the gene’s own values (surrogate = "remap", the default), so a gene detected in a handful of pixels gets a calibrated p-value, and such genes are cheap because their permutations stop early (see How STcompare works). On the whole AKI raster (32,285 genes on 311 pixels), the default tests the 17,902 genes that are detected in both sections and skips the 14,383 that are never detected in one of them. With surrogate = "gaussian", whose surrogates have nearly normal values, the correlation of a rarely detected gene is decided by the few pixels where both samples detect it and its p-value is far too small, so minDetected then defaults to $\sqrt{N}$ of the $N$ shared pixels (18 of 311, 47 of 2170), which skips 4,402 more genes of the AKI raster; genes above that filter but mostly zero can still get p-values that are somewhat too small, so with gaussian surrogates test the genes that have a spatial pattern in both samples, or raise minDetected, for example to 0.1. Genes without a spatial pattern in a sample have little to compare in any mode, which is why the case studies test spatially variable genes.

The smoothing grid: delta

Each surrogate is smoothed with every bandwidth $\delta$ of the grid, and the best match of the variogram, delta star, is kept. The default grid adds 0.01 and 0.05 to 0.1, …, 0.9: with the shorter grid, the smallest value was often the best one in the published analyses, a sign that the best bandwidth was smaller still.

compareSpatial() drops bandwidths whose window would hold fewer than 2 pixels ($\lfloor N\delta \rfloor < 2$), with a message. With 150 pixels, 0.01 is dropped:

pixels <- intersect(colnames(rast$A), colnames(rast$B))[1:150]
small <- compareSpatial(rast$A[, pixels], rast$B[, pixels],
                        nPermutations = 199, nThreads = 2)
#> compareSpatial: delta 0.01 dropped: with 150 shared pixels the smoothing window must hold at least 2 pixels (delta >= 0.0133)
#> compareSpatial: 1 gene on 150 shared pixels in 0.1 s (199 permutations in total, 2 threads)
attr(small, "params")$deltaGrid
#>  [1] 0.05 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90

deltaGridEdge = TRUE flags genes for which more than half of the surrogates, in either direction, chose the smallest or the largest bandwidth, and a message counts them. Their p-values remain valid, but the surrogates may reproduce their autocorrelation less well. If many genes are flagged, add values beyond the edge, for example 0.005 when the smallest value is chosen (it needs $N \ge 400$). Each bandwidth adds to the run time, the small ones on many pixels the most.

Permutations: nPermutations and exceedances

A gene stops when either direction has exceedances null correlations at least as extreme as the observed one, or at nPermutations.

nPermutations sets the resolution of small p-values. A gene that reaches the limit has $p = (b + 1) / (n + 1)$, at least $1 / (n + 1)$. With the Benjamini-Hochberg adjustment of $G$ genes, the gene of rank $k$ is significant at level $\alpha$ only if $p \le \alpha k / G$. If you expect at least $k$ significant genes, a limit of

$$n \ge \frac{G}{\alpha k} - 1$$

lets genes without any exceedance reach significance:

G <- c(100, 1000, 10000)
k <- c(1, 10, 100)
needed <- outer(G, k, function(G, k) pmax(0, ceiling(G / (0.05 * k) - 1)))
dimnames(needed) <- list(genes = G, "significant genes expected" = k)
needed
#>        significant genes expected
#> genes        1    10  100
#>   100     1999   199   19
#>   1000   19999  1999  199
#>   10000 199999 19999 1999

The default, 10,000, suffices for 1000 genes with at least 2 significant genes, or for 10,000 genes with at least 20. When many genes are significant, as in both case studies, 1000 gives the same significant genes 6 to 7 times faster. If genes reach the limit without any exceedance and are still not significant after the adjustment, compareSpatial() says so in a message: a larger nPermutations would resolve them.

exceedances sets the precision of p-values that are not small. The relative standard error of a sequential p-value is about $1 / \sqrt{h}$ for $h$ exceedances, about 30% for the default 10. A gene whose patterns are unrelated costs on average about $h (1 + \ln(n / h))$ permutations, around 80 for the defaults, instead of $n$. A larger exceedances gives more precise large p-values at a proportional cost; exceedances = Inf gives every gene exactly nPermutations permutations and $p = \max(p_X, p_Y)$, that is pValuePermute = max(pValuePermuteX, pValuePermuteY).

The variogram range: maxDistPrctile

The variograms that the surrogates must match are computed up to the maxDistPrctile quantile (25% by default) of the distances between the shared pixels. Pairs of pixels far apart are fewer, their variogram is noisier, and it reflects the shape of the tissue more than the spatial autocorrelation of the gene. A smaller value matches only the short-range autocorrelation; a larger one also the longer-range structure. The default suits most data; compareSpatial() stops with an error that suggests a larger value when too few pixel pairs fall in the range to compute a variogram.

Similarity: foldChange, minQuantile and minPixels

foldChange is on the log2 scale: 1 means within two-fold, log2(1.5) within 1.5-fold. Each gene gets a threshold in each sample, its minQuantile quantile over the shared pixels (the columns t1 and t2), and a pixel counts when it is above the threshold in either sample. If fewer than a share minPixels of the pixels count, the similarity is NA. For a gene expressed in a small region, lower minPixels; to ignore background levels, raise minQuantile.

Threads and run time

The run time is proportional to the total number of permutations (the sum of the nPermutations column) times the cost of one permutation in both directions, divided by the number of threads. All genes share the parts that depend only on the pixel coordinates, and the permutations of all genes are spread over nThreads threads. The cost of a permutation grows with the number of pixels and of bandwidths: one permutation in one direction took about 0.17 ms of one thread for the kidney data (311 pixels, 11 bandwidths) and 1.25 ms for the brain data (2170 pixels, of which the variograms use a random subset of 1000); the remapping of the surrogate onto the gene’s values, one sort of the pixel values, is about 7% of that.

Measured run times of the case studies (Apple M1 Ultra, default settings apart from nPermutations):

Data	Genes	Shared pixels	nPermutations	Threads	Run time
Kidney, control vs AKI	1046	311	1000	8	32 s
Kidney, control vs AKI	1046	311	1000	16	18 s
Kidney, control vs AKI	1046	311	10000	16	118 s
Brain, MERFISH vs Visium	325	2170	1000	8	54 s
Brain, MERFISH vs Visium	325	2170	1000	16	29 s
Brain, MERFISH vs Visium	325	2170	10000	16	189 s

The genes that reach the limit dominate: with the default limit, they used 82-84% of all permutations. To predict the run time of a large analysis, run a random subset of genes with the settings you plan to use and scale the run time (attr(res, "runtime")) by the number of genes.

Threads do not change the results. Each permutation’s random numbers depend only on the seed, the gene and the permutation number, so the results are identical for any number of threads:

one <- compareSpatial(rast[c("A", "C")], nPermutations = 999, nThreads = 1,
                      verbose = FALSE)
two <- compareSpatial(rast[c("A", "C")], nPermutations = 999, nThreads = 2,
                      verbose = FALSE)
identical(as.data.frame(one), as.data.frame(two))
#> [1] TRUE

Set the default for a session with options(STcompare.nThreads = 8). The threads are threads of the compiled code: no R processes are started, so they work the same on every operating system, and they never call R. Use at most the number of cores available, and fewer on a shared computer.

Memory. keepNulls = TRUE keeps two numbers per permutation and direction (the null correlation and the delta star), up to 320 MB for 1000 genes that all run to 10,000 permutations. Use it for a few genes.

Progress. In an interactive session, a progress line shows the share of the work done, the genes finished, the permutations so far, and an estimate of the remaining time:

compareSpatial:  61% | genes done 640/1046 | 1.9M permutations | elapsed 0:12 | ETA 0:08

In a console narrower than the line, the elapsed time is left out. At the end, the line shows the total number of permutations, the sum of the nPermutations column. progress = FALSE turns it off and verbose = FALSE the other messages.

Reproducibility

Results depend only on the data, the settings and seed. Each gene and direction has its own random stream, determined by the seed, the gene name, the direction and the permutation number. A gene therefore gets the same permutations whatever the number of threads, the order of the genes, or the other genes in the call; only pValuePermuteAdj, which adjusts across the genes of the call, depends on the other genes. The permutations are drawn over the shared pixels in the column order of the first sample, as in the legacy functions: reordering its pixels, or swapping the two samples, gives other permutations. A different seed does too: the null correlations are other draws, and the p-values agree within Monte Carlo error (here two unrelated patterns, with keepNulls = TRUE to see the first null correlations of each run):

data(simRanPatternRasts)
unrelated <- function(seed) {
  compareSpatial(simRanPatternRasts[[1]], simRanPatternRasts[[2]], seed = seed,
                 tests = "correlation", nThreads = 2, verbose = FALSE,
                 keepNulls = TRUE)
}
runs <- list(seed0 = unrelated(0), seed1 = unrelated(1))
t(sapply(runs, function(res) {
  c(pValuePermute = res$pValuePermute, nPermutations = res$nPermutations,
    firstNulls = round(attr(res, "details")[[1]]$nullX[1:3], 3))
}))
#>       pValuePermute nPermutations firstNulls1 firstNulls2 firstNulls3
#> seed0     0.9090909            11      -0.177      -0.288      -0.054
#> seed1     0.9090909            11      -0.051       0.038       0.131

The session’s random number generator is left unchanged. compareSpatial() draws its permutations and noise in compiled code. With more than 1000 shared pixels, it also draws the pixels of the variogram subsample from seed with R’s default generator, in a local state that it restores afterwards. It does not affect the random numbers of the rest of an analysis:

set.seed(42)
before <- runif(1)
set.seed(42)
invisible(compareSpatial(rast[c("A", "B")], nPermutations = 99,
                         nThreads = 2, verbose = FALSE))
identical(before, runif(1))
#> [1] TRUE

Platforms. The results are deterministic on a given platform. Between platforms (for example macOS and Linux), the Gaussian noise of the surrogates can differ in the last bits, and on regular grids of pixels, pixel pairs that lie exactly at the largest variogram distance can be binned differently. The surrogates then differ slightly, and p-values can differ within Monte Carlo error.

Record the settings. The result keeps every setting, including the bandwidths actually used, the seed and the assays, in attr(res, "params"), along with the call and the run time.

The legacy functions

The functions of earlier versions remain available with their interface and semantics. They are computed by the same compiled engine as compareSpatial(), much faster than before, and reproduce the published analyses: the same delta stars, the same null correlations up to rounding, and the same p-values, except that p-values are now $(b + 1) / (B + 1)$, never 0, and that the gene-wise functions return them raw in pValuePermuteX and pValuePermuteY and adjusted across genes in pValuePermuteXAdj and pValuePermuteYAdj (earlier versions had no Adj columns: the iterative function returned the adjusted values under the first two names, and spatialCorrelationGeneExp(), which adjusted one gene at a time, returned unadjusted values under them).

Function	What it does	Use it to
compareSpatial()	Both tests for every gene; adaptive permutations; one p-value per gene	Analyse new data
spatialCorrelationGeneExp()	The correlation test for every gene with a fixed number of permutations; p-values of both directions	Reproduce published results
spatialCorrelationGeneExpIterPermutations()	The same, with a second round of more permutations for genes that look significant	Reproduce published results
spatialCorrelationGeneExpWithinSample()	The correlation test for pairs of genes within one sample	Compare genes with each other
spatialCorrelation(), viladomatCorrelation()	The correlation test for vectors of values and coordinates	Use the test outside of SpatialExperiment objects
spatialSimilarity()	The similarity, with the class of every pixel	Get the classes as a table; linearRegression() and pixelClass() plot them
plotGene(), plotOverview()	The standard figures of a compareSpatial() result: one gene, and all genes	Inspect and publish results

Re-running the six published analyses with the legacy functions on macOS (arm64), where they were computed, gives identical delta stars and null correlations within $10^{-13}$, and takes seconds instead of hours (16 threads of an Apple M1 Ultra, against the authors’ runs on 20 to 22 workers):

Analysis	Genes	Shared pixels	Run time	Authors’ run time
Kidney, control vs AKI (100 permutations)	1046	311	1.9 s	not recorded
Kidney, control vs AKI (100, then 1000 permutations)	1046	311	14 s	not recorded
MERFISH replicates, affine alignment	483	1299	40 s	7.65 h
MERFISH replicates, STalign alignment	483	1371	40 s	16.8 h
Brain, MERFISH vs Visium	325	2170	13 s	1.78 h
Brain cell types	16	2174	0.9 s	6.95 min

The legacy functions share one sequence of permutations across genes and directions (the published results depend on it), use a fixed number of permutations per round, and report a p-value per direction. compareSpatial() uses independent permutations for each gene and direction, stops early for genes that are clearly not significant, gives its surrogates the gene’s own values, and reports one p-value per gene. On the published analyses, with its default settings (no gene is skipped; nPermutations = 1000 gives the same significant genes), the two agree on the significance of 98% of the kidney genes and 86% of the brain genes, counting a gene as significant with the legacy function when both of its adjusted p-values (pValuePermuteXAdj and pValuePermuteYAdj) are below 0.05:

Analysis	Significant with both	Legacy function only	compareSpatial() only
Kidney, 1046 genes	729	2	23
Brain, 325 genes	128	0	45

The kidney genes classified differently are near the significance threshold (pValuePermuteAdj between 0.02 and 0.08 with compareSpatial(), and a larger of the two adjusted legacy p-values between 0.04 and 0.14), where Monte Carlo error and the different adjustments (one p-value per gene rather than one per direction) decide. The 45 brain genes that only compareSpatial() finds are all positively correlated ($r$ = 0.06 to 0.24), and they come from the remapped surrogates: with surrogate = "gaussian", the surrogates of the legacy function, and its $\sqrt{N}$ filter (which skips 28 of the brain genes), compareSpatial() agrees with the legacy function on 95% of the brain genes (122 significant with both, 6 with the legacy function only, 10 with compareSpatial() only). In the calibration study, gaussian surrogates were conservative for skewed, zero-inflated genes with a spatial pattern, as the log-normalized brain pixels are, so the additional genes are consistent with power that the gaussian surrogates lose (bench/calibration-results.md in the source repository).