Analysis of Image-Based Spatial Data in Seurat: Vizgen MERFISH Mouse Brain (S2R1)

1 Overview

This document reproduces the Vizgen MERSCOPE/MERFISH section of the Seurat v5 tutorial “Analysis of Image-based Spatial Data in Seurat” (Satija Lab), applied to the Mouse Brain Receptor Map, Slice 2 Replicate 1 (S2R1) dataset.

The dataset contains:

• A targeted 483-gene MERFISH panel (canonical cell-type markers, non-sensory GPCRs, and receptor tyrosine kinases).
• A cell-by-gene count matrix.
• Cell metadata with spatial centroid coordinates.
• Cell segmentation boundaries (HDF5 format, one file per spatial tile).
• Detected-transcript (molecule-level) coordinates.

The workflow performed here is:

1. Validate the local data directory.
2. Load the data with LoadVizgen().
3. Inspect centroids, segmentation boundaries, and molecules.
4. Run quality control.
5. Normalize with SCTransform().
6. Run PCA, UMAP, neighbor-graph construction, and clustering.
7. Visualize clusters and gene expression in UMAP space and in tissue coordinates.
8. Visualize individual detected molecules.
9. Crop a region of interest and visualize cell segmentation boundaries.
10. Identify cluster marker genes.
11. Save the processed object.

The large raw microscopy images/ directory (~23 GB locally, considerably larger on the Vizgen server) is not required for any step shown here and is intentionally excluded.

---

2 load libraries

This dataset was produced using the Vizgen MERSCOPE system, which utilizes the MERFISH technology. The total dataset is available for public download, and contains nine samples (three full coronal slices of the mouse brain, with three biological replicates per slice). The gene panel consists of 483 gene targets, representing known anonical cell type markers, nonsensory G-Protein coupled receptors (GPCRs), and Receptor Tyrosine Kinases (RTKs). In this vignette, we analyze one of the samples - slice 2, replicate 1. The median number of transcripts detected in each cell is 206.

First, we read in the dataset and create a Seurat object.

We use the LoadVizgen() function, which we have written to read in the output of the Vizgen analysis pipeline. The resulting Seurat object contains the following information:

A count matrix, indicating the number of observed molecules for each of the 483 transcripts in each cell. This matrix is analogous to a count matrix in scRNA-seq, and is stored by default in the RNA assay of the Seurat object

3 Loading segmentations

4 Loading segmentations is a slow process and multi processing with the future pacakge is

5 recommended

library(future)

plan("multisession", workers = 4)
options(future.globals.maxSize = 10 * 1024^3)

data_dir <- file.path(getwd(), "data", "vizgen", "s2r1")

vizgen.obj <- LoadVizgen(
  data.dir = data_dir,
  fov = "s2r1"
)

vizgen.obj

An object of class Seurat 
483 features across 83546 samples within 1 assay 
Active assay: Vizgen (483 features, 0 variable features)
 1 layer present: counts
 1 spatial field of view present: s2r1

Cell Centroids: The spatial coordinates marking the centroid for each cell being profiled

# Get the center position of each centroid. There is one row per cell in this dataframe.
head(GetTissueCoordinates(vizgen.obj[["s2r1"]][["centroids"]]))

Cell Segmentation Boundaries: The spatial coordinates that describe the polygon segmentation of each single cell

# Get the coordinates for each segmentation vertice. Each cell will have a variable number of
# vertices describing its shape.
head(GetTissueCoordinates(vizgen.obj[["s2r1"]][["segmentation"]]))

Molecule positions: The spatial coordinates for each individual molecule that was detected during the multiplexed smFISH experiment.

# Fetch molecules positions for Chrm1
head(FetchData(vizgen.obj[["s2r1"]][["molecules"]], vars = "Chrm1"))

6 Preprocessing and unsupervised analysis

We start by performing a standard unsupervised clustering analysis, essentially first treating the dataset as an scRNA-seq experiment. We use SCTransform-based normalization, though we slightly modify the default clipping parameters to mitigate the effect of outliers that we occasionally observe in smFISH experiments. After normalization, we can run dimensional reduction and clustering.

vizgen.obj <- SCTransform(vizgen.obj, assay = "Vizgen", clip.range = c(-10, 10))
vizgen.obj <- RunPCA(vizgen.obj, npcs = 30, features = rownames(vizgen.obj))
vizgen.obj <- RunUMAP(vizgen.obj, dims = 1:30)

Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

vizgen.obj <- FindNeighbors(vizgen.obj, reduction = "pca", dims = 1:30)
vizgen.obj <- FindClusters(vizgen.obj, resolution = 0.3)

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 83546
Number of edges: 3056263

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9595
Number of communities: 26
Elapsed time: 19 seconds

We can then visualize the results of the clustering either in UMAP space (with DimPlot()) or overlaid on the image with ImageDimPlot().

DimPlot(vizgen.obj, reduction = "umap")

ImageDimPlot(vizgen.obj, fov = "s2r1", cols = "polychrome", axes = TRUE)

You can also customize multiple aspect of the plot, including the color scheme, cell border widths, and size (see below).

7 Customizing spatial plots in Seurat

The ImageDimPlot() and ImageFeaturePlot() functions have a few parameters which you can customize individual visualizations. These include:

alpha: Ranges from 0 to 1. Sets the transparency of within-cell coloring. size: determines the size of points representing cells, if centroids are being plotted cols: Sets the color scheme for the internal shading of each cell. Examples settings are polychrome, glasbey, Paired, Set3, and parade. Default is the ggplot2 color palette shuffle.cols: In some cases the selection of cols is more effective when the same colors are assigned to different clusters. Set shuffle.cols = TRUE to randomly shuffle the colors in the palette. border.size: Sets the width of the cell segmentation borders. By default, segmentations are plotted with a border size of 0.3 and centroids are plotted without border. border.color: Sets the color of the cell segmentation borders dark.background: Sets a black background color (TRUE by default) axes: Display Since it can be difficult to visualize the spatial localization patterns of an individual cluster when viewing them all together, we can highlight all cells that belong to a particular cluster:

p1 <- ImageDimPlot(vizgen.obj, fov = "s2r1", cols = "red", cells = WhichCells(vizgen.obj, idents = 1))
p2 <- ImageDimPlot(vizgen.obj, fov = "s2r1", cols = "red", cells = WhichCells(vizgen.obj, idents = 15))
p1 + p2

We can find markers of individual clusters and visualize their spatial expression pattern. We can color cells based on their quantified expression of an individual gene, using ImageFeaturePlot(), which is analagous to the FeaturePlot() function for visualizing expression on a 2D embedding. Since MERFISH images individual molecules, we can also visualize the location of individual molecules.

p1 <- ImageFeaturePlot(vizgen.obj, features = "Slc17a7")
p2 <- ImageDimPlot(vizgen.obj, molecules = "Slc17a7", nmols = 10000, alpha = 0.3, mols.cols = "red")
p1 + p2

Note that the nmols parameter can be used to reduce the total number of molecules shown to reduce overplotting. You can also use the mols.size, mols.cols, and mols.alpha parameter to further optimize.

Plotting molecules is especially useful for visualizing co-expression of multiple genes on the same plot.

p1 <- ImageDimPlot(vizgen.obj, fov = "s2r1", alpha = 0.3, molecules = c("Slc17a7", "Olig1"), nmols = 10000)
markers.14 <- FindMarkers(vizgen.obj, ident.1 = "14")
p2 <- ImageDimPlot(vizgen.obj, fov = "s2r1", alpha = 0.3, molecules = rownames(markers.14)[1:4],
    nmols = 10000)
p1 + p2

The updated Seurat spatial framework has the option to treat cells as individual points, or also to visualize cell boundaries (segmentations). By default, Seurat ignores cell segmentations and treats each cell as a point (‘centroids’). This speeds up plotting, especially when looking at large areas, where cell boundaries are too small to visualize.

We can zoom into a region of tissue, creating a new field of view. For example, we can zoom into a region that contains the hippocampus. Once zoomed-in, we can set DefaultBoundary() to show cell segmentations. You can also ‘simplify’ the cell segmentations, reducing the number of edges in each polygon to speed up plotting.

# create a Crop
cropped.coords <- Crop(vizgen.obj[["s2r1"]], x = c(1750, 3000), y = c(3750, 5250), coords = "plot")
# set a new field of view (fov)
vizgen.obj[["hippo"]] <- cropped.coords

# visualize FOV using default settings (no cell boundaries)
p1 <- ImageDimPlot(vizgen.obj, fov = "hippo", axes = TRUE, size = 0.7, border.color = "white", cols = "polychrome",
    coord.fixed = FALSE)

# visualize FOV with full cell segmentations
DefaultBoundary(vizgen.obj[["hippo"]]) <- "segmentation"
p2 <- ImageDimPlot(vizgen.obj, fov = "hippo", axes = TRUE, border.color = "white", border.size = 0.1,
    cols = "polychrome", coord.fixed = FALSE)

# simplify cell segmentations
vizgen.obj[["hippo"]][["simplified.segmentations"]] <- Simplify(coords = vizgen.obj[["hippo"]][["segmentation"]],
    tol = 3)
DefaultBoundary(vizgen.obj[["hippo"]]) <- "simplified.segmentations"

# visualize FOV with simplified cell segmentations
DefaultBoundary(vizgen.obj[["hippo"]]) <- "simplified.segmentations"
p3 <- ImageDimPlot(vizgen.obj, fov = "hippo", axes = TRUE, border.color = "white", border.size = 0.1,
    cols = "polychrome", coord.fixed = FALSE)

p1 + p2 + p3

#What is the tol parameter? The tol parameter determines how simplified the resulting segmentations are. A higher value of tol will reduce the number of vertices more drastically which will speed up plotting, but some segmentation detail will be lost. See https://rgeos.r-forge.r-project.org/reference/topo-unary-gSimplify.html for examples using different values for tol.

We can visualize individual molecules plotted at higher resolution after zooming-in

# Since there is nothing behind the segmentations, alpha will slightly mute colors
ImageDimPlot(vizgen.obj, fov = "hippo", molecules = rownames(markers.14)[1:4], cols = "polychrome",
    mols.size = 1, alpha = 0.5, mols.cols = c("red", "blue", "yellow", "green"))

8 Mouse Brain: 10x Genomics Xenium In Situ

In this section we’ll analyze data produced by the Xenium platform. The vignette demonstrates how to load the per-transcript location data, cell x gene matrix, cell segmentation, and cell centroid information available in the Xenium outputs. The resulting Seurat object will contain the gene expression profile of each cell, the centroid and boundary of each cell, and the location of each individual detected transcript. The per-cell gene expression profiles are similar to standard single-cell RNA-seq and can be analyzed using the same tools.

This uses the Tiny subset dataset from 10x Genomics provided in the Fresh Frozen Mouse Brain for Xenium Explorer Demo which can be downloaded as described below. These analysis steps are also compatible with the larger Full coronal section, but will take longer to execute.

system("wget https://cf.10xgenomics.com/samples/xenium/1.0.2/Xenium_V1_FF_Mouse_Brain_Coronal_Subset_CTX_HP/Xenium_V1_FF_Mouse_Brain_Coronal_Subset_CTX_HP_outs.zip")

url <- "https://cf.10xgenomics.com/samples/xenium/1.0.2/Xenium_V1_FF_Mouse_Brain_Coronal_Subset_CTX_HP/Xenium_V1_FF_Mouse_Brain_Coronal_Subset_CTX_HP_outs.zip"

download.file(
  url,
  destfile = "Xenium_V1_FF_Mouse_Brain_Coronal_Subset_CTX_HP_outs.zip",
  mode = "wb"
)

First we read in the dataset and create a Seurat object. Provide the path to the data folder for a Xenium run as the input path. The RNA data is stored in the Xenium assay of the Seurat object. This step should take about a minute (you can improve this by installing arrow and hdf5r).

path <- "~/1-Thesis_Final_Year_2025/2025-Year3_Analysis/1-scRNA_RESULTS-19-11-2025/Spatial_Analysis/Image_Based/data/vizgen/s2r1"
# Load the Xenium data
xenium.obj <- LoadXenium(path, fov = "fov", segmentations = "cell")
# remove cells with 0 counts
xenium.obj <- subset(xenium.obj, subset = nCount_Xenium > 0)

Spatial information is loaded into slots of the Seurat object, labelled by the name of “field of view” (FOV) being loaded. Initially all the data is loaded into the FOV named fov. Later, we will make a cropped FOV that zooms into a region of interest.

Standard QC plots provided by Seurat are available via the Xenium assay. Here are violin plots of genes per cell (nFeature_Xenium) and transcript counts per cell (nCount_Xenium)

VlnPlot(xenium.obj, features = c("nFeature_Xenium", "nCount_Xenium"), ncol = 2, pt.size = 0)

Next, we plot the positions of the pan-inhibitory neuron marker Gad1, inhibitory neuron sub-type markers Pvalb, and Sst, and astrocyte marker Gfap on the tissue using ImageDimPlot().

ImageDimPlot(xenium.obj, fov = "fov", molecules = c("Gad1", "Sst", "Pvalb", "Gfap"), nmols = 20000)

Here we visualize the expression level of some key layer marker genes at the per-cell level using ImageFeaturePlot() which is analogous to the FeaturePlot() function for visualizing expression on a 2D embedding. We manually adjust the max.cutoff for each gene to roughly the 90th percentile (which can be specified with max.cutoff=‘q90’) of it’s count distribution to improve contrast.

ImageFeaturePlot(xenium.obj, features = c("Cux2", "Rorb", "Bcl11b", "Foxp2"), max.cutoff = c(25,
    35, 12, 10), size = 0.75, cols = c("white", "red"))

We can zoom in on a chosen area with the Crop() function. Once zoomed-in, we can visualize cell segmentation boundaries along with individual molecules.

cropped.coords <- Crop(xenium.obj[["fov"]], x = c(1200, 2900), y = c(3750, 4550), coords = "plot")
xenium.obj[["zoom"]] <- cropped.coords
# visualize cropped area with cell segmentations & selected molecules
DefaultBoundary(xenium.obj[["zoom"]]) <- "segmentation"
ImageDimPlot(xenium.obj, fov = "zoom", axes = TRUE, border.color = "white", border.size = 0.1, cols = "polychrome",
    coord.fixed = FALSE, molecules = c("Gad1", "Sst", "Npy2r", "Pvalb", "Nrn1"), nmols = 10000)

Next, we use SCTransform for normalization followed by standard dimensionality reduction and clustering. This step takes about 5 minutes from start to finish.

xenium.obj <- SCTransform(xenium.obj, assay = "Xenium")
xenium.obj <- RunPCA(xenium.obj, npcs = 30, features = rownames(xenium.obj))
xenium.obj <- RunUMAP(xenium.obj, dims = 1:30)

Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

xenium.obj <- FindNeighbors(xenium.obj, reduction = "pca", dims = 1:30)
xenium.obj <- FindClusters(xenium.obj, resolution = 0.3)

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 36553
Number of edges: 1340890

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9587
Number of communities: 28
Elapsed time: 4 seconds

We can then visualize the results of the clustering by coloring each cell according to its cluster either in UMAP space with DimPlot() or overlaid on the image with ImageDimPlot().

DimPlot(xenium.obj)

We can visualize the expression level of the markers we looked at earlier on the UMAP coordinates.

FeaturePlot(xenium.obj, features = c("Cux2", "Bcl11b", "Foxp2", "Gad1", "Sst", "Gfap"))

We can now use ImageDimPlot() to color the cell positions colored by the cluster labels determined in the previous step.

ImageDimPlot(xenium.obj, cols = "polychrome", size = 0.75)

Using the positional information of each cell, we compute spatial niches. We use a cortex reference from the the Allen Brain Institute to annotate cells, so we first crop the dataset to the cortex. The Allen Brain reference can be installed here.

Below, we use Slc17a7 expression to help determine the cortical region.

#To be continued………………..

.............................