Introduction

This vignette covers the core operations for working with MultiAssaySpatialExperiment (MASE) objects:

Construction: Building MASE objects from matrices, existing Bioconductor objects, or vendor output (see MultiAssaySpatialExperiment use cases)
Subsetting: Filtering by specimens, assays, columns and spatial regions
Spatial operations: Annotation and aggregation by spatial regions
Labels ↔︎ shapes: Converting between raster masks and vector polygons

Prerequisite: If you are new to MASE, read the Introduction to MultiAssaySpatialExperiment vignette first to understand the MASE structure and components.

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

Building MASE objects

There are several approaches to constructing a MASE object:

Minimal MASE: Start with assays + metadata, add spatial elements later
From scratch: Build all components manually
From existing objects: Convert SpatialExperiment or wrap SingleCellExperiment objects in an ExperimentList
From vendor output: Use platform readers (readXeniumMASE(), etc.), covered in the MultiAssaySpatialExperiment use cases vignette
With prepMASE: Wrap prepMultiAssay() for name harmonization and spatial FK cleanup (see below)

Throughout, the sampleMap is the three-column DataFrame that says which specimen each assay column belongs to, with one row per column: assay names the element of the ExperimentList, primary is a row name of colData, and colname is the column name within that assay. The Introduction to MultiAssaySpatialExperiment vignette covers sampleMap and spatialMap in more detail.

Optional: lazy Parquet I/O via BiocDuckDB

MultiAssaySpatialExperiment does not import BiocDuckDB. For large on-disk datasets, the BiocDuckDB package registers methods on MASE generics (readParquetForMASE, lazy spatialPoints / spatialShapes, etc.) so you can persist and query MASE objects without loading full layers into memory. See the BiocDuckDB documentation for Parquet read/write workflows.

Minimal MASE

The simplest MASE has one experiment, specimen metadata (colData), and a sampleMap. No spatial elements are required initially. In the example below, each cell is treated as its own specimen for simplicity; in practice, colData rows usually represent biological replicates or tissue sections.

# Create a simple expression matrix
counts <- matrix(rpois(50, lambda = 10), nrow = 10, ncol = 5,
                dimnames = list(paste0("Gene", 1:10),
                               paste0("Cell", 1:5)))

# Specimen metadata
specimens <- DataFrame(
  patient_id = c("P1", "P1", "P1", "P2", "P2"),
  tissue = c("cortex", "cortex", "cortex", "medulla", "medulla"),
  row.names = paste0("Cell", 1:5)
)

# Sample map: link assay columns to specimens
sample_map <- DataFrame(
  assay = factor("rna", "rna"),
  primary = paste0("Cell", 1:5),
  colname = paste0("Cell", 1:5)
)

# Construct MASE
mase_minimal <- MultiAssaySpatialExperiment(
  experiments = ExperimentList(rna = counts),
  colData = specimens,
  sampleMap = sample_map
)

mase_minimal
#> A MultiAssaySpatialExperiment object of 1 listed
#>  experiment with a user-defined name and respective class.
#>  Containing an ExperimentList class object of length 1:
#>  [1] rna: matrix with 10 rows and 5 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: 0 elements
#>   spatialShapes: 0 elements
#>   imgData: NULL
#>   spatialMap: NULL

Adding spatial coordinates

To add point coordinates (e.g., cell centroids), create a PointsLayerList and a spatialMap linking assay columns to spatial instances.

# Suppose we have two assays: RNA and protein
rna_counts <- matrix(rpois(50, 10), nrow = 10, ncol = 5,
                    dimnames = list(paste0("Gene", 1:10), paste0("Cell", 1:5)))
prot_counts <- matrix(rpois(15, 5), nrow = 3, ncol = 5,
                     dimnames = list(paste0("Protein", 1:3), paste0("Cell", 1:5)))

# Both assays share the same specimens
specimens2 <- DataFrame(
  patient_id = c("P1", "P1", "P1", "P2", "P2"),
  row.names = paste0("Cell", 1:5)
)

# Sample map for two assays
sample_map2 <- DataFrame(
  assay = factor(rep(c("rna", "protein"), each = 5), c("rna", "protein")),
  primary = rep(paste0("Cell", 1:5), 2),
  colname = rep(paste0("Cell", 1:5), 2)
)

# Point coordinates (shared by both assays)
coords2 <- DataFrame(
  x = runif(5, 0, 100),
  y = runif(5, 0, 100),
  instance_id = paste0("Cell", 1:5)
)

# Spatial map: link assay columns to spatial points
spatial_map2 <- DataFrame(
  assay = factor(rep(c("rna", "protein"), each = 5), c("rna", "protein")),
  colname = rep(paste0("Cell", 1:5), 2),
  element_type = "points",
  region = factor(rep("coords", 10), "coords"),
  instance_id = rep(paste0("Cell", 1:5), 2)
)

# Construct MASE with spatialMap
mase_with_map <- MultiAssaySpatialExperiment(
  experiments = ExperimentList(rna = rna_counts, protein = prot_counts),
  colData = specimens2,
  sampleMap = sample_map2,
  points = PointsLayerList(coords = coords2),
  spatialMap = spatial_map2
)

mase_with_map
#> A MultiAssaySpatialExperiment object of 2 listed
#>  experiments with user-defined names and respective classes.
#>  Containing an ExperimentList class object of length 2:
#>  [1] rna: matrix with 10 rows and 5 columns
#>  [2] protein: matrix with 3 rows and 5 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: 1 element
#>   spatialShapes: 0 elements
#>   imgData: NULL
#>   spatialMap: present

Adding spatial shapes

Shapes (polygons) typically represent cell boundaries, tissue regions, or annotation areas. Use the sf package to create geometries.

# Create simple circular boundaries around each point
centroids <- data.frame(
  x = runif(5, 0, 100),
  y = runif(5, 0, 100),
  instance_id = paste0("Cell", 1:5)
)

# Convert to sf points, then buffer to circles
centroids_sf <- st_as_sf(centroids, coords = c("x", "y"))
circles <- st_buffer(centroids_sf, dist = 5)

# Convert back to DataFrame with geometry column
boundaries <- DataFrame(
  geometry = st_geometry(circles),
  instance_id = paste0("Cell", 1:5)
)

# Add to MASE
mase_with_shapes <- mase_with_map
spatialShapes(mase_with_shapes) <- ShapesLayerList(boundaries = boundaries)

mase_with_shapes
#> A MultiAssaySpatialExperiment object of 2 listed
#>  experiments with user-defined names and respective classes.
#>  Containing an ExperimentList class object of length 2:
#>  [1] rna: matrix with 10 rows and 5 columns
#>  [2] protein: matrix with 3 rows and 5 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: 1 element
#>   spatialShapes: 1 element
#>   imgData: NULL
#>   spatialMap: present

Building from multiple SingleCellExperiment objects

Combine multiple SingleCellExperiment objects into a MASE:

# Create two SingleCellExperiment objects (simulating different technologies)
sce1 <- SingleCellExperiment(
  assays = list(counts = matrix(rpois(40, 10), nrow = 10, ncol = 4,
                               dimnames = list(paste0("Gene", 1:10),
                                              paste0("Cell", 1:4)))),
  colData = DataFrame(tech = "rna", cell_id = paste0("Cell", 1:4))
)

sce2 <- SingleCellExperiment(
  assays = list(counts = matrix(rpois(12, 5), nrow = 3, ncol = 4,
                               dimnames = list(paste0("Protein", 1:3),
                                              paste0("Cell", 1:4)))),
  colData = DataFrame(tech = "protein", cell_id = paste0("Cell", 1:4))
)

# Specimen metadata (shared)
specimens_sce <- DataFrame(
  patient = rep("P1", 4),
  tissue = rep("cortex", 4),
  row.names = paste0("Cell", 1:4)
)

# Sample maps
samplemap_sce <- DataFrame(
  assay = factor(rep(c("rna", "protein"), each = 4), c("rna", "protein")),
  primary = rep(paste0("Cell", 1:4), 2),
  colname = rep(paste0("Cell", 1:4), 2)
)

# Construct MASE
mase_from_sce <- MultiAssaySpatialExperiment(
  experiments = ExperimentList(rna = sce1, protein = sce2),
  colData = specimens_sce,
  sampleMap = samplemap_sce
)

mase_from_sce
#> A MultiAssaySpatialExperiment object of 2 listed
#>  experiments with user-defined names and respective classes.
#>  Containing an ExperimentList class object of length 2:
#>  [1] rna: SingleCellExperiment with 10 rows and 4 columns
#>  [2] protein: SingleCellExperiment with 3 rows and 4 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: 0 elements
#>   spatialShapes: 0 elements
#>   imgData: NULL
#>   spatialMap: NULL

prepMASE and buildSpatialMap

For manual construction, buildSpatialMap() assembles a spatialMap from sampleMap (analogous to MAE’s listToMap()), and prepMASE() wraps prepMultiAssay() while harmonizing spatial slots:

spmap <- buildSpatialMap(sample_map, region = "cells", element_type = "shapes")
prepared <- prepMASE(
  ExperimentList(rna = counts, protein = protein_counts),
  colData = specimens,
  sample_map,
  shapes = ShapesLayerList(cells = cell_boundaries),
  spatialMap = spmap
)
mase <- do.call(MultiAssaySpatialExperiment, prepared)
drops(prepared$metadata$drops)  # rows removed during harmonization, if any

Validation checklist

Before constructing a MASE, ensure:

ExperimentList:
- All experiments have unique column names
colData:
- Row names are specimen identifiers (primary IDs)
- No duplicate row names
sampleMap:
- assay column matches names in ExperimentList
- primary column matches row names in colData
- colname column matches column names in experiments
- No duplicate (assay, colname) pairs
points/shapes (if provided):
- instance_id column exists
- For points: x, y columns exist
- For shapes: geometry column exists (sfc)
spatialMap (if provided):
- assay + colname exist in sampleMap
- element_type is “points” or “shapes”
- region matches names in points or shapes
- instance_id exists in the corresponding spatial element

Check validity manually:

validObject(mase_with_map)
#> [1] TRUE

Most of these are checked at construction time, not just by validObject(). sampleMap rows that fail the colname/primary checks are dropped with a message rather than an error, since a partial mismatch is common when building a sampleMap by hand:

bad_sample_map <- sample_map2
bad_sample_map$colname[1] <- "NotARealColumn"
mase_dropped <- MultiAssaySpatialExperiment(
  experiments = ExperimentList(rna = rna_counts, protein = prot_counts),
  colData = specimens2,
  sampleMap = bad_sample_map
)
#> harmonizing input:
#>   removing 1 sampleMap rows with 'colname' not in colnames of experiments

A broken spatialMap, by contrast, fails validObject() outright, since there is no reasonable way to silently drop a dangling geometry reference:

bad_spatial_map <- spatial_map2
bad_spatial_map$instance_id[1] <- "NotACell"
MultiAssaySpatialExperiment(
  experiments = ExperimentList(rna = rna_counts, protein = prot_counts),
  colData = specimens2,
  sampleMap = sample_map2,
  points = PointsLayerList(coords = coords2),
  spatialMap = bad_spatial_map
)
#> Error in `validObject()`:
#> ! invalid class "MultiAssaySpatialExperiment" object: 
#>     spatialMap: instance_id value(s) not found in points$coords: NotACell

Subsetting operations

MASE extends MultiAssayExperiment with spatial slots (points, shapes, images, labels, imgData, and spatialMap). When you subset a MASE, these slots are automatically updated to stay consistent with the assays and specimen metadata.

The examples below use a separate MASE object from the construction section above, so you can run them independently.

Overview of subsetting operations

Operation	What it filters	Spatial propagation
`subsetByColData(x, y)`	Specimens (primary identifiers)	`spatialMap`, `imgData`
`subsetByRow(x, y, ...)`	Rows of experiments	Spatial layers unchanged
`subsetByColumn(x, y)`	Assay columns (names, indices, or logical vectors per assay)	`spatialMap`, `imgData`, points, shapes, images, labels
`subsetByAssay(x, y)`	Assays	Points, shapes, images, labels and `spatialMap` rows tied to retained assays
`x[i, j, k, ..., drop]`	Row, column, assay indices	Composes the above
`subsetByBoundingBox(x, xmin, xmax, ymin, ymax, ...)`	Points/shapes by rectangle	Assays via `spatialMap`
`subsetByPolygon(x, polygon, ...)`	Points/shapes by polygon	Assays via `spatialMap`

Note: subsetByAssay() keeps only spatial layers referenced in spatialMap for the retained assays. Auxiliary shape layers added for annotation (but not linked in spatialMap) are dropped; run annotateWithRegions() on the full MASE first if you need those layers.

Basic subsetting

Two assays are enough to show what the different subsetting axes actually do. Here P1 and P2 are the specimens, the rows of colData; S1 to S4 are observations, the assay columns. Both assays measure the same four cells, and both are mapped onto one shared points layer.

make_assay <- function() {
  SummarizedExperiment(
    assays = list(counts = matrix(rnorm(20), nrow = 5, ncol = 4,
      dimnames = list(paste0("G", 1:5), paste0("S", 1:4)))),
    # "zone", not "region": this is an arbitrary per-column covariate on the
    # assay's own colData, unrelated to spatialMap's "region" column below,
    # which names a points/shapes layer
    colData = DataFrame(zone = rep(c("core", "margin"), length.out = 4)))
}

pts <- DataFrame(x = 1:4, y = 1:4, instance_id = paste0("S", 1:4))

mase <- MultiAssaySpatialExperiment(
  experiments = ExperimentList(rna = make_assay(), protein = make_assay()),
  colData = DataFrame(row.names = c("P1", "P2")),
  sampleMap = DataFrame(
    assay = factor(rep(c("rna", "protein"), each = 4)),
    primary = rep(c("P1", "P1", "P2", "P2"), 2),
    colname = rep(paste0("S", 1:4), 2)
  ),
  points = PointsLayerList(coords = pts),
  spatialMap = DataFrame(
    assay = factor(rep(c("rna", "protein"), each = 4)),
    colname = rep(paste0("S", 1:4), 2),
    element_type = "points",
    region = "coords",
    instance_id = rep(paste0("S", 1:4), 2)
  )
)

vapply(experiments(mase), ncol, integer(1))
#>     rna protein 
#>       4       4

Subset by specimen

Naming specimens cuts across the whole object. P1 owns S1 and S2 in both assays, so both come back with two columns:

by_specimen <- mase[, "P1"]
vapply(experiments(by_specimen), ncol, integer(1))
#>     rna protein 
#>       2       2

Subset by assay column

Passing a list instead targets one assay by name. Only rna is cut here; protein keeps all four columns, which is the difference from the previous example:

by_column <- mase[, list(rna = c("S1", "S3"))]
#> harmonizing input:
#>   removing 2 sampleMap rows with 'colname' not in colnames of experiments
vapply(experiments(by_column), ncol, integer(1))
#>     rna protein 
#>       2       4

The same list interface filters on observation metadata, which is the idiomatic replacement for ad hoc “query by metadata” patterns:

cdf <- colData(experiments(mase)[["rna"]])
mase_core <- subsetByColumn(mase, list(rna = cdf$zone == "core"))
#> harmonizing input:
#>   removing 2 sampleMap rows with 'colname' not in colnames of experiments
ncol(experiments(mase_core)[["rna"]])
#> [1] 2
nrow(spatialPoints(mase_core)[["coords"]])
#> [1] 4

Bracket notation

The [i, j, k, drop] operator composes feature, specimen and assay subsetting in one call:

mase[1:3, "P1", "rna", drop = FALSE]
#> Warning: 'experiments' dropped; see 'drops()'
#> harmonizing input:
#>   removing 2 sampleMap rows not in names(experiments)
#> A MultiAssaySpatialExperiment object of 1 listed
#>  experiment with a user-defined name and respective class.
#>  Containing an ExperimentList class object of length 1:
#>  [1] rna: SummarizedExperiment with 3 rows and 2 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: 1 element
#>   spatialShapes: 0 elements
#>   imgData: NULL
#>   spatialMap: present

Spatial subsetting

When you have spatialMap linking assay columns to spatial elements, you can subset by bounding box or polygon. Points and shapes within the region are kept; assays are filtered to retained instance_ids via spatialMap.

Subset by bounding box

Filter by a rectangle (xmin, xmax, ymin, ymax):

m2 <- subsetByBoundingBox(mase, xmin = 1.5, xmax = 4.5, ymin = 1.5, ymax = 4.5)
#> harmonizing input:
#>   removing 2 sampleMap rows with 'colname' not in colnames of experiments
m2
#> A MultiAssaySpatialExperiment object of 2 listed
#>  experiments with user-defined names and respective classes.
#>  Containing an ExperimentList class object of length 2:
#>  [1] rna: SummarizedExperiment with 5 rows and 3 columns
#>  [2] protein: SummarizedExperiment with 5 rows and 3 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: 1 element
#>   spatialShapes: 0 elements
#>   imgData: NULL
#>   spatialMap: present
spatialPoints(m2)[["coords"]]
#> DataFrame with 3 rows and 3 columns
#>           x         y instance_id
#>   <integer> <integer> <character>
#> 1         2         2          S2
#> 2         3         3          S3
#> 3         4         4          S4
ncol(experiments(m2)[["assay1"]])
#> NULL

Subset by polygon

Filter by an sf polygon:

poly <- st_polygon(list(matrix(
  c(1.5, 1.5, 4.5, 1.5, 4.5, 4.5, 1.5, 4.5, 1.5, 1.5),
  ncol = 2, byrow = TRUE)))
m3 <- subsetByPolygon(mase, poly)
#> harmonizing input:
#>   removing 2 sampleMap rows with 'colname' not in colnames of experiments
nrow(spatialPoints(m3)[["coords"]])
#> [1] 3
ncol(experiments(m3)[["assay1"]])
#> NULL

Select data from multiple non-contiguous regions

A region of interest need not be a single connected shape. Subsetting once per region and keeping the results separate lets you compare them directly:

# Define two regions of interest
region1 <- st_polygon(list(matrix(c(1, 1, 2, 1, 2, 2, 1, 2, 1, 1), ncol = 2, byrow = TRUE)))
region2 <- st_polygon(list(matrix(c(4, 4, 5, 4, 5, 5, 4, 5, 4, 4), ncol = 2, byrow = TRUE)))

# Subset by each region
m_r1 <- subsetByPolygon(mase, region1)
#> harmonizing input:
#>   removing 4 sampleMap rows with 'colname' not in colnames of experiments
#>   removing 1 colData rownames not in sampleMap 'primary'
m_r2 <- subsetByPolygon(mase, region2)
#> harmonizing input:
#>   removing 6 sampleMap rows with 'colname' not in colnames of experiments
#>   removing 1 colData rownames not in sampleMap 'primary'

c(region1 = ncol(experiments(m_r1)[[1]]),
  region2 = ncol(experiments(m_r2)[[1]]))
#> region1 region2 
#>       2       1

Buffer-based subsets

Find all points within a distance from a reference point:

# Reference point
ref_point <- st_point(c(3, 3))

# Create buffer (e.g., 1.5 units radius)
buffer_region <- st_buffer(ref_point, dist = 1.5)

# Subset by buffer
m_buffer <- subsetByPolygon(mase, buffer_region)
#> harmonizing input:
#>   removing 2 sampleMap rows with 'colname' not in colnames of experiments
m_buffer
#> A MultiAssaySpatialExperiment object of 2 listed
#>  experiments with user-defined names and respective classes.
#>  Containing an ExperimentList class object of length 2:
#>  [1] rna: SummarizedExperiment with 5 rows and 3 columns
#>  [2] protein: SummarizedExperiment with 5 rows and 3 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: 1 element
#>   spatialShapes: 0 elements
#>   imgData: NULL
#>   spatialMap: present

Spatial annotation and aggregation

Spatial omics often measures signals at discrete locations (spots, cell centroids) that lie within larger regions (cells, tissue segments). A common workflow is to:

Annotate each measurement point with the region it falls in
Aggregate assay values by region (sum, mean, count)

MASE provides annotateWithRegions for step 1 and aggregateByRegion for step 2. Both use the central spatialMap table to link assay columns to spatial elements.

Example setup

We construct a minimal MASE with four spots (transcriptomic measurements), coordinates, two cell polygons, and a gene expression assay. The goal is to annotate each spot with its containing cell and then aggregate expression per cell.

# Four spots at coordinates
pts_agg <- DataFrame(
  x = c(1.5, 2.5, 2.5, 3.5),
  y = c(1.5, 1.5, 2.5, 2.5),
  instance_id = paste0("S", 1:4))

# Three cell polygons
cell1 <- st_polygon(list(matrix(c(1, 1, 2, 1, 2, 2, 1, 2, 1, 1), ncol = 2, byrow = TRUE)))
cell2 <- st_polygon(list(matrix(c(2, 1, 3, 1, 3, 2, 2, 2, 2, 1), ncol = 2, byrow = TRUE)))
cell3 <- st_polygon(list(matrix(c(2, 2, 3, 2, 3, 3, 2, 3, 2, 2), ncol = 2, byrow = TRUE)))

shp_df <- DataFrame(
  instance_id = c("cell1", "cell2", "cell3"),
  geometry = st_sfc(cell1, cell2, cell3))

# Gene expression assay (3 genes × 4 spots)
expr <- matrix(c(10, 20, 5,  15,
                 30, 10, 25, 5,
                 5,  15, 20, 30),
  nrow = 3, ncol = 4,
  dimnames = list(paste0("Gene", 1:3), paste0("S", 1:4)))

# Construct MASE
mase_agg <- MultiAssaySpatialExperiment(
  experiments = ExperimentList(rna = expr),
  colData = DataFrame(row.names = "P1"),
  sampleMap = DataFrame(
    assay = factor("rna"),
    primary = rep("P1", 4),
    colname = paste0("S", 1:4)
  ),
  points = PointsLayerList(centroids = pts_agg),
  shapes = ShapesLayerList(cells = shp_df),
  spatialMap = DataFrame(
    assay = factor("rna"),
    colname = paste0("S", 1:4),
    element_type = "points",
    region = "centroids",
    instance_id = paste0("S", 1:4)
  )
)
mase_agg
#> A MultiAssaySpatialExperiment object of 1 listed
#>  experiment with a user-defined name and respective class.
#>  Containing an ExperimentList class object of length 1:
#>  [1] rna: matrix with 3 rows and 4 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: 1 element
#>   spatialShapes: 1 element
#>   imgData: NULL
#>   spatialMap: present

Annotate with regions

annotateWithRegions performs a spatial join: for each point, it finds which shape (if any) matches and records the shape’s instance_id in a new column in spatialMap.

mase_agg_orig <- mase_agg
mase_agg <- annotateWithRegions(mase_agg, points = "centroids", shapes = "cells")
spatialMap(mase_agg)
#> DataFrame with 4 rows and 6 columns
#>      assay     colname element_type      region instance_id       cells
#>   <factor> <character>  <character> <character> <character> <character>
#> 1      rna          S1       points   centroids          S1       cell1
#> 2      rna          S2       points   centroids          S2       cell2
#> 3      rna          S3       points   centroids          S3       cell3
#> 4      rna          S4       points   centroids          S4          NA

The new column cells shows: S1 → cell1, S2 → cell2, S3 → cell3, S4 → NA (no containing cell). To assign S4 to its nearest cell instead, use join = sf::st_nearest_feature.

Aggregate by region

With the annotation in place, aggregateByRegion groups assay columns by region and applies an aggregation function.

Count points per region

FUN = "count" returns the number of points in each shape:

agg_count <- aggregateByRegion(mase_agg, by = "cells", FUN = "count")
agg_count
#> DataFrame with 3 rows and 2 columns
#>         cells     count
#>   <character> <integer>
#> 1       cell1         1
#> 2       cell2         1
#> 3       cell3         1

Sum expression per region

FUN = "sum" sums each gene’s expression over spots in the same cell:

agg_sum <- aggregateByRegion(mase_agg, by = "cells", FUN = "sum")
agg_sum[["rna"]]
#>       cell1 cell2 cell3
#> Gene1    10    15    25
#> Gene2    20    30     5
#> Gene3     5    10     5

Mean expression per region

FUN = "mean" averages expression across spots in each cell:

agg_mean <- aggregateByRegion(mase_agg, by = "cells", FUN = "mean")
agg_mean[["rna"]]
#>       cell1 cell2 cell3
#> Gene1    10    15    25
#> Gene2    20    30     5
#> Gene3     5    10     5

Lower-level spatial joins

annotateWithRegions() wraps a point-in-polygon join and writes results to spatialMap. For ad hoc layer-to-layer joins, use spatialJoin() on DataFrame layers with geometry columns.

Aggregation strategies

Different strategies serve different purposes: sum preserves total signal (useful for RNA-seq or ATAC-seq counts), mean allows comparing regions with different numbers of points (a normalized comparison), and count supports quality control, density analysis, or use as a normalization factor.

# Compare sum vs mean for Gene1
cells <- colnames(agg_sum[["rna"]])
data.frame(
  cell = cells,
  sum = agg_sum[["rna"]]["Gene1", ],
  mean = agg_mean[["rna"]]["Gene1", ],
  n_spots = agg_count$count[match(cells, agg_count$cells)]
)
#>        cell sum mean n_spots
#> cell1 cell1  10   10       1
#> cell2 cell2  15   15       1
#> cell3 cell3  25   25       1

Visualizing aggregated data

After aggregation, visualize results spatially:

# Get aggregated data for Gene1
gene1_expr <- agg_sum[["rna"]]["Gene1", ]

# Get cell polygons as sf
cells_sf <- st_sf(
  instance_id = shp_df$instance_id,
  geometry = shp_df$geometry
)

# Add expression values
cells_sf$Gene1_sum <- gene1_expr[cells_sf$instance_id]

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

# Original points
plot(pts_agg$x, pts_agg$y, 
     pch = 16, cex = 2,
     col = rainbow(4)[rank(expr["Gene1", ])],
     main = "Gene1: Original spots",
     xlab = "x", ylab = "y")
text(pts_agg$x, pts_agg$y, paste0("S", 1:4), pos = 3, cex = 0.7)

# Aggregated cells
plot(st_geometry(cells_sf), 
     col = rainbow(3)[rank(cells_sf$Gene1_sum)],
     main = "Gene1: Aggregated by cell",
     xlab = "x", ylab = "y")
text(st_coordinates(st_centroid(cells_sf)), 
     labels = round(cells_sf$Gene1_sum, 1),
     cex = 0.8)
#> Warning: st_centroid assumes attributes are constant over geometries

Labels ↔︎ shapes interoperability

MASE supports both labels (raster segmentation masks) and shapes (vector polygons). The package stores both slot types; converting between them, rasterization from shapes to labels, and vectorization from labels back to shapes, is a manual workflow using external tools, covered in the next two sections.

Rasterization: shapes → labels

Convert cell boundary polygons to a segmentation mask:

# Example: 3 cell polygons
cell1_rast <- st_polygon(list(matrix(c(0, 0, 1, 0, 1, 1, 0, 1, 0, 0), ncol = 2, byrow = TRUE)))
cell2_rast <- st_polygon(list(matrix(c(1, 0, 2, 0, 2, 1, 1, 1, 1, 0), ncol = 2, byrow = TRUE)))
cell3_rast <- st_polygon(list(matrix(c(0, 1, 1, 1, 1, 2, 0, 2, 0, 1), ncol = 2, byrow = TRUE)))

shapes_df <- DataFrame(
  instance_id = c("cell1", "cell2", "cell3"),
  geometry = st_sfc(cell1_rast, cell2_rast, cell3_rast)
)

MASE stores the result of rasterization rather than performing it. The steps are to choose a target grid extent and resolution, determine which pixels each polygon covers, write the polygon identifier into those pixels, and store the resulting raster as a RasterLayerList in the labels slot. stars::st_rasterize() and terra::rasterize() both do the middle two steps.

Vectorization: labels → shapes

Extract polygon boundaries from a segmentation mask:

# Example label matrix (3 cells)
label_matrix <- matrix(c(
  1, 1, 2, 2,
  1, 1, 2, 2,
  3, 3, 0, 0,
  3, 3, 0, 0
), nrow = 4, byrow = TRUE)

Vectorization runs the same path in reverse: read the label matrix, trace the boundary of each distinct non-zero value, convert those boundaries to sf polygons carrying the label as instance_id, and store them as a ShapesLayerList in the shapes slot. stars::st_as_sf() and sf::st_polygonize() cover the tracing step.

Use cases

Rasterization (shapes → labels):

Generate training data for deep learning models (requires raster masks)
Perform fast pixel-based operations (convolution, morphology)
Memory-efficient storage for dense segmentations

Vectorization (labels → shapes):

Perform spatial operations (area, perimeter, centroid)
Spatial subsetting (point-in-polygon, overlap)
Visualization with ggplot2 + sf

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] BiocBaseUtils_1.15.1     pkgconfig_2.0.3          Matrix_1.7-6            
#> [13] KernSmooth_2.23-27       RColorBrewer_1.1-3       S7_0.2.2                
#> [16] lifecycle_1.0.5          compiler_4.6.1           farver_2.1.2            
#> [19] codetools_0.2-20         htmltools_0.5.9          sys_3.4.3               
#> [22] buildtools_1.0.0         class_7.3-24             sass_0.4.10             
#> [25] yaml_2.3.12              pillar_1.11.1            jquerylib_0.1.4         
#> [28] classInt_0.4-11          DelayedArray_0.39.6      cachem_1.1.0            
#> [31] magick_2.9.1             abind_1.4-8              tidyselect_1.2.1        
#> [34] digest_0.6.39            dplyr_1.2.1              maketools_1.3.2         
#> [37] fastmap_1.2.0            grid_4.6.1               cli_3.6.6               
#> [40] SparseArray_1.13.2       magrittr_2.0.5           S4Arrays_1.13.0         
#> [43] e1071_1.7-17             withr_3.0.3              scales_1.4.0            
#> [46] rmarkdown_2.32           XVector_0.53.0           otel_0.2.0              
#> [49] SpatialExperiment_1.23.0 evaluate_1.0.5           knitr_1.51              
#> [52] rlang_1.3.0              Rcpp_1.1.2               glue_1.8.1              
#> [55] DBI_1.3.0                BiocManager_1.30.27      jsonlite_2.0.0          
#> [58] R6_2.6.1                 units_1.0-1

Working with MultiAssaySpatialExperiment