fcsumap - gated multi-sample flow-cytometry UMAP

A single R script that takes a folder of FCS files, applies a declarative gating strategy, and produces one shared UMAP embedding across all samples, plus the figures and tables needed to interpret and QC it.

Nothing about any particular study is compiled into the script. Panels, markers, populations, thresholds, cohorts and covariates are all discovered from the data or declared in a YAML config, so the same script runs on a different panel without editing code.

Quick commands

One folder, one batch, one marker panel - the minimum viable run:

Rscript scripts/fcsumap.R --dir /path/to/data --outdir results

docker run --rm -v "/path/to/data:/data:ro" -v "$PWD/results:/results" \
  fcs-umap:1.0 --dir /data --outdir /results

Everything past this point in this section is the same run scaled up - more folders, more subfolders, more marker panels in one go. See §3 for what each flag does and §4 for the Docker-specific env vars (PIPELINE_SCRIPT, --memory).

Adding samples to an existing study? Go straight to §3.2 - which files to edit, what breaks a new batch, and what to check before trusting the numbers. For a routine batch nothing under scripts/ changes.

Two defaults worth knowing up front:

• Other CD45+ is hidden from the UMAP figures. It is the leftover, not a population - every CD45+ cell matching no definition in the spec - and it is routinely large enough to blanket the real populations and dominate every legend. Pass --other to draw it. Tables always include it either way.
• Study and sex colours are fixed, not palette positions: the --reference-group is green, the other studies red then purple; male green, female red. See §5.4.

1. What it does, and the one design decision that matters

Conventional practice embeds each sample separately, or embeds gated populations separately. Both make between-sample comparison meaningless: UMAP axes carry no absolute meaning, so two independently-computed embeddings cannot be overlaid, and a cell's position tells you nothing about a cell in the other plot.

This script computes one embedding per panel, from a balanced subsample of CD45+ events pooled across every sample, and then colours that single embedding by population, sample, cohort and covariates. Because all cells share one coordinate system, "cohort A's monocytes sit where cohort B's do not" is a claim the figure can actually support.

Two consequences worth understanding before you read any output:

• Gates are computed per sample, on that sample's own distributions (a valley between modes, not a fixed number), because instrument gain and staining drift between acquisitions. The embedding is shared; the gating is not.
• Populations are labels on cells, not inputs to the embedding. The UMAP is built from marker intensities alone. If a labelled population lands in one compact island, that is independent evidence the gate is sane - it was never told where to go.

2. Requirements

R ≥ 4.4, and:

Source	Packages
Bioconductor	flowCore (required), flowDensity (optional - better threshold finding)
CRAN	uwot, data.table, ggplot2, patchwork, viridis, scales, yaml, optparse, dplyr, tidyr, hexbin
CRAN, optional	readxl - only needed if --absolute-counts points at an .xlsx file (§5.5); .csv/.tsv need nothing extra

Install them all with:

Rscript scripts/install_deps.R

The installer reads the package list out of fcsumap.R itself rather than repeating it, so it cannot fall out of step with what the script loads. A missing optional package is reported as a note; a missing required one exits non-zero.

Memory. Peak usage is driven by events held in memory, not by file size on disk. Use --max-events-per-file to bound it (see §6).

3. Quick start

# 1. Discover files and write a config with every derived threshold in it
Rscript scripts/fcsumap.R --dir /path/to/data --recursive \
        --write-config my_config.yaml

# 2. Read that config. Edit any threshold you disagree with.

# 3. Run
Rscript scripts/fcsumap.R \
        --dir /path/to/data --recursive \
        --config my_config.yaml \
        --outdir results

That works with no metadata at all - you get populations, frequencies and a UMAP coloured by population and sample. Patient covariates require the two optional metadata files in §5.

3.1 Multiple folders, batches, panels, and sample groups

The examples above already scale past one tidy folder of FCS files - nothing below is a different mode, just the same flags applied to messier input.

One folder, cohorts in subfolders - --recursive searches down:

Rscript scripts/fcsumap.R --dir /path/to/study --recursive \
        --config my_config.yaml --outdir results

Several folders that aren't nested under one parent (separate drives, separate acquisition batches shipped as separate archives, cohorts that live next to unrelated data you don't want scanned) - --files takes a comma-separated list of paths or globs, and can be combined with --dir:

Rscript scripts/fcsumap.R \
        --files "/data/batch_2024_01/*,/data/batch_2024_06/*,/mnt/backup/extra_cohort/*" \
        --config my_config.yaml --outdir results

A folder per sample GROUP, kept as separate runs - useful for a quick look at one arm, or to compare a scoped rerun against the full cohort (§8.1 explains --include-qc-failed; the same reasoning applies to any subset run):

for grp in healthy_controls arm_a arm_b; do
  Rscript scripts/fcsumap.R --dir "/data/$grp" \
          --config my_config.yaml --outdir "results_$grp"
done

Point every run's --outdir at its own folder if you keep more than one - every run writes the same filenames (population_frequencies.png etc.) regardless of which subset produced them.

More than one marker panel, or samples from more than one acquisition batch with a different marker set, in the SAME --dir/--files - no flag needed. Panels are fingerprinted from each file's own marker symbols and handled entirely separately: one shared embedding per panel, and once there is more than one panel every output filename gets a _<panel-name> suffix (population_frequencies_panel_1.png, population_frequencies_panel_2.png, ...) so populations scored under different marker sets are never pooled into one figure. --umap-markers/--umap-markers-all (§6) override which markers drive the embedding, if the built-in lineage-marker preference list is a poor fit for your panel.

As many sample groups as the patient table declares - --group-column (default: the resolved cohort column) drives every between-group figure and the density/UMAP-by-group comparisons (§8). It is not limited to two or three groups; group_comparison.png, umap_density_by_group.png and umap_overview_by_group.png all render one column per group, however many there are.

3.2 Adding a new batch of samples

Start with one question: is the marker panel the same?

	What you do
Same panel (same antibodies, same fluorochromes)	Drop the files in, add rows to the two metadata CSVs, re-run. Nothing in scripts/ changes. Gates re-derive per sample automatically.
Different panel	Nothing extra either - panels are fingerprinted from each file's own marker symbols and embedded separately, with a _<panel-name> suffix on every output. You only touch scripts/config_cohorts.yaml if the new panel needs population definitions the current one doesn't have.

Thresholds are not copied from the last run. They are computed per sample from that sample's own distributions, so a new batch acquired at different gain gets its own gates without you editing anything.

The files to touch

File	When	What to add
data/<any folder>/*.fcs	Always	The new acquisitions. Subfolders are fine - --recursive searches down.
data/sample_map_*.csv	Always	One row per new FCS file: file,sample_id,patient_id,timepoint,is_control
data/patient_table_*.csv	If the batch brings new patients	One row per new patient: patient_id,cohort,age_years,sex
data/absolute cell counts.xlsx	Only if you have haemogram counts	One row per new sample (§5.5). Without it, figures fall back to frequencies.
scripts/config_cohorts.yaml	Only if populations, ratios or colours change	Edit the populations: / ratios: / colors: blocks
Your run command or driver script	Only if a flag changes	e.g. a different --reference-group

Nothing under scripts/ needs editing for a routine batch. If you find yourself changing the R script to onboard data, that is a bug - the panel, populations, cohorts and covariates are all meant to come from the data or the config.

Five things that will actually bite you

1. The file column must match the filename on disk byte for byte. Including a trailing copy, and including oddities like ;2 in Blood samples_HC-13;2.fcs. A mismatch aborts the run before any analysis with "these input files are not in the sample map". This is the single most common way a new batch fails, and it is usually a file that got renamed after the map was written.
2. --pattern has to match the new names. The default study command uses '[.]fcs( copy)?$', which accepts .fcs and .fcs copy and nothing else.
3. patient_id must agree across the two CSVs. Matching ignores case and surrounding whitespace, but an unmatched ID is not an error - you get a WARNING: no patient-table row for patient_id: line and silently NA covariates. Read the log.
4. --reference-group must equal a value in the cohort column exactly, spelling and spacing included. It anchors the statistics and the green colour; if it matches nothing, the reference silently falls back to the first group alphabetically.
5. --memory is coupled to --max-events-per-file. Raise the event cap without raising the memory cap and the container is OOM-killed mid-run: exit 137, no R error, no traceback, the log just stops. It reads as a hang, not a crash.

Check these four things after the run

Check	Where
Every new file was read and none silently dropped	Log: INPUT: N FCS file(s) should equal your file count
Which samples were excluded, and why	staining_qc.csv
Gates landed in the right place on the new batch	recon_diagnostics.png, then gating_qc.png
Population sizes are biologically plausible	gate_counts.csv, population_frequencies.csv

Do not read any comparison figure before those four. A gate that drifted on a new batch produces confident, well-formatted, wrong numbers everywhere downstream.

If the new batch needs its own gating strategy

Regenerate a config from the data, edit it, then run with it:

# 1. Derive thresholds from the new batch and write them out
Rscript scripts/fcsumap.R --dir /path/to/new_batch --recursive \
        --write-config new_config.yaml

# 2. Open new_config.yaml. Compare its thresholds against gating_qc.png from a
#    previous run. Edit any you disagree with - the file is data, not code.

# 3. Run with it
Rscript scripts/fcsumap.R --dir /path/to/new_batch --recursive \
        --config new_config.yaml --outdir results_new_batch

Point each run's --outdir at its own folder. Every run writes the same filenames regardless of which subset produced them, so a shared --outdir overwrites the previous batch's results in place.

4. Running in Docker

The pinned image reproduces the exact environment the pipeline was validated in (R 4.4.3, flowCore 2.18.0, flowDensity 1.40.0, uwot 0.2.4). This matters more than it usually does: uwot is stochastic and flowDensity places gate thresholds, so a different version of either changes numbers you might publish, even with the same --seed.

# Build from the repository ROOT (not from scripts/)
docker build -f scripts/Dockerfile -t fcs-umap:1.0 .

# Run: data mounted read-only, results mounted writable
docker run --rm \
  -v "$PWD:/data:ro" \
  -v "$PWD/results:/results" \
  -e OMP_NUM_THREADS=8 \
  fcs-umap:1.0 \
  --dir /data --recursive \
  --config /data/scripts/my_config.yaml \
  --outdir /results

On Windows PowerShell use ${PWD} instead of $PWD.

4.1 Multiple folders, batches, panels, and sample groups in Docker

Same patterns as §3.1 - the container only ever sees what's mounted, so the one thing that differs from the plain Rscript case is getting the right host paths into /data.

Cohorts in subfolders under one parent - mount the parent once, --recursive searches it same as bare Rscript:

docker run --rm -v "/path/to/study:/data:ro" -v "$PWD/results:/results" \
  fcs-umap:1.0 --dir /data --recursive \
  --config /data/scripts/my_config.yaml --outdir /results

Folders that share a common parent but aren't meant to all be scanned - mount the parent once and use --files with subpath globs, exactly as in §3.1:

docker run --rm -v "/data:/data:ro" -v "$PWD/results:/results" fcs-umap:1.0 \
  --files "/data/batch_2024_01/*,/data/batch_2024_06/*" \
  --config /data/scripts/my_config.yaml --outdir /results

Folders that do NOT share a common parent (a second drive, a separate mount point) - mount each one into its own container path, then reference both in --files:

docker run --rm \
  -v "/data/batch_2024_01:/data1:ro" \
  -v "/mnt/backup/extra_cohort:/data2:ro" \
  -v "$PWD/results:/results" \
  fcs-umap:1.0 --files "/data1/*,/data2/*" \
  --config /data/scripts/my_config.yaml --outdir /results

A folder per sample group, kept as separate runs - give each its own results mount so outputs don't overwrite each other:

foreach ($grp in "healthy_controls","arm_a","arm_b") {
  docker run --rm -v "${PWD}\data\$grp:/data:ro" -v "${PWD}\results_$grp:/results" `
    fcs-umap:1.0 --dir /data --config /data/scripts/my_config.yaml --outdir /results
}

Multiple marker panels or acquisition batches with different marker sets, in the same mount - no flag needed, same as bare Rscript: panels are auto-detected and every output past the first panel gets a _<panel-name> suffix. As many sample groups as the patient table declares - also no Docker-specific change; --group-column and the between-group figures (§8) scale the same way regardless of how the data was mounted.

Editing the script without rebuilding the image

The command above runs the copy of fcsumap.R baked into the image at build time. Add one env var and it runs the mounted copy instead - the same /data mount already used for --config/--sample-map/--patient-table, just pointed at the script too:

docker run --rm \
  -v "$PWD:/data:ro" \
  -v "$PWD/results:/results" \
  -e OMP_NUM_THREADS=8 \
  -e PIPELINE_SCRIPT=/data/scripts/fcsumap.R \
  fcs-umap:1.0 \
  --dir /data --recursive \
  --config /data/scripts/my_config.yaml \
  --outdir /results

The image and the script are then fully independent: routine edits to fcsumap.R take effect on the very next run, no rebuild. A rebuild is still required when the script starts needing a package the image doesn't have

• nothing else. Set PIPELINE_SCRIPT on every run and this is simply how the container behaves. See scripts/entrypoint.sh.

Notes:

• No data is ever copied into the image. .dockerignore excludes every FCS file, cohort directory, spreadsheet and generated patient table. Data is mounted at run time. This keeps the build fast and keeps patient-linked files out of any image that might be pushed to a registry.
• The image is verified at build time - it runs --help and prints resolved package versions - so a broken image fails during docker build rather than an hour into an analysis. That check runs against the image's own baked-in copy of the script, since at build time there is nothing else to run it against.
• Give the container enough memory (--memory=16g) or cap events per file.
• The entrypoint resolves which script to run at container start (baked-in copy, or $PIPELINE_SCRIPT when set) - either way, every flag in §6 works unchanged.

5. Input files

5.1 FCS files (required)

Point --dir at a folder; add --recursive if cohorts live in subdirectories.

The default --exclude drops single-stain compensation controls by filename, because those are instrument-setup acquisitions, not samples - embedding them would place a cloud of single-colour events in the shared space. Pass --exclude '' to keep everything.

If a file on disk fails to match --pattern, the script warns and names it rather than proceeding silently. Analysing an unnoticed subset of a cohort is a worse outcome than a failed run. Archives with duplicated names need --pattern '[.]fcs( copy)?$'.

5.2 --sample-map (optional)

Maps each file to a sample identity. Generate a template with --write-sample-map map.csv; see scripts/examples/sample_map_template.csv.

Column	Meaning
file	filename (the only required column)
sample_id	short label used in figures
patient_id	join key to the patient table
timepoint	e.g. baseline / follow-up
is_control	TRUE for unstained/reference tubes
panel	force a panel assignment instead of fingerprinting

is_control is load-bearing and asymmetric. TRUE excludes the tube from the embedding and lets its distributions serve as the negative reference for unimodal markers. FALSE is a positive assertion that this is a biological sample: if it then fails staining QC, it is excluded but never promoted to reference - one bad tube's failure must not propagate into every threshold in the panel.

5.3 --patient-table (optional)

Per-patient covariates, joined on patient_id. German or English column names are both accepted and output is always English - Geschlecht→sex, Geburtsdatum→date_of_birth, Alter→age_years, Groesse→height_cm, Gewicht→weight_kg, and values (weiblich→female) likewise. Already-English input passes through unchanged rather than warning.

Ages are derived from date of birth against --reference-date. Set this to a fixed study date; left at today's date, ages change between runs and results stop being reproducible.

Any covariate column with more than one distinct value is discovered and plotted automatically - including a cohort or group column, which is usually the central variable of the study.

5.4 --config (optional but recommended)

YAML holding the asinh cofactor, per-marker thresholds and the population definitions. Run with --write-config first, inspect what was derived, then edit and re-run with --config.

Populations are declared per marker as above / below:

populations:
  Classical monocytes:
    CD3: below
    CD14: above
    CD16: below

  # any_of is a disjunction: CD16-positive OR CD56-positive
  NK cells:
    CD3: below
    any_of:
      CD16: above
      CD56: above

  # Scatter channels are gateable, for populations defined by physics rather
  # than by an antibody (granulocytes: dense granules scatter strongly)
  Granulocytes:
    CD45: above
    CD3: below
    CD16: above
    SSC-A: above

Scatter channels can define populations but are excluded from UMAP features - the embedding is built from fluorescence only.

Standard immunology ratios (e.g. CD4:CD8) are declared the same way, as a ratios: block referencing population names from populations: above - there is no built-in default, since a hardcoded population name would be silently wrong on a different panel:

ratios:
  CD4_CD8:
    label: "CD4:CD8 ratio"      # optional; defaults to the block's key
    numerator: CD4 T cells
    denominator: CD8 T cells

Each entry writes one row per sample to population_ratios.csv and one panel to population_ratios.png (§8). A ratio whose numerator or denominator population was not scored in any sample is skipped, with a note in the log.

Colours. Every colour every figure draws with lives in one place - default_colors() in the script - and is entirely re-themeable through the config's colors: block, without editing R. Any key you omit keeps its built-in default; see scripts/examples/colors_template.yaml for the full list (pasting it in unedited changes nothing - it's the actual defaults) or run --write-config to get the currently active set written back out alongside everything else that run derived:

colors:
  population_palette: ["#E8112D", "#FF7A00", "#FFC800"]  # discrete: populations/samples
  study_palette: ["#00A651", "#E8112D", "#8E44E8"]       # reference first, then studies
  sex_palette: {male: "#00A651", female: "#E8112D"}
  other_grey: "#7A5C00"          # the "Other CD45+" catch-all
  gate_highlight: "#D62728"      # recon_diagnostics.png: gate box + CD45 cutoff line
  density_viridis: mako          # hexbin/density figures (viridisLite option name/letter)
  intensity_viridis: magma       # marker-intensity continuous UMAP colouring

Two variables carry fixed meaning across every figure. A colour that means one group in one panel and a different group in the next is worse than no colour at all, because the reader carries the association across the page whether or not it holds. So study and sex are pinned before any palette is issued:

Variable	Assignment
Study / cohort	--reference-group is always green; the remaining studies take red then purple, in sorted order
Sex	male green, female red

Every discrete scale in the pipeline checks these first (semantic_colours()) and falls through to population_palette for anything else - so populations, samples and any other categorical keep the ordinary palette.

population_palette is bright, one clear hue per entry, walking the wheel red → orange → yellow → lime → green → cyan → blue → pink. No black, violet or grey. Verified minimum pairwise CIE Lab distance 28.2 at full opacity, well clear of the 12–18 band where colours start being confused. Violet is deliberately absent, which is why nine is the practical ceiling here rather than twelve - beyond that pop_palette() falls back to hcl.colors().

Key	Used in
population_palette	Discrete population/sample colouring that is not study or sex
study_palette	Study/cohort curves and bars, reference first (green)
sex_palette	Every figure that colours by sex
other_grey	"Other CD45+" catch-all (kept out of population_palette on purpose; hidden from UMAPs unless --other)
reference_fill, bar_outline, bracket	group_comparison.png and everything that shares its layout (population_frequencies.png, functional_markers.png, population_ratios.png, absolute_counts.png)
gate_highlight	recon_diagnostics.png: leukocyte gate box, CD45 cutoff line
threshold_ok, threshold_review	gating_qc.png: valley-derived vs. fallback threshold
tile_border, na_fill, count_viridis	absolute_counts_qc.png heatmap
density_viridis	Hexbin/density figures: recon_diagnostics.png, umap_density*.png
intensity_viridis	Marker-intensity / numeric-covariate continuous UMAP colouring
grid_major, axis_ticks, subtitle_text, caption_text, label_text, empty_panel_text	Shared theme chrome and figure-level captions, throughout

5.5 --absolute-counts (optional)

A directly measured absolute cell count per population - from dual-platform counting, counting beads, or a clinical flow lab's own analysis - as opposed to the wbc_per_ul column in §5.3, which derives an absolute count by multiplying this pipeline's own gated percentages by a single whole-blood WBC number. The two are independent: this one does not depend on this pipeline's gates agreeing with whoever produced the counts.

Accepts .xlsx, .csv, or .tsv; every format is normalized to absolute_counts_raw.csv in --outdir before parsing (one code path regardless of source - see that file if a match or a number looks wrong). An .xlsx input needs the readxl package (Rscript install_deps.R).

Format - wide, one row per sample, one column per population, matching how these are typically exported from Excel or a clinical flow report. See scripts/examples/absolute_counts_template.csv:

Row/column	Meaning
Row 1, column 1	ignored (a label like Sample:)
Row 1, columns 2+	population names - free text, need not match --config's population names
Column 1, rows 2+	sample identifier: matched against --sample-map by either patient_id or the acquisition filename (case/whitespace-insensitive)
A column with a header but no numeric value in any row	treated as a units label (e.g. in 1 mL blood) and skipped as data; its text is scanned for mL/µL to auto-convert per-mL counts to the cells_per_ul convention - assumed already cells/µL, with a logged NOTE, if no unit is found
A row with no sample identifier	tolerated and skipped (cohort block separators)

Requires --sample-map (it supplies the join key). Produces, in this order:

• absolute_counts_qc.png - one tile per matched sample x population, log-scale colour. Inspect this first: it is externally supplied data this pipeline did not measure, so a transcription slip or a stray order-of-magnitude error will not look anomalous in a single boxplot the way it jumps out in a full grid. Controls and QC-failed samples are labelled, not hidden, for the same reason.
• absolute_counts.csv - the matched long-format table.
• absolute_counts.png / absolute_counts_stats.csv - the same grouped bar/point figure and Wilcoxon tests as group_comparison.png (§8.1), with --group-column.

Rows that match no sample, or that match a sample outside the current --dir/--files subset, are dropped with a NOTE naming them rather than silently guessed at or included.

6. Command-line reference

Flag	Default	Notes
--dir, --files	-	input; --files takes paths or a quoted glob
--recursive	off	search --dir recursively
--pattern	[.]fcs$	filename regex
--exclude	compensation controls	regex applied after discovery; '' keeps all
--max-events-per-file	0 (all)	events per file, sampled evenly through the acquisition
--cells-per-sample	20000	balanced draw per sample for the embedding
--max-cells	200000	ceiling per embedding
--config, --write-config	-	read / derive-and-exit
--sample-map, --write-sample-map	-	read / template-and-exit
--patient-table	-	covariates CSV
--absolute-counts	-	directly measured absolute counts (.xlsx/.csv/.tsv, §5.5); requires --sample-map
--reference-date	today	set explicitly for reproducible ages
--group-column	resolved cohort	patient-table column defining the comparison groups
--reference-group	first alphabetically	group all others are tested against; drawn unfilled and leftmost
--p-adjust-display	raw	raw or BH - which p-value the figure brackets show; both always in the CSV
--cofactor	derived per panel	asinh cofactor
--n-neighbors, --min-dist	30, 0.3	UMAP structure
--singlet-mad-k	3	singlet band width in MADs
--min-cd45-pct	5	staining QC floor
--include-qc-failed	off	force declared samples (not controls) that fail staining QC into gating/UMAP/frequencies anyway. Their pct_of_cd45_pos carries no real evidence (§8) - meant for small cohorts where excluding any sample is too costly, not for routine use
--umap-markers	built-in lineage list	comma-separated markers preferred as UMAP features; override when the default list has fewer than 2 matches in your panel
--umap-markers-all	off	use every eligible marker as a UMAP feature instead of preferring lineage markers
--viability-marker	auto	dye name override
--threads	0 (all)	UMAP threads
--seed	42	RNG seed
--keep-exprs, --no-session	off	session-state size control
--outdir	results	output directory
--other	off	also draw the Other CD45+ catch-all on the UMAP figures. Hidden by default: it is not a population but the leftover - every CD45+ cell matching no definition in the spec - and at ~36% of cells here it blankets the real populations and dominates every legend. Display only: those cells stay gated, counted and present in cells_umap.csv, gate_counts.csv and every frequency table either way
--flowjo-export	off	ALSO write UMAP-annotated FCS files for interactive use in the FlowJo GUI (§8.2). Strictly additive - every default output is still produced
--flowjo-outdir	<outdir>/flowjo	where --flowjo-export writes
--flowjo-no-concat	off	skip the pooled _ALL_SAMPLES.fcs
--flowjo-no-groups	off	skip the per-study _GROUP_<cohort>.fcs files

--max-events-per-file samples evenly across the acquisition rather than taking the first N events, because acquisition order is time order: a leading slice biases toward whatever was flowing in the first minute.

7. What the pipeline does, step by step

The console log is organised into the same numbered steps.

1. Reading files and resolving markers. Marker names come from the FCS parameter descriptions, falling back to detector names. Files are grouped into panels by marker fingerprint, so a cohort acquired with a different antibody set is embedded separately instead of being silently merged. Compensation is applied from each file's spillover matrix where usable.
2. Deriving transform and gates. An asinh cofactor is derived per panel, then per sample: singlets from an FSC-A/FSC-H MAD band, live cells by the viability dye, and CD45+ by a valley between modes. Where a marker is unimodal there is no valley, and the script says so rather than inventing a cut.
3. Initial QC diagnostics. recon_diagnostics.png - the gates as drawn, per sample. Written before any population is scored, so a bad gate is visible before you have numbers to be attached to.
4. Scoring populations. The declared gate hierarchy is applied. Populations whose markers are absent from a panel are reported in populations_unavailable.csv instead of being reported as zero - absent and empty are different claims.
5. Patient metadata. Translation, column filtering, age derivation, join.
6. Embedding. Balanced subsample of CD45+ events per sample → one UMAP per panel → coloured overlays.
7. Writing tables. Frequencies, MFI, functional-marker blocks, per-cell table.
8. Between-group comparison. Groups resolve from the patient table via the sample map - declared in study metadata rather than parsed out of filenames. Each population is tested against the reference group and the lettered panel figure is drawn. Skipped, with a logged reason, if no group column resolves.
9. Saving session state (skip with --no-session).

8. Output and how to read it

Multiple marker panels. Every file below is per-panel: with one panel (the common case) filenames are exactly as listed. With more than one - different marker sets in the same run, e.g. two acquisition batches on different panels - every filename below (figures, tables, and the group comparison / functional-marker / ratio stats CSVs) gets a _<panel-name> suffix, one full set per panel, because a population scored under one marker set and threshold is not comparable to the same population name scored under another. umap_overview.png/umap_markers.png already worked this way; §7 step 1 has the full picture of how panels are detected.

Figures

File	What it shows
recon_diagnostics.png	Every gate, per sample, with its threshold drawn on the distribution
gating_qc.png	Gate hierarchy per sample with event counts
umap_overview.png	The shared embedding coloured by population, sample, and each discovered covariate
umap_markers.png	The same embedding, one panel per marker, intensity-coloured
umap_density.png	The shared embedding as 2D density, one small panel per sample
umap_density_by_group.png	Same shared axes, one small panel per group instead of per sample (e.g. per cohort) - the side-by-side view companion to umap_overview.png's single overlaid plot. Written only when --group-column resolves to 2+ groups
umap_overview_by_group.png	Combined-vs-per-group small multiples: for every colouring umap_overview.png draws (population, sample, each covariate), one row of panels - "All samples" pooled, then one panel per group, all on the SAME fixed axes and (population/covariate rows) the SAME colour-to-category mapping, so a group-specific shift is visible directly rather than inferred from overlaid colours. Written only when --group-column resolves to 2+ groups
population_frequencies.png	One lettered panel per population: mean bar, median rule, SD whiskers, every sample as a point (all samples pooled, % of CD45+)
group_comparison.png	The same layout split by group, with Wilcoxon brackets against the reference group. Written only when a group column resolves
functional_markers.png	Same layout, one panel per population × functional marker (e.g. HLA-DR on classical monocytes): % positive, split by group with brackets once a group column resolves
population_ratios.png	Same layout, one panel per ratio declared in the config's ratios: block (e.g. CD4:CD8), split by group with brackets once a group column resolves
absolute_counts_qc.png	One tile per sample × population from --absolute-counts, log-scale colour. Inspect before absolute_counts.png - see §5.5. Written only when --absolute-counts is supplied
absolute_counts.png	Same bar/point/brackets layout as group_comparison.png, for the --absolute-counts values instead of this pipeline's own gated frequencies

| umap_multigraph_overlay.png | One composite image answering "what differs, and where". Top: the shared embedding pooled, then one panel per study, with every cluster-subcluster labelled (4a, 4b, …) at the same position in every panel. Below: one peak panel per cluster-subcluster-marker (e.g. 4-b-CD3), each overlaying one curve per study. A port of FlowJo's Layout Editor Make Multigraph Overlay → Histograms. Written only when --group-column resolves to 2+ groups |

Tables

File	Grain
cells_umap.csv	one row per embedded cell: markers, population_label, umap_1/2, sample and patient covariates
subcluster_marker_shifts.csv	one row per cluster × subcluster × marker × study, sorted by |Cliff's delta|. Makes the overlay figure sortable: occupancy_pct_delta (does this study put a different share of its cells in this compartment) and median_delta / cliffs_delta (do the cells that are there express the marker differently)
population_frequencies.csv	sample × population, count, pct_of_parent, pct_of_cd45_pos, and cells_per_ul when a blood count is supplied
group_comparison_stats.csv	one row per population × group: n, mean, median, SD, fold change, Cliff's delta, raw p, BH-adjusted p
population_marker_mfi.csv	every marker in every population - unscoped
functional_markers.csv	functional markers scoped to the populations where they are interpretable
functional_markers_stats.csv	same test as group_comparison_stats.csv, one row per (population × marker) × group. Written only when a group column resolves
population_ratios.csv	sample × ratio: numerator/denominator population, their pct_of_cd45_pos, and the ratio (from the config's ratios: block, e.g. CD4:CD8)
population_ratios_stats.csv	same test as group_comparison_stats.csv, one row per ratio × group. Written only when a group column resolves
absolute_counts.csv	sample × population, matched cells_per_ul from --absolute-counts, plus is_control/qc_status
absolute_counts_stats.csv	same test as group_comparison_stats.csv, one row per --absolute-counts population × group. Written only when a group column resolves
absolute_counts_raw.csv	the --absolute-counts input, normalized to CSV verbatim - open this first if a match or a number looks wrong
thresholds_used.csv	every threshold, its source (valley / quantile_fallback) and needs_review
staining_qc.csv	per-sample verdict and whether it entered the embedding
gate_counts.csv	events surviving each gate step
populations_unavailable.csv	populations not scorable in a panel, and why
patient_metadata_english.csv	translated, column-filtered patient table

Reading the output in the right order

Start with thresholds_used.csv and recon_diagnostics.png, before any frequency table. Every downstream number inherits the gates; a frequency table computed from a misplaced threshold looks exactly as authoritative as a correct one.

Specifically:

• source = quantile_fallback means no valley was found - the marker was unimodal, and the "threshold" is a quantile cut carrying no evidence of a positive population. Any population depending on it is unreliable.
• needs_review = TRUE marks a threshold with no evidence behind it. It is set when and only when the quantile fallback was used, so it is a per-sample flag about that threshold's derivation - not a comparison against other samples.
• scale_outlier = TRUE is the cross-sample check: for each marker within a panel, thresholds are compared against the median of the other samples, and a threshold sitting far from its peers is flagged. This is the column that catches a tube whose gates are individually plausible but inconsistent with the rest of the cohort - the failure mode that silently makes per-sample frequencies non-comparable. threshold_scale_qc.csv reports the same comparison per marker.
• Read scale_outlier alongside the spread of a marker's thresholds. Within one panel they should be broadly similar. A wide spread - especially thresholds near or below zero alongside thresholds in the normal range - means gates are being placed inconsistently. A threshold at or below zero makes essentially every event "positive".
• Then sanity-check against immunology. In blood, CD4 T cells outnumber CD8 T cells - population_ratios.png turns this into an explicit per-sample number instead of an eyeballed heuristic, and is also where a clinically-inverted ratio would first become visible; monocyte subsets are ordered classical > intermediate ≈ non-classical. An all-or-nothing pattern across samples - one population near 0% wherever another is high - is a threshold artefact, not biology.
• Large Other CD45+ means CD45+ cells matched no declared population. Some is expected. A large compact island of it is a population your gates are missing.
• Check island composition before believing any group difference. If a cluster is dominated by one sample, you are looking at a batch effect. Draws are balanced per sample, so a genuine cluster should contain many samples. umap_density_by_group.png and umap_overview_by_group.png are the direct check: the same shared axes, split one panel per group instead of overlaid, so a cohort-specific shift or dropout is visible without hunting through colours in umap_overview.png. The latter keeps population colouring (and the sample panel) per group, not just density, and puts "All samples" in column 1 of every row for a direct pooled-vs-per-group comparison.

Excluded samples still appear in staining_qc.csv and the QC figures - with the reason - so exclusions are auditable rather than invisible.

--include-qc-failed forces excluded samples back in, for cohorts too small to lose any. A sample can fail for two reasons: no separable CD45+ mode was found (CD45 was unimodal - the fallback threshold that produces a number is a fixed top-quantile cut with no evidence of a real positive population behind it), or its stained fraction fell below --min-cd45-pct. Either way, the flag keeps the sample in gating/UMAP/frequencies using whatever threshold was already computed - it does not derive a better one. staining_qc.csv's verdict column says "included anyway (--include-qc-failed)" for every sample this affects, so you can still tell which numbers came from a real gate versus a forced one before trusting them in a comparison.

8.1 Between-group abundance comparison

group_comparison.png and group_comparison_stats.csv answer "do the groups differ in the size of each population". They are produced whenever a group column resolves from the patient table; with a single group only the pooled population_frequencies.png is written.

functional_markers.png/.csv and population_ratios.png/.csv answer the same question for a different quantity - marker positivity within a population (e.g. does HLA-DR on classical monocytes differ), and a derived ratio (e.g. CD4:CD8) - using the identical bar/median/whisker/points layout, Wilcoxon brackets, and BH correction described below. Everything in this section applies to all three equally; only the y-axis and its caveat differ.

# --group-column defaults to the resolved 'cohort' column
Rscript scripts/fcsumap.R --dir data --recursive \
  --patient-table patients.csv --sample-map samples.csv \
  --reference-group "controls" --outdir results

Which quantity is on the y-axis matters more than the p-value. Three columns in the frequency table look interchangeable and are not:

• count is never comparable between samples. It is an event count - it depends on how long the operator ran the tube and on --max-events-per-file. It is never used as a y-axis here.
• pct_of_cd45_pos is comparable but compositional. Populations are constrained to sum to 100%, so one lineage expanding mathematically forces the others down. A significant fall in "% NK" is equally consistent with NK loss and with granulocyte gain. The figure prints this caveat in its subtitle whenever frequencies are what is being shown.
• cells_per_ul is absolute and independent per population. This is what the y-axis switches to automatically once a wbc_per_ul column (leukocytes/µL from a haemogram) exists in the patient table. Without it, no absolute concentration can be produced from FCS files alone - the instrument records events, not volume.

The tests. Each non-reference group is compared against the reference with a Wilcoxon rank-sum test, plus a Kruskal–Wallis omnibus test when three or more groups are present. Rank-based tests are used rather than t-tests because group sizes here are single- to low-double-digit and abundance distributions are right-skewed with outliers, so normality is not a safe assumption.

Read the adjusted p-value column, not just the brackets. Testing 12 populations against 2 groups is 24 tests; under the null, a raw p < 0.05 somewhere is more likely than not. Both p_value and p_adj_BH (Benjamini–Hochberg FDR) are always written to the CSV, and the run log states how many hits survive correction. Brackets show raw p by default - --p-adjust-display BH switches them. If nothing survives correction, the log says so explicitly: treat those hits as hypotheses, not findings.

Read cliffs_delta alongside p. It is a bounded effect size (−1 to +1) that stays interpretable at n = 6, where a ratio of means can be dominated by a single outlier. A large delta with a marginal p is an underpowered real effect; a small delta with a small p is usually a technical artefact.

Every sample is drawn as a point, deliberately. At these group sizes a bar and a whisker cannot show whether a difference rests on one outlier or on a consistent shift. The bar height is the mean and the black rule is the median - when they diverge, the distribution is skewed and the median is the more honest summary. Points are jittered horizontally only, so no point's vertical position (its actual value) is distorted.

All the gating caveats above still apply. A group difference computed from a quantile_fallback threshold, or from a marker flagged scale_outlier, is a threshold artefact with a p-value attached. Check thresholds_used.csv before believing any bracket.

8.2 Subclusters, the multigraph overlay, and the FlowJo hand-off

Subclusters are fitted on the reference group only. With --reference-group set, each cluster is split by k-means run on the reference cells alone; every other study's cells are then assigned to the nearest reference centroid. Clustering the pooled cells instead would let patient cells help draw the very boundaries they are then tested against - a disease-specific subset would carve out its own subcluster and come back looking like a normal compartment that merely contains patient cells. Fitting on healthy cells makes each subcluster a piece of normal biology and asks the honest question: given the compartments a healthy immune system has, where do patient cells fall, and do they look the same once they get there?

That gives the two ways a study can deviate, both reported:

Mode	Column in subcluster_marker_shifts.csv
Occupancy - the compartment holds a different share of the study's cells	occupancy_pct_delta
Intensity - the cells that are there express a marker differently	median_delta, cliffs_delta

Subcluster letters are ordered by position on the embedding (left to right), not by a marker coordinate, so 4a → 4b → 4c trace the cluster across the map and the peak rows follow the same spatial sequence.

Peaks are mode-normalised (each curve scaled so its own maximum is 100 - FlowJo's "Count (%)"). Study sizes differ several-fold, so on a shared density axis the largest study's curve simply towers over the rest and the figure quietly becomes a cell-count comparison. Read peak position and shape here; abundance belongs to group_comparison.png, where it has statistics. A curve is omitted where a study has fewer than 20 cells in the compartment - a KDE from fewer is noise shaped like a result.

FlowJo hand-off. --flowjo-export writes UMAP-annotated FCS files alongside every normal output. FlowJo's command line cannot compute a UMAP, so the embedding is carried in as ordinary FCS parameters (UMAP-1, UMAP-2, plus integer-coded Population, SampleID, CohortID) - see docs/flowjo_interactive.md for why extra parameters rather than a CLR file, and for the GUI steps.

File	Contents
_ALL_SAMPLES.fcs	every cell pooled - open this first
_GROUP_<cohort>.fcs	one per study, for drag-and-drop cohort overlays
<sample>.fcs	one per sample
population_codes.csv	integer → label legend for the three code channels
samples.txt	sample list the FlowJo CLI accepts as a first argument

Cluster numbers in umap_multigraph_overlay.png are the same codes as population_codes.csv and the Population channel of the exported FCS - one numbering across figures, tables and the GUI.

9. Reproducibility

• --seed fixes the embedding; the pinned Docker image fixes the library versions that the seed is only meaningful relative to.
• --reference-date fixes derived ages.
• --write-config captures every derived threshold, so a run can be repeated exactly, and any manual override is a visible diff rather than a remembered decision.
• Version-control the config, the sample map schema and the script. Do not version-control raw FCS or patient-linked metadata: .gitignore excludes both, and de-identified templates in scripts/examples/ keep a fresh clone runnable.

10. Licence

MIT - see LICENSE.