bio-imaging-mass-cytometry-differential-analysis

What this skill does

## 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.

```python
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.

```python
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.

```python
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.

```python
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