fix(science): stop scoring evidence that was never looked up

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.

ml/rarelens_ml/benchmark.py

@@ -0,0 +1,134 @@
+"""Measure the phenotype ranking against every published case in Phenopacket Store.
+
+One demo case ranking correctly is an anecdote. This asks the only question that matters for a
+phenotype-driven tool: given a real patient's reported terms, where does the gene the authors
+actually diagnosed come in a ranking of every gene HPO annotates?
+
+    python -m rarelens_ml.benchmark --phenopackets all_phenopackets.zip
+
+**Read the result with the contamination in mind.** HPO's gene-to-phenotype annotations are
+themselves curated from published case reports — quite possibly the very ones being scored here.
+The median causal gene already carries every one of its patient's terms, so this measures how well
+the ranking retrieves a gene HPO has already been told about. It is an upper bound. A prospective
+number, on a patient whose disease gene nobody has annotated yet, would be lower; how much lower
+this corpus cannot say.
+
+Ties are the other trap. Scoring by term overlap alone puts many genes on identical scores, so the
+honest report is a range: optimistic counts a tie as a win, pessimistic counts every tied gene as
+ranked ahead of the right answer. The truth is between them.
+"""
+import argparse
+import json
+import statistics
+import sys
+import zipfile
+from collections import defaultdict
+from collections.abc import Iterator
+
+from rarelens_ml.hpo import information_content, phenotype_score
+
+
+def gene_annotations(path: str) -> dict[str, set[str]]:
+    """gene -> HPO terms, from the propagated table scripts/load-hpo.py writes (TSV export)."""
+    genes: dict[str, set[str]] = defaultdict(set)
+    with open(path) as fh:
+        for line in fh:
+            gene, _, term = line.rstrip("\n").partition("\t")
+            if gene and term:
+                genes[gene].add(term)
+    return genes
+
+
+def cases(path: str) -> Iterator[tuple[str, list[str]]]:
+    """(causal gene, observed HPO terms) for each phenopacket with exactly one causal gene."""
+    with zipfile.ZipFile(path) as z:
+        for entry in z.namelist():
+            if not entry.endswith(".json"):
+                continue
+            try:
+                packet = json.loads(z.read(entry))
+            except ValueError:
+                continue
+            causal = {
+                g.get("variantInterpretation", {})
+                .get("variationDescriptor", {})
+                .get("geneContext", {})
+                .get("symbol")
+                for i in packet.get("interpretations", [])
+                for g in i.get("diagnosis", {}).get("genomicInterpretations", [])
+            } - {None}
+            terms = [
+                f["type"]["id"]
+                for f in packet.get("phenotypicFeatures", [])
+                if not f.get("excluded")
+            ]
+            if len(causal) == 1 and terms:
+                yield causal.pop(), terms
+
+
+def rank_of(gene: str, terms: list[str], genes: dict[str, set[str]], ic: dict[str, float],
+            default: float) -> tuple[int, int, float]:
+    """(optimistic rank, pessimistic rank, the causal gene's own score)."""
+    target = phenotype_score(terms, genes[gene], ic, default)
+    better = tied = 0
+    for other, annotated in genes.items():
+        if other == gene:
+            continue
+        value = phenotype_score(terms, annotated, ic, default)
+        if value > target:
+            better += 1
+        elif value == target:
+            tied += 1
+    return better + 1, better + tied + 1, target
+
+
+def report(ranks: list[tuple[int, int, float]], n_genes: int, min_terms: int = 1) -> str:
+    n = len(ranks)
+    median_score = statistics.median(t for *_, t in ranks)
+    lines = [
+        (
+            f"{n} published cases with at least {min_terms} HPO term(s), ranked against "
+            f"{n_genes} genes (random top-1 would be {1 / n_genes:.2%})"
+        ),
+        (
+            f"the causal gene's own phenotype score: median {median_score:.2f}"
+            "  <- 1.00 means HPO already carries every one of the patient's terms for that gene"
+        ),
+    ]
+    for label, column in (("optimistic", 0), ("pessimistic", 1)):
+        r = [row[column] for row in ranks]
+        lines.append(
+            f"  {label:12s} top-1 {sum(x == 1 for x in r) / n:6.1%}   "
+            f"top-10 {sum(x <= 10 for x in r) / n:6.1%}   "
+            f"MRR {sum(1 / x for x in r) / n:.3f}   median rank {statistics.median(r):.0f}"
+        )
+    return "\n".join(lines)
+
+
+def main() -> None:
+    p = argparse.ArgumentParser()
+    p.add_argument("--phenopackets", required=True,
+                   help="all_phenopackets.zip from a phenopacket-store release")
+    p.add_argument("--annotations", required=True,
+                   help="TSV of gene<tab>hpo_id, exported from the gene_phenotypes table")
+    p.add_argument("--limit", type=int, help="benchmark only the first N cases (a smoke run)")
+    p.add_argument("--min-terms", type=int, default=1,
+                   help="skip cases with fewer HPO terms; a one-term case can only tie")
+    a = p.parse_args()
+
+    genes = gene_annotations(a.annotations)
+    ic = information_content(genes)
+    default = max(ic.values(), default=1.0)
+    ranks = []
+    for gene, terms in cases(a.phenopackets):
+        if gene in genes and len(terms) >= a.min_terms:
+            ranks.append(rank_of(gene, terms, genes, ic, default))
+        if a.limit and len(ranks) >= a.limit:
+            break
+    if not ranks:
+        sys.exit("no benchmarkable cases: is --annotations the gene_phenotypes export?")
+    print(report(ranks, len(genes), a.min_terms))
+
+
+if __name__ == "__main__":
+    main()

ml/rarelens_ml/features.py

@@ -2,11 +2,23 @@
 
 Training imports it, and train.log_and_register ships this package inside the logged pyfunc
 (code_paths), so serving runs exactly this code on the raw columns below.
+
+**Allele frequency is deliberately not a feature.** It used to be, and it dominated everything:
+the same missense variant scored 0.887 at AF 0 and 0.0003 at AF 0.01, so the model was largely a
+frequency lookup. That caused two problems. It double-counted, because the ranking already scores
+frequency explicitly and auditably in `triage.rarity_score`, putting ~45% of the rank on one
+measurement. And it was circular, because ClinVar's labels are assigned with ACMG criteria that
+call a variant benign *on frequency* (BA1/BS1), so the model was rediscovering the rule used to
+label its own training data — which is most of why the headline AUROC looked so good.
+
+What is left is the variant's predicted effect: what it does to the protein, and how damaging two
+independent predictors think that is. That is evidence the rest of the ranking does not already
+have, which is the only reason to give the model a weight at all.
 """
 import pandas as pd
 
 # What serving must send: raw values as stored in the variants table / its annotations.
-RAW_COLUMNS = ["impact", "consequence", "gnomad_af", "cadd_phred", "am_pathogenicity"]
+RAW_COLUMNS = ["impact", "consequence", "cadd_phred", "am_pathogenicity"]
 
 IMPACT_ORDER = {"MODIFIER": 0, "LOW": 1, "MODERATE": 2, "HIGH": 3}
 
@@ -14,8 +26,8 @@ IMPACT_ORDER = {"MODIFIER": 0, "LOW": 1, "MODERATE": 2, "HIGH": 3}
 def build(df: pd.DataFrame) -> pd.DataFrame:
     out = pd.DataFrame(index=df.index)
     out["impact_rank"] = df["impact"].map(IMPACT_ORDER).fillna(0).astype(int)
-    # No gnomAD record means the variant was not observed: treat as AF 0.
-    out["gnomad_af"] = pd.to_numeric(df["gnomad_af"], errors="coerce").fillna(0.0)
+    # Left as NaN on purpose: LightGBM handles missing natively, and imputing a number here would
+    # assert a score nobody computed.
     out["cadd_phred"] = pd.to_numeric(df["cadd_phred"], errors="coerce")
     out["am_pathogenicity"] = pd.to_numeric(df["am_pathogenicity"], errors="coerce")
     out["consequence"] = df["consequence"].astype("category")

ml/rarelens_ml/hpo.py

@@ -0,0 +1,125 @@
+"""HPO ontology handling: propagation and information content.
+
+Shared by scripts/load-hpo.py, which writes the tables the API ranks against, and
+rarelens_ml.benchmark, which scores that ranking. It lives in the package rather than in the
+script so the arithmetic underneath the project's main scientific claim is covered by tests.
+
+Two ideas, both standard practice and both absent from the first version of the ranking:
+
+**Propagation.** HPO's gene annotations are direct. A gene linked to "Aortic root aneurysm" is not
+also linked to "Aortic aneurysm", so matching case terms by exact ID missed any patient whose
+description sat one level away from the curator's chosen term. The annotation propagation rule
+says a gene annotated with a term is annotated with all of that term's ancestors; matching then
+works in both directions without the ranking knowing the ontology exists.
+
+**Information content.** IC(term) = -ln(share of genes carrying it). After propagation almost
+every gene carries "Abnormality of the cardiovascular system", so its IC is near zero, while
+"Dilated left subclavian artery" is worth a great deal. Counting terms alike let a patient's
+"Global developmental delay" count as much as a near-pathognomonic sign.
+"""
+import io
+import math
+from collections import defaultdict
+from collections.abc import Iterable
+
+# Terms outside this branch (inheritance, clinical modifiers, frequency) describe how a disease
+# behaves rather than what is wrong with the patient, and must not count towards a match.
+PHENOTYPIC_ABNORMALITY = "HP:0000118"
+
+
+def parse_obo(handle: io.TextIOBase) -> tuple[dict[str, set[str]], dict[str, str]]:
+    """Each term's direct parents and its name, from hp.obo. Obsolete terms are dropped."""
+    parents: dict[str, set[str]] = {}
+    names: dict[str, str] = {}
+    term_id: str | None = None
+    name: str | None = None
+    is_a: set[str] = set()
+    obsolete = in_term = False
+
+    def flush() -> None:
+        if term_id and not obsolete:
+            parents[term_id] = is_a
+            names[term_id] = name or term_id
+
+    for raw in handle:
+        line = raw.rstrip("\n")
+        if line.startswith("["):
+            flush()
+            term_id, name, is_a, obsolete = None, None, set(), False
+            in_term = line == "[Term]"
+        elif not in_term:
+            continue
+        elif line.startswith("id: HP:"):
+            term_id = line[4:].strip()
+        elif line.startswith("name: "):
+            name = line[6:].strip()[:200]
+        elif line.startswith("is_a: HP:"):
+            is_a.add(line[6:].split("!")[0].strip())
+        elif line.startswith("is_obsolete: true"):
+            obsolete = True
+    flush()
+    return parents, names
+
+
+def ancestors_of(parents: dict[str, set[str]]) -> dict[str, set[str]]:
+    """Every term's ancestors, itself included.
+
+    Iterative, because HPO is deep enough to exhaust the recursion limit, and in true post-order:
+    a term is resolved only once every parent is resolved. A pre-order walk read backwards looks
+    like it would do, but on a DAG a term can be visited before one of its parents on another
+    branch, and then it silently inherits that parent alone instead of the parent's whole
+    lineage. That dropped Camptodactyly and Chiari malformation out of the phenotype branch
+    entirely, which is what this shape of bug looks like from the outside.
+    """
+    cache: dict[str, set[str]] = {}
+    for start in parents:
+        if start in cache:
+            continue
+        stack: list[tuple[str, bool]] = [(start, False)]
+        while stack:
+            node, resolved = stack.pop()
+            if node in cache:
+                continue
+            if resolved:
+                found = {node}
+                for parent in parents.get(node, ()):
+                    found |= cache.get(parent, {parent})  # fallback guards against a cycle
+                cache[node] = found
+            else:
+                stack.append((node, True))
+                stack.extend((p, False) for p in parents.get(node, ()) if p not in cache)
+    return cache
+
+
+def propagate(
+    direct: Iterable[tuple[str, str]], ancestors: dict[str, set[str]]
+) -> dict[str, set[str]]:
+    """gene -> its annotated terms plus all their ancestors, within the phenotype branch."""
+    genes: dict[str, set[str]] = defaultdict(set)
+    for gene, term in direct:
+        for node in ancestors.get(term, {term}):
+            if node != PHENOTYPIC_ABNORMALITY and PHENOTYPIC_ABNORMALITY in ancestors.get(node, ()):
+                genes[gene].add(node)
+    return dict(genes)
+
+
+def information_content(genes: dict[str, set[str]]) -> dict[str, float]:
+    """-ln(share of genes carrying the term); 0 for a term every gene has."""
+    if not genes:
+        return {}
+    counts: dict[str, int] = defaultdict(int)
+    for terms in genes.values():
+        for term in terms:
+            counts[term] += 1
+    return {term: -math.log(n / len(genes)) for term, n in counts.items()}
+
+
+def phenotype_score(
+    case_terms: Iterable[str], gene_terms: set[str], ic: dict[str, float], default: float
+) -> float:
+    """Information-content-weighted recall; the same arithmetic as app.services.triage."""
+    terms = list(case_terms)
+    total = sum(ic.get(t, default) for t in terms)
+    if total <= 0:
+        return 0.0
+    return sum(ic.get(t, default) for t in terms if t in gene_terms) / total

ml/rarelens_ml/train.py

@@ -29,10 +29,7 @@ PARAMS = {
 POS = {"pathogenic", "likely_pathogenic"}
 NEG = {"benign", "likely_benign"}
 # VEP --tab column -> raw feature column (am_pathogenicity already matches).
-VEP_TO_RAW = {
-    "IMPACT": "impact", "Consequence": "consequence", "gnomADe_AF": "gnomad_af",
-    "CADD_PHRED": "cadd_phred",
-}
+VEP_TO_RAW = {"IMPACT": "impact", "Consequence": "consequence", "CADD_PHRED": "cadd_phred"}
 
 
 def label(clin_sig: object) -> int | None: