From c45943b1964ae92e65389c64f9511b50f9d5649b Mon Sep 17 00:00:00 2001 From: Aaron Lisman Date: Sun, 18 Oct 2026 23:49:14 +0000 Subject: [PATCH 1/6] Risk model: share the modelled-patient filter, and keep each group's members MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Two internal changes Table 1 will build on, with no change to any result: - keepRiskRows(rows, spec) is the fit's own loop that decides who is modelled and counts who is left out and why, moved out of fitRiskmodel. A second reader of the same rows (Table 1) gets the same patients in the same order through it rather than through a copy. - groupRiskmodel also returns `members`: each group's patients as positions in the fit's subjects (no identifiers). It is not added to the view's result. Verified: test/riskmodel-groups.test.ts passes; the riskmodel result for msk_chord_2024 NSCLC (runs 2, covariate GENDER, stratified by SMOKING_PREDICTIONS_3_CLASSES; 7,809 patients) hashes identically (sha1 44e5f901…) before and after. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01UCpBeF9VMrjXyddFKUfiSX --- lib/views/riskmodel.ts | 54 ++++++++++++++++++++++++++---------------- 1 file changed, 34 insertions(+), 20 deletions(-) diff --git a/lib/views/riskmodel.ts b/lib/views/riskmodel.ts index 2e5d81a..ca98675 100644 --- a/lib/views/riskmodel.ts +++ b/lib/views/riskmodel.ts @@ -113,36 +113,48 @@ export class RiskRefusal extends Error { } } -/** Fetch the cohort's rows and fit the ensemble. Throws RiskRefusal when the data cannot support a model. */ -/** the resolver a fit reads through: a patient-level participant (viewRows) with view aggregates (viewMeta) */ -type RiskResolver = PatientLevelCapable & ViewAggregateCapable; +/** A modelled patient's row, as the fit reads it: by position, without any identifier. */ +export interface RiskKeptRow { t: number; entry: number; ev: boolean; stratum: string; genes: string[]; panels: string[]; cov: string[]; half: number } -export async function fitRiskmodel(r: RiskResolver, args: RiskFitArgs): Promise { - const s = riskFitSpec(args.view, args.compile.oncokbPin); - const mark = r.queries.length; - const [rows, panelRows] = await riskViewRows(r, args); - const logged = () => r.queries.slice(mark).map((q) => ({ label: q.label, sql: q.sql, ...(q.ms != null ? { ms: q.ms } : {}) })); - const panelGenes = new Map(panelRows.map((x) => [String(x.panel), new Set((x.genes as unknown as string[]) || [])] as const)); - // who is modelled, and who is left out and why — every exclusion is counted +/** + * Who of the rows `riskmodelRowsSql` returned is modelled, and who is left out and why — every + * exclusion counted. The fit and Table 1 both read their patients through this, so the same rows + * give the same patients in the same order. + */ +export function keepRiskRows(rows: Record[], s: Pick) { let notProfiled = 0, noStratum = 0, noCovariate = 0, noAnchor = 0, anchorAfterStart = 0, noFollowUp = 0; - const kept: { t: number; entry: number; ev: boolean; stratum: string; genes: string[]; panels: string[]; cov: string[]; half: number }[] = []; - for (const x of rows) { + const kept: RiskKeptRow[] = []; + rows.forEach((x) => { const panels = ((x.panels as unknown as string[]) || []).map(String); - if (!panels.length) { notProfiled++; continue; } + if (!panels.length) { notProfiled++; return; } const st = s.strata ? String(x.stratum ?? "") : ""; - if (s.strata && !st) { noStratum++; continue; } + if (s.strata && !st) { noStratum++; return; } const cov = s.covariates.map((_, k) => String(x["c" + k] ?? "")); - if (cov.some((v) => !v)) { noCovariate++; continue; } + if (cov.some((v) => !v)) { noCovariate++; return; } let t = Number(x.m), entry = 0; if (s.anchor) { - if (Number(x.anchored) !== 1) { noAnchor++; continue; } + if (Number(x.anchored) !== 1) { noAnchor++; return; } const off = Number(x.anchor_months); - if (off > 0) { anchorAfterStart++; continue; } + if (off > 0) { anchorAfterStart++; return; } entry = -off; t = t - off; - if (!(t > entry)) { noFollowUp++; continue; } + if (!(t > entry)) { noFollowUp++; return; } } kept.push({ t, entry, ev: Number(x.ev) === 1, stratum: st, genes: ((x.genes as unknown as string[]) || []).map(String), panels, cov, half: Number(x.half) }); - } + }); + return { kept, left: { notProfiled, noStratum, noCovariate, noAnchor, anchorAfterStart, noFollowUp } }; +} + +/** Fetch the cohort's rows and fit the ensemble. Throws RiskRefusal when the data cannot support a model. */ +/** the resolver a fit reads through: a patient-level participant (viewRows) with view aggregates (viewMeta) */ +type RiskResolver = PatientLevelCapable & ViewAggregateCapable; + +export async function fitRiskmodel(r: RiskResolver, args: RiskFitArgs): Promise { + const s = riskFitSpec(args.view, args.compile.oncokbPin); + const mark = r.queries.length; + const [rows, panelRows] = await riskViewRows(r, args); + const logged = () => r.queries.slice(mark).map((q) => ({ label: q.label, sql: q.sql, ...(q.ms != null ? { ms: q.ms } : {}) })); + const panelGenes = new Map(panelRows.map((x) => [String(x.panel), new Set((x.genes as unknown as string[]) || [])] as const)); + const { kept, left: { notProfiled, noStratum, noCovariate, noAnchor, anchorAfterStart, noFollowUp } } = keepRiskRows(rows, s); if (kept.length < 50) throw new RiskRefusal("has " + kept.length + " patients with a survival time, a mutation panel and every covariate on record; a risk model needs at least 50", logged(), true); // candidates: genes every modelled patient's panels cover (WES covers all), mutated often enough let covered = null as Set | null; @@ -237,6 +249,8 @@ export interface RiskGrouping { neverHeldOut: number; /** distinct held-out scores, and how many groups the cuts made */ distinct: number; groupsMade: number; + /** each group's patients, as positions in the fit's subjects (no identifiers); not part of the result */ + members: number[][]; } /** The quantiles a view's groups are cut at: `n` equal groups, or the quantiles it names. */ @@ -313,7 +327,7 @@ export function groupRiskmodel(fit: RiskFit, qs: number[]): RiskGrouping { return { groups, logrank, refit, histogram: { edges: Array.from({ length: bins + 1 }, (_, k) => lo + k * w), counts }, scale: { min: sorted[0], max: sorted[sorted.length - 1] }, percentiles: Array.from({ length: 101 }, (_, k) => sorted[Math.min(sorted.length - 1, Math.floor(k / 100 * sorted.length))]), - neverHeldOut: fit.pool.length - poolScored.length, distinct, groupsMade }; + neverHeldOut: fit.pool.length - poolScored.length, distinct, groupsMade, members }; } // ---- the view --------------------------------------------------------------------------------- From 979c76df16da68c3c2b2ada41b5ffc07c90117e1 Mon Sep 17 00:00:00 2001 From: Aaron Lisman Date: Mon, 19 Oct 2026 00:02:20 +0000 Subject: [PATCH 2/6] Risk model: Table 1, each risk group's clinical and genomic makeup MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A riskmodel view with `describe: {attributes?, genes?}` adds `table1` to its result: per risk group, a categorical attribute's patients per value, a numeric one's median and IQR, a gene's mutated patients over those whose panels cover it — each tested across groups, with Benjamini–Hochberg q. Where and how it is counted - In the engine, from the fit's own rows query (riskmodelRowsSql) with the attributes added (folded value, its spelling, and how many values the patient's samples hold), kept by keepRiskRows exactly as the fit kept them, and refused unless they match the fit's patients one by one (time, entry, event). Group membership is groupRiskmodel's `members`. Only counts leave; no identifier is projected (boundary suite passes). - Gene rows are panel-aware: profiled = the patient's mutation panels cover the gene (WES covers all); a patient whose panels do not is "not profiled", never unmutated. Mutated follows the model's `measure`. - Tests reuse lib/stats.ts: chiSquareTest (categories; genes as mutated vs not among the profiled), kruskalWallis + quantilesFromHistogram (numbers), benjaminiHochberg. - Rows the groups differ in by construction are marked: genes the model kept in any run and covariates "in the model", the endpoint's columns "outcome". They keep p but get no q and stay out of the BH family. The stratum is marked "stratum" and stays in (the score has no stratum term). - Defaults: attributes = stratum, covariates, then the study's attributes by priority without survival columns; of the first 20, those with no value, one value, or >12 categories are left out (listed in `omitted` with why), first 8 kept. Genes = the 15 kept in the most runs. Named lists: ≤20 each (refused past that); an unknown attribute, a category with >100 values, or a gene on no panel and mutated in nobody is set aside by name. Case variants of one gene are one row. - Not in the fit's cache key (test/riskmodel-groups.test.ts): changing Table 1 refits nothing. Invariants: riskmodelMethod gains Table 1's steps and a script recomputing every p, q and median from the counts (scipy chi2_contingency without correction, kruskal, numpy percentile, statsmodels fdr_bh over the same family); lib/verification.ts names Table 1 under χ² and Kruskal–Wallis (no new statistic); DSL reference + docs/dsl.md regenerated; the gallery's riskmodel example names a Table 1; two planted bugs in check-mutants (denominator ignores the panel; model rows in the FDR family). Verified - npm test: 1258 pass. test/riskmodel-table1.test.ts covers the denominator (panel lacking the gene, WES, a group with nobody covered), roles and the BH family, folding, numbers, defaults, set-aside rows, and the refusal when rows do not match the fit. Planting either mutant by hand fails it. - Against ClickHouse, independently: msk_chord_2024 NSCLC, validate, runs 3 — a hand-written query over the validation half gives 3867 modelled patients, SMARCA2 22/1093, ZFHX3 185/3710, KEAP1 535/3867; Table 1's per-group counts sum to exactly these. - check:methods on the riskmodel example: the script reproduces all 21 numbers. Ad-hoc runs with defaults, named lists, strata/covariates and validation all reproduced (29–32 numbers each); the check caught a duplicate-gene bug, fixed above. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01UCpBeF9VMrjXyddFKUfiSX --- docs/dsl.md | 19 +++- lib/dsl-reference.ts | 15 +-- lib/federation-views.ts | 25 ++++- lib/method.ts | 109 ++++++++++++++++++- lib/types.ts | 50 +++++++++ lib/verification.ts | 4 +- lib/views/riskmodel-table1.ts | 196 ++++++++++++++++++++++++++++++++++ lib/views/riskmodel.ts | 90 +++++++++++++--- scripts/check-mutants.ts | 8 ++ test/riskmodel-groups.test.ts | 3 +- test/riskmodel-table1.test.ts | 157 +++++++++++++++++++++++++++ 11 files changed, 646 insertions(+), 30 deletions(-) create mode 100644 lib/views/riskmodel-table1.ts create mode 100644 test/riskmodel-table1.test.ts diff --git a/docs/dsl.md b/docs/dsl.md index a092dc5..ccb3bd5 100644 --- a/docs/dsl.md +++ b/docs/dsl.md @@ -244,7 +244,7 @@ patient — so the arm counts summing to more than the base cohort is expected, | `paired` | `genes?:[≤20], top?:1-40, pairBy?:{attribute, a, b, studies?:{:{attribute, a, b}}, labels?:{a, b}}, measure?:"driver"|"any", cohorts:[names] (≤3)` | **paired samples within patients**: among patients with at least one sample on each side of `pairBy` (a SAMPLE-level attribute and its values for side a and side b — default `SAMPLE_TYPE` Primary / Metastasis; values fold on case and whitespace like every clinical predicate), per gene, the patients with a driver (or, with `measure:"any"`, any alteration) on both sides, on side a only, on side b only — counted only among patients profiled for that gene on BOTH sides by each side's own panels, so a gene one side never assayed never reads as discordant. A patient with several samples on a side is altered there if any is. Concordance = both ÷ altered anywhere; an exact McNemar test (two-sided binomial on the discordant patients) asks whether discordance runs one way, BH q across the genes of each cohort. When no patient has both sides, the refusal lists the values the attribute actually takes. Sample type records where tissue came from, not when — the caption says "metastasis-only" is not "acquired later". Counts only; a remote peer cannot answer yet and is named. **Without `genes` it scans**: every gene with a driver (or any alteration) in some paired patient is tested, q is Benjamini–Hochberg across all of them, and the `top` (default 20) most one-sided are drawn, ordered by McNemar p — how to ask "what is acquired after treatment" without naming genes first. **`pairBy.studies`** gives each study its own attribute and values (PRE-TX / POST-TX in one, pre-treatment / post-treatment in another) so their patients pool into one count; a study not listed uses the top-level attribute, and `labels` names the two sides in the figure (default: their values). With an attribute other than a sample site, the caption says the sides are as each study labels them, and that a side-only call can also be a subclone one biopsy missed or a lower-purity sample under-calling | | `clonality` | `genes:[≤12], measure?:"driver"|"any", qc?:true, cohorts:[names] (≤3)` | **clonality from FACETS allele-specific copy number** (only studies that carry it — the refusal names them): per gene, the cohort's mutation CALLS (one mutation in one sample; OncoKB drivers by default) that FACETS called clonal, subclonal or indeterminate, LOH at the locus (minor copy number 0, of calls where it was estimated) and the cancer-cell fraction in ten bins; per cohort, FACETS's own whole-genome-doubling, purity and ploidy calls. `qc` (default true) uses only samples whose `FACETS_QC` is Pass, with the pass/fail counts reported. Only the study's current mutation profile is read (the table also holds rows orphaned by superseded loads). Two cohorts add Fisher tests on clonal vs subclonal per gene (BH across genes) and on WGD. Counts only; a remote peer cannot answer yet and is named | | `incidence` | `event:{eventType, attrs?:[{key,value} ×≤3], label?}, anchor?:{eventType,key?,value?,occurrence?}, horizons?:[12,24,60], cohorts:[names] (≤6)` | **cumulative incidence with a competing risk**: time from the anchor (default the study's day 0) to the FIRST `event` after it, with death (`OS_STATUS`) as the competing event and `OS_MONTHS` as follow-up — Aalen–Johansen, with 95% intervals from its delta-method standard error (checked against R's cmprsk: estimates identical to five decimals, SEs within 3%). Not 1 − Kaplan–Meier, which censors the dead as if still at risk and overstates incidence wherever death is common; the view reports both at each of `horizons`, so the difference is on the page. Left truncation as in `survival`: a patient anchored before day 0 enters the risk set at day 0. Patients who cannot be at risk are counted by reason, never read as event-free: already had the event at or before the anchor, had it before observation began, no follow-up, no anchor event. Two or more cohorts: Gray's test (weight 1) by permutation of cohort labels over the pooled records — exact under the null, seeded, and computable from the (entry, time, cause) tallies each site sends, which are counts | -| `riskmodel` | `cohorts:[one name], endpoint?, monthsAttribute?/statusAttribute?, measure?:"any"|"driver", minMutated?:30, strata?:"", covariates?:["" ×≤6], anchor?:{eventType,key?,value?,occurrence?}, validate?:false, alpha?:0.5, runs?:10 (≤100), groups?:3 (2-5) | cuts?:[quantiles ×1–4], seed?:1, labelAbove?:0.5` | **a survival risk model from mutated genes**, in the manner of OncoCast (Martin et al.): an elastic-net Cox model (glmnet's algorithm, checked against scikit-survival; Breslow ties; `alpha` the l1 share of the penalty) of one cohort's survival on every gene that the panels of EVERY patient cover — so no patient counts as wild-type for a gene never tested — and that is mutated in at least `minMutated` patients. Fitted `runs` times, each on a random two-thirds of the patients with the penalty chosen by 3-fold cross-validated partial likelihood, scoring the held-out third; a patient's risk score is the mean of their held-out scores, so no one is scored by a model trained on them (a patient never held out has no score and is counted). `covariates` are clinical attributes always in the model and never penalised (numbers as they are, categories as indicators against the commonest value, at most 10 values); a patient missing one is left out and counted. `anchor` starts each patient's clock at a timeline event, as the survival view does, with the delayed entry that implies (patients without the event, or whose event follows the study's day 0, are counted out). `validate` trains only on the discovery half of the engine's fixed split (the `split` leaf's halves) and reports everything below on the validation half, which no model saw, cut at the discovery half's quantiles. Reports the held-out C-index across runs (and on the validation half), each gene's share of runs kept and mean hazard ratio, the covariates' hazard ratios, Kaplan-Meier curves of `groups` risk groups cut at equal quantiles of the score (no outcome used; a cut moves to the nearer edge of a run of tied scores, so equal scores share a group, and groups merge when the scores cannot fill them) with a log-rank test — or, with `cuts` instead of `groups` (not both), cut at the quantiles it names, increasing, each strictly between 0 and 1, as OncoCast's `cuts` (e.g. `[0.8]`, the top fifth against the rest; tied scores still share a group, and the caveat says the cuts were chosen, since cuts placed by looking at the curves overstate the separation — confirm one on the validation half) — each group's 12- and 36-month survival and hazard ratio against the lowest (95% Wald interval), the hazard ratio per standard deviation of the score, and the score's histogram. The score is shown on OncoCast's 0–10 scale as well as raw: (score − min) / (max − min) × 10, min and max the lowest and highest held-out score the cuts are placed on (with `validate`, the discovery half's, so a validation patient can fall outside 0–10); the scale is linear, so it changes no group. The genes are also drawn as OncoCast's volcano — mean hazard ratio where kept (log scale) against the share of runs kept, sized by patients mutated — naming those kept in at least `labelAbove` of the runs (0–1; drawn in the browser, so changing it re-runs nothing). `strata` gives each value its own baseline hazard — use CANCER_TYPE on a mixed cohort, or a gene marking a type carries that type's prognosis. One cohort in one patient-level study; at most 300,000 patients × runs per request (refused, naming the limit). Seeded: the same document gives the same model | +| `riskmodel` | `cohorts:[one name], endpoint?, monthsAttribute?/statusAttribute?, measure?:"any"|"driver", minMutated?:30, strata?:"", covariates?:["" ×≤6], anchor?:{eventType,key?,value?,occurrence?}, validate?:false, alpha?:0.5, runs?:10 (≤100), groups?:3 (2-5) | cuts?:[quantiles ×1–4], seed?:1, labelAbove?:0.5, describe?:{attributes?:["" ×≤20], genes?:["" ×≤20]}` | **a survival risk model from mutated genes**, in the manner of OncoCast (Martin et al.): an elastic-net Cox model (glmnet's algorithm, checked against scikit-survival; Breslow ties; `alpha` the l1 share of the penalty) of one cohort's survival on every gene that the panels of EVERY patient cover — so no patient counts as wild-type for a gene never tested — and that is mutated in at least `minMutated` patients. Fitted `runs` times, each on a random two-thirds of the patients with the penalty chosen by 3-fold cross-validated partial likelihood, scoring the held-out third; a patient's risk score is the mean of their held-out scores, so no one is scored by a model trained on them (a patient never held out has no score and is counted). `covariates` are clinical attributes always in the model and never penalised (numbers as they are, categories as indicators against the commonest value, at most 10 values); a patient missing one is left out and counted. `anchor` starts each patient's clock at a timeline event, as the survival view does, with the delayed entry that implies (patients without the event, or whose event follows the study's day 0, are counted out). `validate` trains only on the discovery half of the engine's fixed split (the `split` leaf's halves) and reports everything below on the validation half, which no model saw, cut at the discovery half's quantiles. Reports the held-out C-index across runs (and on the validation half), each gene's share of runs kept and mean hazard ratio, the covariates' hazard ratios, Kaplan-Meier curves of `groups` risk groups cut at equal quantiles of the score (no outcome used; a cut moves to the nearer edge of a run of tied scores, so equal scores share a group, and groups merge when the scores cannot fill them) with a log-rank test — or, with `cuts` instead of `groups` (not both), cut at the quantiles it names, increasing, each strictly between 0 and 1, as OncoCast's `cuts` (e.g. `[0.8]`, the top fifth against the rest; tied scores still share a group, and the caveat says the cuts were chosen, since cuts placed by looking at the curves overstate the separation — confirm one on the validation half) — each group's 12- and 36-month survival and hazard ratio against the lowest (95% Wald interval), the hazard ratio per standard deviation of the score, and the score's histogram. The score is shown on OncoCast's 0–10 scale as well as raw: (score − min) / (max − min) × 10, min and max the lowest and highest held-out score the cuts are placed on (with `validate`, the discovery half's, so a validation patient can fall outside 0–10); the scale is linear, so it changes no group. The genes are also drawn as OncoCast's volcano — mean hazard ratio where kept (log scale) against the share of runs kept, sized by patients mutated — naming those kept in at least `labelAbove` of the runs (0–1; drawn in the browser, so changing it re-runs nothing). `strata` gives each value its own baseline hazard — use CANCER_TYPE on a mixed cohort, or a gene marking a type carries that type's prognosis. `describe` adds **Table 1**, each risk group's makeup, counted on the patients the curves draw (with `validate`, the validation half): a categorical attribute's patients per value and group (the patient's value, or their first sample's, case and whitespace folded; those whose samples disagree are counted), with a χ² test; a numeric one's median and IQR per group, with Kruskal–Wallis; a gene's mutated patients over those whose mutation panels cover it (`measure` decides what counts as mutated; patients whose panels do not are "not profiled", never unmutated), with a χ² test; Benjamini–Hochberg q-values across rows. Rows the groups differ in by construction are marked — the model's genes (kept in any run) and covariates "in the model", the endpoint's own columns "outcome" — and get a p-value but no q-value, staying out of the correction; the model's stratum is marked "stratum" and stays in, since the score has no stratum term. A list left out takes its default: `attributes`, the model's stratum and covariates then the study's attributes in priority order, survival columns left out, keeping the first 8 that have a value here, take more than one, and (categorical) at most 12 values — those passed over are listed with why; `genes`, the 15 the model kept in the most runs. A named attribute the study lacks, a categorical one with more than 100 values here, and a gene on none of the patients' panels and mutated in none are set aside by name. `describe: {}` is the defaults. One cohort in one patient-level study; at most 300,000 patients × runs per request (refused, naming the limit). Seeded: the same document gives the same model | | `labs` | `test, anchor:{eventType,key?,value?,occurrence?}, eventType?:"Lab_Test", testKey?:"TEST", resultKey?:"RESULT", baselineDays?:90, fromWeeks?:-12, toWeeks?:52, binWeeks?:4, responseDrop?:50, cohorts:[names] (≤4)` | **a lab value's course from a treatment start**: each patient's `test` results (timeline events of `eventType` whose `testKey` names the test, value in `resultKey`) as % change from their own baseline — the last value in the `baselineDays` up to and including the anchor; a non-positive baseline is left out — per `binWeeks` bin from `fromWeeks` to `toWeeks`, drawn as median and interquartile band (from 2-point histograms, exact to within a point; a patient counts once per bin, at their mean; the baseline reading itself is not drawn). Best response is each patient's lowest value after the anchor; the share falling by at least `responseDrop`% is counted exactly (50 on PSA is PCWG3's PSA50), and two cohorts get a Fisher test on it. A **waterfall** of every patient's best change is drawn only from participants entitled to patient-level values (a method only PatientLevelResolver has, like a plot's points); an aggregate-only participant's patients stay in the histograms and the rate and are counted as withheld from the bars. With no evaluable patient, the refusal lists the tests the study records | | `cancertypes` | `genes:[≤30], groupBy?:"CANCER_TYPE"||"study", measure?:"any"|"driver", minSamples?:10, show?:"frequency"|"count", unit?:"patient"|"sample", cohorts:[names] (≤2)` | **Cancer Types Summary**, cBioPortal's results-view tab: per group of the cohort's PATIENTS (the portal's default; `unit:"sample"` counts samples — a patient counts in each group one of their samples is in) — a sample attribute's value (folded on case and whitespace like every clinical predicate; samples with no value form their own group) or the study — the share altered in at least one of `genes`, split by kind: only mutations, only amplifications, only deep deletions, only structural variants, or more than one kind. A group's denominator is its samples profiled for at least one of the genes (the union of their panels, as an oncoprint's merged row counts it); its total sample count is shown beside it. `measure:"driver"` counts OncoKB drivers only. Groups with fewer than `minSamples` profiled samples are pooled into one bar, named and counted, never dropped. Counts only, so every local participant can answer; a remote peer cannot yet and is named | | `overlap` | `cohorts:[names] (2–6)` | **how cohorts' patients intersect**, cBioPortal's group-comparison Overlap tab: for every combination of the cohorts, the patients — and their samples — in exactly those cohorts and none of the others, with each cohort's own size. Two or three cohorts are drawn as the portal draws them, area-proportional diagrams of samples and of patients side by side; more as an UpSet plot. A cohort is a set of patients, so a patient's samples go wherever the patient does. Each study's database answers with one count per combination (a bitmask of membership) — counts only — and a cohort's size here must equal its own count, or the view says the definition gave two answers. Ids are study-qualified, so patients only ever overlap within a study | @@ -3865,7 +3865,7 @@ Prostate cancer in MSK-CHORD, from the first hormone therapy. PSA falls by a med ### Example 56: A survival risk model from mutated genes in lung cancer -Non-small cell lung cancer in MSK-CHORD. The model keeps KEAP1, STK11, SMARCA4 and TP53 as adverse and EGFR as favourable, as the literature does, and its risk thirds separate clearly; a held-out C-index near 0.61 says mutations alone predict survival only modestly. Five runs keep the example quick; OncoCast recommends fifty. +Non-small cell lung cancer in MSK-CHORD. The model keeps KEAP1, STK11, SMARCA4 and TP53 as adverse and EGFR as favourable, as the literature does, and its risk thirds separate clearly; a held-out C-index near 0.61 says mutations alone predict survival only modestly. Five runs keep the example quick; OncoCast recommends fifty. Its Table 1 (`describe`) shows who is in each third: the high-risk third is more often male, a former or current smoker and at stage 4 than the low-risk third. KEAP1 and TP53 are marked “in the model” — the score is built from them, so their enrichment is by construction and gets no q-value — and ZFHX3, which not every panel covers, is counted over each group's patients whose panel does. ```json { @@ -3890,7 +3890,20 @@ Non-small cell lung cancer in MSK-CHORD. The model keeps KEAP1, STK11, SMARCA4 a "cohorts": [ "NSCLC" ], - "runs": 5 + "runs": 5, + "describe": { + "attributes": [ + "STAGE_HIGHEST_RECORDED", + "GENDER", + "SMOKING_PREDICTIONS_3_CLASSES", + "CURRENT_AGE_DEID" + ], + "genes": [ + "KEAP1", + "TP53", + "ZFHX3" + ] + } } ] } diff --git a/lib/dsl-reference.ts b/lib/dsl-reference.ts index 956622e..a87d3f3 100644 --- a/lib/dsl-reference.ts +++ b/lib/dsl-reference.ts @@ -998,15 +998,17 @@ export const EXAMPLES: DslExample[] = [ }, { title: "A survival risk model from mutated genes in lung cancer", - demonstrates: "A `riskmodel` view: an elastic-net Cox model of survival on every mutated gene that all the cohort's panels cover, fitted on random two-thirds of the patients and scored on the held-out third, repeatedly (OncoCast's ensemble) — so each patient's risk score comes from models that never saw them — with the held-out C-index, the genes kept most often, and Kaplan-Meier curves of the risk groups.", + demonstrates: "A `riskmodel` view: an elastic-net Cox model of survival on every mutated gene that all the cohort's panels cover, fitted on random two-thirds of the patients and scored on the held-out third, repeatedly (OncoCast's ensemble) — so each patient's risk score comes from models that never saw them — with the held-out C-index, the genes kept most often, Kaplan-Meier curves of the risk groups, and each group's clinical and genomic makeup (Table 1).", tags: ["riskmodel", "machine learning", "survival", "cox"], - note: "Non-small cell lung cancer in MSK-CHORD. The model keeps KEAP1, STK11, SMARCA4 and TP53 as adverse and EGFR as favourable, as the literature does, and its risk thirds separate clearly; a held-out C-index near 0.61 says mutations alone predict survival only modestly. Five runs keep the example quick; OncoCast recommends fifty.", + note: "Non-small cell lung cancer in MSK-CHORD. The model keeps KEAP1, STK11, SMARCA4 and TP53 as adverse and EGFR as favourable, as the literature does, and its risk thirds separate clearly; a held-out C-index near 0.61 says mutations alone predict survival only modestly. Five runs keep the example quick; OncoCast recommends fifty. " + + "Its Table 1 (`describe`) shows who is in each third: the high-risk third is more often male, a former or current smoker and at stage 4 than the low-risk third. KEAP1 and TP53 are marked “in the model” — the score is built from them, so their enrichment is by construction and gets no q-value — and ZFHX3, which not every panel covers, is counted over each group's patients whose panel does.", doc: { study: "msk_chord_2024", cohorts: { "NSCLC": { all: [{ kind: "clinical", attribute: "CANCER_TYPE", values: ["Non-Small Cell Lung Cancer"], per: "sample" }] } }, - views: [{ kind: "riskmodel", cohorts: ["NSCLC"], runs: 5 }] + views: [{ kind: "riskmodel", cohorts: ["NSCLC"], runs: 5, + describe: { attributes: ["STAGE_HIGHEST_RECORDED", "GENDER", "SMOKING_PREDICTIONS_3_CLASSES", "CURRENT_AGE_DEID"], genes: ["KEAP1", "TP53", "ZFHX3"] } }] } }, { @@ -1502,7 +1504,7 @@ patient — so the arm counts summing to more than the base cohort is expected, | \`paired\` | \`genes?:[≤20], top?:1-40, pairBy?:{attribute, a, b, studies?:{:{attribute, a, b}}, labels?:{a, b}}, measure?:"driver"\|"any", cohorts:[names] (≤3)\` | **paired samples within patients**: among patients with at least one sample on each side of \`pairBy\` (a SAMPLE-level attribute and its values for side a and side b — default \`SAMPLE_TYPE\` Primary / Metastasis; values fold on case and whitespace like every clinical predicate), per gene, the patients with a driver (or, with \`measure:"any"\`, any alteration) on both sides, on side a only, on side b only — counted only among patients profiled for that gene on BOTH sides by each side's own panels, so a gene one side never assayed never reads as discordant. A patient with several samples on a side is altered there if any is. Concordance = both ÷ altered anywhere; an exact McNemar test (two-sided binomial on the discordant patients) asks whether discordance runs one way, BH q across the genes of each cohort. When no patient has both sides, the refusal lists the values the attribute actually takes. Sample type records where tissue came from, not when — the caption says "metastasis-only" is not "acquired later". Counts only; a remote peer cannot answer yet and is named. **Without \`genes\` it scans**: every gene with a driver (or any alteration) in some paired patient is tested, q is Benjamini–Hochberg across all of them, and the \`top\` (default 20) most one-sided are drawn, ordered by McNemar p — how to ask "what is acquired after treatment" without naming genes first. **\`pairBy.studies\`** gives each study its own attribute and values (PRE-TX / POST-TX in one, pre-treatment / post-treatment in another) so their patients pool into one count; a study not listed uses the top-level attribute, and \`labels\` names the two sides in the figure (default: their values). With an attribute other than a sample site, the caption says the sides are as each study labels them, and that a side-only call can also be a subclone one biopsy missed or a lower-purity sample under-calling | | \`clonality\` | \`genes:[≤12], measure?:"driver"\|"any", qc?:true, cohorts:[names] (≤3)\` | **clonality from FACETS allele-specific copy number** (only studies that carry it — the refusal names them): per gene, the cohort's mutation CALLS (one mutation in one sample; OncoKB drivers by default) that FACETS called clonal, subclonal or indeterminate, LOH at the locus (minor copy number 0, of calls where it was estimated) and the cancer-cell fraction in ten bins; per cohort, FACETS's own whole-genome-doubling, purity and ploidy calls. \`qc\` (default true) uses only samples whose \`FACETS_QC\` is Pass, with the pass/fail counts reported. Only the study's current mutation profile is read (the table also holds rows orphaned by superseded loads). Two cohorts add Fisher tests on clonal vs subclonal per gene (BH across genes) and on WGD. Counts only; a remote peer cannot answer yet and is named | | \`incidence\` | \`event:{eventType, attrs?:[{key,value} ×≤3], label?}, anchor?:{eventType,key?,value?,occurrence?}, horizons?:[12,24,60], cohorts:[names] (≤6)\` | **cumulative incidence with a competing risk**: time from the anchor (default the study's day 0) to the FIRST \`event\` after it, with death (\`OS_STATUS\`) as the competing event and \`OS_MONTHS\` as follow-up — Aalen–Johansen, with 95% intervals from its delta-method standard error (checked against R's cmprsk: estimates identical to five decimals, SEs within 3%). Not 1 − Kaplan–Meier, which censors the dead as if still at risk and overstates incidence wherever death is common; the view reports both at each of \`horizons\`, so the difference is on the page. Left truncation as in \`survival\`: a patient anchored before day 0 enters the risk set at day 0. Patients who cannot be at risk are counted by reason, never read as event-free: already had the event at or before the anchor, had it before observation began, no follow-up, no anchor event. Two or more cohorts: Gray's test (weight 1) by permutation of cohort labels over the pooled records — exact under the null, seeded, and computable from the (entry, time, cause) tallies each site sends, which are counts | -| \`riskmodel\` | \`cohorts:[one name], endpoint?, monthsAttribute?/statusAttribute?, measure?:"any"\|"driver", minMutated?:30, strata?:"", covariates?:["" ×≤6], anchor?:{eventType,key?,value?,occurrence?}, validate?:false, alpha?:0.5, runs?:10 (≤100), groups?:3 (2-5) | cuts?:[quantiles ×1–4], seed?:1, labelAbove?:0.5\` | **a survival risk model from mutated genes**, in the manner of OncoCast (Martin et al.): an elastic-net Cox model (glmnet's algorithm, checked against scikit-survival; Breslow ties; \`alpha\` the l1 share of the penalty) of one cohort's survival on every gene that the panels of EVERY patient cover — so no patient counts as wild-type for a gene never tested — and that is mutated in at least \`minMutated\` patients. Fitted \`runs\` times, each on a random two-thirds of the patients with the penalty chosen by 3-fold cross-validated partial likelihood, scoring the held-out third; a patient's risk score is the mean of their held-out scores, so no one is scored by a model trained on them (a patient never held out has no score and is counted). \`covariates\` are clinical attributes always in the model and never penalised (numbers as they are, categories as indicators against the commonest value, at most 10 values); a patient missing one is left out and counted. \`anchor\` starts each patient's clock at a timeline event, as the survival view does, with the delayed entry that implies (patients without the event, or whose event follows the study's day 0, are counted out). \`validate\` trains only on the discovery half of the engine's fixed split (the \`split\` leaf's halves) and reports everything below on the validation half, which no model saw, cut at the discovery half's quantiles. Reports the held-out C-index across runs (and on the validation half), each gene's share of runs kept and mean hazard ratio, the covariates' hazard ratios, Kaplan-Meier curves of \`groups\` risk groups cut at equal quantiles of the score (no outcome used; a cut moves to the nearer edge of a run of tied scores, so equal scores share a group, and groups merge when the scores cannot fill them) with a log-rank test — or, with \`cuts\` instead of \`groups\` (not both), cut at the quantiles it names, increasing, each strictly between 0 and 1, as OncoCast's \`cuts\` (e.g. \`[0.8]\`, the top fifth against the rest; tied scores still share a group, and the caveat says the cuts were chosen, since cuts placed by looking at the curves overstate the separation — confirm one on the validation half) — each group's 12- and 36-month survival and hazard ratio against the lowest (95% Wald interval), the hazard ratio per standard deviation of the score, and the score's histogram. The score is shown on OncoCast's 0–10 scale as well as raw: (score − min) / (max − min) × 10, min and max the lowest and highest held-out score the cuts are placed on (with \`validate\`, the discovery half's, so a validation patient can fall outside 0–10); the scale is linear, so it changes no group. The genes are also drawn as OncoCast's volcano — mean hazard ratio where kept (log scale) against the share of runs kept, sized by patients mutated — naming those kept in at least \`labelAbove\` of the runs (0–1; drawn in the browser, so changing it re-runs nothing). \`strata\` gives each value its own baseline hazard — use CANCER_TYPE on a mixed cohort, or a gene marking a type carries that type's prognosis. One cohort in one patient-level study; at most 300,000 patients × runs per request (refused, naming the limit). Seeded: the same document gives the same model | +| \`riskmodel\` | \`cohorts:[one name], endpoint?, monthsAttribute?/statusAttribute?, measure?:"any"\|"driver", minMutated?:30, strata?:"", covariates?:["" ×≤6], anchor?:{eventType,key?,value?,occurrence?}, validate?:false, alpha?:0.5, runs?:10 (≤100), groups?:3 (2-5) | cuts?:[quantiles ×1–4], seed?:1, labelAbove?:0.5, describe?:{attributes?:["" ×≤20], genes?:["" ×≤20]}\` | **a survival risk model from mutated genes**, in the manner of OncoCast (Martin et al.): an elastic-net Cox model (glmnet's algorithm, checked against scikit-survival; Breslow ties; \`alpha\` the l1 share of the penalty) of one cohort's survival on every gene that the panels of EVERY patient cover — so no patient counts as wild-type for a gene never tested — and that is mutated in at least \`minMutated\` patients. Fitted \`runs\` times, each on a random two-thirds of the patients with the penalty chosen by 3-fold cross-validated partial likelihood, scoring the held-out third; a patient's risk score is the mean of their held-out scores, so no one is scored by a model trained on them (a patient never held out has no score and is counted). \`covariates\` are clinical attributes always in the model and never penalised (numbers as they are, categories as indicators against the commonest value, at most 10 values); a patient missing one is left out and counted. \`anchor\` starts each patient's clock at a timeline event, as the survival view does, with the delayed entry that implies (patients without the event, or whose event follows the study's day 0, are counted out). \`validate\` trains only on the discovery half of the engine's fixed split (the \`split\` leaf's halves) and reports everything below on the validation half, which no model saw, cut at the discovery half's quantiles. Reports the held-out C-index across runs (and on the validation half), each gene's share of runs kept and mean hazard ratio, the covariates' hazard ratios, Kaplan-Meier curves of \`groups\` risk groups cut at equal quantiles of the score (no outcome used; a cut moves to the nearer edge of a run of tied scores, so equal scores share a group, and groups merge when the scores cannot fill them) with a log-rank test — or, with \`cuts\` instead of \`groups\` (not both), cut at the quantiles it names, increasing, each strictly between 0 and 1, as OncoCast's \`cuts\` (e.g. \`[0.8]\`, the top fifth against the rest; tied scores still share a group, and the caveat says the cuts were chosen, since cuts placed by looking at the curves overstate the separation — confirm one on the validation half) — each group's 12- and 36-month survival and hazard ratio against the lowest (95% Wald interval), the hazard ratio per standard deviation of the score, and the score's histogram. The score is shown on OncoCast's 0–10 scale as well as raw: (score − min) / (max − min) × 10, min and max the lowest and highest held-out score the cuts are placed on (with \`validate\`, the discovery half's, so a validation patient can fall outside 0–10); the scale is linear, so it changes no group. The genes are also drawn as OncoCast's volcano — mean hazard ratio where kept (log scale) against the share of runs kept, sized by patients mutated — naming those kept in at least \`labelAbove\` of the runs (0–1; drawn in the browser, so changing it re-runs nothing). \`strata\` gives each value its own baseline hazard — use CANCER_TYPE on a mixed cohort, or a gene marking a type carries that type's prognosis. \`describe\` adds **Table 1**, each risk group's makeup, counted on the patients the curves draw (with \`validate\`, the validation half): a categorical attribute's patients per value and group (the patient's value, or their first sample's, case and whitespace folded; those whose samples disagree are counted), with a χ² test; a numeric one's median and IQR per group, with Kruskal–Wallis; a gene's mutated patients over those whose mutation panels cover it (\`measure\` decides what counts as mutated; patients whose panels do not are "not profiled", never unmutated), with a χ² test; Benjamini–Hochberg q-values across rows. Rows the groups differ in by construction are marked — the model's genes (kept in any run) and covariates "in the model", the endpoint's own columns "outcome" — and get a p-value but no q-value, staying out of the correction; the model's stratum is marked "stratum" and stays in, since the score has no stratum term. A list left out takes its default: \`attributes\`, the model's stratum and covariates then the study's attributes in priority order, survival columns left out, keeping the first 8 that have a value here, take more than one, and (categorical) at most 12 values — those passed over are listed with why; \`genes\`, the 15 the model kept in the most runs. A named attribute the study lacks, a categorical one with more than 100 values here, and a gene on none of the patients' panels and mutated in none are set aside by name. \`describe: {}\` is the defaults. One cohort in one patient-level study; at most 300,000 patients × runs per request (refused, naming the limit). Seeded: the same document gives the same model | | \`labs\` | \`test, anchor:{eventType,key?,value?,occurrence?}, eventType?:"Lab_Test", testKey?:"TEST", resultKey?:"RESULT", baselineDays?:90, fromWeeks?:-12, toWeeks?:52, binWeeks?:4, responseDrop?:50, cohorts:[names] (≤4)\` | **a lab value's course from a treatment start**: each patient's \`test\` results (timeline events of \`eventType\` whose \`testKey\` names the test, value in \`resultKey\`) as % change from their own baseline — the last value in the \`baselineDays\` up to and including the anchor; a non-positive baseline is left out — per \`binWeeks\` bin from \`fromWeeks\` to \`toWeeks\`, drawn as median and interquartile band (from 2-point histograms, exact to within a point; a patient counts once per bin, at their mean; the baseline reading itself is not drawn). Best response is each patient's lowest value after the anchor; the share falling by at least \`responseDrop\`% is counted exactly (50 on PSA is PCWG3's PSA50), and two cohorts get a Fisher test on it. A **waterfall** of every patient's best change is drawn only from participants entitled to patient-level values (a method only PatientLevelResolver has, like a plot's points); an aggregate-only participant's patients stay in the histograms and the rate and are counted as withheld from the bars. With no evaluable patient, the refusal lists the tests the study records | | \`cancertypes\` | \`genes:[≤30], groupBy?:"CANCER_TYPE"\|\|"study", measure?:"any"\|"driver", minSamples?:10, show?:"frequency"\|"count", unit?:"patient"\|"sample", cohorts:[names] (≤2)\` | **Cancer Types Summary**, cBioPortal's results-view tab: per group of the cohort's PATIENTS (the portal's default; \`unit:"sample"\` counts samples — a patient counts in each group one of their samples is in) — a sample attribute's value (folded on case and whitespace like every clinical predicate; samples with no value form their own group) or the study — the share altered in at least one of \`genes\`, split by kind: only mutations, only amplifications, only deep deletions, only structural variants, or more than one kind. A group's denominator is its samples profiled for at least one of the genes (the union of their panels, as an oncoprint's merged row counts it); its total sample count is shown beside it. \`measure:"driver"\` counts OncoKB drivers only. Groups with fewer than \`minSamples\` profiled samples are pooled into one bar, named and counted, never dropped. Counts only, so every local participant can answer; a remote peer cannot yet and is named | | \`overlap\` | \`cohorts:[names] (2–6)\` | **how cohorts' patients intersect**, cBioPortal's group-comparison Overlap tab: for every combination of the cohorts, the patients — and their samples — in exactly those cohorts and none of the others, with each cohort's own size. Two or three cohorts are drawn as the portal draws them, area-proportional diagrams of samples and of patients side by side; more as an UpSet plot. A cohort is a set of patients, so a patient's samples go wherever the patient does. Each study's database answers with one count per combination (a bitmask of membership) — counts only — and a cohort's size here must equal its own count, or the view says the definition gave two answers. Ids are study-qualified, so patients only ever overlap within a study | @@ -1877,8 +1879,9 @@ export const AGENT_GRAMMAR = `Document shape: "baselineDays":90,"fromWeeks":-12,"toWeeks":52,"binWeeks":4,"responseDrop":50,"cohorts":[names?(max 4)]}, // % change from baseline over time (median/IQR), best response, responders (PSA50), waterfall {"kind":"riskmodel","cohorts":[one name],"endpoint":"OS","measure":"any|driver","strata":""?,"covariates":[""...(max 6)]?, - "anchor":{"eventType":"Treatment","key":"SUBTYPE","value":"Immuno"}?,"validate":true?,"alpha":0.5,"runs":10,"groups":3 | "cuts":[0.8],"minMutated":30,"seed":1,"labelAbove":0.5?}, - // OncoCast-style elastic-net Cox ensemble on mutated genes: held-out C-index, genes kept (and their volcano), KM by risk group, score on 0–10 (patient-level, one study) + "anchor":{"eventType":"Treatment","key":"SUBTYPE","value":"Immuno"}?,"validate":true?,"alpha":0.5,"runs":10,"groups":3 | "cuts":[0.8],"minMutated":30,"seed":1,"labelAbove":0.5?, + "describe":{"attributes":[""...(max 20)]?,"genes":[""...(max 20)]?}?}, + // OncoCast-style elastic-net Cox ensemble on mutated genes: held-out C-index, genes kept (and their volcano), KM by risk group, score on 0–10, and with describe each group's makeup (Table 1; model rows marked, no q) (patient-level, one study) {"kind":"cancertypes","genes":[...(max 30)],"groupBy":"CANCER_TYPE||study","measure":"any|driver","minSamples":10,"show":"frequency|count","unit":"patient|sample"}, // the portal's Cancer Types Summary: % of samples altered per cancer type, split by mutation/amp/deep del/SV/multiple {"kind":"overlap","cohorts":[names?(2-6)]}, diff --git a/lib/federation-views.ts b/lib/federation-views.ts index 84f95a8..1e0a025 100644 --- a/lib/federation-views.ts +++ b/lib/federation-views.ts @@ -949,6 +949,10 @@ export interface RiskmodelRowsOpts { driverPin?: string | null; /** re-anchor the clock at a timeline event, with the delayed entry that implies */ anchor?: { eventType: string; key?: string; value?: string; occurrence?: "FIRST" | "LAST" } | null; + /** Table 1's clinical attributes, returned as d0, d1, … (the folded value: the patient's, or the + * first sample's), r0, r1, … (that value as spelled) and u0, u1, … (how many distinct folded + * values the patient's samples hold) */ + describe?: string[]; } /** @@ -962,9 +966,9 @@ export interface RiskmodelRowsOpts { * patients needs their rows. */ export function riskmodelRowsSql(cohortSql: string, study: string, db: string, o: RiskmodelRowsOpts): string { - const { endpoint, columns, strata = null, covariates = [], driverPin = null, anchor = null } = o; + const { endpoint, columns, strata = null, covariates = [], driverPin = null, anchor = null, describe = [] } = o; if (!/^[A-Z][A-Z0-9_]{0,39}$/.test(endpoint)) throw new Error("bad survival endpoint '" + String(endpoint).slice(0, 40) + "'"); - for (const a of [strata, ...covariates]) if (a && !/^[A-Za-z0-9_:/.-]{1,128}$/.test(a)) throw new Error("bad attribute name '" + String(a).slice(0, 40) + "'"); + for (const a of [strata, ...covariates, ...describe]) if (a && !/^[A-Za-z0-9_:/.-]{1,128}$/.test(a)) throw new Error("bad attribute name '" + String(a).slice(0, 40) + "'"); const months = columns?.months || (endpoint + "_MONTHS"), status = columns?.status || (endpoint + "_STATUS"); const pid = pidOf("patient_unique_id", study); const attrs = [...(strata ? [strata] : []), ...covariates]; @@ -1022,16 +1026,31 @@ export function riskmodelRowsSql(cohortSql: string, study: string, db: string, o " AND " + foldedEqualsSql("ce.event_type", anchor.eventType) + "\n" + anchorAttr + " AND p.stable_id IN (SELECT pid FROM cohort)\n" + " GROUP BY pid\n" + + ")" : "") + + (describe.length ? + ",\ndescr AS (\n" + + " SELECT " + pid + " AS pid,\n" + + describe.map((a, k) => { + const is = "attribute_name = " + sq(a), fold = clinicalValueFoldSql("attribute_value"); + return " argMinIf(" + fold + ", sample_unique_id, " + is + ") AS d" + k + ", argMinIf(attribute_value, sample_unique_id, " + is + ") AS r" + k + + ", uniqExactIf(" + fold + ", " + is + ") AS u" + k; + }).join(",\n") + "\n" + + " FROM " + clinicalTable(db, describe) + "\n" + + " WHERE cancer_study_identifier = " + sq(study) + " AND attribute_name IN (" + describe.map(sq).join(",") + ") AND trimBoth(attribute_value) != ''\n" + + " AND " + pid + " IN (SELECT pid FROM cohort)\n" + + " GROUP BY pid\n" + ")" : "") + "\n" + "SELECT os.months AS m, if(os.status LIKE '1%', 1, 0) AS ev, os.half AS half, panels.panels AS panels, mut.genes AS genes" + (strata ? ", attrs.a0 AS stratum" : "") + covariates.map((_, k) => ", attrs.a" + (k + (strata ? 1 : 0)) + " AS c" + k).join("") + - (anchor ? ", anch.hit AS anchored, round(anch.day / 30.4375, 4) AS anchor_months" : "") + "\n" + + (anchor ? ", anch.hit AS anchored, round(anch.day / 30.4375, 4) AS anchor_months" : "") + + describe.map((_, k) => ", descr.d" + k + " AS d" + k + ", descr.r" + k + " AS r" + k + ", descr.u" + k + " AS u" + k).join("") + "\n" + "FROM os\n" + " LEFT JOIN panels ON panels.pid = os.pid\n" + " LEFT JOIN mut ON mut.pid = os.pid\n" + (attrs.length ? " LEFT JOIN attrs ON attrs.pid = os.pid\n" : "") + (anchor ? " LEFT JOIN anch ON anch.pid = os.pid\n" : "") + + (describe.length ? " LEFT JOIN descr ON descr.pid = os.pid\n" : "") + "WHERE os.months IS NOT NULL AND os.months > 0\n" + // a fixed order: the ensemble's seeded splits are drawn over row positions, so the same // document gives the same model only if the rows come back in the same order. The id orders diff --git a/lib/method.ts b/lib/method.ts index de3ce0f..43b09d6 100644 --- a/lib/method.ts +++ b/lib/method.ts @@ -16,6 +16,7 @@ import { OQL_UNALTERED } from "./plot-oql.ts"; import { BUILD_COMMIT } from "./build-info.ts"; import { COX_MAX_ITERATIONS, LOGISTIC_MAX_ITERATIONS } from "./stats.ts"; import { SYSTEMIC_LINE_CLASSES, CLASS_KEYS } from "./treatment-lines.ts"; +import { table1RowKey, TABLE1_CONSIDERED, TABLE1_DEFAULT_ATTRIBUTES, TABLE1_DEFAULT_MAX_VALUES, TABLE1_MAX_VALUES, type Table1Histogram } from "./views/riskmodel-table1.ts"; export interface MethodStep { kind: "fetch" | "transform" | "statistic"; @@ -2189,6 +2190,8 @@ export interface RiskmodelTrace { study: string; cohort: string; cohortPatients: number; rowsReturned: number; r: import("./types.ts").RiskmodelResult; lambdas: number; folds: number; testShare: number; + /** Table 1's numeric rows: each group's value distribution, for the script */ + table1Histograms: Table1Histogram[]; } /** @@ -2242,8 +2245,110 @@ export function riskmodelMethod(t: RiskmodelTrace): MethodRecord { ...r.groups.map((g) => [g.name, g.curve.n.toLocaleString() + " patients, " + g.curve.events.toLocaleString() + " events, median " + (g.curve.median == null ? "not reached" : g.curve.median.toFixed(1) + " months") + (g.hr != null ? ", HR " + g.hr.toFixed(2) : "")] as [string, string])] }, ]; - return { steps, implementation: { functions: ["coxnetEnsemble", "coxnetPath", "cvLambda", "concordance", "coxFit"], file: "lib/coxnet.ts", commit: BUILD_COMMIT }, - verification: "lib/coxnet.ts's path is checked against scikit-survival's CoxnetSurvivalAnalysis, its delayed-entry, stratified fits (with standard errors) against statsmodels' PHReg and its C-index against lifelines, by test/coxnet-golden.test.ts (fixtures from python/golden/coxnet.py)." }; + const fitVerification = "lib/coxnet.ts's path is checked against scikit-survival's CoxnetSurvivalAnalysis, its delayed-entry, stratified fits (with standard errors) against statsmodels' PHReg and its C-index against lifelines, by test/coxnet-golden.test.ts (fixtures from python/golden/coxnet.py)."; + const implementation = { functions: ["coxnetEnsemble", "coxnetPath", "cvLambda", "concordance", "coxFit"], file: "lib/coxnet.ts", commit: BUILD_COMMIT }; + if (!r.table1) return { steps, implementation, verification: fitVerification }; + const t1 = table1Record(r.table1, t.table1Histograms); + return { steps: [...steps, ...t1.steps], + implementation: { ...implementation, functions: [...implementation.functions, "buildTable1", "chiSquareTest", "kruskalWallis", "quantilesFromHistogram", "benjaminiHochberg"], + file: implementation.file + ", lib/views/riskmodel-table1.ts, lib/stats.ts" }, + verification: fitVerification + " The script below recomputes Table 1's tests and q-values from its counts (the model itself it does not refit); χ² matches scipy " + VERIFIED.chi2.version + + " chi2_contingency, Kruskal–Wallis and quantiles scipy/numpy, Benjamini–Hochberg statsmodels multipletests (test/tables-golden.test.ts, test/moments-golden.test.ts).", + python: withResult(t1.python), check: t1.check }; +} + +/** Table 1's steps, and the script that recomputes its tests from its counts. */ +function table1Record(t: import("./types.ts").RiskTable1, histograms: Table1Histogram[]) { + const fmtP = (p: number | null) => p == null ? "—" : pFmt(p); + const shown = t.rows.filter((x) => x.kind !== "skipped"); + const groups = t.groups.map((g) => g.name); + const steps: MethodStep[] = [ + { kind: "fetch", title: "Table 1: the groups' patients again, with their clinical values", + detail: "The riskmodel:table1 query above is the riskmodel:rows query with Table 1's attributes added: for each, the patient's value (a sample-level attribute: their first sample's), " + + "case and whitespace folded as everywhere in this app and shown in its commonest spelling, and how many values their samples hold. Its rows are kept exactly as the model's were and must match " + + "the model's patients one by one — survival time, entry and event — or the view is refused, so each group's column describes the patients its curve drew" + + (t.reportedOn === "validation half" ? " (the validation half)" : "") + ". riskmodel:table1-attributes gives each attribute's type (NUMBER is numeric), riskmodel:panels each mutation panel's genes. " + + "Patient-level: the rows never leave this server's engine; only the counts below do.", + params: t.groups.map((g) => [g.name, g.n.toLocaleString("en-US") + " patients"] as [string, string]) }, + { kind: "transform", title: "Table 1's rows", + detail: (t.defaults.attributes + ? "Attributes by the default rule: the model's stratum and covariates, then the study's attributes in its priority order, leaving out survival endpoints' columns (…_MONTHS, …_STATUS), " + + "considering the first " + TABLE1_CONSIDERED + "; of those, an attribute is left out if no patient here has a value, if it takes one value among them, or if it is categorical with more than " + + TABLE1_DEFAULT_MAX_VALUES + " values (free text, ids); the first " + TABLE1_DEFAULT_ATTRIBUTES + " that remain are shown." + : "Attributes as the view names them; a categorical one with more than " + TABLE1_MAX_VALUES + " values among these patients is set aside.") + " " + + (t.defaults.genes ? "Genes: the " + t.rows.filter((x) => x.kind === "gene").length + " the model kept in the most runs." : "Genes as the view names them.") + " " + + "A gene is counted as mutated the way the model counts it (" + (t.measure === "driver" ? "OncoKB drivers" : "any mutation") + "), and only among the patients whose mutation panels cover it " + + "(a whole-exome patient covers every gene): a patient whose panels do not is not profiled for it, never unaltered, and is counted as such in each group. " + + "A row the groups differ in by construction is marked: a gene the model kept in any run and a covariate are \"in the model\" (the score is built from them), the endpoint's own columns are \"outcome\" (the groups are cut to differ in survival); the model's stratum is marked \"stratum\" but is a description, since the score has no stratum term.", + params: [...t.omitted.map((o) => ["left out: " + o.name, o.why] as [string, string]), + ...t.rows.filter((x) => x.kind === "skipped").map((x) => ["set aside: " + x.name, (x as { skipped: string }).skipped] as [string, string])] }, + { kind: "statistic", title: "Each row tested across the groups", + detail: "A categorical attribute: Pearson χ² on the group × value table, no continuity correction, an all-zero row or column dropped first, patients with no value left out. " + + "A gene: Pearson χ² on the group × (mutated, not mutated) table of its profiled patients. A numeric attribute: Kruskal–Wallis on each group's values, ranks averaged within ties and tie-corrected; medians and quartiles interpolating linearly between ranks. " + + "The smallest expected count is reported, since χ² is poor below about 5.", + params: shown.map((x) => [table1RowKey(x).replace(/^attribute:/, ""), x.p == null ? "not tested: " + (x as { untested?: string }).untested + : (x.test === "chi-squared" ? "χ² " : "H ") + x.statistic + " on " + x.df + " df · p " + fmtP(x.p) + (x.role ? " · " + x.role : "")] as [string, string]) }, + { kind: "statistic", title: "Benjamini–Hochberg across Table 1's open rows", + detail: "q-values over the " + t.tested + " tested rows that are not in the model and not its outcome. Those rows differ between the groups by construction, so they get a p-value but no q-value and do not enter the correction, where their p-values would shrink every other row's q." }, + ]; + const tested = shown.filter((x) => x.p != null); + const cat = tested.filter((x) => x.kind === "cat") as Extract[]; + const genes = tested.filter((x) => x.kind === "gene") as Extract[]; + const numKeys = new Set(tested.filter((x) => x.kind === "num").map(table1RowKey)); + const hist = histograms.filter((h) => numKeys.has(h.row)); + const family = tested.filter((x) => x.q != null).map(table1RowKey); + const python = [ + "# Recompute Table 1's tests from its counts: scipy " + VERIFIED.chi2.version + " and statsmodels " + VERIFIED.bh.version + " are the versions our code is tested against.", + "import io", + "import numpy as np", + "import pandas as pd", + "from scipy.stats import chi2_contingency, kruskal", + "from statsmodels.stats.multitest import multipletests", + "", + "groups = " + JSON.stringify(groups), + RESULT_INIT, + "", + "def chi2_p(table):", + " table = table[table.sum(axis=1) > 0][:, table.sum(axis=0) > 0] # drop an all-zero row or column", + " return chi2_contingency(table, correction=False).pvalue", + ...(cat.length ? ["", + ...pyTable("categorical", "Categorical rows: patients per group (by index into groups) and value (by index into the row's values)", ["row", "group", "level", "n"], + cat.flatMap((x) => x.levels.flatMap((l, li) => l.counts.map((n, gi) => [table1RowKey(x), gi, li, n]))), "table1_categorical.csv"), + "for row, g in categorical.groupby(\"row\", sort=False):", + " table = g.pivot_table(index=\"group\", columns=\"level\", values=\"n\", aggfunc=\"sum\", fill_value=0).to_numpy()", + " RESULT[f\"table1:{row}:p\"] = chi2_p(table)", + " print(f\"{row}: chi-squared p={RESULT[f'table1:{row}:p']:.4g}\")"] : []), + ...(genes.length ? ["", + ...pyTable("genes", "Gene rows: per group (by index into groups), patients with the gene mutated and patients whose panels cover it", ["row", "group", "altered", "profiled"], + genes.flatMap((x) => x.perGroup.map((c, gi) => [table1RowKey(x), gi, c.altered, c.profiled])), "table1_genes.csv"), + "for row, g in genes.groupby(\"row\", sort=False):", + " table = np.column_stack([g[\"altered\"], g[\"profiled\"] - g[\"altered\"]])", + " RESULT[f\"table1:{row}:p\"] = chi2_p(table)", + " print(f\"{row}: \" + \", \".join(f\"{a}/{p}\" for a, p in zip(g[\"altered\"], g[\"profiled\"])) + f\" mutated; chi-squared p={RESULT[f'table1:{row}:p']:.4g}\")"] : []), + ...(hist.length ? ["", + ...pyTable("numeric", "Numeric rows: each group's values, as (value, patients)", ["row", "group", "value", "n"], hist.map((h) => [h.row, h.group, h.value, h.n]), "table1_numeric.csv"), + "for row, g in numeric.groupby(\"row\", sort=False):", + " values = {grp: np.repeat(h[\"value\"].to_numpy(), h[\"n\"].to_numpy()) for grp, h in g.groupby(\"group\", sort=False)}", + " for grp, v in values.items():", + " for name, q in ((\"median\", 50), (\"q1\", 25), (\"q3\", 75)):", + " RESULT[f\"table1:{row}:{grp}:{name}\"] = np.percentile(v, q)", + " RESULT[f\"table1:{row}:p\"] = kruskal(*values.values()).pvalue", + " print(f\"{row}: Kruskal-Wallis p={RESULT[f'table1:{row}:p']:.4g}\")"] : []), + "", + "# q-values over the open rows only: those in the model or its outcome differ by construction", + "family = " + JSON.stringify(family), + "if family:", + " for row, qv in zip(family, multipletests([RESULT[f\"table1:{r}:p\"] for r in family], method=\"fdr_bh\")[1]):", + " RESULT[f\"table1:{row}:q\"] = qv", + ].join("\n") + "\n"; + const medians = tested.filter((x) => x.kind === "num").flatMap((x) => (x as Extract).perGroup + .map((c, gi) => ({ key: table1RowKey(x) + ":" + groups[gi], ...c })).filter((c) => c.n > 0)); + const check = { + ...checks("table1:", tested.map((x) => ({ key: table1RowKey(x), p: x.p, ...(x.q != null ? { q: x.q } : {}) })), (x) => x.key, [["p", "p"], ["q", "q"]]), + // medians and quartiles are reported to 3 decimals + ...checks("table1:", medians, (m) => m.key, [["median", "median", { abs: 5e-4 + 1e-9 }], ["q1", "q1", { abs: 5e-4 + 1e-9 }], ["q3", "q3", { abs: 5e-4 + 1e-9 }]]), + }; + return { steps, python, check }; } export function timelineMethod(t: { cohorts: string[]; queryText: string; sortName: string | null; desc: boolean; shown: number; diff --git a/lib/types.ts b/lib/types.ts index 141d5b7..ba95ebd 100644 --- a/lib/types.ts +++ b/lib/types.ts @@ -963,6 +963,54 @@ export interface RiskmodelView { /** the volcano names the genes kept in at least this share of runs (0 to 1, default 0.5); drawn in * the browser from the whole result, so changing it re-runs nothing */ labelAbove?: number; + /** Table 1: the risk groups' clinical and genomic makeup, tested across groups. Present, the view + * computes it; a list left out is the default (the study's leading attributes; the genes the + * model kept most often) */ + describe?: { attributes?: string[]; genes?: string[] }; +} + +/** Why a Table 1 row differs between the risk groups by construction, or is the model's stratum. */ +export type RiskTable1Role = "in the model" | "outcome" | "stratum" | null; + +interface RiskTable1Test { + /** null when the row cannot be tested (the reason in `untested`) */ + test: "chi-squared" | "Kruskal-Wallis" | null; statistic: number | null; df: number | null; p: number | null; + /** Benjamini–Hochberg across the tested rows that are not in the model or the outcome; null for those */ + q: number | null; + untested?: string; +} + +export type RiskTable1Row = + | ({ kind: "cat"; name: string; per: "patient" | "sample"; role: RiskTable1Role; + /** every value recorded among the groups' patients, commonest first, with its count per group */ + levels: { value: string; counts: number[] }[]; + /** per group, patients with no value recorded */ + missing: number[]; + /** patients whose samples hold more than one value (counted under their first sample's) */ + multiValued: number; minExpected: number | null } & RiskTable1Test) + | ({ kind: "num"; name: string; per: "patient" | "sample"; role: RiskTable1Role; + perGroup: { n: number; median: number | null; q1: number | null; q3: number | null; missing: number }[]; + multiValued: number } & RiskTable1Test) + | ({ kind: "gene"; name: string; role: RiskTable1Role; + /** per group: patients with the gene mutated, patients whose panels cover it, and those whose do not */ + perGroup: { altered: number; profiled: number; notProfiled: number }[]; minExpected: number | null } & RiskTable1Test) + | { kind: "skipped"; name: string; what: "attribute" | "gene"; skipped: string }; + +export interface RiskTable1 { + /** the patients described: the same as the groups' curves */ + reportedOn: "held-out scores" | "validation half"; + groups: { name: string; n: number }[]; + /** the gene rows count mutations as the model does */ + measure: "any" | "driver"; + /** what the table holds, as a document would name it — what "Add to document" writes */ + describe: { attributes: string[]; genes: string[] }; + /** the default lists were used (no list named in the view) */ + defaults: { attributes: boolean; genes: boolean }; + rows: RiskTable1Row[]; + /** attributes the default rule considered and left out, with why */ + omitted: { name: string; why: string }[]; + /** rows in the Benjamini–Hochberg family */ + tested: number; } export interface RiskmodelResult { @@ -1009,6 +1057,8 @@ export interface RiskmodelResult { /** the held-out score (raw) at each whole percentile 0…100 of the scores the cuts are placed on, * the rank a cut at that quantile starts from — what the risk-groups dialog places its handles by */ percentiles: number[]; + /** with `describe`: Table 1, the groups' makeup */ + table1?: RiskTable1; } export interface LabsView { diff --git a/lib/verification.ts b/lib/verification.ts index b5510a5..f4b984c 100644 --- a/lib/verification.ts +++ b/lib/verification.ts @@ -35,7 +35,7 @@ export const STAT_CHECKS: StatCheck[] = [ { method: "Logistic regression", usedFor: "logistic view, odds ratios", keys: ["logistic"], agreement: "coefficients, standard errors and p-values to 6 significant digits", test: "test/regression-golden.test.ts", fixture: "regression", caseKeys: ["logistic"] }, - { method: "Fisher's exact test, χ² test", usedFor: "gene comparisons, co-occurrence, clinical tables", keys: ["fisher", "chi2"], + { method: "Fisher's exact test, χ² test", usedFor: "gene comparisons, co-occurrence, clinical tables, the risk model's Table 1", keys: ["fisher", "chi2"], agreement: "p-values to 6 significant digits, including tables with zeros and very large counts", test: "test/tables-golden.test.ts", fixture: "tables", caseKeys: ["fisher", "chi2"] }, { method: "Benjamini–Hochberg false discovery rate", usedFor: "every q-value", keys: ["bh"], @@ -44,7 +44,7 @@ export const STAT_CHECKS: StatCheck[] = [ { method: "Cochran–Mantel–Haenszel test, McNemar's test", usedFor: "stratified and paired comparisons", keys: ["cmh", "mcnemar"], agreement: "p-values and odds ratios to 6 significant digits", test: "test/tables-golden.test.ts", fixture: "tables", caseKeys: ["cmh", "mcnemar"] }, - { method: "Welch's t-test, Welch's ANOVA, Mann–Whitney, Kruskal–Wallis, medians", usedFor: "comparing a clinical value between groups; each part of a composition between cohorts", keys: ["welch", "anova", "ranks"], + { method: "Welch's t-test, Welch's ANOVA, Mann–Whitney, Kruskal–Wallis, medians", usedFor: "comparing a clinical value between groups (the risk model's Table 1 among them); each part of a composition between cohorts", keys: ["welch", "anova", "ranks"], agreement: "p-values to 6 significant digits", test: "test/moments-golden.test.ts", fixture: "moments", caseKeys: ["welch", "anova", "ranks"] }, { method: "Pearson correlation and its meta-analysis across studies", usedFor: "co-expression, scatter plots", keys: ["pearson", "meta"], diff --git a/lib/views/riskmodel-table1.ts b/lib/views/riskmodel-table1.ts new file mode 100644 index 0000000..7c569d6 --- /dev/null +++ b/lib/views/riskmodel-table1.ts @@ -0,0 +1,196 @@ +// Table 1 for the risk model: each risk group's clinical and genomic makeup, tested across groups. +// +// Pure: lib/views/riskmodel.ts fetches the rows (the fit's own query, with Table 1's attributes +// added), matches them to the fit's patients and hands them here by position, without identifiers; +// only counts come out. Gene rates are panel-aware (docs/PERFORMANCE.md): a group's rate for a gene +// is over its patients whose panels cover that gene, and a patient whose panels do not is "not +// profiled", never unaltered. Rows the groups differ in by construction — the model's genes and +// covariates, its outcome — are marked and kept out of the q-value family. +import { benjaminiHochberg, chiSquareTest, kruskalWallis, quantilesFromHistogram, type Histogram } from "../stats.ts"; +import type { RiskTable1, RiskTable1Row, RiskTable1Role } from "../types.ts"; + +export const TABLE1_MAX_ATTRIBUTES = 20, TABLE1_MAX_GENES = 20; +/** the default attributes: at most this many, chosen from the first TABLE1_CONSIDERED in the study's priority order */ +export const TABLE1_DEFAULT_ATTRIBUTES = 8; +export const TABLE1_CONSIDERED = 20; +/** a default attribute takes at most this many values among the patients (free text and ids take more) */ +export const TABLE1_DEFAULT_MAX_VALUES = 12; +/** past this many values a categorical attribute's chi-squared table is too sparse to say anything */ +export const TABLE1_MAX_VALUES = 100; +/** a survival endpoint's columns: never a default, since the groups are cut to differ in survival */ +export const OUTCOME_COLUMN = /_(MONTHS|STATUS)$/; + +/** A clinical attribute as the study's metadata describes it. */ +export interface Table1AttrMeta { id: string; numeric: boolean; patient: boolean } + +/** One modelled patient, by position in the fit: their mutation panels and mutated genes, and per Table 1 attribute its folded value, its spelling and how many values their samples hold. */ +export interface Table1Patient { panels: string[]; genes: string[]; d: string[]; r: string[]; u: number[] } + +export interface Table1Input { + study: string; + reportedOn: RiskTable1["reportedOn"]; measure: RiskTable1["measure"]; + /** each group's name and its patients, as positions in `patients` */ + groups: { name: string; members: number[] }[]; + /** every modelled patient, by the fit's position */ + patients: Table1Patient[]; + /** the attributes fetched, in the order of each patient's d/r/u; `meta` null when the study has no such attribute */ + attributes: { name: string; meta: Table1AttrMeta | null }[]; + /** the attributes came from the default rule (filtered and capped here), rather than being named */ + attributesDefault: boolean; + genes: string[]; genesDefault: boolean; + /** the genes each mutation panel covers (WES, which covers every gene, is not listed) */ + panelGenes: Map>; + model: { genes: string[]; covariates: string[]; strata: string | null; outcome: string[] }; +} + +/** A numeric row's value distribution in one group, for the method record's script. */ +export interface Table1Histogram { row: string; group: string; value: number; n: number } + +const roleOfAttribute = (name: string, m: Table1Input["model"]): RiskTable1Role => + m.covariates.includes(name) ? "in the model" : m.outcome.includes(name) ? "outcome" : name === m.strata ? "stratum" : null; + +/** The row's key in the method record's checks. */ +export const table1RowKey = (row: { kind: string; name: string }) => (row.kind === "gene" ? "gene:" : "attribute:") + row.name; + +const untested = (why: string) => ({ test: null, statistic: null, df: null, p: null, q: null, untested: why }); + +export function buildTable1(t: Table1Input): { table: RiskTable1; histograms: Table1Histogram[] } { + const reported = t.groups.flatMap((g) => g.members); + const rows: RiskTable1Row[] = []; + const omitted: RiskTable1["omitted"] = []; + const histograms: Table1Histogram[] = []; + let defaultsKept = 0; + + t.attributes.forEach(({ name, meta }, k) => { + const isDefault = t.attributesDefault; + if (isDefault && defaultsKept >= TABLE1_DEFAULT_ATTRIBUTES) return; + const skip = (why: string) => { if (isDefault) omitted.push({ name, why }); else rows.push({ kind: "skipped", name, what: "attribute", skipped: why }); }; + if (!meta) { skip("not a clinical attribute of " + t.study); return; } + const role = roleOfAttribute(name, t.model), per = meta.patient ? "patient" as const : "sample" as const; + const multiValued = reported.filter((i) => t.patients[i].u[k] > 1).length; + if (meta.numeric) { + const hists = t.groups.map((g) => { + const counts = new Map(); + for (const i of g.members) { + const d = t.patients[i].d[k]; + const v = d === "" ? NaN : Number(d); + if (Number.isFinite(v)) counts.set(v, (counts.get(v) || 0) + 1); + } + return [...counts].map(([v, n]) => ({ v, n })).sort((a, b) => a.v - b.v) as Histogram[]; + }); + const distinct = new Set(hists.flatMap((h) => h.map((x) => x.v))).size; + if (!distinct) { skip("no value recorded for these patients"); return; } + if (isDefault && distinct < 2) { skip("takes one value among these patients"); return; } + if (isDefault) defaultsKept++; + const perGroup = t.groups.map((g, gi) => { + const n = hists[gi].reduce((s, x) => s + x.n, 0); + const [median, q1, q3] = quantilesFromHistogram(hists[gi], [0.5, 0.25, 0.75]); + return { n, median, q1, q3, missing: g.members.length - n }; + }); + t.groups.forEach((g, gi) => { for (const h of hists[gi]) histograms.push({ row: table1RowKey({ kind: "num", name }), group: g.name, value: h.v, n: h.n }); }); + const filled = perGroup.filter((p) => p.n > 0).length; + const test = distinct < 2 ? untested("takes one value among these patients") + : filled < 2 ? untested("fewer than two groups have a value") + : (() => { const kw = kruskalWallis(hists); return { test: "Kruskal-Wallis" as const, statistic: kw.h, df: kw.df, p: kw.p, q: null }; })(); + rows.push({ kind: "num", name, per, role, perGroup, multiValued, ...test }); + return; + } + // categories, by the folded value; each shown as its commonest spelling + const spelled = new Map>(); + const byGroup = t.groups.map((g) => { + const counts = new Map(); + for (const i of g.members) { + const p = t.patients[i], d = p.d[k]; + if (d === "") continue; + counts.set(d, (counts.get(d) || 0) + 1); + let sp = spelled.get(d); if (!sp) spelled.set(d, sp = new Map()); + sp.set(p.r[k], (sp.get(p.r[k]) || 0) + 1); + } + return counts; + }); + const keys = [...spelled.keys()]; + if (!keys.length) { skip("no value recorded for these patients"); return; } + if (isDefault && keys.length < 2) { skip("takes one value among these patients"); return; } + const limit = isDefault ? TABLE1_DEFAULT_MAX_VALUES : TABLE1_MAX_VALUES; + if (keys.length > limit) { skip("takes " + keys.length + " values here (" + (isDefault ? "a default takes at most " + limit : "Table 1 describes at most " + limit) + ")"); return; } + if (isDefault) defaultsKept++; + const label = (key: string) => [...spelled.get(key)!].sort((a, b) => b[1] - a[1] || a[0].localeCompare(b[0]))[0][0]; + const levels = keys.map((key) => ({ value: label(key), counts: byGroup.map((c) => c.get(key) || 0) })) + .sort((a, b) => b.counts.reduce((s, x) => s + x, 0) - a.counts.reduce((s, x) => s + x, 0) || a.value.localeCompare(b.value)); + const missing = t.groups.map((g, gi) => g.members.length - levels.reduce((s, l) => s + l.counts[gi], 0)); + const filled = byGroup.filter((c) => c.size > 0).length; + let test, minExpected: number | null = null; + if (keys.length < 2) test = untested("takes one value among these patients"); + else if (filled < 2) test = untested("fewer than two groups have a value"); + else { + const x = chiSquareTest(t.groups.map((_, gi) => levels.map((l) => l.counts[gi]))); + test = { test: "chi-squared" as const, statistic: x.chi2, df: x.df, p: x.p, q: null }; + minExpected = x.minExpected; + } + rows.push({ kind: "cat", name, per, role, levels, missing, multiValued, minExpected, ...test }); + }); + + // genes: a symbol resolves against the panels' genes and the genes mutated in these patients + const known = new Set(); + for (const s of t.panelGenes.values()) for (const g of s) known.add(g); + for (const p of t.patients) for (const g of p.genes) known.add(g); + const byUpper = new Map(); + for (const g of known) if (!byUpper.has(g.toUpperCase())) byUpper.set(g.toUpperCase(), g); + for (const asked of t.genes) { + const gene = known.has(asked) ? asked : byUpper.get(asked.toUpperCase()); + if (!gene) { rows.push({ kind: "skipped", name: asked, what: "gene", skipped: "on none of the mutation panels these patients were sequenced on, and mutated in none of them" }); continue; } + // two spellings of one gene ("keap1", "KEAP1") are one row + if (rows.some((x) => x.kind === "gene" && x.name === gene)) continue; + const covers = (p: Table1Patient) => p.panels.some((pn) => pn === "WES" || !!t.panelGenes.get(pn)?.has(gene)); + const perGroup = t.groups.map((g) => { + let altered = 0, profiled = 0; + for (const i of g.members) { + const p = t.patients[i]; + if (!covers(p)) continue; + profiled++; + if (p.genes.includes(gene)) altered++; + } + return { altered, profiled, notProfiled: g.members.length - profiled }; + }); + const role: RiskTable1Role = t.model.genes.includes(gene) ? "in the model" : null; + let test, minExpected: number | null = null; + if (perGroup.filter((x) => x.profiled > 0).length < 2) test = untested("fewer than two groups have a patient sequenced on a panel covering it"); + else if (perGroup.every((x) => x.altered === 0) || perGroup.every((x) => x.altered === x.profiled)) test = untested(perGroup.every((x) => x.altered === 0) ? "mutated in none of these patients" : "mutated in every profiled patient"); + else { + const x = chiSquareTest(perGroup.map((c) => [c.altered, c.profiled - c.altered])); + test = { test: "chi-squared" as const, statistic: x.chi2, df: x.df, p: x.p, q: null }; + minExpected = x.minExpected; + } + rows.push({ kind: "gene", name: gene, role, perGroup, minExpected, ...test }); + } + + // Benjamini–Hochberg over the rows that are open questions; the model's own rows and its outcome differ by construction + const family = rows.filter((r): r is Exclude => r.kind !== "skipped" && r.p != null && r.role !== "in the model" && r.role !== "outcome"); + benjaminiHochberg(family.map((r) => r.p!)).forEach((q, i) => { family[i].q = q; }); + + const named = (what: "attribute" | "gene") => rows.filter((r) => r.kind === "skipped" ? r.what === what : (r.kind === "gene") === (what === "gene")).map((r) => r.name); + return { + table: { + reportedOn: t.reportedOn, measure: t.measure, + groups: t.groups.map((g) => ({ name: g.name, n: g.members.length })), + describe: { attributes: named("attribute"), genes: named("gene") }, + defaults: { attributes: t.attributesDefault, genes: t.genesDefault }, + rows, omitted, tested: family.length, + }, + histograms, + }; +} + +/** + * The attributes the default rule considers, in order: the model's stratum and covariates, then the + * study's attributes by priority, leaving out survival endpoints' columns — at most TABLE1_CONSIDERED. + * buildTable1 then drops those that cannot describe these patients and keeps TABLE1_DEFAULT_ATTRIBUTES. + */ +export function table1DefaultCandidates(metaInOrder: string[], model: { covariates: string[]; strata: string | null }): string[] { + const out: string[] = []; + for (const a of [...(model.strata ? [model.strata] : []), ...model.covariates, ...metaInOrder]) { + if (out.length >= TABLE1_CONSIDERED) break; + if (!out.includes(a) && !OUTCOME_COLUMN.test(a)) out.push(a); + } + return out; +} diff --git a/lib/views/riskmodel.ts b/lib/views/riskmodel.ts index ca98675..feae8d9 100644 --- a/lib/views/riskmodel.ts +++ b/lib/views/riskmodel.ts @@ -14,7 +14,9 @@ import { kaplanMeier, logRankK } from "../stats.ts"; import { riskmodelMethod } from "../method.ts"; import type { ViewRun } from "./run.ts"; import { overLimit, refuse } from "./shared.ts"; -import { RISK_LABEL_ABOVE, riskGroupNames } from "../figures/riskmodel.ts"; +import { RISK_LABEL_ABOVE, RISK_FIGURE_GENES, riskGroupNames } from "../figures/riskmodel.ts"; +import { tableAttributesSql } from "../clinical-table.ts"; +import { buildTable1, table1DefaultCandidates, TABLE1_MAX_ATTRIBUTES, TABLE1_MAX_GENES, type Table1Patient } from "./riskmodel-table1.ts"; const RISK_MAX_RUNS = 100; /** patients x runs past which one request would not finish inside the analysis route's time limit */ @@ -41,13 +43,13 @@ export interface RiskFitSpec { } /** - * The view's settings that decide only how its groups are cut and drawn, never the model: the one + * The view's settings that decide only how its groups are cut, drawn and described (Table 1), never the model: the one * list a setting must be on to be left out of the fit's cache key. Everything else in the view — * including any setting added later — is in the key, so a setting nobody thought about costs at * worst a needless refit, never a fit served for a different model. A setting goes here only if * riskFitSpec does not read it (test/riskmodel-groups.test.ts checks each one). */ -export const RISK_PRESENTATION_KEYS = ["groups", "cuts", "labelAbove", "note", "cohorts"] as const; +export const RISK_PRESENTATION_KEYS = ["groups", "cuts", "labelAbove", "note", "cohorts", "describe"] as const; /** The view without its presentation settings, keys sorted: the part of the fit's key the view gives. */ export function riskModelView(view: RiskmodelView): RiskmodelView { @@ -119,12 +121,12 @@ export interface RiskKeptRow { t: number; entry: number; ev: boolean; stratum: s /** * Who of the rows `riskmodelRowsSql` returned is modelled, and who is left out and why — every * exclusion counted. The fit and Table 1 both read their patients through this, so the same rows - * give the same patients in the same order. + * give the same patients in the same order. `at[k]` is kept patient k's position in `rows`. */ export function keepRiskRows(rows: Record[], s: Pick) { let notProfiled = 0, noStratum = 0, noCovariate = 0, noAnchor = 0, anchorAfterStart = 0, noFollowUp = 0; - const kept: RiskKeptRow[] = []; - rows.forEach((x) => { + const kept: RiskKeptRow[] = [], at: number[] = []; + rows.forEach((x, row) => { const panels = ((x.panels as unknown as string[]) || []).map(String); if (!panels.length) { notProfiled++; return; } const st = s.strata ? String(x.stratum ?? "") : ""; @@ -140,8 +142,9 @@ export function keepRiskRows(rows: Record[], s: Pick entry)) { noFollowUp++; return; } } kept.push({ t, entry, ev: Number(x.ev) === 1, stratum: st, genes: ((x.genes as unknown as string[]) || []).map(String), panels, cov, half: Number(x.half) }); + at.push(row); }); - return { kept, left: { notProfiled, noStratum, noCovariate, noAnchor, anchorAfterStart, noFollowUp } }; + return { kept, at, left: { notProfiled, noStratum, noCovariate, noAnchor, anchorAfterStart, noFollowUp } }; } /** Fetch the cohort's rows and fit the ensemble. Throws RiskRefusal when the data cannot support a model. */ @@ -207,12 +210,15 @@ export async function fitRiskmodel(r: RiskResolver, args: RiskFitArgs): Promise< }; } -/** The fit's two queries, through the patient-level resolver alone (test/boundary.test.ts sweeps their SQL). */ -function riskViewRows(r: RiskResolver, args: RiskFitArgs) { +/** + * The fit's two queries — and Table 1's, which are the same with its attributes added — through the + * patient-level resolver alone (test/boundary.test.ts sweeps their SQL). + */ +function riskViewRows(r: RiskResolver, args: RiskFitArgs, table1?: { describe: string[] }) { const s = riskFitSpec(args.view, args.compile.oncokbPin); return Promise.all([ - r.viewRows("riskmodel:rows", args.def, (c, st, d) => riskmodelRowsSql(c, st, d, { endpoint: s.endpoint, columns: s.columns ?? undefined, strata: s.strata, - covariates: s.covariates, driverPin: s.driverPin, anchor: s.anchor })), + r.viewRows(table1 ? "riskmodel:table1" : "riskmodel:rows", args.def, (c, st, d) => riskmodelRowsSql(c, st, d, { endpoint: s.endpoint, columns: s.columns ?? undefined, strata: s.strata, + covariates: s.covariates, driverPin: s.driverPin, anchor: s.anchor, ...(table1 ? { describe: table1.describe } : {}) })), r.viewMeta("riskmodel:panels", (st, d) => riskmodelPanelGenesSql(st, d)), ]); } @@ -330,6 +336,44 @@ export function groupRiskmodel(fit: RiskFit, qs: number[]): RiskGrouping { neverHeldOut: fit.pool.length - poolScored.length, distinct, groupsMade, members }; } +// ---- Table 1 ---------------------------------------------------------------------------------- + +/** + * Table 1, the groups' clinical and genomic makeup (counted and tested by lib/views/riskmodel-table1.ts). + * Reads the fit's own rows again with Table 1's attributes added and keeps them as the fit did, and + * refuses unless they are the fit's patients in the fit's order — every survival time, entry and + * event the same — so each group's members are exactly the patients its curve drew. + */ +export async function riskTable1(r: RiskResolver, args: RiskFitArgs, fit: RiskFit, g: RiskGrouping, describe: NonNullable) { + const s = riskFitSpec(args.view, args.compile.oncokbPin), mark = r.queries.length; + const logged = () => r.queries.slice(mark).map((q) => ({ label: q.label, sql: q.sql, ...(q.ms != null ? { ms: q.ms } : {}) })); + const metaRows = await r.viewMeta("riskmodel:table1-attributes", (st, d) => tableAttributesSql(st, d)); + const meta = new Map(metaRows.map((m) => [String(m.id), { id: String(m.id), numeric: String(m.datatype).toUpperCase() === "NUMBER", patient: Number(m.patient) === 1 }] as const)); + const attrs = describe.attributes ? [...new Set(describe.attributes)] + : table1DefaultCandidates(metaRows.map((m) => String(m.id)), { covariates: s.covariates, strata: s.strata }); + const fetched = attrs.filter((a) => meta.has(a)); + const [rows, panelRows] = await riskViewRows(r, args, { describe: fetched }); + const { kept, at } = keepRiskRows(rows, s); + const same = kept.length === fit.n && kept.every((k, i) => k.t === fit.subjects.time[i] && k.entry === fit.subjects.entry[i] && (k.ev ? 1 : 0) === fit.subjects.event[i]); + if (!same) throw new RiskRefusal("Table 1 could not be matched to the model's patients: the cohort's rows changed after the model was fitted (" + kept.length + " patients now, " + fit.n + " then); run the view again", logged()); + const col = new Map(fetched.map((a, k) => [a, k] as const)); + const cell = (x: Record, p: string, a: string) => col.has(a) ? x[p + col.get(a)] : undefined; + const patients: Table1Patient[] = kept.map((k, i) => { + const x = rows[at[i]]; + return { panels: k.panels, genes: k.genes, d: attrs.map((a) => String(cell(x, "d", a) ?? "")), r: attrs.map((a) => String(cell(x, "r", a) ?? "")), u: attrs.map((a) => Number(cell(x, "u", a) ?? 0)) }; + }); + const months = s.columns?.months ?? s.endpoint + "_MONTHS", status = s.columns?.status ?? s.endpoint + "_STATUS"; + const { table, histograms } = buildTable1({ + study: args.study, reportedOn: fit.validated ? "validation half" : "held-out scores", measure: s.measure, + groups: g.groups.map((x, gi) => ({ name: x.name, members: g.members[gi] })), patients, + attributes: attrs.map((name) => ({ name, meta: meta.get(name) ?? null })), attributesDefault: !describe.attributes, + genes: describe.genes ? [...new Set(describe.genes)] : fit.features.slice(0, RISK_FIGURE_GENES).map((f) => f.gene), genesDefault: !describe.genes, + panelGenes: new Map(panelRows.map((x) => [String(x.panel), new Set((x.genes as unknown as string[]) || [])] as const)), + model: { genes: fit.features.map((f) => f.gene), covariates: s.covariates, strata: s.strata, outcome: [months, status] }, + }); + return { table, histograms, queries: logged() }; +} + // ---- the view --------------------------------------------------------------------------------- /** 0.8 → "80th", 0.333 → "33.3rd" */ @@ -395,6 +439,16 @@ export async function runRiskmodel(view: RiskmodelView, ctx: ViewRun): Promise x == null || (Array.isArray(x) && x.every((v) => typeof v === "string")); + if (typeof describe !== "object" || Array.isArray(describe) || !isList(describe.attributes) || !isList(describe.genes)) + return refuse("riskmodel", ctx, "describe is { attributes?: [names], genes?: [symbols] }, each list optional (left out, Table 1 uses its default)"); + for (const a of describe.attributes || []) if (!/^[A-Za-z0-9_:/.-]{1,128}$/.test(a)) return refuse("riskmodel", ctx, "“" + a.slice(0, 40) + "” is not a clinical attribute's name"); + for (const x of describe.genes || []) if (!/^[A-Za-z0-9][A-Za-z0-9_.-]{0,39}$/.test(x)) return refuse("riskmodel", ctx, "“" + x.slice(0, 40) + "” is not a gene symbol"); + const e = overLimit("Table 1", TABLE1_MAX_ATTRIBUTES, (describe.attributes || []).length) || overLimit("Table 1", TABLE1_MAX_GENES, (describe.genes || []).length); + if (e) return refuse("riskmodel", ctx, e.replace("at most " + TABLE1_MAX_ATTRIBUTES, "at most " + TABLE1_MAX_ATTRIBUTES + ((describe.attributes || []).length > TABLE1_MAX_ATTRIBUTES ? " attributes" : " genes"))); + } if (cohorts[name].count * runs > RISK_WORK_LIMIT) return refuse("riskmodel", ctx, "a risk model fits at most " + RISK_WORK_LIMIT.toLocaleString("en-US") + " patients × runs in one request; “" + name + "” has " + cohorts[name].count.toLocaleString("en-US") + " patients, so at most " + Math.max(1, Math.floor(RISK_WORK_LIMIT / Math.max(1, cohorts[name].count))) + " runs (or narrow the cohort)"); @@ -409,6 +463,14 @@ export async function runRiskmodel(view: RiskmodelView, ctx: ViewRun): Promise ({ study, ...q })) } : out; } const g = groupRiskmodel(fit, groupQuantiles(cutsChosen ? { cuts: cutsChosen } : { groups: nGroups })); + let t1: Awaited> | null = null; + if (describe) { + try { t1 = await riskTable1(r, args, fit, g, describe); } + catch (e) { + if (!(e instanceof RiskRefusal)) throw e; + return { ...refuse("riskmodel", ctx, e.message), queries: [...fit.queries, ...e.queries].map((q) => ({ study, ...q })) }; + } + } const cs = [...fit.cIndexRuns].filter((c) => Number.isFinite(c)).sort((x, y) => x - y); const pct = (f: number) => cs.length ? cs[Math.min(cs.length - 1, Math.round(f * (cs.length - 1)))] : NaN; @@ -423,11 +485,12 @@ export async function runRiskmodel(view: RiskmodelView, ctx: ViewRun): Promise fit.subjects.event[i] === 1).length, cIndex: fit.validationCIndex! } : null, cIndex: { median: pct(0.5), q1: pct(0.25), q3: pct(0.75), values: fit.cIndexRuns }, features: fit.features, groups: g.groups, logrank: g.logrank, refit: g.refit, histogram: g.histogram, scale: g.scale, percentiles: g.percentiles, labelAbove, cutsChosen, + ...(t1 ? { table1: t1.table } : {}), }; return { - kind: "riskmodel", cohorts: vNames, riskmodel: result, queries: fit.queries.map((q) => ({ study, ...q })), + kind: "riskmodel", cohorts: vNames, riskmodel: result, queries: [...fit.queries, ...(t1 ? t1.queries : [])].map((q) => ({ study, ...q })), method: riskmodelMethod({ study, cohort: name, cohortPatients: cohorts[name].count, rowsReturned: fit.rowsReturned, r: result, - lambdas: RISK_LAMBDAS, folds: RISK_FOLDS, testShare: RISK_TEST_SHARE }), + lambdas: RISK_LAMBDAS, folds: RISK_FOLDS, testShare: RISK_TEST_SHARE, table1Histograms: t1 ? t1.histograms : [] }), participation: ctx.participationFor([study], true), caveat: "Association, not cause: a gene's weight is its link to " + endpoint + " with the other candidates" + (fit.covariates.length ? " and the covariates" : "") + " held fixed, in this cohort" + (strata ? ", within each " + strata : "") + ". " + @@ -441,6 +504,7 @@ export async function runRiskmodel(view: RiskmodelView, ctx: ViewRun): Promise { test("the settings left out only cut and draw the groups: the model's settings do not read them", () => { for (const k of RISK_PRESENTATION_KEYS) { - const changed = { ...base, [k]: k === "cohorts" ? ["B"] : k === "cuts" ? [0.8] : k === "note" ? "other" : k === "labelAbove" ? 0.9 : 5 } as RiskmodelView; + const changed = { ...base, [k]: k === "cohorts" ? ["B"] : k === "cuts" ? [0.8] : k === "note" ? "other" : k === "labelAbove" ? 0.9 + : k === "describe" ? { attributes: ["STAGE"], genes: ["KEAP1"] } : 5 } as RiskmodelView; assert.equal(key(changed), key(base), k + " is left out of the key"); assert.deepEqual(riskFitSpec(changed, "v7"), riskFitSpec(base, "v7"), k + " must not change what the fit reads"); } diff --git a/test/riskmodel-table1.test.ts b/test/riskmodel-table1.test.ts new file mode 100644 index 0000000..9dd99fd --- /dev/null +++ b/test/riskmodel-table1.test.ts @@ -0,0 +1,157 @@ +// Table 1 of the risk model (lib/views/riskmodel-table1.ts buildTable1, lib/views/riskmodel.ts riskTable1): +// each risk group's makeup, counted from patients by position. Synthetic patients, so every count +// can be read off — above all the panel-aware denominator: a patient whose panels do not cover a +// gene is not profiled for it, never unaltered. +import { test, describe } from "node:test"; +import assert from "node:assert/strict"; +import { buildTable1, table1DefaultCandidates, TABLE1_DEFAULT_ATTRIBUTES, TABLE1_DEFAULT_MAX_VALUES, type Table1Input, type Table1Patient } from "../lib/views/riskmodel-table1.ts"; +import { riskTable1, RiskRefusal, type RiskFit, type RiskGrouping, type RiskFitArgs } from "../lib/views/riskmodel.ts"; +import { chiSquareTest } from "../lib/stats.ts"; +import type { RiskTable1Row } from "../lib/types.ts"; + +const patient = (o: Partial = {}): Table1Patient => ({ panels: ["IMPACT468"], genes: [], d: [], r: [], u: [], ...o }); + +/** Two groups of patients; attributes and genes as given. */ +function input(patients: Table1Patient[], o: Partial = {}): Table1Input { + const half = Math.ceil(patients.length / 2); + return { + study: "s", reportedOn: "held-out scores", measure: "any", + groups: [{ name: "Low risk", members: patients.map((_, i) => i).slice(0, half) }, { name: "High risk", members: patients.map((_, i) => i).slice(half) }], + patients, attributes: [], attributesDefault: false, genes: [], genesDefault: false, + panelGenes: new Map([["IMPACT468", new Set(["KEAP1", "TP53"])], ["IMPACT341", new Set(["TP53"])]]), + model: { genes: [], covariates: [], strata: null, outcome: ["OS_MONTHS", "OS_STATUS"] }, + ...o, + }; +} +const row = (rows: RiskTable1Row[], name: string) => rows.find((x) => x.name === name)!; + +describe("gene rows: the panel-aware denominator", () => { + test("a patient whose panels omit the gene is not profiled, not unaltered; whole-exome covers every gene", () => { + const ps = [ + patient({ genes: ["KEAP1"] }), patient(), patient({ panels: ["IMPACT341"] }), patient({ panels: ["IMPACT341"] }), // low: 1 of 2 profiled + patient({ panels: ["WES"], genes: ["KEAP1"] }), patient({ genes: ["KEAP1"] }), patient({ panels: ["IMPACT341"] }), patient(), // high: 2 of 3 profiled + ]; + const t = buildTable1(input(ps, { genes: ["KEAP1"] })).table; + const k = row(t.rows, "KEAP1") as Extract; + assert.deepEqual(k.perGroup, [{ altered: 1, profiled: 2, notProfiled: 2 }, { altered: 2, profiled: 3, notProfiled: 1 }]); + assert.equal(k.p, chiSquareTest([[1, 1], [2, 1]]).p, "tested on mutated vs not among the profiled only"); + }); + + test("a gene on one panel among several is counted over that panel's patients only, not the group", () => { + // the dangerous edit: dividing by the group's size. Every patient here is on IMPACT341, which lacks KEAP1 + const ps = Array.from({ length: 6 }, () => patient({ panels: ["IMPACT341"] })); + const t = buildTable1(input(ps, { genes: ["KEAP1"] })).table; + const k = row(t.rows, "KEAP1") as Extract; + assert.deepEqual(k.perGroup.map((c) => c.profiled), [0, 0]); + assert.equal(k.p, null); + assert.match(String(k.untested), /sequenced on a panel covering it/); + }); + + test("a symbol resolves case-insensitively, once; one on no panel and mutated in nobody is set aside by name", () => { + const ps = [patient({ genes: ["KEAP1"] }), patient(), patient(), patient({ genes: ["KEAP1"] })]; + const t = buildTable1(input(ps, { genes: ["keap1", "KEAP1", "NOTAGENE"] })).table; + assert.deepEqual(t.rows.map((x) => x.kind + ":" + x.name), ["gene:KEAP1", "skipped:NOTAGENE"]); + assert.deepEqual(t.describe.genes, ["KEAP1", "NOTAGENE"]); + }); +}); + +describe("rows the groups differ in by construction", () => { + const ps = Array.from({ length: 40 }, (_, i) => patient({ + genes: i >= 20 && i % 2 ? ["KEAP1"] : i % 5 === 0 ? ["TP53"] : [], + d: [i < 20 ? "never" : i % 3 ? "former" : "never", String(i % 4 ? "female" : "male"), i < 20 ? "0:living" : "1:deceased"], + r: [i < 20 ? "Never" : i % 3 ? "Former" : "Never", String(i % 4 ? "Female" : "Male"), i < 20 ? "0:LIVING" : "1:DECEASED"], u: [1, 1, 1], + })); + const meta = { numeric: false, patient: true }; + const t = buildTable1(input(ps, { + attributes: [{ name: "SMOKING", meta: { id: "SMOKING", ...meta } }, { name: "SEX", meta: { id: "SEX", ...meta } }, { name: "OS_STATUS", meta: { id: "OS_STATUS", ...meta } }], + genes: ["KEAP1", "TP53"], model: { genes: ["KEAP1"], covariates: ["SMOKING"], strata: "SEX", outcome: ["OS_MONTHS", "OS_STATUS"] }, + })).table; + + test("the model's genes and covariates are 'in the model', its endpoint 'outcome', its stratum 'stratum'", () => { + assert.deepEqual(t.rows.map((x) => x.name + "=" + ("role" in x ? x.role : "")), ["SMOKING=in the model", "SEX=stratum", "OS_STATUS=outcome", "KEAP1=in the model", "TP53=null"]); + }); + + test("they keep a p-value but get no q-value and stay out of the Benjamini–Hochberg family; the stratum stays in", () => { + for (const n of ["SMOKING", "OS_STATUS", "KEAP1"]) { + const x = row(t.rows, n) as Exclude; + assert.ok(x.p != null && x.q == null, n); + } + assert.equal(t.tested, 2, "SEX and TP53"); + // with two rows in the family, the larger p is its own q + const sex = row(t.rows, "SEX") as Exclude, tp53 = row(t.rows, "TP53") as Exclude; + assert.equal(Math.max(sex.q!, tp53.q!), Math.max(sex.p!, tp53.p!)); + }); + + test("categories fold, and show their commonest spelling", () => { + const smoking = row(t.rows, "SMOKING") as Extract; + assert.deepEqual(smoking.levels.map((l) => l.value).sort(), ["Former", "Never"]); + assert.deepEqual(smoking.missing, [0, 0]); + }); +}); + +describe("clinical rows", () => { + test("numbers: medians per group and Kruskal–Wallis; a value that is not a number is missing", () => { + const ps = [10, 20, 30, 40, 50, 60, 70, 80].map((v, i) => patient({ d: [i === 1 ? "na" : String(v)], r: [String(v)], u: [i === 7 ? 2 : 1] })); + const t = buildTable1(input(ps, { attributes: [{ name: "AGE", meta: { id: "AGE", numeric: true, patient: false } }] })); + const age = row(t.table.rows, "AGE") as Extract; + assert.deepEqual(age.perGroup.map((c) => [c.n, c.median, c.missing]), [[3, 30, 1], [4, 65, 0]]); + assert.equal(age.test, "Kruskal-Wallis"); + assert.equal(age.multiValued, 1); + assert.equal(age.per, "sample"); + assert.equal(t.histograms.filter((h) => h.row === "attribute:AGE").reduce((s, h) => s + h.n, 0), 7); + }); + + test("an attribute the study does not have is set aside by name, and one with no values among these patients too", () => { + const ps = [patient({ d: ["", ""], r: ["", ""], u: [0, 0] }), patient({ d: ["", ""], r: ["", ""], u: [0, 0] })]; + const t = buildTable1(input(ps, { attributes: [{ name: "NOPE", meta: null }, { name: "EMPTY", meta: { id: "EMPTY", numeric: false, patient: true } }] })).table; + assert.deepEqual(t.rows.map((x) => x.kind === "skipped" ? x.skipped : ""), ["not a clinical attribute of s", "no value recorded for these patients"]); + }); +}); + +describe("the default attributes", () => { + test("the stratum and covariates come first, survival endpoints' columns never", () => { + assert.deepEqual(table1DefaultCandidates(["OS_STATUS", "STAGE", "AGE", "PFS_MONTHS", "SEX"], { covariates: ["AGE"], strata: "SEX" }), ["SEX", "AGE", "STAGE"]); + }); + + test("one with a single value, none, or too many categories is left out with why; at most the default count are kept", () => { + const names = ["CONST", "WIDE", "EMPTY", ...Array.from({ length: TABLE1_DEFAULT_ATTRIBUTES + 2 }, (_, k) => "A" + k)]; + const ps = Array.from({ length: 30 }, (_, i) => patient({ + d: names.map((n) => n === "CONST" ? "x" : n === "WIDE" ? "v" + i : n === "EMPTY" ? "" : String(i % 2)), + r: names.map(() => "v"), u: names.map(() => 1), + })); + const t = buildTable1(input(ps, { attributesDefault: true, attributes: names.map((name) => ({ name, meta: { id: name, numeric: false, patient: true } })) })).table; + assert.deepEqual(t.omitted.map((o) => o.name), ["CONST", "WIDE", "EMPTY"]); + assert.match(t.omitted[1].why, new RegExp("takes 30 values here \\(a default takes at most " + TABLE1_DEFAULT_MAX_VALUES + "\\)")); + assert.deepEqual(t.describe.attributes, names.slice(3, 3 + TABLE1_DEFAULT_ATTRIBUTES)); + assert.equal(t.defaults.attributes, true); + }); +}); + +describe("Table 1 reads the fit's own patients", () => { + // a resolver that answers the rows query and the metadata, as a patient-level participant would + const rowsOf = (months: number[]) => months.map((m, i) => ({ m, ev: i % 2, half: 0, panels: ["IMPACT468"], genes: i % 3 ? [] : ["KEAP1"], d0: i % 2 ? "a" : "b", r0: i % 2 ? "A" : "B", u0: 1 })); + const resolver = (rows: Record[]) => ({ + queries: [] as { label: string; sql: string }[], + viewRows: async () => rows, + viewMeta: async (label: string) => label === "riskmodel:panels" ? [{ panel: "IMPACT468", genes: ["KEAP1"] }] : [{ id: "X", datatype: "STRING", patient: 1, priority: 1 }], + }); + const fit = (months: number[]) => ({ + n: months.length, validated: false, features: [{ gene: "KEAP1", selected: 1, meanCoef: 1, hr: 2.7, mutated: 2 }], + subjects: { time: months, entry: months.map(() => 0), event: months.map((_, i) => (i % 2) as 0 | 1), stratum: months.map(() => 0) }, + }) as unknown as RiskFit; + const grouping = { groups: [{ name: "Low risk" }, { name: "High risk" }], members: [[0, 1, 2], [3, 4, 5]] } as unknown as RiskGrouping; + const args = { study: "s", def: { all: [] }, compile: { oncokbPin: null, profiles: [], gaProfiles: [] }, + view: { kind: "riskmodel", endpoint: "OS", minMutated: 5, runs: 1 } } as unknown as RiskFitArgs; + const months = [1, 2, 3, 4, 5, 6]; + + test("the same rows give each group its members' values", async () => { + const t = await riskTable1(resolver(rowsOf(months)) as never, args, fit(months), grouping, { attributes: ["X"] }); + const x = row(t.table.rows, "X") as Extract; + assert.deepEqual(x.levels.map((l) => [l.value, ...l.counts]), [["A", 1, 2], ["B", 2, 1]]); + }); + + test("rows that are not the fit's patients, in the fit's order, are refused rather than described", async () => { + await assert.rejects(riskTable1(resolver(rowsOf([1, 2, 3, 4, 6, 5])) as never, args, fit(months), grouping, {}), (e) => e instanceof RiskRefusal && /matched to the model's patients/.test(e.message)); + await assert.rejects(riskTable1(resolver(rowsOf(months.slice(1))) as never, args, fit(months), grouping, {}), RiskRefusal); + }); +}); From 8bc536f4dcda615d7d4936913d39a4c9205b425b Mon Sep 17 00:00:00 2001 From: Aaron Lisman Date: Mon, 19 Oct 2026 00:06:01 +0000 Subject: [PATCH 3/6] Risk model: the figure's caption states Table 1 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The figure does not draw Table 1 (the page shows it in a dialog), so the caption the server SVG carries — MCP, exports — says it in words when the view has `describe`: the groups and their sizes, how many rows and how many are in the q-value family, then each row's headline per group (a category's commonest value as a share of patients with a value, a number's median, a gene's share mutated among patients whose panels cover it) with its q, in three parts: differing at q < 0.05, not, and different by construction (in the model / outcome, with p and no q). Defaults passed over and set-aside rows are named with why. A stress document (riskmodel-table1-full: 20 attributes, 20 genes, five groups, a stratum and a covariate, including outcome columns and an unknown gene) joins scripts/figure-stress.json. Verified: FIGURE_CHECK_ONLY=risk check:figures — the gallery example and the stress document lint clean at 560, 936 and 1352px (6 figures, 0 issues); the 936px caption was rendered and read; npm test passes. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01UCpBeF9VMrjXyddFKUfiSX --- lib/figures/riskmodel.ts | 33 +++++++++++++++- scripts/figure-stress.json | 77 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 108 insertions(+), 2 deletions(-) diff --git a/lib/figures/riskmodel.ts b/lib/figures/riskmodel.ts index bc66187..60035ca 100644 --- a/lib/figures/riskmodel.ts +++ b/lib/figures/riskmodel.ts @@ -1,9 +1,9 @@ // Risk model figure: the risk groups' survival curves and the genes the model kept. import { palette, frame, caption, linear, log, axisLeft, axisBottom, leftAxisWidth, bottomAxisHeight, label, textWidth, - seriesStyle, SERIES, esc, fmtInt, fmtNum, fmtPct, fmtShare, fmtRatio, fmtP, fmtMonths, FS, SW, type SvgOpts + seriesStyle, SERIES, esc, fmtInt, fmtNum, fmtPct, fmtShare, fmtRatio, fmtP, fmtQ, fmtMonths, FS, SW, type SvgOpts } from "../svg-kit.ts"; -import type { ViewResult } from "../types.ts"; +import type { ViewResult, RiskTable1 } from "../types.ts"; import { stepPath } from "../chart/marks.ts"; import { plotHeight } from "./survival.ts"; import { timeScale, timeTitle } from "../chart/time.ts"; @@ -170,6 +170,34 @@ function riskVolcano(m: Risk, p: ReturnType, x0: number, y0: num * The risk model: Kaplan–Meier curves of the risk groups, then the genes the model kept most * often with their mean hazard ratio (a bar on a log scale about 1), then what the numbers rest on. */ +/** + * Table 1 in words, for the figure's caption: the figure does not draw it (the page shows it in a + * dialog), so an image that travels alone still says what it holds and what it found. + */ +export function table1Words(t: RiskTable1): string { + const groups = t.groups.map((g) => g.name.replace(/ risk$/, "")).join(" / "); + const shown = t.rows.filter((r) => r.kind !== "skipped"); + const share = (n: number, of: number) => of ? fmtPct(100 * n / of) : "–"; + const summary = (r: typeof shown[number]): string => { + if (r.kind === "gene") return r.name + " mutated in " + r.perGroup.map((c) => share(c.altered, c.profiled)).join(" / ") + " of patients whose panels cover it"; + if (r.kind === "num") return r.name + " median " + r.perGroup.map((c) => c.median == null ? "–" : fmtNum(c.median)).join(" / "); + const top = r.levels[0]; + return r.name + " " + top.value + " in " + top.counts.map((n, gi) => share(n, top.counts[gi] + r.levels.slice(1).reduce((s, l) => s + l.counts[gi], 0))).join(" / "); + }; + const stat = (r: typeof shown[number]) => r.q != null ? " (" + fmtQ(r.q) + ")" : r.p != null ? " (" + fmtP(r.p) + ")" : r.untested ? " (not tested: " + r.untested + ")" : ""; + const open = shown.filter((r) => r.role !== "in the model" && r.role !== "outcome"); + const differ = open.filter((r) => r.q != null && r.q < 0.05), rest = open.filter((r) => !differ.includes(r)); + const built = shown.filter((r) => r.role === "in the model" || r.role === "outcome"); + const set = t.rows.filter((r) => r.kind === "skipped"); + return "Table 1 (not drawn; in the result) compares the groups (" + groups + ", " + t.groups.map((g) => fmtInt(g.n)).join(" / ") + " patients" + + (t.reportedOn === "validation half" ? ", the validation half" : "") + ") on " + shown.length + " rows, each tested across groups (χ² or Kruskal–Wallis), q-values Benjamini–Hochberg over " + t.tested + " of them; a category's share is of the patients with a value. " + + (differ.length ? "Differ at q < 0.05: " + differ.map((r) => summary(r) + stat(r)).join("; ") + ". " : "No open row differs at q < 0.05. ") + + (rest.length ? "Not at q < 0.05: " + rest.map((r) => summary(r) + stat(r)).join("; ") + ". " : "") + + (built.length ? "Different by construction, so no q-value: " + built.map((r) => summary(r) + " (" + r.role + (r.p != null ? ", " + fmtP(r.p) : "") + ")").join("; ") + ". " : "") + + (t.omitted.length ? "Left out of the default attributes: " + t.omitted.map((o) => o.name + " (" + o.why + ")").join("; ") + ". " : "") + + (set.length ? "Set aside: " + set.map((r) => r.name + " (" + (r as { skipped: string }).skipped + ")").join("; ") + ". " : ""); +} + export function riskmodelSvg(v: ViewResult, o: SvgOpts): string { const m = v.riskmodel!; const p = palette(o.theme); @@ -294,6 +322,7 @@ export function riskmodelSvg(v: ViewResult, o: SvgOpts): string { (m.candidates > m.features.length ? " " + fmtInt(m.candidates - m.features.length) + " of the " + fmtInt(m.candidates) + " candidate genes were never kept, so are not on the volcano." : "") + (left.length ? " Left out: " + left.join("; ") + "." : "") + (m.features.length > shown.length ? " " + (m.features.length - shown.length) + " more genes were kept in at least one run; they are in the result, not drawn." : "") + + (m.table1 ? " " + table1Words(m.table1) : "") + " " + (v.caveat || ""); const cap = caption(text, W, yy + 4, p.muted); return frame(W, yy + 4 + cap.height + 8, out.join("") + cap.svg, p.bg, grown, diff --git a/scripts/figure-stress.json b/scripts/figure-stress.json index 0ab3de1..e37d50c 100644 --- a/scripts/figure-stress.json +++ b/scripts/figure-stress.json @@ -2161,5 +2161,82 @@ "sort": "mutational_signatures_contribution_v2_Signature7" } ] + }, + "riskmodel-table1-full": { + "study": "msk_chord_2024", + "cohorts": { + "Non-small cell lung cancer, sequenced, with overall survival on record": { + "all": [ + { + "kind": "clinical", + "attribute": "CANCER_TYPE", + "values": [ + "Non-Small Cell Lung Cancer" + ], + "per": "sample" + } + ] + } + }, + "views": [ + { + "kind": "riskmodel", + "cohorts": [ + "Non-small cell lung cancer, sequenced, with overall survival on record" + ], + "runs": 2, + "groups": 5, + "strata": "SMOKING_PREDICTIONS_3_CLASSES", + "covariates": [ + "GENDER" + ], + "describe": { + "attributes": [ + "SMOKING_PREDICTIONS_3_CLASSES", + "GENDER", + "STAGE_HIGHEST_RECORDED", + "RACE", + "ETHNICITY", + "SAMPLE_TYPE", + "MSI_TYPE", + "CLINICAL_GROUP", + "CLINICAL_SUMMARY", + "PATHOLOGICAL_GROUP", + "CANCER_TYPE_DETAILED", + "ICD_O_HISTOLOGY_DESCRIPTION", + "CURRENT_AGE_DEID", + "TMB_NONSYNONYMOUS", + "FRACTION_GENOME_ALTERED", + "MUTATION_COUNT", + "OS_STATUS", + "OS_MONTHS", + "PRIOR_MED_TO_MSK", + "HISTORY_OF_PDL1" + ], + "genes": [ + "KEAP1", + "TP53", + "STK11", + "EGFR", + "KRAS", + "SMARCA4", + "NFE2L2", + "ZFHX3", + "SMARCA2", + "PTPRD", + "EPHA3", + "TGFBR2", + "PTEN", + "CDKN2A", + "RB1", + "ARID1A", + "NF1", + "BRAF", + "MET", + "NOTAGENE1" + ] + } + } + ] } } From 07a6f51fc87e27c02215c0175dc90e61af24bf20 Mon Sep 17 00:00:00 2001 From: Aaron Lisman Date: Mon, 19 Oct 2026 00:21:02 +0000 Subject: [PATCH 4/6] Risk model: a Table 1 dialog, with rows to add and "Add to document" MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A "Table 1" button in the riskmodel card's toolbar (its own zone, HTML beside the figure) opens every risk group's makeup side by side: a category's patients per value (top 8, the rest pooled as a disclosed "other"; not recorded shown), a number's median [IQR], a gene's mutated over profiled with the not-profiled count; p, q, and the row's role ("in the model", "outcome", "stratum", each titled with what it means) and a "sparse" flag below an expected count of 5. Clicking a group's header marks its column. The dialog owns its draft `describe` (attributes from the study's catalog in a PickMenu, genes typed in and removed as chips, "defaults" to go back): a change re-runs only this request — the view's one-view document with the draft, through /api/analysis, so the model is the cached fit — and the page re-renders nothing. While it runs the table is covered ("updating Table 1…") and the controls stay usable; a failure says why with retry. When the draft is what the document holds, the view's own result is shown with no request. "Add to document" writes the resolved lists into the view (patchView), so a link and the caption carry it; "Remove from document" takes them out. e2e (riskmodel Table 1 block): one button, aligned (misalignedControls on the toolbar and the dialog); groups side by side, a header click marks a column; the model's genes marked with no q; adding SMARCA2 with the request held covers the table and leaves the controls usable, then is one request and shows not-profiled counts; "Add to document" puts describe in the URL and "Table 1 (not drawn; in the result)" in the caption (caption lines joined as read, since it wraps); reopening is served from the view's result with no request; no JS errors. Verified: the block passes against this folder's dev server on cold seeds (13; and 11, which first exposed the caption-wrap flaw in the check); screenshots of the toolbar and dialog at 2x read and compared with the risk-groups dialog; tsc and eslint clean. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01UCpBeF9VMrjXyddFKUfiSX --- app/lab/kind-cards.tsx | 2 + app/lab/risk-table1-dialog.tsx | 230 +++++++++++++++++++++++++++++++++ e2e.cjs | 70 ++++++++++ 3 files changed, 302 insertions(+) create mode 100644 app/lab/risk-table1-dialog.tsx diff --git a/app/lab/kind-cards.tsx b/app/lab/kind-cards.tsx index 4304f3c..980927d 100644 --- a/app/lab/kind-cards.tsx +++ b/app/lab/kind-cards.tsx @@ -14,6 +14,7 @@ import { Input } from "../_ui/input"; import { OptionSelect } from "../_ui/choice.tsx"; import { FIELD } from "../_ui/toolbar.ts"; import { RISK_LABEL_ABOVE } from "../../lib/figures/riskmodel.ts"; +import { RiskTable1Button } from "./risk-table1-dialog.tsx"; import { RiskGroupsButton } from "./risk-groups-dialog.tsx"; /** @@ -458,6 +459,7 @@ export function KindCardControls({ c }: { c: KindCardCtx }) { +
view +
-
view -
-
adjust a.attribute !== spec.strata).map((a) => ({ value: a.attribute, label: a.attribute, hint: a.kind === "num" ? "number" : "category" }))} chosen={spec.covariates || []} max={6} busy={busy} menuId={"risk-cov:" + c.vi} menuStore={RISK_MENUS} @@ -444,19 +475,7 @@ export function KindCardControls({ c }: { c: KindCardCtx }) { }
model - - - +
diff --git a/e2e.cjs b/e2e.cjs index 850f585..a649ef5 100644 --- a/e2e.cjs +++ b/e2e.cjs @@ -1073,6 +1073,15 @@ async function main() { ok("risk model: draws the 0–10 score and the volcano", /score \(0–10\)/.test(await text()) && /Selection frequency against mean hazard ratio/.test(await text())); const askew = await misalignedControls(page, ".riskmodel-tb"); ok("risk model: its controls share one height and centre line", askew.length === 0, askew.join(" | ")); + // the fitting settings sit in the view's settings menu; a value the menu does not offer still shows + await page.locator(".risksettingsbtn").click(); + const settings = page.locator("[data-slot=popover-content]"); + await settings.waitFor({ timeout: 10000 }); + const runsShown = (await settings.locator(".setrow", { hasText: "runs" }).locator("[role=combobox]").innerText()).trim(); + const settingsText = (await settings.innerText()).replace(/\s+/g, " "); + await page.keyboard.press("Escape"); + ok("risk model: penalty, runs and validation are in the settings menu, showing the document's runs", + runsShown === "2" && /penalty/.test(settingsText) && /validate on held-out half/.test(settingsText), "runs " + JSON.stringify(runsShown)); let requests = 0; page.on("request", (r) => { if (r.url().includes("/api/analysis")) requests++; }); await page.locator(".riskmodel-tb label.tbfield", { hasText: "label genes" }).locator("[role=combobox]").click();