Third Party Tools

This page walks through reading the per-sample outputs from the results-folder into the two most common single-cell / spatial transcriptomics ecosystems:

• Scanpy (Python, AnnData)
• Seurat (R, Seurat object)

Both walkthroughs follow the same structure:

1. Read the Matrix Market expression matrix (matrix.mtx.gz + features.tsv.gz + barcodes.tsv.gz).
2. Parse spatial coordinates (X, Y) out of the barcode strings.
3. Attach the registered OME-TIFF image and align coordinates to the image's local space.
4. (Optional) Load cell / nuclei contour CSVs.

Resource requirements: DRAGEN StrataMap outputs are large. A single sample's raw (per-barcode) MTX can contain tens of millions of barcodes and hundreds of millions of non-zero counts, resulting in files that are several gigabytes on disk and expand to 20–60 GB in memory when loaded. Cell-binned and grid-binned matrices are substantially smaller, but still routinely require 8–32 GB RAM depending on tissue area and cell density. Plan for at least 32 GB RAM if you intend to load raw data, run dimensionality reduction, or hold multiple samples in memory simultaneously.

Coordinate note: Transcripts in barcodes.tsv.gz are in global substrate coordinates (µm). Registered OME-TIFF images and TIFF masks are in local sample coordinates (µm). To overlay points on the image you must subtract Global_left / Global_top (read from the OME-TIFF metadata) from the spatial coordinates. See Global and Local Coordinates.

The barcode format is <type>:Y:X, where <type> is one of:

• SBC — raw / un-binned (per spatial barcode), e.g. SBC:3084.702:29810.565
• bin<N> — grid-binned at N µm, e.g. bin10:10005:27955
• cell<ID> — cell- or nuclei-binned, e.g. cell10119:5238:29404

Throughout the examples below, replace SAMPLEID with your actual sample name. The sample name is the same as the folder name inside results/.

Per-sample inputs used in the examples:

Variable	Path
mtx_path (cell-bin)	results/SAMPLEID/SAMPLEID_cell_binned/
mtx_path (grid)	results/SAMPLEID/SAMPLEID_grid_binned_10um/
mtx_path (raw)	results/SAMPLEID/SAMPLEID_raw/
ome_tiff_path	results/SAMPLEID/SAMPLEID.ome.tiff
cell_contour_path	results/SAMPLEID/SAMPLEID_Expanded_5um_cell_contour_coords.csv
nuclei_contour_path	results/SAMPLEID/SAMPLEID_nuclei_contour_coords.csv

The OME-TIFF contains three image planes:

Index (Python tif.pages[i])	Index (R readTIFF(..., all=i))	Contents
0	1	Original-resolution H&E (RGB, variable µm/pixel)
1	2	Tissue mask (1 channel, 1 µm/pixel)
2	3	1 µm/pixel H&E (RGB) — used in the examples below

---

Scanpy (Python)

Tested with Scanpy 1.10.4.

Requirements

pip install scanpy==1.10.4 anndata==0.11.4 pandas==2.3.3 numpy==2.4.2 matplotlib==3.10.8 tifffile==2026.2.24 ome-types==0.5.2

1. Read the Matrix Market files

matrix.mtx.gz is stored as genes × barcodes, so transpose to cells × genes before assigning metadata.

import scanpy as sc
import pandas as pd
import numpy as np

mtx_path = "results/SAMPLEID/SAMPLEID_cell_binned"

adata = sc.read_mtx(f"{mtx_path}/matrix.mtx.gz").T

adata.obs_names = pd.read_csv(
    f"{mtx_path}/barcodes.tsv.gz", sep="\t", header=None, compression="gzip"
)[0].values

genes = pd.read_csv(
    f"{mtx_path}/features.tsv.gz", sep="\t", header=None, compression="gzip"
)
adata.var_names = genes[1].values
adata.var["gene_ids"] = genes[0].values
adata.var["feature_types"] = genes[2].values
adata.var_names_make_unique()

2. Parse spatial coordinates from the barcodes

Each barcode encodes its spatial position as <type>:Y:X (µm, global coordinates).

coords = pd.read_csv(f"{mtx_path}/barcodes.tsv.gz", delimiter=":", header=None)
coords.columns = ["barcode", "Y", "X"]

sp_coords = coords[["X", "Y"]].to_numpy(dtype=np.float32)
adata.obsm = {"spatial": sp_coords}

At this point adata.obsm["spatial"] is in global substrate space and is enough for any analysis that does not need the H&E image.

3. Attach the registered image and convert to local coordinates

To overlay points on the registered OME-TIFF, read Global_top / Global_left from the OME metadata and subtract them from the spatial coordinates so they live in the image's local pixel/µm space.

The 1 µm/pixel H&E is at page index 2 (the third page), which means a scalefactor of 1 (1 pixel = 1 µm).

import json
from tifffile import TiffFile
from ome_types import from_xml

sample_id = "SAMPLEID"
ome_tiff_path = "results/SAMPLEID/SAMPLEID.ome.tiff"

with TiffFile(ome_tiff_path) as tif:
    ome = from_xml(tif.ome_metadata)
    ome_obj = json.loads(ome.json())
    value_dict = ome_obj["structured_annotations"]["map_annotations"][0]["value"]
    Global_top  = int(round(float(value_dict["GlobalPos_top"])))
    Global_left = int(round(float(value_dict["GlobalPos_left"])))

    image = tif.pages[2].asarray()  # 1 um/pixel H&E

adata.uns["spatial"] = {
    sample_id: {
        "images": {"hires": image},
        "scalefactors": {
            "tissue_hires_scalef": 1,
            "spot_diameter_fullres": 1,
        },
    }
}

adata.obsm["spatial"][:, 0] -= Global_left
adata.obsm["spatial"][:, 1] -= Global_top

If you want full-resolution H&E instead, load tif.pages[0] and set tissue_hires_scalef = 1 / pixel_size_um, where pixel_size_um comes from ome_obj["images"][0]["pixels"]["physical_size_x"] / 1000.

4. (Optional) Load cell / nuclei contours

Contour CSVs are also in global coordinates and need the same offset.

def load_contours(path, Global_left, Global_top):
    df = pd.read_csv(path)
    contours = {}
    for cell_id, group in df.groupby("cell_id"):
        coords = group[["vertex_x", "vertex_y"]].values.astype(np.float32)
        coords[:, 0] -= Global_left
        coords[:, 1] -= Global_top
        contours[str(cell_id)] = coords
    return contours

adata.uns["cell_contours"]   = load_contours("results/SAMPLEID/SAMPLEID_Expanded_5um_cell_contour_coords.csv", Global_left, Global_top)
adata.uns["nuclei_contours"] = load_contours("results/SAMPLEID/SAMPLEID_nuclei_contour_coords.csv",            Global_left, Global_top)

5. Quick sanity-check plot

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(8, 8))
ax.imshow(adata.uns["spatial"][sample_id]["images"]["hires"])
ax.scatter(adata.obsm["spatial"][:, 0], adata.obsm["spatial"][:, 1], s=0.5, c="red")
ax.set_aspect("equal")
plt.show()

---

Seurat (R)

Tested with Seurat 5.3.0. The API examples use the Seurat v5 Assay5 class; they will not work on Seurat v4 without modification.

Requirements

The easiest way to get a working environment is via conda / mamba:

mamba create -n seurat_doc -c conda-forge -y \
    r-base=4.3.3 r-seurat=5.3.0 r-matrix=1.6_5 r-jsonlite=2.0.0 r-tiff=0.1_12 r-xml2=1.4.0 r-ggplot2=3.5.2
mamba activate seurat_doc

1. Read the Matrix Market files

Use gzfile() + Matrix::readMM so you don't need any extra .gz reader package.

suppressPackageStartupMessages({
  library(Seurat)
  library(Matrix)
})

mtx_path <- "results/SAMPLEID/SAMPLEID_cell_binned"

mat       <- Matrix::readMM(gzfile(file.path(mtx_path, "matrix.mtx.gz")))
barcodes  <- readLines(gzfile(file.path(mtx_path, "barcodes.tsv.gz")))
genes_raw <- readLines(gzfile(file.path(mtx_path, "features.tsv.gz")))
genes     <- as.data.frame(do.call(rbind, strsplit(genes_raw, "\t")), stringsAsFactors = FALSE)
colnames(genes) <- c("gene_id", "gene_name", "feature_type")

rownames(mat) <- make.unique(genes$gene_name)
colnames(mat) <- barcodes

seurat_obj <- CreateSeuratObject(counts = mat)

# Seurat v5: attach feature-level metadata via AddMetaData on the assay.
feat_meta <- data.frame(
  gene_ids      = genes$gene_id,
  feature_types = genes$feature_type,
  row.names     = rownames(mat)
)
seurat_obj[["RNA"]] <- AddMetaData(seurat_obj[["RNA"]], metadata = feat_meta)

Seurat may warn Feature names cannot have underscores ('_'), replacing with dashes ('-'). That is expected and only renames features in the assay; the gene_ids column above still holds the original Ensembl IDs.

2. Parse spatial coordinates from the barcodes

coords_split <- strsplit(barcodes, ":")
coords_df <- data.frame(
  Y = as.numeric(sapply(coords_split, `[`, 2)),
  X = as.numeric(sapply(coords_split, `[`, 3)),
  row.names = barcodes
)

seurat_obj$spatial_x <- coords_df$X
seurat_obj$spatial_y <- coords_df$Y

sp_coords <- as.matrix(coords_df[, c("X", "Y")])
colnames(sp_coords) <- c("spatial_1", "spatial_2")
seurat_obj[["spatial"]] <- CreateDimReducObject(
  embeddings = sp_coords, key = "spatial_", assay = "RNA"
)

At this point the spatial reduction is in global substrate space.

3. Attach the registered image and convert to local coordinates

tiff::readTIFF(path, payload = FALSE) returns just the metadata (no pixels) for the first IFD, including the OME-XML in description. We pull Global_top / Global_left straight out of that, then read only the page we actually want with all = N.

⚠️ Do not use tiff::readTIFF(path, all = TRUE) on these OME-TIFFs — the original-resolution H&E plane is multiple gigabytes uncompressed and will OOM or segfault.

library(tiff)

sample_id     <- "SAMPLEID"
ome_tiff_path <- "results/SAMPLEID/SAMPLEID.ome.tiff"

# Pull the OME-XML out of the first IFD's ImageDescription tag (no pixels read).
meta <- tiff::readTIFF(ome_tiff_path, payload = FALSE)
desc <- meta$description[1]

# Tags look like: <M K="GlobalPos_top">2467.07</M>
m_top  <- regmatches(desc, regexpr('GlobalPos_top">[^<]+',  desc))
m_left <- regmatches(desc, regexpr('GlobalPos_left">[^<]+', desc))
Global_top  <- as.integer(round(as.numeric(sub('GlobalPos_top">',  "", m_top))))
Global_left <- as.integer(round(as.numeric(sub('GlobalPos_left">', "", m_left))))

# Read only the 1 um/pixel H&E (the 3rd page, 1-indexed).
image       <- tiff::readTIFF(ome_tiff_path, all = 3, as.is = TRUE)[[1]]
scalefactor <- 1

seurat_obj@misc$spatial <- list()
seurat_obj@misc$spatial[[sample_id]] <- list(
  images       = list(hires = image),
  scalefactors = list(
    tissue_hires_scalef   = scalefactor,
    spot_diameter_fullres = 1
  )
)

# Shift global -> local so points overlay the image.
seurat_obj$spatial_x <- seurat_obj$spatial_x - Global_left
seurat_obj$spatial_y <- seurat_obj$spatial_y - Global_top

emb <- seurat_obj[["spatial"]]@cell.embeddings
emb[, 1] <- emb[, 1] - Global_left
emb[, 2] <- emb[, 2] - Global_top
seurat_obj[["spatial"]]@cell.embeddings <- emb

4. (Optional) Load cell / nuclei contours

load_contours <- function(path, Global_left, Global_top) {
  df  <- read.csv(path)
  out <- list()
  for (cid in unique(df$cell_id)) {
    coords <- as.matrix(df[df$cell_id == cid, c("vertex_x", "vertex_y")])
    coords[, 1] <- coords[, 1] - Global_left
    coords[, 2] <- coords[, 2] - Global_top
    out[[as.character(cid)]] <- coords
  }
  out
}

seurat_obj@misc$contours <- list(
  cell_contours   = load_contours("results/SAMPLEID/SAMPLEID_Expanded_5um_cell_contour_coords.csv", Global_left, Global_top),
  nuclei_contours = load_contours("results/SAMPLEID/SAMPLEID_nuclei_contour_coords.csv",            Global_left, Global_top)
)

5. Quick sanity-check plot

The image needs to be loaded as a numeric [0, 1] array for annotation_raster (i.e. drop as.is = TRUE if you used it above when computing offsets — or read it once more without it for plotting).

library(ggplot2)

image_plot <- tiff::readTIFF(ome_tiff_path, all = 3)[[1]]  # numeric in [0, 1]

ggplot() +
  annotation_raster(as.raster(image_plot),
                    xmin = 0, xmax = ncol(image_plot),
                    ymin = 0, ymax = nrow(image_plot)) +
  geom_point(
    data = data.frame(x = seurat_obj$spatial_x, y = seurat_obj$spatial_y),
    aes(x, y), size = 0.5, color = "red"
  ) +
  coord_fixed(xlim = c(0, ncol(image_plot)),
              ylim = c(nrow(image_plot), 0),
              expand = FALSE)

Additional Resources

• Matrix Market format specification: https://math.nist.gov/MatrixMarket/formats.html
• Scanpy documentation: https://scanpy.readthedocs.io/
• AnnData documentation: https://anndata.readthedocs.io/
• Seurat documentation: https://satijalab.org/seurat/