SVG-Rotation-Bench

Benchmarking Rotation Invariance and Performance of Kernel-Based Methods for Spatially Variable Genes (SVGs) Detection

Overview

Spatially resolved transcriptomics (SRT) measures gene expression while preserving the spatial coordinates of each observation. The detection of genes whose expression shows a non-random pattern across the tissue, also known as Spatially Variable Genes (SVGs), is the entry point of any SRT analysis. Therefore, a mis-specified SVG list propagates errors through the entire pipeline.

The field now offers dozens of SVG detection methods, but different methods produce markedly different SVG lists on the same data.

This repository provides a systematic benchmarking framework that evaluates five representative SVG detection methods along two independent dimensions of quality:

1. Classification performance — the ability to distinguish spatially variable genes from non-variable ones (e.g. false positive rate and sensitivity).
2. Rotation invariance — the consistency of results when the spatial coordinates are rotated.

The framework runs in two modes: a simulated mode with a known ground truth (scDesign3-generated synthetic data with a controlled signal gradient) and a whole mode that runs the rotation benchmark directly on the full, real geneset of each slice.

Background: why rotation invariance matters

The absolute coordinate system (the x and y axes of the slide) carries no biological meaning: the biology of the tissue is the same regardless of the angle at which the section was positioned. Yet methods whose statistical formulation depends on the absolute coordinate system can, in principle, produce different results for the same tissue when the spatial coordinates are rotated.

The central hypothesis of this work is that methods operating on relative distances between spots (neighbor graphs, Euclidean-distance kernels) are rotation-invariant by construction, whereas methods that project coordinates onto absolute axes are vulnerable to rotation.

Su & Cui (2025, Nature Communications) illustrated this problem on the mouse olfactory bulb: SPARK-X identified 2,321 SVGs in the original orientation, a number that dropped to 548 under a 60° rotation, with only seven genes detected consistently between the 30° and 45° rotations. This repository extends that investigation by:

1. Quantifying rotation invariance under controlled simulation with a known ground truth (scDesign3), rather than on real data without a reference.
2. Evaluating both dimensions jointly (invariance and classification calibration), making their independence explicit.
3. Extending the method panel to the random-effects regression family (SpatialDE, nnSVG), previously treated only by theoretical argument.
4. Scaling the evaluation across four tissue slices of the 10x Visium mouse brain dataset, in both simulated and whole-transcriptome modes.

Scope and Included Methods

We benchmark 5 representative methods across two main statistical approaches:

1. Dependence Tests (Evaluating global independence between expression and location):

• SPARK-X: Non-parametric covariance test.
• SMASH: Generalized non-parametric method incorporating distance-based kernels.

2. Random-Effect Regression (Modeling spatial variation as a Gaussian Process):

• SpatialDE: Linear mixed model using standard Gaussian Processes (GP).
• nnSVG: Nearest-Neighbor Gaussian Process (NNGP) for linear scalability.

3. Baseline (Rotation invariant by design, graph based)

• Moran's I: Spatial autocorrelation statistic (via squidpy).

Pipeline

The benchmark is organized into four sequential stages. Each stage iterates over 4 tissue slices of the 10x Visium mouse brain dataset (stxBrain.SeuratData: anterior1, anterior2, posterior1, posterior2) and two modes:

• simulated — synthetic counts generated by scDesign3 from a fitted reference model, with signal strength controlled by a mixing parameter alpha.
• whole — the full, real geneset of the slice, run through the rotation benchmark directly so that rotation-invariance metrics (Jaccard, Venn) can be computed on the real data, in parallel to the simulated runs.

Control what runs by editing the SLICES and MODES vectors at the top of each script.

src/
├── 01_simulation/
│   ├── scripts/
│   │   ├── run_scdesign3_simulation.R  # sim mode: fit scDesign3 + alpha sweep
│   │   ├── export_whole_dataset.R       # whole mode: export full real counts
│   │   └── sanity_simulation_plots.R   # optional sim sanity diagnostics
│   └── outputs/scDesign3/{slice}/{simulated|whole}/
│       ├── data/                       # counts.csv, location.csv
│       └── figures/                    # Diagnostic plots (.png, .pdf)
├── 02_rotation/
│   ├── scripts/
│   │   ├── generate_rotated_locations.R # rotates coordinates
│   │   ├── convert_to_anndata.py       # converts data to .h5ad
│   │   └── sanity_rotation_plots.py    # optional rotation sanity diagnostics
│   └── outputs/{slice}/{mode}/
│       ├── locations/                  # rotated_locations_{angle}.csv
│       └── anndata/
│           ├── data/                   # scdesign3_angle{angle}.h5ad
│           └── figures/                # Verification plots (.png)
├── 03_benchmark/                       # Run SVG detection methods across all rotations
│   ├── scripts/                        # Per-method benchmark scripts
│   ├── tools/                          # External tools (SMASH)
│   └── outputs/{slice}/{mode}/{method}/  # Raw results per slice/mode/method/angle
└── 04_metrics/                         # Collect and compare results per method/rotation
    ├── scripts/
    └── outputs/{slice}/{mode}/
        ├── module1_rotation/           # Jaccard + Venn (both modes)
        └── module2_performance/        # FPR/sensitivity/confusion/runtime (simulated only)

Setup

1. Python environment (conda)

The Python-based methods and data conversion scripts require the following environment:

conda env create -f environment.yml
conda activate svg-rotation-bench

This installs: numpy, scipy, pandas, matplotlib, h5py, scikit-learn, scikit-image, geopandas, scanpy, squidpy, anndata, statsmodels, tqdm, plus the PyPI packages SpatialDE and NaiveDE.

2. R dependencies

The R-based methods and metrics require packages from both CRAN and Bioconductor.

CRAN packages

install.packages(c(
  "ggplot2", "reshape2", "patchwork", "Seurat",
  "cowplot", "dplyr", "scales", "MASS", "ggvenn"
))

Bioconductor packages

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install(c(
  "SingleCellExperiment", "SpatialExperiment", "scran",
  "nnSVG", "anndata"
))

SeuratData (CRAN)

The reference dataset is distributed through SeuratData. The first run of run_scdesign3_simulation.R and export_whole_dataset.R will automatically install the stxBrain dataset if it is missing; you can also do it manually:

install.packages("SeuratData")
library(SeuratData)
InstallData("stxBrain")

The dataset provides four slices (anterior1, anterior2, posterior1, posterior2) loadable via LoadData("stxBrain", type = <slice>). Spatial coordinates are retrieved with GetTissueCoordinates(seu, scale = "lowres").

SPARK-X (GitHub)

SPARK-X is not available on CRAN or Bioconductor. Install from the authors' GitHub repository:

if (!requireNamespace("devtools", quietly = TRUE))
    install.packages("devtools")

devtools::install_github("xzhoulab/SPARK")

scDesign3 (GitHub)

scDesign3 is not available on CRAN or Bioconductor. Install from the authors' GitHub repository:

if (!require("devtools", quietly = TRUE))
    install.packages("devtools")
devtools::install_github("SONGDONGYUAN1994/scDesign3")

3. External tools (manual download)

Two methods require manual cloning/downloading into src/03_benchmark/tools/:

SMASH

Clone the SMASH repository into the tools directory:

cd src/03_benchmark/tools/
git clone https://github.com/sealx017/SMASH-package.git 

4. SpatialDE compatibility patch

SpatialDE is unmaintained and incompatible with scipy >= 1.12 (removed scipy.misc.derivative and scipy.arange). The benchmark script (run_spatialde.py) includes runtime compatibility patches to address these issues. No manual intervention is required.

Stage 1: Simulation (src/01_simulation/)

Synthetic spatial transcriptomics data is generated using scDesign3, fitted per-slice on the 10x Visium mouse brain reference dataset (stxBrain.SeuratData, four slices: anterior1, anterior2, posterior1, posterior2). Each slice is processed independently because slices have distinct coordinate systems and spatial structure.

The full pipeline runs in two modes:

• Marginal Modeling: Learns the true spatial mean $\mu_s(s)$ using a GAM with Gaussian Process smoother.
• Joint Modeling: Preserves gene-gene correlation using a Gaussian Copula.
• Ground Truth Generation: True SVGs are defined by mixing the spatial signal with a randomized null signal $\mu_{ns}(s)$ using a mixing parameter $\alpha$:
  $\mu(s) = \alpha \cdot \mu_s(s) + (1 - \alpha) \cdot \mu_{ns}(s)$

Scripts:

• scripts/run_scdesign3_simulation.R — Per slice: loads stxBrain via SeuratData, pre-filters to top 200 SVGs (Moran's I), fits scDesign3 (GP marginal models + Gaussian copula), selects the top 50 genes by deviance explained, runs the alpha sweep (0.00-1.00 in 0.05 steps), and writes counts.csv/location.csv to outputs/scDesign3/{slice}/simulated/data/. Also produces one real-vs-sim sanity plot per slice.
• scripts/export_whole_dataset.R — Per slice: exports the full stxBrain counts and spatial coordinates to outputs/scDesign3/{slice}/whole/data/ so the benchmark can be run on the complete real geneset.
• scripts/sanity_simulation_plots.R — Optional standalone diagnostic: plots one selected gene (default Mbp) across alpha levels {1, 0.6, 0.2, 0} for each slice.

Stage 2: Rotation (src/02_rotation/)

The generated count matrix is frozen. To test rotation invariance, the original spatial coordinates matrix $S \in \mathbb{R}^{n \times 2}$ is multiplied by a rotation matrix $R$ at specific angles $\theta \in {0^\circ, 30^\circ, 45^\circ, 60^\circ}$:

$$S^* = S R^T$$

Where $R$ is the standard rotation matrix:

$$R = \begin{bmatrix} \cos(\frac{\theta}{180}\pi) & -\sin(\frac{\theta}{180}\pi) \\ \sin(\frac{\theta}{180}\pi) & \cos(\frac{\theta}{180}\pi) \end{bmatrix}$$

Scripts:

• scripts/generate_rotated_locations.R — Reads outputs/scDesign3/{slice}/{mode}/data/location.csv and outputs per-angle rotated coordinate files to outputs/{slice}/{mode}/locations/.
• scripts/convert_to_anndata.py — Builds separate AnnData .h5ad files for each rotation angle, combining the fixed counts.csv with the rotated spatial coordinates. .h5ad files saved to outputs/{slice}/{mode}/anndata/data/. For simulated mode the gene_alpha feature metadata is parsed normally; for whole mode the var metadata simply stores the plain gene symbol.
• scripts/sanity_rotation_plots.py — Optional standalone diagnostic: plots one (slice, mode) at 0/30/45/60 deg side-by-side; can also produce a combined all-slices x all-angles figure when ALL_SLICES_FIGURE = True.

Stage 3: Benchmark (src/03_benchmark/)

Each SVG detection method is implemented as a standalone script under scripts/. Each script iterates over the selected slices and modes and runs the method across all four rotation angles. Raw outputs (p-values, test statistics per gene) are saved to outputs/{slice}/{mode}/{method}/.

Scripts:

• scripts/run_sparkx.R — SPARK-X non-parametric covariance test per slice/mode/angle.
• scripts/run_nnsvg.R — nnSVG nearest-neighbor Gaussian process per slice/mode/angle.
• scripts/run_spatialde.py — SpatialDE Gaussian process regression per slice/mode/angle.
• scripts/run_moransi.py — Moran's I spatial autocorrelation per slice/mode/angle.
• scripts/run_smash.py — SMASH non-parametric kernel test per slice/mode/angle.

Stage 4: Metrics (src/04_metrics/)

Raw benchmark outputs from Stage 3 are evaluated per (slice, mode). Metrics are kept focused on a small, decision-relevant set; per-slice only (no cross-slice aggregation for now).

Module 1 — Rotation Invariance (both simulated and whole modes):

• Jaccard index of significant SVG sets across all angle pairs (0/30/45/60), with a cross-method heatmap.
• Venn diagrams (4-way) of significant SVG sets per method.

Module 2 — Statistical Performance (simulated mode only, since alpha ground truth is required):

• FPR at alpha = 0, 0 deg only.
• Sensitivity by alpha curve at 0 deg.
• Confusion matrix heatmap at 0 deg (adj p < 0.05).

Runtime barplot (both modes).

Scripts:

• scripts/compute_all_metrics.R — Iterates over the selected slices and modes, reads outputs/{slice}/{mode}/{method}/scdesign3_angle{angle}_results.rds, computes the metrics above, and writes CSVs + figures (PNG 300dpi + PDF) to outputs/{slice}/{mode}/{module1_rotation|module2_performance}/. Module 2 is skipped automatically when mode == "whole".