A review of the ranking's arithmetic found four things wrong, all of which
made the score look better informed than it was. Measurements below are from
this repo, not estimates.
**Components now abstain instead of inventing a number.** A run without a VEP
cache returns no allele frequencies, and rarity_score(None) read that as
"absent from gnomAD, therefore maximally rare" and awarded every variant a
free 0.25. jobs.has_frequencies / has_effect_scores record what the run
actually produced, absent components are dropped from the weighted mean, and
the remaining weights are renormalised so the score keeps its meaning. The UI
shows "not looked up" rather than a bar, and the funnel stops calling a step
"rare" when nothing was filtered.
**Allele frequency is no longer a model feature.** It dominated: the same
missense variant scored 0.887 at AF 0 and 0.0003 at AF 0.01. That double-
counted, because the ranking already scores frequency explicitly, putting
~45% of every rank on one measurement; and it was circular, because ACMG
assigns ClinVar's benign labels using frequency (BA1/BS1). Retraining without
it moves missense AUROC from 0.872 to 0.500 — exactly random. The old figure
was allele frequency, not variant-effect knowledge. The model therefore
abstains unless CADD or AlphaMissense is present, since otherwise it only
restates the consequence class.
**Phenotype matching is weighted by information content** and HPO annotations
are propagated up the ontology. Counting terms alike let "global
developmental delay" (IC 0.93) count as much as "dilated left subclavian
artery" (IC 7.88).
**A real bug in the propagation, found by checking it.** The ancestor walk
read a pre-order DFS backwards, which on a DAG lets a term resolve before one
of its parents and inherit that parent alone instead of its lineage. It
dropped 399 terms out of the phenotype branch, Camptodactyly and Chiari
malformation among them. Now a true post-order, tested against a reference
transitive closure.
The ontology arithmetic moved to rarelens_ml.hpo so it is covered by tests,
and rarelens_ml.benchmark measures the whole thing: across 10,178 published
cases the causal gene ranks first 45.9-81.0% of the time against 5,269 genes,
versus 0.02% for chance. docs/data.md reports that with its contamination
(HPO's annotations come from these same case reports), and includes the
measurement showing information-content weighting earns its place while
propagation does not - kept anyway, for a reason the docs argue rather than
assume.
Five things had drifted:
- The tech table listed GKE and Cloud SQL as the cloud, but the default
deployment track has been serverless Cloud Run + Google Batch since
docs/cloud.md; GKE and Cloud SQL are behind flags. Cloud Run was missing
entirely, including as the Nextflow driver.
- The button is "Analyse case", not "Run VEP annotation" (twice).
- samples were renamed to cases, so LOCAL_DATA_ROOT holds a case's vcf_uri.
- The Data section omitted HPO gene-to-phenotype, which is half the ranking,
and the phenopackets behind the published case.
- The status list described the sample list and filtered variant table the UI
no longer has, and read as pending work when all of it ships. Replaced with
what actually works, plus the two known gaps: no allele frequencies without
the VEP cache, and no authentication on the deployed demo.
uv generated it while running ruff and pytest. The repo already tracks
web/package-lock.json, so pinning the resolved Python dependencies the same
way keeps CI and local runs reproducible instead of re-resolving each time.
`make published-case` reads a GA4GH phenopacket from Monarch's Phenopacket
Store and takes two things from it verbatim: the HPO terms the authors
reported and the variant they called causal. The default is the TGFBR2
proband from Loeys et al., Nat Genet 2005 (doi:10.1038/ng1511), the paper
that first defined Loeys-Dietz syndrome -- 30 reported terms and
NM_003242.6:c.1069G>T p.(Gly357Trp).
The rest of that patient's genome is not public, so background variants come
from GIAB HG002 around the locus. They are drawn from coding exons where
possible, via Ensembl's REST API: of ~4,000 HG002 variants in the window only
9 are coding, so a random sample is entirely intronic, the consequence filter
discards all of it, and the causal variant is left as the only candidate --
a funnel that proves nothing.
The real run ranks TGFBR2 first at 0.897 against an OSBPL10 missense at
0.547. Both are rare missense variants the model scores identically (0.887);
only the phenotype separates them, which is the argument for phenotype-driven
triage in one table.
Documented with three caveats rather than left implicit: the phenotype match
is partly circular because HPO's gene annotations are themselves curated from
published cases; rarity contributes nothing without the VEP cache
(--af_gnomade is rejected with --database, and plain --af returns nothing
even for rs429358 at ~15% global frequency); and one healthy genome is not a
diagnostic exome.
A case analysed before a model existed stayed unscored forever: the only control on the page
was "Re-analyse case", which re-runs five minutes of VEP to obtain a score that takes a second.
When any candidate is unscored the page now offers "Score variants" on its own.
Tests: web 34; svelte-check clean.
"Variants are unscored" was accurate: nothing was ever trained, so a quarter of every rank was
dead weight and the UI leaked a connection error at the reader.
- scripts/make-training-set.sh derives a training table from ClinVar directly. ClinVar already
carries the molecular consequence, the gene and an allele frequency, which is the feature set
serving sends, so this avoids running VEP over hundreds of thousands of variants. 2-star
records only.
- train.py now holds out whole genes (GroupShuffleSplit). docs/data.md had said to do this since
the data pass; the code was still doing a random split, which is the leak Grimm 2015 describes.
- evaluate() reports missense on its own. On the last run: AUROC 0.986 over 74,239 held-out
variants, but 0.872 over the 13,553 missense ones, and the docs say plainly why even that is
flattered — within missense the only live feature is allele frequency, and ClinVar's benign
calls often use allele frequency as evidence (ACMG BA1/BS1), so the feature partly caused the
label.
- the 503 now names what is missing (model@alias via tracking URI) and leaves the exception in
the server log instead of the UI.
- make training-set / make train; the 58 MB table is gitignored.
Verified end to end: model registered as v2, the simulated NF2 case scores 0.999 on the planted
variant, and it now ranks 1.00 with all four components live.
Tests: api 77, ml 22, loader 16, web 32; ruff, mypy, svelte-check clean.
From clicking through the redesigned UI:
- scoring a case without a model registry painted a red failure across a case that had in
fact analysed fine. It is now a quiet note saying the model term contributes 0, because
scoring is an optional fourth of the rank, not the analysis.
- the MLflow default moves to port 5001. On macOS, AirPlay Receiver owns 5000, which is why
the registry answered "403" rather than refusing the connection; docker-compose publishes
5001 to match.
- a funnel step that kept nothing drew a visible bar. Zero now draws zero.
- "1 candidates".
- the funnel's fixed grid columns forced a horizontal scrollbar on the report.
The report also lists the top undecided candidates now: the first thing anyone opens has no
decisions in it, and "Shortlisted (0)" alone said nothing about what the tool found.
Tests: api 77, web 32; ruff, mypy, svelte-check clean.
A table with filters made the user do the work. Rare disease triage is a different task:
which few variants could explain *this* patient's phenotype, and why. The app now answers
that, and lets a reviewer act on the answer.
Domain
- a case is a proband: a VCF plus the HPO terms observed in the patient (samples -> cases)
- HPO's gene-to-phenotype annotations are loaded as reference data (scripts/load-hpo.py)
- each candidate can be shortlisted or dismissed with a reason and a note
Ranking (app/services/triage.py, 21 tests)
- weighted sum of phenotype match, rarity, consequence severity and the model's score,
with every component shown next to the candidate
- rarity and consequence filter; phenotype only ranks, because a real diagnosis can sit in
a gene nobody has annotated yet and filtering on it would hide exactly that case
- ClinVar is deliberately not an input: it appears beside the result as independent
confirmation, so nothing ranks highly merely because ClinVar already said pathogenic
UI
- the funnel is the headline: variants called -> rare -> coding candidates -> phenotype-matched
- ranked candidates with evidence chips, not a grid of everything; filters are demoted
- a variant panel showing the score breakdown, the matched HPO terms, the raw VEP record and
links out to Ensembl/gnomAD/ClinVar, with the decision controls
- a printable case report: phenotype, funnel, shortlisted variants with reasons, provenance
API: /cases with phenotypes, /cases/{id}/candidates (funnel + ranked + weights),
/variants/{id}, /variants/{id}/decision, /cases/{id}/report, /phenotypes for the picker.
Scoring moved under the case and now answers 503 with the reason when no model registry is
reachable, instead of a 500.
Verified end to end on a simulated proband (scripts/make-demo-case.sh: real GIAB HG002
background + one real ClinVar 2-star pathogenic NF2 variant). 13 variants called -> 1 coding
candidate, and the planted variant ranks first at 0.80 on phenotype 1.00, rarity 1.00 and
consequence 1.00, with ClinVar agreeing afterwards.
Tests: api 75, ml 18, loader 16, web 27; ruff, mypy, svelte-check, terraform validate, both
kustomize overlays and the Nextflow stub run all clean.
"running" for two and a half minutes tells the user nothing. The job page now shows a
spinner, the elapsed time, and the pipeline step Nextflow is actually on.
- the API streams the Nextflow output into jobs.log as it arrives, instead of keeping it
only when the run dies. Writes are throttled to one every 3s, or immediately when a new
process starts, and are skipped for a job that has already finished, so a late line
cannot overwrite the loader's result.
- web/src/lib/progress.ts formats the elapsed time and picks the latest [PROCESS] line.
No percentage: the pipeline cannot honestly estimate one.
- the spinner animates only under prefers-reduced-motion: no-preference.
Verified on a live run: the page showed "VEP (tiny)" for the duration, then the variant
table replaced it on success.
Tests: api 53, web 20; ruff, mypy, svelte-check clean.
Makes a real annotation runnable locally without the 25 GB VEP cache, which is what
the demo needs and what a reviewer can reproduce in minutes.
- params.vep_database (VEP_DATABASE=true) queries Ensembl's public database instead of
a local cache. Slower per variant and fewer fields, so --everything is swapped for the
flags the loader actually stores. Its cache placeholder is NO_CACHE, not NO_FILE:
Nextflow rejects two staged inputs sharing a filename.
- PIPELINE_DATABASE_URL is handed to the pipeline when set. The loader runs inside a
container, where the API's own localhost URL would point at the container itself.
- README: how to run the UI's annotate button locally against host Nextflow + Docker.
Verified end to end on pipeline/tests/data/tiny.vcf: bcftools norm split the multiallelic
record, VEP 113 annotated 4 variants live, the loader wrote them and marked the job
succeeded, and the UI shows them. The deletion came back as 22:42126611 CT>C with exact
VCF alleles, which is the case the audit's ID-tagging fix exists for.
Tests: api 51, loader 16, stub run 3/3; ruff, mypy clean.
An end-to-end audit found the repo could not build, test or run as shipped. This
fixes every finding, then adds a Cloud Run track so the demo costs about £1/month
idle instead of ~£150.
CI (red on its first run)
- api: setuptools could not build the package (flat layout with app/ and alembic/)
- web: missing @types/node; `vitest run` exited 1 with no test files
- pipeline: the stub run needed a gitignored VCF, and no process had a stub block
- ruff pinned, mypy configured, DB tests on real Postgres (pgserver locally, service in CI)
ML serving (scores were meaningless)
- the registered model now carries its own feature engineering and returns predict_proba,
so serving sends raw columns and cannot drift from training
- resolve by registry alias (stages are deprecated in MLflow 3) and record the real
version; re-scoring upserts instead of failing on the unique constraint
- ClinVar labels parsed from VEP's lowercase terms
Pipeline
- exact ref/alt recovered from a CHROM_POS_REF_ALT VCF ID; loading is idempotent
- job status reaches running/failed/succeeded, so the UI stops polling dead jobs
- DATABASE_URL travels in the environment or a Nextflow secret, never on a command line
- VEP cache and plugins staged as inputs; the gcp profile runs tasks on Google Batch
Deployment
- the API serves /api (matching the ingress); the web app reads its API URL at runtime
- migrations run in an init container under a Postgres advisory lock
- terraform: custom VPC shared with Batch, private Cloud SQL, API enablement, Workload
Identity bindings, Secret Manager, deletion protection
- serverless track, now the default: Cloud Run services scaling to zero, a Cloud Run job
for the Nextflow driver, and Neon or Cloud SQL behind one DATABASE_URL secret. GKE and
Argo remain, behind -var deploy_kubernetes=true. See docs/cloud.md.
Correctness and security
- 409 on duplicate sample names, 422 on bad paging, natural chromosome ordering, wider
VEP text columns, enum dropped on downgrade, the sample's assembly actually used
- vcf_uri restricted to gs:// objects or files under the data root, blocking option injection
- CORS restricted to configured origins; `make down` no longer deletes volumes
Data
- docs/data.md records the peer-reviewed, openly licensed sources (GIAB HG002, ClinVar,
gnomAD) with citations and an honest evaluation plan; `make data` fetches a chr22 slice
Verified: api 50 tests, ml 18, loader 16, web 12; ruff, mypy, svelte-check, terraform
validate and both kustomize overlays clean.
Kemal Yaylali
kemal