Introduction

This vignette demonstrates realistic workflows using MultiAssaySpatialExperiment (MASE) for multi-assay spatial transcriptomics analysis. It is divided into two parts:

Part 1: Data import: Reading spatial omics data from vendor outputs (Xenium, Visium, Visium HD, CosMx, MERSCOPE) and converting from other Bioconductor spatial formats.

Part 2: Real-world workflows: Complete analysis examples integrating multiple technologies, spatial transformations, and cross-platform label transfer.

Scientific question: How do cell-level and tissue-level spatial patterns complement each other, and can high-resolution Xenium data help interpret lower-resolution Visium spots?

library(MultiAssaySpatialExperiment)
library(SummarizedExperiment)
library(SingleCellExperiment)
library(S4Vectors)
library(sf)

# For visualization
library(ggplot2)
set.seed(123)

Part 1: Data import

MASE provides readers for major spatial omics platforms and coercion methods for interoperability with other Bioconductor spatial classes.

Reading Xenium data

10x Genomics Xenium provides single-cell resolution subcellular imaging. The Xenium Explorer output directory contains:

cell_feature_matrix.h5: Cell × gene counts
cells.parquet or cells.csv.gz: Cell coordinates and metadata
nucleus_boundaries.parquet: Cell segmentation polygons (optional)
transcripts.parquet: Transcript-level coordinates (optional, large)

Basic usage

Two arguments matter more than the rest. segmentations chooses whether to load the cell boundaries, the nucleus boundaries, or neither, and it is what decides whether the object gets a shapes layer at all. add_transcripts controls the transcript table, which is by far the largest part of a Xenium run and is worth skipping until you need molecule-level detail.

mase_xenium <- readXeniumMASE(
  data_dir = "path/to/xenium_output",
  sample_id = "sample1",
  segmentations = "cell",
  add_transcripts = FALSE  # Skip transcripts for faster loading
)

The result holds the cell-by-gene matrix as a SingleCellExperiment in experiments(mase)$rna, the cell centroids in spatialPoints(), and, when segmentations is set, the cell or nucleus boundaries in spatialShapes(). spatialMap() records which spatial element each assay column sits at.

Options

Parameter	Default	Purpose
`data_dir`	required	Path to Xenium output directory
`sample_id`	`NULL`	Sample identifier (uses directory name if NULL)
`segmentations`	`"cell"`	Segmentation to load: `"cell"`, `"nucleus"`, or `"both"`
`add_transcripts`	`FALSE`	Load transcript-level coordinates (large!)
`images`	`TRUE`	Load image metadata
`load_images`	`FALSE`	Load image raster data

When to load transcripts

Transcript-level data (millions of points) is useful for visualizing subcellular localization, quality control (transcript counts per cell), or custom aggregation strategies.

Skip transcripts (add_transcripts = FALSE) for faster loading and smaller memory footprint when working with cell-level data only.

Reading Visium data

10x Genomics Visium provides spot-based spatial transcriptomics. Space Ranger output contains:

filtered_feature_bc_matrix.h5: Spot × gene counts
spatial/tissue_positions.csv: Spot coordinates
spatial/scalefactors_json.json: Image scale factors
spatial/tissue_hires_image.png: H&E image (optional)

Basic usage

data selects the filtered matrix (spots under tissue) or the raw one (every spot on the slide); filtered is what most analyses want. images = TRUE loads the H&E image, which is what makes imgData() and the image-aware plotting useful later, at the cost of reading the image file.

mase_visium <- readVisiumMASE(
  data_dir = "path/to/spaceranger_output",
  sample_id = "sample1",
  images = TRUE,
  data = "filtered"
)

The result holds the spot-by-gene matrix in experiments() and the spot coordinates in spatialPoints(). With images = TRUE the H&E image is available through spatialImages(), and imgData() records which image belongs to which specimen along with its scale factor.

Options

Parameter	Default	Purpose
`data_dir`	required	Path to Space Ranger output directory
`sample_id`	`NULL`	Sample identifier (uses directory name if NULL)
`type`	`"HDF5"`	Matrix format: `"HDF5"` or `"sparse"`
`data`	`"filtered"`	Matrix to read: `"filtered"` (spots under tissue) or `"raw"`
`images`	`TRUE`	Load image metadata
`load_images`	`FALSE`	Load image raster data
`unit`	`"pixel"`	Coordinate units: `"pixel"` or `"micron"`
`min_area`	`NULL`	Minimum polygon area for filtering
`block_size`	100 MB	Block size for chunked reading

Visium HD support

Visium HD provides higher resolution with multiple binning levels (2 µm, 8 µm, 16 µm). Use the dedicated HD reader and specify bin size codes "002", "008", or "016":

# Read Visium HD at 8 µm resolution
mase_visium_hd <- readVisiumHDMASE(
  data_dir = "path/to/spaceranger_hd_output",
  sample_id = "sample_hd",
  bin_size = "008",
  images = TRUE
)

Reading CosMx data

NanoString CosMx provides subcellular resolution with compartment-specific segmentation (nucleus, membrane, cytoplasm). Output files include:

*_exprMat_file.csv: Expression matrix
*_metadata_file.csv: Cell metadata
*_tx_file.csv: Transcript coordinates
CellComposite/: Composite images

Basic usage

CosMx runs are organized into fields of view (FOVs), and a whole slide can be more than you want in memory. fov_ids restricts the read to named FOVs; leaving it NULL loads all of them.

mase_cosmx <- readCosMxMASE(
  data_dir = "path/to/cosmx_output",
  sample_id = "cosmx_sample",
  fov_ids = NULL  # NULL loads all FOVs; or c("1", "3", "5")
)

As with the other readers, the counts land in experiments(), the cell coordinates in spatialPoints(), and the segmentation in spatialShapes(), which for CosMx includes subcellular compartments when the run provides them.

Options

Parameter	Default	Purpose
`data_dir`	required	Path to CosMx output directory
`sample_id`	`NULL`	Sample identifier (uses directory name if NULL)
`fov_ids`	`NULL`	FOVs to load (`NULL` loads all)
`load_transcripts`	`FALSE`	Load transcript-level coordinates
`images`	`TRUE`	Load image metadata
`load_images`	`FALSE`	Load image raster data
`min_area`	`NULL`	Minimum polygon area for filtering

FOV (Field of View) handling

CosMx data is organized by FOVs. Options:

fov_ids = NULL: Load all FOVs (combines into single MASE)
fov_ids = c("1", "3", "5"): Load specific FOVs
fov_ids = "1": Load a single FOV

Reading MERSCOPE data (MERFISH)

Vizgen MERSCOPE provides high-throughput spatial transcriptomics (often referred to as MERFISH). Output directory contains:

cell_by_gene.csv: Cell × gene matrix
cell_metadata.csv: Cell coordinates and metadata
detected_transcripts.csv: Transcript coordinates
cell_boundaries/: Cell segmentation polygons (Parquet)

Basic usage

segmentation names which segmentation run to load, since a MERSCOPE output directory can hold more than one; "cellpose" is the usual default. As with Xenium, load_transcripts gates the largest table in the run.

mase_merscope <- readMERSCOPEMASE(
  data_dir = "path/to/vizgen_output",
  sample_id = "merscope_sample",
  segmentation = "cellpose",
  load_transcripts = FALSE
)

MERSCOPE runs are often three-dimensional, so spatialPoints() carries a z column when the data provide one; the segmentation polygons are in spatialShapes().

Options

Parameter	Default	Purpose
`data_dir`	required	Path to MERSCOPE output directory
`sample_id`	`NULL`	Sample identifier (uses directory name if NULL)
`fov_ids`	`NULL`	FOVs to load (`NULL` loads all)
`segmentation`	`"cellpose"`	Segmentation to load: `"cellpose"`, `"watershed"`, or `"both"`
`load_transcripts`	`FALSE`	Load transcript-level coordinates
`images`	`TRUE`	Load image metadata
`load_images`	`FALSE`	Load image raster data
`min_area`	`NULL`	Minimum polygon area for filtering
`min_qv`	20	Minimum quality value for transcripts

Choosing a reader

Select the reader that matches your platform output directory:

Platform	Reader
10x Xenium	`readXeniumMASE()`
10x Visium	`readVisiumMASE()`
10x Visium HD	`readVisiumHDMASE()`
NanoString CosMx	`readCosMxMASE()`
Vizgen MERSCOPE	`readMERSCOPEMASE()`

Converting from other Bioconductor classes

If the data are already in a single-assay spatial class, there is no need to go back to the vendor output. Coercion is the usual way into MASE when you want to add a second assay to an analysis that started with just one.

From SpatialExperiment

library(SpatialExperiment)

# spe <- read10xVisium("path/to/visium")
mase_from_spe <- as(spe, "MultiAssaySpatialExperiment")

The assays become the single element of the ExperimentList, colData carries over as specimen metadata, and spatialCoords becomes a points layer. The sampleMap and spatialMap are generated for you, so the result is immediately a valid MASE with one assay: the point of the exercise is that you can now add a second one alongside it.

From SpatialFeatureExperiment

library(SpatialFeatureExperiment)

# sfe <- read10xVisiumSFE("path/to/visium")
mase_from_sfe <- as(sfe, "MultiAssaySpatialExperiment")

SpatialFeatureExperiment already carries geometry, so more survives the trip: every assay is preserved, each colGeometry becomes a points or shapes layer depending on whether it holds centroids or polygons, and the annotGeometries are kept as further shapes layers. The generated spatialMap is what ties those geometries back to the assay columns they describe.

Reading multiple samples

For multi-sample experiments, read each sample separately and combine:

mase1 <- readXeniumMASE(data_dir = "path/to/sample1", sample_id = "P001_xenium")
mase2 <- readXeniumMASE(data_dir = "path/to/sample2", sample_id = "P002_xenium")

mase_multi <- c(mase1, mase2)

Each sample_id becomes a specimen in the combined object’s colData, and both sets of spatial layers are carried across under their own names. c() also rebuilds sampleMap and spatialMap over the union, which is what makes the result subsettable by specimen afterwards: asking for one patient’s data returns their assay columns and their geometry, and nothing from the other sample.

Part 2: Real-world workflows

Part 1 got data into the object. The workflows below are what you do once it is there, and each is chosen to exercise a different reason for using a multi-assay container: combining two platforms over one tissue, reconciling coordinate systems that do not agree, keeping many specimens straight, and the shape a typical analysis session takes end to end.

Unlike Part 1, these examples run on simulated data, so they execute when the vignette is built and you can modify them without vendor output to hand.

Workflow 1: Xenium + Visium integration

Background

A common experimental design combines Xenium, which gives single-cell resolution and subcellular localization but covers a limited gene panel (~300-500 genes), with Visium, which gives spot-level resolution (55 µm spots) across the whole transcriptome, at the cost of multiple cells per spot. Integrating both makes it possible to: 1. Use Xenium to identify cell types and their spatial organization 2. Use Visium to capture broader transcriptomic programs 3. Transfer cell type labels from Xenium to Visium spots 4. Deconvolve Visium spots using Xenium cell proportions

Data preparation

This example uses synthetic data that mimics a real experiment. In practice, you would load your own Xenium cell-by-gene matrix and cell centroids, and your own Visium spot-by-gene matrix and spot centroids.

# Xenium data: single-cell resolution
# Simulate 500 cells, 50 genes (subset of full panel)
n_cells <- 500
n_genes_xenium <- 50

xenium_counts <- matrix(
  rpois(n_cells * n_genes_xenium, lambda = 5),
  nrow = n_genes_xenium,
  ncol = n_cells,
  dimnames = list(
    paste0("Gene", 1:n_genes_xenium),
    paste0("Cell", 1:n_cells)
  )
)

# Cell centroids (tissue coordinates in micrometers)
xenium_coords <- DataFrame(
  x = runif(n_cells, 0, 1000),
  y = runif(n_cells, 0, 1000),
  instance_id = paste0("Cell", 1:n_cells)
)

# Cell type annotations (from clustering or known markers)
cell_types <- sample(
  c("Neuron", "Astrocyte", "Microglia", "Oligodendrocyte"),
  n_cells,
  replace = TRUE,
  prob = c(0.5, 0.2, 0.15, 0.15)
)

# Create SingleCellExperiment
xenium_sce <- SingleCellExperiment(
  assays = list(counts = xenium_counts),
  colData = DataFrame(
    cell_type = cell_types,
    tech = "xenium"
  )
)

# Visium data: spot-level resolution
# Simulate 100 spots, 1000 genes (whole transcriptome, downsampled)
n_spots <- 100
n_genes_visium <- 1000

# Create gene names: first 50 overlap with Xenium, rest are Visium-only
visium_gene_names <- c(
  paste0("Gene", 1:n_genes_xenium),  # Overlapping genes (50)
  paste0("GeneV", 1:(n_genes_visium - n_genes_xenium))  # Visium-only (950)
)

visium_counts <- matrix(
  rpois(n_spots * n_genes_visium, lambda = 10),
  nrow = n_genes_visium,
  ncol = n_spots,
  dimnames = list(
    visium_gene_names,
    paste0("Spot", 1:n_spots)
  )
)

# Spot centroids (same coordinate system as Xenium)
visium_coords <- DataFrame(
  x = runif(n_spots, 0, 1000),
  y = runif(n_spots, 0, 1000),
  instance_id = paste0("Spot", 1:n_spots)
)

# Create SummarizedExperiment
visium_se <- SummarizedExperiment(
  assays = list(counts = visium_counts),
  colData = DataFrame(
    tech = rep("visium", n_spots),
    row.names = paste0("Spot", 1:n_spots)
  )
)

Build MASE object

Combine Xenium and Visium from the same tissue section into one MASE object. Both assays share a single specimen (colData row); individual cells and spots are observations linked through sampleMap and spatialMap.

# Specimen metadata: one row per biological unit (here, one tissue section)
specimen_id <- "section_01"
specimens <- DataFrame(
  patient_id = "P001",
  tissue = "mouse_brain",
  section_id = specimen_id,
  row.names = specimen_id
)

n_total <- n_cells + n_spots

# Sample map: link each assay column (cell or spot) to the shared specimen
sample_map <- DataFrame(
  assay = factor(
    c(rep("xenium", n_cells), rep("visium", n_spots)),
    c("xenium", "visium")
  ),
  primary = rep(specimen_id, n_total),
  colname = c(paste0("Cell", 1:n_cells), paste0("Spot", 1:n_spots))
)

# Spatial map: link assay columns to spatial coordinates
spatial_map <- DataFrame(
  assay = factor(
    c(rep("xenium", n_cells), rep("visium", n_spots)),
    c("xenium", "visium")
  ),
  colname = c(paste0("Cell", 1:n_cells), paste0("Spot", 1:n_spots)),
  element_type = rep("points", n_total),
  region = factor(
    c(rep("xenium_coords", n_cells), rep("visium_coords", n_spots)),
    c("xenium_coords", "visium_coords")
  ),
  instance_id = c(paste0("Cell", 1:n_cells), paste0("Spot", 1:n_spots))
)

# Construct MASE
mase <- MultiAssaySpatialExperiment(
  experiments = ExperimentList(
    xenium = xenium_sce,
    visium = visium_se
  ),
  colData = specimens,
  sampleMap = sample_map,
  points = PointsLayerList(
    xenium_coords = xenium_coords,
    visium_coords = visium_coords
  ),
  spatialMap = spatial_map
)

mase
#> A MultiAssaySpatialExperiment object of 2 listed
#>  experiments with user-defined names and respective classes.
#>  Containing an ExperimentList class object of length 2:
#>  [1] xenium: SingleCellExperiment with 50 rows and 500 columns
#>  [2] visium: SummarizedExperiment with 1000 rows and 100 columns
#> Functionality:
#>  experiments() - obtain the ExperimentList instance
#>  colData() - the primary/phenotype DataFrame
#>  sampleMap() - the sample coordination DataFrame
#>  `$`, `[`, `[[` - extract colData columns, subset, or experiment
#>  *Format() - convert into a long or wide DataFrame
#>  assays() - convert ExperimentList to a SimpleList of matrices
#>  exportClass() - save data to flat files
#> Spatial elements:
#>   spatialImages: 0 elements
#>   spatialLabels: 0 elements
#>   spatialPoints: 2 elements
#>   spatialShapes: 0 elements
#>   imgData: NULL
#>   spatialMap: present

Exploratory analysis

With both assays in one object, the first thing worth checking is that they describe the same piece of tissue. Everything downstream, transferring labels, deconvolving spots, comparing expression, assumes the two coordinate systems agree, so it is cheap to confirm that before relying on it.

Do the two assays overlap in space?

Plotting the two side by side on the same axes shows both the resolution difference, thousands of cells against a coarse grid of spots, and whether they cover the same region:

# Extract coordinates
xenium_pts <- spatialPoint(mase, "xenium_coords")
visium_pts <- spatialPoint(mase, "visium_coords")

# Plot both assays
par(mfrow = c(1, 2))

# Xenium cells colored by type
# widen the x range so the legend sits beside the points rather than on top of them
plot(xenium_pts$x, xenium_pts$y,
     col = as.factor(cell_types),
     pch = 16, cex = 0.5,
     main = "Xenium: Single cells",
     xlab = "x (µm)", ylab = "y (µm)",
     xlim = c(min(xenium_pts$x), max(xenium_pts$x) * 1.35))
legend("topright", legend = levels(as.factor(cell_types)),
       col = seq_along(levels(as.factor(cell_types))), pch = 16,
       cex = 0.7, bg = "white", box.col = "grey70")

# Visium spots
plot(visium_pts$x, visium_pts$y,
     pch = 21, bg = "lightblue", cex = 2,
     main = "Visium: Spots (55 µm)",
     xlab = "x (µm)", ylab = "y (µm)")

Which genes are shared between assays

The two platforms target different panels, so the first thing to establish is how much of the gene space they have in common.

xenium_genes <- rownames(experiments(mase)$xenium)
visium_genes <- rownames(experiments(mase)$visium)

shared_genes <- intersect(xenium_genes, visium_genes)

c(shared = length(shared_genes),
  xenium_only = length(setdiff(xenium_genes, visium_genes)),
  visium_only = length(setdiff(visium_genes, xenium_genes)))
#>      shared xenium_only visium_only 
#>          50           0         950

Spatial annotation: Cell types in Visium spots

Use Xenium cell types to annotate Visium spots. For each spot, find which cells fall within it and calculate cell type proportions. The loop below is explicit for teaching; when spot polygons are stored in spatialShapes(), you can also use annotateWithRegions() (see Working with MultiAssaySpatialExperiment).

# Create circular regions for Visium spots (55 µm diameter = 27.5 µm radius)
visium_circles <- st_buffer(
  st_sf(visium_pts,
        geometry = st_as_sfc(paste0("POINT(", visium_pts$x, " ", visium_pts$y, ")"))),
  dist = 27.5
)

visium_shapes <- DataFrame(
  geometry = st_geometry(visium_circles),
  instance_id = visium_pts$instance_id
)

# Add shapes to MASE
spatialShapes(mase) <- ShapesLayerList(visium_spots = visium_shapes)

# For each Visium spot, find which Xenium cells are inside
xenium_sf <- st_sf(xenium_pts,
                   geometry = st_as_sfc(paste0("POINT(", xenium_pts$x, " ", xenium_pts$y, ")")))
visium_sf <- st_sf(
  instance_id = visium_shapes$instance_id,
  geometry = visium_shapes$geometry
)

# Spatial intersection
intersections <- st_intersects(visium_sf, xenium_sf)

# Calculate cell type proportions per spot
spot_cell_types <- lapply(seq_along(intersections), function(i) {
  cell_indices <- intersections[[i]]
  if (length(cell_indices) == 0) {
    return(data.frame(
      spot = visium_pts$instance_id[i],
      Neuron = 0, Astrocyte = 0, Microglia = 0, Oligodendrocyte = 0,
      n_cells = 0
    ))
  }
  types <- cell_types[cell_indices]
  props <- base::table(types) / length(types)
  data.frame(
    spot = visium_pts$instance_id[i],
    Neuron = as.numeric(props["Neuron"]),
    Astrocyte = as.numeric(props["Astrocyte"]),
    Microglia = as.numeric(props["Microglia"]),
    Oligodendrocyte = as.numeric(props["Oligodendrocyte"]),
    n_cells = length(cell_indices),
    stringsAsFactors = FALSE
  )
})

spot_annotations <- do.call(rbind, spot_cell_types)
spot_annotations[is.na(spot_annotations)] <- 0

head(spot_annotations)
#>    spot Neuron Astrocyte Microglia Oligodendrocyte n_cells
#> 1 Spot1    1.0         0       0.0               0       1
#> 2 Spot2    0.0         0       0.0               0       0
#> 3 Spot3    0.0         1       0.0               0       1
#> 4 Spot4    0.0         0       0.0               0       0
#> 5 Spot5    0.0         0       0.0               0       0
#> 6 Spot6    0.5         0       0.5               0       2

Reading the composition result

Each spot now has a cell-type composition rather than a single label. The plots below show, per spot, the proportion of each Xenium cell type it contains. This is the quantity a deconvolution method would otherwise have to estimate; here it is counted directly, because the two assays share a coordinate system and the cells inside each spot are known.

# Plot cell type proportions per spot
par(mfrow = c(2, 2))

for (cell_type in c("Neuron", "Astrocyte", "Microglia", "Oligodendrocyte")) {
  plot(visium_pts$x, visium_pts$y,
       pch = 21, cex = 3,
       bg = rgb(spot_annotations[[cell_type]], 0, 0, alpha = 0.7),
       main = paste(cell_type, "proportion"),
       xlab = "x (µm)", ylab = "y (µm)")
}

Spatial aggregation: Summarize Xenium by regions

Aggregate Xenium gene expression by tissue regions (e.g., cortical layers). First, define regions, then use aggregateByRegion().

# Define two regions: "cortex" (y > 500) and "subcortex" (y <= 500)
# In practice, these would be anatomical annotations
region_polygons <- list(
  cortex = st_polygon(list(cbind(
    c(0, 1000, 1000, 0, 0),
    c(500, 500, 1000, 1000, 500)
  ))),
  subcortex = st_polygon(list(cbind(
    c(0, 1000, 1000, 0, 0),
    c(0, 0, 500, 500, 0)
  )))
)

regions_sf <- st_sf(
  region_id = c("cortex", "subcortex"),
  geometry = st_sfc(region_polygons)
)

# Convert to DataFrame for MASE
anatomical_regions_df <- DataFrame(
  geometry = st_geometry(regions_sf),
  instance_id = c("cortex", "subcortex")
)

# Add to MASE, as a shapes layer named "anatomical_regions" (not "regions": that
# name is reserved for spatialMap's own `region` column, which names the points
# layer each row's instance_id comes from, "xenium_coords" here)
spatialShapes(mase) <- ShapesLayerList(
  visium_spots = visium_shapes,
  anatomical_regions = anatomical_regions_df
)

# Annotate Xenium centroids with anatomical regions, then aggregate expression.
# Run on the full MASE so the anatomical_regions layer is retained
# (subsetByAssay drops shapes not referenced in spatialMap).
mase <- annotateWithRegions(mase,
    points = "xenium_coords", shapes = "anatomical_regions")
xenium_by_region <- aggregateByRegion(mase, by = "anatomical_regions",
    assays = "xenium", FUN = "sum")
xenium_by_region[["xenium"]][1:5, ]
#>       cortex subcortex
#> Gene1   1256      1285
#> Gene2   1215      1220
#> Gene3   1311      1248
#> Gene4   1194      1276
#> Gene5   1171      1286

Comparing expression between the two platforms

With the shared genes identified, the question is whether the two platforms agree on the genes they both measure.

# Extract counts for shared genes
xenium_shared <- assay(experiments(mase)$xenium, "counts")[shared_genes, ]
visium_shared <- assay(experiments(mase)$visium, "counts")[shared_genes, ]

# Calculate mean expression per gene
xenium_means <- rowMeans(xenium_shared)
visium_means <- rowMeans(visium_shared)

# Plot correlation
plot(log1p(xenium_means), log1p(visium_means),
     pch = 16, col = rgb(0, 0, 0, 0.5),
     xlab = "log1p(mean Xenium counts)",
     ylab = "log1p(mean Visium counts)",
     main = "Gene expression correlation")
abline(0, 1, col = "red", lty = 2)

# Correlation
cor_val <- cor(xenium_means, visium_means, method = "spearman")
text(0.5, max(log1p(visium_means)) * 0.9,
     paste("Spearman rho =", round(cor_val, 3)),
     pos = 4)

Integration insights

Key findings from this analysis:

Spatial resolution: Xenium provides single-cell resolution, while Visium captures tissue-level patterns with multiple cells per spot.
Gene overlap: Shared genes allow direct comparison and validation between technologies.
Cell type transfer: Xenium cell type information can be transferred to Visium spots, revealing the cellular composition that a deconvolution method would otherwise have to estimate.
Regional patterns: Both assays can be aggregated by anatomical regions for comparative analysis.
Complementary data: Xenium’s high resolution complements Visium’s broad transcriptomic coverage.

Workflow 2: Coordinate transformation and alignment

When integrating data from different spatial platforms, coordinate systems often differ. This workflow demonstrates how to align Xenium and Visium data using spatial transformations.

Background

Xenium and Visium may use different coordinate systems: Xenium reports microns (µm) from a reference origin, while Visium reports grid indices (row/column) or pixel coordinates. To overlay them, you need to:

Convert Visium grid indices to physical coordinates (microns)
Apply transformations (translation, scaling, rotation) to align modalities
Store transformed coordinates in MASE

Converting Visium grid to microns

Visium spots are arranged in a hexagonal grid. Each spot is ~55 µm diameter, with center-to-center spacing of ~100 µm.

# Visium spot coordinates (grid indices from Space Ranger)
visium_grid <- DataFrame(
  row_idx = c(0, 0, 1, 1, 2, 2),
  col_idx = c(0, 1, 0, 1, 0, 1),
  instance_id = paste0("Spot", 1:6)
)

# Convert to microns (approximate)
# Spot spacing: ~100 µm horizontally, ~87 µm vertically (hexagonal)
visium_microns <- DataFrame(
  x = visium_grid$col_idx * 100,  # Column → x
  y = visium_grid$row_idx * 87,   # Row → y (hexagonal offset)
  instance_id = visium_grid$instance_id
)

visium_microns
#> DataFrame with 6 rows and 3 columns
#>           x         y instance_id
#>   <numeric> <numeric> <character>
#> 1         0         0       Spot1
#> 2       100         0       Spot2
#> 3         0        87       Spot3
#> 4       100        87       Spot4
#> 5         0       174       Spot5
#> 6       100       174       Spot6

Applying affine transformations

MASE allows storing the original and transformed coordinates side by side as named layers of one object that stay attached to the assays through spatialMap; that step comes right after this one, in Storing transformed coordinates. Computing the transform itself is plain R: an affine transform is a matrix multiply, and MASE does not impose a particular package or framework for computing it.

An affine transform composes translation, rotation and scaling into a single matrix:

# Define affine matrix: translate by (50, 30), scale by 1.2, no rotation
# Format: 2x3 matrix [a b c; d e f] where:
#   x' = a*x + b*y + c
#   y' = d*x + e*y + f

affine_matrix <- matrix(c(
  1.2, 0.0, 50,   # x' = 1.2*x + 0*y + 50 (scale + translate)
  0.0, 1.2, 30    # y' = 0*x + 1.2*y + 30 (scale + translate)
), nrow = 2, byrow = TRUE)

# Apply transformation
apply_affine <- function(pts, mat) {
  x_new <- mat[1, 1] * pts$x + mat[1, 2] * pts$y + mat[1, 3]
  y_new <- mat[2, 1] * pts$x + mat[2, 2] * pts$y + mat[2, 3]
  DataFrame(x = x_new, y = y_new, instance_id = pts$instance_id)
}

# Simulate Xenium coordinates
xenium_original <- DataFrame(
  x = c(100, 150, 200, 250, 300),
  y = c(100, 150, 200, 150, 100),
  instance_id = paste0("Cell", 1:5)
)

# Transform Xenium to align with Visium
xenium_aligned <- apply_affine(xenium_original, affine_matrix)

xenium_aligned
#> DataFrame with 5 rows and 3 columns
#>           x         y instance_id
#>   <numeric> <numeric> <character>
#> 1       170       150       Cell1
#> 2       230       210       Cell2
#> 3       290       270       Cell3
#> 4       350       210       Cell4
#> 5       410       150       Cell5

Landmark-based alignment

When you have known landmarks (tissue features visible in both modalities), use them to compute an optimal transformation:

# Example: three landmarks in both Xenium and Visium
# In practice, these would be manually identified features (blood vessels, etc.)

xenium_landmarks <- matrix(c(
  100, 100,
  200, 150,
  300, 200
), ncol = 2, byrow = TRUE)

visium_landmarks <- matrix(c(
  50, 50,
  150, 100,
  250, 150
), ncol = 2, byrow = TRUE)

The workflow around those landmarks is the same whichever package computes the transformation. Identify at least three corresponding points in both modalities, then fit a transformation to them: Procrustes analysis for a rigid fit, least squares for a full affine one, or a dedicated registration package such as RNiftyReg, imager or EBImage. Apply the result to every point, check it visually, and store the transformed coordinates back into the object as their own layer.

Storing transformed coordinates

Store both original and aligned coordinates as separate layers for flexibility:

# Create MASE with both coordinate systems
# Original Xenium coordinates
pts_original <- PointsLayerList(xenium_original = xenium_original)

# Aligned Xenium coordinates
pts_aligned <- PointsLayerList(xenium_aligned = xenium_aligned)

# Visium coordinates
pts_visium <- PointsLayerList(visium_coords = visium_microns)

Keeping the original and aligned coordinates as separate layers of one object is deliberate: the two can be compared directly, the transformation can be refined and re-applied without re-reading anything, and downstream code chooses a coordinate system by naming a layer rather than by tracking which object holds which version.

Checking the alignment visually

Plotting the two modalities before and after the transformation is the quickest way to tell whether the alignment is sound.

par(mfrow = c(1, 2))

# Before alignment
plot(xenium_original$x, xenium_original$y,
     pch = 16, col = "red", cex = 1.5,
     main = "Before alignment",
     xlab = "x (µm)", ylab = "y (µm)",
     xlim = c(0, 400), ylim = c(0, 350))
points(visium_microns$x, visium_microns$y,
       pch = 21, bg = "blue", cex = 3)
legend("topright", legend = c("Xenium (cells)", "Visium (spots)"),
       col = c("red", "blue"), pch = c(16, 21), pt.cex = c(1.5, 2))

# After alignment
plot(xenium_aligned$x, xenium_aligned$y,
     pch = 16, col = "red", cex = 1.5,
     main = "After alignment",
     xlab = "x (µm)", ylab = "y (µm)",
     xlim = c(0, 400), ylim = c(0, 350))
points(visium_microns$x, visium_microns$y,
       pch = 21, bg = "blue", cex = 3)
legend("topright", legend = c("Xenium (aligned)", "Visium (spots)"),
       col = c("red", "blue"), pch = c(16, 21), pt.cex = c(1.5, 2))

Workflow 3: Multi-sample analysis

Once several patients or tissue sections are involved, the thing that usually goes wrong is bookkeeping: which columns of which assay belong to which patient, and which geometry belongs with them. Handling that by hand across a list of objects is where orphaned coordinates and mismatched subsets come from.

In MASE the answer is that specimens are first-class. Each patient is a row of colData, and every assay column and every spatial element is tied back to one of those rows through sampleMap and spatialMap.

specimens_multi <- DataFrame(
  patient_id = c(rep("P001", 150), rep("P002", 120), rep("P003", 180)),
  tissue = "mouse_brain",
  condition = c(rep("control", 150), rep("treated", 120), rep("treated", 180)),
  row.names = paste0("Sample", 1:450)
)

# sampleMap is built exactly as in the single-sample example above, with one
# row per assay column across all three patients.

Because those maps exist, the questions you actually want to ask are ordinary subsetting operations, and each returns a complete, self-consistent object rather than a matrix that has lost track of its coordinates:

mase[, colData(mase)$patient_id == "P001", ]   # one patient, every assay
mase[, colData(mase)$condition == "control", ] # one arm of the study
subsetByAssay(mase, "xenium")                  # one modality, every patient

Comparing conditions then means running the same per-patient analysis over those subsets, with the confidence that each carries its own geometry.

Workflow 4: Iterative analysis pattern

The three workflows above each made a specific point. This one is the shape most sessions actually take, and the reason to show it is the ordering: every step narrows the object, and doing them in this order means each subsequent step runs on less data while still being correct.

Start by reading the vendor output, then drop the specimens that failed quality control. Doing QC first is a habit worth keeping, because everything after it is more expensive per observation:

mase <- readXeniumMASE(data_dir = "path/to/xenium_output", sample_id = "sample1")
mase_qc <- mase[, colData(mase)$qc_pass, ]

Restricting to a region of interest is the next cut. It is a spatial operation rather than an assay one, but because the links are maintained it narrows the assays too, so the annotation step that follows sees only the cells you care about:

library(sf)

roi <- st_polygon(list(matrix(
  c(100, 200, 500, 200, 500, 600, 100, 600, 100, 200),
  ncol = 2, byrow = TRUE)))
mase_roi <- subsetByPolygon(mase_qc, roi)

Then annotate and aggregate. These two belong together and in this order: aggregateByRegion() summarises over an annotation, so it needs annotateWithRegions() to have recorded one on spatialMap first. Running them on the narrowed object rather than the original is the payoff for the two cuts above:

mase_annotated <- annotateWithRegions(mase_roi,
    points = "centroids", shapes = "cells")

cell_expr <- aggregateByRegion(mase_annotated, by = "cells", FUN = "sum")

cell_expr is an ordinary list of matrices, so differential expression, clustering and the rest proceed with whatever tools you normally use. The object is the thing that kept the coordinates, the specimens and the counts aligned up to that point; it does not need to be involved afterwards.

Summary

Across these workflows the same property did the work each time. Two platforms could be compared because both were mapped onto one coordinate system; a region of interest could be cut without losing track of which cells it contained; three patients could be held in one object and separated again on demand. In each case what made it straightforward was that the links between assay columns, specimens and spatial elements are recorded once, in sampleMap and spatialMap, and are maintained by the operations rather than by the analyst.

That is also the limit of what the container claims. It does not implement deconvolution, registration or differential expression; Workflow 2 computed its coordinate transform in plain R for that reason, with MASE only storing the result. What it offers is that the inputs to those methods, and their outputs, stay consistent with each other and with the tissue they came from, and that a single-assay slice can be handed to SpatialExperiment or SpatialFeatureExperiment whenever a tool expects one.

Session info

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats4    stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#>  [1] ggplot2_4.0.3                       sf_1.1-2                           
#>  [3] SingleCellExperiment_1.35.2         MultiAssaySpatialExperiment_0.99.12
#>  [5] MultiAssayExperiment_1.39.1         SummarizedExperiment_1.43.0        
#>  [7] Biobase_2.73.2                      GenomicRanges_1.65.4               
#>  [9] Seqinfo_1.3.2                       IRanges_2.47.5                     
#> [11] S4Vectors_0.51.9                    BiocGenerics_0.59.12               
#> [13] generics_0.1.4                      MatrixGenerics_1.25.0              
#> [15] matrixStats_1.5.0                   BiocStyle_2.41.0                   
#> 
#> loaded via a namespace (and not attached):
#>  [1] gtable_0.3.6             rjson_0.2.23             xfun_0.60               
#>  [4] bslib_0.12.0             lattice_0.23-1           vctrs_0.7.3             
#>  [7] tools_4.6.1              tibble_3.3.1             proxy_0.4-29            
#> [10] pkgconfig_2.0.3          Matrix_1.7-6             KernSmooth_2.23-27      
#> [13] RColorBrewer_1.1-3       S7_0.2.2                 lifecycle_1.0.5         
#> [16] compiler_4.6.1           farver_2.1.2             codetools_0.2-20        
#> [19] htmltools_0.5.9          sys_3.4.3                buildtools_1.0.0        
#> [22] class_7.3-24             sass_0.4.10              yaml_2.3.12             
#> [25] pillar_1.11.1            jquerylib_0.1.4          classInt_0.4-11         
#> [28] DelayedArray_0.39.6      cachem_1.1.0             magick_2.9.1            
#> [31] abind_1.4-8              tidyselect_1.2.1         digest_0.6.39           
#> [34] dplyr_1.2.1              maketools_1.3.2          fastmap_1.2.0           
#> [37] grid_4.6.1               cli_3.6.6                SparseArray_1.13.2      
#> [40] magrittr_2.0.5           S4Arrays_1.13.0          e1071_1.7-17            
#> [43] withr_3.0.3              scales_1.4.0             rmarkdown_2.32          
#> [46] XVector_0.53.0           otel_0.2.0               SpatialExperiment_1.23.0
#> [49] evaluate_1.0.5           knitr_1.51               rlang_1.3.0             
#> [52] Rcpp_1.1.2               glue_1.8.1               DBI_1.3.0               
#> [55] BiocManager_1.30.27      jsonlite_2.0.0           R6_2.6.1                
#> [58] units_1.0-1

MultiAssaySpatialExperiment use cases

Introduction

Part 1: Data import

Reading Xenium data

Basic usage

Options

When to load transcripts

Reading Visium data

Basic usage

Options

Visium HD support

Reading CosMx data

Basic usage

Options

FOV (Field of View) handling

Reading MERSCOPE data (MERFISH)

Basic usage

Options

Choosing a reader

Converting from other Bioconductor classes

From SpatialExperiment

From SpatialFeatureExperiment

Reading multiple samples

Part 2: Real-world workflows

Workflow 1: Xenium + Visium integration

Background

Data preparation

Build MASE object

Exploratory analysis

Do the two assays overlap in space?

Which genes are shared between assays

Spatial annotation: Cell types in Visium spots

Reading the composition result

Spatial aggregation: Summarize Xenium by regions

Comparing expression between the two platforms

Integration insights

Workflow 2: Coordinate transformation and alignment

Background

Converting Visium grid to microns

Applying affine transformations

Landmark-based alignment

Storing transformed coordinates

Checking the alignment visually

Workflow 3: Multi-sample analysis

Workflow 4: Iterative analysis pattern

Summary

See also

Session info