diff --git a/.changeset/integer-count-bin-edges.md b/.changeset/integer-count-bin-edges.md new file mode 100644 index 0000000..239ac95 --- /dev/null +++ b/.changeset/integer-count-bin-edges.md @@ -0,0 +1,27 @@ +--- +'@platforma-open/milaboratories.feature-integration.per-cell-metrics': patch +'@platforma-open/milaboratories.feature-integration.ui': patch +--- + +The fitted background draws the distribution the run holds + +Two defects, one picture. The fitted-background grid is the only way to see whether a tag's counts +separated into two populations, so a hump it invents is the one error this surface cannot carry. + +**Bins now sit on whole numbers.** `np.geomspace` puts edges between whole numbers, and a UMI count is +a whole number, so a bin could fall strictly between two counts and stand empty at every weight a run +could produce. On a 15-tag run topping out at 5,155 counts the bin at [2.039, 2.911) held nothing on +all 15 panels and read as a missing bar. Edges are now whole and strictly increasing, so every bin +holds at least one count. The low end steps by 1. Above that the step is geometric as before. The last +edge is one past the top count, which makes every bin half-open. A run now takes at most 24 bins +instead of always 24. + +**Bars are now density, cells per count.** Bin width in counts rises across the edge set — 1, 1, 2, 3, +4, 6, 9, 14, 21 on one real run — so a bin spanning 4 counts stood about four times a neighbour +spanning 1 at equal density. That step read as a second hump on tags whose counts only decay. Dividing +each weight by the counts its bin spans removes it. The y axis is labelled "Cells per count", and each +panel's caption carries the cell count the plot is drawn from. + +The hover readout on that grid now reports the density rather than the cell count, because +`PlChartHistogram` prints the number it is handed under a fixed `count:` label. The caption carries the +magnitude instead. 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/model/src/index.ts b/model/src/index.ts index d30a3dc..157c915 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"; @@ -41,6 +41,12 @@ 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, + // tag_distribution.py DEFAULT_INITIAL_SIGNAL_WEIGHT, the source paper's own initial weight. A default + // only -- unlike the line above this is not a floor, because it states what the EXPERIMENT holds + // rather than what the method licenses, and a sorted or synthetic population can sit far from it. + expectedBinderFraction: 0.1, // 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 +70,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 +221,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 +235,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 +524,10 @@ type BlockDataV2 = Omit< | "distributionMinCells" | "countFloor" | "boundCutoff" + | "boundProbability" | "minVotingCells" | "minAgreement" + | "expectedBinderFraction" | "gateThreshold" | "grouping" | "contendingGroups" @@ -710,6 +731,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 +750,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 +890,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/software/per-cell-metrics/src/emit_verdicts.py b/software/per-cell-metrics/src/emit_verdicts.py index 568b63f..1eae731 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 ( @@ -122,13 +123,16 @@ ) from tag_distribution import ( DEFAULT_DISTRIBUTION_MIN_CELLS, + DEFAULT_INITIAL_SIGNAL_WEIGHT, 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 +219,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 +244,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=DEFAULT_INITIAL_SIGNAL_WEIGHT) 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 +294,20 @@ 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. + if 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 +497,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 +551,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 +568,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 +720,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 +1040,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 +1057,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 +1542,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 +1641,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..269b765 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,9 @@ 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. +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 +170,9 @@ class _Mixture(NamedTuple): signal: int -def _fit_two_component_nb(counts: np.ndarray) -> _Mixture | None: +def _fit_two_component_nb( + counts: np.ndarray, initial_signal_weight: float = DEFAULT_INITIAL_SIGNAL_WEIGHT +) -> _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 +192,13 @@ 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)) + # Split the counts so that `initial_signal_weight` of them sit above the pivot. Each side then + # seeds one component, and the two side sizes are the starting weights. + pivot = float(np.quantile(x, 1.0 - 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 +282,7 @@ def fit_tag_probabilities( counts: np.ndarray, min_cells: int = DEFAULT_DISTRIBUTION_MIN_CELLS, scored: np.ndarray | None = None, + initial_signal_weight: float = DEFAULT_INITIAL_SIGNAL_WEIGHT, ) -> TagFit: """One tag's counts across one sample's cells, as a probability of binding per cell. @@ -308,7 +322,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 +348,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 +368,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 = DEFAULT_INITIAL_SIGNAL_WEIGHT, ) -> TagFits: """One fit per (sample, tag) the panel declares, scored per cell. @@ -394,7 +414,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 +430,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 +458,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..4f64c19 100644 --- a/test/src/qcDefaults.test.ts +++ b/test/src/qcDefaults.test.ts @@ -85,6 +85,18 @@ 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", + ], + [ + "expectedBinderFraction", + "DEFAULT_EXPECTED_BINDER_FRACTION", + "tag_distribution.py", + "DEFAULT_INITIAL_SIGNAL_WEIGHT", + ], ["minVotingCells", "DEFAULT_MIN_VOTING_CELLS", "combine.py", "DEFAULT_MIN_VOTERS"], [ "panelReferenceMinMembers", @@ -108,7 +120,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|EXPECTED_BINDER_FRACTION|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..2a530ed 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..65b5b92 100644 --- a/ui/src/pages/MainPage.vue +++ b/ui/src/pages/MainPage.vue @@ -35,6 +35,7 @@ import { AGGREGATE_DETECTION_DEFAULTS, groupingColumns, QC_LINE_DEFAULTS, + VERDICT_DEFAULTS, } from "@platforma-open/milaboratories.feature-integration.model"; import type { ImportFileHandle } from "@platforma-sdk/model"; import { parseTagCsvMeta } from "../csvMeta"; @@ -227,6 +228,19 @@ const allSources = computed(() => referenceSources.value?.options ?? []); // // The data keeps the share, so this is a display conversion and not a migration. Rounded on the way in, // because a percentage entered as an integer must come back as the same integer. +// The expected binder fraction as a PERCENTAGE, where the data holds a share from 0 to 1. Same display +// conversion as the agreement limit below, and for the same reason: a scientist states "about 30% of these +// cells bound", not "0.3". Not rounded on the way out, so a typed value comes back as the number typed. +const binderPercent = computed({ + get: () => { + const share = app.model.data.expectedBinderFraction; + return typeof share === "number" ? share * 100 : undefined; + }, + set: (percent: number | undefined) => { + app.model.data.expectedBinderFraction = typeof percent === "number" ? percent / 100 : undefined; + }, +}); + const agreementPercent = computed({ get: () => { const share = app.model.data.minAgreement; @@ -290,7 +304,7 @@ const groupingOptions = computed(() => [ value: TAG_GROUPING_VALUE, // Labelled so it cannot be mistaken for one of the panel's own columns, which is what naming it // after the barcode column did. - label: "Each barcode on its own — one identity per barcode", + label: "One identity per barcode", }, ...panelPropertyOptions.value, ]); @@ -320,10 +334,10 @@ const combineColumnError = computed(() => { if (!c) return undefined; if (c === app.model.data.barcodeSeqColumn || c === app.model.data.featureNameColumn) return ( - `The Combine-mode column must be a column of its own — it holds each feature's mode ` + - `("sum" or "all"), not barcodes or feature names. It's currently set to "${c}", the same ` + - `column used for the ${c === app.model.data.barcodeSeqColumn ? "barcode sequence" : "feature name"}. ` + - `Pick a different column, or clear it to sum all co-barcodes.` + `The stored Combine-mode column is "${c}", which this panel already uses for the ` + + `${c === app.model.data.barcodeSeqColumn ? "barcode sequence" : "feature name"}. ` + + `A combine-mode column must be a column of its own. It holds each feature's mode, "sum" or "all". ` + + `Upload the panel file again to clear it. Each feature then sums the counts of its barcodes.` ); return undefined; }); @@ -836,8 +850,8 @@ const gridOptions = { @@ -875,8 +889,8 @@ const gridOptions = { > @@ -968,8 +982,56 @@ const gridOptions = { A cell reads bound where its score reaches this number, from 0 to 100. The score is how certain it is that the antigen makes up more than 92.5% of the antigen and baseline counts.

- Certainty, not strength — two counts against zero score low. Cell Ranger says - this score does not measure binding strength. + Certainty, not strength: two counts against zero score low. The score does not + measure binding strength. + + + + + + + + + @@ -1057,11 +1119,11 @@ const gridOptions = { > @@ -1083,8 +1145,8 @@ const gridOptions = { @@ -1114,7 +1176,7 @@ const gridOptions = { Reads matched per cell in the cell list. The measurement warns below this count. It has no alert line, because the vendor published one boundary.

Default 5000. The vendor recommends this minimum for this assay type. Nothing calibrates - it against your own data, and no test asserts it. + it against your own data. @@ -1130,10 +1192,10 @@ const gridOptions = { @@ -1164,10 +1225,10 @@ const gridOptions = { label="Undeclared reads warn" > @@ -1200,8 +1260,8 @@ const gridOptions = { diff --git a/workflow/src/column-specs.lib.tengo b/workflow/src/column-specs.lib.tengo index 7b64e35..6cd0dde 100644 --- a/workflow/src/column-specs.lib.tengo +++ b/workflow/src/column-specs.lib.tengo @@ -35,9 +35,9 @@ guardNoScore := func(annotations) { } /* A NEW column must never key on pl7.app/feature/featureId while carrying barcode sequences: the axis - identity would be unchanged and its value space inverted, so no query would fail and joins would return - wrong rows. Every antigen axis is minted under its own name. The block test asserts the same over - EMITTED columns, which needs a running workflow. */ +identity would be unchanged and its value space inverted, so no query would fail and joins would return +wrong rows. Every antigen axis is minted under its own name. The block test asserts the same over +EMITTED columns, which needs a running workflow. */ // Split from the assertion below so it can be tested both ways. Tengo has no try/catch, so a predicate is // the only part of a guard a test can exercise. usesLegacyFeatureAxis := func(axesSpec) { @@ -208,7 +208,7 @@ perCellSummaryOutput := func(blockId, sampleAxisName) { valueType: "String", annotations: a(84000, true, { "pl7.app/label": "Feature breakdown", - "pl7.app/description": "Every feature this cell has signal for, as 'feature (fraction%, UMI count)', sorted by descending fraction (dominant feature first).", + "pl7.app/description": "Every feature this cell has signal for, as 'feature (share%, N UMI)', comma-separated, largest share first. A share under 1% prints as '<1%'.", "pl7.app/isSummary": "true" }) } @@ -310,9 +310,9 @@ qcSummaryColumnsSpec := func(sampleAxisSpec) { numCol("featuresDetected", "featuresDetected", "Features detected", 86000, undefined), numCol("totalUniqueUmis", "totalUniqueUmis", "Total distinct UMIs", 85000, undefined), numCol("medianUmisPerCell", "medianUmisPerCell", "Median UMIs / cell", 84000, ".1f"), - numCol("panelAssignedFraction", "panelAssignedFraction", "Panel-assigned fraction", 83000, ".2p", "Fraction of feature-barcode reads kept after correcting their barcode against the panel; reads too far from any panel entry are dropped. Its complement is the share landing in barcodes the panel never declared, which is where this measurement's line comes from. A low value flags a panel or read-geometry mismatch."), - numCol("cellBarcodeValidFraction", "cellBarcodeValidFraction", "Valid cell-barcode fraction", 82000, ".2p", "Fraction of reads whose cell barcode corrects onto the whitelist this chemistry produces. A low value points at the wrong whitelist or the wrong read geometry rather than at the panel."), - numCol("aggregateBarcodeFraction", "aggregateBarcodeFraction", "Aggregate-barcode read fraction", 81000, ".2p", "Reads in barcodes flagged as aggregates by the top-100 IQR rule, over readsTotal."), + numCol("panelAssignedFraction", "panelAssignedFraction", "Panel-assigned fraction", 83000, ".2p"), + numCol("cellBarcodeValidFraction", "cellBarcodeValidFraction", "Valid cell-barcode fraction", 82000, ".2p"), + numCol("aggregateBarcodeFraction", "aggregateBarcodeFraction", "Aggregate-barcode read fraction", 81000, ".2p"), numCol("aggregateBarcodesFlagged", "aggregateBarcodesFlagged", "Aggregate barcodes flagged", 80000, undefined), numCol("aggregateBarcodeThreshold", "aggregateBarcodeThreshold", "Aggregate-barcode UMI threshold", 79000, ".1f") ], @@ -394,15 +394,17 @@ identityAxis := func(blockId, groupingId) { } // panelAxis: one distinct declared tag set. No panel file names its panel, so the id is a hash of the -// sorted tag list (emit_verdicts.py `_panel_id`), stable across re-runs of the same declaration. Where one -// panel covers every sample the axis takes a single value and still renders, holding that one value on -// every row: the SDK suppresses no axis for being constant. +// sorted tag list (emit_verdicts.py `_panel_id`), stable across re-runs of the same declaration. + panelAxis := func(blockId) { return { name: "pl7.app/antigen/panelId", type: "String", domain: { "pl7.app/blockId": blockId }, - annotations: { "pl7.app/label": "Panel" } + annotations: { + "pl7.app/label": "Panel", + "pl7.app/table/visibility": "hidden" + } } } @@ -551,7 +553,7 @@ verdictsImportSpec := func(setAxisSpec, identityAxisSpec, served) { "pl7.app/label": "Why unsettled", "pl7.app/isDiscreteFilter": "true", "pl7.app/discreteValues": UNRELIABLE_REASONS, - "pl7.app/description": "A statement left unsettled by a position the experiment never asked calls for a panel change; one left unsettled by a reading that did not survive calls for a re-run. A bare 'unreliable' cannot tell them apart." + "pl7.app/description": "Why this verdict is not settled. 'never-offered' pairs with 'never asked'. The other five explain 'unreliable': no comparator, all cells gated, tie, below the agreement limit, too few voters." }) }), // The two sets `four-state-verdict` names, and they are NOT the same set. The cells a question @@ -563,7 +565,7 @@ verdictsImportSpec := func(setAxisSpec, identityAxisSpec, served) { name: "pl7.app/antigen/cellsAsked", valueType: "Int", annotations: a(99000, true, { - "pl7.app/label": "Cells the question was put to", + "pl7.app/label": "Cells asked", "pl7.app/min": "0", "pl7.app/description": "How many of this clonotype's cells sat in a sample whose panel offered this antigen and whose reads carried at least one of its tags. Cells a gate set aside, or with no baseline to read against, are in this number and did not vote." }) @@ -597,10 +599,10 @@ verdictsImportSpec := func(setAxisSpec, identityAxisSpec, served) { name: "pl7.app/antigen/wasCompeted", valueType: "String", annotations: a(96500, true, { - "pl7.app/label": "Reading was competed", + "pl7.app/label": "Competed by a bound antigen", "pl7.app/isDiscreteFilter": "true", "pl7.app/discreteValues": BOOL_VALUES, - "pl7.app/description": "True where this antigen read 'not bound' and something it was declared to compete with read 'bound' for the same clonotype. A statement can test this; the state itself is unchanged." + "pl7.app/description": "True where this antigen read 'not bound' and something it was declared to compete with read 'bound' for the same clonotype. A statement can test this. The state itself is unchanged." }) }), verdictCol("competedWith", "competedWith", { @@ -654,16 +656,16 @@ setCountsImportSpec := func(setAxisSpec, served) { // print the same word, so the count is what tells them apart. It does not vary by identity, which is // why it belongs here rather than in the card. countCol("cellCount", "pl7.app/antigen/cellCount", "Cells", 96000, true, - "How many cells this clonotype has. The verdicts rest on these cells; how many of them could answer at a given identity travels with that identity's verdict."), + "How many cells this clonotype has. The verdicts rest on these cells. How many of them could answer at a given identity travels with that identity's verdict."), // Off by default. It qualifies the cell count directly above it -- forty cells of which thirty-eight // read nothing is a different clonotype from forty that all read something -- so it is ordered next to // it rather than among the identity counts below. countCol("cellsReadingNothing", "pl7.app/antigen/cellsReadingNothing", "Cells that read nothing", 95500, false, - "How many of this clonotype's cells were left with no count on any tag, the baseline included, once the minimum count had run. It changes no verdict — those cells vote not bound like any other — and it separates a negative resting on cells that read something from one resting on cells that read nothing."), + "How many of this clonotype's cells were left with no count on any tag, the baseline included, once the minimum count had run. It changes no verdict. Those cells vote not bound like any other. It separates a negative resting on cells that read something from one resting on cells that read nothing."), countCol("boundCount", "pl7.app/antigen/boundCount", "Identities bound", 95000, true, "How many distinct antigen identities this clonotype bound. Breadth, never strength: nothing here says how well it bound any of them."), countCol("offeredCount", "pl7.app/antigen/offeredCount", "Identities offered", 94000, true, - "How many identities this clonotype's cells were actually stained with — the denominator a rate must use, since a clone offered eight of ten and binding all eight failed nothing."), + "How many identities this clonotype's cells were offered: declared on the sample's panel and present in its reads. A rate must use this denominator. A clone offered eight of ten and binding all eight failed nothing."), countCol("settledCount", "pl7.app/antigen/settledCount", "Identities settled", 93500, true, "Of the identities offered, how many the data could settle either way."), countCol("unsettledCount", "pl7.app/antigen/unsettledCount", "Identities unsettled", 93000, true, @@ -753,11 +755,10 @@ identityPivotColumns := func(identities, groupingId, served, colName, valueType, if !is_undefined(conflicts[identity]) { // One barcode, several names, and the identity is unaffected. The KEY is the barcode, and the // readings under this column are that one barcode's. - note = { "pl7.app/description": "The panel gives this one barcode more than one name — " + + note = { "pl7.app/description": "The panel gives this one barcode more than one name: " + text.join(conflicts[identity], ", ") + - ". Different samples name it differently, so every name it was given is shown. This is still a " + - "single barcode and a single identity; only the naming disagrees. Check the panel file if the " + - "names should have matched." } + ". Different samples name it differently, so every name is shown. It is one barcode and one " + + "identity. Only the naming disagrees. Check the panel file if the names should match." } } cols = append(cols, { column: identity, @@ -809,7 +810,6 @@ identityPivotImportSpec := func(setAxisSpec, identities, groupingId, served, col } } - // --- The punchcard's pivot: result_identity_punch.csv, keyed (setId) -------------------- // // One column per identity whose value carries the state AND everything behind it, as @@ -832,7 +832,7 @@ identitySummaryImportSpec := func(setAxisSpec, identities, groupingId, served) { "pl7.app/antigen/identityVerdict", "String", 92000, "identity_", false, { "pl7.app/isDiscreteFilter": "true", "pl7.app/discreteValues": VERDICT_STATES - }) + }) } // --- Re-derivation material ------------------------------------------------------------- @@ -883,7 +883,7 @@ cellPunchImportSpec := func(sampleAxisSpec, cellAxisSpec, identities, groupingId // subset of the larger. annotations: a(95000, true, { "pl7.app/label": "Identities this cell bound", - "pl7.app/description": "How many identities this cell read as bound, over the identities its sample was stained for. A silent position counts as the not-bound it resolves to; a cell that could not be compared counts none.", + "pl7.app/description": "How many identities this cell read as bound, over the identities its sample was stained for. A silent position counts as the not-bound it resolves to. A cell that could not be compared counts none.", "pl7.app/min": "0" }) } @@ -923,7 +923,7 @@ cellTagCountsImportSpec := func(sampleAxisSpec, cellAxisSpec, tagAxisSpec) { valueType: "Int", annotations: a(80000, true, { "pl7.app/label": "Tag UMI count", - "pl7.app/description": "Molecules seen for this tag in this cell, after the count floor. Sparse: a cell with no row for a tag saw nothing for it, which is a reading, not a gap.", + "pl7.app/description": "Molecules seen for this tag in this cell, after the count floor. Sparse: a cell with no row for a tag saw nothing for it, which is a reading, not a gap. The baseline tag has no row here.", "pl7.app/min": "0", "pl7.app/isAbundance": "true", "pl7.app/abundance/unit": "molecules", @@ -936,6 +936,35 @@ cellTagCountsImportSpec := func(sampleAxisSpec, cellAxisSpec, tagAxisSpec) { } } +// result_cell_raw_counts.csv, keyed (sampleId, cellId, tag). +cellRawTagCountsImportSpec := func(sampleAxisSpec, cellAxisSpec, tagAxisSpec) { + return { + axes: [ + { column: "sampleId", spec: sampleAxisSpec }, + { column: "cellId", spec: cellAxisSpec }, + { column: "tag", spec: tagAxisSpec } + ], + columns: [{ + column: "umiCount", + id: "tagRawUmiCount", + spec: { + name: "pl7.app/antigen/rawUmiCount", + valueType: "Int", + annotations: a(80000, true, { + "pl7.app/label": "Tag UMI count (pre-floor)", + "pl7.app/description": "Molecules seen for this tag in this cell, BEFORE the count floor: no count has been zeroed and the comparator's exemption from the floor does not apply, so a cell's tags sum to what that cell held.", + "pl7.app/min": "0", + "pl7.app/isAbundance": "true", + "pl7.app/abundance/unit": "molecules", + "pl7.app/abundance/normalized": "false" + }) + } + }], + storageFormat: "Parquet", + partitionKeyLength: 1 + } +} + // result_cell_scalars.csv, keyed (sampleId, cellId). // // The comparator is in the domain for the same reason it is on the verdicts: referenceCount and @@ -962,22 +991,22 @@ cellScalarsImportSpec := func(sampleAxisSpec, cellAxisSpec, served) { scalarCol("referenceCount", "pl7.app/antigen/referenceCount", "Int", a(79000, true, { "pl7.app/label": "Baseline reading", - "pl7.app/description": "What this cell's reading was compared against. Empty where the run had no baseline for the cell.", + "pl7.app/description": "The declared baseline tag's count in this cell. Under a declared baseline every reading is compared against it. Under a fitted distribution it feeds only the admissibility gate. Empty where the cell has none.", "pl7.app/min": "0" - })), + })), scalarCol("admissibility", "pl7.app/antigen/admissibility", "String", a(78000, true, { "pl7.app/label": "Admissibility", "pl7.app/isDiscreteFilter": "true", "pl7.app/discreteValues": ADMISSIBILITY_VALUES, "pl7.app/description": "Whether this cell's readings could be compared at all, and if not, why. An inadmissible cell casts no vote; it does not vote 'not bound'." - })), + })), scalarCol("inCellList", "pl7.app/antigen/inCellList", "String", a(77000, false, { "pl7.app/label": "In cell list", "pl7.app/isDiscreteFilter": "true", "pl7.app/discreteValues": BOOL_VALUES - })) + })) ], storageFormat: "Parquet", partitionKeyLength: 0 @@ -1157,7 +1186,7 @@ panelLabelsImportSpec := func(panelAxisSpec) { spec: { name: "pl7.app/label", valueType: "String", - annotations: a(0, true, { "pl7.app/label": "Panel" }) + annotations: a(0, false, { "pl7.app/label": "Panel" }) } }], storageFormat: "Parquet", @@ -1240,8 +1269,9 @@ panelMismatchImportSpec := func(panelAxisSpec, tagAxisSpec) { // // Barcodes the reads carried that no panel declares get their own table, keyed by sequence, because they // have no row above -- they are not in the panel. It is the one thing on this surface that carries a -// status: the share of a sample's reads landing in undeclared barcodes. That status stays the barcode's -// and never the sample's, so it sits here rather than on any row of qcImportSpec's per-sample list. +// status, and that status reads the ROW's own share of its sample's pre-refine reads. It stays the +// barcode's and never the sample's, so it sits here rather than on any row of qcImportSpec's per-sample +// list. The sample-level share keeps its own column beside it and carries no status. // // Usually this table has no row for a sample at all, and that is the wanted outcome. The description // below says so, since an empty grid otherwise reads as a check that never ran. @@ -1261,7 +1291,7 @@ undeclaredBarcodeImportSpec := func(sampleAxisSpec, tagAxisSpec) { annotations: a(79000, true, { "pl7.app/label": "Reads", "pl7.app/min": "0", - "pl7.app/description": "Reads carrying this sequence, from the pre-refine pass -- before refine-tags would have snapped or dropped it." + "pl7.app/description": "Reads that carry this sequence. The pre-refine pass supplies the count. It is taken before refine-tags snapped the sequence onto a panel barcode, or dropped it." }) } }, @@ -1276,7 +1306,7 @@ undeclaredBarcodeImportSpec := func(sampleAxisSpec, tagAxisSpec) { "pl7.app/format": ".2p", "pl7.app/min": "0", "pl7.app/max": "1", - "pl7.app/description": "This one sequence's reads over every read the pre-refine pass saw for this sample. Reads here are not reads the run lost -- refine-tags snaps a sequence close enough to a panel barcode onto it, and the sample's own quality page carries the share it rescued that way." + "pl7.app/description": "The reads of this one sequence, divided by every read the pre-refine pass saw for this sample. The Status column reads this number. A read here is not always a read the run lost: the Reads correction rescued column gives the share refine-tags snapped onto the panel." }) } }, @@ -1287,11 +1317,11 @@ undeclaredBarcodeImportSpec := func(sampleAxisSpec, tagAxisSpec) { name: "pl7.app/antigen/undeclaredBarcodeShare", valueType: "Double", annotations: a(78000, true, { - "pl7.app/label": "Undeclared share (whole sample)", + "pl7.app/label": "Sample's undeclared share", "pl7.app/format": ".2p", "pl7.app/min": "0", "pl7.app/max": "1", - "pl7.app/description": "The share of this SAMPLE's reads landing in barcodes nobody declared -- every undeclared sequence together, whether or not it has a row here. Repeated on every one of that sample's rows, and the number the Status column reads. For this one sequence's own share, read the column beside it." + "pl7.app/description": "The share of this sample's reads that land in barcodes nobody declared. It counts every undeclared sequence together. It includes sequences that have no row in this table. The same number shows on every row of the sample. It carries no status. For the share of one sequence, read the column beside this one." }) } }, @@ -1305,7 +1335,7 @@ undeclaredBarcodeImportSpec := func(sampleAxisSpec, tagAxisSpec) { "pl7.app/label": "Status", "pl7.app/isDiscreteFilter": "true", "pl7.app/discreteValues": QC_STATUSES, - "pl7.app/description": "Reads the undeclared share of the whole sample, not this one sequence's: warns above a 0.50 share and alerts at 1.0 (inherited from the field), so it is the same on every row of a sample. This status is the barcode's, never a sample's: it never rolls into any sample's own status. No row here is the outcome to want -- an empty table means every barcode the reads carried was declared." + "pl7.app/description": "How this one sequence reads against its own share of the sample's reads. By default it warns above 1% and alerts above 5%. This status belongs to the barcode. It never rolls up into the status of a sample. No row in this table is the outcome to want. An empty table means every barcode the reads carried was declared." }) } } @@ -1537,19 +1567,19 @@ qcReagentImportSpec := func(panelAxisSpec, tagAxisSpec, identityAxisSpec) { // Seen in prints the ratio, because the count alone does not say out of how many. The two counts it is // built from stay as their own columns, optional, so the figure can be sorted and filtered as a number. col("seenIn", "pl7.app/antigen/reagentSeenIn", "String", "Seen in", 79000, true, - "Samples of this panel carrying any count of the barcode, over the samples it declares. 0/4 means no read carried it."), + "Samples of this panel where a listed cell carries any count of the barcode, over the samples the panel declares. 0/4 means none did."), col("samplesSeenIn", "pl7.app/antigen/reagentSamplesSeenIn", "Int", "Seen in (count)", 78950, false, - "Samples of this panel carrying any count of the barcode. Zero means no read carried it."), + "Samples of this panel where a listed cell carries any count of the barcode. Zero means none did."), col("samplesInPanel", "pl7.app/antigen/reagentSamplesInPanel", "Int", "Samples in panel", 78500, false, "The denominator for Seen in: the samples this panel declares."), col("samplesSeenInNames", "pl7.app/antigen/reagentSamplesSeenInNames", "String", "Seen in (samples)", 78900, false, - "Samples of this panel carrying any count of the barcode, named. Empty means none did."), + "Samples of this panel where a listed cell carries any count of the barcode, named. Empty means none did."), col("samplesInPanelNames", "pl7.app/antigen/reagentSamplesInPanelNames", "String", "Declared in (samples)", 78700, false, "The samples this panel declares the tag for, named. Reading it apart from Seen in (samples) tells a sample whose panel never declared the tag from one that declared it and saw nothing."), col("cellsWithCount", "pl7.app/antigen/reagentCellsWithCount", "Int", "Cells with any count", 78000, true, "Cells holding any count of the barcode, taken before the minimum."), col("cellsAboveTheLine", "pl7.app/antigen/reagentCellsAboveTheLine", "Double", "Cells called bound", 77000, true, - "Cells the verdict read called bound. Empty for a tag supplying the baseline; the reason column names that case."), + "Listed cells whose count of this one tag read bound, after the minimum. Empty for the baseline tag."), col("medianCountPerCell", "pl7.app/antigen/reagentMedianCount", "Double", "Median count per cell", 76000, true, "Median count over the cells holding any count, taken before the minimum."), // The two disagreement rates print a rate or, where none exists, the words for why. A single-tag @@ -1606,12 +1636,12 @@ qcSampleSummaryImportSpec := func(sampleAxisSpec) { "pl7.app/label": "Status", "pl7.app/isDiscreteFilter": "true", "pl7.app/discreteValues": QC_STATUSES, - "pl7.app/description": "The worst status among this sample's own measurements that carry one, from result_qc.csv's rollup row. Empty where nothing at this sample carried a status." + "pl7.app/description": "The worst status among this sample's own measurements that carry one. Empty where none does." }) } }, num("readsTotal", "readsTotal", "Reads parsed", 95000, undefined, - "Every read the parser saw, and the share matching the tag pattern."), + "Every read the parser saw, before pattern matching."), num("readsMatched", "readsMatched", "Reads matched", 94000, undefined, undefined), num("matchedFraction", "matchedFraction", "Matched fraction", 93000, ".2p", undefined), num("cellsDetected", "cellsDetected", "Cell barcodes detected", 92000, undefined, undefined), @@ -1633,11 +1663,11 @@ qcSampleSummaryImportSpec := func(sampleAxisSpec) { num("aggregateBarcodeFraction", "aggregateBarcodeFraction", "Aggregate-barcode read fraction", 83000, ".2p", "Reads in barcodes flagged as aggregates by the top-100 IQR rule, over readsTotal."), num("floorRemoved", "floorRemoved", "Counts removed by the minimum", 82000, undefined, - "Readings the minimum zeroed, and cells whose every non-reference reading was removed."), - num("uniqueCountsPerCell", "uniqueCountsPerCell", "Unique counts per cell", 81000, ".1f", - "Reads and distinct UMIs per cell barcode."), + "Readings the minimum count set to zero, summed over the sample. The baseline tag is exempt."), + num("uniqueCountsPerCell", "uniqueCountsPerCell", "Median unique counts per listed cell", 81000, ".1f", + "The median of each listed cell's total unique count across all tags. Cells outside the cell list are excluded."), num("highReferenceCells", "highReferenceCells", "Sticky reference cells", 80000, undefined, - "Cells whose reference reading exceeded the declared gate, or the spread of those readings.") + "With a gate set, how many cells read above it on the baseline tag. With no gate, the median baseline reading.") ], storageFormat: "Parquet", partitionKeyLength: 0 @@ -1669,6 +1699,7 @@ export { identityPunchImportSpec: identityPunchImportSpec, cellPunchImportSpec: cellPunchImportSpec, cellTagCountsImportSpec: cellTagCountsImportSpec, + cellRawTagCountsImportSpec: cellRawTagCountsImportSpec, cellScalarsImportSpec: cellScalarsImportSpec, offeredImportSpec: offeredImportSpec, tagIdentityLinkerImportSpec: tagIdentityLinkerImportSpec, diff --git a/workflow/src/main.tpl.tengo b/workflow/src/main.tpl.tengo index 326dfd0..4debb48 100644 --- a/workflow/src/main.tpl.tengo +++ b/workflow/src/main.tpl.tengo @@ -412,6 +412,8 @@ wf.body(func(args) { setParam("distributionMinCells", args.distributionMinCells) setParam("countFloor", args.countFloor) setParam("boundCutoff", args.boundCutoff) + setParam("boundProbability", args.boundProbability) + setParam("expectedBinderFraction", args.expectedBinderFraction) setParam("minVotingCells", args.minVotingCells) setParam("minAgreement", args.minAgreement) setParam("gateThreshold", args.gateThreshold) @@ -457,9 +459,9 @@ wf.body(func(args) { identitySummary: verdictRun.output("identitySummary"), identityPunch: verdictRun.output("identityPunch"), cellPunch: verdictRun.output("cellPunch"), - // cellCounts is deliberately NOT passed. The per-cell per-tag counts are the run's largest table and no - // reader exists for them on either side of the block boundary, so the import template does not build - // them into columns. An input it would ignore reads as a consumer that is not there. + // cellCounts, the POST-floor per-cell per-tag counts, is still deliberately NOT passed. + // cellRawCounts, the same grain BEFORE the floor, IS passed. + cellRawCounts: verdictRun.output("cellRawCounts"), cellScalars: verdictRun.output("cellScalars"), offered: verdictRun.output("offered"), tagIdentity: verdictRun.output("tagIdentity"), @@ -530,6 +532,8 @@ wf.body(func(args) { // per-cell reference readings are still reported, as an OUTPUT below. The per-cell per-tag counts are not // built at any grain. blockExports.antigenVerdicts = verdictImport.output("antigenVerdicts") + // Per-cell per-tag capture, pre-floor. + blockExports.antigenCellTagCounts = verdictImport.output("cellTagCounts") // The same frame as an OUTPUT, because a block's own exports are not in its own result pool. Without // this the block that produced the verdicts would be the one place that cannot read them. // diff --git a/workflow/src/verdict-args.lib.tengo b/workflow/src/verdict-args.lib.tengo index 1eef173..d2cf0da 100644 --- a/workflow/src/verdict-args.lib.tengo +++ b/workflow/src/verdict-args.lib.tengo @@ -25,6 +25,8 @@ PARAMETER_NAMES := [ "aggregateBarcodeError", "aggregateBarcodeWarn", "boundCutoff", + "boundProbability", + "expectedBinderFraction", "captureMap", "cellBarcodeValidError", "cellBarcodeValidWarn", @@ -52,6 +54,10 @@ PARAMETER_NAMES := [ // always on the command line, where the run record and a re-run both see it. DEFAULT_COUNT_FLOOR := 4 DEFAULT_BOUND_CUTOFF := 75 +// verdict.py DISTRIBUTION_BOUND_PROBABILITY. Both the default and the floor of what a run may ask for. +DEFAULT_BOUND_PROBABILITY := 0.9 +// The source paper's own initial weight. A statement about the experiment, so the form may move it. +DEFAULT_EXPECTED_BINDER_FRACTION := 0.1 DEFAULT_MIN_VOTING_CELLS := 1 DEFAULT_PANEL_MIN_MEMBERS := 25 DEFAULT_DISTRIBUTION_MIN_CELLS := 300 @@ -63,8 +69,8 @@ DEFAULT_CELL_BARCODE_VALID_ERROR := 0.50 DEFAULT_READS_PER_CELL_WARN := 5000 DEFAULT_AGGREGATE_BARCODE_WARN := 0.05 DEFAULT_AGGREGATE_BARCODE_ERROR := 1.0 -DEFAULT_UNDECLARED_BARCODE_WARN := 0.5 -DEFAULT_UNDECLARED_BARCODE_ERROR := 1.0 +DEFAULT_UNDECLARED_BARCODE_WARN := 0.01 +DEFAULT_UNDECLARED_BARCODE_ERROR := 0.05 DEFAULT_USABLE_READ_WARN := 0.20 DEFAULT_USABLE_READ_ERROR := 0.0 @@ -131,6 +137,8 @@ build := func(params) { add("--distribution-min-cells", _intArg(_num(params.distributionMinCells, DEFAULT_DISTRIBUTION_MIN_CELLS))) add("--floor", _intArg(_num(params.countFloor, DEFAULT_COUNT_FLOOR))) add("--cutoff", string(_num(params.boundCutoff, DEFAULT_BOUND_CUTOFF))) + add("--bound-probability", string(_num(params.boundProbability, DEFAULT_BOUND_PROBABILITY))) + add("--initial-signal-weight", string(_num(params.expectedBinderFraction, DEFAULT_EXPECTED_BINDER_FRACTION))) add("--min-voters", _intArg(_num(params.minVotingCells, DEFAULT_MIN_VOTING_CELLS))) // A switch, and it carries its value rather than standing alone. The vector is flag/value pairs: `valueOf` diff --git a/workflow/src/verdict-args.test.tengo b/workflow/src/verdict-args.test.tengo index 87f5bac..95d14ab 100644 --- a/workflow/src/verdict-args.test.tengo +++ b/workflow/src/verdict-args.test.tengo @@ -162,8 +162,8 @@ Test_line_thresholds_are_always_present := func() { test.isEqual("5000", va.valueOf(args, "--reads-per-cell-warn"), "the depth line is stated") test.isEqual("0.05", va.valueOf(args, "--aggregate-barcode-warn"), "the aggregate-barcode warn line is stated") test.isEqual("1", va.valueOf(args, "--aggregate-barcode-error"), "the aggregate-barcode alert line is stated") - test.isEqual("0.5", va.valueOf(args, "--undeclared-barcode-warn"), "the undeclared-barcode warn line is stated") - test.isEqual("1", va.valueOf(args, "--undeclared-barcode-error"), "the undeclared-barcode alert line is stated") + test.isEqual("0.01", va.valueOf(args, "--undeclared-barcode-warn"), "the undeclared-barcode warn line is stated") + test.isEqual("0.05", va.valueOf(args, "--undeclared-barcode-error"), "the undeclared-barcode alert line is stated") test.isEqual("0.2", va.valueOf(args, "--usable-read-warn"), "the usable-read warn line is stated") test.isEqual("0", va.valueOf(args, "--usable-read-error"), "the usable-read alert line is stated") } diff --git a/workflow/src/verdict-import.tpl.tengo b/workflow/src/verdict-import.tpl.tengo index fbe7c3a..96f80d3 100644 --- a/workflow/src/verdict-import.tpl.tengo +++ b/workflow/src/verdict-import.tpl.tengo @@ -254,6 +254,15 @@ self.body(func(inputs) { // (sampleId, tag) rather than any axis the tables above use, and because the retired check's // "undeclared-in-panel" direction was structurally unreachable there -- refine-tags has already snapped // every barcode onto the panel before that check's counts are built. + // Per-cell per-tag CAPTURE, exported for a downstream per-cell composition plot. Its own frame: the + // grain is (sample, cell, tag) and nothing else here shares it -- `cellRefFb` above is per cell, and the + // reagent table is per (tag, identity, panel). The tag label travels WITH it + cellTagFb := pframes.pFrameBuilder() + addTo(cellTagFb, "cellRawCounts", inputs.cellRawCounts, + columnSpecs.cellRawTagCountsImportSpec(sampleAxis, cellAxis, tagAxis)) + addTo(cellTagFb, "cellTagLabels", inputs.tagLabels, + columnSpecs.tagLabelsImportSpec(tagAxis)) + undeclaredFb := pframes.pFrameBuilder() addTo(undeclaredFb, "undeclaredBarcodes", inputs.undeclaredBarcodes, columnSpecs.undeclaredBarcodeImportSpec(sampleAxis, tagAxis)) @@ -266,6 +275,7 @@ self.body(func(inputs) { qcSummaryTable: qcSummaryFb.build(), qcDistributions: distributionFb.build(), reagentTable: reagentFb.build(), - undeclaredBarcodesTable: undeclaredFb.build() + undeclaredBarcodesTable: undeclaredFb.build(), + cellTagCounts: cellTagFb.build() } }) diff --git a/workflow/src/verdict-run.tpl.tengo b/workflow/src/verdict-run.tpl.tengo index ecd1ae7..8e744f7 100644 --- a/workflow/src/verdict-run.tpl.tengo +++ b/workflow/src/verdict-run.tpl.tengo @@ -30,6 +30,7 @@ RESULT_TABLES := [ { out: "identityPunch", file: "result_identity_punch.csv" }, { out: "cellPunch", file: "result_cell_punch.csv" }, { out: "cellCounts", file: "result_cell_counts.csv" }, + { out: "cellRawCounts", file: "result_cell_raw_counts.csv" }, { out: "cellScalars", file: "result_cell_scalars.csv" }, { out: "offered", file: "result_offered.csv" }, { out: "tagIdentity", file: "result_tag_identity.csv" },