rarelens: building a rare-disease variant triage platform, and the four things I got wrong

rarelens takes a patient’s genome and their clinical phenotype and narrows thousands of variants to a handful a scientist can actually review, showing the evidence behind every rank. It is a working monorepo: a Nextflow DSL2 pipeline that annotates variants with Ensembl VEP, a FastAPI service that ranks them against Human Phenotype Ontology annotations, a LightGBM model served from an MLflow registry, a SvelteKit interface built around triage decisions rather than table filtering, and Terraform for two deployment tracks on Google Cloud — a serverless one that idles at roughly £1 a month, and a Kubernetes one with Argo Workflows, Argo Events and ArgoCD behind a feature flag.

Thank you for reading this post, don't forget to subscribe!

https://gitea.yaylali.uk/kemal/rarelens

Everything runs on public, openly licensed data. No patient data is used, accepted, or possible to load.

This is a technical tour. It is also, deliberately, an account of the things that turned out to be wrong — including a model metric that collapsed from 0.872 to exactly 0.500 the moment I took away the feature that was doing all the work, and an ontology bug that silently deleted 399 terms.

The problem

A rare-disease proband’s exome contains something like twenty thousand coding variants. Perhaps a few hundred are rare and protein-altering. Exactly one, usually, explains why the patient is ill.

Finding it is not a filtering problem so much as an evidence problem: the same variant is uninteresting in one patient and diagnostic in another, and what changes between them is the clinical phenotype. That observation drives the whole design. The interface is not a variant table with filters — it is a ranked shortlist where each candidate carries the four pieces of evidence that put it there, and where a reviewer shortlists or dismisses with a reason that ends up in a case report.

Architecture diagram of the rarelens system — One repository holds the pipeline, the API, the interface, model serving and the infrastructure for two deployment tracks. The scientist never sees any of it except the middle column.

The Nextflow pipeline

Annotation is a three-process DSL2 workflow, each process a pinned container: bcftools normalises, Ensembl VEP 113 annotates, and a loader writes results into PostgreSQL and marks the job succeeded.

The Nextflow DSL2 pipeline — The same workflow file runs on a laptop, on Google Batch and under Argo Workflows; the executor is a profile, not a rewrite.

Variant identity has to survive annotation

VEP’s Location and Allele output columns trim indel alleles and shift positions by one, which is correct for VEP’s own purposes and wrong as a primary key. A deletion written by the caller as CT>C comes back as -> at a different coordinate, and the loader then stores a variant that does not exist in the input file.

The fix is to carry identity through a field VEP does not touch. The normalise step sets the VCF ID to CHROM_POS_REF_ALT, and the loader parses identity back out of it:

bcftools norm -m -any -f $genome | \
  bcftools annotate --set-id '%CHROM\_%POS\_%REF\_%FIRST_ALT'

Verified on a real run: 22:42126611 CT>C round-trips with its alleles intact.

Stub blocks make continuous integration possible at all

A VEP cache is 25 GB. No pull request is going to download one. Every process therefore carries a stub: block, and CI runs the workflow with -stub-run against a three-record fixture, checking channel wiring and process contracts in seconds without a container or a cache in sight. It catches the failure that actually happens in practice — a renamed output, a channel that emits the wrong cardinality — while leaving scientific correctness to the tests that can afford to be slow.

Event-driven execution, three ways

POST /cases/{id}/annotate writes a job row and hands off. What happens next is a configuration decision, not a code path the caller knows about.

One API call, three execution backends — One entry point, three backends, chosen from settings alone.

Cloud Run job. The Nextflow driver runs as a Cloud Run job started through the Jobs API with argument overrides. It scales to zero between runs, and its service account holds run.jobsExecutorWithOverrides on exactly one job — not project-wide.
Pub/Sub. The job is published as an event. An Argo Events sensor subscribes and triggers an Argo Workflow on GKE. This is the decoupled path: retries, ordering and back-pressure become the queue’s problem rather than the API’s, and the API can be restarted mid-pipeline without losing work.
Local subprocess. Nextflow runs directly and its stdout is streamed into the job log, so the interface shows live progress. This is what a developer gets with nothing configured, and it is the same code path the other two wrap.

All three converge on the same job row, so the interface polls one endpoint regardless. The database URL is passed by environment or as a Nextflow secret and never appears on a command line, keeping it out of .command.sh and the workflow logs.

Kubernetes, Argo and GitOps

The Kubernetes track is Kustomize bases with a local overlay (kind, an in-cluster Postgres) and a GCP overlay (Cloud SQL, Workload Identity). Argo Workflows runs the annotation WorkflowTemplate; Argo Events holds the Pub/Sub EventSource and the sensor that triggers it; ArgoCD reconciles the cluster from the repository, and a green CI run on main bumps image tags in the GCP overlay so that deployment is a commit rather than a command.

One bug from this area is worth recording because it is invisible until it bites. Kustomize generates hashed ConfigMap names so a configuration change forces a rollout, but the hashed name is only substituted into workloads Kustomize believes are in scope. The overlays were missing namespace: rarelens, so the substitution silently did not happen and pods mounted a ConfigMap name that no longer existed. The symptom was a pod stuck in CreateContainerConfigError with nothing wrong in the manifests as written.

Infrastructure as code, with cost as a constraint

Terraform provisions both tracks from one root module, with the expensive half behind flags:

terraform apply -var project=<id>                        # serverless: Cloud Run + Batch
terraform apply -var project=<id> -var deploy_kubernetes=true \
                -var deploy_cloud_sql=true               # adds GKE, Argo, Cloud SQL

A portfolio platform is idle more than 99% of the time, which makes idle cost the only cost that matters. The serverless track is built around that: Cloud Run at min-instances=0, a driver that exists only while a pipeline runs, and Batch on Spot VMs. Idle cost lands near £1 a month, almost all of it the database. CI authenticates through Workload Identity Federation, so there is no service account key anywhere in the repository.

Integrating public biological infrastructure

Almost nothing here is self-contained, and integrating public bioinformatics services is most of the work: Ensembl VEP (cache mode, or against Ensembl’s public database server when 25 GB is not available), gnomAD v4.1 allele frequencies, the Human Phenotype Ontology, ClinVar for training labels filtered to two-star review status, Monarch’s Phenopacket Store, and MLflow as a model registry.

One of those deserves a note. I had written off allele frequencies as impossible without the full VEP cache, having tested only VEP’s own flags. Testing the assumption disproved it: gnomAD’s public bucket files are tabix-indexed, so a range request returns a few hundred records without fetching the file.

# real gnomAD frequencies with no bulk download, verified end to end
bcftools view -r chr3:30672000-30673000 \
  https://storage.googleapis.com/gcp-public-data--gnomad/release/4.1/...chr3.vcf.bgz
# -> 392 records, 5 KB, then fed to VEP with --custom alongside --database

The ranking, and what it refuses to claim

The rank is a weighted mean of four components a reviewer can audit. ClinVar is deliberately not one of them — it sits beside the result as independent confirmation, so nothing ranks highly merely because ClinVar already called it pathogenic.

Evidence that was never looked up must abstain — Rarity and consequence filter; phenotype only ranks, because a real diagnosis can sit in a gene nobody has annotated yet.

The abstention rule replaced a genuine bug. Run without a VEP cache there are no allele frequencies — and the code read a missing frequency as absent from gnomAD, therefore maximally rare, handing every variant a free quarter of its score. The model, separately, was returning 0.887 for every variant on features it had never been given. Two of four components were fiction, and the total looked fully informed.

The fix is structural rather than cosmetic: the loader records what the annotation run actually produced, components without evidence return null instead of a number, and the weights renormalise over whatever is left. The interface prints “not looked up” where it would otherwise have drawn a bar.

Measuring it against 10,178 published cases

One demonstration case ranking correctly is an anecdote. The benchmark asks the only question that matters for a phenotype-driven tool: given a real patient’s reported terms, where does the gene their authors diagnosed rank among all 5,269 HPO-annotated genes?

Benchmarking the phenotype ranking — Retrieval against 5,269 candidate genes, and an honest report of what the measurement cannot tell you.

Two things there matter more than the headline. The first is contamination: HPO’s gene annotations are curated from these same case reports, so the median causal gene already carries every one of its patient’s terms. The number is an upper bound and is labelled as one.

The second is that when I measured my own improvements, one of them did not work. Information-content weighting helped; ontology propagation cost about as much as weighting gained. That result is published with the table rather than quietly dropped, because a benchmark you only consult when it agrees with you is not a benchmark.

The number that changed what the model is allowed to do

The pathogenicity model is LightGBM trained on ClinVar labels, held out by gene rather than by variant. That distinction is not pedantry: a random split puts variants of the same gene on both sides, and the model then scores the gene instead of the variant — exactly the inflation Grimm and colleagues documented for this class of tool in 2015.

Missense AUROC falls to 0.500 without allele frequency — Removing one feature moved missense AUROC from 0.872 to 0.500.

The model originally took allele frequency as a feature and looked respectable. Probing it showed frequency dominating everything — the same missense variant scored 0.887 at frequency zero and 0.0003 at one per cent. That is two separate problems. It double-counted, because the ranking already scores frequency explicitly, putting roughly 45% of every rank on one measurement. And it was circular, because ACMG’s BA1/BS1 criteria assign ClinVar’s benign labels using allele frequency, so the model was rediscovering the rule that had generated its own labels.

Retraining without it returned missense AUROC 0.500 — exactly random. With frequency gone and no CADD or AlphaMissense scores in the training table, nothing is left but the consequence class, so every missense variant scores identically. The conclusion is unambiguous and slightly uncomfortable: the model never had variant-effect knowledge. It now abstains from the ranking unless it has a predictor the other components do not already provide, and 0.500 is the measurement that justifies the abstention.

The application

The demonstration case is a real published patient: the TGFBR2 proband from Loeys and colleagues, Nature Genetics 2005 — the paper that first defined Loeys–Dietz syndrome. Their thirty reported phenotype terms and the causal variant come straight from the publication, read from a GA4GH phenopacket; the background variants come from a public reference genome, because the rest of that patient’s genome is not public.

A published case in the rarelens interface — The funnel across the top is the story. Each candidate carries its evidence as chips, and the note under the header states plainly which components did not score and why.

The evidence panel for a candidate variant — Selecting a candidate opens the arithmetic. Every component is shown with its weight and its contribution, and the two that had nothing to go on say so rather than displaying a zero.

The case report: shortlisted and dismissed variants with the reviewer’s reasons, the funnel counts, and provenance — so a result can be reproduced or challenged later.

Engineering practice

Tests where the risk is. 99 Python tests plus 51 in the front end, concentrated on the ranking arithmetic, the ontology handling, the loader’s idempotency and the security boundary. Database tests run against a real PostgreSQL — an embedded server locally, a service container in CI — because the schema is part of the behaviour.
Migrations are tested both ways. Alembic upgrade and downgrade are exercised against a live database, which is how a downgrade that left a stale enum type behind was caught.
CI runs the whole estate: lint and types for the API, unit tests for the model code, a real Postgres for the loader, a Nextflow stub run, and terraform validate.
Input validation as a security boundary. A case’s VCF URI must be a gs:// object or an absolute path beneath a configured data root, with a VCF suffix. That is what stops a crafted path becoming a Nextflow option or reading an arbitrary file, and it has its own test file.

What is still wrong

An honest introduction should end with the open problems, not the achievements.

The model has no features in its training data. Its only two remaining inputs, CADD and AlphaMissense, are absent from all 688,362 training rows. It abstains today, so nothing is broken — but installing the VEP plugins would flip it on while it still knows nothing. Serving should refuse unless the model was trained on the features it is handed.
The benchmark is contaminated and no amount of code fixes it. A leave-one-publication-out rebuild is possible — HPO’s annotation file carries the source PMID — and is the honest next step.
No baseline comparison. Scoring the same held-out rows with CADD and AlphaMissense would give the model something to beat. AlphaMissense is a 0.64 GB download and tractable; CADD’s whole-genome file is 87.5 GB and is not.
Rarity contributes nothing without the cache in the default demonstration, though the gnomAD streaming route above now makes that solvable.

rarelens is a self-training project built in the open on public data, licensed AGPL-3.0. It is not a clinical tool and makes no diagnostic claim: ACMG/AMP treats computational predictions as supporting evidence only, never sufficient alone.