feat: redesign around phenotype-driven triage, not variant filtering

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.
2026-09-12 08:30:44 +01:00
parent abde5ec6e4
commit 07a01715fd
47 changed files with 2159 additions and 539 deletions
@@ -0,0 +1,65 @@
+#!/usr/bin/env python3
+"""Load HPO's gene-to-phenotype annotations into the gene_phenotypes table.
+
+This is the reference data the ranking matches a case's phenotype against. Source file:
+https://purl.obolibrary.org/obo/hp/hpoa/genes_to_phenotype.txt (HPO release, ~20 MB).
+
+Cite the Human Phenotype Ontology when showing results; see docs/data.md.
+"""
+import argparse
+import csv
+import io
+import os
+import sys
+import urllib.request
+
+from sqlalchemy import create_engine, text
+from sqlalchemy.engine import make_url
+
+URL = "https://purl.obolibrary.org/obo/hp/hpoa/genes_to_phenotype.txt"
+
+
+def rows(handle: io.TextIOBase) -> list[tuple[str, str, str]]:
+    """Unique (gene, term) pairs; the file repeats them once per associated disease."""
+    seen: set[tuple[str, str]] = set()
+    out: list[tuple[str, str, str]] = []
+    for row in csv.DictReader(handle, delimiter="\t"):
+        gene, hpo_id, name = row["gene_symbol"], row["hpo_id"], row["hpo_name"]
+        if not gene or not hpo_id or (gene, hpo_id) in seen:
+            continue
+        seen.add((gene, hpo_id))
+        out.append((gene[:60], hpo_id[:20], name[:200]))
+    return out
+
+
+def main() -> None:
+    p = argparse.ArgumentParser()
+    p.add_argument("--url", default=URL)
+    p.add_argument("--file", help="use a local copy instead of downloading")
+    a = p.parse_args()
+
+    url = os.environ.get("DATABASE_URL")
+    if not url:
+        sys.exit("DATABASE_URL is not set")
+
+    if a.file:
+        with open(a.file) as fh:
+            annotations = rows(fh)
+    else:
+        print(f"downloading {a.url}", file=sys.stderr)
+        with urllib.request.urlopen(a.url) as response:  # noqa: S310 - fixed HPO release URL
+            annotations = rows(io.TextIOWrapper(response, encoding="utf-8"))
+    print(f"{len(annotations)} gene/term pairs", file=sys.stderr)
+
+    engine = create_engine(make_url(url).set(drivername="postgresql+psycopg"))
+    with engine.begin() as conn:
+        conn.execute(text("TRUNCATE gene_phenotypes RESTART IDENTITY"))
+        cursor = conn.connection.cursor()
+        with cursor.copy("COPY gene_phenotypes (gene_symbol, hpo_id, hpo_name) FROM STDIN") as copy:
+            for row in annotations:
+                copy.write_row(row)
+    print(f"loaded {len(annotations)} annotations", file=sys.stderr)
+
+
+if __name__ == "__main__":
+    main()
@@ -0,0 +1,52 @@
+#!/usr/bin/env bash
+# Build a SIMULATED proband VCF: real GIAB HG002 variants as background, plus one real ClinVar
+# pathogenic variant in a disease gene, which becomes the diagnosis the triage should find.
+#
+# Spiking a known variant into a public genome is how phenotype-driven triage tools are
+# benchmarked. Every input here is public and openly licensed, and this is not a real patient.
+# Provenance and citations: docs/data.md
+set -euo pipefail
+
+IMAGE=${BCFTOOLS_IMAGE:-quay.io/biocontainers/bcftools:1.20--h8b25389_0}
+GENE=${GENE:-NF2}
+REGION=${REGION:-chr22:20000000-31000000}
+BACKGROUND=${BACKGROUND:-12}          # kept small: VEP's database mode is slow per variant
+OUT_DIR=${OUT_DIR:-data}
+CLINVAR=${CLINVAR:-https://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh38/clinvar.vcf.gz}
+GIAB=${GIAB:-https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/release/AshkenazimTrio/HG002_NA24385_son/NISTv4.2.1/GRCh38/HG002_GRCh38_1_22_v4.2.1_benchmark.vcf.gz}
+
+mkdir -p "$OUT_DIR"
+docker run --rm -v "$PWD/$OUT_DIR:/out" "$IMAGE" bash -eu -c "
+  echo '==> background: GIAB HG002 $REGION' >&2
+  bcftools view -H -r '$REGION' '$GIAB' \
+    | awk '\$5 !~ /,/ && length(\$4) < 20 && length(\$5) < 20' \
+    | awk 'NR % 149 == 0' \
+    | head -n $BACKGROUND \
+    | awk -v OFS='\t' '{ sub(/^chr/, \"\", \$1); print \$1, \$2, \".\", \$4, \$5, \".\", \"PASS\", \".\" }' \
+    > /out/_background.tsv
+
+  echo '==> diagnosis: a ClinVar pathogenic variant in $GENE' >&2
+  bcftools view -H -r 22 -i 'INFO/CLNSIG ~ \"Pathogenic\" && INFO/GENEINFO ~ \"${GENE}:\" && INFO/CLNREVSTAT ~ \"multiple_submitters\"' '$CLINVAR' \
+    | awk 'length(\$4) < 20 && length(\$5) < 20' \
+    | head -n 1 \
+    | awk -v OFS='\t' '{ print \$1, \$2, \".\", \$4, \$5, \".\", \"PASS\", \".\" }' \
+    > /out/_spike.tsv
+
+  test -s /out/_background.tsv || { echo 'no background variants found' >&2; exit 1; }
+  test -s /out/_spike.tsv || { echo 'no 2-star pathogenic ClinVar variant found for $GENE' >&2; exit 1; }
+
+  {
+    printf '##fileformat=VCFv4.2\n'
+    printf '##source=rarelens SIMULATED proband: GIAB HG002 background + one ClinVar pathogenic %s variant\n' '$GENE'
+    printf '##contig=<ID=22,length=50818468>\n'
+    printf '#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n'
+    sort -k2,2n /out/_background.tsv /out/_spike.tsv
+  } | bgzip > /out/proband-simulated.vcf.gz
+  tabix -f -p vcf /out/proband-simulated.vcf.gz
+  rm -f /out/_background.tsv
+"
+echo
+echo "wrote $OUT_DIR/proband-simulated.vcf.gz"
+echo "the planted diagnosis (chrom pos ref alt):"
+awk '{print "  " $1, $2, $4, $5}' "$OUT_DIR/_spike.tsv"
+rm -f "$OUT_DIR/_spike.tsv"