Kemal Yaylali
e76ae847a1
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.
171 lines
7.6 KiB
Python
171 lines
7.6 KiB
Python
"""Ranking is the scientific claim this app makes, so it is tested as pure logic."""
|
|
import pytest
|
|
|
|
from app.models import Prediction, Variant
|
|
from app.services import triage
|
|
|
|
# A run with a VEP cache and plugins: every line of evidence was looked up.
|
|
FULL = triage.Evidence(frequencies=True, effect_scores=True)
|
|
# VEP's database mode: no frequencies, no CADD/AlphaMissense.
|
|
DATABASE_ONLY = triage.Evidence(frequencies=False, effect_scores=False)
|
|
|
|
|
|
def variant(**kw: object) -> Variant:
|
|
fields: dict = {
|
|
"id": 1, "chrom": "22", "pos": 100, "ref": "A", "alt": "G",
|
|
"gene": "NF2", "impact": "HIGH", "consequence": "frameshift_variant",
|
|
"gnomad_af": None, "clinvar_sig": None, "annotations": {},
|
|
}
|
|
fields.update(kw)
|
|
score = fields.pop("score", None)
|
|
v = Variant(**fields)
|
|
if score is not None:
|
|
v.prediction = Prediction(model_name="m", model_version="1", score=float(score))
|
|
return v
|
|
|
|
|
|
def ontology(gene_terms: dict[str, set[str]], ic: dict[str, float] | None = None) -> triage.Ontology:
|
|
"""Equal information content unless a test is specifically about specificity."""
|
|
terms = {t for terms in gene_terms.values() for t in terms}
|
|
return triage.Ontology(gene_terms=gene_terms, ic=ic or dict.fromkeys(terms, 1.0))
|
|
|
|
|
|
def test_weights_sum_to_one() -> None:
|
|
assert sum(triage.WEIGHTS.values()) == pytest.approx(1.0)
|
|
|
|
|
|
@pytest.mark.parametrize(
|
|
("af", "expected"),
|
|
[(None, 1.0), (0.0, 1.0), (0.00005, 0.8), (0.0005, 0.5), (0.005, 0.2), (0.05, 0.0)],
|
|
)
|
|
def test_rarity_rewards_absence_from_gnomad(af: float | None, expected: float) -> None:
|
|
assert triage.rarity_score(af) == expected
|
|
|
|
|
|
@pytest.mark.parametrize(
|
|
("impact", "expected"),
|
|
[("HIGH", 1.0), ("MODERATE", 0.6), ("LOW", 0.2), ("MODIFIER", 0.0), (None, 0.0), ("?", 0.0)],
|
|
)
|
|
def test_consequence_severity(impact: str | None, expected: float) -> None:
|
|
assert triage.consequence_score(impact) == expected
|
|
|
|
|
|
def test_phenotype_match_is_the_share_of_the_patients_terms() -> None:
|
|
case_terms = ["HP:0000365", "HP:0009592", "HP:0002321", "HP:0000598"]
|
|
o = triage.Ontology(
|
|
gene_terms={"NF2": {"HP:0000365", "HP:0009592"}}, ic=dict.fromkeys(case_terms, 1.0)
|
|
)
|
|
score, matched = triage.phenotype_score("NF2", case_terms, o)
|
|
assert score == pytest.approx(0.5) # 2 of 4 terms, all equally specific
|
|
assert matched == ["HP:0000365", "HP:0009592"]
|
|
|
|
|
|
def test_a_specific_term_outweighs_a_common_one() -> None:
|
|
"""Counting terms alike let 'global developmental delay' rival a near-pathognomonic sign."""
|
|
ic = {"HP:0001263": 0.1, "HP:0000193": 6.0} # developmental delay vs bifid uvula
|
|
case_terms = ["HP:0001263", "HP:0000193"]
|
|
common = triage.phenotype_score("A", case_terms, triage.Ontology({"A": {"HP:0001263"}}, ic))
|
|
specific = triage.phenotype_score("B", case_terms, triage.Ontology({"B": {"HP:0000193"}}, ic))
|
|
assert common[0] == pytest.approx(0.1 / 6.1)
|
|
assert specific[0] == pytest.approx(6.0 / 6.1)
|
|
assert specific[0] > common[0] * 10
|
|
|
|
|
|
def test_an_unknown_term_is_treated_as_maximally_specific() -> None:
|
|
"""It can never match, so it must depress every gene equally rather than be ignored."""
|
|
o = triage.Ontology({"NF2": {"HP:0000365"}}, {"HP:0000365": triage.DEFAULT_IC})
|
|
score, _ = triage.phenotype_score("NF2", ["HP:0000365", "HP:9999999"], o)
|
|
assert score == pytest.approx(0.5)
|
|
|
|
|
|
def test_phenotype_match_is_zero_for_genes_hpo_has_never_annotated() -> None:
|
|
assert triage.phenotype_score("NOVEL1", ["HP:0000365"], ontology({})) == (0.0, [])
|
|
|
|
|
|
def test_phenotype_match_is_zero_when_no_phenotype_was_entered() -> None:
|
|
assert triage.phenotype_score("NF2", [], ontology({"NF2": {"HP:0000365"}})) == (0.0, [])
|
|
|
|
|
|
def test_the_funnel_counts_each_narrowing_step() -> None:
|
|
variants = [
|
|
variant(id=1, gnomad_af=None, impact="HIGH", gene="NF2"), # rare, coding, matched
|
|
variant(id=2, gnomad_af=0.0002, impact="MODERATE", gene="CHEK2"), # rare, coding
|
|
variant(id=3, gnomad_af=0.3, impact="HIGH", gene="NF2"), # common
|
|
variant(id=4, gnomad_af=None, impact="MODIFIER", gene="NF2"), # rare, non-coding
|
|
]
|
|
funnel = triage.funnel(variants, ["HP:0000365"], ontology({"NF2": {"HP:0000365"}}), FULL)
|
|
assert (funnel.total, funnel.rare, funnel.candidates, funnel.phenotype_matched) == (4, 3, 2, 1)
|
|
assert funnel.frequencies is True
|
|
|
|
|
|
def test_the_funnel_admits_when_the_rare_step_filtered_nothing() -> None:
|
|
variants = [variant(id=1, gnomad_af=None, impact="HIGH")]
|
|
assert triage.funnel(variants, [], ontology({}), DATABASE_ONLY).frequencies is False
|
|
|
|
|
|
def test_the_diagnosis_outranks_the_noise() -> None:
|
|
o = ontology({"NF2": {"HP:0000365", "HP:0009592"}})
|
|
case_terms = ["HP:0000365", "HP:0009592"]
|
|
diagnosis = variant(id=1, gene="NF2", impact="HIGH", gnomad_af=None, score=0.94)
|
|
plausible = variant(id=2, gene="CHEK2", impact="MODERATE", gnomad_af=0.0004, score=0.55)
|
|
noise = variant(id=3, gene="TTN", impact="MODERATE", gnomad_af=0.0009, score=0.10)
|
|
|
|
ranked = triage.rank([noise, plausible, diagnosis], case_terms, o, FULL)
|
|
assert [c.variant.id for c in ranked] == [1, 2, 3]
|
|
top = ranked[0]
|
|
assert top.matched_terms == case_terms
|
|
assert top.components["phenotype"] == 1.0
|
|
assert top.score == pytest.approx(0.35 + 0.25 + 0.20 + 0.20 * 0.94)
|
|
|
|
|
|
def test_evidence_that_was_never_looked_up_abstains_instead_of_scoring_full_marks() -> None:
|
|
"""The bug this guards: a database-mode run gave every variant rarity 1.0 for free."""
|
|
v = variant(id=1, gene="NF2", impact="HIGH", gnomad_af=None, score=0.94)
|
|
o = ontology({"NF2": {"HP:0000365"}})
|
|
[candidate] = triage.rank([v], ["HP:0000365"], o, DATABASE_ONLY)
|
|
|
|
assert candidate.components["rarity"] is None
|
|
assert candidate.components["model"] is None
|
|
# Only phenotype (0.35) and consequence (0.20) had evidence, renormalised over 0.55.
|
|
assert candidate.score == pytest.approx((0.35 * 1.0 + 0.20 * 1.0) / 0.55)
|
|
|
|
|
|
def test_the_model_abstains_when_it_has_no_feature_the_ranking_lacks() -> None:
|
|
"""Without CADD or AlphaMissense the model only restates the consequence class."""
|
|
v = variant(id=1, gnomad_af=None, score=0.89)
|
|
evidence = triage.Evidence(frequencies=True, effect_scores=False)
|
|
[candidate] = triage.rank([v], [], ontology({}), evidence)
|
|
assert candidate.components["model"] is None
|
|
assert candidate.scored is True # a prediction exists; it just does not earn a weight
|
|
|
|
|
|
def test_weights_in_use_renormalise_to_one() -> None:
|
|
full = triage.weights_in_use(FULL)
|
|
assert sum(full.values()) == pytest.approx(1.0)
|
|
assert full == {name: pytest.approx(w) for name, w in triage.WEIGHTS.items()}
|
|
|
|
partial = triage.weights_in_use(DATABASE_ONLY)
|
|
assert set(partial) == {"phenotype", "consequence"}
|
|
assert sum(partial.values()) == pytest.approx(1.0)
|
|
assert partial["phenotype"] == pytest.approx(0.35 / 0.55, abs=1e-4)
|
|
|
|
|
|
def test_an_unscored_variant_still_ranks_and_says_so() -> None:
|
|
[candidate] = triage.rank([variant(id=1, gnomad_af=None)], [], ontology({}), FULL)
|
|
assert candidate.components["model"] is None
|
|
assert candidate.scored is False
|
|
|
|
|
|
def test_common_and_non_coding_variants_are_not_candidates() -> None:
|
|
variants = [
|
|
variant(id=1, gnomad_af=0.2, impact="HIGH"),
|
|
variant(id=2, gnomad_af=None, impact="MODIFIER"),
|
|
]
|
|
assert triage.rank(variants, [], ontology({}), FULL) == []
|
|
|
|
|
|
def test_ranking_is_deterministic_for_equal_scores() -> None:
|
|
a = variant(id=7, gene="AAA", chrom="1", pos=10, gnomad_af=None)
|
|
b = variant(id=3, gene="BBB", chrom="1", pos=10, gnomad_af=None)
|
|
assert [c.variant.id for c in triage.rank([a, b], [], ontology({}), FULL)] == [3, 7]
|