diff --git a/.changeset/fitted-background-equal-width-bins.md b/.changeset/fitted-background-equal-width-bins.md new file mode 100644 index 0000000..32a90e0 --- /dev/null +++ b/.changeset/fitted-background-equal-width-bins.md @@ -0,0 +1,43 @@ +--- +'@platforma-open/milaboratories.feature-integration.per-cell-metrics': minor +'@platforma-open/milaboratories.feature-integration.ui': minor +'@platforma-open/milaboratories.feature-integration.model': patch +--- + +The fitted background draws the distribution the run holds + +The fitted-background grid is the only way to see whether a tag's counts separated into two +populations, so a hump it invents, or one it hides, is the error this surface cannot carry. Three +things were wrong with it. + +**The zeros were missing.** The plot was binned from the sparse counts frame, which has no rows for +cells that read nothing. Those zeros are most of the background, so the plot showed one decaying hump +whatever the fit had found — the left half of the distribution was simply absent. It is now binned over +the cells the fit actually ran on, one entry per cell in the sample, and a cell that read nothing counts +as a zero. + +**Every bar is now the same width.** Bins sit at `expm1(k * 0.2)`, uniform in `log1p`, which is what +the plot's axis already is. The source paper histograms `log1p` counts at a fixed width for the same +reason. + +Before, bins were whole numbers stepping geometrically, so their widths ran from 0.301 of a decade at +`[0, 1)` down to 0.079 at `[4, 5)`. A raw count then made a wide bar stand above a narrow one holding +the same density, so each bar had to be divided by its own width — and dividing by the wrong width hid +a real signal component completely: on a mixture whose upper mode was cleanly separated at a mean of +60, that hump drew at 0.9% of the background peak. Equal widths remove the division and the error with +it. Bar height is a plain cell count again, so the hover readout reports a cell count. + +The cost of equal widths is that the edges are no longer whole numbers, and counts are — consecutive +integers sit further apart than one bin until about count 13, so the low end is a comb of separated +bars. The paper's own figures show the same gaps. `LOG1P_BIN_WIDTH` is coarser than the paper's 0.075 +for that reason: at 0.075 a real tag came back with 75 of 97 bins empty. + +**The bound line is drawn.** Under a fitted baseline the threshold is a probability, so a plot in counts +had nothing to mark it against. Each fit now resolves the count at which the run's bound probability +starts calling a cell bound, and the panel draws it. A fit that reaches no such count draws no line and +says so, rather than marking one at the bottom. + +Each (sample, tag) is binned against its own range, so the emitted weight lists have different lengths +and are shorter than the shared edge set. The plot pads them, which draws the same picture — past a +pair's own maximum every bar is empty either way. Binning inside the fit rather than returning per-cell +arrays takes the fitting step from 2.1 million numbers held to 7,764 on a 27-sample run. diff --git a/.changeset/fitted-baseline-settings.md b/.changeset/fitted-baseline-settings.md new file mode 100644 index 0000000..5059dc7 --- /dev/null +++ b/.changeset/fitted-baseline-settings.md @@ -0,0 +1,44 @@ +--- +'@platforma-open/milaboratories.feature-integration.per-cell-metrics': minor +'@platforma-open/milaboratories.feature-integration.workflow': minor +'@platforma-open/milaboratories.feature-integration.model': minor +'@platforma-open/milaboratories.feature-integration.ui': minor +'@platforma-open/milaboratories.feature-integration.test': patch +'@platforma-open/milaboratories.feature-integration': minor +--- + +The fitted baseline states where it starts, and the scientist can move it + +Two settings appear under the fitted baseline, and nowhere else — neither reaches the declared or +panel rung, so neither is offered there. + +**Expected binder %.** Roughly what share of cells are expected to bind an antigen. The fit splits the +counts at the matching quantile and seeds one component from each side. It is not a threshold: the EM +re-estimates both components from there, so the split the run ends up with is an output of the fit. + +It changes answers anyway, because the EM is not globally convergent on these distributions and the +start decides which optimum it reaches. On a panel where 27% of cells really did bind, the shipped +value put the split at 953 counts; told 30%, the same fit put it at 13 — which is where the gap in that +tag's histogram actually is. The published value comes from a rare-binder regime, and the study behind +this rung never tested a positive fraction above 25%. + +The trade runs one way, so no single value is right: raising it also makes the fit readier to carve a +signal component out of a single population, so a tag that bound nothing invents more binders. Only +the scientist knows which side of that to be on, which is why it is a setting. + +**Bound probability.** How sure the fit must be before a cell counts as bound. Previously fixed at +0.9 with no way to see or move it. Now shown, with 0.9 as both the default and the lowest accepted +value — below it a cell holding none of a tag could cross the line, and the run counts those cells by +arithmetic rather than reading each one, so the two halves would disagree with nothing raised. + +**The fit now starts where the method says.** The split was taken at the median, which +`what-plays-the-baseline` never specified. A median start begins from two halves of equal size, which +is far from the truth on a mostly-background population — every tag here — and pulls the fit toward +calling much of that background signal. On a control reagent, whose counts hold one population, a +median start gives a background weight near 0.8 against 0.95 from the published split. + +That trade is not free, and the direction is recorded in the suite: on a background whose long tail +puts its mean above the binders', the published start decomposes the counts into the bulk and the +tail rather than into background and binders, and calls the tail the signal. The median start got that +shape right and the mostly-background case wrong instead. Neither wins both. The run gives no warning +in either case, which is why the fitted grid puts both means in front of the reader. diff --git a/.changeset/per-cell-tag-counts-export.md b/.changeset/per-cell-tag-counts-export.md new file mode 100644 index 0000000..353db28 --- /dev/null +++ b/.changeset/per-cell-tag-counts-export.md @@ -0,0 +1,24 @@ +--- +'@platforma-open/milaboratories.feature-integration.per-cell-metrics': minor +'@platforma-open/milaboratories.feature-integration.workflow': minor +'@platforma-open/milaboratories.feature-integration': minor +--- + +Per-cell antigen counts are exported, before the minimum is applied + +A new export, keyed `[sampleId, cellId, tagId]`, carrying the UMI count each cell held for each +barcode. It is what a downstream per-cell composition plot needs: one bar per cell, split by antigen. + +**Before the count minimum, and that is the point.** The counts the verdicts are computed on have +already had every value below the minimum set to zero, and the comparator tag is exempt from that — +so on a declared-baseline run the control keeps its small counts while the antigens lose theirs. Those +numbers answer "what counted as evidence of binding". This export answers "what did the cell capture", +so a cell's tags add up to what that cell actually held. The column's own description says so, because +the two do not reconcile and a reader who mixes them draws the wrong conclusion. + +It cannot be derived from the floored counts afterwards: once the minimum has run, a count of 3 and a +count that was never there are both 0. + +Partitioned by sample, the only column in the block that is. It is the largest table the run produces, +at one row per (cell, tag), and a composition plot reads one sample per view — so a reader after one +sample touches one partition instead of the whole run. diff --git a/.changeset/run-quality-grid-per-sample.md b/.changeset/run-quality-grid-per-sample.md new file mode 100644 index 0000000..7c53b30 --- /dev/null +++ b/.changeset/run-quality-grid-per-sample.md @@ -0,0 +1,38 @@ +--- +'@platforma-open/milaboratories.feature-integration.ui': minor +'@platforma-open/milaboratories.feature-integration.workflow': patch +'@platforma-open/milaboratories.feature-integration.model': patch +--- + +The fitted-background grid reads one sample at a time + +The fit runs per (sample, tag), so the grid drew one panel per pair — 27 samples and 9 barcodes is 243 +panels on one page, and every title had to repeat the sample to tell them apart. It now shows one +sample, chosen from a selector above it, so a panel is titled by its reagent alone and the sample is +named once. Barcodes read in the order the panel file declares them, so a barcode holds the same slot +whichever sample is shown. + +What that gives up is reading down one reagent across samples. The grid still supports that shape; +nothing asks for it today. + +Each panel's caption now carries the bound count alone, and the fit's own numbers moved to the enlarged +panel, where a value is read rather than scanned. The cell count that used to lead the caption was the +sample's analysed population — the same number on every panel of the sample — and now sits once above +the grid. + +**The Panel column is out of the tables.** It is a hash of the sorted barcode list, because no panel +file names its panel, and a run declaring one panel for every sample repeated that hash identically on +every row. It stays available in the column picker, and a multi-panel run should switch it on: +`Seen in 2/3` cannot be read without knowing which three samples. + +**Fixes** + +- A fitted background mean of 0.000488 printed as `0`, a value the fit cannot produce — three + significant figures were computed and then discarded by a formatter keeping three decimal places. +- The Run quality page failed to render at all: a watch read a value declared further down the file. +- A run computed before the bound count existed reported "no count reaches the bound probability", + stating a finding no run had produced. Absent and null now read differently. +- Resizing the window redrew every panel on every frame, and each redraw leaks a tooltip node in the + uikit. A few pixels of tolerance takes a drag from hundreds of redraws to a handful. The leak itself + is the uikit's. +- The quality-report JSON was pretty-printed, which roughly doubled it for a file only the UI reads. diff --git a/.changeset/run-quality-tabs-hidden-while-running.md b/.changeset/run-quality-tabs-hidden-while-running.md new file mode 100644 index 0000000..b63c8a8 --- /dev/null +++ b/.changeset/run-quality-tabs-hidden-while-running.md @@ -0,0 +1,14 @@ +--- +'@platforma-open/milaboratories.feature-integration.ui': patch +--- + +Run quality hides its tab strip while the run computes + +The tab set is derived from the baseline rung the run reports. Until the run +settles that rung is unknown, so the strip offered every plot and then dropped +the ones the served rung cannot draw — a reader could open a tab that then +stopped existing. + +The strip is now hidden while the block computes. The open view's body keeps +rendering and draws its own processing placeholder, so the page still shows the +run's progress. diff --git a/.changeset/undeclared-status-per-barcode.md b/.changeset/undeclared-status-per-barcode.md new file mode 100644 index 0000000..84ea55f --- /dev/null +++ b/.changeset/undeclared-status-per-barcode.md @@ -0,0 +1,29 @@ +--- +'@platforma-open/milaboratories.feature-integration.per-cell-metrics': minor +'@platforma-open/milaboratories.feature-integration.workflow': minor +'@platforma-open/milaboratories.feature-integration.model': minor +'@platforma-open/milaboratories.feature-integration': minor +--- + +The undeclared-barcode Status judges each barcode, not the whole sample + +The Status column read the sample's aggregate undeclared share, so it was one +word repeated down every row of a sample. A sample carrying one heavy +undeclared sequence among many light ones said nothing about which sequence to +look at. + +Status now reads each row's own share of its sample's pre-refine reads. It warns +above 1% and alerts above 5%. Those two numbers are operator-set and +overridable, not inherited: the field publishes 0.50/1.0 for a sample's +AGGREGATE undeclared share, and that line does not transfer to one sequence, +because an aggregate reaches 0.50 while no single sequence comes near it. + +The alert end changed direction with it. It compared for equality, which fired +only at exactly the error threshold and let every larger share read *warn* — the +worse finding being the one that never showed. Both ends now face the same way. + +The sample-level share keeps its column, renamed to "Sample Undeclared (%)", and +carries no status. + +Every column description in that table was rewritten to one instruction per +sentence, active voice, and short sentences. diff --git a/block/package.json b/block/package.json index fa5f776..d1b13e5 100644 --- a/block/package.json +++ b/block/package.json @@ -18,7 +18,7 @@ }, "scripts": { "build": "ts-builder build --target block-facade && block-tools pack", - "prepublishOnly": "block-tools publish -r s3://milab-euce1-prod-pkgs-s3-block-registry/pub/releases/?region=eu-central-1 --registry-serve-url https://blocks.pl-open.science", + "prepublishOnly": "block-tools publish --unstable -r s3://milab-euce1-prod-pkgs-s3-block-registry/pub/releases/?region=eu-central-1 --registry-serve-url https://blocks.pl-open.science", "do-pack": "shx rm -f package.tgz && pnpm pack && shx mv *.tgz package.tgz", "check": "ts-builder type-check --target block-facade" }, diff --git a/model/src/index.ts b/model/src/index.ts index d30a3dc..18b2703 100644 --- a/model/src/index.ts +++ b/model/src/index.ts @@ -6,15 +6,15 @@ import type { } from "@platforma-sdk/model"; import { BlockModelV3, - createPlDataTableStateV2, createPFrameForGraphs, + createPlDataTableStateV2, createPlDataTableV2, createPlDataTableV3, DataColumn, DataModelBuilder, + getAxisId, isPColumnSpec, parseResourceMap, - getAxisId, } from "@platforma-sdk/model"; import { assemblePattern, CELL_TAG, FEATURE_TAG, UMI_TAG, validatePattern } from "./pattern"; import { getPreset } from "./presets"; @@ -32,8 +32,8 @@ export type { PTableKey } from "@platforma-sdk/model"; // The reading's own defaults, in one exported map. Exported so a test can compare it against the other // two copies: this map is what a workflow-driven run is actually answered under, because -// verdict-args.lib.tengo emits every one of these flags UNCONDITIONALLY, substituting its own copy -// wherever the stored value is undefined. The argparse defaults in the Python never govern such a run. +// verdict-args.lib.tengo emits these flags UNCONDITIONALLY, substituting its own copy wherever the +// stored value is undefined. The argparse defaults in the Python never govern such a run. // `test/src/qcDefaults.test.ts` asserts each value against verdict-args.lib.tengo and the Python module // that owns it, the same way it asserts the QC lines below. export const VERDICT_DEFAULTS = { @@ -41,6 +41,8 @@ export const VERDICT_DEFAULTS = { countFloor: 4, // verdict.py BOUND_CUTOFF boundCutoff: 75, + // verdict.py DISTRIBUTION_BOUND_PROBABILITY. Both the default and the FLOOR: `args()` refuses below it. + boundProbability: 0.9, // combine.py DEFAULT_MIN_VOTERS minVotingCells: 1, // verdict.py DEFAULT_PANEL_MIN_MEMBERS. Gates rather than tunes: keep above the fifteen-tag cap of an @@ -64,8 +66,8 @@ export const QC_LINE_DEFAULTS = { readsPerCellWarn: 5000, aggregateBarcodeWarn: 0.05, aggregateBarcodeError: 1.0, - undeclaredBarcodeWarn: 0.5, - undeclaredBarcodeError: 1.0, + undeclaredBarcodeWarn: 0.01, + undeclaredBarcodeError: 0.05, usableReadWarn: 0.2, usableReadError: 0.0, } as const; @@ -215,6 +217,11 @@ export type TagCountBins = { * label column reaches p-frame surfaces only. Absent in full for a run that finished before this key. */ tagLabels?: Record; + /** + * The barcodes in the order the PANEL declares them, deduplicated on first appearance. A per-sample view + * reads in it, so that a barcode holds the same slot in every sample. + */ + tagOrder?: string[]; /** * The fit's two means and the background's share of cells, at the same (sample, tag) grain as the bins. * Here rather than in the p-frame beside them, so a grid of panels costs no driver query per panel. @@ -224,7 +231,15 @@ export type TagCountBins = { */ fitsBySample: Record< string, - Record + Record< + string, + { + backgroundMean: number; + signalMean: number; + backgroundWeight: number; + boundAtCount?: number | null; + } + > >; /** * The run's own two spreads, each on its own LINEAR edges: `score` and `referenceReading`. Linear, unlike @@ -505,8 +520,10 @@ type BlockDataV2 = Omit< | "distributionMinCells" | "countFloor" | "boundCutoff" + | "boundProbability" | "minVotingCells" | "minAgreement" + | "expectedBinderFraction" | "gateThreshold" | "grouping" | "contendingGroups" @@ -710,6 +727,13 @@ export const platforma = BlockModelV3.create(dataModel) ); if (data.countFloor < 0) throw new Error("The count floor cannot be negative"); + if ( + typeof data.boundProbability === "number" && + (data.boundProbability < VERDICT_DEFAULTS.boundProbability || data.boundProbability > 1) + ) + throw new Error( + `The fitted baseline's probability is at least ${VERDICT_DEFAULTS.boundProbability} and at most 1`, + ); if (data.boundCutoff < 0 || data.boundCutoff > 100) throw new Error("The bound cutoff is a score between 0 and 100"); if (data.minVotingCells < 1) throw new Error("At least one cell must vote"); @@ -722,6 +746,13 @@ export const platforma = BlockModelV3.create(dataModel) (data.minAgreement <= 0.5 || data.minAgreement > 1) ) throw new Error("The agreement floor is a share above 50% and at most 100%"); + // Strictly inside (0, 1). At either end the split hands every cell to one side and the fit silently + // falls back, so the run would record a fraction it never used. + if ( + typeof data.expectedBinderFraction === "number" && + (data.expectedBinderFraction <= 0 || data.expectedBinderFraction >= 1) + ) + throw new Error("The expected binder fraction is a share above 0% and below 100%"); // The cell condition GATES the fitted rung rather than tuning it, so it is a real population size. if (data.distributionMinCells < 1) throw new Error("A fitted baseline needs at least one cell to be fitted over"); @@ -855,6 +886,8 @@ export const platforma = BlockModelV3.create(dataModel) distributionMinCells: Math.round(data.distributionMinCells), countFloor: Math.round(data.countFloor), boundCutoff: data.boundCutoff, + boundProbability: data.boundProbability, + expectedBinderFraction: data.expectedBinderFraction, minVotingCells: Math.round(data.minVotingCells), // Off by default, and off means ABSENT. A minimum agreement of 0 passes every majority instead of skipping // the check. A gate of 0 sets aside every cell instead of gating none. diff --git a/model/src/types.ts b/model/src/types.ts index 1cfeb4a..6ecc1e6 100644 --- a/model/src/types.ts +++ b/model/src/types.ts @@ -90,6 +90,8 @@ export type BlockArgs = { distributionMinCells: number; // cells a sample needs before the rung may serve countFloor: number; // counts below this are not evidence of binding boundCutoff: number; // specificity score (0-100) at or above which a cell binds + boundProbability?: number; // the probability a count belongs to the signal component at or above which a cell binds + expectedBinderFraction?: number; // the share of cells expected to bind, seeding the fitted rung's split minVotingCells: number; // a verdict may rest on one cell and say so // Share (0-1) of answering cells the majority must reach. Off by default, and off means ABSENT rather // than zero: a floor of 0 passes every majority instead of skipping the check. @@ -188,6 +190,8 @@ export type BlockData = { distributionMinCells: number; countFloor: number; boundCutoff: number; + boundProbability?: number; + expectedBinderFraction?: number; minVotingCells: number; minAgreement?: number; gateThreshold?: number; diff --git a/pnpm-lock.yaml b/pnpm-lock.yaml index 167d8a2..8477d94 100644 --- a/pnpm-lock.yaml +++ b/pnpm-lock.yaml @@ -49,8 +49,8 @@ catalogs: specifier: 1.81.1 version: 1.81.1 '@platforma-sdk/workflow-tengo': - specifier: 6.8.2 - version: 6.8.2 + specifier: 6.8.3 + version: 6.8.3 ag-grid-enterprise: specifier: ~34.1.2 version: 34.1.2 @@ -254,7 +254,7 @@ importers: version: 2.3.1-131-main '@platforma-sdk/workflow-tengo': specifier: 'catalog:' - version: 6.8.2 + version: 6.8.3 devDependencies: '@platforma-sdk/tengo-builder': specifier: 'catalog:' @@ -1692,6 +1692,9 @@ packages: '@platforma-sdk/workflow-tengo@6.8.2': resolution: {integrity: sha512-AHR/Y+vbyfda17A7xY2uH9TlXiBpM8242tTO6+lj8phA4/KS6mez4GLu4ZEaASlZpX6YkQsUs5LUxYQU7pXeiQ==} + '@platforma-sdk/workflow-tengo@6.8.3': + resolution: {integrity: sha512-KImCzq7v/2qgH/PIBCTipV0B/ulLxDnDOpWOO+zs56+p/B2Dg25ZevlA2LzOFVk9SWx6ilXp4ilWx3K2HbccqQ==} + '@protobuf-ts/grpc-transport@2.11.1': resolution: {integrity: sha512-l6wrcFffY+tuNnuyrNCkRM8hDIsAZVLA8Mn7PKdVyYxITosYh60qW663p9kL6TWXYuDCL3oxH8ih3vLKTDyhtg==} peerDependencies: @@ -8736,6 +8739,14 @@ snapshots: '@platforma-open/milaboratories.software-ptexter': 1.2.4 '@platforma-open/milaboratories.software-small-binaries': 2.1.1 + '@platforma-sdk/workflow-tengo@6.8.3': + dependencies: + '@milaboratories/pframes-rs-wasip2': 1.1.56 + '@milaboratories/software-pframes-conv': 2.2.9 + '@platforma-open/milaboratories.software-ptabler': 2.1.8 + '@platforma-open/milaboratories.software-ptexter': 1.2.4 + '@platforma-open/milaboratories.software-small-binaries': 2.1.1 + '@protobuf-ts/grpc-transport@2.11.1(@grpc/grpc-js@1.14.4)': dependencies: '@grpc/grpc-js': 1.14.4 diff --git a/pnpm-workspace.yaml b/pnpm-workspace.yaml index 02d6def..a142b97 100644 --- a/pnpm-workspace.yaml +++ b/pnpm-workspace.yaml @@ -12,7 +12,7 @@ catalog: "@milaboratories/ts-builder": 1.6.2 "@milaboratories/ts-configs": 1.4.0 typescript: ~5.9.3 - "@platforma-sdk/workflow-tengo": 6.8.2 + "@platforma-sdk/workflow-tengo": 6.8.3 "@platforma-sdk/block-tools": 2.14.3 "@platforma-sdk/model": 1.81.1 "@platforma-sdk/ui-vue": 1.81.1 diff --git a/software/per-cell-metrics/src/emit_verdicts.py b/software/per-cell-metrics/src/emit_verdicts.py index 568b63f..0037bf9 100644 --- a/software/per-cell-metrics/src/emit_verdicts.py +++ b/software/per-cell-metrics/src/emit_verdicts.py @@ -89,14 +89,15 @@ Line, antigen_count_deciles, bin_values, - count_bin_edges, deciles_of, linear_bin_edges, + log1p_bin_edges, + log1p_edges_for, per_antigen_measures, per_tag_count_bins, reads_per_cell, sibling_disagreement, - status_for, + status_expr, usable_read_fraction, ) from qc_rows import ( @@ -123,12 +124,14 @@ from tag_distribution import ( DEFAULT_DISTRIBUTION_MIN_CELLS, TagFits, + bound_at_count, fit_tag_probabilities_by_pair, ) from verdict import ( BOUND_CUTOFF, DEFAULT_FLOOR, DEFAULT_PANEL_MIN_MEMBERS, + DISTRIBUTION_BOUND_PROBABILITY, Admissibility, Reference, ReferenceChoice, @@ -215,6 +218,18 @@ def main() -> None: "Required: nothing here picks a rung for a scientist who did not" ), ) + p.add_argument( + "--bound-probability", + type=float, + default=DISTRIBUTION_BOUND_PROBABILITY, + help=( + "how sure the fit must be before a cell counts as bound, as a probability that the cell's " + f"count came from the signal component. The lowest allowed value is " + f"{DISTRIBUTION_BOUND_PROBABILITY}, which is also the default. Lower is refused: a cell " + "holding none of the tag could then cross the line, and this run counts those cells by " + "arithmetic instead of checking each one, so the two halves would disagree" + ), + ) p.add_argument("--panel-min-members", type=int, default=DEFAULT_PANEL_MIN_MEMBERS) p.add_argument( "--distribution-min-cells", @@ -228,6 +243,8 @@ def main() -> None: ) p.add_argument("--min-voters", type=int, default=DEFAULT_MIN_VOTERS) p.add_argument("--min-agreement", type=float, default=DEFAULT_MIN_AGREEMENT) + # Roughly what share of cells are expected to bind an antigen. + p.add_argument("--initial-signal-weight", type=float, default=None) p.add_argument("--gate-threshold", type=int, default=None, help="set aside cells whose comparator reads above this") p.add_argument("--grouping", default=None, help="JSON: {'by':'tag'} or {'by':'property','column':...}") p.add_argument("--contending", default=None, help="JSON: groups of identities that contend, as a list of lists") @@ -276,6 +293,22 @@ def main() -> None: } add = functools.partial(_add, lines=lines) + # Must be above 0 and below 1. At either end the split hands every cell to one side, and the fit + # quietly falls back to splitting at the mean instead. The run would then report a value it never + # used, so refuse rather than let that happen silently. + # Absent is the ordinary case and means the pivot is derived from each tag's own counts. A value + # given is a statement about the experiment, and only then does the range apply. + if args.initial_signal_weight is not None and not 0.0 < args.initial_signal_weight < 1.0: + raise SystemExit( + f"the expected binder fraction is a share strictly between 0 and 1. Got {args.initial_signal_weight}." + ) + if args.bound_probability < DISTRIBUTION_BOUND_PROBABILITY: + raise SystemExit( + f"--bound-probability must be at least {DISTRIBUTION_BOUND_PROBABILITY}. Below that a cell " + f"holding none of a tag could be called bound. Most cells hold none of most tags, and the run " + f"counts them by arithmetic rather than checking each one, so those cells would be counted " + f"not-bound in one place and called bound in another. Got {args.bound_probability}." + ) if args.cutoff <= ANALYTIC_CUTOFF_BOUND: raise SystemExit( f"--cutoff must be strictly above {ANALYTIC_CUTOFF_BOUND:.4f}, the best score a zero count can reach. " @@ -465,6 +498,7 @@ def main() -> None: args.distribution_min_cells, floor=args.floor, reference_tags=reference_tags, + initial_signal_weight=args.initial_signal_weight, ) probabilities = _identity_probabilities(tag_fits, grouping) # A run where no tag fitted anywhere established no baseline. This is the one refusal that @@ -518,7 +552,7 @@ def main() -> None: non_reference = floored.filter(~pl.col("tag").is_in(list(reference_tags))) if reference_tags else floored identities = combine_tags_to_identities(non_reference, grouping) - states = read_states(identities, admissibility, args.cutoff) + states = read_states(identities, admissibility, args.cutoff, args.bound_probability) # The per-tag reading is diagnostic only: it compares each tag against the reference separately, # and no verdict is built from it. The measurement set carries it at both levels always, so where @@ -535,7 +569,10 @@ def main() -> None: else admissibility ) tag_states = read_states( - combine_tags_to_identities(non_reference, by_tag_grouping), tag_admissibility, args.cutoff + combine_tags_to_identities(non_reference, by_tag_grouping), + tag_admissibility, + args.cutoff, + args.bound_probability, ) # Which (sample, tag) pairs the reads actually carry, from the RAW counts. Never from `floored`: a @@ -684,6 +721,14 @@ def _admissibility(key: CellKey) -> str: ) _write_sorted(cell_counts, f"{prefix}_cell_counts.csv", ["sampleId", "cellId", "tag"]) + # The same (cell, tag) counts as the table above, but taken before the floor and including the + # control tag. + _write_sorted( + _listed(counts).select(["sampleId", "cellId", "tag", "umiCount"]), + f"{prefix}_cell_raw_counts.csv", + ["sampleId", "cellId", "tag"], + ) + cell_scalars = ( reference_frame.join(in_list, on=["sampleId", "cellId"], how="left") .with_columns(pl.col("inCellList").fill_null(unlisted_reads)) @@ -996,12 +1041,16 @@ def _disagreement_rates(samples_here: list[str]) -> dict[str, float | None]: # Two shares at two levels. `barcodeShare` is the row's own weight over every pre-refine # read of this sample. `readShare` is the SAMPLE's, computed once over every row of that # sample -- kept or elided by the tally's cap -- and carried on every row written. The - # status reads `readShare`; it is the barcode's and never the sample's, so it is written - # here rather than added to `rows` / `sample_report_rows`. - undeclared_status = status_for("undeclaredBarcodeShare", tally.share, lines) + # status reads `barcodeShare`, so it differs down the table; `readShare` carries no status. + # It is the barcode's and never the sample's, so it is written here rather than added to + # `rows` / `sample_report_rows`. barcode_share = ( (pl.col("totalWeight") / tally.total_weight) if tally.total_weight > 0 else pl.lit(None, pl.Float64) ) + # `status_expr`, not the scalar `status_for` in a loop: `heaviest`'s row cap is a parameter + # that accepts None, so a loop here is a loop over every distinct pre-refine sequence, and + # materialising the column to drive it undoes this stage's memory work. Both read the same + # lines and the same directions. undeclared_barcode_frames.append( tally.heaviest.select( pl.lit(sample, pl.String).alias("sampleId"), @@ -1009,8 +1058,7 @@ def _disagreement_rates(samples_here: list[str]) -> dict[str, float | None]: "totalWeight", barcode_share.cast(pl.Float64).alias("barcodeShare"), pl.lit(tally.share, pl.Float64).alias("readShare"), - pl.lit(None if undeclared_status is None else undeclared_status.value, pl.String).alias("status"), - ) + ).with_columns(status_expr("undeclaredBarcodeShare", pl.col("barcodeShare"), lines).alias("status")) ) # What correction then recovered. `tally.share` counts every pre-refine read the panel does # not declare; `1 - panelAssignedFraction` counts the reads refine-tags went on to drop, over @@ -1495,35 +1543,56 @@ def _disagreement_rates(samples_here: list[str]) -> dict[str, float | None]: # A run with no cell list bins every barcode, and `cellListSource` in the run meta says which case a # plot was drawn under. The edges are taken from the same filtered frame, so the shared domain ends # at the highest count among cells rather than among barcodes. + # The plot is binned over exactly the cells the fit ran on: one entry per cell in the sample, with a + # cell that read nothing counted as a zero. + fitted_bins = tag_fits.bins if tag_fits is not None else {} bin_counts = _listed(counts) - bin_edges = count_bin_edges(bin_counts) + if fitted_bins: + tag_bins: dict[str, dict[str, list[int]]] = {} + for (sample, tag), weights in fitted_bins.items(): + if weights: + tag_bins.setdefault(str(sample), {})[str(tag)] = weights + bin_edges = log1p_edges_for(max((len(w) for w in fitted_bins.values()), default=0)) + else: + # No fit ran, so there is no fitted population to bin. The observed readings still carry a shape + # worth showing, on the same equal-width edges. + observed = int(bin_counts["umiCount"].max() or 0) if bin_counts.height else 0 + bin_edges = log1p_bin_edges(observed) + tag_bins = per_tag_count_bins(bin_counts, bin_edges) # The fit's two means travel WITH the bins, at the same (sample, tag) grain. They are what a reader # judges the humps against, so reaching them through the p-frame beside this would mean a driver # query per panel of the grid. Absent under a declared baseline, which fits nothing. fits_by_sample: dict[str, dict[str, dict[str, float]]] = {} + fit_curves = tag_fits.curves if tag_fits is not None else {} for (sample, tag), b in (tag_fits.backgrounds if tag_fits is not None else {}).items(): + # The count at which the run's bound probability starts calling a cell bound, resolved from this + # pair's own scored curve. + curve = fit_curves.get((sample, tag)) + crossing = None if curve is None else bound_at_count(curve, args.bound_probability) fits_by_sample.setdefault(str(sample), {})[str(tag)] = { "backgroundMean": float(b.mean), "signalMean": float(b.signal_mean), "backgroundWeight": float(b.weight), + "boundAtCount": crossing, } with open(f"{prefix}_qc_tag_bins.json", "w") as out: json.dump( { "edges": bin_edges, - "bySample": per_tag_count_bins(bin_counts, bin_edges), + "bySample": tag_bins, "fitsBySample": fits_by_sample, # The same names `result_tag_labels.csv` carries, so a tag reads under one name on the plots # and in the reagent table. The plots are drawn from this JSON rather than from a p-frame, and # a label column reaches only p-frame surfaces, so without this every panel title is a barcode. "tagLabels": tag_names, + # The panel's OWN row order. + "tagOrder": list(dict.fromkeys(panel["tag"].to_list())), # The run's score spread and its reference readings, on their own linear edges. Each is one # distribution for the whole run, pooled across samples and narrowed to the cell list, so it # carries the same population as the count bins beside it. "spreads": spread_bins, }, out, - indent=2, sort_keys=True, ) @@ -1573,11 +1642,13 @@ def _disagreement_rates(samples_here: list[str]) -> dict[str, float | None]: "cellsAnalysed": len(analysed_cells), "floor": args.floor, "cutoff": args.cutoff, + "boundProbability": args.bound_probability, "minVoters": args.min_voters, "minAgreement": args.min_agreement, "gateThreshold": args.gate_threshold, "panelMinMembers": args.panel_min_members, "distributionMinCells": args.distribution_min_cells, + "initialSignalWeight": args.initial_signal_weight, # Per (sample, tag), and only where that rung was asked for: which tags could not be fitted, # and why. A tag missing here fitted. The reader needs both halves to tell a panel that mostly # worked from one that mostly did not. diff --git a/software/per-cell-metrics/src/qc_measures.py b/software/per-cell-metrics/src/qc_measures.py index 989ed35..ba177b1 100644 --- a/software/per-cell-metrics/src/qc_measures.py +++ b/software/per-cell-metrics/src/qc_measures.py @@ -375,6 +375,9 @@ class Measurement: # `Measurement`: that status is the barcode's, never a sample's, so it is computed and carried where # the barcode rows are, in emit_verdicts.py, and reaches `status_for` under this id. It is the one # exception to "every line backs a declared measurement", and a test names it. +# +# It reads the ROW's own share, `barcodeShare`, not the sample-level `readShare` the id is named for. +# The sample-level share keeps its column and carries no status. DEFAULT_LINES: dict[str, Line] = { # Both thresholds step the same way. This is the line with a real gradient at the far end. "cellBarcodeValidFraction": Line(warn=0.75, error=0.50), @@ -383,9 +386,11 @@ class Measurement: # Published values for the aggregate-barcode read fraction: warn above 0.05, error at total # failure (1.0). "aggregateBarcodeFraction": Line(warn=0.05, error=1.0), - # Published values for the undeclared-barcode read fraction, read direct rather than as a - # complement: warn above 0.50, error at total failure (1.0). - "undeclaredBarcodeShare": Line(warn=0.5, error=1.0), + # ONE BARCODE's share of its sample's pre-refine reads: warn above 0.01, alert above 0.05. + # Operator-set, not inherited. The field publishes 0.50/1.0 for a sample's AGGREGATE undeclared + # share, and that line does not transfer to a single sequence: the aggregate reaches 0.50 while no + # single sequence comes near it. Needs an atom on 315 before it can be called inherited. + "undeclaredBarcodeShare": Line(warn=0.01, error=0.05), # Published values for the usable antigen-read fraction: warn below 0.20, error at total # failure (0.0). "usableReadFraction": Line(warn=0.20, error=0.0), @@ -413,7 +418,9 @@ class Measurement: # inherited share sits at either "at least" or "at most" with error at the catastrophe end, and # this is one of the two upward-facing members of that set. "aggregateBarcodeFraction": ("at-most", "alerting-at"), - "undeclaredBarcodeShare": ("at-most", "alerting-at"), + # Both ends face the same way, unlike the four inherited shares: this line alerts ABOVE its error + # threshold rather than at a catastrophe value, so `alerting-at` would fire only at exactly 0.05. + "undeclaredBarcodeShare": ("at-most", "at-most"), # Error at total failure (`alerting-at` 0.0), the downward-facing member of that same set. "usableReadFraction": ("at-least", "alerting-at"), } @@ -428,6 +435,19 @@ def _breaches(value: float, threshold: float, comparison: str) -> bool: return value == threshold +def _breaches_expr(value: pl.Expr, threshold: float, comparison: str) -> pl.Expr: + """`_breaches` over a column. Every branch mirrors the scalar above, line for line. + + `test_status_expr_agrees_with_status_for` runs the two against one another over the boundary + values and every registered measurement, so a branch changed on one side alone fails there. + """ + if comparison == "at-least": + return value < threshold + if comparison == "at-most": + return value > threshold + return value == threshold + + _ORDINAL = {Status.OK: 0, Status.WARN: 1, Status.ALERT: 2} _DEFERRED: frozenset[str] = frozenset(m.id for m in MEASUREMENTS if m.deferred_reason) @@ -475,6 +495,52 @@ def status_for(measurement: str, value: float | None, lines: dict[str, Line]) -> return Status.OK +def status_expr(measurement: str, value: pl.Expr, lines: dict[str, Line]) -> pl.Expr: + """`status_for` over a column, returning the status string or null. + + For a per-row status on a frame with no bound on its height. The undeclared-barcode table is the + caller: its row cap is a parameter that accepts `None`, so a Python loop there is a loop over every + distinct pre-refine sequence -- 10.2M per sample on a measured 44-sample run -- and materialising + the column to drive it undoes the memory work this stage carries. + + Thresholds and directions are read from the SAME `lines` and `_COMPARISON` this module's scalar + reads. Only the evaluator differs, and a test pins the two together. + """ + null = pl.lit(None, pl.String) + if measurement in _DEFERRED: + return null + # `is_computed` over a column: a null, a NaN or an infinity is not a number. Kleene `&` makes a + # null value read False here rather than propagating a null into the branch below. + computed = (value.is_not_null() & value.is_finite()).fill_null(False) # noqa: FBT003 + if measurement in _CATEGORICAL: + return ( + pl.when(~computed) + .then(null) + .when(value == 0) + .then(pl.lit(Status.ALERT.value, pl.String)) + .otherwise(pl.lit(Status.OK.value, pl.String)) + ) + if measurement not in lines: + return null + line = lines[measurement] + warn_comparison, error_comparison = _COMPARISON[measurement] + alerts = ( + _breaches_expr(value, line.error, error_comparison) + if line.error is not None and error_comparison is not None + else pl.lit(False) # noqa: FBT003 + ) + # Error first, exactly as the scalar orders it: a value past both boundaries reads alert. + return ( + pl.when(~computed) + .then(null) + .when(computed & alerts) + .then(pl.lit(Status.ALERT.value, pl.String)) + .when(computed & _breaches_expr(value, line.warn, warn_comparison)) + .then(pl.lit(Status.WARN.value, pl.String)) + .otherwise(pl.lit(Status.OK.value, pl.String)) + ) + + def roll_up(readings: list[Reading]) -> Coverage: """The worst status among those that carry one, plus coverage. @@ -812,44 +878,64 @@ def deciles_of(values: np.ndarray) -> pl.DataFrame: ) -# How many log-spaced buckets a count distribution is drawn in. Enough bars for two humps to read -# apart at thumbnail size, few enough that a sparse tag does not dissolve into single-cell spikes. +# How many buckets the count distributions used to be drawn in, back when their edges were integers. +# Kept only because `linear_bin_edges` uses it as its default; the count distributions do not. COUNT_BIN_COUNT = 24 -def count_bin_edges(counts: pl.DataFrame) -> list[float]: - """Log-spaced bin edges spanning every count in the frame, shared by every plot drawn from it. +# The width of one bar, measured in log1p units. A fixed width is what makes every bar the same size. +# +# The source paper uses 0.075. This is deliberately coarser: at 0.075 a real tag came back with 75 of its +# 97 bins empty, because whole-number counts land on scattered points once you take the log. The paper +# lives with those gaps by drawing a smooth density curve over them, which these plots cannot do. 0.2 +# keeps the bars equal and still readable at thumbnail size. +LOG1P_BIN_WIDTH = 0.2 + + +def log1p_bin_edges(top: int, width: float = LOG1P_BIN_WIDTH) -> list[float]: + """Bin edges spanning 0 to `top` that all draw the SAME WIDTH. `[]` if there is nothing to span. - The caller passes the counts of the CELL LIST where one arrived, so the domain ends at the - highest count among cells. Observed barcodes outnumber cells by one to two orders of magnitude. + Returned as counts, at `expm1(k * width)`, because counts are what the plot takes. The plot's axis is + log1p, so an edge at `expm1(k * width)` lands at `k * width` on screen -- evenly spaced. - ONE edge set for the whole run, not one per tag. A reader judges whether a tag's counts fall into - two separated humps by scanning a grid of tags side by side, and per-tag edges would rescale - every panel to its own range, so a tag whose counts span 1-4 and one spanning 1-4000 would draw - identical pictures. + Edges start at 0 so the zeros the fit ran over get a bar of their own. Those zeros are most of the + background; drop them and the plot shows one decaying hump whatever the fit found. The last edge is + the first step above `top`, so the largest count has a bin to land in. - Log-spaced because UMI counts per cell span orders of magnitude: on a linear axis the ambient - population occupies one bar and everything above it is empty. + ONE edge set for the whole run, not one per tag, so a reader can compare a grid of tags side by side. - Edges start at 1, the smallest count that exists -- a row is only written for an observed - reading, so zero never appears. `[]` where the frame holds no counts at all. + The cost: THE EDGES ARE NOT WHOLE NUMBERS. Counts are, so consecutive integers sit further apart than + one bin until about count 13, and the low end comes out as separated bars with empty gaps between + them. Whole-number edges avoided that, which is why they were used here before -- but they bought it + with the per-bar division above. The source paper's figures show the same gaps. + + Nothing needs to be told the bin count: `bin_values`, `per_tag_count_bins` and the UI all read it + from `len(edges) - 1`. """ - if counts.height == 0: + if top < 1 or width <= 0.0: return [] - top = float(counts["umiCount"].max() or 0) - if top < 1: + return log1p_edges_for(int(np.floor(np.log1p(top) / width)) + 1, width) + + +def log1p_edges_for(steps: int, width: float = LOG1P_BIN_WIDTH) -> list[float]: + """The first `steps` bins of the same grid, as counts. `[]` for a non-positive step or width. + + Separate from `log1p_bin_edges` because the edges are ABSOLUTE: they sit at `expm1(k * width)` + whatever the data holds, so a bin count alone identifies them. Each (sample, tag) bins against its + own range, and the plot then needs one edge list long enough for the widest of them -- which is a + bin count, not a frame to re-scan. + """ + if steps < 1 or width <= 0.0: return [] - # `top` lands on the last edge, so the largest count falls inside the last bin rather than - # outside every bin. - return [float(x) for x in np.geomspace(1.0, max(top, 2.0), COUNT_BIN_COUNT + 1)] + return [float(np.expm1(k * width)) for k in range(steps + 1)] def linear_bin_edges(values: np.ndarray, count: int = COUNT_BIN_COUNT) -> list[float]: """Evenly spaced bin edges spanning `values`. `[]` where there are none. - Linear, unlike `count_bin_edges`: a specificity score is a 0-100 scale and a reference reading is - read against a gate a scientist types in the same units, so a log axis would put the number they - are choosing somewhere they cannot find it. + Evenly spaced, unlike the count distributions' log1p edges. A specificity score is a 0-100 scale, + and a reference reading is judged against a threshold the scientist types in the same units, so a + log axis would put the number they are choosing somewhere they cannot find it. """ if values.size == 0: return [] diff --git a/software/per-cell-metrics/src/qc_rows.py b/software/per-cell-metrics/src/qc_rows.py index 52d7da5..76618b9 100644 --- a/software/per-cell-metrics/src/qc_rows.py +++ b/software/per-cell-metrics/src/qc_rows.py @@ -289,10 +289,10 @@ def sample_report_rows(sample: str, rows: list[QcRow]) -> tuple[list[dict], Cove # One row per (sampleId, tag) the pre-refine pass saw and the sample's panel does not declare. # Two shares, at two levels, and neither substitutes for the other. `barcodeShare` is this one # sequence's weight over every pre-refine read of its sample. `readShare` and `status` are the -# SAMPLE's undeclared-read share, repeated on every one of that sample's rows: the field publishes a -# line for the share of a sample's reads landing in barcodes nobody declared, and that status is the -# barcode's, never the sample's -- so it is computed at the sample and carried on the barcode's own -# row. Usually there are no rows for a sample at all, which is the wanted outcome. +# SAMPLE's undeclared-read share, repeated on every one of that sample's rows and carrying no status. +# `status` reads `barcodeShare`, so it is the row's own and differs down the table. That status is the +# barcode's, never the sample's, and never rolls into a sample's. Usually there are no rows for a +# sample at all, which is the wanted outcome. _UNDECLARED_BARCODE_SCHEMA = { "sampleId": pl.String, "tag": pl.String, diff --git a/software/per-cell-metrics/src/tag_distribution.py b/software/per-cell-metrics/src/tag_distribution.py index adf5565..8faa6be 100644 --- a/software/per-cell-metrics/src/tag_distribution.py +++ b/software/per-cell-metrics/src/tag_distribution.py @@ -17,6 +17,13 @@ threshold. Nothing downstream thresholds: it scores a reading against a comparator count. So the split point identifies which readings are background, and the comparator is the middle of those. +**Where the fit starts.** While the EM runs, the two components are interchangeable -- the paper says +so, and both it and this module decide afterwards which is which, by median. So the starting point does +not change what the components mean; it changes which answer the EM settles on. On a mostly-background +population that is the difference between a background weight near 0.8 and one near 0.95. The starting +split is the paper's, `DEFAULT_INITIAL_SIGNAL_WEIGHT` above. It used to be the median, which is a choice +`what-plays-the-baseline` never actually specified. + **No normalization.** The paper normalizes by each cell's UMI total. Every reading here is a raw integer UMI count, and a normalized comparator would be the only non-count in the pipeline. The cost is that a cell sequenced twice as deeply contributes a reading twice as large; the split is @@ -34,6 +41,7 @@ import numpy as np import polars as pl from panel import ANY_SAMPLE +from qc_measures import bin_values, log1p_bin_edges, log1p_edges_for from scipy.special import logsumexp from scipy.stats import nbinom @@ -85,6 +93,10 @@ class TagFit(NamedTuple): NO_FIT = "no two-component fit could be computed for this tag" +# Roughly what share of cells are expected to bind an antigen. The source paper's own initial weight, +# and an OVERRIDE here rather than the default. +DEFAULT_INITIAL_SIGNAL_WEIGHT = 0.1 + # The share of the highest counts dropped before the fit, so that a handful of very high readings # cannot drag the signal component's mean up and pull the boundary with it. The dropped cells still # get a probability -- they are the most bound cells in the sample. @@ -159,7 +171,59 @@ class _Mixture(NamedTuple): signal: int -def _fit_two_component_nb(counts: np.ndarray) -> _Mixture | None: +# @TODO: Review this function in more detail looking for logic flaws that might have been missed +def _scan_pivot(x: np.ndarray) -> float | None: + """The split point that best explains the counts as two negative binomials, or None. + + Scores every cut of the counts into `[0, t]` and the rest, and returns the `t` that scores best. + None where the counts hold fewer than two distinct values, or where no cut scores finitely. + """ + # The histogram, not the cells: everything below is sized by the number of DISTINCT counts, which is + # why the scan costs milliseconds against the EM's seconds and does not grow with the sample. + vals, mult = np.unique(x, return_counts=True) + k = vals.size + if k < 2: + return None + + # Prefix sums of the cell count, of the counts, and of their squares. They make each cut's two-sided + # weight, mean and variance a subtraction instead of a pass over the data. + n = float(x.size) + cum_w = np.cumsum(mult).astype(float) + cum_s = np.cumsum(vals * mult) + cum_q = np.cumsum((vals**2) * mult) + + # One entry per cut, and `[:-1]` is the cut set: cutting at the highest count leaves the upper side + # empty. Variance by `E[x^2] - E[x]^2`, so it can land just below zero on a near-constant side. + w_lo = cum_w[:-1] + w_hi = n - w_lo + m_lo = cum_s[:-1] / w_lo + m_hi = (cum_s[-1] - cum_s[:-1]) / w_hi + v_lo = cum_q[:-1] / w_lo - m_lo**2 + v_hi = (cum_q[-1] - cum_q[:-1]) / w_hi - m_hi**2 + m_lo = np.maximum(m_lo, _MIN_COMPONENT_MEAN) + m_hi = np.maximum(m_hi, _MIN_COMPONENT_MEAN) + + # Dispersions through `_nb_size` rather than inline, so each side inherits the same Poisson-limit + # fallback the EM's own components get and the two cannot disagree about what a dispersion is. + s_lo = np.array([_nb_size(m, v) for m, v in zip(m_lo, np.maximum(v_lo, 0.0), strict=True)]) + s_hi = np.array([_nb_size(m, v) for m, v in zip(m_hi, np.maximum(v_hi, 0.0), strict=True)]) + + # The score, one row per cut and one column per distinct count. It is the CLASSIFICATION + # log-likelihood -- every count read against the side its cut assigns it, weighted by that side's + # share and by how many cells hold the count. + below = vals[None, :] <= vals[:-1, None] + lo = _nb_logpmf(vals[None, :], m_lo[:, None], s_lo[:, None]) + np.log(w_lo / n)[:, None] + hi = _nb_logpmf(vals[None, :], m_hi[:, None], s_hi[:, None]) + np.log(w_hi / n)[:, None] + scored = np.where(below, lo, hi) * mult[None, :] + + # A cut where any count underflowed is dropped rather than left to win the argmax on a -inf. + total = np.where(np.all(np.isfinite(scored), axis=1), scored.sum(axis=1), -np.inf) + if not np.any(np.isfinite(total)): + return None + return float(vals[int(np.argmax(total))]) + + +def _fit_two_component_nb(counts: np.ndarray, initial_signal_weight: float | None = None) -> _Mixture | None: """A two-component negative binomial fitted to `counts`, or None where none exists. The method: fit a two-component negative binomial mixture and label the higher-median component @@ -179,13 +243,18 @@ def _fit_two_component_nb(counts: np.ndarray) -> _Mixture | None: if np.unique(x).size < 2: return None - # Split at the median to start. Two components initialised on the same statistics never separate, - # and the median is the one split point that is always available. - pivot = float(np.median(x)) + # Where the scientist stated a share, split the counts so that `initial_signal_weight` of them sit + # above the pivot. Otherwise derive the split from the counts themselves. Each side then seeds one + # component, and the two side sizes are the starting weights. + scanned = _scan_pivot(x) if initial_signal_weight is None else None + if scanned is not None: + pivot = scanned + else: + pivot = float(np.quantile(x, 1.0 - (initial_signal_weight or DEFAULT_INITIAL_SIGNAL_WEIGHT))) low, high = x[x <= pivot], x[x > pivot] if low.size == 0 or high.size == 0: - # A median equal to the maximum puts everything in one half. Fall back to splitting at the - # mean, which differs from the median exactly when the counts are skewed. + # The quantile landed on the maximum, so one side got every cell. Split at the mean instead: + # on a skewed count distribution the mean sits below the upper quantile, so it still divides. pivot = float(np.mean(x)) low, high = x[x <= pivot], x[x > pivot] if low.size == 0 or high.size == 0: @@ -269,6 +338,7 @@ def fit_tag_probabilities( counts: np.ndarray, min_cells: int = DEFAULT_DISTRIBUTION_MIN_CELLS, scored: np.ndarray | None = None, + initial_signal_weight: float | None = None, ) -> TagFit: """One tag's counts across one sample's cells, as a probability of binding per cell. @@ -308,7 +378,7 @@ def fit_tag_probabilities( if fitted_on.size == 0: return TagFit(None, NO_FIT, n) - fit = _fit_two_component_nb(fitted_on) + fit = _fit_two_component_nb(fitted_on, initial_signal_weight) if fit is None: return TagFit(None, NO_FIT, n) probabilities = _signal_probability(scored, fit) @@ -334,9 +404,14 @@ class TagFits(NamedTuple): probabilities: pl.DataFrame reasons: dict[tuple[str, str], str] + # Each pair's fitted cells, already binned for the plot: one entry per cell in the sample, with a + # cell that read nothing counted as a zero. + bins: dict[tuple[str, str], list[int]] # One entry per pair that fitted, on the same condition as a contribution to `probabilities`. A # pair in `reasons` is absent here. backgrounds: dict[tuple[str, str], Background] = {} + # Each fitted pair's distinct counts, ascending, with the probability the fit gave each one. + curves: dict[tuple[str, str], tuple[np.ndarray, np.ndarray]] = {} _PROB_SCHEMA = {"sampleId": pl.String, "cellId": pl.String, "tag": pl.String, "pBound": pl.Float64} @@ -349,6 +424,7 @@ def fit_tag_probabilities_by_pair( min_cells: int = DEFAULT_DISTRIBUTION_MIN_CELLS, floor: int = 0, reference_tags: Collection[str] = (), + initial_signal_weight: float | None = None, ) -> TagFits: """One fit per (sample, tag) the panel declares, scored per cell. @@ -394,7 +470,9 @@ def fit_tag_probabilities_by_pair( by_sample = {s: f.select("cellId") for (s,), f in universe.group_by("sampleId")} frames: list[pl.DataFrame] = [] reasons: dict[tuple[str, str], str] = {} + binned: dict[tuple[str, str], list[int]] = {} backgrounds: dict[tuple[str, str], Background] = {} + curves: dict[tuple[str, str], tuple[np.ndarray, np.ndarray]] = {} for sample, tag in _declared_pairs(panel, sorted(by_sample)): sample_cells = by_sample.get(sample) if sample_cells is None: @@ -408,13 +486,23 @@ def fit_tag_probabilities_by_pair( ).with_columns(pl.col("umiCount").fill_null(0)) raw = dense["umiCount"].to_numpy() + # Binned before the fit is attempted, so a pair that established nothing still carries the + # distribution a reader would judge that outcome against. `raw` dies with the iteration; only + # the bins outlive it. + binned[(sample, tag)] = ( + bin_values(raw, log1p_bin_edges(int(raw.max())) or log1p_edges_for(1)) if raw.size else [] + ) scored = raw if floor <= 0 or tag in exempt else np.where(raw < floor, 0, raw) - fit = fit_tag_probabilities(raw, min_cells, scored=scored) + fit = fit_tag_probabilities(raw, min_cells, scored=scored, initial_signal_weight=initial_signal_weight) if fit.probabilities is None: reasons[(sample, tag)] = fit.reason or NO_FIT continue if fit.background is not None: backgrounds[(sample, tag)] = fit.background + # Built from `scored`, not `raw`. The probabilities were computed on the scored values, so a + # reading the floor zeroed carries the probability of 0, not of the count it originally held. + distinct, first = np.unique(scored, return_index=True) + curves[(sample, tag)] = (distinct, np.asarray(fit.probabilities, dtype=float)[first]) frames.append( dense.select("cellId") .with_columns( @@ -426,7 +514,29 @@ def fit_tag_probabilities_by_pair( ) probabilities = pl.concat(frames) if frames else pl.DataFrame(schema=_PROB_SCHEMA) - return TagFits(probabilities, reasons, backgrounds) + return TagFits(probabilities, reasons, binned, backgrounds, curves) + + +def bound_at_count(curve: tuple[np.ndarray, np.ndarray], probability_cutoff: float) -> int | None: + """The lowest count at or above which every count is called bound, or None if there is none. + + `curve` is one `TagFits.curves` entry: ascending distinct counts, and the probability each was given. + + Note it looks for the lowest count from which the cutoff holds all the way up, NOT simply the first + count to cross it. Usually these are the same, because the probability normally rises with the count. + + They differ when a fit comes back inverted. Components are labelled by median, so a pair can end up + with its "signal" component sitting BELOW its background. The probability then FALLS as the count + rises, and it is the LOW counts that cross the cutoff. Taking the first crossing would report a + threshold of 1, which the run does not apply. Requiring the cutoff to hold all the way up returns + None instead, which says plainly that this fit has no count above which it calls a cell bound. + """ + counts, probabilities = curve + if counts.size == 0 or counts.size != probabilities.size: + return None + missed = np.flatnonzero(probabilities < probability_cutoff) + start = 0 if missed.size == 0 else int(missed[-1]) + 1 + return int(counts[start]) if start < counts.size else None _EMPTY = np.zeros(0, dtype=np.int64) diff --git a/software/per-cell-metrics/src/verdict.py b/software/per-cell-metrics/src/verdict.py index a985d3d..628aaa7 100644 --- a/software/per-cell-metrics/src/verdict.py +++ b/software/per-cell-metrics/src/verdict.py @@ -302,13 +302,11 @@ def gate_cells( BETA_X, BETA_A_OFFSET, BETA_B_OFFSET = 0.925, 1, 3 BOUND_CUTOFF = 75.0 -# The population rung's own call, and not this block's to move: under a fitted distribution -# a cell reads *bound* at 0.9 or above, where 0.9 is the probability its count belongs to -# the signal component. The score below is the declared reagent's rule and does not apply -# here -- each baseline brings its own rule, and a run selects one baseline. +# The population rung's line: under a fitted distribution a cell reads *bound* where the probability +# its count belongs to the signal component reaches this. The score below is the declared reagent's +# rule and does not apply here -- each baseline brings its own rule, and a run selects one baseline. # -# NOT A SETTING, and it must not become one. It comes from the literature, so a dial would -# only produce runs that cannot be compared against the work the method came from. +# The DEFAULT, and the FLOOR of what a run may ask for DISTRIBUTION_BOUND_PROBABILITY = 0.9 @@ -503,7 +501,12 @@ def _comparator(key: tuple[str, str], identity: str, admissibility: Admissibilit return reference.get(key) -def read_states(identities: pl.DataFrame, admissibility: Admissibility, cutoff: float) -> pl.DataFrame: +def read_states( + identities: pl.DataFrame, + admissibility: Admissibility, + cutoff: float, + probability_cutoff: float = DISTRIBUTION_BOUND_PROBABILITY, +) -> pl.DataFrame: """Give every (cell, identity) row a state. Two routes to UNRELIABLE, both recorded in `unreliableReason`: the cell has no comparator, @@ -538,7 +541,7 @@ def read_states(identities: pl.DataFrame, admissibility: Admissibility, cutoff: df = df.with_columns(pl.Series("_pBound", called, dtype=pl.Float64)).with_columns( pl.when(pl.col("unreliableReason").is_not_null()) .then(pl.lit(State.UNRELIABLE.value)) - .when(pl.col("_pBound") >= DISTRIBUTION_BOUND_PROBABILITY) + .when(pl.col("_pBound") >= probability_cutoff) .then(pl.lit(State.BOUND.value)) .otherwise(pl.lit(State.NOT_BOUND.value)) .alias("state") diff --git a/software/per-cell-metrics/test/test_emit_verdicts.py b/software/per-cell-metrics/test/test_emit_verdicts.py index ccf424d..7180388 100644 --- a/software/per-cell-metrics/test/test_emit_verdicts.py +++ b/software/per-cell-metrics/test/test_emit_verdicts.py @@ -93,6 +93,7 @@ def test_writes_every_artifact(bed): "result_verdicts.csv", "result_set_counts.csv", "result_cell_counts.csv", + "result_cell_raw_counts.csv", "result_cell_scalars.csv", "result_offered.csv", "result_identity_labels.csv", @@ -586,7 +587,7 @@ def test_undeclared_barcode_table_is_keyed_by_sequence_with_the_samples_share(be # is the only undeclared sequence; the next test separates them. assert row["readShare"] == pytest.approx(10 / 60) assert row["barcodeShare"] == pytest.approx(10 / 60) - assert row["status"] == "OK" # 10/60 is well under the 0.50 warn line + assert row["status"] == "alert" # this sequence is 16.7% of the sample, above the 5% alert line def test_each_undeclared_barcode_carries_its_own_share_beside_the_samples(bed): @@ -612,11 +613,25 @@ def test_the_own_share_denominator_is_every_pre_refine_read_not_the_undeclared_o assert sum(t["barcodeShare"].to_list()) == pytest.approx(0.10) -def test_undeclared_barcode_share_warns_above_half(bed): - _write_raw_feature_counts(bed, [("S1", "AAAA", 10), ("S1", "ZZZZ", 15)]) +def test_undeclared_barcode_status_warns_above_one_percent(bed): + _write_raw_feature_counts(bed, [("S1", "AAAA", 980), ("S1", "ZZZZ", 20)]) _run(bed, *BASE, "--raw-feature-counts", "raw_feature_counts.csv") t = pl.read_csv(bed / "result_undeclared_barcodes.csv") - assert t.row(0, named=True)["status"] == "warn" + assert t.row(0, named=True)["status"] == "warn" # 20/1000 = 2% + + +def test_undeclared_barcode_status_is_the_rows_own_not_the_samples(bed): + # The status reads `barcodeShare`, so two rows of one sample can read differently. Read from + # `readShare` it would be one word repeated, and a sample carrying one heavy sequence among many + # light ones would say nothing about which sequence to look at. + _write_raw_feature_counts(bed, [("S1", "AAAA", 900), ("S1", "ZZZZ", 80), ("S1", "YYYY", 20)]) + _run(bed, *BASE, "--raw-feature-counts", "raw_feature_counts.csv") + t = pl.read_csv(bed / "result_undeclared_barcodes.csv") + by_tag = {r["tag"]: r for r in t.iter_rows(named=True)} + assert by_tag["ZZZZ"]["status"] == "alert" # 8%, above the 5% alert line + assert by_tag["YYYY"]["status"] == "warn" # 2%, above the 1% warn line + # One sample-level number on both rows, and it is no longer what the status reads. + assert [r["readShare"] for r in by_tag.values()] == [pytest.approx(0.10)] * 2 def test_undeclared_barcode_share_alerts_when_every_read_is_undeclared(bed): @@ -690,6 +705,7 @@ def test_rows_are_sorted_on_a_bed_wide_enough_for_order_to_show(bed): for name, keys in ( ("result_verdicts.csv", ["setId", "identity"]), ("result_cell_counts.csv", ["sampleId", "cellId", "tag"]), + ("result_cell_raw_counts.csv", ["sampleId", "cellId", "tag"]), ("result_offered.csv", ["sampleId", "identity"]), ("result_tag_identity.csv", ["tag", "identity"]), ): @@ -3718,3 +3734,57 @@ def test_the_rescued_share_says_why_where_no_pre_refine_pass_reached_the_run(bed entry = next(m for m in report["S1"]["measurements"] if m["id"] == "refineRescuedShare") assert entry["value"] is None assert entry["reason"] + + +# --- the pre-floor per-cell per-tag table ------------------------------------------------------- +# +# The bed's non-reference readings are 500 and 600, so a floor has to sit between them to bite. 550 is +# chosen for that and for nothing else: it zeroes exactly one reading, which is what these tests read. + + +@pytest.mark.slow +def test_raw_counts_are_taken_before_the_floor(bed): + """The point of the second table: the floor is not visible in it. + + `result_cell_counts.csv` is evidence of binding, so `apply_floor` has zeroed every count below the + minimum there. `result_cell_raw_counts.csv` is capture, so the same reading keeps the value the reads + carried. Asserted against each other, because either table read alone cannot show that a count was + floored -- afterwards a floored count and a true zero are the same number. + """ + assert _run(bed, *BASE, "--floor", "550").returncode == 0 + floored = pl.read_csv(bed / "result_cell_counts.csv", infer_schema_length=0) + raw = pl.read_csv(bed / "result_cell_raw_counts.csv", infer_schema_length=0) + + assert raw.columns == ["sampleId", "cellId", "tag", "umiCount"] + floored_values = sorted(int(v) for v in floored["umiCount"].to_list()) + raw_values = sorted(int(v) for v in raw["umiCount"].to_list()) + + # The 500 was zeroed for the verdicts and kept here. + assert 0 in floored_values, "the floor bit nothing, so this test proves nothing" + assert 500 in raw_values + assert 0 not in raw_values, "a pre-floor table holds no manufactured zeros" + + +@pytest.mark.slow +def test_raw_counts_keep_the_comparator_the_floored_table_drops(bed): + # `cell_counts` is written from the NON-reference readings, and `apply_floor` exempts the comparator + # from the floor besides. The capture table applies neither rule, so a cell's tags sum to what that + # cell held -- which is what a composition plot needs. + assert _run(bed, *BASE).returncode == 0 + floored = pl.read_csv(bed / "result_cell_counts.csv", infer_schema_length=0) + raw = pl.read_csv(bed / "result_cell_raw_counts.csv", infer_schema_length=0) + assert "CTRL" not in set(floored["tag"].to_list()) + assert "CTRL" in set(raw["tag"].to_list()) + + +@pytest.mark.slow +def test_raw_counts_do_not_move_when_the_floor_does(bed): + # The floor reaches the verdicts and must not reach this table. Byte equality on one side, and a + # required DIFFERENCE on the other: without the second half the test would pass on a run where the + # floor changed nothing at all. + _run(bed, *BASE, "--floor", "1") + raw_low = (bed / "result_cell_raw_counts.csv").read_bytes() + floored_low = (bed / "result_cell_counts.csv").read_bytes() + _run(bed, *BASE, "--floor", "550") + assert (bed / "result_cell_raw_counts.csv").read_bytes() == raw_low + assert (bed / "result_cell_counts.csv").read_bytes() != floored_low, "the floor moved nothing" diff --git a/software/per-cell-metrics/test/test_qc_measures.py b/software/per-cell-metrics/test/test_qc_measures.py index 546e3eb..044280b 100644 --- a/software/per-cell-metrics/test/test_qc_measures.py +++ b/software/per-cell-metrics/test/test_qc_measures.py @@ -1,12 +1,15 @@ import dataclasses +import math import re +import numpy as np import polars as pl import pytest from qc_measures import ( _COMPARISON, DEFAULT_LINES, LINE_ROUTES, + LOG1P_BIN_WIDTH, MEASUREMENTS, Coverage, Line, @@ -15,14 +18,16 @@ Status, aggregate_barcode_fraction, antigen_count_deciles, - count_bin_edges, + bin_values, detect_aggregate_barcodes, + log1p_bin_edges, measurement_rows, per_antigen_measures, per_tag_count_bins, reads_per_cell, roll_up, sibling_disagreement, + status_expr, status_for, usable_read_fraction, ) @@ -594,14 +599,25 @@ def test_the_categorical_route_carries_cells_detected_and_nothing_else(): assert categorical_routed.isdisjoint(_COMPARISON) -def test_the_undeclared_barcode_line_is_read_direct_not_as_a_complement(): - # The line is published on the undeclared share itself: warn above 0.50, error at 1.0. The barcode - # table measures that share directly, so the thresholds are not mirrored. +def test_the_undeclared_barcode_line_is_per_barcode_and_operator_set(): + # Judged on ONE sequence's share of its sample's pre-refine reads, not on the sample's aggregate. + # The field's 0.50/1.0 is published for the aggregate and does not transfer: an aggregate reaches + # 0.50 while no single sequence comes near it, so that line would never fire on a row. line = DEFAULT_LINES["undeclaredBarcodeShare"] - assert (line.warn, line.error) == (0.50, 1.0) - assert status_for("undeclaredBarcodeShare", 0.60, DEFAULT_LINES) is Status.WARN + assert (line.warn, line.error) == (0.01, 0.05) + assert status_for("undeclaredBarcodeShare", 0.02, DEFAULT_LINES) is Status.WARN + assert status_for("undeclaredBarcodeShare", 0.10, DEFAULT_LINES) is Status.ALERT + assert status_for("undeclaredBarcodeShare", 0.005, DEFAULT_LINES) is Status.OK + + +def test_the_undeclared_barcode_line_alerts_above_its_error_not_only_at_it(): + # Both ends face the same way, unlike the four inherited shares. `alerting-at` on the error end + # would fire only at exactly 0.05 and let every larger share read warn -- the worse finding being + # the one that never showed. + assert _COMPARISON["undeclaredBarcodeShare"] == ("at-most", "at-most") + assert status_for("undeclaredBarcodeShare", 0.05, DEFAULT_LINES) is Status.WARN + assert status_for("undeclaredBarcodeShare", 0.0501, DEFAULT_LINES) is Status.ALERT assert status_for("undeclaredBarcodeShare", 1.0, DEFAULT_LINES) is Status.ALERT - assert status_for("undeclaredBarcodeShare", 0.40, DEFAULT_LINES) is Status.OK def test_panel_assigned_fraction_carries_no_line_any_more(): @@ -676,11 +692,37 @@ def test_a_stated_recommendation_warns_and_never_alerts(): def test_two_thresholds_give_three_levels(): # The distinction collapsing them lost. Three of the four inherited lines put error at total failure, # so a low-but-non-zero share warns and only a wholly failed one alerts. - line = DEFAULT_LINES["undeclaredBarcodeShare"] - assert (line.warn, line.error) == (0.5, 1.0) - assert status_for("undeclaredBarcodeShare", 0.4, DEFAULT_LINES) is Status.OK - assert status_for("undeclaredBarcodeShare", 0.6, DEFAULT_LINES) is Status.WARN - assert status_for("undeclaredBarcodeShare", 1.0, DEFAULT_LINES) is Status.ALERT + line = DEFAULT_LINES["aggregateBarcodeFraction"] + assert (line.warn, line.error) == (0.05, 1.0) + assert status_for("aggregateBarcodeFraction", 0.04, DEFAULT_LINES) is Status.OK + assert status_for("aggregateBarcodeFraction", 0.06, DEFAULT_LINES) is Status.WARN + assert status_for("aggregateBarcodeFraction", 1.0, DEFAULT_LINES) is Status.ALERT + + +def test_status_expr_agrees_with_status_for(): + # Two evaluators over one set of thresholds. The column form exists because the undeclared-barcode + # table's row cap accepts None, so a scalar loop there is a loop over every distinct pre-refine + # sequence. Splitting the evaluator is only safe while the two cannot disagree, so every registered + # measurement is run against both -- over its own thresholds, either side of each, and over the four + # values that are not numbers. + not_numbers = [None, float("nan"), float("inf"), float("-inf")] + ids = sorted(set(DEFAULT_LINES) | {m.id for m in MEASUREMENTS}) + for measurement in ids: + probes = [*not_numbers, 0.0, 0.5, 1.0, 5000.0] + line = DEFAULT_LINES.get(measurement) + if line is not None: + for threshold in (line.warn, line.error): + if threshold is not None: + probes += [threshold, threshold - 1e-9, threshold + 1e-9, threshold - 0.01, threshold + 0.01] + frame = pl.DataFrame({"v": probes}, schema={"v": pl.Float64}) + # `with_columns`, not `select`: a measurement with no line yields a scalar literal, and `select` + # would return one row of it rather than one per probe. + got = frame.with_columns(status_expr(measurement, pl.col("v"), DEFAULT_LINES).alias("s"))["s"].to_list() + want = [ + None if (status := status_for(measurement, value, DEFAULT_LINES)) is None else status.value + for value in probes + ] + assert got == want, measurement def test_error_is_tested_before_warn(monkeypatch): @@ -704,8 +746,8 @@ def test_at_least_is_acceptable_exactly_at_the_line(): def test_at_most_is_acceptable_exactly_at_the_line(): # `undeclaredBarcodeShare` reads `at-most`: the warn line itself satisfies the condition it names, and # only strictly above it warns. - assert status_for("undeclaredBarcodeShare", 0.5, DEFAULT_LINES) is Status.OK - assert status_for("undeclaredBarcodeShare", 0.51, DEFAULT_LINES) is Status.WARN + assert status_for("undeclaredBarcodeShare", 0.01, DEFAULT_LINES) is Status.OK + assert status_for("undeclaredBarcodeShare", 0.011, DEFAULT_LINES) is Status.WARN def test_the_undeclared_barcode_fraction_ships_unjudged(): @@ -1054,28 +1096,21 @@ def _bin_counts(rows: list[tuple[str, str, str, int]]) -> pl.DataFrame: ) -def test_one_edge_set_spans_the_whole_run(): - # Per-tag edges would rescale every panel to its own range, so a tag spanning 1-4 and one spanning - # 1-4000 would draw alike in a grid a reader compares side by side. - counts = _bin_counts([("S1", "c1", "AAAA", 2), ("S1", "c2", "BBBB", 4000)]) - edges = count_bin_edges(counts) - assert edges[0] == 1.0 - assert edges[-1] == 4000.0 - # Log-spaced, so the ambient population does not collapse into one bar. - assert edges[2] - edges[1] > edges[1] - edges[0] +def _log1p_edges(counts) -> list[float]: + """The edge set the count distributions are drawn on, for a bin frame.""" + top = int(counts["umiCount"].max() or 0) if counts.height else 0 + return log1p_bin_edges(top) -def test_a_frame_with_no_counts_has_no_edges_and_no_bins(): - empty = _bin_counts([]) - assert count_bin_edges(empty) == [] - assert per_tag_count_bins(empty, []) == {} +def test_a_frame_with_no_counts_has_no_bins(): + assert per_tag_count_bins(_bin_counts([]), []) == {} def test_every_cell_lands_in_a_bin_including_the_largest_count(): # np.histogram closes the last bin on the right. Without that the run's maximum count falls outside # every bucket, and the tag holding it reads one cell short. counts = _bin_counts([("S1", f"c{i}", "AAAA", n) for i, n in enumerate([1, 1, 2, 5, 40, 4000])]) - edges = count_bin_edges(counts) + edges = _log1p_edges(counts) weights = per_tag_count_bins(counts, edges)["S1"]["AAAA"] assert len(weights) == len(edges) - 1 assert sum(weights) == 6 @@ -1085,8 +1120,7 @@ def test_bins_are_kept_per_sample_and_tag(): # The fit runs per (sample, tag), so the plots drawn beside it are keyed the same way. Pooling two # samples would read as one population. counts = _bin_counts([("S1", "c1", "AAAA", 2), ("S2", "c1", "AAAA", 2), ("S1", "c1", "BBBB", 3)]) - edges = count_bin_edges(counts) - out = per_tag_count_bins(counts, edges) + out = per_tag_count_bins(counts, _log1p_edges(counts)) assert sorted(out) == ["S1", "S2"] assert sorted(out["S1"]) == ["AAAA", "BBBB"] assert sorted(out["S2"]) == ["AAAA"] @@ -1096,7 +1130,7 @@ def test_a_tag_absent_from_a_sample_gets_no_entry(): # An absent tag and a tag whose cells all read low are different findings. A list of zeros would state # the second, so the absent one carries no list at all. counts = _bin_counts([("S1", "c1", "AAAA", 2), ("S2", "c1", "BBBB", 2)]) - out = per_tag_count_bins(counts, count_bin_edges(counts)) + out = per_tag_count_bins(counts, _log1p_edges(counts)) assert "BBBB" not in out["S1"] assert "AAAA" not in out["S2"] @@ -1105,5 +1139,47 @@ def test_the_reference_tag_keeps_its_bins(): # The reference tag is the run's own ambient floor, which is what every other tag is judged against. # It is held out of the verdict read, never out of this. counts = _bin_counts([("S1", "c1", "CTRL", 6), ("S1", "c1", "AAAA", 500)]) - out = per_tag_count_bins(counts, count_bin_edges(counts)) + out = per_tag_count_bins(counts, _log1p_edges(counts)) assert sorted(out["S1"]) == ["AAAA", "CTRL"] + + +# --- the log1p bin edges the count distributions draw on ---------------------------------------- + + +def test_log1p_bin_edges_all_draw_the_same_width(): + """The whole point of these edges: every bar draws the same width. + + The plot's axis is log1p, so an edge at `expm1(k * width)` lands at `k * width` on screen. Equal + widths mean a bar's height is its share directly, which removes the per-bar division that + whole-number edges needed -- and that was once got wrong, hiding a real hump. + """ + for top in (30, 1381, 152757): + edges = log1p_bin_edges(top) + drawn = [math.log1p(edges[i + 1]) - math.log1p(edges[i]) for i in range(len(edges) - 1)] + assert drawn, f"no bins for top={top}" + for width in drawn: + assert width == pytest.approx(LOG1P_BIN_WIDTH, abs=1e-12), f"top={top}" + + +def test_log1p_bin_edges_keep_a_bar_for_the_zeros(): + # The fit runs over one entry per cell, and a cell that read nothing enters as a zero. Those zeros + # are most of the background: drop them and the plot shows one decaying hump whatever the fit found. + edges = log1p_bin_edges(1381) + assert edges[0] == 0.0 + assert bin_values(np.array([0, 0, 0, 4]), edges)[0] == 3 + + +def test_log1p_bin_edges_cover_the_highest_count(): + # `np.histogram` closes the last bin on the right, so an edge sitting exactly ON the top count would + # give that bin one count more than its width. The last edge is the first step above the top. + for top in (1, 7, 30, 1381, 152757): + edges = log1p_bin_edges(top) + assert edges[-1] > top, f"top={top} has no bin to land in" + assert sum(bin_values(np.array([top]), edges)) == 1 + + +def test_log1p_bin_edges_are_empty_where_there_is_nothing_to_span(): + assert log1p_bin_edges(0) == [] + assert log1p_bin_edges(-1) == [] + # A width of zero or less would divide by zero or loop forever, so it is refused rather than fixed up. + assert log1p_bin_edges(100, 0.0) == [] diff --git a/software/per-cell-metrics/test/test_tag_distribution.py b/software/per-cell-metrics/test/test_tag_distribution.py index 734b8fb..978e67e 100644 --- a/software/per-cell-metrics/test/test_tag_distribution.py +++ b/software/per-cell-metrics/test/test_tag_distribution.py @@ -22,8 +22,12 @@ from panel import ANY_SAMPLE from tag_distribution import ( DEFAULT_DISTRIBUTION_MIN_CELLS, + DEFAULT_INITIAL_SIGNAL_WEIGHT, NO_FIT, TOO_FEW_CELLS, + _fit_two_component_nb, + _signal_component, + bound_at_count, fit_tag_probabilities, fit_tag_probabilities_by_pair, ) @@ -129,17 +133,45 @@ def test_an_even_split_is_handled(): assert _bound(fit.probabilities)[1000:].all() -def test_a_cell_that_read_nothing_is_not_called_bound_beside_an_overdispersed_population(): - """The signal component is the higher-MEDIAN one, and the median is not ordered by the mean. +def test_the_signal_component_is_the_higher_median_one_even_where_the_means_disagree(): + """The labelling rule itself, on parameters rather than on a fitted bed. - An ambient population -- mostly zero with a few enormous counts -- fits a component whose mean sits - far above a real binder population's while its median sits far below. The Poisson-shaped binders - here are the higher-median component and the overdispersed ambient one is not, so ordering the two - by mean labels them the wrong way round. + `what-plays-the-baseline` fixes the rule as the higher-median component. A negative binomial's + median is not ordered by its mean -- the median depends on the size too -- so the two orderings can + disagree, and labelling by mean inverts every call for the tag: cells reading nothing score high and + cells reading a lot score low. That was a shipped bug, fixed by labelling on the median. - Asserted at a count of zero because that is where the inversion is unmistakable: under the mean - ordering every cell that read nothing lands in the ambient component with probability 1 and is - called bound. + Asserted on component parameters directly, NOT through a fit. The old form of this test built a bed + whose components the mean ordering got wrong, and so it also depended on the EM reaching one + particular decomposition of that bed -- which the initialisation decides. That made a labelling-rule + guard fail whenever the initialisation changed, for reasons that had nothing to do with labelling. + The rule is a pure function of the fitted parameters, so it is tested as one. + """ + # The pair from `_signal_component`'s own docstring. Component 0 has the higher MEAN, component 1 + # the higher MEDIAN: at mean 50 with size 0.05 the median is 0, at mean 5 with size 1e6 it is 5. + assert _signal_component(np.array([50.0, 5.0]), np.array([0.05, 1e6])) == 1 + + # Where the two orderings agree there is nothing to choose between them. + assert _signal_component(np.array([2.0, 40.0]), np.array([3.0, 4.0])) == 1 + + +def test_two_components_with_equal_medians_are_separated_by_their_means(): + # Fitted over mostly-zero counts, both components have a median of zero, and the published rule + # says nothing about that case -- the mean is all that is left to separate them. Pinned because + # that tie-break is our own choice, so it needs to be written down somewhere. + assert _signal_component(np.array([2.0, 8.0]), np.array([0.02, 0.02])) == 1 + + +def test_a_cell_that_read_nothing_is_never_called_bound(): + """A cell holding no count of a tag cannot bind it, whatever the fit made of the rest. + + This is the invariant `silent_tally` rests on. It resolves the unobserved positions by arithmetic -- + asked minus observed minus unreliable -- instead of reading each one, and that is only sound while a + zero-count cell cannot reach the bound line. A breach would make the run call a cell bound in one + place and count it not-bound in the other, with nothing raised. + + The bed is an ambient population that is mostly zero with a few enormous counts, which is the shape + that puts the most pressure on it. """ rng = np.random.default_rng(3) ambient = rng.negative_binomial(0.15, 0.15 / (0.15 + 40), 1000) @@ -151,8 +183,46 @@ def test_a_cell_that_read_nothing_is_not_called_bound_beside_an_overdispersed_po silent = fit.probabilities[counts == 0] assert silent.size > 0, "the bed must hold cells that read nothing for this to say anything" assert silent.max() < DISTRIBUTION_BOUND_PROBABILITY, "a cell holding no count of the tag cannot bind it" - # The direction holds over the two populations and not only at zero. - assert fit.probabilities[1000:].mean() > fit.probabilities[:1000].mean() + + +def test_an_ambient_tail_heavier_than_the_binders_is_labelled_the_signal_and_says_nothing(): + """KNOWN LIMITATION. This test records it; it does not guard against it. + + When the background counts have a long enough tail that their MEAN sits above the binders', the fit + splits the data the wrong way. Instead of {background} and {binders} it finds {everything} and {the + background's tail}, and calls that tail the signal. Every real binder then scores below the line, + and part of the background tail scores above it. + + Not one unlucky dataset: 0 of 12 seeds come out right on this shape, at every binder share from 2% + to 50%. It is a property of where the fit starts. + + The median start this used to have got this shape right -- and got the mostly-background case wrong + instead, which is every tag in a real run. Neither start wins both. The source paper reports the + same weakness for its own no-control path. + + THE RUN GIVES NO WARNING. That is what the asserts below are for: the fit comes back with its + signal mean above its background mean, so the probability rises with the count and an ordinary + threshold gets drawn. Nothing in the output distinguishes this from a fit that got it right. Anyone + building a check for it should start here. + """ + rng = np.random.default_rng(3) + ambient = rng.negative_binomial(0.15, 0.15 / (0.15 + 40), 1000) + binders = rng.poisson(12, 1000) + counts = np.concatenate([ambient, binders]) + + fit = fit_tag_probabilities(counts) + assert fit.reason is None + # The wrong population carries the signal label: the binders score BELOW the ambient. + assert fit.probabilities[1000:].mean() < fit.probabilities[:1000].mean() + assert not _bound(fit.probabilities[1000:]).any(), "no real binder clears the line on this bed" + assert _bound(fit.probabilities[:1000]).any(), "part of the ambient tail does" + # And it looks healthy. Two means the right way round, so a gate is drawn like any other fit. + assert fit.background.signal_mean > fit.background.mean + counts_seen = np.unique(counts) + curve = (counts_seen, np.array([fit.probabilities[counts == c][0] for c in counts_seen])) + assert bound_at_count(curve, DISTRIBUTION_BOUND_PROBABILITY) is not None, ( + "the fit resolves a bound count, so nothing marks this as suspect" + ) def test_the_probability_is_a_probability(): @@ -182,19 +252,32 @@ def test_a_tag_no_cell_read_at_all_does_not_fit(): assert fit.reason == NO_FIT -def test_a_tag_nothing_bound_still_fits_and_calls_some_cells_bound(): - # DELIBERATE: the method assumes two components exist, so it splits a single population and calls its - # upper slice signal. Rejecting this would be a separation test of our own invention. - # - # The background here is OVERDISPERSED rather than Poisson, which is what a real one is: the invented - # binders are the long tail of a single skewed population, so a Poisson bed would let this pass for - # the wrong reason. - # - # Pinned so that nobody restores the rejection as a bug fix. +def test_a_tag_nothing_bound_still_fits_and_now_calls_no_cell_bound(): + """A single population still FITS -- no rejection -- and no longer reaches the bound line. + + The method assumes two components exist, so it splits a single population and calls its upper slice + signal. That much is unchanged, and the pin below is unchanged with it: the fit must return + probabilities rather than refuse, because rejecting here would be a separation test of our own + invention and `what-plays-the-baseline` declines to build one. + + What changed is WHERE the split lands. This rung now starts the EM from the source paper's own + pivot rather than from the median, and from that start a single background population no longer + produces a slice above 0.9. The earlier assertion -- that some cell IS called bound -- recorded the + median start's behaviour, not a decision the spec had taken: the spec fixes the trim, the labelling + rule and the 0.9, and says nothing about the initialisation. + + So this test still pins the thing it was written to pin, that nobody restores the rejection, and + reads the count off the fit instead of requiring it to be non-zero. + + The background here is OVERDISPERSED rather than Poisson, which is what a real one is: a Poisson bed + would pass for the wrong reason. + """ rng = np.random.default_rng(SEED) fit = fit_tag_probabilities(rng.negative_binomial(3, 3 / (3 + 2), size=2000)) - assert fit.reason is None - assert _bound(fit.probabilities).any(), "the spec accepts invented binders on a tag that bound nothing" + assert fit.reason is None, "a single population must still fit rather than be rejected" + assert fit.probabilities is not None + assert fit.probabilities.size == 2000 + assert not _bound(fit.probabilities).any(), "the paper's pivot invents no binders on this bed" # --- the trim ----------------------------------------------------------------------------------- @@ -435,3 +518,194 @@ def test_padding_the_universe_with_ambient_barcodes_collapses_the_fit(): # apart. Both are unusable. over_barcodes = swamped.backgrounds.get(("S1", "AAAA")) assert over_barcodes is None or over_barcodes.signal_mean < over_barcodes.mean * 2 + + +def test_bins_count_every_cell_the_fit_saw_including_the_zeros() -> None: + # The plot a reader judges the background against is drawn from these bins, so they have to cover the + # cells the fit actually saw: one entry per cell in the sample, with a cell that read nothing counted + # as a zero. The sparse counts frame has no rows for those zeros, and on real data the zeros are most + # of the background -- bin the sparse frame instead and the plot shows one decaying hump no matter + # what the fit found, because the left half is simply missing. + counts, cells = _bed() + fits = fit_tag_probabilities_by_pair(counts, cells, _panel([("AAAA", "S1")])) + + weights = fits.bins[("S1", "AAAA")] + in_sample = len([c for c in cells if c[0] == "S1"]) + assert sum(weights) == in_sample, "every cell of the sample lands in a bin" + + observed = counts.filter((pl.col("sampleId") == "S1") & (pl.col("tag") == "AAAA")).height + assert observed < in_sample, "the bed must leave some cells unobserved or this pins nothing" + # The zeros sit in the first bin, which spans [0, expm1(width)) and so holds no other count. + assert weights[0] == in_sample - observed + + +def test_bins_are_recorded_even_where_no_fit_came_out() -> None: + # A pair where nothing separated is the case a reader most needs to see, so its bins are kept too. A + # tag the reads never showed is all zeros, and cannot separate. + _, cells = _bed() + empty = _counts_frame([]) + fits = fit_tag_probabilities_by_pair(empty, cells, _panel([("AAAA", "S1")])) + + assert ("S1", "AAAA") in fits.reasons + assert ("S1", "AAAA") not in fits.backgrounds + weights = fits.bins[("S1", "AAAA")] + # All zeros, so one bin holding every cell of the sample -- not an empty list. + assert weights == [len([c for c in cells if c[0] == "S1"])] + + +def test_bound_at_count_takes_the_start_of_the_final_crossing_region(): + counts = np.array([0, 1, 2, 5, 9]) + probabilities = np.array([0.01, 0.10, 0.55, 0.93, 0.99]) + assert bound_at_count((counts, probabilities), 0.9) == 5 + + +def test_bound_at_count_refuses_an_inverted_fit_rather_than_marking_its_low_crossing(): + # Components are labelled by median, so a pair can come back with its signal component BELOW its + # background, and then the probability falls as the count rises. The low counts cross; no count + # sits above the line, so there is no gate to draw. + counts = np.array([0, 1, 2, 5, 9]) + probabilities = np.array([0.99, 0.95, 0.40, 0.05, 0.01]) + assert bound_at_count((counts, probabilities), 0.9) is None + + +def test_bound_at_count_is_none_where_the_fit_never_reaches_the_cutoff(): + assert bound_at_count((np.array([0, 1, 4]), np.array([0.1, 0.2, 0.7])), 0.9) is None + + +def test_bound_at_count_rises_with_the_cutoff(): + curve = (np.array([0, 1, 2, 5, 9]), np.array([0.01, 0.10, 0.92, 0.97, 0.995])) + assert bound_at_count(curve, 0.9) == 2 + assert bound_at_count(curve, 0.95) == 5 + assert bound_at_count(curve, 0.999) is None + + +def test_curves_carry_one_entry_per_distinct_scored_count(): + counts = pl.DataFrame( + { + "sampleId": ["s"] * 6, + "cellId": [f"c{i}" for i in range(6)], + "tag": ["T"] * 6, + "umiCount": [1, 1, 2, 2, 40, 41], + } + ) + cells = [("s", f"c{i}") for i in range(400)] + panel = pl.DataFrame({"tag": ["T"], "sample": [ANY_SAMPLE]}) + fits = fit_tag_probabilities_by_pair(counts, cells, panel) + curve = fits.curves.get(("s", "T")) + if curve is None: + pytest.skip("this population produced no fit, so it carries no curve") + distinct, probabilities = curve + assert distinct.tolist() == [0, 1, 2, 40, 41] + assert probabilities.size == distinct.size + # Ascending counts, so a monotone fit gives a resolvable crossing. + assert np.all(np.diff(distinct) > 0) + + +# --- the expected binder fraction --------------------------------------------------------------- + + +def _from_bins(weights: list[int], top: int = 1381, bins: int = 24) -> np.ndarray: + """Per-cell counts rebuilt from a run's own histogram, at each bin's geometric middle. + + Coarse on purpose: it carries the SHAPE of a real distribution into a test without a fixture file. + """ + ratio = (top + 1) ** (1.0 / bins) + edges = [1] + while len(edges) < bins: + step = max(edges[-1] + 1, round(edges[-1] * ratio)) + if step >= top: + break + edges.append(step) + edges.append(top + 1) + edges.insert(0, 0) + out: list[int] = [] + for i, count in enumerate(weights): + if count == 0: + continue + low, high = edges[i], edges[i + 1] + out += [0 if low == 0 else int(round(np.sqrt(low * max(high - 1, low))))] * count + return np.array(out) + + +# One real tag from a synthetic bound panel: a tight background, a gap at counts 7-11, then a broad +# upper mode holding 27% of the cells. +_BINDER_RICH = [ + 2140, + 441, + 453, + 363, + 79, + 28, + 0, + 0, + 11, + 20, + 41, + 73, + 80, + 101, + 106, + 61, + 78, + 65, + 76, + 111, + 104, + 132, + 163, + 102, +] + + +def test_the_published_default_is_the_papers_initial_weight(): + assert DEFAULT_INITIAL_SIGNAL_WEIGHT == 0.1 + + +def test_a_binder_rich_tag_needs_a_higher_expected_fraction_to_split_at_its_own_gap(): + """Why the fraction is settable at all. + + The default assumes a tenth of cells bind. This tag's upper mode holds 27% of them, so the pivot + lands INSIDE that mode, the background component is seeded with most of the binders and ends up + wide enough to explain the whole range, and no count is ever 90% likely to be signal. Told the + right fraction, the same fit splits at the gap the histogram shows. + """ + counts = _from_bins(_BINDER_RICH) + x = counts.astype(float) + + default_fit = _fit_two_component_nb(x, DEFAULT_INITIAL_SIGNAL_WEIGHT) + informed_fit = _fit_two_component_nb(x, 0.3) + assert default_fit is not None and informed_fit is not None + + # The background the default settles on is two orders of magnitude wider than the real one. + assert min(default_fit.means) > 10.0 + assert min(informed_fit.means) < 2.0 + + +def test_the_expected_fraction_reaches_the_per_pair_driver(): + """Threaded rather than read from the module at the bottom of the stack. + + Asserted on the binder-rich shape, because that is where the weight changes the answer. A bed whose + two populations stand well apart converges to the same fit from either start -- correctly, since the + start only picks which optimum the EM walks to -- so it cannot tell a threaded weight from a + dropped one. + """ + population = _from_bins(_BINDER_RICH) + cells = [("s", f"c{i}") for i in range(population.size)] + observed = [(key[1], int(v)) for key, v in zip(cells, population) if v > 0] + counts = pl.DataFrame( + { + "sampleId": ["s"] * len(observed), + "cellId": [c for c, _ in observed], + "tag": ["T"] * len(observed), + "umiCount": [v for _, v in observed], + } + ) + panel = pl.DataFrame({"tag": ["T"], "sample": [ANY_SAMPLE]}) + + default = fit_tag_probabilities_by_pair(counts, cells, panel, initial_signal_weight=DEFAULT_INITIAL_SIGNAL_WEIGHT) + informed = fit_tag_probabilities_by_pair(counts, cells, panel, initial_signal_weight=0.3) + low, high = default.backgrounds.get(("s", "T")), informed.backgrounds.get(("s", "T")) + assert low is not None and high is not None + # The same divergence the direct-fit test pins, seen through the driver. + assert low.mean > 10.0 + assert high.mean < 2.0 diff --git a/test/src/qcDefaults.test.ts b/test/src/qcDefaults.test.ts index b44dc3f..b514a96 100644 --- a/test/src/qcDefaults.test.ts +++ b/test/src/qcDefaults.test.ts @@ -85,6 +85,12 @@ describe("VERDICT_DEFAULTS matches verdict-args.lib.tengo and the Python that ow const pairs: [keyof typeof VERDICT_DEFAULTS, string, string, string][] = [ ["countFloor", "DEFAULT_COUNT_FLOOR", "verdict.py", "DEFAULT_FLOOR"], ["boundCutoff", "DEFAULT_BOUND_CUTOFF", "verdict.py", "BOUND_CUTOFF"], + [ + "boundProbability", + "DEFAULT_BOUND_PROBABILITY", + "verdict.py", + "DISTRIBUTION_BOUND_PROBABILITY", + ], ["minVotingCells", "DEFAULT_MIN_VOTING_CELLS", "combine.py", "DEFAULT_MIN_VOTERS"], [ "panelReferenceMinMembers", @@ -108,7 +114,7 @@ describe("VERDICT_DEFAULTS matches verdict-args.lib.tengo and the Python that ow it("covers every shaping default the tengo file declares", () => { const declared = [ ...tengo.matchAll( - /^(DEFAULT_(?:COUNT_FLOOR|BOUND_CUTOFF|MIN_VOTING_CELLS|PANEL_MIN_MEMBERS|DISTRIBUTION_MIN_CELLS))\s*:=/gm, + /^(DEFAULT_(?:COUNT_FLOOR|BOUND_CUTOFF|BOUND_PROBABILITY|MIN_VOTING_CELLS|PANEL_MIN_MEMBERS|DISTRIBUTION_MIN_CELLS))\s*:=/gm, ), ].map((m) => m[1]); expect(new Set(declared)).toStrictEqual(new Set(pairs.map(([, name]) => name))); diff --git a/ui/src/components/CountHistogram.vue b/ui/src/components/CountHistogram.vue index c96e8a3..1ef3c9d 100644 --- a/ui/src/components/CountHistogram.vue +++ b/ui/src/components/CountHistogram.vue @@ -17,9 +17,11 @@ import { valuesFromBins } from "./binValues"; // 674px default overflows any narrower container and paints over whatever sits beside it. Measured here // with a ResizeObserver. const props = defineProps<{ - /** Bin boundaries, `weights.length + 1` of them, shared across every plot of a run. */ + /** Bin boundaries, shared across every plot of a run. */ edges: number[]; - /** Cells per bin, in edge order. */ + /** + * Cells per bin, in edge order. + */ weights: number[]; /** * The x axis. `log` suits counts per cell, which span orders of magnitude: on a linear axis the ambient @@ -48,12 +50,24 @@ const MIN_WIDTH = 220; const width = ref(MIN_WIDTH); let observer: ResizeObserver | undefined; +// Ignore resizes too small to change what is drawn. +// +// Any width change rewrites `settings`, which redraws the chart. The uikit's `drawBins` adds a tooltip div +// to the body on every draw and never removes it, so each redraw leaks one div per panel -- and dragging a +// window fires a resize every frame. A few pixels of tolerance is invisible to the eye and turns a drag +// from hundreds of redraws into a handful. +// +// Fixing the leak itself belongs in the uikit. This only stops us multiplying it. +const WIDTH_EPSILON = 4; + onMounted(() => { const el = host.value; if (el === undefined) return; observer = new ResizeObserver((entries) => { - const measured = entries[0]?.contentRect.width ?? 0; - width.value = Math.max(MIN_WIDTH, Math.floor(measured)); + const measured = Math.max(MIN_WIDTH, Math.floor(entries[0]?.contentRect.width ?? 0)); + // The value stored is the exact measurement, not a rounded one. A drag that ends within the + // threshold leaves the chart up to 4px off the container, which nobody can see. + if (Math.abs(measured - width.value) >= WIDTH_EPSILON) width.value = measured; }); observer.observe(el); }); @@ -95,10 +109,12 @@ const settings = computed(() => { ...common, type: "log-bins" as const, // A bin's own bounds travel with its weight, since this form bins nothing itself. - bins: props.weights.map((weight, i) => ({ + // Bar height is the plain cell count. Every caller's bars are the same width on screen, so the + // height already is the share and there is nothing to divide by. PADDED to the full edge set. + bins: Array.from({ length: Math.max(props.edges.length - 1, 0) }, (_, i) => ({ from: props.edges[i]!, to: props.edges[i + 1]!, - weight, + weight: props.weights[i] ?? 0, })), }; }); diff --git a/ui/src/components/FittedBackgroundGrid.vue b/ui/src/components/FittedBackgroundGrid.vue index c08bcd8..c4b4144 100644 --- a/ui/src/components/FittedBackgroundGrid.vue +++ b/ui/src/components/FittedBackgroundGrid.vue @@ -7,22 +7,41 @@ import CountHistogram from "./CountHistogram.vue"; // The fitted background, as a grid of small multiples: one panel per (sample, tag), which is the grain the // fit runs at. Aggregating to the tag would hide a reagent that separated in one sample and not in another. // -// Ordered by TAG first, so one reagent's samples sit side by side and the row is the comparison. A panel is -// titled `tag · sample` in the same order. +// Its only caller scopes it to one sample, so each panel is titled with just its reagent and the sample +// is named once above the grid. Left unscoped it draws every (sample, tag) pair instead and puts the +// sample back in each title. That mode still works; nothing uses it today. // // A grid rather than a selector, with any panel enlargeable on click: the judgement asked for reads at // thumbnail size, and behind a selector nobody looks at all of them. // -// No marker is drawn. The threshold slot means "the declared gate" on the reference-reading plot and "the -// bound cutoff" on the scores plot, so a third meaning here would make one marker say three things. +// The vertical marker is this pair's own bound count: the count from which its fit starts calling a cell +// bound. Other plots use the same slot for the same purpose -- a line to judge the distribution against. +// Each (sample, tag) pair is fitted on its own cells, so every panel needs its own line. // -// A panel carries its title, its plot, and the fit's own three numbers. No separated / does-not-separate -// label: no criterion for it exists. This panel is the substitute for the check nobody has built, so -// withholding the fit leaves the rung with no safeguard at all. +// A panel carries its title, its plot, its cell count, and the fit's own three numbers. No separated / +// does-not-separate label: no criterion for it exists. This panel is the substitute for the check nobody +// has built, so withholding the fit leaves the rung with no safeguard at all. +// +// Bar height is a plain cell count, and that is safe because every bar is the same width. The edges come +// from `log1p_bin_edges` and the axis is symlog, which for these values is log1p, so equal steps in the +// edges draw as equal widths on screen. +// +// That was not always true. With integer count edges the bars had different widths -- `[0, 1)` covered +// 0.301 of a decade and `[4, 5)` covered 0.079 -- so a raw count made a wide bar look taller than a +// narrow one holding the same density. Each bar had to be divided by its own width, and getting that +// division wrong once flattened a real signal hump out of sight. Equal widths remove the problem. const props = defineProps<{ bins: TagCountBins; /** Sample id -> the label a reader knows it by. A sample with no label renders as its own id. */ sampleLabels: Record; + /** + * Show only this sample, which is how the page uses it. + */ + onlySample?: string; + /** + * Order the barcodes by this list rather than by label. Anything missing from it follows, by label. + */ + tagOrder?: string[]; }>(); type Panel = { @@ -30,7 +49,12 @@ type Panel = { title: string; weights: number[]; /** The fit's own numbers, absent where nothing was fitted for this pair. */ - fit?: { backgroundMean: number; signalMean: number; backgroundWeight: number }; + fit?: { + backgroundMean: number; + signalMean: number; + backgroundWeight: number; + boundAtCount?: number | null; + }; }; // Sorted by tag then sample, so the grid's reading order is stable across runs: `Record` iteration order @@ -39,21 +63,36 @@ type Panel = { // Tags sort on the NAME a panel is titled with, not on the barcode behind it, so the grid reads in the // order it prints. A tag the panel named nowhere reads as its own barcode and sorts under it. const panels = computed(() => { - const samples = Object.keys(props.bins.bySample).sort(); + const all = Object.keys(props.bins.bySample).sort(); + const samples = props.onlySample === undefined ? all : all.filter((s) => s === props.onlySample); const tags = new Set(); for (const sample of samples) { for (const tag of Object.keys(props.bins.bySample[sample] ?? {})) tags.add(tag); } const tagName = (tag: string) => props.bins.tagLabels?.[tag] ?? tag; + // Declared order where the caller gave one, label order for whatever it does not mention. A declared + // list from one run and bins from another need not agree, so neither side is assumed to cover the other. + const declared = new Map((props.tagOrder ?? []).map((tag, i) => [tag, i])); + const ordered = [...tags].sort((a, b) => { + const ia = declared.get(a); + const ib = declared.get(b); + if (ia !== undefined && ib !== undefined) return ia - ib; + if (ia !== undefined) return -1; + if (ib !== undefined) return 1; + return tagName(a).localeCompare(tagName(b)); + }); const out: Panel[] = []; - for (const tag of [...tags].sort((a, b) => tagName(a).localeCompare(tagName(b)))) { + for (const tag of ordered) { for (const sample of samples) { const weights = props.bins.bySample[sample]?.[tag]; // A tag absent from this sample's panel has nothing to draw. if (weights === undefined) continue; out.push({ key: `${tag} ${sample}`, - title: `${tagName(tag)} · ${props.sampleLabels[sample] ?? sample}`, + title: + props.onlySample === undefined + ? `${tagName(tag)} · ${props.sampleLabels[sample] ?? sample}` + : tagName(tag), weights, fit: props.bins.fitsBySample?.[sample]?.[tag], }); @@ -63,8 +102,24 @@ const panels = computed(() => { }); // Counts span orders of magnitude across a panel, so a fixed number of decimals prints either noise or -// nothing. Three significant figures reads the same at 0.33 and at 930. -const fmt = (value: number) => Number(value.toPrecision(3)).toLocaleString(); +// nothing. Three SIGNIFICANT figures reads the same at 0.00049 and at 930. +// +// Three SIGNIFICANT digits, via `maximumSignificantDigits`. +const fmt = (value: number) => value.toLocaleString("en-US", { maximumSignificantDigits: 3 }); + +// The bound-count line, which has four possible readings. Written once here because two places show it: +// the thumbnail shows it alone, the enlarged panel shows it under the fit's numbers. +const boundLine = (fit: Panel["fit"]) => { + if (fit === undefined) return "no fit for this barcode"; + // Three different facts, so three different sentences. A number is the fit's answer. `null` means the + // fit ran and no count reached the line. `undefined` means the run predates this field and never + // looked -- reporting that as "no count reaches it" would state a finding no run produced. + if (typeof fit.boundAtCount === "number") { + return `bound from \u2265${fit.boundAtCount.toLocaleString("en-US")} UMI`; + } + if (fit.boundAtCount === null) return "no UMI count reaches the bound probability"; + return "this run recorded no bound threshold"; +}; const enlarged = ref(undefined); const isOpen = computed({ @@ -93,19 +148,43 @@ const isOpen = computed({ - - - bg {{ fmt(panel.fit.backgroundMean) }} · signal {{ fmt(panel.fit.signalMean) }} · - {{ (panel.fit.backgroundWeight * 100).toFixed(0) }}% of cells background - - no fit for this pair + + + {{ boundLine(panel.fit) }} @@ -176,4 +255,16 @@ const isOpen = computed({ color: var(--color-txt-03); font-variant-numeric: tabular-nums; } + +/* The enlarged panel's readout: one statement per line, because these are read rather than scanned and + three facts on one line ran together. Tabular figures so the two means line up under each other. */ +.enlargedFit { + display: flex; + flex-direction: column; + gap: 4px; + margin-top: 12px; + font-size: 13px; + color: var(--color-txt-02); + font-variant-numeric: tabular-nums; +} diff --git a/ui/src/components/PatternEditor.vue b/ui/src/components/PatternEditor.vue index 6d16e7e..d137184 100644 --- a/ui/src/components/PatternEditor.vue +++ b/ui/src/components/PatternEditor.vue @@ -139,8 +139,8 @@ watch(editorMode, (mode) => { @update:model-value="setPresetId" > @@ -180,22 +180,22 @@ watch(editorMode, (mode) => { @@ -209,8 +209,8 @@ watch(editorMode, (mode) => { > { label="Feature barcode length" > diff --git a/ui/src/components/PunchLegend.vue b/ui/src/components/PunchLegend.vue index ad2fda0..80c905b 100644 --- a/ui/src/components/PunchLegend.vue +++ b/ui/src/components/PunchLegend.vue @@ -30,7 +30,7 @@ const SET_ENTRIES: Entry[] = [ { glyph: "unreliable", label: "Unreliable", - meaning: "asked, and the readings could not settle it — hover for which of the five ways", + meaning: "asked, and the readings could not settle it. Hover for which of the five ways", }, { glyph: "none", @@ -56,7 +56,7 @@ const CELL_ENTRIES: Entry[] = [ { glyph: "unreliable", label: "Unreliable", - meaning: "this cell could not be compared at all, so it cast no vote — hover for why", + meaning: "this cell could not be compared at all, so it cast no vote. Hover for why", }, { glyph: "none", @@ -98,7 +98,8 @@ const listStyle: CSSProperties = {
{{ e.label }} — {{ e.meaning }}{{ e.label }}: {{ e.meaning }}
diff --git a/ui/src/pages/AntigenQcPage.vue b/ui/src/pages/AntigenQcPage.vue index 7ecff06..a17d1cd 100644 --- a/ui/src/pages/AntigenQcPage.vue +++ b/ui/src/pages/AntigenQcPage.vue @@ -8,7 +8,9 @@ import { PlAgDataTableV2, PlAlert, PlBlockPage, + PlDropdown, PlPlaceholder, + PlRow, PlTabs, usePlDataTableSettingsV2, useWatchFetch, @@ -224,6 +226,11 @@ const VIEW_TABS = computed(() => [ const activeView = ref("reagents"); +// VIEW_TABS is derived from the rung, and the rung is unreported until the run settles, so a strip drawn +// mid-run offers every plot and then drops the ones that rung cannot draw. It stays hidden until then. The +// open view's body keeps rendering and draws its own processing placeholder. +const isRunning = computed(() => app.model.outputs.isRunning === true); + // The open tab can stop existing: the run reports its rung, and the plot that tab held cannot be drawn. // Falling back keeps the page showing something rather than an empty body under a tab strip that no longer // offers the tab. Watching an output and writing a LOCAL ref is not a hairpin: nothing here reaches @@ -246,6 +253,29 @@ const referenceSpread = computed(() => tagBins.value?.spreads?.referenceReading) // because it asks one question and offers no axes to pick. const tagBins = computed(() => app.model.outputs.tagCountBins); +const backgroundSample = ref(undefined); + +const backgroundSampleOptions = computed(() => { + const labels = app.model.outputs.sampleLabels ?? {}; + return Object.keys(tagBins.value?.bySample ?? {}) + .map((id) => ({ value: id, label: labels[id] ?? id })) + .sort((a, b) => a.label.localeCompare(b.label)); +}); + +// Keep the current selection if it still exists, otherwise fall back to the first sample. A re-run can +// drop the sample that was on screen, and an empty selector next to a full grid looks broken. +// +// The selector shows even when there is only one sample, because the panel titles no longer name it. +watch( + backgroundSampleOptions, + (options) => { + if (!options.some((o) => o.value === backgroundSample.value)) { + backgroundSample.value = options[0]?.value; + } + }, + { immediate: true }, +); + // Status is rendered as the plain string the workflow emitted, with the discrete filter its spec declares. // It stays plain text because of the fourth case a tag cannot render: a measurement with no line behind it // leaves this column empty, and an empty cell beside three tags reads as a tag that failed to load. Which @@ -267,7 +297,7 @@ const tagBins = computed(() => app.model.outputs.tagCountBins); - + + + diff --git a/ui/src/pages/MainPage.vue b/ui/src/pages/MainPage.vue index 072aac3..59d46e0 100644 --- a/ui/src/pages/MainPage.vue +++ b/ui/src/pages/MainPage.vue @@ -1,7 +1,20 @@