Spatial transcriptomics with R and Python: a Visium workflow

Spatial assays retain a relationship between measured RNA and a position in a tissue section. This Visium-oriented walkthrough follows the same analysis questions in R and Python: load a Space Ranger output, inspect expression on the tissue image, and compare expression-based clusters with anatomy.

Scope and input: both readers below explicitly use filtered_feature_bc_matrix.h5 beside a spatial/ directory. A matrix-market folder alone is not the H5 file these calls expect. For your own conventional Visium v1/v2 output, point both readers at the same outs/ directory with those assets. The Space Ranger output guide documents the files. These calls do not cover Xenium, MERFISH, Visium HD binning, segmentation, deconvolution, or a complete QC workflow.

Prepare one public sample for both languages

Download the standard-library input preparer into a new project. With Python 3.10 or later, run the command below from that project root. It downloads approximately 30 MB from 10x Genomics, checks recorded SHA256 fingerprints and stages the mouse-brain H5 matrix and images together. If an ordinary request receives HTTP 403 or a transient network error, it retries the same official files with size-bounded, validated byte-range requests; the final fingerprints must still match. It refuses an existing destination; do not remove checksum checks if an upstream source changes. No R or Squidpy installation is needed for this preparation step.

python3 prepare_inputs.py --output data/visium/V1_Adult_Mouse_Brain

data/visium/V1_Adult_Mouse_Brain/
  filtered_feature_bc_matrix.h5
  source-manifest.json
  spatial/
    tissue_positions_list.csv
    scalefactors_json.json
    tissue_lowres_image.png
    tissue_hires_image.png

Verified boundary: the source downloads, fingerprints, HDF5 signature, archive contents and shared-directory preparation were checked. Both branches then loaded 2,702 spots and 32,285 input features, preserved raw-count library-size QC, completed expression-neighbor clustering, and rendered the actual image-backed QC, marker and cluster plots. See the recorded source fingerprints and input-preparation notes. Successful execution and inspected overlays do not establish tissue identity, biological cluster labels or a valid scientific conclusion. Before analysis, record versions in the active environments: Python python -c "import scanpy,squidpy; print(scanpy.__version__, squidpy.__version__)"; R R -q -e 'sessionInfo()'. Check official package documentation for compatibility.

1. Load a Visium library

Keep the expression matrix and tissue metadata together. Squidpy's Visium reader loads the count matrix and, when available in the expected output, spot coordinates, image, and scale factors. Seurat's loader similarly expects the matrix file and image data under the specified directory. A missing or mismatched spatial/ directory is a common reason an image-backed plot cannot be made.

R: Seurat

The H5 reader requires hdf5r in the active R library. Install it with install.packages("hdf5r") if missing, then confirm it loads before this branch. On the tested Ubuntu source-install route, HDF5 development files (libhdf5-dev) and an existing R build toolchain were also needed. This is an additional reader dependency, not a complete fresh R/Seurat installation guide.

library(Seurat)
library(ggplot2)
if (!requireNamespace("hdf5r", quietly = TRUE)) stop("Install hdf5r to read the H5 matrix.")

shared_input <- "data/visium/V1_Adult_Mouse_Brain"
spatial_data <- Load10X_Spatial(data.dir = shared_input, filename = "filtered_feature_bc_matrix.h5")
spatial_data <- SCTransform(spatial_data, assay = "Spatial", verbose = FALSE)

SCTransform can take several minutes on this sample without the optional glmGamPoi backend. The tested run used Seurat's slower native fallback; inspect its messages and let the computation finish rather than changing packages midway through a run.

Python: Squidpy and Scanpy

import scanpy as sc
import squidpy as sq

shared_input = "data/visium/V1_Adult_Mouse_Brain"
adata = sq.read.visium(shared_input, counts_file="filtered_feature_bc_matrix.h5")
adata.var_names_make_unique()
sc.pp.calculate_qc_metrics(adata, inplace=True)
sc.pp.normalize_total(adata, inplace=True)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, flavor="seurat", n_top_genes=2000)

Python's normalization and highly-variable-gene selection and Seurat's SCTransform are alternative preprocessing choices; they are not identical transformations or numerically interchangeable. Inspect object and assay/layer names when adapting to another release. If loading fails, confirm the selected directory directly contains filtered_feature_bc_matrix.h5 and spatial/. For your own Space Ranger run this is usually outs/; the prepared public sample uses the named directory above. If counts load without an image, inspect the spatial assets and library metadata before interpreting plots.

2. Plot a feature against tissue context

First plot a basic library-size QC measure, then choose a feature that exists in the object and plot it over the image. The code uses Snap25 and Mbp as example labels only; their presence, biological relevance, and pattern depend on species, assay, tissue, and feature naming. Inspect the image key and color scale before interpreting a spatial gradient.

R

SpatialFeaturePlot(spatial_data, features = c("nCount_Spatial", "Snap25")) &
  theme(legend.position = "right")

The returned R figure combines panels with patchwork. Its & operator applies the legend theme to every panel; using + here only styles the last panel and can leave the QC colorbar cramped above the first plot.

Python

sq.pl.spatial_scatter(adata, color=["total_counts", "Snap25"], img_alpha=0.5, cmap="magma")

Calculate total_counts before normalization and log transformation, as above. It then records raw counts per spot; recomputing it from log-transformed expression would no longer measure library size. The marker panel uses the processed expression matrix, so its values are on a different scale from the raw-count QC panel.

Squidpy's spatial_scatter is the current image-overlay plotting API. Older Scanpy examples use sc.pl.spatial, but current Scanpy documentation marks it deprecated in favor of Squidpy's plot. A quick success check is: nonzero observations and features; nonempty, finite library-size values; a tissue image visible beneath spatially aligned spots; and a readable color scale for QC and the selected marker. Library size is a QC view, not a universal pass/fail threshold. If the image is absent or spots are displaced, verify the input's spatial/ files, library key, and image/coordinate scale; if the marker is missing, choose an observed feature and the intended assay/layer. Do not continue to biological interpretation until these checks pass.

3. Cluster expression profiles, then inspect their locations

Expression-neighbor clustering groups observations by transcriptional similarity. It does not by itself define spatial domains. Build the expression graph, cluster, and then display cluster labels in tissue coordinates. A spatial-neighbor graph is a separate choice for spatial statistics and should be constructed only when the question and coordinate system justify it.

R

spatial_data <- RunPCA(spatial_data, assay = "SCT", verbose = FALSE)
spatial_data <- FindNeighbors(spatial_data, reduction = "pca", dims = 1:30)
spatial_data <- FindClusters(spatial_data, verbose = FALSE)
SpatialDimPlot(spatial_data, label = TRUE, label.size = 3)

Python

sc.pp.pca(adata)
sc.pp.neighbors(adata)
sc.tl.leiden(adata, key_added="clusters")
sq.pl.spatial_scatter(adata, color="clusters", img_alpha=0.4)

These are compact workflow fragments, not a fully specified analysis: package versions, feature filtering, dimensionality choices, clustering resolution, and reproducibility settings need to be chosen and recorded for the dataset. A cluster boundary can reflect cell mixture, tissue composition, depth, or technical structure. Compare the map with morphology and QC before assigning a tissue-domain label.

4. Interpret a pattern without overclaiming

A feature appearing in a region supports a spatial association in this section. Co-location alone does not establish direct cell-cell interaction, lineage, or a causal tissue mechanism. Standard spot-level comparisons also need to account for dependence between nearby observations when making inferential claims. Spatially variable-gene methods and deconvolution address different questions; neither substitutes for checking their assumptions, reference compatibility, and uncertainty.

For spot-based assays such as conventional Visium, each observation can contain RNA from multiple cells. A region-level change may therefore arise from a shift in cell-type mixture, a within-cell expression change, or both. Deconvolution methods estimate mixture from a reference and add assumptions; their estimates are not direct cell counts. Imaging-based assays and newer Visium formats have different resolution and input contracts, so do not transfer spot-level interpretations automatically.

Practical check

For one conventional Visium library, record the input directory and versions using the commands above, run one language's loading and plotting fragments, and compare a library-size QC measure and an observed marker with the tissue image. Your check passes only when counts/features are nonempty, the image and spots are visibly aligned, and both color panels are interpretable; investigate the path, spatial metadata, image selection, feature name, or assay/layer if not. In a short note, state the coordinate/image choice, whether the pattern is qualitative or tested statistically, at least one plausible technical or composition-related alternative, and what additional evidence your interpretation would need.