Explore a Thyra Dataset¶

This notebook shows how to read Mass Spectrometry Imaging (MSI) data from a SpatialData store produced by Thyra.

For installation and the conversion step that produces such a store, see the Tutorial. This notebook picks up from its output.

Topics covered:

• Loading a SpatialData zarr dataset
• Accessing TIC (Total Ion Current) images
• Accessing optical images and MSI-optical overlay
• Accessing the average mass spectrum
• Pixel counts and data representation
• Accessing raw MSI data and pixel coordinates
• Extracting ion images at specific m/z values
• Working with 3D data and z-slices
• Ion mobility: the axis, the mass-mobility heatmap, the grid table
• Tandem MS: the isolation schedule and the demultiplexed sibling table

Runs with no configuration. If you have no dataset yet, the setup cell below generates a small synthetic one and converts it, so every cell can be executed as-is. To use your own data instead, set the THYRA_TUTORIAL_ZARR environment variable before starting Jupyter.

Setup¶

The cell below resolves which dataset to open:

1. THYRA_TUTORIAL_ZARR, if that environment variable is set
2. otherwise example_data/synthetic_brain.zarr, generating and converting it on first run

Generating and converting the synthetic dataset takes a few seconds and needs no vendor SDK, so it works on Windows, macOS and Linux alike.

Note that convert_msi() called without a resampling_config keeps the original m/z axis, whereas the thyra command-line tool resamples by default. Everything below works either way -- only the number of columns in the table differs.

In [1]:

import warnings

# Silence noisy third-party import warnings so the output below stays
# readable. These come from dask and pkg_resources, not from Thyra.
warnings.filterwarnings("ignore", category=FutureWarning)
warnings.filterwarnings("ignore", message=".*pkg_resources is deprecated.*")

import json
import os
import shutil
import sqlite3
from pathlib import Path

import numpy as np
import matplotlib.pyplot as plt
import scipy.sparse as sparse
import spatialdata as sd

# Use your own dataset by setting THYRA_TUTORIAL_ZARR:
#   Windows (cmd)   set THYRA_TUTORIAL_ZARR=D:\data\my_dataset.zarr
#   macOS / Linux   export THYRA_TUTORIAL_ZARR=/data/my_dataset.zarr
# Otherwise a small synthetic dataset is generated and converted here.
DEFAULT_ZARR = Path("example_data/synthetic_brain.zarr")
ZARR_PATH = Path(os.environ.get("THYRA_TUTORIAL_ZARR") or DEFAULT_ZARR)


def build_example_dataset(zarr_path):
    """Generate the synthetic example dataset and convert it."""
    from thyra import convert_msi
    from thyra.tools.make_example_data import generate_example_imzml

    imzml = zarr_path.with_suffix(".imzML")
    if not imzml.exists():
        print("Generating the synthetic example dataset...")
        generate_example_imzml(imzml)

    print("Converting...")
    if not convert_msi(imzml, zarr_path):
        shutil.rmtree(zarr_path, ignore_errors=True)
        raise RuntimeError("Conversion failed -- see the log output above.")


if not ZARR_PATH.exists():
    if ZARR_PATH != DEFAULT_ZARR:
        raise FileNotFoundError(
            f"THYRA_TUTORIAL_ZARR points at a path that does not exist: {ZARR_PATH}"
        )
    build_example_dataset(DEFAULT_ZARR)

sdata = sd.read_zarr(ZARR_PATH)
print(f"Loaded: {ZARR_PATH}")
print("Tables:", list(sdata.tables.keys()))
print("Images:", list(sdata.images.keys()))
print("Shapes:", list(sdata.shapes.keys()))

import warnings # Silence noisy third-party import warnings so the output below stays # readable. These come from dask and pkg_resources, not from Thyra. warnings.filterwarnings("ignore", category=FutureWarning) warnings.filterwarnings("ignore", message=".*pkg_resources is deprecated.*") import json import os import shutil import sqlite3 from pathlib import Path import numpy as np import matplotlib.pyplot as plt import scipy.sparse as sparse import spatialdata as sd # Use your own dataset by setting THYRA_TUTORIAL_ZARR: # Windows (cmd) set THYRA_TUTORIAL_ZARR=D:\data\my_dataset.zarr # macOS / Linux export THYRA_TUTORIAL_ZARR=/data/my_dataset.zarr # Otherwise a small synthetic dataset is generated and converted here. DEFAULT_ZARR = Path("example_data/synthetic_brain.zarr") ZARR_PATH = Path(os.environ.get("THYRA_TUTORIAL_ZARR") or DEFAULT_ZARR) def build_example_dataset(zarr_path): """Generate the synthetic example dataset and convert it.""" from thyra import convert_msi from thyra.tools.make_example_data import generate_example_imzml imzml = zarr_path.with_suffix(".imzML") if not imzml.exists(): print("Generating the synthetic example dataset...") generate_example_imzml(imzml) print("Converting...") if not convert_msi(imzml, zarr_path): shutil.rmtree(zarr_path, ignore_errors=True) raise RuntimeError("Conversion failed -- see the log output above.") if not ZARR_PATH.exists(): if ZARR_PATH != DEFAULT_ZARR: raise FileNotFoundError( f"THYRA_TUTORIAL_ZARR points at a path that does not exist: {ZARR_PATH}" ) build_example_dataset(DEFAULT_ZARR) sdata = sd.read_zarr(ZARR_PATH) print(f"Loaded: {ZARR_PATH}") print("Tables:", list(sdata.tables.keys())) print("Images:", list(sdata.images.keys())) print("Shapes:", list(sdata.shapes.keys()))

Loaded: example_data\synthetic_brain.zarr
Tables: ['msi_dataset_z0']
Images: ['msi_dataset_z0_tic']
Shapes: ['msi_dataset_z0_pixels']

SpatialData Structure Overview¶

A Thyra-converted dataset contains:

Element	Key pattern	Description
Table	{dataset_id}_z{i}	AnnData table with intensity matrix (pixels x m/z), coordinates in .obs, m/z axis in .var, average spectrum in .uns
TIC Image	{dataset_id}_z{i}_tic	2D image of total ion current per pixel, shape (c, y, x) with c=1
Pixel Shapes	{dataset_id}_z{i}_pixels	GeoDataFrame with pixel geometries
Optical Image	{dataset_id}_optical_{suffix}	Microscopy image (if available)

{i} is the slice number, so a 2D dataset has _z0: with the default dataset_id the keys are msi_dataset_z0, msi_dataset_z0_tic and msi_dataset_z0_pixels. A dataset converted with --handle-3d drops the _z{i} part and holds all slices in one table.

Accessing TIC Images¶

The TIC image represents the sum of all intensities across the mass spectrum for each pixel.

In [2]:

# Find TIC image key -- match "_tic" suffix to avoid matching "optical"
tic_keys = [key for key in sdata.images.keys() if key.endswith("_tic")]
print("TIC image keys:", tic_keys)

tic_key = tic_keys[0]
tic_image = sdata.images[tic_key]

# Convert to numpy -- the stored shape is (c, y, x) with c=1
# Use np.asarray() for robust access regardless of backend (xarray, dask, zarr)
tic_array = np.asarray(tic_image)[0]  # drop channel dim -> (y, x)

print(f"TIC image shape: {tic_array.shape}")
print(f"Intensity range: {tic_array.min():.2f} -- {tic_array.max():.2f}")

# Find TIC image key -- match "_tic" suffix to avoid matching "optical" tic_keys = [key for key in sdata.images.keys() if key.endswith("_tic")] print("TIC image keys:", tic_keys) tic_key = tic_keys[0] tic_image = sdata.images[tic_key] # Convert to numpy -- the stored shape is (c, y, x) with c=1 # Use np.asarray() for robust access regardless of backend (xarray, dask, zarr) tic_array = np.asarray(tic_image)[0] # drop channel dim -> (y, x) print(f"TIC image shape: {tic_array.shape}") print(f"Intensity range: {tic_array.min():.2f} -- {tic_array.max():.2f}")

TIC image keys: ['msi_dataset_z0_tic']

TIC image shape: (36, 48)
Intensity range: 16176.91 -- 146471.64

In [3]:

# Visualize the TIC image
plt.figure(figsize=(10, 8))
plt.imshow(tic_array, cmap="viridis")
plt.colorbar(label="TIC Intensity")
plt.xlabel("X (pixels)")
plt.ylabel("Y (pixels)")
plt.title("Total Ion Current (TIC) Image")
plt.tight_layout()
plt.show()

# Visualize the TIC image plt.figure(figsize=(10, 8)) plt.imshow(tic_array, cmap="viridis") plt.colorbar(label="TIC Intensity") plt.xlabel("X (pixels)") plt.ylabel("Y (pixels)") plt.title("Total Ion Current (TIC) Image") plt.tight_layout() plt.show()

Accessing Optical Images¶

When converted with include_optical=True, Thyra stores microscopy images alongside the MSI data. Keys follow the pattern {dataset_id}_optical_{suffix} where suffix is highres, derived, or derived from the filename.

Optical images and MSI-to-optical alignment come from vendor metadata, so they are present for Bruker acquisitions and absent for plain imzML input -- including the synthetic example dataset. The next few cells therefore print a notice instead of a figure unless you are pointing at a Bruker conversion. See Tutorial, Part 2 for a dataset that exercises this path.

In [4]:

# Find optical image keys (exclude TIC images)
optical_keys = [key for key in sdata.images.keys() if "optical" in key]
print("Optical image keys:", optical_keys)

if optical_keys:
    # The primary optical image (highres) has an Identity transform --
    # it defines the coordinate system that the TIC affine maps into.
    # Other optical images may have Scale transforms to align with it.
    from spatialdata.transformations import get_transformation, Identity

    primary_key = None
    for key in optical_keys:
        transforms = get_transformation(sdata.images[key], get_all=True)
        for t in transforms.values():
            if isinstance(t, Identity):
                primary_key = key
                break
        if primary_key:
            break

    # Fallback: pick the one with "highres" in the name, or the largest
    if primary_key is None:
        highres = [k for k in optical_keys if "highres" in k]
        primary_key = highres[0] if highres else optical_keys[0]

    print(f"\nPrimary optical image: {primary_key}")
    opt_image = np.asarray(sdata.images[primary_key])  # shape: (c, y, x)
    print(f"Shape: {opt_image.shape} (channels, y, x), dtype: {opt_image.dtype}")

    # Display -- transpose (c, y, x) -> (y, x, c) for RGB
    if opt_image.shape[0] == 3:
        opt_rgb = np.moveaxis(opt_image, 0, -1)  # (y, x, 3)
    elif opt_image.shape[0] == 1:
        opt_rgb = opt_image[0]  # grayscale
    else:
        opt_rgb = np.moveaxis(opt_image[:3], 0, -1)  # take first 3 channels

    plt.figure(figsize=(10, 8))
    plt.imshow(opt_rgb)
    plt.title(f"Optical Image -- {primary_key}")
    plt.axis("off")
    plt.tight_layout()
    plt.show()
else:
    print("No optical images found (convert with include_optical=True)")

# Find optical image keys (exclude TIC images) optical_keys = [key for key in sdata.images.keys() if "optical" in key] print("Optical image keys:", optical_keys) if optical_keys: # The primary optical image (highres) has an Identity transform -- # it defines the coordinate system that the TIC affine maps into. # Other optical images may have Scale transforms to align with it. from spatialdata.transformations import get_transformation, Identity primary_key = None for key in optical_keys: transforms = get_transformation(sdata.images[key], get_all=True) for t in transforms.values(): if isinstance(t, Identity): primary_key = key break if primary_key: break # Fallback: pick the one with "highres" in the name, or the largest if primary_key is None: highres = [k for k in optical_keys if "highres" in k] primary_key = highres[0] if highres else optical_keys[0] print(f"\nPrimary optical image: {primary_key}") opt_image = np.asarray(sdata.images[primary_key]) # shape: (c, y, x) print(f"Shape: {opt_image.shape} (channels, y, x), dtype: {opt_image.dtype}") # Display -- transpose (c, y, x) -> (y, x, c) for RGB if opt_image.shape[0] == 3: opt_rgb = np.moveaxis(opt_image, 0, -1) # (y, x, 3) elif opt_image.shape[0] == 1: opt_rgb = opt_image[0] # grayscale else: opt_rgb = np.moveaxis(opt_image[:3], 0, -1) # take first 3 channels plt.figure(figsize=(10, 8)) plt.imshow(opt_rgb) plt.title(f"Optical Image -- {primary_key}") plt.axis("off") plt.tight_layout() plt.show() else: print("No optical images found (convert with include_optical=True)")

Optical image keys: []
No optical images found (convert with include_optical=True)

Overlaying MSI Data on Optical Image¶

The optical image typically shows the full slide (e.g. multiple tissue sections), while the TIC only covers the measured region. Thyra stores an Affine transform on the TIC image that maps its raster indices into the optical image's pixel coordinate space. We extract this transform and use it to place the TIC at the correct position on the optical.

In [5]:

from scipy.ndimage import zoom as ndi_zoom
from spatialdata.transformations import get_transformation

if optical_keys:
    opt_h, opt_w = opt_rgb.shape[:2]
    tic_h, tic_w = tic_array.shape

    # Extract the affine transform stored on the TIC image.
    # This is a 3x3 matrix: [[sx, 0, tx], [0, sy, ty], [0, 0, 1]]
    # mapping TIC raster indices (x, y) to optical pixel coordinates.
    tic_element = sdata.images[tic_key]
    transform = get_transformation(tic_element, get_all=True)

    # Pick the dataset-specific or global transform
    t = list(transform.values())[0]
    affine_matrix = np.array(t.to_affine_matrix(input_axes=("x", "y"), output_axes=("x", "y")))
    print(f"TIC affine transform:\n{affine_matrix}")

    sx, sy = affine_matrix[0, 0], affine_matrix[1, 1]
    tx, ty = affine_matrix[0, 2], affine_matrix[1, 2]

    # Compute where the TIC lands on the optical image
    # TIC cell (0,0) spans [0, 1) x [0, 1), so its corner maps to (tx, ty)
    # and the far corner (w,h) to (tx + w*sx, ty + h*sy)
    opt_x_start = tx
    opt_y_start = ty
    opt_x_end = opt_x_start + tic_w * sx
    opt_y_end = opt_y_start + tic_h * sy

    print(f"\nTIC region on optical: x=[{opt_x_start:.0f}, {opt_x_end:.0f}], "
          f"y=[{opt_y_start:.0f}, {opt_y_end:.0f}]")
    print(f"Optical size: {opt_w}x{opt_h}, TIC size: {tic_w}x{tic_h}")
    print(f"Scale: {sx:.2f}x, {sy:.2f}x per TIC pixel")

    # Resize TIC to match optical pixel scale
    tic_resized = ndi_zoom(tic_array, (abs(sy), abs(sx)), order=1)

    # Create an optical-sized canvas and place the TIC at the correct position
    tic_on_optical = np.zeros((opt_h, opt_w))
    y0 = int(round(opt_y_start))
    x0 = int(round(opt_x_start))
    # Clip to stay within optical bounds
    src_y0 = max(0, -y0)
    src_x0 = max(0, -x0)
    dst_y0 = max(0, y0)
    dst_x0 = max(0, x0)
    copy_h = min(tic_resized.shape[0] - src_y0, opt_h - dst_y0)
    copy_w = min(tic_resized.shape[1] - src_x0, opt_w - dst_x0)

    if copy_h > 0 and copy_w > 0:
        tic_on_optical[dst_y0:dst_y0 + copy_h, dst_x0:dst_x0 + copy_w] = \
            tic_resized[src_y0:src_y0 + copy_h, src_x0:src_x0 + copy_w]

    # Normalize for overlay
    tic_max = tic_on_optical.max()
    tic_norm = tic_on_optical / tic_max if tic_max > 0 else tic_on_optical

    # Plot: optical, TIC placed, overlay
    fig, axes = plt.subplots(1, 3, figsize=(20, 7))

    axes[0].imshow(opt_rgb)
    axes[0].set_title("Optical")
    axes[0].axis("off")

    axes[1].imshow(opt_rgb)
    mask = tic_norm > 0.01
    tic_overlay = np.ma.masked_where(~mask, tic_norm)
    axes[1].imshow(tic_overlay, cmap="hot", alpha=0.6, vmin=0, vmax=1)
    axes[1].set_title("Optical + TIC Overlay")
    axes[1].axis("off")

    # Zoomed view of just the TIC region
    pad = 50  # pixels of padding around the TIC region
    axes[2].imshow(opt_rgb)
    axes[2].imshow(tic_overlay, cmap="hot", alpha=0.6, vmin=0, vmax=1)
    axes[2].set_xlim(max(0, x0 - pad), min(opt_w, x0 + tic_resized.shape[1] + pad))
    axes[2].set_ylim(min(opt_h, y0 + tic_resized.shape[0] + pad), max(0, y0 - pad))
    axes[2].set_title("Zoomed Overlay")
    axes[2].axis("off")

    plt.tight_layout()
    plt.show()
else:
    print("No optical image available for overlay")

from scipy.ndimage import zoom as ndi_zoom from spatialdata.transformations import get_transformation if optical_keys: opt_h, opt_w = opt_rgb.shape[:2] tic_h, tic_w = tic_array.shape # Extract the affine transform stored on the TIC image. # This is a 3x3 matrix: [[sx, 0, tx], [0, sy, ty], [0, 0, 1]] # mapping TIC raster indices (x, y) to optical pixel coordinates. tic_element = sdata.images[tic_key] transform = get_transformation(tic_element, get_all=True) # Pick the dataset-specific or global transform t = list(transform.values())[0] affine_matrix = np.array(t.to_affine_matrix(input_axes=("x", "y"), output_axes=("x", "y"))) print(f"TIC affine transform:\n{affine_matrix}") sx, sy = affine_matrix[0, 0], affine_matrix[1, 1] tx, ty = affine_matrix[0, 2], affine_matrix[1, 2] # Compute where the TIC lands on the optical image # TIC cell (0,0) spans [0, 1) x [0, 1), so its corner maps to (tx, ty) # and the far corner (w,h) to (tx + w*sx, ty + h*sy) opt_x_start = tx opt_y_start = ty opt_x_end = opt_x_start + tic_w * sx opt_y_end = opt_y_start + tic_h * sy print(f"\nTIC region on optical: x=[{opt_x_start:.0f}, {opt_x_end:.0f}], " f"y=[{opt_y_start:.0f}, {opt_y_end:.0f}]") print(f"Optical size: {opt_w}x{opt_h}, TIC size: {tic_w}x{tic_h}") print(f"Scale: {sx:.2f}x, {sy:.2f}x per TIC pixel") # Resize TIC to match optical pixel scale tic_resized = ndi_zoom(tic_array, (abs(sy), abs(sx)), order=1) # Create an optical-sized canvas and place the TIC at the correct position tic_on_optical = np.zeros((opt_h, opt_w)) y0 = int(round(opt_y_start)) x0 = int(round(opt_x_start)) # Clip to stay within optical bounds src_y0 = max(0, -y0) src_x0 = max(0, -x0) dst_y0 = max(0, y0) dst_x0 = max(0, x0) copy_h = min(tic_resized.shape[0] - src_y0, opt_h - dst_y0) copy_w = min(tic_resized.shape[1] - src_x0, opt_w - dst_x0) if copy_h > 0 and copy_w > 0: tic_on_optical[dst_y0:dst_y0 + copy_h, dst_x0:dst_x0 + copy_w] = \ tic_resized[src_y0:src_y0 + copy_h, src_x0:src_x0 + copy_w] # Normalize for overlay tic_max = tic_on_optical.max() tic_norm = tic_on_optical / tic_max if tic_max > 0 else tic_on_optical # Plot: optical, TIC placed, overlay fig, axes = plt.subplots(1, 3, figsize=(20, 7)) axes[0].imshow(opt_rgb) axes[0].set_title("Optical") axes[0].axis("off") axes[1].imshow(opt_rgb) mask = tic_norm > 0.01 tic_overlay = np.ma.masked_where(~mask, tic_norm) axes[1].imshow(tic_overlay, cmap="hot", alpha=0.6, vmin=0, vmax=1) axes[1].set_title("Optical + TIC Overlay") axes[1].axis("off") # Zoomed view of just the TIC region pad = 50 # pixels of padding around the TIC region axes[2].imshow(opt_rgb) axes[2].imshow(tic_overlay, cmap="hot", alpha=0.6, vmin=0, vmax=1) axes[2].set_xlim(max(0, x0 - pad), min(opt_w, x0 + tic_resized.shape[1] + pad)) axes[2].set_ylim(min(opt_h, y0 + tic_resized.shape[0] + pad), max(0, y0 - pad)) axes[2].set_title("Zoomed Overlay") axes[2].axis("off") plt.tight_layout() plt.show() else: print("No optical image available for overlay")

No optical image available for overlay

Accessing the Average Mass Spectrum¶

The average spectrum is the mean intensity across all non-empty pixels for each m/z value. It is stored in the .uns attribute of the main data table.

In [6]:

# Get the main MSI table
table_key = list(sdata.tables.keys())[0]
msi_table = sdata.tables[table_key]
print(f"Table: {table_key}, shape: {msi_table.shape} (pixels x m/z values)")

# Average spectrum from .uns
avg_spectrum = msi_table.uns["average_spectrum"]

# m/z axis from .var
mz_values = msi_table.var["mz"].values

print(f"m/z range: {mz_values.min():.4f} -- {mz_values.max():.4f}")
print(f"Number of m/z bins: {len(mz_values)}")

# Get the main MSI table table_key = list(sdata.tables.keys())[0] msi_table = sdata.tables[table_key] print(f"Table: {table_key}, shape: {msi_table.shape} (pixels x m/z values)") # Average spectrum from .uns avg_spectrum = msi_table.uns["average_spectrum"] # m/z axis from .var mz_values = msi_table.var["mz"].values print(f"m/z range: {mz_values.min():.4f} -- {mz_values.max():.4f}") print(f"Number of m/z bins: {len(mz_values)}")

Table: msi_dataset_z0, shape: (1728, 4000) (pixels x m/z values)
m/z range: 250.0000 -- 1200.0000
Number of m/z bins: 4000

In [7]:

# Plot the average mass spectrum
plt.figure(figsize=(12, 5))
plt.plot(mz_values, avg_spectrum, linewidth=0.5)
plt.xlabel("m/z")
plt.ylabel("Average Intensity")
plt.title("Average Mass Spectrum")
plt.xlim(mz_values.min(), mz_values.max())
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# Plot the average mass spectrum plt.figure(figsize=(12, 5)) plt.plot(mz_values, avg_spectrum, linewidth=0.5) plt.xlabel("m/z") plt.ylabel("Average Intensity") plt.title("Average Mass Spectrum") plt.xlim(mz_values.min(), mz_values.max()) plt.grid(True, alpha=0.3) plt.tight_layout() plt.show()

Pixel Count and Data Representation¶

MSI data is stored in a rectangular grid, but not all grid positions contain measurements. The TIC image makes it easy to distinguish data pixels from empty ones.

In [8]:

# Count pixels
total_pixels = tic_array.size
non_empty_pixels = np.count_nonzero(tic_array)

print(f"Total pixels in grid:          {total_pixels:,}")
print(f"Non-empty pixels (with data):  {non_empty_pixels:,}")
print(f"Coverage:                      {non_empty_pixels / total_pixels * 100:.1f}%")

# Count pixels total_pixels = tic_array.size non_empty_pixels = np.count_nonzero(tic_array) print(f"Total pixels in grid: {total_pixels:,}") print(f"Non-empty pixels (with data): {non_empty_pixels:,}") print(f"Coverage: {non_empty_pixels / total_pixels * 100:.1f}%")

Total pixels in grid:          1,728
Non-empty pixels (with data):  1,728
Coverage:                      100.0%

Accessing Raw MSI Data¶

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

In [9]:

X = msi_table.X  # sparse matrix (pixels x m/z)
print(f"Data matrix type: {type(X).__name__}")
print(f"Shape: {X.shape} (pixels x m/z values)")

if sparse.issparse(X):
    nnz = X.nnz
    total = X.shape[0] * X.shape[1]
    print(f"Non-zero entries: {nnz:,} / {total:,} ({nnz / total * 100:.2f}%)")

X = msi_table.X # sparse matrix (pixels x m/z) print(f"Data matrix type: {type(X).__name__}") print(f"Shape: {X.shape} (pixels x m/z values)") if sparse.issparse(X): nnz = X.nnz total = X.shape[0] * X.shape[1] print(f"Non-zero entries: {nnz:,} / {total:,} ({nnz / total * 100:.2f}%)")

Data matrix type: csc_matrix
Shape: (1728, 4000) (pixels x m/z values)
Non-zero entries: 6,885,838 / 6,912,000 (99.62%)

Pixel Coordinates¶

In [10]:

# Coordinates are stored in .obs
print("Available obs columns:", list(msi_table.obs.columns))

# Access spatial coordinates
x_coords = msi_table.obs["x"].values.astype(float)
y_coords = msi_table.obs["y"].values.astype(float)

print(f"X range: {x_coords.min():.1f} -- {x_coords.max():.1f}")
print(f"Y range: {y_coords.min():.1f} -- {y_coords.max():.1f}")
print(f"Number of pixels with coordinates: {len(x_coords):,}")

# Coordinates are stored in .obs print("Available obs columns:", list(msi_table.obs.columns)) # Access spatial coordinates x_coords = msi_table.obs["x"].values.astype(float) y_coords = msi_table.obs["y"].values.astype(float) print(f"X range: {x_coords.min():.1f} -- {x_coords.max():.1f}") print(f"Y range: {y_coords.min():.1f} -- {y_coords.max():.1f}") print(f"Number of pixels with coordinates: {len(x_coords):,}")

Available obs columns: ['y', 'x', 'region', 'spatial_x', 'spatial_y', 'region_number', 'instance_key']
X range: 0.0 -- 47.0
Y range: 0.0 -- 35.0
Number of pixels with coordinates: 1,728

Extracting Individual Spectra¶

In [11]:

# Get spectrum for a specific pixel
pixel_idx = 0
if sparse.issparse(X):
    pixel_spectrum = X[pixel_idx].toarray().flatten()
else:
    pixel_spectrum = np.asarray(X[pixel_idx]).flatten()

plt.figure(figsize=(12, 5))
plt.plot(mz_values, pixel_spectrum, linewidth=0.5)
plt.xlabel("m/z")
plt.ylabel("Intensity")
plt.title(f"Mass Spectrum -- Pixel {pixel_idx}")
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# Get spectrum for a specific pixel pixel_idx = 0 if sparse.issparse(X): pixel_spectrum = X[pixel_idx].toarray().flatten() else: pixel_spectrum = np.asarray(X[pixel_idx]).flatten() plt.figure(figsize=(12, 5)) plt.plot(mz_values, pixel_spectrum, linewidth=0.5) plt.xlabel("m/z") plt.ylabel("Intensity") plt.title(f"Mass Spectrum -- Pixel {pixel_idx}") plt.grid(True, alpha=0.3) plt.tight_layout() plt.show()

Extracting Ion Images¶

To visualise the spatial distribution of a specific m/z value, extract the corresponding intensities from the matrix and reshape them into image coordinates.

The cell below starts from the most intense peak in the average spectrum, which is guaranteed to be a populated bin. To map an arbitrary target m/z, integrate over a tolerance window instead -- see Searching for a Specific m/z Value below.

In [12]:

# Find the highest peak in the average spectrum
peak_idx = np.argmax(avg_spectrum)
target_mz = mz_values[peak_idx]
print(f"Highest average peak: m/z {target_mz:.4f}")

# Extract intensities for this m/z across all pixels
if sparse.issparse(X):
    ion_intensities = np.asarray(X[:, peak_idx].toarray()).flatten()
else:
    ion_intensities = np.asarray(X[:, peak_idx]).flatten()

print(f"Intensity range: {ion_intensities.min():.2f} -- {ion_intensities.max():.2f}")

# Find the highest peak in the average spectrum peak_idx = np.argmax(avg_spectrum) target_mz = mz_values[peak_idx] print(f"Highest average peak: m/z {target_mz:.4f}") # Extract intensities for this m/z across all pixels if sparse.issparse(X): ion_intensities = np.asarray(X[:, peak_idx].toarray()).flatten() else: ion_intensities = np.asarray(X[:, peak_idx]).flatten() print(f"Intensity range: {ion_intensities.min():.2f} -- {ion_intensities.max():.2f}")

Highest average peak: m/z 760.5151
Intensity range: 0.00 -- 5468.85

In [13]:

# Reconstruct the ion image from pixel coordinates
x_unique = np.sort(np.unique(x_coords))
y_unique = np.sort(np.unique(y_coords))

ion_image = np.zeros((len(y_unique), len(x_unique)))

# Map coordinates to grid indices
x_step = x_unique[1] - x_unique[0] if len(x_unique) > 1 else 1.0
y_step = y_unique[1] - y_unique[0] if len(y_unique) > 1 else 1.0
x_indices = np.round((x_coords - x_unique.min()) / x_step).astype(int)
y_indices = np.round((y_coords - y_unique.min()) / y_step).astype(int)

for xi, yi, val in zip(x_indices, y_indices, ion_intensities):
    if 0 <= xi < len(x_unique) and 0 <= yi < len(y_unique):
        ion_image[yi, xi] = val

plt.figure(figsize=(10, 8))
plt.imshow(ion_image, cmap="hot")
plt.colorbar(label="Intensity")
plt.title(f"Ion Image -- m/z {target_mz:.4f}")
plt.xlabel("X (pixels)")
plt.ylabel("Y (pixels)")
plt.tight_layout()
plt.show()

# Reconstruct the ion image from pixel coordinates x_unique = np.sort(np.unique(x_coords)) y_unique = np.sort(np.unique(y_coords)) ion_image = np.zeros((len(y_unique), len(x_unique))) # Map coordinates to grid indices x_step = x_unique[1] - x_unique[0] if len(x_unique) > 1 else 1.0 y_step = y_unique[1] - y_unique[0] if len(y_unique) > 1 else 1.0 x_indices = np.round((x_coords - x_unique.min()) / x_step).astype(int) y_indices = np.round((y_coords - y_unique.min()) / y_step).astype(int) for xi, yi, val in zip(x_indices, y_indices, ion_intensities): if 0 <= xi < len(x_unique) and 0 <= yi < len(y_unique): ion_image[yi, xi] = val plt.figure(figsize=(10, 8)) plt.imshow(ion_image, cmap="hot") plt.colorbar(label="Intensity") plt.title(f"Ion Image -- m/z {target_mz:.4f}") plt.xlabel("X (pixels)") plt.ylabel("Y (pixels)") plt.tight_layout() plt.show()

Searching for a Specific m/z Value¶

In [14]:

# Look up a specific m/z by integrating over a tolerance window.
#
# Do NOT use a single nearest bin (np.abs(mz_values - target).argmin()):
# on a resampled axis most bins are empty by construction, so the nearest
# bin to an arbitrary target is very often zero even when the peak is
# present a few bins away.
target_mz = 760.5  # change this to your m/z of interest
mz_tolerance = 0.25

lo, hi = np.searchsorted(mz_values, [target_mz - mz_tolerance, target_mz + mz_tolerance])

if hi > lo:
    window_values = np.asarray(X[:, lo:hi].sum(axis=1)).ravel()
    populated = np.count_nonzero(np.asarray(avg_spectrum)[lo:hi])
    print(f"Window m/z {target_mz} +/- {mz_tolerance}: {hi - lo} bins, "
          f"{populated} of them populated")
    print(f"Summed intensity range: {window_values.min():.2f} -- {window_values.max():.2f}")
else:
    print(f"No m/z bins within {mz_tolerance} of {target_mz} "
          f"(axis covers {mz_values.min():.2f} -- {mz_values.max():.2f})")

# Look up a specific m/z by integrating over a tolerance window. # # Do NOT use a single nearest bin (np.abs(mz_values - target).argmin()): # on a resampled axis most bins are empty by construction, so the nearest # bin to an arbitrary target is very often zero even when the peak is # present a few bins away. target_mz = 760.5 # change this to your m/z of interest mz_tolerance = 0.25 lo, hi = np.searchsorted(mz_values, [target_mz - mz_tolerance, target_mz + mz_tolerance]) if hi > lo: window_values = np.asarray(X[:, lo:hi].sum(axis=1)).ravel() populated = np.count_nonzero(np.asarray(avg_spectrum)[lo:hi]) print(f"Window m/z {target_mz} +/- {mz_tolerance}: {hi - lo} bins, " f"{populated} of them populated") print(f"Summed intensity range: {window_values.min():.2f} -- {window_values.max():.2f}") else: print(f"No m/z bins within {mz_tolerance} of {target_mz} " f"(axis covers {mz_values.min():.2f} -- {mz_values.max():.2f})")

Window m/z 760.5 +/- 0.25: 2 bins, 2 of them populated
Summed intensity range: 1.78 -- 9976.33

Recalculating TIC and Average Spectrum¶

You can recompute TIC and average spectrum directly from the intensity matrix to verify the stored values.

In [15]:

# Recalculate TIC per pixel
if sparse.issparse(X):
    recalc_tic = np.asarray(X.sum(axis=1)).flatten()
else:
    recalc_tic = X.sum(axis=1)

print(f"Recalculated TIC range: {recalc_tic.min():.2f} -- {recalc_tic.max():.2f}")

# Recalculate TIC per pixel if sparse.issparse(X): recalc_tic = np.asarray(X.sum(axis=1)).flatten() else: recalc_tic = X.sum(axis=1) print(f"Recalculated TIC range: {recalc_tic.min():.2f} -- {recalc_tic.max():.2f}")

Recalculated TIC range: 16176.91 -- 146471.63

In [16]:

# Recalculate average spectrum (non-empty pixels only)
if sparse.issparse(X):
    nonzero_rows = np.diff(X.tocsr().indptr) > 0
    X_nonzero = X[nonzero_rows]
    recalc_avg = np.asarray(X_nonzero.mean(axis=0)).flatten()
else:
    nonzero_rows = np.any(X > 0, axis=1)
    recalc_avg = X[nonzero_rows].mean(axis=0)

# Compare with stored average
plt.figure(figsize=(12, 5))
plt.plot(mz_values, avg_spectrum, label="Stored", alpha=0.7)
plt.plot(mz_values, recalc_avg, label="Recalculated", alpha=0.7, linestyle="--")
plt.xlabel("m/z")
plt.ylabel("Intensity")
plt.title("Average Spectrum: Stored vs Recalculated")
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# Recalculate average spectrum (non-empty pixels only) if sparse.issparse(X): nonzero_rows = np.diff(X.tocsr().indptr) > 0 X_nonzero = X[nonzero_rows] recalc_avg = np.asarray(X_nonzero.mean(axis=0)).flatten() else: nonzero_rows = np.any(X > 0, axis=1) recalc_avg = X[nonzero_rows].mean(axis=0) # Compare with stored average plt.figure(figsize=(12, 5)) plt.plot(mz_values, avg_spectrum, label="Stored", alpha=0.7) plt.plot(mz_values, recalc_avg, label="Recalculated", alpha=0.7, linestyle="--") plt.xlabel("m/z") plt.ylabel("Intensity") plt.title("Average Spectrum: Stored vs Recalculated") plt.legend() plt.grid(True, alpha=0.3) plt.tight_layout() plt.show()

Working with 3D Data and Z-Slices¶

For 3D datasets, Thyra stores each z-slice as a separate table and TIC image with _z{i} suffixes.

In [17]:

# Every Thyra table carries a _z{i} suffix, so the presence of "_z" alone
# does not mean the dataset is 3D -- a plain 2D acquisition is just _z0.
# A genuine 3D dataset has more than one slice.
slice_table_keys = sorted(key for key in sdata.tables.keys() if "_z" in key)
slice_tic_keys = sorted(
    key for key in sdata.images.keys() if "_z" in key and key.endswith("_tic")
)

print(f"Slices found: {len(slice_table_keys)} -> {slice_table_keys}")
print()

if len(slice_table_keys) > 1:
    print(f"3D dataset with {len(slice_table_keys)} z-slices")
    print(f"  TIC keys: {slice_tic_keys}")

    # Access a single z-slice
    z_table = sdata.tables[slice_table_keys[0]]
    print(f"  Slice 0 shape: {z_table.shape}")

    if slice_tic_keys:
        z_tic = np.asarray(sdata.images[slice_tic_keys[0]])[0]
        plt.figure(figsize=(8, 6))
        plt.imshow(z_tic, cmap="viridis")
        plt.colorbar(label="TIC")
        plt.title(f"TIC -- {slice_tic_keys[0]}")
        plt.tight_layout()
        plt.show()
else:
    print("Single 2D dataset (one z-slice).")
    print("Convert with --handle-3d to combine a volume into one table instead.")

# Every Thyra table carries a _z{i} suffix, so the presence of "_z" alone # does not mean the dataset is 3D -- a plain 2D acquisition is just _z0. # A genuine 3D dataset has more than one slice. slice_table_keys = sorted(key for key in sdata.tables.keys() if "_z" in key) slice_tic_keys = sorted( key for key in sdata.images.keys() if "_z" in key and key.endswith("_tic") ) print(f"Slices found: {len(slice_table_keys)} -> {slice_table_keys}") print() if len(slice_table_keys) > 1: print(f"3D dataset with {len(slice_table_keys)} z-slices") print(f" TIC keys: {slice_tic_keys}") # Access a single z-slice z_table = sdata.tables[slice_table_keys[0]] print(f" Slice 0 shape: {z_table.shape}") if slice_tic_keys: z_tic = np.asarray(sdata.images[slice_tic_keys[0]])[0] plt.figure(figsize=(8, 6)) plt.imshow(z_tic, cmap="viridis") plt.colorbar(label="TIC") plt.title(f"TIC -- {slice_tic_keys[0]}") plt.tight_layout() plt.show() else: print("Single 2D dataset (one z-slice).") print("Convert with --handle-3d to combine a volume into one table instead.")

Slices found: 1 -> ['msi_dataset_z0']

Single 2D dataset (one z-slice).
Convert with --handle-3d to combine a volume into one table instead.

Pixel Size and Dataset Metadata¶

Pixel size and conversion metadata are stored in sdata.attrs (SpatialData-level attributes). Additional metadata (instrument info, acquisition params, etc.) is in the table's .uns.

In [18]:

# Pixel size from SpatialData attrs
if "pixel_size_x_um" in sdata.attrs:
    print(f"Pixel size X: {sdata.attrs['pixel_size_x_um']} um")
    print(f"Pixel size Y: {sdata.attrs['pixel_size_y_um']} um")
else:
    print("Pixel size not found in sdata.attrs")

# Dataset info
if "msi_dataset_info" in sdata.attrs:
    info = sdata.attrs["msi_dataset_info"]
    print(f"\nDataset ID: {info.get('dataset_id')}")
    print(f"Dimensions (x, y, z): {info.get('dimensions_xyz')}")
    print(f"Non-empty pixels: {info.get('non_empty_pixels'):,}")
    print(f"Total grid pixels: {info.get('total_grid_pixels'):,}")

# Table-level metadata
print(f"\nTable .uns keys: {list(msi_table.uns.keys())}")

# Pixel size from SpatialData attrs if "pixel_size_x_um" in sdata.attrs: print(f"Pixel size X: {sdata.attrs['pixel_size_x_um']} um") print(f"Pixel size Y: {sdata.attrs['pixel_size_y_um']} um") else: print("Pixel size not found in sdata.attrs") # Dataset info if "msi_dataset_info" in sdata.attrs: info = sdata.attrs["msi_dataset_info"] print(f"\nDataset ID: {info.get('dataset_id')}") print(f"Dimensions (x, y, z): {info.get('dimensions_xyz')}") print(f"Non-empty pixels: {info.get('non_empty_pixels'):,}") print(f"Total grid pixels: {info.get('total_grid_pixels'):,}") # Table-level metadata print(f"\nTable .uns keys: {list(msi_table.uns.keys())}")

Pixel size X: 25.0 um
Pixel size Y: 25.0 um

Dataset ID: msi_dataset
Dimensions (x, y, z): [48, 36, 1]
Non-empty pixels: 1,728
Total grid pixels: 1,728

Table .uns keys: ['format_specific', 'essential_metadata', 'average_spectrum', 'spatialdata_attrs', 'msi_metadata', 'raw_metadata', 'regions']

Ion Mobility¶

A Bruker timsTOF acquisition sweeps a mobility ramp at every pixel, and the MSI table above is always summed over it. The store says so, and a conversion asked for --mobility-grid also keeps the ramp:

Where	What
uns["mobility_axis"]	the axis that was summed away: the 1/K0 of every TIMS scan, its unit and the vendor calibration
uns["mobility_heatmap"]	the dataset's mean mass-mobility frame, about 4,000 m/z bins x 256 mobility channels
{table}_mobility	a sibling table with one column per occupied (m/z bin, mobility channel), written with --mobility-grid

An imzML export without a mobility array -- the synthetic dataset above -- has none of them, so the cells below open a second store. Point THYRA_TUTORIAL_TIMS_ZARR at a timsTOF conversion of your own and THYRA_TUTORIAL_MSMS_ZARR at a PASEF MS/MS one; with neither set, the repository's synthetic TIMS fixture is converted here instead.

The outputs saved in this notebook come from two real acquisitions: a 20 um MALDI-2 timsTOF measurement of a biofilm (18,838 pixels) for mobility, and a 713-pixel PASEF MS/MS acquisition of a tissue section, 15 scheduled precursors, for the tandem-MS section below. The field-by-field contract is in Output Format.

In [19]:

# With neither environment variable set, the repository's synthetic TIMS fixture
# is converted instead: that needs a checkout (the fixture lives under tests/)
# and the bundled Bruker library, which loads on Windows and Linux.
TIMS_ZARR = os.environ.get("THYRA_TUTORIAL_TIMS_ZARR")
MSMS_ZARR = os.environ.get("THYRA_TUTORIAL_MSMS_ZARR") or TIMS_ZARR
FIXTURE = next((p for p in (Path("../tests/data/fixtures/synthetic_tims.d"),
                            Path("tests/data/fixtures/synthetic_tims.d")) if p.exists()), None)

# Three precursors in disjoint slices of the fixture's 240-scan ramp, the shape a
# scheduled MALDI PASEF method has: writing them into a copy of the fixture is
# what turns it into an MS/MS acquisition.
WINDOWS = [(936.578, 49.6, 0, 80), (353.320, 35.2, 80, 160), (313.275, 34.3, 160, 240)]
PASEF_TABLE = """CREATE TABLE PasefFrameMsMsInfo (
    Frame INTEGER, ScanNumBegin INTEGER, ScanNumEnd INTEGER, IsolationMz REAL,
    IsolationWidth REAL, CollisionEnergy REAL, Precursor INTEGER)"""


def build_synthetic_tims(zarr_path):
    """Convert the test fixture, as a PASEF acquisition, with a mobility grid."""
    from thyra.convert import convert_msi

    source = zarr_path.with_name("synthetic_tims_pasef.d")
    if not source.exists():
        shutil.copytree(FIXTURE, source)
        with sqlite3.connect(source / "analysis.tdf") as db:
            db.execute("UPDATE Frames SET MsMsType = 8")
            db.execute(PASEF_TABLE)
            db.executemany(
                "INSERT INTO PasefFrameMsMsInfo VALUES (?, ?, ?, ?, 1.0, ?, NULL)",
                [(frame, begin, end, mz, energy)
                 for frame, in db.execute("SELECT Id FROM Frames").fetchall()
                 for mz, energy, begin, end in WINDOWS],
            )
    if not convert_msi(source, zarr_path, dataset_id="tims", pixel_size_um=20.0,
                       mobility_grid=True):
        shutil.rmtree(zarr_path, ignore_errors=True)
        raise RuntimeError("Conversion failed -- see the log output above.")


if TIMS_ZARR is None and FIXTURE is not None:
    TIMS_ZARR = Path("example_data/synthetic_tims.zarr")
    MSMS_ZARR = MSMS_ZARR or TIMS_ZARR
    if not TIMS_ZARR.exists():
        print("Converting the synthetic TIMS fixture...")
        build_synthetic_tims(TIMS_ZARR)

# With neither environment variable set, the repository's synthetic TIMS fixture # is converted instead: that needs a checkout (the fixture lives under tests/) # and the bundled Bruker library, which loads on Windows and Linux. TIMS_ZARR = os.environ.get("THYRA_TUTORIAL_TIMS_ZARR") MSMS_ZARR = os.environ.get("THYRA_TUTORIAL_MSMS_ZARR") or TIMS_ZARR FIXTURE = next((p for p in (Path("../tests/data/fixtures/synthetic_tims.d"), Path("tests/data/fixtures/synthetic_tims.d")) if p.exists()), None) # Three precursors in disjoint slices of the fixture's 240-scan ramp, the shape a # scheduled MALDI PASEF method has: writing them into a copy of the fixture is # what turns it into an MS/MS acquisition. WINDOWS = [(936.578, 49.6, 0, 80), (353.320, 35.2, 80, 160), (313.275, 34.3, 160, 240)] PASEF_TABLE = """CREATE TABLE PasefFrameMsMsInfo ( Frame INTEGER, ScanNumBegin INTEGER, ScanNumEnd INTEGER, IsolationMz REAL, IsolationWidth REAL, CollisionEnergy REAL, Precursor INTEGER)""" def build_synthetic_tims(zarr_path): """Convert the test fixture, as a PASEF acquisition, with a mobility grid.""" from thyra.convert import convert_msi source = zarr_path.with_name("synthetic_tims_pasef.d") if not source.exists(): shutil.copytree(FIXTURE, source) with sqlite3.connect(source / "analysis.tdf") as db: db.execute("UPDATE Frames SET MsMsType = 8") db.execute(PASEF_TABLE) db.executemany( "INSERT INTO PasefFrameMsMsInfo VALUES (?, ?, ?, ?, 1.0, ?, NULL)", [(frame, begin, end, mz, energy) for frame, in db.execute("SELECT Id FROM Frames").fetchall() for mz, energy, begin, end in WINDOWS], ) if not convert_msi(source, zarr_path, dataset_id="tims", pixel_size_um=20.0, mobility_grid=True): shutil.rmtree(zarr_path, ignore_errors=True) raise RuntimeError("Conversion failed -- see the log output above.") if TIMS_ZARR is None and FIXTURE is not None: TIMS_ZARR = Path("example_data/synthetic_tims.zarr") MSMS_ZARR = MSMS_ZARR or TIMS_ZARR if not TIMS_ZARR.exists(): print("Converting the synthetic TIMS fixture...") build_synthetic_tims(TIMS_ZARR)

In [20]:

def summed_table(store):
    """The MSI table of a store, as opposed to its sibling tables."""
    key = next(k for k in store.tables if not k.endswith(("_mobility", "_msms")))
    return key, store.tables[key]


def ion_image(table, values):
    """Lay one value per pixel out on the acquisition grid."""
    x = table.obs["x"].to_numpy().astype(int)
    y = table.obs["y"].to_numpy().astype(int)
    image = np.full((y.max() + 1, x.max() + 1), np.nan)
    image[y, x] = values
    return image


tims = sd.read_zarr(TIMS_ZARR) if TIMS_ZARR else None
if MSMS_ZARR and str(MSMS_ZARR) == str(TIMS_ZARR):
    msms = tims                                          # one store in both roles
else:
    msms = sd.read_zarr(MSMS_ZARR) if MSMS_ZARR else None
for role, store in (("mobility", tims), ("MS/MS", msms)):
    if store is None:
        print(f"{role:>8}: no store -- set the environment variable above")
        continue
    key, table = summed_table(store)
    print(f"{role:>8}: {key} {table.shape}, siblings {[k for k in store.tables if k != key]}")

def summed_table(store): """The MSI table of a store, as opposed to its sibling tables.""" key = next(k for k in store.tables if not k.endswith(("_mobility", "_msms"))) return key, store.tables[key] def ion_image(table, values): """Lay one value per pixel out on the acquisition grid.""" x = table.obs["x"].to_numpy().astype(int) y = table.obs["y"].to_numpy().astype(int) image = np.full((y.max() + 1, x.max() + 1), np.nan) image[y, x] = values return image tims = sd.read_zarr(TIMS_ZARR) if TIMS_ZARR else None if MSMS_ZARR and str(MSMS_ZARR) == str(TIMS_ZARR): msms = tims # one store in both roles else: msms = sd.read_zarr(MSMS_ZARR) if MSMS_ZARR else None for role, store in (("mobility", tims), ("MS/MS", msms)): if store is None: print(f"{role:>8}: no store -- set the environment variable above") continue key, table = summed_table(store) print(f"{role:>8}: {key} {table.shape}, siblings {[k for k in store.tables if k != key]}")

mobility: biofilm_z0 (18838, 92861), siblings ['biofilm_z0_mobility']
   MS/MS: brain1_z0 (713, 599146), siblings ['brain1_z0_mobility', 'brain1_z0_msms']

The axis that was summed away¶

uns["mobility_axis"] is the ramp itself: one 1/K0 value per TIMS scan, from the vendor calibration, decreasing with scan number. It is a feature coordinate -- it never enters obs, a coordinate system or a transform. A source with no mobility dimension writes no block at all, so "summed over mobility" and "never had any" stay distinguishable.

In [21]:

tims_key, tims_table = summed_table(tims)
axis = tims_table.uns["mobility_axis"]
ramp = np.asarray(axis["values"])

print(f"{axis['type_name']} in {axis['unit_name']}, read from {axis['source']}")
print(f"{axis['n_scans']} scans, {ramp[0]:.4f} down to {ramp[-1]:.4f},"
      f" {abs(np.diff(ramp)).mean():.5f} per scan")
low, high = np.asarray(axis["acq_range"])
print(f"declared acquisition range: {low:.4f} to {high:.4f}")
print(f"mobility-resolved sibling:  {axis.get('resolved_table', 'none written')}")

tims_key, tims_table = summed_table(tims) axis = tims_table.uns["mobility_axis"] ramp = np.asarray(axis["values"]) print(f"{axis['type_name']} in {axis['unit_name']}, read from {axis['source']}") print(f"{axis['n_scans']} scans, {ramp[0]:.4f} down to {ramp[-1]:.4f}," f" {abs(np.diff(ramp)).mean():.5f} per scan") low, high = np.asarray(axis["acq_range"]) print(f"declared acquisition range: {low:.4f} to {high:.4f}") print(f"mobility-resolved sibling: {axis.get('resolved_table', 'none written')}")

inverse reduced ion mobility in volt-second per square centimeter, read from bruker_tdf
1525 scans, 1.0658 down to 0.6656, 0.00026 per scan
declared acquisition range: 0.6451 to 1.0246
mobility-resolved sibling:  biofilm_z0_mobility

The mean mass-mobility frame¶

uns["mobility_heatmap"] is every pixel's raw (m/z, mobility, intensity) points binned and averaged over pixels. It is the discovery surface: look at it to see whether mobility separates anything in this dataset before paying for a mobility-resolved table. The 256 mobility channels are the same 256 the grid table below uses, over the same edges, so a box drawn here maps onto grid channels by integer index.

In [22]:

heat = tims_table.uns["mobility_heatmap"]
counts = np.asarray(heat["counts"])          # (m/z bins, 256), mean intensity per pixel
mz_edges = np.asarray(heat["mz_edges"])
k0_edges = np.asarray(heat["mobility_edges"])
print(f"{counts.shape[0]} m/z bins x {counts.shape[1]} mobility channels,"
      f" {(counts > 0).sum():,} of them occupied")

plt.figure(figsize=(10, 4.5), dpi=90)
plt.pcolormesh(mz_edges, k0_edges, np.log1p(counts).T, cmap="magma", shading="flat")
plt.colorbar(label="log(1 + mean intensity)")
plt.xlabel("m/z")
plt.ylabel("1/K0 (V s cm$^{-2}$)")
plt.title(f"Mass-mobility heatmap -- {tims_key}")
plt.tight_layout()
plt.show()

heat = tims_table.uns["mobility_heatmap"] counts = np.asarray(heat["counts"]) # (m/z bins, 256), mean intensity per pixel mz_edges = np.asarray(heat["mz_edges"]) k0_edges = np.asarray(heat["mobility_edges"]) print(f"{counts.shape[0]} m/z bins x {counts.shape[1]} mobility channels," f" {(counts > 0).sum():,} of them occupied") plt.figure(figsize=(10, 4.5), dpi=90) plt.pcolormesh(mz_edges, k0_edges, np.log1p(counts).T, cmap="magma", shading="flat") plt.colorbar(label="log(1 + mean intensity)") plt.xlabel("m/z") plt.ylabel("1/K0 (V s cm$^{-2}$)") plt.title(f"Mass-mobility heatmap -- {tims_key}") plt.tight_layout() plt.show()

3870 m/z bins x 256 mobility channels, 965,393 of them occupied

The heatmap and the mean spectrum are the same measurement¶

Summing the heatmap over its mobility channels must give back the average spectrum, coarsened onto mz_edges. That holds exactly under the default --tdf-spectrum scan_sum, where both are built from the same raw scans. Under --tdf-spectrum vendor_centroid it does not: the vendor peak picker keeps only the current inside the peaks it picks, while the heatmap bins every raw point (Design Decisions).

In [23]:

avg_spectrum = np.asarray(tims_table.uns["average_spectrum"])
mz_axis = tims_table.var["mz"].to_numpy()

# Coarsen the average spectrum with the heatmap's own edges, then compare.
bin_of_mz = np.clip(np.searchsorted(mz_edges, mz_axis, "right") - 1, 0, counts.shape[0] - 1)
coarsened = np.bincount(bin_of_mz, weights=avg_spectrum, minlength=counts.shape[0])
marginal = counts.sum(axis=1)

print(f"largest disagreement {np.abs(marginal - coarsened).max():.2e}"
      f" on a spectrum peaking at {coarsened.max():,.0f}"
      f" ({np.abs(marginal - coarsened).max() / coarsened.max():.1e} relative)")

plt.figure(figsize=(11, 3.5))
centres = (mz_edges[:-1] + mz_edges[1:]) / 2
plt.plot(centres, coarsened, lw=0.8, label="average_spectrum, coarsened to mz_edges")
plt.plot(centres, marginal, lw=0.8, ls="--", label="mobility_heatmap counts.sum(axis=1)")
plt.xlabel("m/z")
plt.ylabel("mean intensity")
plt.legend()
plt.tight_layout()
plt.show()

avg_spectrum = np.asarray(tims_table.uns["average_spectrum"]) mz_axis = tims_table.var["mz"].to_numpy() # Coarsen the average spectrum with the heatmap's own edges, then compare. bin_of_mz = np.clip(np.searchsorted(mz_edges, mz_axis, "right") - 1, 0, counts.shape[0] - 1) coarsened = np.bincount(bin_of_mz, weights=avg_spectrum, minlength=counts.shape[0]) marginal = counts.sum(axis=1) print(f"largest disagreement {np.abs(marginal - coarsened).max():.2e}" f" on a spectrum peaking at {coarsened.max():,.0f}" f" ({np.abs(marginal - coarsened).max() / coarsened.max():.1e} relative)") plt.figure(figsize=(11, 3.5)) centres = (mz_edges[:-1] + mz_edges[1:]) / 2 plt.plot(centres, coarsened, lw=0.8, label="average_spectrum, coarsened to mz_edges") plt.plot(centres, marginal, lw=0.8, ls="--", label="mobility_heatmap counts.sum(axis=1)") plt.xlabel("m/z") plt.ylabel("mean intensity") plt.legend() plt.tight_layout() plt.show()

largest disagreement 2.20e-03 on a spectrum peaking at 41,409 (5.3e-08 relative)

The mobility-resolved table¶

--mobility-grid bins every pixel's points onto one set of mobility channels shared across the conversion and writes {table}_mobility, where a column is an occupied (m/z bin, channel) pair. Its var["mz"] therefore repeats, so consumers must discriminate on "mobility" in var.columns and never on the element name; var["mz_index"] points back at the column of the MSI table the feature was summed into.

Two species at one m/z that separate in mobility are one column there and two here, which is the whole reason to ask for the table (Design Decisions says why it is opt in). The query below looks for such a case: among the strongest m/z bins, one whose channels hold two features at least 12 channels apart -- further apart than a TIMS peak is wide, so they are two features and not one peak's own tails.

In [24]:

grid = tims.tables[axis["resolved_table"]]
print(f"{axis['resolved_table']}: {grid.shape},"
      f" mobility in var: {'mobility' in grid.var.columns}")
print(grid.var.head(3), "\n")
print({k: v for k, v in grid.uns["mobility_grid"].items() if k != "edges"})

mz_index = grid.var["mz_index"].to_numpy()
channel = grid.var["mobility_index"].to_numpy()
mobility = grid.var["mobility"].to_numpy()
X_grid = grid.X.tocsc()
current = np.asarray(X_grid.sum(axis=0)).ravel()
pixels = np.diff(X_grid.indptr)

per_bin = np.bincount(mz_index, weights=current, minlength=tims_table.n_vars)
candidates = np.argsort(per_bin)[-300:]                     # the 300 strongest m/z bins
starts = np.searchsorted(mz_index, candidates, "left")      # var is sorted by (mz, mobility)
ends = np.searchsorted(mz_index, candidates, "right")

pair = None
for column, start, end in zip(candidates, starts, ends):
    ranked = start + np.argsort(-current[start:end])
    ranked = ranked[pixels[ranked] >= 0.2 * grid.n_obs]     # carried by more than a few pixels
    if ranked.size < 2:
        continue
    apex = ranked[0]
    far = ranked[np.abs(channel[ranked] - channel[apex]) >= 12]
    if far.size and current[far[0]] >= current[apex] / 5:   # and within 5x in current
        score = current[apex] + current[far[0]]
        if pair is None or score > pair[0]:
            pair = (score, column, apex, far[0])

if pair is None:
    print("\nNo m/z bin on this store holds two features that far apart.")
else:
    _, column, first, second = pair
    block = np.flatnonzero(mz_index == column)
    print(f"\nm/z {mz_axis[column]:.4f}: one column of {tims_key},"
          f" {block.size} in {axis['resolved_table']}")
    for feature in (first, second):
        print(f"   1/K0 {mobility[feature]:.4f}: {current[feature]:,.0f} counts"
              f" in {pixels[feature]:,} pixels")

grid = tims.tables[axis["resolved_table"]] print(f"{axis['resolved_table']}: {grid.shape}," f" mobility in var: {'mobility' in grid.var.columns}") print(grid.var.head(3), "\n") print({k: v for k, v in grid.uns["mobility_grid"].items() if k != "edges"}) mz_index = grid.var["mz_index"].to_numpy() channel = grid.var["mobility_index"].to_numpy() mobility = grid.var["mobility"].to_numpy() X_grid = grid.X.tocsc() current = np.asarray(X_grid.sum(axis=0)).ravel() pixels = np.diff(X_grid.indptr) per_bin = np.bincount(mz_index, weights=current, minlength=tims_table.n_vars) candidates = np.argsort(per_bin)[-300:] # the 300 strongest m/z bins starts = np.searchsorted(mz_index, candidates, "left") # var is sorted by (mz, mobility) ends = np.searchsorted(mz_index, candidates, "right") pair = None for column, start, end in zip(candidates, starts, ends): ranked = start + np.argsort(-current[start:end]) ranked = ranked[pixels[ranked] >= 0.2 * grid.n_obs] # carried by more than a few pixels if ranked.size < 2: continue apex = ranked[0] far = ranked[np.abs(channel[ranked] - channel[apex]) >= 12] if far.size and current[far[0]] >= current[apex] / 5: # and within 5x in current score = current[apex] + current[far[0]] if pair is None or score > pair[0]: pair = (score, column, apex, far[0]) if pair is None: print("\nNo m/z bin on this store holds two features that far apart.") else: _, column, first, second = pair block = np.flatnonzero(mz_index == column) print(f"\nm/z {mz_axis[column]:.4f}: one column of {tims_key}," f" {block.size} in {axis['resolved_table']}") for feature in (first, second): print(f" 1/K0 {mobility[feature]:.4f}: {current[feature]:,.0f} counts" f" in {pixels[feature]:,} pixels")

biofilm_z0_mobility: (18838, 13502098), mobility in var: True
                mz  mobility  mz_index  mobility_index
mz1_im1  220.00165  0.667901         1               1
mz1_im3  220.00165  0.671028         1               3
mz1_im5  220.00165  0.674155         1               5 

{'channel_width': 0.0015633925461622541, 'law': 'linear', 'upper': 1.0657848887235417, 'lower': 0.6655563969060047, 'n_channels': 256}

m/z 272.1968: one column of biofilm_z0, 251 in biofilm_z0_mobility
   1/K0 0.8618: 14,554,084 counts in 15,712 pixels
   1/K0 0.8430: 3,250,400 counts in 15,417 pixels

In [25]:

if pair is None:
    print("Nothing to plot on this store.")
else:
    order = block[np.argsort(mobility[block])]
    plt.figure(figsize=(11, 3))
    plt.plot(mobility[order], current[order], lw=1)
    for feature in (first, second):
        plt.axvline(mobility[feature], color="tab:red", ls="--", lw=1)
    plt.xlabel("1/K0 (V s cm$^{-2}$)")
    plt.ylabel("total intensity")
    plt.title(f"the {block.size} channels of m/z {mz_axis[column]:.4f}")
    plt.tight_layout()
    plt.show()

if pair is None: print("Nothing to plot on this store.") else: order = block[np.argsort(mobility[block])] plt.figure(figsize=(11, 3)) plt.plot(mobility[order], current[order], lw=1) for feature in (first, second): plt.axvline(mobility[feature], color="tab:red", ls="--", lw=1) plt.xlabel("1/K0 (V s cm$^{-2}$)") plt.ylabel("total intensity") plt.title(f"the {block.size} channels of m/z {mz_axis[column]:.4f}") plt.tight_layout() plt.show()

A feature is a mobility peak, not a single channel, so each is imaged over the channels around its apex -- the same reason the m/z lookup above integrates a window rather than trusting one bin. Summing all the bin's channels instead gives back the merged column exactly, per pixel: that is the marginal invariant uns["mobility_marginal"] records.

In [26]:

if pair is None:
    print("Nothing to image on this store.")
else:
    summed_column = np.asarray(tims_table.X.tocsc()[:, column].todense()).ravel()
    all_channels = np.asarray(X_grid[:, block].sum(axis=1)).ravel()
    print(f"every channel of the bin adds back up to the merged column:"
          f" {np.abs(all_channels - summed_column).max():.1e} largest difference")

    half = max(1, abs(int(channel[first]) - int(channel[second])) // 3)
    images = [(f"summed m/z {mz_axis[column]:.4f}", summed_column)]
    for feature in (first, second):
        window = block[np.abs(channel[block] - channel[feature]) <= half]
        values = np.asarray(X_grid[:, window].sum(axis=1)).ravel()
        images.append((f"1/K0 {mobility[feature]:.4f} +/- {half} channels", values))
    print(f"the two features image at r ="
          f" {np.corrcoef(images[1][1], images[2][1])[0, 1]:.3f}")

    fig, axes = plt.subplots(1, 3, figsize=(9, 3.2), dpi=72)
    for ax, (title, values) in zip(axes, images):
        image = ion_image(tims_table, values)
        ax.imshow(image, cmap="inferno", vmax=np.nanpercentile(image, 99))
        ax.set_title(title, fontsize=10)
        ax.axis("off")
    plt.tight_layout()
    plt.show()

if pair is None: print("Nothing to image on this store.") else: summed_column = np.asarray(tims_table.X.tocsc()[:, column].todense()).ravel() all_channels = np.asarray(X_grid[:, block].sum(axis=1)).ravel() print(f"every channel of the bin adds back up to the merged column:" f" {np.abs(all_channels - summed_column).max():.1e} largest difference") half = max(1, abs(int(channel[first]) - int(channel[second])) // 3) images = [(f"summed m/z {mz_axis[column]:.4f}", summed_column)] for feature in (first, second): window = block[np.abs(channel[block] - channel[feature]) <= half] values = np.asarray(X_grid[:, window].sum(axis=1)).ravel() images.append((f"1/K0 {mobility[feature]:.4f} +/- {half} channels", values)) print(f"the two features image at r =" f" {np.corrcoef(images[1][1], images[2][1])[0, 1]:.3f}") fig, axes = plt.subplots(1, 3, figsize=(9, 3.2), dpi=72) for ax, (title, values) in zip(axes, images): image = ion_image(tims_table, values) ax.imshow(image, cmap="inferno", vmax=np.nanpercentile(image, 99)) ax.set_title(title, fontsize=10) ax.axis("off") plt.tight_layout() plt.show()

every channel of the bin adds back up to the merged column: 0.0e+00 largest difference
the two features image at r = 0.593

Tandem MS¶

An MS/MS imaging run measures fragments, so var["mz"] is fragment m/z and nothing about the axis says so. uns["msms_schedule"] is what says it: the MS level, the isolation windows, their collision energies, and -- on Bruker PASEF, where each precursor owns a slice of the mobility ramp -- the scan range each window occupied.

When a pixel's frame isolates several precursors, the stored spectrum sums all of them: it is a chimera and must not be read as any one precursor's fragment spectrum. Thyra says so at WARNING during conversion and, by default, also writes the split apart as {table}_msms (Design Decisions).

In [27]:

msms_key, msms_table = summed_table(msms)
schedule = msms_table.uns["msms_schedule"]
msms_ramp = np.asarray(msms_table.uns["mobility_axis"]["values"])

print(f"MS level {schedule['ms_level']}, {schedule['n_windows']} precursors per pixel,"
      f" merged into one spectrum: {schedule['merges_precursors']},"
      f" same schedule at every pixel: {schedule['constant_across_pixels']}")
print(f"demultiplexed sibling: {schedule.get('resolved_table', 'none written')}\n")

target = np.asarray(schedule["isolation_window_target"])
begin = np.asarray(schedule["scan_begin"])
end = np.asarray(schedule["scan_end"])
print(" precursor       isolation window      energy      scans          1/K0")
for i in np.argsort(target):
    low = target[i] - np.asarray(schedule["isolation_window_lower_offset"])[i]
    high = target[i] + np.asarray(schedule["isolation_window_upper_offset"])[i]
    print(f"  {target[i]:9.3f}   {low:8.3f} - {high:<8.3f}"
          f" {np.asarray(schedule['collision_energy'])[i]:5.1f} eV"
          f"  {begin[i]:5d} - {end[i]:<5d}"
          f"  {msms_ramp[end[i] - 1]:.3f} - {msms_ramp[begin[i]]:.3f}")

msms_key, msms_table = summed_table(msms) schedule = msms_table.uns["msms_schedule"] msms_ramp = np.asarray(msms_table.uns["mobility_axis"]["values"]) print(f"MS level {schedule['ms_level']}, {schedule['n_windows']} precursors per pixel," f" merged into one spectrum: {schedule['merges_precursors']}," f" same schedule at every pixel: {schedule['constant_across_pixels']}") print(f"demultiplexed sibling: {schedule.get('resolved_table', 'none written')}\n") target = np.asarray(schedule["isolation_window_target"]) begin = np.asarray(schedule["scan_begin"]) end = np.asarray(schedule["scan_end"]) print(" precursor isolation window energy scans 1/K0") for i in np.argsort(target): low = target[i] - np.asarray(schedule["isolation_window_lower_offset"])[i] high = target[i] + np.asarray(schedule["isolation_window_upper_offset"])[i] print(f" {target[i]:9.3f} {low:8.3f} - {high:<8.3f}" f" {np.asarray(schedule['collision_energy'])[i]:5.1f} eV" f" {begin[i]:5d} - {end[i]:<5d}" f" {msms_ramp[end[i] - 1]:.3f} - {msms_ramp[begin[i]]:.3f}")

MS level 2, 15 precursors per pixel, merged into one spectrum: True, same schedule at every pixel: True
demultiplexed sibling: brain1_z0_msms

 precursor       isolation window      energy      scans          1/K0
    313.275    312.775 - 313.775   34.3 eV   2933 - 3015   0.889 - 0.916
    353.320    352.820 - 353.820   35.2 eV   2832 - 2901   0.927 - 0.949
    367.335    366.835 - 367.835   36.2 eV   2678 - 2798   0.961 - 0.999
    393.300    392.800 - 393.800   40.3 eV   2099 - 2351   1.105 - 1.185
    456.395    455.895 - 456.895   38.1 eV   2463 - 2539   1.045 - 1.069
    521.361    520.861 - 521.861   38.9 eV   2382 - 2429   1.080 - 1.095
    615.533    615.033 - 616.033   43.6 eV   1783 - 1855   1.262 - 1.285
    621.486    620.986 - 621.986   42.5 eV   1889 - 2004   1.215 - 1.251
    701.510    701.010 - 702.010   44.4 eV   1683 - 1747   1.296 - 1.316
    723.490    722.990 - 723.990   45.2 eV   1582 - 1647   1.327 - 1.347
    751.524    751.024 - 752.024   46.0 eV   1485 - 1549   1.358 - 1.377
    784.577    784.077 - 785.077   46.8 eV   1387 - 1450   1.388 - 1.407
    857.527    857.027 - 858.027   47.5 eV   1292 - 1353   1.418 - 1.436
    923.500    923.000 - 924.000   48.4 eV   1177 - 1259   1.447 - 1.471
    936.578    936.078 - 937.078   49.6 eV    994 - 1127   1.487 - 1.527

One precursor's block¶

Each precursor's columns are a contiguous block, sorted by (precursor_mz, precursor_mobility, mz), and the block's row sums are its precursor ion image. Two precursors can legitimately share an m/z -- a method isolates one mass at two mobility positions to keep an isomer pair apart -- so the block identity is precursor_index, never a float comparison on precursor_mz, and precursor_index means nothing outside this one store: two samples align on (precursor_mz, precursor_mobility) (Output Format).

A column inside the block is one fragment's image. In the outputs saved here the strongest of them lies inside the isolation window itself -- an ion that reached the collision cell and survived it, which is ordinary in a fragment spectrum -- and the fragments below the window are far weaker.

In [28]:

demux = msms.tables[schedule["resolved_table"]]
var = demux.var
print(f"{schedule['resolved_table']}: {demux.shape},"
      f" precursor_mz in var: {'precursor_mz' in var.columns}")
print(var.head(3), "\n")

X_msms = demux.X.tocsc()
precursor_index = var["precursor_index"].to_numpy()
fragment_current = np.asarray(X_msms.sum(axis=0)).ravel()
block_current = np.bincount(precursor_index, weights=fragment_current)
block_mz = var["precursor_mz"].to_numpy()[np.searchsorted(precursor_index,
                                                          np.arange(block_current.size))]

# Look the precursor up by nearest value, then slice on its block index.
wanted = float(block_mz[np.argmax(block_current)])   # the one carrying the most current
index = int(precursor_index[np.argmin(np.abs(var["precursor_mz"].to_numpy() - wanted))])
block_columns = np.flatnonzero(precursor_index == index)

print(f"precursor {var['precursor_mz'].iloc[block_columns[0]]:.3f}"
      f" at 1/K0 {var['precursor_mobility'].iloc[block_columns[0]]:.4f}:"
      f" block {index}, {block_columns.size:,} fragment columns,"
      f" {block_current[index] / block_current.sum() * 100:.0f}% of the table's current")
print(f"the {block_current.size} blocks carry"
      f" {block_current.min():,.0f} to {block_current.max():,.0f} counts")

demux = msms.tables[schedule["resolved_table"]] var = demux.var print(f"{schedule['resolved_table']}: {demux.shape}," f" precursor_mz in var: {'precursor_mz' in var.columns}") print(var.head(3), "\n") X_msms = demux.X.tocsc() precursor_index = var["precursor_index"].to_numpy() fragment_current = np.asarray(X_msms.sum(axis=0)).ravel() block_current = np.bincount(precursor_index, weights=fragment_current) block_mz = var["precursor_mz"].to_numpy()[np.searchsorted(precursor_index, np.arange(block_current.size))] # Look the precursor up by nearest value, then slice on its block index. wanted = float(block_mz[np.argmax(block_current)]) # the one carrying the most current index = int(precursor_index[np.argmin(np.abs(var["precursor_mz"].to_numpy() - wanted))]) block_columns = np.flatnonzero(precursor_index == index) print(f"precursor {var['precursor_mz'].iloc[block_columns[0]]:.3f}" f" at 1/K0 {var['precursor_mobility'].iloc[block_columns[0]]:.4f}:" f" block {index}, {block_columns.size:,} fragment columns," f" {block_current[index] / block_current.sum() * 100:.0f}% of the table's current") print(f"the {block_current.size} blocks carry" f" {block_current.min():,.0f} to {block_current.max():,.0f} counts")

brain1_z0_msms: (713, 264006), precursor_mz in var: True
               precursor_mz         mz  precursor_index  mz_index  \
p313.275_mz5        313.275  50.001375                0         5   
p313.275_mz18       313.275  50.004625                0        18   
p313.275_mz40       313.275  50.010126                0        40   

               precursor_mobility  
p313.275_mz5             0.902907  
p313.275_mz18            0.902907  
p313.275_mz40            0.902907   

precursor 751.524 at 1/K0 1.3674: block 10, 13,535 fragment columns, 14% of the table's current
the 15 blocks carry 391,808 to 1,562,162 counts

In [29]:

precursor_image = np.asarray(X_msms[:, block_columns].sum(axis=1)).ravel()

# One fragment, integrated over a 0.05 window of the borrowed MS1 axis: on a
# resampled axis a peak covers several bins, exactly as for an MS1 ion image.
fragment_mz = var["mz"].to_numpy()[block_columns]
strongest = fragment_mz[np.argmax(fragment_current[block_columns])]
window = block_columns[np.abs(fragment_mz - strongest) <= 0.05]
fragment_image = np.asarray(X_msms[:, window].sum(axis=1)).ravel()

fig, axes = plt.subplots(1, 2, figsize=(9, 4.5))
for ax, values, title in ((axes[0], precursor_image, f"precursor {wanted:.3f}, all fragments"),
                          (axes[1], fragment_image, f"fragment m/z {strongest:.4f}")):
    image = ion_image(msms_table, values)
    ax.imshow(image, cmap="inferno", vmax=np.nanpercentile(image, 99))
    ax.set_title(title, fontsize=10)
    ax.axis("off")
plt.tight_layout()
plt.show()

# The split is a filter on the scan number, so the blocks add back up: under the
# default scan_sum every point of a frame falls in exactly one isolation window.
conserved = demux.uns["demultiplexed_current"]
print(f"the split holds {conserved['current_ratio']:.4f}x the current of"
      f" {conserved['summed_table']}, per pixel"
      f" {conserved['current_ratio_pixel_min']:.4f} to"
      f" {conserved['current_ratio_pixel_max']:.4f}")

precursor_image = np.asarray(X_msms[:, block_columns].sum(axis=1)).ravel() # One fragment, integrated over a 0.05 window of the borrowed MS1 axis: on a # resampled axis a peak covers several bins, exactly as for an MS1 ion image. fragment_mz = var["mz"].to_numpy()[block_columns] strongest = fragment_mz[np.argmax(fragment_current[block_columns])] window = block_columns[np.abs(fragment_mz - strongest) <= 0.05] fragment_image = np.asarray(X_msms[:, window].sum(axis=1)).ravel() fig, axes = plt.subplots(1, 2, figsize=(9, 4.5)) for ax, values, title in ((axes[0], precursor_image, f"precursor {wanted:.3f}, all fragments"), (axes[1], fragment_image, f"fragment m/z {strongest:.4f}")): image = ion_image(msms_table, values) ax.imshow(image, cmap="inferno", vmax=np.nanpercentile(image, 99)) ax.set_title(title, fontsize=10) ax.axis("off") plt.tight_layout() plt.show() # The split is a filter on the scan number, so the blocks add back up: under the # default scan_sum every point of a frame falls in exactly one isolation window. conserved = demux.uns["demultiplexed_current"] print(f"the split holds {conserved['current_ratio']:.4f}x the current of" f" {conserved['summed_table']}, per pixel" f" {conserved['current_ratio_pixel_min']:.4f} to" f" {conserved['current_ratio_pixel_max']:.4f}")

the split holds 1.0000x the current of brain1_z0, per pixel 1.0000 to 1.0000

Which defaults produced this store¶

The reading that produced the numbers above is recorded, not implied: processing names the spectrum the TDF reader took, and the versioned metadata block names both sibling tables. Every entry is argued on the Design Decisions page.

In [30]:

for store in ([tims] if msms is tims else [tims, msms]):
    key, table = summed_table(store)
    meta = table.uns["msi_metadata"]
    analysis = meta["ms_analysis"]
    conversion = json.loads(meta["processing"])[0]          # a JSON string in uns
    print(f"{key} -- schema {meta['schema_version']},"
          f" thyra {conversion['software']['version']}")
    print(f"   tdf_spectrum                   {conversion['parameters'].get('tdf_spectrum', '-')}")
    print(f"   ion_mobility.resolved_table    "
          f"{analysis.get('ion_mobility', {}).get('resolved_table', '-')}")
    print(f"   fragmentation.resolved_table   "
          f"{analysis.get('fragmentation', {}).get('resolved_table', '-')}")

for store in ([tims] if msms is tims else [tims, msms]): key, table = summed_table(store) meta = table.uns["msi_metadata"] analysis = meta["ms_analysis"] conversion = json.loads(meta["processing"])[0] # a JSON string in uns print(f"{key} -- schema {meta['schema_version']}," f" thyra {conversion['software']['version']}") print(f" tdf_spectrum {conversion['parameters'].get('tdf_spectrum', '-')}") print(f" ion_mobility.resolved_table " f"{analysis.get('ion_mobility', {}).get('resolved_table', '-')}") print(f" fragmentation.resolved_table " f"{analysis.get('fragmentation', {}).get('resolved_table', '-')}")

biofilm_z0 -- schema 0.5.0, thyra 3.19.0
   tdf_spectrum                   scan_sum
   ion_mobility.resolved_table    biofilm_z0_mobility
   fragmentation.resolved_table   -
brain1_z0 -- schema 0.5.0, thyra 3.19.0
   tdf_spectrum                   scan_sum
   ion_mobility.resolved_table    brain1_z0_mobility
   fragmentation.resolved_table   brain1_z0_msms

Implementation Notes¶

Average spectrum is calculated using only non-empty pixels (those with at least one non-zero intensity). This prevents artificial dilution from empty grid positions.

TIC image orientation follows standard imaging convention: first axis is Y (rows, top to bottom), second is X (columns, left to right).

Sparse storage: MSI data is highly sparse. Keep data in sparse format when possible and avoid .toarray() on the full matrix for large datasets.

Loading: Use spatialdata.read_zarr() to load datasets. This returns a SpatialData object with lazy-loaded elements.