henhaobio
Sign in

bioRxivBioinformaticsPreprint

Statistical tests for bivariate spatial association across multi-omics data with disjoint coordinates

Stijn HawinkeliD, Wangjun HuiD, Britta VelteniD, Steven MaereiD

PDFDOICode
In-depth full-text digest available · Jump to digest ↓

Abstract

Spatial biology has entered a new era of multimodal profiling, with multiple, high-dimensional spatial omics types being measured on consecutive tissue slices, or co-assayed on the same slice. Interest then lies in statistical testing for spatial association between the features of the different modalities, to gain insight in biological processes. One major challenge is the multitude of bivariate combinations, leading to high computational demands. Another difficulty is the difference in spatial resolution between technologies, implying no one-to-one matching between the measurement spots of the different modalities, even after alignment. As a result, common statistical measures such as joint distributions and correlations are not defined, and tests need to rely on spatial vicinity only. Moreover, we argue that many existing bivariate association tests address an inappropriate null hypothesis, or make inappropriate assumptions, both implying absence of spatial autocorrelation in any of the features and leading to misleading conclusions. As a remedy, we modify tests for the detection of spatially variable genes (Morans I, Gaussian processes and generalized additive models) to derive bivariate spatial association tests across modalities with non-overlapping coordinate sets, and provide variance estimators that do account for spatial autocorrelation. We develop inference methods for single sections as well as for replicated experiments with multiple sections, and compare their performance in nonparametric and parametric simulations. Finally, we apply the newly developed methods to two co-assayed spatial transcriptomics and metabolomics datasets from mouse and human. The full suite of tests is available from github.com/sthawinke/sbivar as the R-package sbivar. Author summarySpatial biology investigates molecules and cells in their spatial context in living tissues. For this purpose, recent technologies measure different classes of biomolecules on the same or adjacent tissue sections. Here we focus on methods to detect colocalisation of molecule pairs, which may provide clues on their biological function, e.g. involvement in the same pathways. Many existing tests for identifying such colocalized pairs rely on the two technologies being measured on the exact same locations, which is often not the case due to differences in resolution. Hence they first require location matching by granulating the technology with the best resolution to that of the worst, discarding information and failing to exploit the full spatial resolution of the data. Moreover, existing methods ignore spatial patterning or autocorrelation naturally present in spatial omics data, leading to many false findings. As a remedy, we present a suite of scalable methods that work directly on disjoint coordinate sets and do account for spatial autocorrelation, and confirm their validity in simulation studies. Next we apply the new methods to metabolite and gene expression measurements of mouse and human datasets, uncovering interesting colocalized gene-metabolite pairs. The methods presented are available in the sbivar R-package from github.com/sthawinke/sbivar.

In-depth summary

Claude Opus 5.5 · xhigh effort · 8 Nov 2026 · from the full text of version 2 (23 Jul 2026) · CC BY 4.0

Autocorrelation-aware bivariate Moran's I, bivariate Gaussian process tests and generalized additive models (R package sbivar) test cross-modality colocalisation on disjoint coordinates; many existing permutation tests inflate false positives.

Cross-modality colocalisation tests in spatial omics usually coarsen the finer-resolution modality to match spots, and most ignore spatial autocorrelation (SAC), so they effectively test a randomisation null and report univariate patterns as associations. SAC-aware options the paper reviews need joint coordinates (modified t-test, SpaceBF) or computationally demanding random-shift permutations with discretised p-values; sbivar tests directly on disjoint coordinate sets.

Key results

  1. ShownWith one modality's coordinates permuted on real data (so the other keeps its spatial autocorrelation), Lee's L, bivariate Geary's C, SpaGene and bivariate Moran's I with permutation or SpatialDM nulls had inflated type I error; sbivar's Moran's I and GAMs stayed valid.

    Fig 2, averaged over all Godfrey and Vicari samples, pairs among the 30 most abundant features per modality, 300 permutations for permutation methods. Pearson, modified t-test and random-shift Moran's I also stayed valid; randomisation-null Moran's I was very conservative. Inflation peaked when the metabolome (highest SAC) stayed unpermuted; permuting both modalities kept all methods valid. · Nonparametric simulations; Nonparametric simulations

  2. SuggestedIn simulations with SAC, Pearson, GAMs and Moran's I under randomisation, permutation or SpatialDM nulls were liberal; sbivar's independence-null Moran's I and bivariate Gaussian process (GP) tests, plus the modified t-test, controlled type I error with good power.

    Fig 3, Figs S13-S15: 50 Monte-Carlo runs per scenario with ten features per modality on grids of n=m=625 or n=625, m=900; rates are shares of raw p-values below 0.05. Generalized additive mixed models (GAMMs) also had good power but were slightly liberal in some scenarios; random-shift Moran's I and Kendall's correlation had low power; inflation was worst when both features had SAC. · Parametric simulations; Parametric simulations

  3. SuggestedOn disjoint grids, tests using the original coordinates (Moran's I, bivariate GPs, GAMMs) mostly outpowered tests run after nearest-neighbour matching (modified t-test, Kendall's correlation); simulated alignment error cut the power of many methods but type I error control held.

    Fig 3 disjoint rows (n=625, m=900; matching with RANN nn2). Alignment test: landmark noise with standard deviations 0.5, 1, 1.5, 2 and 3 on Vicari coordinates, 10 repeats per sample and level (Fig S17). Bivariate GPs and especially GAMMs were slowest and scaled worst with spot number; GAMs, modified t-test and sbivar Moran's I were faster (Fig S16). · Parametric simulations; Parametric simulations

  4. SuggestedFor multi-image designs, per-image association measures entered into a linear model kept type I error nominal for Pearson, Moran's I and GAMs; inverse-variance weighting inflated type I error on real-data nulls and did not consistently add power, whereas scaling Moran's I by its maximum raised power.

    Permuted Godfrey data per cancer type and Vicari data across mice (Fig S12); simulations with nine replicates per pattern on disjoint grids with negative-binomial spot numbers (Figs S18-S19). Pearson after spot matching had low power overall; without variance weighting, GAMs were slowest and Pearson fastest. · Generalized additive models (GAMs); Parametric simulations; Nonparametric simulations; Parametric simulations

  5. SuggestedPer-section tests on human lung and breast cancer (Godfrey) found more gene-metabolite pairs than the original absolute-correlation 0.4 cutoff, e.g. KLK10, KLK11, KLK12 and KLK13 with PI 34:1 in LC 091; in mouse striatum (Vicari), CSRP1 tracked palmitoyl-L-carnitine by GAMs and Moran's I.

    Benjamini-Hochberg corrected; KLK pairs recurred more weakly in BC 515 Sections 1 and 2 (Fig S26, Table S1); CSRP1 pair in V11L12-038 A1 and B1 (Fig S35). Some dopamine-gene pairs from the original spot-matched Vicari analysis were not confirmed by GAMs. Replicated analyses (per cancer type; all striatum and substantia nigra samples) gave no significant pairs after correction. · Real data analysis; Real data analysis

Methods, data and reproducibility

Three tests of the independence null (SAC allowed in each variable) on disjoint coordinates: bivariate Moran's I with a variance built from variogram-estimated autocorrelation; a score test for a bivariate Gaussian process (GP) cross-covariance term; and a z-test on the correlation of GAM spline surfaces on a shared grid (delta-method variance). Multi-image data use per-image measures in linear models.

  • Spatial transcriptomics and metabolomics on the same human lung and breast cancer tissue sections (Godfrey data); spots pre-matched by the original authors by aggregating metabolomics spots · human
  • Spatial transcriptomics and metabolomics on the same sections of 6-hydroxydopamine (6-OHDA) unilaterally lesioned mouse brains (Parkinson's model; Vicari data), aligned with MAGPIE · mouse · n = three mouse brains, measured in triplicate
  • Parametric Gaussian simulations on regular joint and disjoint grids (null, dependent GP, streak and zone patterns), plus nine-replicate designs and alignment-error simulations · n = 50 runs per scenario, 10 features per modality; 625 and 900 spots
Methods & tools
sbivar, gstat, mgcv, MAGPIE, SpatialDM, SpaGene, spdep, bispdep, SpatialPack, NTSS, RANN, stats
Reproducibility gaps
No accession numbers or download links are given for the Godfrey and Vicari datasets; Vicari alignment landmarks and resulting coordinates are provided as S1 Data. No software versions or hardware are stated.
Limitations (6)
  • GAMs can still be liberal under short-range spatial autocorrelation, since splines miss fine-grained structure and spline-coefficient variance is underestimated; GAMMs correct this only at a major computational cost. · stated by the authors
  • Bivariate GPs are mathematically and computationally demanding, scale poorly with spot number and give only a p-value without an association-strength estimate, so they are not used for replicated designs. · stated by the authors
  • Being nonparametric, bivariate Moran's I struggles with the heteroskedasticity typical of omics data; it and the GP test do not estimate the interaction range, so several prespecified ranges must be tested and combined. · stated by the authors
  • Power comparisons are inherently incomplete because bivariate spatial association can take infinitely many forms; the authors say additional independent benchmarks are needed. · stated by the authors
  • Replicated real-data analyses were too small to separate real from spurious associations: no pair was significant after multiplicity correction in either dataset. · stated by the authors
  • The highlighted pairs (KLK transcripts with PI 34:1; CSRP1 with palmitoyl-L-carnitine) come from per-section analyses, are not significant in replicated analyses, and are explained via literature on PI3K/AKT signalling and astrocyte lipid handling, with no orthogonal validation; real-data testing covers two datasets. · noted in this summary

Read this if you test colocalisation between transcripts and metabolites (or other modality pairs) measured on mismatched spot grids and need false-positive control under spatial autocorrelation.

Citations
0

Reproducibility ledger

Nobody has reported trying to reproduce this paper yet. Did you replicate an experiment or run the code? Your report saves the next lab a week.

Sign in to file a reproducibility report.

Related papers

Similarity Map
SelectedSame fieldCross-field
THISLi '26Jarrahi '26Li '22Wang '26Liu '22Hao '26
60% semantic matchBioinformatics
Jinpu Li, Mauminah Raina, Yiqing Wang+23 more·bioRxiv·Aug 10, 2026