bio-workflows-cytometry-pipeline

What this skill does

## Version Compatibility

Reference examples tested with: CATALYST 1.26+, diffcyt 1.22+, FlowSOM 2.10+, flowCore 2.14+, flowWorkspace 4.14+, flowStats 4.14+, edgeR 4.0+, limma 3.58+, ggplot2 3.5+; Python (partial alt) flowkit 1.1+.

Before using code patterns, verify installed versions match. If versions differ:
- R: `packageVersion('<pkg>')` then `?function_name` to verify parameters
- Python: `pip show <package>` then `help(module.function)` to check signatures

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt rather than retrying. Each stage defers depth to its component skill.

# Flow Cytometry Pipeline

**"Process my cytometry data from FCS to differential populations"** -> read raw -> compensate/unmix -> transform -> QC -> remove doublets -> cluster (or gate) -> annotate -> test DA/DS, with the sample as the unit of inference.
- R: `flowCore` + `CATALYST::prepData/cluster/runDR` + `diffcyt::diffcyt()`

## The Single Most Important Modern Insight -- A Pipeline Is a Chain of Irreversible Decisions, and the Unit of Inference Is the Sample

Each early choice silently gates the validity of the final test: reading raw (not log-linearized), compensating BEFORE transforming, removing margin events before density QC, assigning type-vs-state markers correctly, and removing doublets before clustering. None of these is recoverable downstream - a doublet clustered as a "double-positive," a state marker used for clustering, or an uncompensated channel becomes a false population that the differential test then "confirms." The second load-bearing thread is that the SAMPLE/subject, not the cell, is the experimental unit: diffcyt aggregates cells to per-sample-per-cluster counts (DA) and medians (DS) before testing, so biological replication (>= 2-3 per group) is mandatory and a per-cell test is invalid. Two normalization layers sit at OPPOSITE ends of the pipeline - EQ-bead drift correction on raw counts at the very front (CyTOF), and CytoNorm cross-batch harmonization after clustering - and conflating them is a classic error.

## Decision Tree: Which Path

| Situation | Path | Why |
|-----------|------|-----|
| Conventional fluorescence flow | compensate (`$SPILLOVER`/flowStats) -> logicle -> ... | optical spillover; logicle handles negatives |
| Spectral cytometer (Aurora/ID7000) | UNMIX (not compensate) -> arcsinh ~150 | overdetermined system; fluorescence-scale |
| Mass cytometry (CyTOF) | EQ-bead normalize (raw) -> arcsinh cofactor 5 -> `compCytof` if needed | metals barely spill (~1-4%); drift correction first |
| High-dim discovery, no prior gates | cluster (FlowSOM via CATALYST) | scales; finds unexpected populations |
| Well-defined populations / rare events (MRD) | hierarchical gating (openCyto) | interpretable; clustering fails for ultra-rare |
| Multi-batch / multi-day | anchor sample per batch -> CytoNorm (after clustering) | model batch in the design for inference |

## Pipeline Overview

```
FCS -> compensate/unmix -> transform -> QC (margins, time, dead) -> doublets
     -> [ cluster (FlowSOM) | gate (openCyto) ] -> annotate -> diffcyt DA/DS -> report
EQ-bead drift normalization (CyTOF) runs on raw counts BEFORE everything; CytoNorm runs LAST.
```

## 1. Panel, Metadata, and Load

**Goal:** Define the type/state panel and sample metadata, then load FCS.

**Approach:** Panel `marker_class` drives everything downstream (type clusters, state is tested); metadata keys samples to condition/subject. See flow-cytometry/fcs-handling.

```r
library(CATALYST); library(diffcyt); library(flowCore); library(ggplot2)

panel <- data.frame(
  fcs_colname = c('FSC-A','SSC-A','CD45','CD3','CD4','CD8','CD19','CD14','Ki67','IFNg'),
  antigen     = c('FSC','SSC','CD45','CD3','CD4','CD8','CD19','CD14','Ki67','IFNg'),
  marker_class = c('none','none','type','type','type','type','type','type','state','state'))
md <- data.frame(file_name = list.files('data', pattern = '\\.fcs$'),
                 sample_id = paste0('S', 1:8),
                 condition = rep(c('Control','Treatment'), each = 4),
                 patient_id = rep(paste0('P', 1:4), 2))
fs <- read.flowSet(file.path('data', md$file_name), transformation = FALSE, truncate_max_range = FALSE)
```

## 2. Compensate / Unmix, then Transform

**Goal:** Remove spillover on linear data, then variance-stabilize.

**Approach:** Conventional flow compensates (matrix before transform); CyTOF skips fluorescence compensation and uses cofactor 5; spectral unmixes then uses ~150. See flow-cytometry/compensation-transformation.

```r
fs_comp <- compensate(fs, spillover(fs[[1]])[[1]])      # conventional flow; CyTOF: omit or use compCytof
COFACTOR <- 150                                          # 5 for CyTOF, ~150 for fluorescence/spectral
sce <- prepData(fs_comp, panel, md, transform = TRUE, cofactor = COFACTOR, FACS = TRUE)
```

## 3. QC (order matters)

**Goal:** Remove margin/boundary events and time anomalies before any density step.

**Approach:** Margins first, then time-based cleaning; on CyTOF, EQ-bead drift correction happens upstream on raw counts. See flow-cytometry/cytometry-qc and flow-cytometry/bead-normalization.

```r
# per-sample sanity + sample-similarity MDS (flag outlier samples)
plotExprs(sce, color_by = 'condition'); plotMDS(sce, color_by = 'condition')
# event-level cleaning runs per-FCS upstream: PeacoQC::RemoveMargins() -> PeacoQC()/flowAI on transformed data
```

## 4. Remove Doublets

**Goal:** Drop aggregates before clustering so they don't form phantom double-positives.

**Approach:** Flow uses the FSC-A vs FSC-H diagonal; CyTOF uses DNA intercalator + Gaussian/Event_length. See flow-cytometry/doublet-detection.

```r
# CyTOF (FACS=TRUE retained Event_length on the arcsinh scale):
e <- assay(sce, 'exprs')
if (all(c('DNA1','Event_length') %in% rownames(sce))) {
  keep <- e['DNA1', ] > quantile(e['DNA1', ], 0.05) &
          e['Event_length', ] <= quantile(e['Event_length', ], 0.99)
  sce <- sce[, keep]
}
```

## 5. Cluster (FlowSOM) or Gate

**Goal:** Define populations by unsupervised clustering on TYPE markers (discovery) or hierarchical gating (defined/rare).

**Approach:** `cluster()` wraps FlowSOM+ConsensusClusterPlus; over-provision the grid, set a seed. See flow-cytometry/clustering-phenotyping (clustering) and flow-cytometry/gating-analysis (gating).

```r
sce <- cluster(sce, features = 'type', xdim = 10, ydim = 10, maxK = 20, seed = 42)
```

## 6. Annotate and Visualize Structure

**Goal:** Label metaclusters from marker medians; embed for display only.

**Approach:** Median heatmap drives annotation; UMAP colors by cluster but is never used to define or quantify populations.

```r
plotExprHeatmap(sce, features = 'type', by = 'cluster_id', k = 'meta20', scale = 'last')
sce <- runDR(sce, dr = 'UMAP', features = 'type', cells = 2000)
plotDR(sce, dr = 'UMAP', color_by = 'meta20')
```

## 7. Differential Abundance and State

**Goal:** Test which populations change in frequency (DA) or state-marker expression (DS) between conditions.

**Approach:** The `diffcyt()` wrapper aggregates to the sample level; results live in `res$res`. See flow-cytometry/differential-analysis.

```r
design   <- createDesignMatrix(ei(sce), cols_design = 'condition')
contrast <- createContrast(c(0, 1))                        # Treatment vs Control
res_DA <- diffcyt(sce, clustering_to_use = 'meta20', analysis_type = 'DA',
                  method_DA = 'diffcyt-DA-edgeR', design = design, contrast = contrast)
res_DS <- diffcyt(sce, clustering_to_use = 'meta20', analysis_type = 'DS',
                  method_DS = 'diffcyt-DS-limma', design = design, contrast = contrast)
da <- as.data.frame(SummarizedExperiment::rowData(res_DA$res))   # cluster_id, logFC, p_val, p_adj
```

## 8. Visualize Results and Export

**Goal:** Summarize significant populations and persist results.

**Approach:** Pass the inner result object (`res$res`) to plotting; export tables and the SCE.

`