rarelens/docs/data.md

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; 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
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
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; Karczewski et al., Nature 581:434–443, 2020. 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
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
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
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

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
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; Schubach et al., Nucleic Acids Res 52(D1), 2024. 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

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 (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 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). 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