Output Format¶

In plain words

Look at the result explains what the result holds and how to open it, in plain words. This page has every detail.

Thyra converts MSI data into SpatialData objects stored as Zarr directories. This page describes what the output contains and how to work with it.

Loading a Dataset¶

import spatialdata as sd

sdata = sd.read_zarr("output.zarr")

print("Tables:", list(sdata.tables.keys()))
print("Images:", list(sdata.images.keys()))
print("Shapes:", list(sdata.shapes.keys()))

Structure Overview¶

A converted dataset contains the following elements:

Element	Key Pattern	Description
Table	`{dataset_id}_z{z}`	AnnData with intensity matrix (pixels x m/z), coordinates in `.obs`, m/z axis in `.var`
TIC Image	`{dataset_id}_z{z}_tic`	2D total ion current image, shape `(1, y, x)`
Pixel Shapes	`{dataset_id}_z{z}_pixels`	GeoDataFrame with pixel box geometries
Optical Images	`{dataset_id}_optical_{name}`	Microscopy images (when available)

3D mode

When converted with --handle-3d, the _z{z} suffix is dropped and all slices are merged into a single table with x, y, z coordinates in .obs. With more than one slice, the TIC image becomes a single volume of shape (c, z, y, x) rather than one 2D image per slice — see 3D Data / Z-Slices.

The pixel shapes stay two-dimensional: one flat square per pixel per slice, carrying no depth. Read a pixel's depth from the table's spatial_z instead. See Pixel footprints in a volume for why, and Coordinate systems for how the three elements line up.

Default dataset ID

The default dataset_id is msi_dataset, so typical keys look like msi_dataset_z0, msi_dataset_z0_tic, etc. Change it with --dataset-id.

Coordinate systems

Every element above carries a transform to the "global" coordinate system (and the same transform to a second one named after the dataset id). In a micrometer store, elements also carry "acquisition": "global" moved to where the source put its raster, so two stores of one slide line up. Thyra writes a self-describing coordinate_systems metadata attr at the zarr top level so consumers know what "global" is in (micrometers or pixels). See Coordinate Systems for the contract and how to read it.

TIC Images¶

The TIC (Total Ion Current) image shows the sum of all intensities per pixel.

import numpy as np
import matplotlib.pyplot as plt

tic_key = [k for k in sdata.images if k.endswith("_tic")][0]
tic_array = np.asarray(sdata.images[tic_key])[0]  # drop channel dim -> (y, x)

plt.imshow(tic_array, cmap="viridis")
plt.colorbar(label="TIC Intensity")
plt.title("Total Ion Current")
plt.show()

Optical Images¶

Bruker acquisitions can carry microscopy images. With --include-optical, the default, they are stored alongside the MSI data.

optical_keys = [k for k in sdata.images if "optical" in k]
print("Optical images:", optical_keys)

if optical_keys:
    opt_image = np.asarray(sdata.images[optical_keys[0]])
    # Shape is (channels, y, x) -- transpose for display
    opt_rgb = np.moveaxis(opt_image[:3], 0, -1)
    plt.imshow(opt_rgb)
    plt.title("Optical Image")
    plt.show()

TIC-to-Optical Overlay¶

The TIC image carries an affine transform that maps it into the optical image's coordinate space. This means overlays work automatically in tools like napari.

To inspect the transform:

from spatialdata.transformations import get_transformation

tic_element = sdata.images[tic_key]
transforms = get_transformation(tic_element, get_all=True)
affine = list(transforms.values())[0]

matrix = np.array(affine.to_affine_matrix(
    input_axes=("x", "y"), output_axes=("x", "y")
))
print(f"Scale: {matrix[0,0]:.2f}x, {matrix[1,1]:.2f}x")
print(f"Offset: ({matrix[0,2]:.0f}, {matrix[1,2]:.0f})")

How alignment works

The alignment image defines the reference coordinate system. It has an Identity transform, or a translation by the crop's origin when the conversion covers one region and the image was cropped to it (see Supported Formats). The TIC image has an Affine transform (scale + offset) that positions it in the optical coordinate space. This comes from the teaching points in the .mis file (Bruker data): the spots sit on a lattice one raster step apart, and a TIC cell's centre lands on its spot. Without enough teaching points, each region is stretched over its Area (see D31).

Which image is which: `attrs["optical_images"]`¶

An element name is derived, not the filename -- sample_0000.tif and sample_0000.jpg both become {dataset_id}_optical_highres -- and the alignment image is otherwise only implied, by its transform, and only when the alignment was applied. The store says both outright:

optical = sdata.attrs.get("optical_images")
if optical:
    for key, info in optical["elements"].items():
        print(key, "<-", info["source_file"])
    print("alignment image:", optical["alignment_element"])

Field	Meaning
`elements`	One entry per optical element in the store, keyed by element name. `source_file` is the name of the file it was read from (no directory: the acquiring machine's path does not travel with the store). An element cropped to one region also has `crop`: the crop's `origin` and `size` and the whole scan's `full_size`, in image pixels.
`alignment_element`	The element the `.mis` `<ImageFile>` designates -- the one the teaching points are stated in, so the one the other images are scaled into. `None` when the source designates none, or when the designated file is not in this store.

The whole section is absent when the conversion put no optical image in the store, and an image whose pixels could not be read is taken back out of it, so a name here always resolves in sdata.images.

Not on the element

This is in the store's root attrs and not on the images/* groups because per-element attributes do not survive the write: an element group carries exactly ome and spatialdata_attrs, whatever the producer put in the element's .attrs.

Mass Spectrum Data¶

Average Spectrum¶

Each table stores a pre-computed average spectrum in uns:

table_key = list(sdata.tables.keys())[0]
msi_table = sdata.tables[table_key]

mz_values = msi_table.var["mz"].values
avg_spectrum = msi_table.uns["average_spectrum"]

plt.plot(mz_values, avg_spectrum, linewidth=0.5)
plt.xlabel("m/z")
plt.ylabel("Average Intensity")
plt.title("Average Mass Spectrum")
plt.show()

It is that table's own X.mean(axis=0) -- the intensity summed down the m/z axis and divided by the rows the table holds, not by the grid positions it covers and not by the spectra the reader handed in. Read it off the matrix and you get the same numbers back:

import numpy as np

X = msi_table.X.toarray() if hasattr(msi_table.X, "toarray") else msi_table.X
np.allclose(msi_table.uns["average_spectrum"], X.mean(axis=0))   # True

Two consequences of "rows, not spectra". A grid position the acquisition never reached has no row and so is in neither the sum nor the count. A position the source measured twice has one row holding the sum of both spectra, so it contributes two spectra of ion current over one row -- which is what that row really holds.

On a multi-slice source converted as 2D each plane's table carries its own mean, so a plane brighter or dimmer than the average is described by its own spectrum.

Per-Region Average Spectrum¶

For multi-region datasets, Thyra also stores a mean spectrum per acquisition region in uns["average_spectrum_per_region"]. Each key is the region number (as a string), and the value is a 1-D array matching the m/z axis.

if "average_spectrum_per_region" in msi_table.uns:
    per_region = msi_table.uns["average_spectrum_per_region"]
    for region_id, spectrum in per_region.items():
        plt.plot(mz_values, spectrum, label=f"Region {region_id}", linewidth=0.5)
    plt.xlabel("m/z")
    plt.ylabel("Average Intensity")
    plt.legend()
    plt.title("Average Spectrum per Region")
    plt.show()

Each region's vector is the mean of the rows in that region, on the same "rows, not spectra" rule as average_spectrum.

Note

This key is present whenever the reader records regions: on every multi-region dataset, and on every solariX dataset, whose reader records a region for each spectrum even when there is only one. Other datasets only have the global average_spectrum.

Regions span the slices

A region is an in-plane footprint: get_region_map() is keyed on (x, y) and has no z component, so on a multi-slice store one region covers the same area on every plane. average_spectrum_per_region is therefore dataset-wide and identical in every plane's table, while average_spectrum beside it is that plane's alone. On the single-table stores that regions actually occur on, the two are over the same rows.

Intensity Matrix¶

The intensity matrix is stored as a sparse matrix. Each row is one pixel, each column is one m/z bin:

X = msi_table.X  # sparse (pixels x m/z)
print(f"Shape: {X.shape}")
print(f"Non-zero: {X.nnz:,} ({X.nnz / (X.shape[0] * X.shape[1]) * 100:.2f}%)")

Sparse format

The storage is CSC (Compressed Sparse Column) on every route, which is fast for extracting ion images (column = one m/z across all pixels). If you need fast per-pixel access, call X.tocsr() on the matrix you read back -- that is one conversion in memory, and it is why there is no format choice at write time.

The stored matrix is canonical: within every column the row indices are in ascending order and each appears once, whatever order the source delivered its pixels in (a multi-area Bruker acquisition comes back area by area, for instance). A consumer can binary-search a column's indices directly, and scipy reports has_sorted_indices as true without a sort_indices() pass.

One pixel is one row. A source that measures the same coordinate twice -- a processed imzML with a repeated coordinate -- has both spectra summed into that row, which is what the pixel's TIC image and its non_empty_pixels count then describe as well. The conversion says so in the log, naming how many positions it applied to. Reading the store is not how you would find out otherwise: scipy and dask both merge repeated entries as they read, so read_zarr would show the sum whether or not the stored arrays held one entry or two.

Every stored intensity is a real, non-negative number

Non-finite and negative intensities are dropped from a spectrum before it is resampled, and the conversion says how many at WARNING the first time it happens.

A NaN is not a measurement, and one of them makes every aggregate over its column NaN -- the average spectrum, the TIC image, any ion image. A negative value is a baseline subtraction that overshot rather than a smaller measurement, and nothing that reads a store treats it as signal: a TIC, an ion image and a mean spectrum all read it as removing current that was never there.

Dropping them before either resampling method runs is also what makes the two agree. Given a negative value, nearest_neighbor stored it (leaving the pixel with a TIC of 0) while tic_preserving found a non-positive total for the spectrum and zeroed the whole thing, so the same file converted two ways gave two different stores with nothing in either saying why.

The sibling tables built from the raw scans -- the mobility grid, the demultiplexed MS/MS split -- drop the same points by the same rule, so a grid's marginal over channels still reproduces the summed table's column on a source that carries them.

Zero is not affected by this; the sparse write drops explicit zeros as it always did. A non-finite m/z is handled by the axis instead. With --no-resample the source's m/z values are var["mz"], so such a source is refused outright: dropping them would desynchronise every index paired with them. A resampled axis is laid over the finite m/z only, and a point at a non-finite m/z is dropped before either resampling method sees the spectrum, and counted with the intensities above.

Ion Images¶

To visualise the spatial distribution of a specific m/z value, sum the columns inside a window around it. On a resampled axis most bins are empty, so the single bin nearest the target usually gives an all-zero image:

target_mz, tol = 760.5, 0.25
lo, hi = np.searchsorted(mz_values, [target_mz - tol, target_mz + tol])

# Sum the window's columns of the sparse matrix, one value per pixel
ion_values = np.asarray(X[:, lo:hi].sum(axis=1)).ravel()

# Reconstruct image from pixel coordinates
x_coords = msi_table.obs["x"].values.astype(int)
y_coords = msi_table.obs["y"].values.astype(int)

ion_image = np.zeros((y_coords.max() + 1, x_coords.max() + 1))
ion_image[y_coords, x_coords] = ion_values

plt.imshow(ion_image, cmap="hot")
plt.colorbar(label="Intensity")
plt.title(f"m/z {target_mz} +/- {tol}")
plt.show()

Individual Spectra¶

pixel_idx = 0
spectrum = X[pixel_idx].toarray().flatten()

plt.plot(mz_values, spectrum, linewidth=0.5)
plt.xlabel("m/z")
plt.ylabel("Intensity")
plt.title(f"Pixel {pixel_idx}")
plt.show()

var columns¶

var always carries mz (numeric, finite, strictly increasing). Readers whose native axis is not m/z keep that axis alongside it, and annotation tools add formula, adduct, annotation_source and fdr -- names fixed by the metadata schema so they mean the same thing in every store. thyra validate checks the mz contract on every table. msi_metadata.provenance.mass_axis holds a SHA-256 of this axis, and thyra validate checks it.

Ion mobility¶

Mobility is a feature coordinate, next to m/z. It never enters obs, a coordinate system or a transform: those say where a pixel is, mobility says what was measured there. Two things follow.

The MSI table above is always summed over mobility. When the source had a mobility dimension (Bruker TDF with TIMS engaged, or an imzML export that carries a mobility array), the table's uns gains two blocks, and uns["msi_metadata"]["ms_analysis"]["ion_mobility"] carries the schema-level summary. A source without mobility writes none of them, so a consumer can tell "summed over mobility" from "never had any".

uns["mobility_axis"] describes the axis the table was summed over. Plain-name keys (never a CV accession: a colon is not a legal zarr key on Windows), numpy arrays rather than lists:

key	value
`present`	`True`
`type_name`, `type_accession`	the quantity: `inverse reduced ion mobility` / `MS:1002815` (TIMS), or drift time / `MS:1002476`
`unit_name`, `unit_accession`	e.g. `volt-second per square centimeter` / `MS:1002814` for 1/K0
`n_scans`	length of `values` (the TIMS ramp length on a Bruker source, the drift bin count on a Waters one)
`values`	`float64[n_scans]`: the native axis when it is shared across pixels -- the 1/K0 of every TIMS scan (from the vendor calibration, decreasing with scan number) or of every feature of a continuous imzML export, or the drift time in ms of every Waters drift bin (ascending, from 0); absent when mobility is per pixel
`acq_range`	`float64[2]`, `[lower, upper]`: the acquired range as declared, else the span of `values`
`calibration`	`{model_type, coefficients}`: the vendor `TimsCalibration` row, provenance only
`source`	`"bruker_tdf"`, `"imzml"` or `"waters_raw"`
`bin_width`	`float`, only where the source defines its own equal-width bins (Waters drift bins): the width of one, in the axis unit. Its presence is what gives the heatmap one channel per bin (D53)
`resolved_table`	element key of the mobility-resolved sibling table, when one was written

uns["mobility_heatmap"] is the table's mean mass-mobility frame: the raw (m/z, mobility, intensity) points of every pixel, binned and averaged over pixel positions (two spectra on one position count once, as they make one row), accumulated from the raw scan read during conversion. It is the discovery surface -- look at it to see whether mobility separates anything before asking for a mobility-resolved table -- and what a viewer's mobility panel draws. It covers the same pixels as uns["average_spectrum"] beside it: on a multi-slice source converted as 2D, each plane's table holds its own plane's heatmap.

key	value
`mz_edges`	`float64[m + 1]`: bin edges on the common m/z axis, which is coarsened by an integer factor to about 4,000 bins (`m` is the axis length itself when it is shorter)
`mobility_edges`	`float64[k + 1]`, `k = 256`: equal-width bins in the axis unit, ascending, spanning the smallest to the largest mobility value the source holds (the axis `values` when shared; every pixel's own values on a processed imzML, read once for this). A source whose axis is fewer than 256 equal-width bins of its own (Waters drift bins, `uns["mobility_axis"]["bin_width"]` set) gets one channel per bin instead, edges half a bin outside the first and last value (D53)
`counts`	`float32[m, k]`: mean intensity per bin over pixel positions, the rows `uns["average_spectrum"]` divides by
`current_ratio`	`float`: the heatmap's ion current over `uns["average_spectrum"]`'s -- 1.0 under `scan_sum`; above 1 under `vendor_centroid` (about 1.01 to 1.23), whose picked peaks keep only part of the raw current the heatmap holds

k is not tuned: the mobility grid table below defaults to the same channels over the same edges -- literally the same rule and the same generator -- so a box drawn on the heatmap maps onto grid channels by integer index in both directions. Read k off mobility_edges or the shape of counts, never assume 256: on a Waters run it is the drift bin count, 200 on every run measured, because 256 channels over 200 bins would leave 56 of them empty in a fixed comb. The heatmap keeps the m/z points the summed table keeps, edge bins included. The m/z binning is the converter's own nearest-bin rule, coarsened, so under a lossless summed spectrum (--tdf-spectrum scan_sum) the heatmap summed over mobility, counts.sum(axis=1), equals uns["average_spectrum"] coarsened to mz_edges. Under vendor_centroid it does not: the centroid keeps only the current inside the peaks its picker assigns (79 to 99 percent on the imaging acquisitions measured, see Design Decisions) and merges bins, while the heatmap is built from every raw point. Which case a given store is in is recorded rather than left to be found by subtraction: current_ratio is the heatmap's total over the stored mean spectrum's, the same number uns["mobility_marginal"] carries for the grid table, and the conversion says it at WARNING when it is not 1. On a Bruker source the heatmap costs one extra library call per frame (about a millisecond); --no-mobility-heatmap skips it.

heat = msi_table.uns["mobility_heatmap"]
counts = np.asarray(heat["counts"])              # (m, k)
mz_centres = np.asarray(heat["mz_edges"])
mz_centres = (mz_centres[:-1] + mz_centres[1:]) / 2
k0_centres = np.asarray(heat["mobility_edges"])
k0_centres = (k0_centres[:-1] + k0_centres[1:]) / 2
plt.pcolormesh(k0_centres, mz_centres, np.log1p(counts))
plt.xlabel("1/K0 (V s cm^-2)"); plt.ylabel("m/z")

Both blocks are summaries of the whole dataset. Nothing in them is per pixel, and nothing in them enters obs or a coordinate system.

When every pixel shares one set of (m/z, mobility) feature pairs -- a continuous imzML export with a mobility array, as TIMSImaging writes -- Thyra also writes a mobility-resolved sibling table, {table}_mobility, unless --no-mobility-table is given:

	MSI table `{id}_z0`	mobility table `{id}_z0_mobility`
rows	pixels	the same pixels, same `obs`, same `region`
`var["mz"]`	strictly increasing, unique	non-decreasing with duplicates
`var["mobility"]`	absent	the feature's 1/K0 (or drift time)
sort	by `mz`	lexicographic `(mz, mobility)`
`var["ccs"]`	absent	the feature's collision cross section in square angstroms, for an assumed charge of 1; on a 1/K0 axis only (see below)
also		`mz_index` (column on the MSI axis), `mobility_index` (rank of the mobility value, or the grid channel), `uns["feature_axis"]`, `uns["mobility_axis"]`, `uns["ccs"]` with `var["ccs"]`

Two isomers at one m/z that separate in mobility are one column in the MSI table and two in the mobility table. The (mz, mobility) sort means an m/z window is one contiguous column block whose mobility-summed image equals the MSI table's, and a mobility window is a mask inside that block.

Consumers must discriminate on "mobility" in var.columns, never on the element name: the mobility table is not a second MSI dataset, and code that assumes a unique, strictly increasing var["mz"] must not be pointed at it. thyra validate applies the pair contract (non-decreasing mz, unique sorted pairs) to tables carrying mobility and the strict contract to all others.

Collision cross section¶

When the mobility axis is 1/K0 in V s/cm2 (uns["mobility_axis"]["type_accession"] is MS:1002815, its unit_accession MS:1002814), the mobility table also carries var["ccs"]: each feature's collision cross section (CCS), in square angstroms, worked out from its own mz and mobility by the Mason-Schamp equation. Every TIMS route qualifies: a Bruker TDF with --mobility-grid, and a continuous imzML whose mobility array states 1/K0 in that unit. A drift-time axis gets no column: turning a drift time into a CCS needs the drift cell's length, field and pressure, which no source Thyra reads states. Neither does an axis whose quantity or unit the source leaves unstated.

The value rests on three assumptions, and uns["ccs"] on the same table records them:

key	value
`quantity_accession`, `quantity_name`	`MS:1002954`, collisional cross sectional area
`unit_accession`, `unit_name`	`UO:0000324`, square angstrom
`method`	`"Mason-Schamp"`
`charge_assumed`	`1`. MALDI ions are almost always singly charged; no source states a feature's charge
`drift_gas`, `drift_gas_mass_da`	`"nitrogen"`, `28.0134`
`temperature_k`	`305.0`
`conditions_source`	`"convention"`: the gas and temperature are the ones Bruker's own conversion uses, not conditions read from the run

For an ion of charge z the true value is between sqrt(z) and z times the stored one, because the charge also changes the ion's mass in the reduced mass. Recompute it rather than rescale it:

from thyra.core.mobility import mason_schamp_ccs

mob = sdata.tables["msi_dataset_z0_mobility"]
ccs_2plus = mason_schamp_ccs(mob.var["mobility"], mob.var["mz"], 2)

The same formula fills the column on every route and every platform, so a store does not depend on which vendor library the machine has. Bruker's own library gives the same values to 1e-6 relative from m/z 100 up, and to 2.2e-6 at m/z 50. On a grid table, mz is the mass bin and mobility the channel centre, so the value is that cell's, not one peak's. Two features that share an m/z keep the order of their mobility: the larger 1/K0 has the larger CCS.

The same table from a common mobility grid¶

A source whose mobility values differ per pixel -- a Bruker TDF, where a frame is a point cloud and no two pixels are promised the same (m/z, 1/K0) pairs, or a processed imzML export with a mobility array, as TIMSCONVERT writes it -- has no shared feature axis to read off. --mobility-grid fills the same {table}_mobility element for it by binning: every pixel's points go onto one set of mobility channels shared across the conversion, and (m/z bin, channel) becomes the feature axis.

The two mechanisms produce the same kind of table -- same element key, same var columns, same sort, same discriminator -- and a consumer does not need to tell them apart to read either. What says which one filled it is uns["mobility_grid"], present only on a binned table:

key	value
`law`	how the channel edges are spaced; `"linear"` (equal width in the axis unit) is the only law today
`lower`, `upper`	the range the channels span, in the axis unit
`n_channels`	how many channels: by default the heatmap's own, 256, or one per bin on a source with fewer equal-width bins (D53)
`channel_width`	`(upper - lower) / n_channels`, recorded so a reader can see it without arithmetic
`edges`	`float64[n_channels + 1]`, the channel edges themselves

msi_metadata.ms_analysis.ion_mobility.grid carries law, lower, upper and n_channels too, so a consumer reading only the summed table finds them.

Three things are worth knowing before asking for one:

The edges come from the values, never a declared range. A real TDF's per-scan 1/K0 overhangs its declared OneOverK0AcqRange by a few scans (1.00003 to 1.29133 against a declared 1.0 to 1.29 on one measured acquisition), and uns["mobility_heatmap"] already bins over the values. A processed imzML has no shared axis, so its values are read from every pixel, once, before the heatmap's pass (D46). --mobility-min / --mobility-max override them, at the cost of the alignment below.
The default channels are an alignment anchor, not a tuning knob. They are the heatmap's own channels over the heatmap's own edges (256, or one per Waters drift bin), so a box drawn on the heatmap selects grid channels by integer index with no resampling and no edge off-by-one. --mobility-bins changes them and gives that up, even when it says 256. The width the anchor realizes on a typical 0.29 1/K0 span is 0.0011, finer than the 0.002 to 0.02 band TIMS resolving power supports; that is said at INFO and the count is not moved for it, because the alignment is worth more than the size.
Its marginal reproduces the summed table under the default --tdf-spectrum scan_sum. A grid table is built from raw scans, so its marginal reproduces the summed table only when the summed table was built from the same scans, which the default is. On an imzML that means each point placed on its nearest mass bin: no resampling, or nearest_neighbor. An explicit vendor_centroid, a peak-picked spectrum over the same ramp, is kept and said at WARNING; uns["mobility_marginal"] then records by how much the two differ.

The marginal invariant. Summing a grid table's channels within one m/z bin reproduces that bin's column of the summed table, per pixel. That is what uns["mobility_marginal"] records rather than merely asserting:

key	value
`summed_table`	element key of the table the marginal is compared against
`current_ratio`	total ion current of the grid table over the summed table's; exactly `1.0` under `scan_sum`
`current_ratio_pixel_min` / `_max`	the same ratio across pixels

The comparison is per pixel, one bounded pass over each table's memmaps. A per-cell deviation is not stored: it would need the marginal and its difference from the summed table in memory, each as large as the summed table (see Design Decisions). The snippet below computes it from the two stored matrices when it is wanted.

grid = sdata.tables["msi_dataset_z0_mobility"]
mz_index = grid.var["mz_index"].to_numpy()
marginal = np.zeros(sdata.tables["msi_dataset_z0"].n_vars)
np.add.at(marginal, mz_index, np.asarray(grid.X[0].todense()).ravel())
# equals the summed table's row 0, to floating point

A grid is refused, at INFO or WARNING and never as an exception, when:

refused when	because
the source has no mobility dimension	there is nothing to bin
`--mobility-grid` was not given	binning is opt in: it costs a pass over the source and a much larger table
the source gives no mobility values to bin over	a Bruker reader opened without its vendor library cannot supply them, and the declared range is not a substitute
the grid spans more pairs than the counting pass can hold (a raw, unresampled axis of millions of bins)	the count array is `4 bytes x bins x channels`, capped at 1 GB; resample to fewer mass bins
the `var` frame of the occupied `(m/z bin, channel)` pairs is projected to take more than half of the machine's free memory (338 bytes per pair: 330 measured, plus the `ccs` column), or the pairs pass an absolute cap of 100,000,000	the count, the projection and the free memory are printed; resample to fewer mass bins or ask for fewer channels. A projection past a quarter of free memory is attempted with a `WARNING`. See Design Decisions

How it is built, and why its size does not matter. The table is built the way the summed table is built, in two passes over the raw scans. The first pass counts, per (m/z bin, channel) cell of the grid, how many pixels occupy it -- a dense count over the grid's span, 142 MB on a default-resampled timsTOF axis, sized before the first pixel is read and independent of how many pixels there are. The second pass scatters each pixel's cells straight into memmapped CSC arrays in a scratch directory next to the output (.thyra_mobility_*, removed once the table is written), 12 bytes of disk per stored non-zero. Nothing is ever held for the whole image: memory is the count array plus one frame. On a Bruker TDF both passes are the summed table's own: the converter reads each frame once per pass and derives the summed spectrum, the heatmap's points, the grid's cells and the MS/MS split from that one read, so none of the siblings adds a read of the source (see Design Decisions). On a processed imzML the two passes are the heatmap's pass, which also counts, and one more for the scatter, after a read of the mobility arrays alone for the range. Every point is mapped onto the mass axis once per pass and shared between the heatmap and the grid, since the mapping costs more than the vendor read.

The size ceiling is on the pairs that carry signal, not on the pairs the grid spans. The two differ by an order of magnitude -- 200 frames of a measured timsTOF acquisition occupied 3.9M of a possible 35.5M -- so refusing on the span would turn away conversions that fit ninefold over. The span is said at INFO when it passes the ceiling; the count is what refuses, checked the moment the first pass ends and before a single value has been scattered, so a refusal costs one read and no memory.

Bruker TDF and processed imzML exports with a mobility array get a grid. A continuous imzML export with a mobility array already has a shared feature axis and is read off it; --mobility-grid is ignored there, with a WARNING. A point whose mobility value is not a finite number has no channel: the heatmap and the grid leave it out, the summed table keeps it, and a WARNING gives the count.

Fragmentation (MS/MS)¶

An MS/MS imaging run measures fragments, so var["mz"] is fragment m/z. Nothing about the axis says so -- a converted MS/MS store is otherwise shaped exactly like an MS1 one -- so when the source reports fragmentation the table carries uns["msms_schedule"]:

key	value
`ms_level`	`2` for a fragment spectrum; the block is absent for MS1
`n_windows`	number of distinct isolation windows in the schedule: isolation events, of which one precursor can own several (the same mass at two collision energies is two windows); when `constant_across_pixels` is false, not every pixel has all of them
`merges_precursors`	`True` when the windows isolate more than one distinct precursor (by target m/z and mobility scan range), so the stored spectrum sums them; one precursor fragmented twice does not merge
`constant_across_pixels`	whether every pixel was fragmented on the same schedule
`isolation_window_target`	`float64[n]`: the isolated m/z of each window
`isolation_window_lower_offset` / `_upper_offset`	`float64[n]`: the window spans `target - lower` to `target + upper`
`collision_energy`	`float64[n]`: in electronvolts
`scan_begin` / `scan_end`	`int64[n]`: the mobility scans each window occupies, when the source separates them that way (Bruker PASEF and DIA frames)
`resolved_table`	element key of the demultiplexed sibling table, when one was written

Field names and CV terms follow mzPeak's spectra_metadata_precursors: ms_level is MS:1000511, the isolation terms are MS:1000827 / 828 / 829, and activation is MS:1000133 with MS:1000045 collision energy in UO:0000266. The accessions travel as values (..._accession keys), never as dict keys -- a colon is not a legal Windows path character and zarr writes a key as a directory name. A source reporting a single full isolation width has it halved into two equal offsets.

The same facts appear in the versioned schema block as uns["msi_metadata"]["ms_analysis"]["fragmentation"], where the precursor list is a JSON string (a list of objects does not round-trip through AnnData/zarr); read_msi_metadata_blocks and thyra validate decode it.

A multi-precursor pixel is a chimera

The MSI table sums a frame into one spectrum per pixel. When merges_precursors is True, that spectrum holds fragments of every precursor the frame isolated, with nothing marking which came from which -- so it must not be read as the fragment spectrum of any one of them. Conversion says so at WARNING, and by default also writes the split apart as a second table, below (--no-msms-table turns that off); the MSI table itself is unchanged either way.

if "msms_schedule" in msi_table.uns:
    sched = msi_table.uns["msms_schedule"]
    print("MS level:", sched["ms_level"])
    for mz, ce in zip(sched["isolation_window_target"], sched["collision_energy"]):
        print(f"  precursor {mz:.3f} at {ce:.1f} eV")

Demultiplexed MS/MS table¶

When the source isolates several precursors per pixel in disjoint mobility scan ranges (Bruker PASEF frames, or type-9 frames of one DIA window group, the layout MALDI iprm-PASEF is reported to use: a scheduled list of precursors serially fragmented at every pixel), Thyra also writes them split apart as a sibling table, {table}_msms, by default (--no-msms-table opts out; see Design Decisions). Each precursor's block is its precursor ion image, and each column inside the block is one fragment's image:

	MSI table `{id}_z0`	MS/MS table `{id}_z0_msms`
rows	pixels	the same pixels, same `obs`, same `region`
a column	one m/z bin, all precursors summed	one m/z bin of one precursor
`var["mz"]`	strictly increasing, unique	the fragment m/z, restarting at every precursor
`var["precursor_mz"]`	absent	the isolated m/z the fragments came from
sort	by `mz`	lexicographic `(precursor_mz, precursor_mobility, mz)`
also		`precursor_mobility` (the 1/K0 it was isolated at), `precursor_index` (its position in this store's precursor axis), `mz_index` (column on the MSI axis), `uns["feature_axis"]`, `uns["msms_schedule"]`, `uns["demultiplexed_current"]`

The fragment axis is the MSI table's own mass axis: var["mz"] is msi.var["mz"][var["mz_index"]], so a column of this table and the corresponding column of the summed table are the same m/z bin. One precursor's fragments are therefore a contiguous column block whose row sums are its ion image, and the blocks add back up: summing all fragment columns of every precursor reproduces the summed table's TIC per pixel, when every scan that carries current lies in an isolation window. Exactly, under the default --tdf-spectrum scan_sum: the split is built from the raw scans, and only a summed spectrum built from the same scans can add back up to it. Two things break the identity, and both are said at WARNING and recorded in uns["demultiplexed_current"]. An explicit vendor_centroid keeps only the current inside the peaks it picks (79 to 99 percent on the imaging acquisitions measured), so the split holds more than the summed table. A schedule whose windows leave out scans that carry current puts that current in no block, so the split holds less. Every acquired PASEF schedule measured so far covers all of it (see Design Decisions).

msms = sdata.tables["msi_dataset_z0_msms"]
# Look the precursor up once, then slice on its block index: never
# compare precursor_mz with == , and never assume an m/z is unique.
var = msms.var
index = var["precursor_index"][np.argmin(np.abs(var["precursor_mz"] - 936.578))]
block = (var["precursor_index"] == index).to_numpy()
image = np.asarray(msms.X[:, block].sum(axis=1)).ravel()

Consumers must discriminate on "precursor_mz" in var.columns, never on the element name. thyra validate checks that precursor_mz is non-decreasing, that each precursor_index owns one contiguous block, and that mz increases strictly inside it; the strict single-column contract applies to every other table.

Two precursors can share an m/z

A method may isolate the same mass at two mobility positions -- that is how an isomer pair is targeted, and separating them is what the mobility dimension is for. Thyra keeps them as two column blocks and never sums them back together. So:

the block identity is precursor_index, not precursor_mz;
(precursor_mz, mz) may legitimately repeat, which is why validation does not use that pair;
precursor_mobility is what tells the two apart (the 1/K0 at the middle of the window's scan range -- a window spans a slice of the ramp, and a scheduled method reports no apex).

Aligning two datasets

precursor_index is a position in one store's precursor axis and means nothing outside it: two samples whose schedules differ in length give the same index to different precursors. Align on (precursor_mz, precursor_mobility) -- by matching the pair, never by comparing it with ==: two runs acquired back to back agreed on the 1/K0 bitwise, while two of the same method on different days came out up to 4.0e-3 apart, because each file carries its own TIMS calibration. The var index labels are named after the precursor's m/z for the same reason -- p936.578_mz1732, never p14_mz1732 -- so anndata.concat cannot silently merge two unrelated precursors. Measured on two real acquisitions in opposite polarities, 15 precursors against 13 with none in common: zero labels shared, against 5,262 a rank-named scheme would have shared, every one of them naming two different precursors. Differing mass axes do not save you there -- both started at m/z 50, so the low mz_index values line up. Concatenate with join="inner" to keep a var at all: on an outer join every column missing from one side makes merge="unique" drop the whole column, while on an inner join it keeps precursor_mz, mz and mz_index and drops exactly precursor_index and precursor_mobility -- anndata finding, on its own, which of the four travel between stores. The table carries no mobility

column: the scan range is how the precursors are separated, not what they are indexed by, and a table matching both discriminators would tell a consumer nothing about which kind it holds.

The split is a filter on the scan number, never an estimate, so Thyra refuses rather than approximates and says at INFO which condition failed. They are asked in this order, which is by how much each says about the acquisition rather than by how cheap it is to check:

refused when	because
the schedule varies from pixel to pixel	the precursors are not a global feature axis
the source records no precursor for the fragment frames	there is nothing to separate them by
there is one precursor	the summed table already is its fragment spectrum
the isolation windows overlap, or carry no scan range	they cannot be separated by mobility alone

The order matters because the last three are facts about what the source recorded and the first is a fact about the method. A diaPASEF file gives each fragment frame one of several window groups, which take turns frame by frame, and keeps survey frames among the fragment ones. Its windows overlap across groups, but only because the groups take turns: the varying schedule is the design, and the overlap its consequence. A data-dependent (DDA) PASEF run trips the first row and the last at once -- one measured holds 14,952 distinct windows, each in one to eight frames, overlapping on the ramp -- and the varying schedule is the design while the overlap is only its consequence.

A refusal writes no sibling table, is never an exception and never touches the summed table. Bruker TDF is the only source that reports what the split needs today.

uns["demultiplexed_current"] records how much of the summed table's ion current the split holds: current_ratio over the whole image, and current_ratio_pixel_min / _max across pixels. Under the default --tdf-spectrum scan_sum it is exactly 1.0 when the windows cover every scan that carries current, and below 1 when they do not: the synthetic type-9 test fixture, whose three windows span scans 60 to 180 of a 240-scan ramp, gives 0.92. Under an explicit --tdf-spectrum vendor_centroid it is above 1 -- the vendor peak picker drops the index bins it assigns to no peak while the split reads raw scans, which on a real acquisition is about 1.5% overall and up to 1.14x on a single pixel. The two tables genuinely do not add up in that mode, and this block is where the store says so. It compares the two tables' per-pixel ion current, one bounded pass over each.

Deliberate limits

The feature axis depends on the data. Only (precursor, bin) pairs that carry signal become columns, so two datasets converted with identical settings get different var. The alternative is precursors x the whole mass axis -- millions of empty columns -- and is not worth it. Align on the intrinsic columns, as above.
The fragment axis is borrowed from MS1. It is the summed table's mass axis, so a resampling grid chosen for intact ions also sets fragment resolution. That coupling is what makes mz_index meaningful and the conservation check exact; it is a choice, not a necessity.
Do not select precursors by float equality. Look the precursor up once and slice on precursor_index within that store.
Built out of core, like the mobility grid. Two passes over the precursor spectra -- count the occupied (precursor, bin) pairs, then scatter into memmapped CSC arrays in a scratch directory next to the output (.thyra_msms_*) -- so memory is the count array (4 bytes x precursors x mass bins) plus one frame, whatever the pixel count. Those are the summed table's own two passes, fed from the same frame read, so the split adds no read of the source. The largest acquisition this has run on is still 713 pixels with 15 precursors; var grows with the occupied pairs.
Two acquired schedules have been tested, from one instrument. 15 precursors in positive mode and 13 in negative sharing none of them, over four images of 487 to 713 pixels, plus a third shape made by deleting one window from a copy, the single-precursor refusal on a real MsMsType = 2 acquisition and the varying-schedule refusal on real DDA-PASEF (see Design Decisions). All four come from the same instrument and method family, and none exceeds 713 pixels.
Type-9 (DIA-layout) frames are tested on a synthetic file only. No MALDI iprm-PASEF run has been converted yet, so which frame type it writes is not confirmed. The split of a one-group schedule is checked point by point on the synthetic fixture. A run whose survey frames share spots with fragment frames is refused, and those pixels sum both MS levels (said at WARNING). Windows in PrmFrameMsMsInfo (prm-PASEF) are not read.
A label means the same bin only when the axes match. mz_index is a column of that store's mass axis, and the default resampling grid follows the acquisition range, so two runs share an axis only when they share that range -- 599,146 bins to m/z 1000 against 635,610 to 1200 on the two measured. Check var["mz"], not just the label, before reading across stores.

Relation to other MS/MS imaging representations

The open formats solve this at the raw layer by never merging: an imzML or mzML spectrum carries its own precursor, and mzPeak links a spectrum row to a precursor table by index. The chimera is created by the analysis layer's one-spectrum-per-pixel model, so this table is that per-spectrum precursor reference translated onto a feature axis.

Feature-based workflows (MZmine's SIMSEF, for instance) instead attach one representative MS2 spectrum to each MS1 feature, which identifies the feature but keeps no spatial information about the fragments. This table is the stronger form of the same data: summing a precursor's block collapses it to that per-precursor ion image, while keeping the block gives every fragment its own image -- which is what a spatial co-localisation check between a fragment and its precursor needs.

The (precursor_mz, mz) feature-pair layout is Thyra's own; no open analysis-layer convention for per-precursor ion images exists to follow. The vocabulary is not: ms_level, the isolation window terms and the activation terms are PSI-MS, spelled as mzPeak spells them.

Pixel Coordinates¶

Coordinates are stored in the table's .obs DataFrame:

print("Columns:", list(msi_table.obs.columns))

Column	Type	Description
`x`, `y`	int	Raster grid coordinates (pixel indices)
`spatial_x`, `spatial_y`	float	Physical coordinates in micrometers
`region`	categorical	SpatialData annotation target: the key of the pixel shapes element, the same on every row
`region_number`	int	The acquisition region the pixel was measured in
`acquisition_order`	int64	The order the pixels were acquired in, when the source records it (see below)
`instance_key`	str	A copy of the index, which SpatialData uses to link each row to its pixel shape

The DataFrame index is instance_id (a string pixel identifier).

region and region_number are unrelated. region tells SpatialData which element the table annotates; to split pixels by acquisition region, use region_number.

Acquisition order¶

The rows follow the raster: left to right, top to bottom. The instrument need not have measured them in that order. A serpentine scan comes back along every second row, and a slide with several regions is measured one region after another. So a drift plotted against the row number is plotted against the raster, not against time. acquisition_order gives the time order:

obs = msi_table.obs
as_acquired = obs.sort_values("acquisition_order")

It is an int64 that grows with acquisition time and differs from row to row. Where the source numbers its spectra, it is that number; otherwise it is the spectrum's position in the order the source lists them. It need not start at 0 or run without gaps, so sort by it or use it to find a spectrum in the source, but do not use it as a row index.

Source	`acquisition_order` is
Bruker timsTOF (`.tsf`, `.tdf`)	`Frames.Id`, from 1, as the file stores it. A `--region` conversion keeps the region's own ids
Bruker solariX	`Spectra.Id` in `peaks.sqlite`
imzML	The spectrum's position in the file's spectrum list, from 0 (its mzML `index`). imzML does not say whether the writer listed spectra as they were acquired; the list is the only order the file has
mzPeak	The spectrum's `spectrum_index`
Waters	The scan's number in the run, from 0. MassLynx counts scans per function, and the count carries on through the functions a long raster is split into
PHI	Absent. Each frame passes over the whole raster and a pixel's spectrum sums all of them, so every pixel was measured throughout the acquisition
Bruker rapifleX	Absent. The reader walks the raster, not the acquisition's position log

A missing column means the order is unknown, not that it is the raster order. A position measured twice is one row holding the sum of both spectra, in every table, and its acquisition_order is the earlier one. The sibling tables (_mobility, _msms) carry the same column, row for row. Why the column looks like this, and what else was considered, is in design decision D26.

Regions¶

Datasets acquired from multi-region slides (e.g., multiple tissue sections on one slide) store region information in two places. The timsTOF, solariX, Rapiflex and mzPeak readers record regions; every other format writes region 1 for each pixel. thyra.list_regions gives the same summary from the raw input, before converting.

Per-Pixel Region Number¶

The region_number column in .obs indicates which acquisition region each pixel belongs to:

print(msi_table.obs["region_number"].value_counts())

Region Summary¶

Region metadata -- including human-readable names from the instrument's Area definitions and each region's bounds (x_min, y_min, x_max, y_max, 0-based, in the store's pixel grid) -- is stored as JSON in uns:

import json

regions = json.loads(msi_table.uns["regions"])
for r in regions:
    print(f"Region {r['region_number']}: {r.get('name', 'unnamed')} "
          f"({r['n_spectra']:,} spectra)")

Example output:

Region 0: E2506 (104,321 spectra)
Region 1: Matrix (1,053 spectra)

Filtering by region

To work with only one region:

mask = msi_table.obs["region_number"] == 0
tissue_table = msi_table[mask]

3D Data / Z-Slices¶

By default, Thyra stores each z-slice as a separate table and TIC image:

slice_tables = sorted(k for k in sdata.tables if "_z" in k)
print(f"{len(slice_tables)} z-slices: {slice_tables}")

# Access a single slice
z0_table = sdata.tables[slice_tables[0]]
print(f"Slice 0: {z0_table.shape}")

3D mode (`--handle-3d`)¶

When converted with --handle-3d, all slices are combined into a single table with x, y, z coordinates in .obs, and the per-slice TIC images are replaced by one volume:

Element	Key	Shape / dims
Table	`{dataset_id}`	pixels x m/z, with `x`, `y`, `z` and `spatial_x`, `spatial_y`, `spatial_z` in `.obs`
TIC volume	`{dataset_id}_tic`	`(c, z, y, x)` — one channel, then the three spatial axes
Pixel Shapes	`{dataset_id}_pixels`	GeoDataFrame of 2D pixel boxes — one per pixel per slice, no z

dataset_id = "msi_dataset"  # the default; see --dataset-id
volume = sdata.images[f"{dataset_id}_tic"]
print(volume.dims)                      # ('c', 'z', 'y', 'x')

arr = np.asarray(volume.data)[0]        # drop channel -> (z, y, x)
plt.imshow(arr[0], cmap="viridis")      # first slice

Note the axis order is (c, z, y, x), not (c, x, y, z): index a slice with arr[z], not arr[..., z].

Voxel depth¶

The volume's transform to "global" scales by two distinct numbers — the in-plane pixel pitch for x and y, and the slice spacing for z. It first moves each voxel back half a step, so a voxel's centre lands on its row's spatial_x, spatial_y and spatial_z:

from spatialdata.transformations import get_transformation

axes = ("c", "z", "y", "x")
matrix = get_transformation(volume, to_coordinate_system="global").to_affine_matrix(
    input_axes=axes, output_axes=axes
)
print(matrix[1, 1])   # um per slice step
print(matrix[3, 3])   # um per pixel in x
print(matrix[1, 4])   # minus half a slice step

The slice spacing comes from --z-spacing. When nothing supplies one, Thyra falls back to the in-plane pitch and records that it did:

cs = sdata.attrs["coordinate_systems"]["global"]
cs["z_spacing_um"]      # the number used
cs["z_spacing_source"]  # "manual" | "automatic" | "assumed_isotropic"

A z_spacing_source of "assumed_isotropic" means nobody supplied a spacing and the in-plane pitch was reused — treat the depth as unknown rather than as measured. Section thickness is set by the microtome, not by the raster, so the two agree only by coincidence. See --z-spacing.

These keys only appear on volumes

z_spacing_um and z_spacing_source are written only when the store actually holds a multi-slice volume. Their absence is how a 2D store says it has no z axis, which is why convention_version stays at 1 — the keys are additive and a consumer that never reads 3D sees the schema it already knows.

Pixel footprints in a volume¶

The pixel polygons are two-dimensional, on every route including a volume. A slice's depth is not on the geometry; read it from the table's spatial_z, or from the volume's own transform:

dataset_id = "msi_dataset"  # the default; see --dataset-id
shapes = sdata.shapes[f"{dataset_id}_pixels"]
print(shapes.geometry.iloc[0].has_z)    # False, on 2D and 3D alike

# Depth per pixel, in micrometres, from the table:
obs = sdata.tables[dataset_id].obs
print(obs.loc["0", "spatial_z"])        # z_index * z_spacing_um

Every footprint therefore resolves to the correct place in x and y, and carries no claim about z. Overlaying the shapes on the volume is exact in-plane and needs spatial_z to pick the slice.

Why the footprints are flat

They briefly were not. v3.2.0 made them POLYGON Z at the depth of their slice, which is geometrically the more honest representation, and it broke spatialdata's spatial queries: a bounding box enclosing an entire test volume returned 26 of 30 footprints and 26 of 30 table rows, with no exception and no warning. A z-restricted query returned the same rows whether z was inside or far outside the data.

spatialdata asks for 2D here: ShapesModel.validate warns that a 3-dimensional geometry column "could led to unexpected behaviors" and names force_2d() as the remedy. 3D shapes are not on its roadmap (#109 has been idle since June 2023), and the live 2.5D discussion (#961) covers points, images and labels only. Serial-section MSI is 2.5D in that sense.

So this is a deliberate trade: a documented gap in z, in exchange for queries that return every pixel. Expressing the depth as a Translation on flat geometry does not close the gap either — the transform silently drops z. Both negative results are pinned by tests/unit/converters/test_3d_pixel_shapes_z.py, which will fail if upstream changes. See Coordinate systems.

Dataset Metadata¶

Global metadata¶

Stored in sdata.attrs:

if "pixel_size_x_um" in sdata.attrs:
    print(
        f"Pixel size: {sdata.attrs['pixel_size_x_um']} x "
        f"{sdata.attrs['pixel_size_y_um']} um"
    )

if "msi_dataset_info" in sdata.attrs:
    info = sdata.attrs["msi_dataset_info"]
    print(f"Dimensions: {info.get('dimensions_xyz')}")
    print(f"Non-empty pixels: {info.get('non_empty_pixels'):,}")

Table-level metadata¶

Instrument info, acquisition parameters, and resampling config are in msi_table.uns:

print("uns keys:", list(msi_table.uns.keys()))

Provenance¶

uns["essential_metadata"] records what the store was made from, and how the source was interpreted. It is the only place a converted dataset says where it came from, so it is written the same way by every converter path:

provenance = msi_table.uns["essential_metadata"]

print(provenance["source_path"])     # the file this was converted from
print(provenance["dimensions"])      # source grid, [x, y, z]
print(provenance["mass_range"])      # SOURCE m/z range, not the target axis
print(provenance["spectrum_type"])   # "centroid spectrum" / "profile spectrum"
print(provenance["thyra_version"])   # the Thyra that wrote the store

mass_range describes the source, not the resampled axis -- for the axis the data actually sits on, read msi_table.var["mz"].

spectrum_type is the representation the reader delivered. It is also in msi_metadata as ms_analysis.spectrum_representation (see below), except when Thyra only guessed it.

Beside it, when the source format provides them:

Key	Contents
`format_specific`	Vendor metadata (imzML file mode and UUID; flexImaging areas and teaching points for rapiflex; `instrument_source_type` for Bruker tsf/tdf, Bruker's own code for what ionised the sample, recorded raw because its labels are not documented in the file; for Bruker tsf/tdf also `instrument_calibration`, the m/z calibration the run started with as the analysis database's `CalibrationInfo` states it, and `calibration`, when the acquisition has a `calibration.sqlite`: how many states it holds and which one the conversion applied -- with `--no-recalibrated` none, so those keys are empty and the newest state is under `latest_calibration_*`)
`acquisition_params`	Polarity, scan range, laser settings, timestamp and method name, in the vendor's spelling and unit
`instrument_info`	Instrument model, serial, software version
`raw_metadata`	Source metadata as read, for round-trip fidelity; for a Bruker acquisition with a `.mis` file, `mis_metadata` is that file as parsed
`regions`	Acquisition region summary, as a JSON string (see Regions)

A section the source format has nothing for is omitted rather than written empty, so "instrument_info" not in uns means "this format does not carry it" rather than "it was carried and lost".

These sections hold what the source states, with two exceptions that apply to every format alike. A field that names a person -- who calibrated the instrument, who operated it, a contact's name or the details that reach them (address, e-mail, phone) -- is not copied. A path the vendor recorded on the acquisition PC keeps only its last component, the file name. Neither does anything in a store, and the source still holds both. The paths Thyra records about where it read the source stay whole: provenance.source_path above, and in format_specific the data_path of a Bruker tsf/tdf or Waters source (for Bruker also database_path and binary_file), the ibd_file of an imzML and the mis_file of a solariX acquisition. The store's root attributes repeat three of these sections, as format_specific_metadata, acquisition_parameters and instrument_information, and follow the same rule.

Within these sections, a list that holds anything besides numbers -- imzML cvParams (a list of objects) is the main case -- is stored as a JSON string; decode it with json.loads:

import json

cv_params = json.loads(msi_table.uns["raw_metadata"]["cvParams"])
print(cv_params[0])   # {"name": ..., "accession": "MS:...", "value": ...}

Purely numeric lists stay plain arrays. The JSON encoding exists because AnnData/zarr cannot round-trip such lists faithfully (a list of objects comes back as repr strings), and because a stored list of strings reads back as a numpy string array -- and on numpy 2.1-2.2, deepcopying such an array segfaults the Python process outright (numpy#28609, fixed in numpy 2.3). Every table copy deepcopies uns (AnnData.copy, spatialdata.polygon_query, joins), so a store carrying one would crash those readers with no traceback.

Stores written by Thyra <= 3.7.2 on numpy 2.1-2.2

Older stores still carry cvParams as a string array. On numpy 2.1.x-2.2.x (Google Colab ships 2.1.3), convert those arrays to plain lists right after loading, before anything copies the table:

from thyra.metadata import sanitize_uns_string_arrays

for table in sdata.tables.values():
    table.uns = sanitize_uns_string_arrays(table.uns)

Environments on numpy 2.0 or >= 2.3 are unaffected either way.

Structured metadata: `uns["msi_metadata"]`¶

The sections above preserve what the source said, in the source's own vocabulary. uns["msi_metadata"] is the normalised, versioned view of the same facts: fixed field names, PSI-MS/NCBITaxon/UBERON/CHEBI ontology terms, and a schema a validator can hold it to.

block = msi_table.uns["msi_metadata"]

print(block["schema_version"])                    # "0.14.0"
print(block["ms_analysis"]["pixel_size_um"])      # {"x": 20.0, "y": 20.0}
print(block["ms_analysis"].get("spectrum_representation"))  # "profile" / "centroid"
print(block.get("sample", {}).get("sample_id"))   # only when you named the sample
print(block["provenance"]["source_format"])       # "imzml"
print(block.get("acquisition"))                   # see below; absent for imzML
print(block.get("calibration"))                   # see below; absent for imzML
print(block.get("alignment"))                     # see below; absent for imzML

block["acquisition"] is the normalised view of the run itself, filled from acquisition_params by every reader that has the facts and absent (not empty) when a reader has none of them. One line per field, with its unit:

Field	Meaning and unit
`acquisition_datetime`	Start of the acquisition, ISO 8601 `YYYY-MM-DDThh:mm:ss[.fff]`; carries a UTC offset only when the source recorded one (Bruker tsf/tdf and solariX do, PHI and Waters do not)
`laser_power_percent`	Laser power as the percentage of the laser's range the vendor software shows; Bruker tsf/tdf and solariX
`laser_frequency_hz`	Laser repetition rate in hertz; Bruker tsf/tdf and solariX
`shots_per_pixel`	Laser shots summed into one pixel's spectrum; Bruker tsf/tdf, solariX and rapiflex
`method_file`	File name of the acquisition method (`.m` for Bruker, `.EXP` for Waters, `*.par` for rapiflex), never a path

The raw vendor values stay in acquisition_params under their own spellings (laser_frequency on Bruker tsf/tdf, laser_rep_rate on solariX, acquisition_date on PHI and Waters, and so on), so nothing is lost by the normalisation.

block["calibration"] is the same kind of view of how the source's m/z values were calibrated, filled by Bruker tsf/tdf, PHI and Waters and absent for the other formats. "The calibration" is the one the values rest on: the most recent recalibration when there is one, otherwise the calibration the acquisition ran under.

Field	Meaning and unit
`calibration_datetime`	When the calibration was made, ISO 8601 at the precision the source records (Waters to the minute), with a UTC offset only when the source recorded one (Bruker does)
`recalibrated`	Whether a calibration made after the acquisition replaced the one it ran under; Bruker tsf/tdf with a `calibration.sqlite`, and PHI
`original_calibration_datetime`	For recalibrated data, when the calibration the recalibration replaced was made
`software`, `software_version`	What made the calibration, in the source's words; Bruker tsf/tdf
`n_reference_peaks`	How many reference peaks the calibration was fitted to; Bruker tsf/tdf and PHI
`mz_standard_deviation_ppm`	The reference peaks' residual after the calibration: the square root of their summed squared m/z errors in ppm over one less than their number. Never zero and never over fewer than two peaks: neither is a measurement, and vendors write both where they have nothing to state
`lock_mass_corrected`	Whether the m/z values were corrected against a lock mass; Waters, as MassLynx reports it

Which calibration a conversion applied is not in the section, because a conversion can be told to apply one the source does not consider current; it is the m/z calibration step of processing (PSI-MS MS:1001485). See Metadata Schema for where each vendor states each field and what is deliberately left out.

block["alignment"] says which optical image the source registers its raster onto, and how. It comes from the .mis file beside a Bruker tsf/tdf, rapiflex or solariX acquisition, and is absent for everything else:

Field	Meaning and unit
`optical_image_file`	File name of the optical image, the `.mis` `<ImageFile>`; never a path
`method`	`"teaching points"` when the `.mis` has the three that register the image
`teaching_points`	One `{image_x_px, image_y_px, stage_x_um, stage_y_um}` per point: a pixel of the image, from its top-left corner, and the stage position of the same feature in micrometres. Stored as a JSON string, like `processing`

When the store holds that image, attrs["optical_images"]["alignment_element"] names it.

It is written by every converter path identically, validated by thyra validate, and exported to a METASPACE submission by thyra export-metaspace. See Metadata Schema for the full contract.

Recipes¶

Side-by-side TIC and ion image¶

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

ax1.imshow(tic_array, cmap="viridis")
ax1.set_title("TIC")

target_mz, tol = 760.5, 0.25
lo, hi = np.searchsorted(mz_values, [target_mz - tol, target_mz + tol])
ion_values = np.asarray(X[:, lo:hi].sum(axis=1)).ravel()

x_coords = msi_table.obs["x"].values.astype(int)
y_coords = msi_table.obs["y"].values.astype(int)
ion_image = np.zeros((y_coords.max() + 1, x_coords.max() + 1))
ion_image[y_coords, x_coords] = ion_values

ax2.imshow(ion_image, cmap="hot")
ax2.set_title(f"m/z {target_mz} +/- {tol}")

plt.tight_layout()
plt.show()

Export ion image to TIFF¶

from PIL import Image

# Normalise to 0-255
ion_norm = (ion_image / ion_image.max() * 255).astype(np.uint8)
Image.fromarray(ion_norm).save("ion_image.tiff")

Top N most intense m/z values¶

avg = msi_table.uns["average_spectrum"]
top_n = 10
top_indices = np.argsort(avg)[-top_n:][::-1]

for idx in top_indices:
    print(f"  m/z {mz_values[idx]:.4f}  avg intensity {avg[idx]:.1f}")

Summary statistics¶

print(f"Dataset: {table_key}")
print(f"  Pixels: {msi_table.n_obs:,}")
print(f"  m/z bins: {msi_table.n_vars:,}")
print(f"  m/z range: {mz_values.min():.2f} -- {mz_values.max():.2f}")
print(f"  Sparsity: {(1 - X.nnz / (X.shape[0] * X.shape[1])) * 100:.1f}%")
if "pixel_size_x_um" in sdata.attrs:
    print(
        f"  Pixel size: {sdata.attrs['pixel_size_x_um']} x "
        f"{sdata.attrs['pixel_size_y_um']} um"
    )

Read both axes

pixel_size_x_um and pixel_size_y_um are equal on a square raster, which is nearly every acquisition -- but a DESI method with DesiXStep != DesiYStep is not square, and the store carries the real pitch on each axis in every block that states one: these attrs, coordinate_systems.global and its raster_to_global_affine, the image and shapes transforms, obs["spatial_x"]/["spatial_y"], the pixel footprints and msi_metadata.ms_analysis.pixel_size_um. Reading only pixel_size_x_um as the pixel size renders such a raster squashed by y/x. Pass --pixel-size to declare the raster square instead, which applies one number to both axes.