Bio imaging mass cytometry differential analysis

Version Compatibility

Reference examples tested with: diffcyt 1.22+ (R), lme4 1.1+ (R), statsmodels 0.14+, scanpy 1.10+, sccoda 0.1.9+, numpy 1.26+, pandas 2.2+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures
• R: packageVersion('<pkg>') then ?function_name to verify parameters

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Notes specific to this skill: cell-type proportions are compositional (they sum to 1), so a real increase in one type forces apparent decreases in others -- a Dirichlet/CLR-aware method (scCODA) or a reference cell type is needed, not independent per-type tests. statsmodels mixedlm fits patient as a random effect. diffcyt operates on per-sample cluster counts (DA) and per-sample median marker expression (DS).

IMC Differential Analysis

"Compare cell types and spatial structure between my conditions" -> Aggregate to the patient, then test across patients -- never across cells.

• Python: statsmodels.formula.api.mixedlm, sccoda, scanpy
• R: diffcyt, lme4::lmer, SpaceANOVA for spatial differential testing

The Single Most Important Modern Insight -- the replicate is the patient, not the cell, and the million-cell count is a red herring

After phenotyping and spatial analysis, every interesting claim ("disease has more Tregs", "responders have more CD8-tumor contact") is a comparison BETWEEN groups, and the single most common fatal error is testing it at the cell level. Hundreds of thousands of cells from one patient are not independent replicates -- they are correlated reads of one biological sample, and cells within an image are massively spatially autocorrelated. Testing at the cell level inflates n by orders of magnitude and manufactures significance: a per-cell test over 50,000 cells reports p~0 for trivial effects because the effective sample size is the number of PATIENTS (often 10-40), not cells (Squair 2021 Nat Commun 12:5692). In imaging this is worse than in scRNA because slide and ROI add nesting levels and ROIs are not random samples of the tissue. The correct spine is invariant across every differential question: compute a per-image (or per-ROI) summary, aggregate to ONE value per patient, then test across patients with the patient as the unit -- a mixed model with patient as a random effect (image nested within patient), a pseudobulk-style per-patient summary, or a cell-count-weighted average of per-image statistics (Samorodnitsky and Wu 2024 Brief Bioinform 25:bbae522). Two riders complete the picture: cell-type proportions are COMPOSITIONAL (they are constrained to sum to 1, so a real rise in one type mechanically depresses the others, and independent per-type tests double-count this), and acquisition BATCH drifts by day/run and can align with clinical group, so batch must be a covariate and acquisition order randomized against condition. Phenotyping-method choice and statistical-unit choice are orthogonal: getting the cell types right does not excuse testing them wrong.

Differential Question Taxonomy

Question	Per-image summary	Patient-level test
Cell-type abundance differs between groups	per-image cell-type proportions	mixed model on proportions; scCODA (compositional); diffcyt-DA
A functional/state marker differs within a type	per-image median marker per type	pseudobulk per patient + limma/edgeR; diffcyt-DS
A spatial interaction/niche differs between groups	per-image enrichment z / Ripley's K / CN abundance	mixed model / cell-count-weighted; SpaceANOVA (FANOVA on cross-K)

Decision Tree by Scenario

Scenario	Recommended	Why
Compare cell-type proportions, several types shift	scCODA (or CLR + mixed model)	handles the compositional constraint and the reference-type problem
Standard cytometry-style DA on clusters	diffcyt-DA (edgeR on per-sample counts)	designed for per-sample cluster counts; established
Multiple ROIs per patient	mixed model, patient random effect (image nested)	respects nesting; or cell-count-weighted aggregate
Few patients (n < ~10)	simple per-patient test; report low power honestly	do NOT rescue power with cell count
Differential functional-marker expression within a type	pseudobulk per patient + limma/edgeR (diffcyt-DS)	aggregates out cell-level pseudoreplication
Spatial interaction/niche across groups	per-image spatial stat -> patient aggregate; SpaceANOVA	the spatial statistic is the summary; the unit is still the patient
Acquisition batch aligns with group	include batch covariate; if confounded, the contrast is unrescuable	randomize acquisition order against condition

Aggregate to the Patient (the spine)

Goal: Collapse millions of cells to one summary per patient before any test.

Approach: Compute per-image cell-type proportions, then aggregate images to their patient. Every downstream test consumes this patient-level table, not the cell table.

import pandas as pd

# obs has one row per cell with image_id, patient, condition, cell_type
counts = obs.groupby(['patient', 'condition', 'image_id', 'cell_type']).size().unstack(fill_value=0)
image_prop = counts.div(counts.sum(axis=1), axis=0)            # per-image proportions
patient_prop = image_prop.groupby(['patient', 'condition']).mean()   # one row per patient

Differential Abundance with a Mixed Model

Goal: Test a cell type's proportion across groups while respecting patient/ROI nesting.

Approach: Fit a mixed model with patient as a random effect when multiple ROIs per patient exist; this absorbs within-patient correlation that a fixed-effect test would treat as independent replication.

import statsmodels.formula.api as smf

# one row per image; proportion of the target type; patient random intercept
df = image_prop.reset_index().rename(columns={'Treg': 'prop'})   # 'Treg' = an actual cell_type column
model = smf.mixedlm('prop ~ condition + batch', df, groups=df['patient'])   # batch as covariate
res = model.fit()
print(res.summary())   # the condition coefficient is tested with patient as the unit

Compositional Differential Abundance (scCODA)

Goal: Avoid the false "everything changed" artifact when proportions are constrained to sum to 1.

Approach: Model the counts as compositional against a reference cell type; a change is interpreted relative to that reference rather than as an independent per-type shift.

import sccoda.util.cell_composition_data as dat
from sccoda.util import comp_ana as mod

# patient-level cell-type COUNTS (not proportions); pick a biologically stable reference type
data = dat.from_pandas(patient_counts, covariate_columns=['condition'])
analysis = mod.CompositionalAnalysis(data, formula='condition', reference_cell_type='Epithelial')
result = analysis.sample_hmc()
result.summary()

Differential Spatial Feature

Goal: Test whether a spatial interaction or niche differs between groups, at the patient unit.

Approach: Treat the per-image spatial statistic (a neighborhood-enrichment z, a Ripley's cross-K curve, a CN abundance) as the summary, aggregate to patient, and test across patients with FDR over cell-type pairs. SpaceANOVA does this as a functional ANOVA on per-image cross-K with subject structure.

import statsmodels.formula.api as smf
from statsmodels.stats.multitest import multipletests

# per_image_pair: rows = (image_id, patient, condition, batch, enrichment_z) for ONE type pair
pvals = {}
for pair, sub in per_image_pair.groupby('pair'):
    res = smf.mixedlm('enrichment_z ~ condition + batch', sub, groups=sub['patient']).fit()
    pvals[pair] = res.pvalues['condition[T.responder]']
padj = dict(zip(pvals, multipletests(list(pvals.values()), method='fdr_bh')[1]))   # FDR across pairs

Differential State Within a Type

Goal: Test whether a functional/state marker (Ki67, PD-1) differs within a cell type between groups.

Approach: Pseudobulk to one value per patient per cell type (median marker expression among that type's cells), then test across patients. diffcyt-DS (R) formalizes this on per-sample medians; the Python pseudobulk path is below.

import numpy as np

t = adata[adata.obs['cell_type'] == 'T cell']
ki67 = t[:, 'Ki67'].X
ki67 = ki67.toarray().ravel() if hasattr(ki67, 'toarray') else np.asarray(ki67).ravel()
pb = t.obs.assign(ki67=ki67).groupby(['patient', 'condition'])['ki67'].median().reset_index()
print(smf.ols('ki67 ~ condition', pb).fit().pvalues['condition[T.responder]'])   # one value per patient

Per-Method Failure Modes

Cell-level testing

Trigger: a t-test/Wilcoxon/regression over individual cells. Mechanism: correlated cells from one sample are pseudoreplicates; effective n is the patient count. Symptom: p~0 for trivial effects; "significant" findings that do not replicate. Fix: aggregate to per-patient summaries; test across patients.

Independent per-type proportion tests

Trigger: a separate test per cell type on proportions. Mechanism: proportions sum to 1, so a real rise in one type depresses others mechanically. Symptom: many types appear to change in opposite directions. Fix: compositional model (scCODA) or CLR transform with a reference type.

Over-correcting batch

Trigger: aggressive integration to make clusters patient-agnostic, then testing on corrected data. Mechanism: correction can treat real between-patient biology as batch. Symptom: the disease signal disappears. Fix: correct minimally, validate that invariant types align while variable types stay separate, and keep integration out of the across-patient inference path.

Rescuing power with cell count

Trigger: claiming significance from n=4 patients because millions of cells were imaged. Mechanism: cell count is not replication. Symptom: confident claims from few patients. Fix: report the patient n and the true power honestly; collect more patients.

Quantitative Thresholds

Threshold	Source	Rationale
Unit of replication = patient count (often 10-40)	Squair 2021 Nat Commun 12:5692	cells/ROIs are pseudoreplicates
Cell-count-weighted per-image aggregation	Samorodnitsky and Wu 2024 Brief Bioinform 25:bbae522	controls type-I error with high power (vs unweighted ROI averaging)
Reference cell type for compositional DA	Buttner 2021 Nat Commun 12:6876	proportions are not independent
BH-FDR across cell-type pairs and radii	multiplicity	~200 pairs x radii guarantees false positives
Batch covariate; randomize acquisition order	spatial dossier	run drift can align with clinical group

Common Errors

Error / symptom	Cause	Solution
p~0 with a trivial effect	cell-level test	per-patient aggregation; mixed model
All cell types "changed"	compositional constraint ignored	scCODA / CLR with reference type
Disease signal vanished after integration	batch over-correction	minimal correction; keep it out of the test
Significant with n=4 patients	power rescued by cell count	report patient n; do not over-claim
Many significant pairs	no FDR across pairs/radii	BH-FDR over the full grid

References

• Squair JW, Gautier M, Kathe C, et al. 2021. Confronting false discoveries in single-cell differential expression. Nat Commun 12:5692. — pseudoreplication; aggregate to sample.
• Samorodnitsky S, Wu MC. 2024. Statistical analysis of multiple regions-of-interest in multiplexed spatial proteomics data. Brief Bioinform 25(6):bbae522. — ROI aggregation, cell-count-weighted averaging, SPOT omnibus.
• Seal S, Neelon B, Angel PM, et al. 2024. SpaceANOVA: Spatial Co-occurrence Analysis of Cell Types in Multiplex Imaging Data Using Point Process and Functional ANOVA. J Proteome Res 23(4):1131-1143. — FANOVA on per-image cross-K with subject structure.
• Weber LM, Nowicka M, Soneson C, Robinson MD. 2019. diffcyt: Differential discovery in high-dimensional cytometry via high-resolution clustering. Commun Biol 2:183. — diffcyt-DA/DS.
• Buttner M, Ostner J, Muller CL, Theis FJ, Schubert B. 2021. scCODA is a Bayesian model for compositional single-cell data analysis. Nat Commun 12:6876. — compositional differential abundance.
• Schurch CM, Bhate SS, Barlow GL, et al. 2020. Coordinated Cellular Neighborhoods Orchestrate Antitumoral Immunity at the Colorectal Cancer Invasive Front. Cell 182(5):1341-1359.e19. — cellular neighborhoods compared across patients.

Related Skills

• phenotyping - supplies the cell-type labels whose proportions are compared
• spatial-analysis - supplies the per-image spatial statistics that become patient-level summaries
• quality-metrics - batch must be diagnosed and entered as a covariate
• experimental-design/randomization-blocking - the experimental-unit and pseudoreplication foundation
• clinical-biostatistics/subgroup-analysis - multiplicity and effect estimation in clinical cohorts
• flow-cytometry/differential-analysis - diffcyt-DA/DS for suspension cytometry