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.
244 lines
17 KiB
Markdown
244 lines
17 KiB
Markdown
# Data: what rarelens actually runs on
|
||
|
||
Everything below is public, peer-reviewed and consented for open redistribution. No patient data,
|
||
no data access agreement, nothing that needs an application. These are the references to quote
|
||
when showing the platform to someone.
|
||
|
||
Citations were verified against [PubMed](https://pubmed.ncbi.nlm.nih.gov/); each row links its DOI.
|
||
|
||
## The demo slice
|
||
|
||
`make data` fetches two real files, chromosome 22 only (roughly 100 MB, minutes rather than hours):
|
||
|
||
| File | What it is | Role |
|
||
|---|---|---|
|
||
| `data/example.vcf.gz` | GIAB HG002 (NA24385) v4.2.1 benchmark calls, GRCh38, chr22 | the sample a scientist annotates |
|
||
| `data/clinvar.chr22.vcf.gz` | ClinVar, GRCh38, chr22 | training labels, and the ClinVar column in the UI |
|
||
|
||
HG002 is the NIST Genome in a Bottle Ashkenazi son, recruited through the Personal Genome Project,
|
||
which consents participants to unrestricted public release. It is the reference genome the field
|
||
benchmarks variant callers against, so it is both realistic and unambiguously shareable.
|
||
|
||
## Datasets
|
||
|
||
| Dataset | Used for | Access | Terms | Citation |
|
||
|---|---|---|---|---|
|
||
| **ClinVar** (GRCh38) | pathogenic/benign labels, ClinVar column | `ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh38/` | NCBI public domain | Landrum et al., *Nucleic Acids Res* 48(D1):D835–D844, 2020. [10.1093/nar/gkz972](https://doi.org/10.1093/nar/gkz972) |
|
||
| **Genome in a Bottle** HG002 v4.2.1 | the demo sample | `ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/`, `s3://giab` | open, no use restriction | Zook et al., *Nat Biotechnol* 37:561–566, 2019. [10.1038/s41587-019-0074-6](https://doi.org/10.1038/s41587-019-0074-6) |
|
||
| **gnomAD** v4 | allele frequency feature and filter | `gs://gcp-public-data--gnomad`, `s3://gnomad-public-us-east-1` | free use, no restriction | Chen et al., *Nature* 625:92–100, 2024. [10.1038/s41586-023-06045-0](https://doi.org/10.1038/s41586-023-06045-0); Karczewski et al., *Nature* 581:434–443, 2020. [10.1038/s41586-020-2308-7](https://doi.org/10.1038/s41586-020-2308-7) |
|
||
| **1000 Genomes** 30x | optional cohort/trio data | EBI FTP, `s3://1000genomes` | fully open, no access restriction | Byrska-Bishop et al., *Cell* 185(18):3426–3440.e19, 2022. [10.1016/j.cell.2022.08.004](https://doi.org/10.1016/j.cell.2022.08.004) |
|
||
| **MANE Select** | one transcript per gene, if transcript choice ever matters | Ensembl/RefSeq | open | Morales et al., *Nature* 604:310–315, 2022. [10.1038/s41586-022-04558-8](https://doi.org/10.1038/s41586-022-04558-8) |
|
||
| **Human Phenotype Ontology** gene-to-phenotype | what the phenotype half of the ranking matches against (`make hpo`) | `purl.obolibrary.org/obo/hp/hpoa/genes_to_phenotype.txt` | free to use with attribution | Gargano et al., *Nucleic Acids Res* 52(D1):D1333–D1346, 2024. [10.1093/nar/gkad1005](https://doi.org/10.1093/nar/gkad1005) |
|
||
| **Phenopacket Store** | published patients: the reported phenotype and causal variant of a real case (`make published-case`) | `github.com/monarch-initiative/phenopacket-store` | BSD-3-Clause | Danis et al., *HGG Adv* 6(1):100371, 2025. [10.1016/j.xhgg.2024.100371](https://doi.org/10.1016/j.xhgg.2024.100371) |
|
||
|
||
## Tools and scores
|
||
|
||
| Tool | Role | Terms | Citation |
|
||
|---|---|---|---|
|
||
| **Ensembl VEP** 113 | annotation (`pipeline/modules/vep.nf`) | Apache 2.0 | McLaren et al., *Genome Biol* 17:122, 2016. [10.1186/s13059-016-0974-4](https://doi.org/10.1186/s13059-016-0974-4) |
|
||
| **CADD** | `cadd_phred` feature | free for non-commercial use; commercial licence required | Rentzsch et al., *Nucleic Acids Res* 47(D1):D886–D894, 2019. [10.1093/nar/gky1016](https://doi.org/10.1093/nar/gky1016); Schubach et al., *Nucleic Acids Res* 52(D1), 2024. [10.1093/nar/gkad989](https://doi.org/10.1093/nar/gkad989) |
|
||
| **AlphaMissense** | `am_pathogenicity` feature | predictions moved to CC BY 4.0 in March 2024 (originally CC BY-NC-SA) | Cheng et al., *Science* 381:eadg7492, 2023. [10.1126/science.adg7492](https://doi.org/10.1126/science.adg7492) |
|
||
|
||
Neither score is required: `rarelens_ml.features` treats a missing CADD or AlphaMissense value as
|
||
NaN and LightGBM handles it, so the pipeline runs without the plugin data.
|
||
|
||
## The simulated proband
|
||
|
||
`make demo-case` builds `data/proband-simulated.vcf.gz`: real GIAB HG002 variants as background
|
||
plus one real ClinVar 2-star pathogenic variant in a disease gene (NF2 by default, giving
|
||
neurofibromatosis type 2). Spiking a known variant into a public genome is how phenotype-driven
|
||
triage tools are benchmarked, every input is public, and the file's header says SIMULATED. It is
|
||
not a patient, and no part of it is invented: both the background and the planted variant are real
|
||
published records.
|
||
|
||
The point of it is that the case has a right answer, so the ranking can be checked rather than
|
||
admired. Give the case the phenotype of the planted disease and the planted variant should rank
|
||
first — on phenotype, rarity and consequence, with ClinVar agreeing only afterwards.
|
||
|
||
**Caveat when running without a VEP cache.** `VEP_DATABASE=true` queries Ensembl's public database
|
||
instead of the 25 GB cache. It returns no gnomAD frequencies and no plugin scores, so the rarity
|
||
and model components have nothing to work with. They abstain: the job records what the run looked
|
||
up (`jobs.has_frequencies`, `jobs.has_effect_scores`), those components are dropped from the score
|
||
rather than given a default, and the remaining weights are renormalised. Fine for showing the
|
||
mechanics; use the cache for anything you would quote.
|
||
|
||
## A published case
|
||
|
||
`make published-case` builds a case around a patient who actually exists in the literature. It
|
||
reads a GA4GH phenopacket from Monarch's Phenopacket Store, which curates published case reports
|
||
into machine-readable records and keeps the PMID on each one, and takes two things from it
|
||
verbatim: the phenotype terms the authors reported, and the variant they called causal.
|
||
|
||
The default is **Loeys-Dietz syndrome**, from the paper that first defined it — Loeys et al.,
|
||
*Nat Genet* 37:275–281, 2005, [10.1038/ng1511](https://doi.org/10.1038/ng1511) (PMID 15731757).
|
||
Family 4, individual II-1: 30 reported HPO terms, from hypertelorism and a bifid uvula to arterial
|
||
tortuosity and an aortic root aneurysm, and a heterozygous *TGFBR2* missense variant,
|
||
`NM_003242.6:c.1069G>T` `p.(Gly357Trp)`, at GRCh38 chr3:30672252 G>T.
|
||
|
||
The phenotype and the answer are real. The rest of that patient's genome is not public, and
|
||
triage means nothing if the causal variant is the only variant in the file, so background variants
|
||
come from GIAB HG002 in a 3 Mb window around the locus. The file is therefore a published
|
||
diagnosis inside a public background genome — not anyone's exome. That is the standard
|
||
construction for benchmarking phenotype-driven triage, and it is the reason this case can be
|
||
redistributed at all.
|
||
|
||
The background is taken from coding exons wherever possible, using coding-exon coordinates from
|
||
Ensembl's public REST API. This matters more than it sounds: of the ~4,000 HG002 variants in that
|
||
window only 9 fall in coding exons, so a random sample is entirely intronic, the consequence
|
||
filter throws all of it away, and the causal variant ends up the only candidate left — a funnel
|
||
that proves nothing. Real coding variants give the ranking something it has to rank *against*.
|
||
|
||
What that run looks like: 21 variants in, **2** surviving the consequence filter, **1** matching
|
||
the phenotype. Run without a VEP cache, so two of the four components have nothing to go on and
|
||
abstain; the remaining weights are renormalised over 0.55.
|
||
|
||
| | score | phenotype (0.64) | rarity | consequence (0.36) | model |
|
||
|---|---|---|---|---|---|
|
||
| *TGFBR2* 3:30672252 missense | **0.855** | 1.00 (30/30 terms) | not looked up | 0.60 | not looked up |
|
||
| *OSBPL10* 3:31748090 missense | 0.218 | 0.00 | not looked up | 0.60 | not looked up |
|
||
|
||
This is the whole argument for phenotype-driven triage in one table. Both are rare missense
|
||
variants, identical on every piece of evidence this run has except one. What separates the
|
||
published diagnosis from an incidental variant in a lipid-transport gene is the patient's
|
||
phenotype, and nothing else. ClinVar's "pathogenic" on the first row is shown afterwards as
|
||
independent confirmation — it is not an input to the rank.
|
||
|
||
Three things to say out loud when showing it:
|
||
|
||
- **The phenotype match is partly circular.** HPO's gene-to-phenotype annotations are themselves
|
||
curated from published cases, quite possibly including this one. A 30/30 term match against
|
||
*TGFBR2* is evidence the plumbing works, not evidence the ranking would find a novel gene.
|
||
- **Two of the four components abstained, and that is the honest outcome.** In database mode VEP
|
||
returns no frequencies and no plugin scores, so rarity and the model have nothing to say.
|
||
`--af_gnomade` is rejected outright with `--database`, and plain `--af` returns nothing even for
|
||
common variants — checked against rs429358, roughly 15% globally. An earlier version of this
|
||
table read 1.00 for rarity on both rows, which was not a measurement: it was the absence of one.
|
||
- **The background is one healthy genome, not a diagnostic exome.** A real case would have
|
||
thousands of rare coding variants to discard, not a handful.
|
||
|
||
Other cases work the same way — any phenopacket with GRCh38 coordinates will do:
|
||
|
||
```
|
||
scripts/make-published-case.py --phenopacket <raw phenopacket-store JSON URL>
|
||
```
|
||
|
||
## The model, and why it currently abstains
|
||
|
||
`make training-set` builds a training table straight from ClinVar rather than running VEP over
|
||
hundreds of thousands of variants: ClinVar already carries the molecular consequence (`MC`) and the
|
||
gene (`GENEINFO`). Only 2-star-and-above records are kept. `make train` then fits LightGBM and
|
||
points the `production` alias at the new version. 312,025 training and 74,239 held-out variants
|
||
across 7,728 and 1,932 genes, with no gene on both sides.
|
||
|
||
The model used to take allele frequency as a feature. Removing it is the single most informative
|
||
thing in this document, because of what happened to the numbers:
|
||
|
||
| | AUROC | AUPRC | missense AUROC | missense AUPRC |
|
||
|---|---|---|---|---|
|
||
| v2, with gnomAD allele frequency | 0.986 | 0.954 | 0.872 | 0.725 |
|
||
| v3, allele frequency removed | 0.966 | 0.881 | **0.500** | 0.398 |
|
||
|
||
**0.500 on missense is exactly random.** Strip frequency out and the model cannot tell one missense
|
||
variant from another at all, because nothing is left but the consequence class — every missense row
|
||
gets the identical score. So the respectable-looking 0.872 was not variant-effect knowledge. It was
|
||
allele frequency, and ClinVar's benign calls are *made with* allele frequency under ACMG's BA1/BS1
|
||
criteria. The feature had partly caused the label; the model had rediscovered the labelling rule.
|
||
|
||
The remaining 0.966 is the same trick one level up: ClinVar's pathogenic set is largely loss of
|
||
function and its benign set largely is not, so a model handed the consequence class separates them
|
||
without knowing anything hard.
|
||
|
||
Two consequences, both deliberate:
|
||
|
||
1. **Allele frequency is no longer a feature.** The ranking already scores frequency explicitly, as
|
||
a step function a reviewer can read (`triage.rarity_score`). Feeding it to the model as well put
|
||
roughly 45% of every rank on one measurement counted twice.
|
||
2. **The model abstains unless it has CADD or AlphaMissense.** Without them it can only restate the
|
||
consequence class, which the ranking already scores — and 0.500 is the measurement saying so.
|
||
Install the plugin data (see data/README.md) and the model earns its 0.20 back; until then it
|
||
contributes nothing, and the UI says so rather than showing a number.
|
||
|
||
It is also not comparable to published CADD or AlphaMissense figures, which are trained and
|
||
evaluated on different data. A fair comparison scores the same held-out rows with all three, which
|
||
needs the plugin data and is the obvious next step.
|
||
|
||
## How well does the phenotype ranking actually work?
|
||
|
||
`make benchmark` runs the ranking against every case in Phenopacket Store: given a real patient's
|
||
reported terms, where does the gene their authors diagnosed come in a ranking of all 5,269 genes
|
||
HPO annotates? Ties are reported as a range, because term-overlap scoring puts many genes on
|
||
identical scores — optimistic counts a tie as a win, pessimistic counts every tied gene as ahead.
|
||
|
||
| cases | | top-1 | top-10 | MRR |
|
||
|---|---|---|---|---|
|
||
| all 10,178 | optimistic | 81.0% | 87.2% | 0.830 |
|
||
| | pessimistic | 45.9% | 69.7% | 0.541 |
|
||
| the 6,485 with ≥6 terms | optimistic | 77.1% | 85.4% | 0.797 |
|
||
| | pessimistic | 59.5% | 81.0% | 0.670 |
|
||
|
||
Random guessing would put the right gene first 0.02% of the time, so the phenotype term is doing
|
||
real work. Two caveats, and the first is severe.
|
||
|
||
**The benchmark is contaminated.** The median causal gene already carries *every one* of its
|
||
patient's terms, because HPO's gene annotations are curated from these same case reports. This
|
||
measures how well the ranking retrieves a gene HPO has already been told about — an upper bound.
|
||
A prospective number, on a patient whose gene nobody has annotated yet, would be lower, and this
|
||
corpus cannot say by how much.
|
||
|
||
**Information content earns its place; propagation does not, measurably.** Both were added to
|
||
replace plain term counting, and the corpus was asked whether they helped. On the 6,485 cases with
|
||
at least six terms, pessimistic figures (the honest end of the range):
|
||
|
||
| scoring | top-1 | top-10 | MRR |
|
||
|---|---|---|---|
|
||
| direct annotations, count terms (the original) | 61.8% | 80.1% | 0.682 |
|
||
| direct annotations, weight by information content | **63.6%** | **83.5%** | **0.706** |
|
||
| propagated annotations, count terms | 58.2% | 77.7% | 0.653 |
|
||
| propagated + information content (shipped) | 59.5% | 81.0% | 0.670 |
|
||
|
||
Weighting by specificity helps: it breaks ties, which is exactly what it is for. Propagation costs
|
||
about as much as weighting gains, and the shipped combination is a wash against the original — a
|
||
point or two either way inside a contaminated benchmark.
|
||
|
||
Propagation is kept anyway, and the reason is worth stating plainly rather than hiding behind the
|
||
numbers. This corpus cannot show what propagation is for: its terms were chosen by the same
|
||
curators whose choices HPO records, so the IDs already line up and exact matching is flattered.
|
||
The app's users pick their own terms from a search box and will not line up that neatly. The
|
||
measurement is here so anyone who disagrees can act on it — the four rows above are one flag and
|
||
one argument to `make benchmark` apart.
|
||
|
||
## Evaluating the model honestly
|
||
|
||
The model trains on ClinVar labels and is scored on ClinVar-labelled variants, which is exactly
|
||
where published benchmarks go wrong. What to do about it:
|
||
|
||
1. **Never let the label into the features.** `CLIN_SIG` is excluded by construction; `clinvar_sig`
|
||
is stored for display only (`rarelens_ml/features.py` lists the feature columns). Allele
|
||
frequency was removed for a related reason: ACMG uses it to *assign* the benign label.
|
||
2. **Split by gene, not by variant.** Random splits put variants from the same gene on both sides,
|
||
and a model can then score a gene rather than a variant. Grimm et al. showed this inflates
|
||
reported accuracy for exactly this class of tool: *Hum Mutat* 36:513–523, 2015.
|
||
[10.1002/humu.22768](https://doi.org/10.1002/humu.22768) *Implemented*:
|
||
`rarelens_ml.train.split_by_gene` holds out whole genes.
|
||
3. **Prefer a time-based holdout.** Train on an older ClinVar release (monthly archives live under
|
||
`vcf_GRCh38/archive_2.0/`) and test only on variants classified after that date. This is the
|
||
closest thing to a prospective evaluation available without new patients.
|
||
4. **Filter labels by review status.** ClinVar's `CLNREVSTAT` marks how much evidence backs a
|
||
classification; two-star and above ("multiple submitters, no conflicts") is the usual bar.
|
||
*Known gap*: VEP's `CLIN_SIG` does not carry review status, so this needs ClinVar annotated as a
|
||
custom field before it can be enforced.
|
||
5. **Report against published baselines on the same rows.** CADD PHRED and AlphaMissense are
|
||
already columns in the variant table, so AUROC and AUPRC for the model next to those two, with
|
||
the variant count, is a fair comparison rather than a number with nothing to beat.
|
||
6. **Report AUPRC, not just AUROC.** Pathogenic variants are the minority class; AUROC flatters.
|
||
|
||
## What must not be claimed
|
||
|
||
ACMG/AMP treats computational predictions as *supporting* evidence only, never sufficient on their
|
||
own for classifying a variant (Richards et al., *Genet Med* 17:405–424, 2015.
|
||
[10.1038/gim.2015.30](https://doi.org/10.1038/gim.2015.30)). rarelens is a learning platform on
|
||
public data: it makes no diagnostic claim, and the UI shows a score next to the evidence rather
|
||
than a verdict. For what a real diagnostic pipeline looks like end to end, see the 100,000 Genomes
|
||
Project rare-disease pilot: Smedley et al., *N Engl J Med* 385:1868–1880, 2021.
|
||
[10.1056/NEJMoa2035790](https://doi.org/10.1056/NEJMoa2035790)
|