pab1s
2026.09.14 13 min read

From five million variants to ten candidates

Part two of a series on the 2026 MVA hackathon: what the prioritization pipeline does to the VCF, one invented row followed through every gate, and why the filters set the ceiling.

· Pablo Olivares

The problem fits in one sentence: I have a table with 5,012,204 rows, and at most one of them is the reason a child is sick. Find it.

If you have built a search engine or a recommender, the shape of that task should already feel familiar. By the end of this post I want it to feel like work you have done before, with different column names. One note on scope first: this series covers only what the hackathon organizers published and methods that are already public, nothing about the specific child and nothing about variants I have found. The disease background, the data formats and the scoring live in the first post, so I skip them here.

Where the table comes from

A sequencer cannot read a 3.1-billion-letter genome in one go. It reads millions of short fragments, around 150 letters each, every letter carrying its own error estimate, and it has no idea where any fragment belongs. An upstream pipeline fixes that, passing from the full raw FASTQ to the VCF in five steps: check each fragment’s quality and trim the junk, place every fragment on the reference genome, remove duplicates and correct quality scores that drifted for boring instrument reasons, propose every position where this person differs from the reference, and throw out the proposals that look like machine error rather than biology. What survives is a single table, the VCF, with one row per real difference between this person and the reference. The hackathon hands us that table already built, and everything below starts from it.

What the pipeline does with those rows is the rest of this post (Fig. 1):

From a variant table to a ranked shortlist Two inputs, VCF and HPO terms, flow through annotation, frequency and inheritance filtering, separate deleteriousness and phenotype scoring, then aggregation to genes and one fused ordering. VCF: 5,012,204 rows variants plus caller diagnostics HPO terms coded clinical findings annotate VEP/Jannovar, gnomAD, dbNSFP filter: 5M to hundreds frequency, consequence, inheritance damage scores CADD, REVEL, splice phenotype scores Exomiser, LIRICAL, Phen2Gene aggregate, fuse, order: top 10
Fig. 1: the interpretation pipeline. Annotation adds columns, filtering removes rows that cannot explain rare disease, two scorers judge damage and fit separately, and fusion produces one order with one number per candidate.

Making the table smaller

Nobody scores five million variants one at a time. The first move is to throw most of the table away with rules that say, in effect, that a given row cannot be the cause of a rare disease. As we walk the stages I’ll carry one invented row through them: a heterozygous missense change in a placeholder gene, visible in every figure from Fig. 2 onward. Every value attached to it is invented for illustration, and nothing there comes from the hackathon case.

Before the cutting starts, the rows get wider. Each variant is looked up in a few outside catalogs: VEP or Jannovar for what it does to nearby genes, gnomAD for how common it is in a large population of presumably healthy people, dbNSFP for a bundle of precomputed pathogenicity predictions, and RefSeq or Ensembl for the gene models underneath all of it. In database terms this is a join, and in ML terms it is feature enrichment keyed on position and allele: rows gain columns and none disappear. That is exactly what happens to the example row in Fig. 2, which picks up its consequence, its population frequency, and two pathogenicity scores without moving an inch.

The example row gains four annotation columns A row, chr12:1,000,000 A greater than G, genotype 0/1, points at four added columns: missense, allele frequency 0.004 percent, CADD 28, REVEL 0.91. annotation adds columns, removes nothing chr12:1,000,000 A>G · 0/1 missense AF 0.004% CADD 28 REVEL 0.91
Fig. 2: the example row, annotated. Four columns arrive from the catalogs. The row is neither moved nor removed.

Three rules then do the cutting, and no model is involved in any of them:

  • Frequency drops anything more than about 1% of healthy people carry, since a variant that common cannot be causing a rare disease. In Fig. 3 the example row survives that gate at 0.004%, while the neighbour below it, carried by 12% of healthy people, is dropped.

    The frequency gate keeps the rare row, drops the common one Two rows. The example row at allele frequency 0.004 percent is kept. A neighbour at 12 percent is dropped. chr12:1,000,000 A>G · gnomAD AF 0.004% kept chr4:222,222,222 C>T · gnomAD AF 12% dropped
    Fig. 3: the frequency gate. The example row is rare enough to pass. The common neighbour goes, and does not come back.

  • Consequence keeps only variants that change a protein or disturb splicing, because the rest rarely explain symptoms. In Fig. 4 that keeps our missense row and drops the synonymous neighbour sitting beside it.

    The consequence gate keeps the missense row, drops the synonymous one Two rows. The example row is missense and is kept. A synonymous neighbour is dropped. chr12:1,000,000 A>G · missense kept chr12:1,000,200 C>T · synonymous dropped
    Fig. 4: the consequence gate. Missense survives. The synonymous row beside it rarely explains symptoms, so it goes.

  • Inheritance asks whether the variant fits the way the disease runs in this family, which includes the compound heterozygous case: two different damaging variants in the same gene, one inherited from each parent. For the example row a second damaging hit in the same gene arrives from the other parent, so both rows fit the pattern and stay, as Fig. 5 shows. A lone damaging hit with nothing behind it does not.

    The inheritance gate keeps the compound heterozygous pair The example row, compound heterozygous with a second hit in the same gene in trans, is kept. A damaging row with no second hit is dropped. chr12:1,000,000 A>G · compound het with chr12:1,000,480 kept chr7:888,888,888 A>C · damaging, but no second hit dropped
    Fig. 5: the inheritance gate. Two damaging hits in one gene, one from each parent, fit the family pattern. A lone hit does not.

In rough numbers, frequency takes the 5,012,204 rows down to the order of 10510^5, consequence takes that to the order of 10410^4, and inheritance leaves a few hundred.

What deserves attention is what this stage is. It is the cheap first stage of a two-stage recommender, implemented as deterministic rules rather than learned, and every row it drops is gone for good. The dashed rows in Figs. 3 to 5 do not come back at the scoring stage or at any other point, so the recall of this stage is a ceiling on everything after it. If the causal variant is dropped here, nothing downstream brings it back, however clever the scorer turns out to be. That single fact shapes the whole analysis more than any choice made later, which is why the auditing effort belongs on the filter thresholds rather than on the choice of scorer.

Judging what survives

A few hundred candidates remain, and each one can now be examined properly. Two different questions get asked of them.

The first is whether the variant is damaging on its own molecular merits. CADD, REVEL, SIFT and PolyPhen, together with various splice predictors, answer exactly that, one variant at a time, as frozen pointwise features: pretrained somewhere else, applied here, with the usual transfer risk that implies. None of them has ever heard of this patient.

The second is whether the gene matches what the clinicians observe. The symptoms have been coded as HPO terms, the Human Phenotype Ontology, and three public tools compare them against what is known about genes and diseases. This is retrieval and ranking: the symptoms are the query, genes and diseases are the corpus, and the three tools are three rankers over it. Exomiser measures the similarity between the patient’s symptoms and known gene-to-phenotype relationships in humans and in model organisms, using the structure of the ontology, and returns one combined score. LIRICAL treats every symptom as a diagnostic test: each symptom earns a likelihood ratio, and the ratios are multiplied together with a genotype-based ratio mixed in, symptoms assumed independent. The result is a posttest probability per candidate disease, with each symptom’s contribution still visible. Phen2Gene looks the symptoms up in a precomputed gene index that weights specific findings above vague ones (“seizures” tells you little, a rare specific finding tells you a lot), roughly the role inverse document frequency plays in search. It answers in about a second, leaving variant filtering to a separate companion step.

The example row scores well on both questions, as Fig. 6 puts side by side: the damage predictors put it near the top of their range, and the gene it sits in matches the coded symptoms closely enough to matter. The two panels judge the row independently, which is exactly what makes the fusion step afterwards non-trivial.

The example row under both scorers Two panels. Damage: CADD 28 and REVEL 0.91, both high. Phenotype fit: 0.86, the gene matches the coded symptoms. is it damaging? CADD 28 · REVEL 0.91 both score it high does the gene match? phenotype fit 0.86 the gene matches
Fig. 6: the two scores. Damage and phenotype fit judge the example row independently, and neither knows about the other.

Merging into one ranking

The last move is fusion. Variant-level scores roll up to genes and diseases, combinations the inheritance pattern cannot support are pruned, and damage is merged with phenotype fit into one number and one rank per candidate. The output is the top ten. For the example row the two scores merge to 0.88, which lands it at rank two, the highlighted slot in Fig. 7.

The example row fused to one number, ranked second Damage and phenotype fit merge to 0.88, and the row lands at rank two of the final top ten, shown as ten slots with the second highlighted. damage and phenotype fit merge to 0.88 the final top ten rank 2
Fig. 7: fusion. Two scores become one number, and the example row lands at rank two of the final top ten.

The same path as Figs. 2 to 7, compressed into a few lines:

# VCF + HPO terms -> top 10
annotated = annotate(vcf)                 # VEP/Jannovar, gnomAD, dbNSFP: add columns
rare      = filter_frequency(annotated)   # common in healthy people -> out
coding    = filter_consequence(rare)      # protein-altering or splice only
fitting   = filter_inheritance(coding)    # family pattern, incl. compound het
damage    = score_damage(fitting)         # CADD, REVEL: fixed, pretrained functions
relevance = score_phenotype(fitting, hpo) # Exomiser / LIRICAL / Phen2Gene
ranked    = fuse(damage, relevance)       # roll up to gene, merge, calibrate
return top(ranked, 10)

Where fusion gets subtle

The final merge has to satisfy two metrics at once, and they want different things. Rank Points cares only about the order of the top ten. F-max thresholds the raw scores across the whole list. Take a monotone transform of the scores and the top-ten order is untouched, so the rank-2 slot in Fig. 7 stays exactly where it is and Rank Points does not move at all, while the rescaled scores can put the F-max threshold in a very different place. That is why LIRICAL’s calibrated probabilities and a rescaled Exomiser score can produce the same top ten and still behave completely differently under F-max. The last step is therefore a ranking problem and a calibration problem at the same time, and the pipeline has to get both right.

Put together, the shape is the familiar one: a brutal rules-based first stage whose recall caps everything after it, then a ranker that blends frozen pretrained scores with a retrieval-style match to the symptoms, evaluated on ordering and on calibration at once. That is variant prioritization, and it is the system I am building on.