Bone-marrow cell populations

A patient's bone marrow is measured one cell at a time. The quantity to report is the fraction of cells in each of six populations. This page reduces every cell to one of 8 categories chosen from the estimation problem, freezes the categories, and asks how well the fractions of patients the categories never saw can be recovered from 8 integer counts.

1. The cells and a patient's composition

Flow cytometry passes cells one by one through a laser and records, for each cell, a forward-scatter and a side-scatter intensity and one fluorescence intensity per marker channel: twelve numbers per cell here. In the FlowCyt benchmark an expert has labelled every cell of thirty bone-marrow aspirates as one of six populations: T cells, B cells, monocytes, mast cells, the CD34-positive progenitors (HSPCs), and everything else. The clinical quantity is the composition, the fraction of a patient's cells in each population, and the differences between patients are what a report is about.

• T cells
• B cells
• monocytes
• mast cells
• HSPCs

▸ marks the ten held-out patients, frozen before any fitting.

Each bar is one patient's composition over the five target populations, drawn at their true fractions; the dominant remaining compartment, about nine tenths of every patient, is left out so the variation is visible. The spread between patients is the thing being measured, and it is much wider than any of the estimation errors below.

The cohort was split once, before anything was fitted, into twenty reference patients the categories are built from and 10 held-out patients that nothing downstream has seen. Every number on this page that is called held out was computed on those 10 patients.

2. Why compress the cells

A patient's cells can be fitted directly, cell by cell, and the result section below does exactly that as the reference. The reason to reduce them to a handful of integer counts is that counts are what can be archived, transmitted, audited and compared across instruments: a category rule plus 8 integers is a summary a laboratory can publish, whereas a hundred thousand continuous classifier outputs is not. The cost of the compression is the question.

There is a floor on how few categories can work. Six fractions that sum to one have 5 free coordinates, and a fixed number of cells spread over $B$ categories gives only $B - 1$ independent counts. Five categories therefore cannot identify the composition, whatever they are: with 4 independent counts for 5 unknowns the fit has a flat direction. The bin-budget section shows that this is exactly what happens.

3. From marker channels to a score

FS INT

SS INT

FL1 INT_CD14-FITC

FL2 INT_CD19-PE

FL3 INT_CD13-ECD

FL4 INT_CD33-PC5.5

FL5 INT_CD34-PC7

FL6 INT_CD117-APC

FL7 INT_CD7-APC700

FL8 INT_CD16-APC750

FL9 INT_HLA-PB

FL10 INT_CD45-KO

• T cells
• B cells
• monocytes
• mast cells
• HSPCs

Every channel, drawn as a density per population, with the dominant compartment hidden so the rest are visible at all. The horizontal axis is the robust arcsinh (cofactor 150), the study's own transform: on a linear intensity axis every panel collapses into one spike against the origin. The two scatter channels separate the broad compartments, CD34 is what makes the smallest compartment visible, and CD45, the pan-leukocyte backbone, is where everything overlaps. No single channel separates the populations. What separates them is the joint distribution of all twelve, and that is what the score captures.

ScoreQuant works with the score of a cell: the derivative of the log likelihood of the cell with respect to the parameters being estimated, one coordinate per parameter. For a mixture of six populations with fractions $\theta$ and the last population absorbed as the simplex reference, the score has five coordinates, and each is a function of the ratios between the population densities at the cell. Nobody has those densities, so they are estimated: a classifier cross-fitted across the reference patients gives each cell a calibrated posterior over the populations, dividing by the training proportions turns posteriors into density ratios, and the mixture parameterization turns ratios into the score.

• T cells
• B cells
• monocytes
• mast cells
• HSPCs

Two of the boxes in that chain are deliberately not ScoreQuant: the classifier is not, and the mixture likelihood that finally turns counts into fractions is not. What the library receives is a table of score rows with their weights, declared as an estimated ratio rather than an exact score, so nothing downstream can describe the result in exact-Fisher language. The classifier, its calibration and its cross-fitting protocol are on the score-model page (opens in a new tab) of the reference study.

4. The reference patients and the call

The code on this page runs the whole chain on a committed fixture-scale table: 7,038 reference cells the rule is fitted on, 3,518 further reference cells kept apart to calibrate the categories, and 10,240 cells of the 10 held-out patients. The numbers displayed on the page come from the full study on 600,000 cells, and the closing section says what the fixture-scale run gives instead.

import numpy as np

import scorequant as sq
from examples.cell_population import fit_binned_mixture, fit_unbinned_mixture
from examples.cell_population.likelihood import estimate_bin_templates

table = np.load("examples/data/flowcyt_walkthrough.npz")
parameters = tuple(str(name) for name in table["schema"])
theta0 = table["theta0"]

reference = sq.ScoreSample(
    table["partition_scores"].astype(float),
    table["partition_weights"],
    schema=sq.ScoreSchema(parameters),
    provenance=sq.ScoreProvenance(
        kind="estimated_ratio",
        description="patient-cross-fitted calibrated FlowCyt mixture score",
    ),
)
assert reference.provenance.exact_fisher is False

The weights make the reference measure: every reference patient counts equally within a population, and the populations are weighted by the reference composition, so no single patient's abundance sets the categories. The call takes the source, a bin budget, a criterion and the configuration of the solver that criterion selects.

quantizer = sq.fit_quantizer(
    reference,
    n_bins=8,
    criterion=sq.DOptimality(),
    config=sq.SoftVoronoiConfig(seed=2026, initializer_restarts=4, max_steps=160),
)
assert quantizer.n_bins == 8
assert quantizer.information_kind == "supplied_score_surrogate"

DOptimality treats every population fraction as equally worth keeping: it maximises the determinant of the information the 8 labels carry about the five fractions, which is the same as minimising the volume of their joint uncertainty ellipsoid. The soft Voronoi solver relaxes the hard assignment into a temperature-annealed soft one, optimises the cell centres, and hardens at the end. The result is a rule: a function from a cell's five score coordinates to a category, which is what makes it applicable to patients it never saw. This is the distinction the whole page rests on. A partition of the reference rows would label those rows and stop; fit_quantizer returns something that labels the next patient.

5. What the eight categories are

A five-dimensional partition cannot be drawn, but it can be described without projecting it: for each category, the probability that a cell landing there belongs to each population, at the reference composition.

Reference population composition of each learned bin, as probabilities that sum to one across a row.

Bin	T cells	B cells	monocytes	mast cells	HSPCs	other
bin 6	0.94	0.00	0.00	0.00	0.00	0.06
bin 1	0.01	0.86	0.01	0.00	0.01	0.11
bin 4	0.16	0.33	0.17	0.03	0.19	0.12
bin 2	0.15	0.31	0.15	0.08	0.18	0.13
bin 5	0.01	0.03	0.84	0.00	0.01	0.10
bin 3	0.02	0.04	0.02	0.00	0.78	0.14
bin 0	0.00	0.00	0.00	0.00	0.00	1.00
bin 7	0.06	0.12	0.06	0.01	0.31	0.45

Each row is one of the 8 learned categories; each entry is the probability that a cell in that category belongs to the column population, at the reference composition.

Six of the 8 categories act as gates: each is dominated by one population, including one for the rare progenitors and one for the mast cells, which is the category a single marker could not have found. The remaining two are mixed. They are not a failure of the fit. A gate per population would spend six categories and identify the composition only through the gates' purities; the mixed categories are where the classifier is uncertain, and keeping them as their own counts is what lets the fit separate the populations along the directions the gates confuse. The criterion chose them because they carry information about the fractions, not because they look like anything on a two-channel plot.

6. Applying the frozen rule to held-out patients

The rule is frozen. It is applied to the calibration cells, whose expert labels give each category's probability under each population, $P(\text{bin} \mid \text{population})$; and to the held-out cells, whose predicted categories are counted per patient. The count vector and the calibration table are the whole input to the mixture fit. The held-out labels are used only afterwards, to score the estimate.

template_bins = quantizer.predict_scores(table["template_scores"].astype(float))
templates = estimate_bin_templates(
    table["template_labels"], template_bins, table["template_patients"], n_bins=8
)

held_out_scores = table["held_out_scores"].astype(float)
held_out_patients = table["held_out_patients"]
held_out_bins = quantizer.predict_scores(held_out_scores)

patients = np.unique(held_out_patients)
counts = np.array([np.bincount(held_out_bins[held_out_patients == p], minlength=8) for p in patients])
binned = np.array([fit_binned_mixture(row, templates).fractions for row in counts])
unbinned = np.array(
    [
        fit_unbinned_mixture(table["held_out_ratios"][held_out_patients == p].astype(float)).fractions
        for p in patients
    ]
)
expert = np.array(
    [
        np.bincount(table["held_out_labels"][held_out_patients == p], minlength=6)
        / np.count_nonzero(held_out_patients == p)
        for p in patients
    ]
)
assert counts.shape == (len(patients), 8)
assert binned.shape == unbinned.shape == expert.shape

The unbinned fit is the reference the compression is measured against: the same classifier's continuous ratios, straight into the mixture likelihood, with no compression at all. The macro RMSE below is the root-mean-square error of the fitted fractions over the 10 held-out patients, averaged over the five target populations.

def macro_rmse(estimate):
    per_population = np.sqrt(np.mean((estimate - expert) ** 2, axis=0))
    return float(np.mean(per_population[:5]))


binned_rmse = macro_rmse(binned)
unbinned_rmse = macro_rmse(unbinned)

7. The result

Held-out macro RMSE, from the counts

0.001929

Error in the recovered population fractions from the eight counts, averaged over the five target populations and the held-out patients. The quantity this page estimates.

Held-out macro RMSE, unbinned

0.001728

The same estimate with no compression: the continuous classifier ratios straight into the mixture fit. The floor the counts are measured against.

Held-out D-efficiency

0.9846

Geometric mean of the eigenvalues of the binned information relative to the unbinned information on the held-out cells: the fraction of the score's information the eight labels keep.

Information kind

supplied_score_surrogate

What the library calls this result. The scores were estimated by a classifier, so it is a surrogate and never exact Fisher information.

The full study on 600,000 cells, at 8 categories. Every number is measured on the 10 held-out patients.

• T cells
• B cells
• monocytes
• mast cells
• HSPCs

Filled: fitted from the eight counts. Hollow: fitted from the unbinned classifier ratios. An estimate of exactly zero is drawn on the lower edge.

Every held-out patient and every target population sits on the diagonal at every abundance, from the mast cells near one part in ten thousand to the T cells near one in ten, and the filled marks from the 8 counts sit where the hollow unbinned marks sit. Compression costs, as it must: the information carried by a label is never more than the information carried by the thing it replaced, and the binned error is above the unbinned one. The finding is the size of the cost. The two errors are of the same order, and both are far below the spread between patients in the first chart, so replacing a patient's cells with 8 integers does not change what the report says.

The exact finite-D solver, run on the same reference rows and compiled into its Voronoi rule, reaches 0.002094 on the same patients with a D-efficiency of 0.9853: slightly worse on the error, slightly better on the information. The two disagree in opposite directions because the downstream fit is a fixed-total multinomial likelihood on the counts, not the local Fisher quadratic the criterion optimises, so any sentence about this study has to say which method produced the number in it. Both are on the solver-comparison page (opens in a new tab).

8. Why it works: the retained information

The D-efficiency is the explanation, not the headline. The 8 labels keep 0.9846 of the information the estimated score carries about the five fractions, on cells the rule never saw; the criterion optimised that quantity on the reference cells, and it transferred. A high retained information is a local, quadratic statement about the likelihood near the reference composition, and the held-out patients differ from the reference composition, so the retention does not guarantee the estimation error; it makes it unsurprising.

What the retention measures is how much of the estimated score the labels keep, a number that stays high whether or not the estimate itself is good. The density-ratio page measures the gap between those two things on a problem where the exact score is known; the unbinned row of the table above is this page's own check, because a poor estimate would have shown up there first.

9. What simpler reductions retain

Every method below sees the same reference cells, gets the same budget, is frozen the same way and is scored the same way: fit the fractions from the counts on the held-out patients, and measure the result against the expert composition.

Held-out D-efficiency at 8 bins, compared across 5 binning methods

Method	Held-out D-efficiency at 8 bins
Grid on two marker PCs	0.0223
One score direction	0.0378
k-means on the markers	0.0704
Soft Voronoi (ScoreQuant)	0.9846
D exchange (ScoreQuant)	0.9853

Three ordinary reductions against the two information-optimal ones, at the same bin budget on the same held-out cells.

A rectangular grid on the two leading principal components of the marker channels is a two-dimensional gate with straight edges. Clustering the twelve-dimensional marker vectors groups cells that look alike and counts the groups. Cutting a single score direction into intervals keeps the classifier output but summarises it in one dimension. All three land an order of magnitude away on estimation error as well, at every budget:

• D exchange (compiled)
• k-means on raw markers
• One score direction
• Normalized trace (k-means)
• Soft Voronoi
• Two-dimensional grid
• Unbinned classifier ratio

10. The bin budget

• D exchange (compiled)
• k-means on raw markers
• One score direction
• Normalized trace (k-means)
• Soft Voronoi
• Two-dimensional grid

Five bins fail, as the second section said they must: the held-out D-efficiency is 0.3911 and the macro RMSE 0.030002, because 4 independent counts cannot resolve 5 fractions and the fit has a flat direction. Eight is the knee: 0.9846 retained and 0.001929 error. Ten bins retain 0.9903 and reach 0.002023; thirty retain 0.9984 and reach 0.002600. More categories keep more of the local information, and the fraction estimates get no better, because the calibration table has more entries to estimate from the same reference cells and the smallest categories run empty on some patients. The simpler reductions do improve with more bins, the grid from 0.083325 at eight to 0.062011 at thirty, and remain an order of magnitude away.

11. Limits, and how to reproduce this

The fixture-scale run is not the full study. The code on this page, run on the committed table, gives a held-out macro RMSE of 0.0090 from the counts against 0.0016 unbinned, with a held-out D-efficiency of 0.9774. The information retention is close to the full study's; the estimation error is several times worse, because the calibration table $P(\text{bin} \mid \text{population})$ is estimated from 3,518 cells instead of hundreds of thousands, and at that size the calibration, not the categories, limits the fit. The full study's numbers, the ones displayed above, come from 600,000 cells drawn from the complete benchmark.

The score is a surrogate. The classifier's ratios are estimated, so the result is recorded as supplied_score_surrogate, never as exact Fisher information.

Aiming at one fraction does not help here. The study also asked for the categories that best resolve the progenitor fraction alone, with the rest of the composition profiled out. On the committed fixture the plain D partition already retains 0.9954 of that profiled information against a certified ceiling of 0.9997, so there was almost nothing left for a profiled criterion to win. That negative result is on the profiled-criterion page (opens in a new tab).

The data. The benchmark is © Bini, Nassajian Mojarrad, Liarou, Matthes, Marchand-Maillet (CHIL 2024), licensed CC-BY-NC-SA-4.0, and lives in its source repository; this repository does not redistribute it. What is committed is a small deterministic fixture, the fixture-scale walkthrough table this page runs on, and the evidence file of the full study. The acquisition protocol, the per-patient estimates, the solver comparison and the uncertainty audit are in the FlowCyt case study (opens in a new tab).

JAX_ENABLE_X64=1 MPLBACKEND=Agg uv run python -m examples.cell_population --full
uv run python website/scripts/generate_showcase.py --force