Imaging Mass Cytometry Pipeline

Pipeline Overview

Raw MCD/TIFF Files ──> Image Processing ──> Cell Masks │ ▼ ┌─────────────────────────────────────────────┐ │ imc-pipeline │ ├─────────────────────────────────────────────┤ │ 1. Data Preprocessing (spillover, hot px) │ │ 2. Cell Segmentation (Cellpose/Mesmer) │ │ 3. Single-cell Quantification │ │ 4. Clustering & Phenotyping │ │ 5. Spatial Analysis │ │ 6. Visualization │ └─────────────────────────────────────────────┘ │ ▼ Cell Types + Spatial Neighborhoods

Complete steinbock Workflow

Step 1: Setup and Preprocessing

Initialize steinbock project

steinbock preprocess imc
--mcd data/*.mcd
--panel panel.csv
--output raw/

Hot pixel filtering

steinbock preprocess imc hotpixel
--input raw/
--output img/
--threshold 50

Create nuclear and membrane channels

steinbock preprocess mosaic
--input img/
--channels panel.csv
--output mosaics/

Step 2: Cell Segmentation

Using Cellpose

steinbock segment cellpose
--input img/
--panel panel.csv
--channel DNA1 DNA2
--output masks/
--diameter 20

Alternative: Using Mesmer

steinbock segment mesmer
--input img/
--panel panel.csv
--nuclear DNA1 DNA2
--membrane CD45
--output masks/

Step 3: Single-cell Quantification

Extract intensities

steinbock measure intensities
--input img/
--masks masks/
--panel panel.csv
--output intensities/

Measure cell properties (area, etc.)

steinbock measure regionprops
--masks masks/
--output regionprops/

Extract neighbor relationships

steinbock measure neighbors
--masks masks/
--output neighbors/
--distance 15

Complete Python Workflow

import pandas as pd import numpy as np import anndata as ad import scanpy as sc import squidpy as sq from pathlib import Path

=== 1. LOAD DATA ===

data_dir = Path('steinbock_output')

intensities = pd.read_csv(data_dir / 'intensities.csv', index_col=0) regionprops = pd.read_csv(data_dir / 'regionprops.csv', index_col=0) neighbors = pd.read_csv(data_dir / 'neighbors.csv')

print(f'Loaded {len(intensities)} cells')

=== 2. CREATE ANNDATA ===

adata = ad.AnnData(X=intensities.values, obs=regionprops, var=pd.DataFrame(index=intensities.columns)) adata.obs['image_id'] = [idx.split('_')[0] for idx in intensities.index] adata.obs['cell_id'] = intensities.index

Add spatial coordinates

adata.obsm['spatial'] = regionprops[['centroid_y', 'centroid_x']].values

=== 3. PREPROCESSING ===

Arcsinh transform (cofactor 5 for IMC)

adata.X = np.arcsinh(adata.X / 5)

Scale for clustering

sc.pp.scale(adata, max_value=10) adata.raw = adata.copy()

=== 4. DIMENSIONALITY REDUCTION ===

sc.pp.pca(adata, n_comps=20) sc.pp.neighbors(adata, n_neighbors=15) sc.tl.umap(adata)

=== 5. CLUSTERING ===

sc.tl.leiden(adata, resolution=0.8) print(f'Found {adata.obs["leiden"].nunique()} clusters')

=== 6. PHENOTYPING ===

Marker expression per cluster

sc.tl.rank_genes_groups(adata, 'leiden', method='wilcoxon') marker_genes = sc.get.rank_genes_groups_df(adata, group=None)

Annotate clusters based on markers

cluster_annotations = { '0': 'T cells', '1': 'Macrophages', '2': 'Tumor', '3': 'B cells', '4': 'Stromal' } adata.obs['cell_type'] = adata.obs['leiden'].map(cluster_annotations)

=== 7. SPATIAL ANALYSIS ===

Build spatial graph

sq.gr.spatial_neighbors(adata, coord_type='generic', delaunay=True)

Neighborhood enrichment

sq.gr.nhood_enrichment(adata, cluster_key='cell_type')

Co-occurrence analysis

sq.gr.co_occurrence(adata, cluster_key='cell_type')

Ripley's statistics

sq.gr.ripley(adata, cluster_key='cell_type', mode='L')

=== 8. VISUALIZATION ===

import matplotlib.pyplot as plt

UMAP by cell type

fig, axes = plt.subplots(1, 2, figsize=(14, 5)) sc.pl.umap(adata, color='cell_type', ax=axes[0], show=False) sc.pl.umap(adata, color='leiden', ax=axes[1], show=False) plt.savefig('umap_celltypes.png', dpi=150, bbox_inches='tight')

Spatial plot

fig, ax = plt.subplots(figsize=(10, 10)) sq.pl.spatial_scatter(adata[adata.obs['image_id'] == 'image1'], color='cell_type', shape=None, size=10, ax=ax) plt.savefig('spatial_celltypes.png', dpi=150, bbox_inches='tight')

Neighborhood enrichment heatmap

sq.pl.nhood_enrichment(adata, cluster_key='cell_type') plt.savefig('neighborhood_enrichment.png', dpi=150, bbox_inches='tight')

=== 9. DIFFERENTIAL ANALYSIS ===

Compare conditions

adata.obs['condition'] = adata.obs['image_id'].map({ 'image1': 'Control', 'image2': 'Control', 'image3': 'Treatment', 'image4': 'Treatment' })

Cell type proportions

proportions = adata.obs.groupby(['image_id', 'condition', 'cell_type']).size().unstack(fill_value=0) proportions = proportions.div(proportions.sum(axis=1), axis=0)

Save results

adata.write('imc_analysis.h5ad') proportions.to_csv('cell_type_proportions.csv') print('Analysis complete!')

R Alternative (imcRtools)

library(imcRtools) library(cytomapper) library(CATALYST)

Read steinbock output

spe <- read_steinbock('steinbock_output/')

Transform

assay(spe, 'exprs') <- asinh(counts(spe) / 5)

Cluster

spe <- runDR(spe, features = rownames(spe), exprs_values = 'exprs', dr = 'UMAP') spe <- cluster(spe, features = rownames(spe), exprs_values = 'exprs', xdim = 10, ydim = 10, maxK = 20)

Spatial analysis

spe <- buildSpatialGraph(spe, img_id = 'image_id', type = 'expansion', threshold = 20) spe <- aggregateNeighbors(spe, colPairName = 'neighborhood', by = 'cluster_id')

Spatial context

cn <- detectCommunity(spe, colPairName = 'neighborhood', size_threshold = 10, group_by = 'image_id')

Plot

plotSpatial(spe, img_id = 'image1', node_color_by = 'cluster_id')

QC Checkpoints

Stage Check Action if Failed

Preprocessing No hot pixel streaks Lower threshold

Segmentation

80% cells detected Adjust diameter

Quantification All markers extracted Check panel.csv

Clustering 5-20 clusters Adjust resolution

Spatial Neighbors detected Check distance

Workflow Variants

High-plex Panels (40+ markers)

Use batch-aware clustering

import scvi

scvi.model.SCVI.setup_anndata(adata, batch_key='image_id') model = scvi.model.SCVI(adata) model.train() adata.obsm['X_scvi'] = model.get_latent_representation() sc.pp.neighbors(adata, use_rep='X_scvi')

Tumor Microenvironment Analysis

Spatial interactions with tumor

tumor_cells = adata[adata.obs['cell_type'] == 'Tumor'].obs_names sq.gr.ligrec(adata, cluster_key='cell_type', source_groups=['Tumor'], target_groups=['T cells', 'Macrophages'])

Related Skills

• imaging-mass-cytometry/data-preprocessing - Hot pixel, spillover

• imaging-mass-cytometry/cell-segmentation - Cellpose/Mesmer details

• imaging-mass-cytometry/phenotyping - Cluster annotation

• imaging-mass-cytometry/spatial-analysis - Spatial statistics

• imaging-mass-cytometry/interactive-annotation - Manual cell labeling

• imaging-mass-cytometry/quality-metrics - QC metrics

• single-cell/clustering - Clustering methods

• spatial-transcriptomics/spatial-statistics - Related spatial methods