diff --git a/Architecture.md b/Architecture.md new file mode 100644 index 0000000..680db87 --- /dev/null +++ b/Architecture.md @@ -0,0 +1,57 @@ +# Architecture + +One monorepo holds the scientific pipeline, the API, the interface, model serving and the +infrastructure for two deployment tracks. + +``` + scientist ──▶ SvelteKit UI ──▶ FastAPI ──▶ PostgreSQL + (/api proxy) │ ▲ + │ │ + ├──▶ MLflow registry (model) + │ + └──▶ Nextflow pipeline ──▶ bcftools ──▶ VEP ──▶ loader +``` + +The interface proxies `/api` from its own origin, so there is no CORS anywhere and one public +address serves both. + +## Why the pieces are where they are + +**The pipeline is separate from the API on purpose.** Annotation takes minutes to hours and needs +containers, reference data and a lot of memory; the API needs to answer in milliseconds. They +communicate through a `jobs` row and nothing else, which is what lets the same pipeline run as a +local subprocess, a Cloud Run job or an Argo Workflow without the API knowing. +See [Execution backends](Execution-backends). + +**Feature engineering travels with the model.** The registered MLflow artifact carries +`rarelens_ml.features` as model code, so serving cannot drift from training. The API sends raw +columns and gets a probability back; it does not know what the features are. + +**The ontology work happens at load time, not query time.** HPO annotations are propagated up the +ontology and each term's information content is computed once by `scripts/load-hpo.py`, so the +ranking is a set lookup and a weighted sum. See [Ranking](Ranking). + +## Data model + +``` +cases ──┬── case_phenotypes the patient's HPO terms + └── jobs ──── variants ──┬── predictions model score + └── variant_decisions shortlist / dismiss + reason +gene_phenotypes HPO gene→term, propagated (reference data) +hpo_terms term → information content (reference data) +``` + +`jobs` also records **what the annotation run produced** — `has_frequencies`, +`has_effect_scores`. That is not bookkeeping: the ranking uses it to decide which components are +allowed to score at all. See [Ranking](Ranking). + +## Deployment tracks + +| Track | What runs | When to use it | +|---|---|---| +| **Serverless** (default) | Cloud Run for api + web, Cloud Run job for the Nextflow driver, Google Batch for pipeline tasks | almost always; idles near £1/month | +| **Kubernetes** (flagged off) | GKE Autopilot, Argo Workflows, Argo Events, ArgoCD | to demonstrate the GitOps path, then destroy | +| **Railway** | Postgres + api + web, pre-seeded, no pipeline | putting it in front of people quickly | + +Both cloud tracks run identical pipeline code; the executor is a Nextflow profile. +See [Deployments](Deployments) and [Infrastructure](Infrastructure). diff --git a/Benchmarks.md b/Benchmarks.md new file mode 100644 index 0000000..1c4e29b --- /dev/null +++ b/Benchmarks.md @@ -0,0 +1,77 @@ +# Benchmarks + +Two separate measurements: how well the **phenotype ranking** retrieves the right gene, and how +well the **model** classifies pathogenicity. They are not the same thing and should never be +quoted as one. + +## Phenotype ranking + +```bash +make benchmark # downloads phenopacket-store, exports gene_phenotypes, runs the ranking +``` + +Given a real patient's reported HPO terms, where does the gene their authors diagnosed rank among +all 5,269 HPO-annotated genes? Ties are reported as a range: 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 be 0.02%. + +### The benchmark is contaminated + +HPO's gene-to-phenotype annotations are curated **from these same case reports**. The median +causal gene already carries every one of its patient's terms. So 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. + +Do not quote the headline without this paragraph. + +### A change that did not work + +Information-content weighting and ontology propagation both replaced plain term counting. Asked +whether they helped — pessimistic figures, the 6,485 cases with six or more terms: + +| scoring | top-1 | top-10 | MRR | +|---|---|---|---| +| count terms (the original) | 61.8% | 80.1% | 0.682 | +| **+ information content** | **63.6%** | **83.5%** | **0.706** | +| + propagation | 58.2% | 77.7% | 0.653 | +| + both (shipped) | 59.5% | 81.0% | 0.670 | + +Weighting earns its place. Propagation costs about what weighting gains. It was kept anyway, +because this corpus **cannot show what propagation is for** — its term IDs were chosen by the same +curators HPO records, so exact matching is flattered, while real users pick their own terms from a +search box. That is an argument, not a measurement, and it is labelled as one. + +## The model + +```bash +make training-set && make train +``` + +Held out by gene, never by variant: a random split puts variants of the same gene on both sides and +the model scores the gene instead of the variant (Grimm et al., *Hum Mutat* 2015). +312,025 training and 74,239 test variants; 7,728 and 1,932 genes; no gene on both sides. + +| | 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 is exactly random.** Strip frequency out and nothing is left but the consequence class, so +every missense variant scores identically. The respectable-looking 0.872 was allele frequency, not +variant-effect knowledge — and ACMG's BA1/BS1 criteria assign ClinVar's *benign* labels using +frequency, so the feature had partly caused the label. + +The overall 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. + +This is why the model abstains unless it has CADD or AlphaMissense. See [Roadmap](Roadmap) for why +that is not yet a complete fix. diff --git a/Data-sources.md b/Data-sources.md new file mode 100644 index 0000000..65b8541 --- /dev/null +++ b/Data-sources.md @@ -0,0 +1,55 @@ +# Data sources + +Everything is public, peer-reviewed and consented for open redistribution. **No patient data, no +data access agreement, nothing that needs an application.** Full citations with DOIs are in +`docs/data.md` in the repository — quote from there, not from here. + +| Source | Used for | Licence | +|---|---|---| +| ClinVar | training labels, the ClinVar column | NCBI public domain | +| Genome in a Bottle HG002 | the background genome | open, no restriction | +| gnomAD v4 | allele frequency | free use | +| Human Phenotype Ontology | what the phenotype half matches against | free with attribution | +| Phenopacket Store (Monarch) | the published demonstration case, and the benchmark | BSD-3-Clause | +| Ensembl VEP 113 | annotation | Apache 2.0 | +| CADD, AlphaMissense | optional model features | see `docs/data.md` — CADD needs a commercial licence | + +## Loading it + +```bash +make data # GIAB + ClinVar, chr22 only, ~100 MB +make hpo # HPO annotations, propagated and IC-weighted +make demo-case # a simulated proband: GIAB background + one ClinVar pathogenic variant +make published-case # a real published patient, from a GA4GH phenopacket +make training-set # ~370k labelled variants from ClinVar, 2-star and above +``` + +## The demonstration cases + +**The simulated proband** spikes one real ClinVar pathogenic variant into a real public genome. +It has a right answer, so the ranking can be checked rather than admired. + +**The published case** invents nothing. It reads a GA4GH phenopacket curated from a peer-reviewed +case report and takes the patient's reported HPO terms and the authors' causal variant verbatim. +The default is the *TGFBR2* proband from Loeys et al., *Nat Genet* 2005 — the paper that first +defined Loeys–Dietz syndrome — with 30 reported terms and `NM_003242.6:c.1069G>T`. + +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 around the locus, drawn +from **coding exons** where possible. That last detail matters more than it sounds: of roughly +4,000 HG002 variants in a 3 Mb window only 9 are coding, so a random sample is entirely intronic, +the consequence filter discards all of it, and the causal variant is left as the only candidate — +a funnel that proves nothing. + +## Allele frequencies without the 25 GB cache + +VEP's `--af_gnomade` is rejected outright with `--database`, and plain `--af` returns nothing even +for common variants. But gnomAD's public bucket is tabix-indexed, so a range request works: + +```bash +bcftools view -r chr3:30672000-30673000 \ + https://storage.googleapis.com/gcp-public-data--gnomad/release/4.1/vcf/exomes/...chr3.vcf.bgz +# 392 records, 5 KB -> feed to VEP with --custom alongside --database +``` + +Verified end to end. Not yet wired into the pipeline — see [Roadmap](Roadmap). diff --git a/Deployments.md b/Deployments.md new file mode 100644 index 0000000..227b7a6 --- /dev/null +++ b/Deployments.md @@ -0,0 +1,64 @@ +# Deployments + +Three ways to run it. Pick on cost and on whether you need the pipeline. + +| | Pipeline runs? | Idle cost | Effort | +|---|---|---|---| +| [GCP serverless](#gcp-serverless-default) | yes | ~£1/month | terraform apply | +| [GCP Kubernetes](Kubernetes-and-GitOps) | yes | control plane + Cloud SQL | flag, then destroy | +| [Railway](#railway) | **no** | $5/month flat | `railway up` | + +## GCP serverless (default) + +```bash +cd infra/terraform +terraform init -backend-config="bucket=" +export TF_VAR_database_url='postgresql+asyncpg://user:pass@host/db?sslmode=require' +terraform apply -var project= +cd ../.. && make serverless-deploy PROJECT= TAG= +``` + +`terraform output web_url` is the address to share: it serves the UI and proxies `/api`, so there +is one public address and no CORS. Upload the VEP cache to `gs://-rarelens-data/refs/vep` +before a real annotation, and set `-var model_uri=gs://.../models/pathogenicity/1` to score without +an MLflow server. + +**Set a billing budget first.** The demo has no authentication. + +## Railway + +Three services — Postgres, api, web — with the API on the private network only, so the UI's `/api` +proxy is the single public entry point. + +```bash +railway link --project --environment production --service api +cd api && railway up --service api # the repo is on Gitea, so deploy from the working copy +cd ../web && railway up --service web +``` + +Variables: on the API, `DATABASE_URL` pointing at `postgres.railway.internal` and +`DATABASE_IDLE_CONNECTIONS=false`; on the UI, `API_INTERNAL_URL=http://api.railway.internal:8000`, +`PUBLIC_PIPELINE_ENABLED=false`, `BASIC_AUTH_USER` and `BASIC_AUTH_PASSWORD`. + +Three things to know: + +- **The pipeline cannot run there.** Nextflow shells out to `docker run` for VEP and bcftools, and + a container platform gives you a container, not a Docker daemon. Cases are annotated elsewhere + and copied up with `scripts/seed-remote.sh`; `PUBLIC_PIPELINE_ENABLED=false` hides the buttons + that would otherwise be left to fail. +- **`DATABASE_IDLE_CONNECTIONS=false` is what makes sleeping work.** Railway decides a service is + idle from its *outbound* traffic, and a pooled database connection is outbound traffic, so the + default pool keeps the service awake and billable for ever. It uses `NullPool` instead — a + connection per request, which is nothing at demo traffic and the wrong trade under real load. +- **Serverless must be enabled per service and only applies on the next deploy.** Leave it off for + Postgres, which holds the volume. Expect the documented 502 on the first request to a slept + service; it answers on retry. + +## Seeding a remote database + +```bash +scripts/seed-remote.sh "postgresql://user:pass@host:port/db" +``` + +Copies the analysed cases and the HPO reference data. `--data-only`: the target's schema comes from +Alembic, so the two cannot disagree about it. Run `alembic upgrade head` against the target first. diff --git a/Execution-backends.md b/Execution-backends.md new file mode 100644 index 0000000..c75feda --- /dev/null +++ b/Execution-backends.md @@ -0,0 +1,34 @@ +# Execution backends + +`POST /cases/{id}/annotate` writes a `jobs` row and hands off. `app/services/events.launch()` +then picks a backend from configuration alone — the pipeline code is identical in all three cases. + +| Setting present | Backend | Behaviour | +|---|---|---| +| `cloudrun_job` | Cloud Run job | executes the Nextflow driver with argument overrides via the Jobs API | +| `pubsub_topic` | Pub/Sub | publishes the job; Argo Events triggers an Argo Workflow | +| neither | local subprocess | runs Nextflow directly, streaming stdout into the job log | + +They are checked in that order. + +## Why three + +**Cloud Run job** is the default because it scales to zero: nothing runs between pipelines. The +service account holds `run.jobsExecutorWithOverrides` on exactly one job, not project-wide. + +**Pub/Sub** is the decoupled path. Retries, ordering and back-pressure become the queue's problem +instead of the API's, and the API can restart mid-pipeline without losing work. This is the right +shape once more than one thing produces annotation requests. + +**The local subprocess** is what a developer gets with nothing configured. `_watch` streams the +last `PROGRESS_LINES` of Nextflow's output into the job row every `PROGRESS_INTERVAL_S` seconds so +the UI can show live progress, and `_record_progress` only writes while the job is still `running`, +so a late write cannot resurrect a finished job. + +All three converge on the same `jobs` row. The interface polls one endpoint and does not know or +care which ran. + +## When the API has no Nextflow + +If `shutil.which("nextflow")` finds nothing, the job is marked failed with the command to run by +hand. That is the normal state in docker-compose, and on Railway — see [Deployments](Deployments). diff --git a/Gotchas.md b/Gotchas.md new file mode 100644 index 0000000..f676140 --- /dev/null +++ b/Gotchas.md @@ -0,0 +1,70 @@ +# Gotchas + +Things that cost real debugging time here. Each one is a bug that was found and fixed; they are +recorded because the symptom rarely points at the cause. + +## Science + +**A missing measurement is not a measurement of zero.** `rarity_score(None)` read "this run has no +allele frequencies" as "absent from gnomAD, therefore maximally rare" and handed every variant a +free 0.25. The fix is structural: jobs record what they looked up, and a component with no +evidence abstains. If you add a component, decide what makes it abstain first. See +[Ranking](Ranking). + +**A DAG walked pre-order and read backwards loses ancestors.** Computing each HPO term's ancestors +with a pre-order DFS reversed *looks* like a post-order, but on a DAG a term can be reached before +one of its parents, and then inherits that parent alone instead of the parent's whole lineage. It +silently dropped **399 terms** out of the phenotype branch — including Camptodactyly and Chiari +malformation — so cases mentioning them scored lower for no visible reason. Now a true post-order, +tested against a reference transitive closure. + +**A feature can cause its own label.** The model took allele frequency as a feature, while ACMG's +BA1/BS1 criteria assign ClinVar's *benign* labels using allele frequency. It was rediscovering the +rule that generated its training data, and the impressive missense AUROC of 0.872 fell to exactly +0.500 when the feature was removed. See [Benchmarks](Benchmarks). + +**Split by gene, not by variant.** A random split puts variants of the same gene on both sides and +the model learns the gene. Grimm et al. documented this inflation for exactly this class of tool. + +## Pipeline + +**VEP's coordinates are not a primary key.** `Location` and `Allele` trim indel alleles and shift +positions; `CT>C` comes back as `->` elsewhere. Identity travels in the VCF ID instead. See +[Pipeline](Pipeline). + +**`--af_gnomade` is rejected with `--database`,** and plain `--af` silently returns nothing even +for a variant at 15% global frequency. Frequencies genuinely need the cache — or a range request +against gnomAD's public bucket, see [Data sources](Data-sources). + +**Nextflow rejects two staged inputs sharing a filename.** The cache placeholder and the plugin +placeholder must be different files (`assets/NO_CACHE` and `assets/NO_FILE`), or the run fails with +an input collision that names neither. + +**`head` plus `pipefail` kills a script silently.** `head` closes the pipe, the upstream process +takes SIGPIPE, and with `pipefail` the whole script aborts with no message. + +## Infrastructure + +**Kustomize only substitutes hashed ConfigMap names in workloads it considers in scope.** Missing +`namespace:` in an overlay means the hash is generated but never substituted, and pods mount a name +that does not exist. The symptom is `CreateContainerConfigError` with correct-looking manifests. + +**A connection pool stops a container ever sleeping.** Platforms that sleep idle containers measure +*outbound* traffic; a pooled database connection is outbound traffic. See [Deployments](Deployments). + +## Local environment + +**macOS AirPlay Receiver owns port 5000** and answers 403, which looks exactly like an auth failure +from MLflow. MLflow is on 5001 here. + +**`*.rlwy.net` may resolve to 0.0.0.0** on a machine with DNS filtering, so Railway's TCP proxy +looks refused when it is fine. Check with a public resolver and connect by IP. + +## Security boundaries, not conveniences + +**`vcf_uri` validation is a boundary.** A case's URI must be a `gs://` object or an absolute path +beneath `LOCAL_DATA_ROOT`, ending in a VCF suffix. Without it, a crafted value becomes a Nextflow +option or reads an arbitrary file. It has its own tests; do not relax it for convenience. + +**The database URL never goes on a command line.** Environment or Nextflow secret only, so it stays +out of `.command.sh` and the workflow logs. diff --git a/Home.md b/Home.md index 27b0ee6..bb8ed9f 100644 --- a/Home.md +++ b/Home.md @@ -1,3 +1,48 @@ # rarelens -Initialising the wiki. +Rare-disease variant triage on public data: a proband's variants narrowed against their phenotype, +each candidate carrying the evidence for its rank. + +**This is a research and self-training platform. It is not a clinical tool and makes no diagnostic +claim. No patient data is used, accepted, or possible to load.** + +## Start here + +| If you want to… | Read | +|---|---| +| understand how the pieces fit together | [Architecture](Architecture) | +| run it on your machine | [Local development](Local-development) | +| know what the ranking actually does | [Ranking](Ranking) | +| know how well it works, and how well it doesn't | [Benchmarks](Benchmarks) | +| change the pipeline | [Pipeline](Pipeline) | +| deploy it | [Deployments](Deployments) | +| avoid a trap someone already hit | [Gotchas](Gotchas) | + +## What it does + +1. A **case** is a proband: a VCF plus the patient's phenotype as HPO terms. +2. A **Nextflow** pipeline normalises the VCF, annotates it with Ensembl VEP, and loads the + variants into PostgreSQL. +3. The **API** ranks each variant on four auditable components — phenotype fit, rarity, + consequence severity, and a model score — and refuses to score any component the annotation run + did not actually measure. +4. The **interface** shows the narrowing as a funnel and the evidence behind each candidate. A + reviewer shortlists or dismisses with a reason, and that trail becomes a case report. + +## The shape of the repository + +``` +api/ FastAPI, Pydantic v2, SQLAlchemy 2.0 async, Alembic (~1,380 lines) +ml/ LightGBM model, HPO ontology handling, benchmark (~500 lines) +pipeline/ Nextflow DSL2, bcftools + VEP + loader (~620 lines) +web/ SvelteKit 5, adapter-node (~1,460 lines) +infra/ Terraform (GCP), Kustomize, Argo Workflows, ArgoCD (~890 lines) +scripts/ data preparation and seeding (~505 lines) +docs/ architecture, cloud choice, data provenance, write-up +``` + +## A note on how this wiki is written + +Pages record what is true now, including what is broken. Where a number appears it came from a run +in this repository and can be reproduced by the command next to it. If you change something that +makes a number here wrong, change the number too — a wiki nobody trusts is worse than no wiki. diff --git a/Infrastructure.md b/Infrastructure.md new file mode 100644 index 0000000..8297838 --- /dev/null +++ b/Infrastructure.md @@ -0,0 +1,54 @@ +# Infrastructure + +`infra/terraform` provisions both deployment tracks from one root module, with the expensive half +behind flags. + +```bash +terraform apply -var project= # serverless: Cloud Run + Batch +terraform apply -var project= -var deploy_kubernetes=true \ + -var deploy_cloud_sql=true # adds GKE, Argo, Cloud SQL +``` + +## What the default track creates + +| Resource | Why | +|---|---| +| Cloud Run services (api, web) | `min-instances=0`, so nothing runs when nobody is looking | +| Cloud Run **job** (nextflow driver) | exists only while a pipeline runs | +| Google Batch | VEP tasks on Spot VMs | +| GCS bucket | VCFs, the VEP cache, model artifacts | +| Secret Manager | the database URL, and any Nextflow secrets | +| Artifact Registry | images CI pushes | +| Service accounts + IAM | one per component, least privilege | + +GKE Autopilot and Cloud SQL are opt-in because a Kubernetes control plane and a managed database +are most of what a demonstration estate costs. + +## Cost as a design constraint + +A portfolio platform is idle more than 99% of the time, so idle cost is the only cost that matters. +The reasoning, including why Google Cloud rather than AWS, is in `docs/cloud.md`. The short version: +GKE's free tier covers one control plane where EKS charges about $73/month, and Cloud Run's +scale-to-zero plus Batch on Spot puts the idle bill near £1/month — almost all of it the database. + +The Kubernetes track is meant to be turned on, demonstrated, and destroyed: + +```bash +make serverless-destroy PROJECT= +``` + +## Authentication from CI + +Workload Identity Federation. **No service account key exists in the repository or in repository +secrets.** CI exchanges its OIDC token for short-lived credentials. + +## Configuration that matters + +| Variable | Effect | +|---|---| +| `deploy_kubernetes` | GKE, Argo, the whole second track | +| `deploy_cloud_sql` | managed Postgres instead of an external URL | +| `database_url` | an external Postgres (a free tier, say) for the serverless track | +| `model_uri` | score from an artifact, with no MLflow server running | + +Terraform is `fmt -check`ed and validated on every pull request — see [Testing and CI](Testing-and-CI). diff --git a/Kubernetes-and-GitOps.md b/Kubernetes-and-GitOps.md new file mode 100644 index 0000000..c5e2c0b --- /dev/null +++ b/Kubernetes-and-GitOps.md @@ -0,0 +1,46 @@ +# Kubernetes and GitOps + +Off by default. `terraform apply -var deploy_kubernetes=true` turns it on; it exists to +demonstrate the GitOps path and is meant to be destroyed afterwards. + +## Layout + +``` +infra/k8s/base/ api, web, postgres, ingress +infra/k8s/overlays/local/ kind: in-cluster Postgres, NodePort +infra/k8s/overlays/gcp/ Cloud SQL, Workload Identity, image tags +infra/argo-workflows/ annotate WorkflowTemplate, EventSource, Sensor, RBAC +infra/argocd/app.yaml the Application that reconciles the cluster +``` + +## Locally + +```bash +make kind # builds the images, loads them into a kind cluster, applies the local overlay +``` + +## The GitOps loop + +1. A pull request merges to `main`. +2. CI builds and pushes images to Artifact Registry. +3. CI bumps the image tags in `infra/k8s/overlays/gcp`. +4. ArgoCD sees the commit and reconciles the cluster. + +Deployment is therefore a commit, and rollback is a revert. + +## Argo Workflows and Argo Events + +`annotate.yaml` is a `WorkflowTemplate` running the same containers the Nextflow pipeline uses. +`events.yaml` holds the Pub/Sub `EventSource` and the `Sensor` that triggers the template when the +API publishes a job. `rbac.yaml` grants the sensor permission to create workflows in the +`rarelens` namespace and nothing else. + +## The ConfigMap trap + +Kustomize appends a content hash to generated ConfigMap names so that a config change forces a +rollout. That substitution only happens for workloads Kustomize considers in scope. Both overlays +were missing `namespace: rarelens`, so the hashed name was generated but never substituted, and +pods mounted a ConfigMap name that no longer existed. + +The symptom is a pod stuck in `CreateContainerConfigError` while every manifest reads correctly. +See [Gotchas](Gotchas). diff --git a/Local-development.md b/Local-development.md new file mode 100644 index 0000000..5868f03 --- /dev/null +++ b/Local-development.md @@ -0,0 +1,51 @@ +# Local development + +```bash +make up # postgres + api + web + mlflow via docker-compose +make migrate # alembic upgrade head +make hpo # HPO annotations: what the ranking matches against +make published-case # a real published patient +make test # api, ml, loader and web tests +``` + +Then open . + +## Making the "Analyse case" button work + +docker-compose has no Nextflow, so the button marks the job failed with the command to run instead. +To get the real thing, run the API **on the host**, where Nextflow and Docker are: + +```bash +docker compose up -d db +cd api && DATABASE_URL=postgresql+asyncpg://rarelens:rarelens@localhost:5432/rarelens \ + PIPELINE_DATABASE_URL=postgresql+asyncpg://rarelens:rarelens@host.docker.internal:5432/rarelens \ + LOCAL_DATA_ROOT=$PWD/.. VEP_DATABASE=true \ + uv run --extra dev uvicorn app.main:app --port 8000 +``` + +- `PIPELINE_DATABASE_URL` is what the loader container gets. Inside it, the API's own `localhost` + would be the container itself. +- `LOCAL_DATA_ROOT` is the directory a case's `vcf_uri` must sit under. This is a security + boundary, not a convenience — see [Gotchas](Gotchas). + +## Ports + +| Service | Port | Note | +|---|---|---| +| web | 5173 | | +| api | 8000 | | +| postgres | 5432 | | +| MLflow | **5001** | not 5000: macOS AirPlay Receiver owns 5000 and answers 403 | + +## Training a model + +```bash +make training-set # ClinVar-derived table, ~370k labelled variants +make train # fits, reports held-out metrics by gene split, moves the production alias +``` + +The API scores with `models:/rarelens-pathogenicity@production`. Setting `MODEL_URI` to an artifact +path skips the registry entirely, which is how the serverless deployment scores with no MLflow +server running. + +What the metrics do and do not mean is in [Benchmarks](Benchmarks). Read that before quoting them. diff --git a/Pipeline.md b/Pipeline.md new file mode 100644 index 0000000..0bdab8a --- /dev/null +++ b/Pipeline.md @@ -0,0 +1,69 @@ +# Pipeline + +`pipeline/main.nf` is a three-process Nextflow DSL2 workflow. Each process runs in a pinned +container. + +| Process | Container | Does | +|---|---|---| +| `NORMALISE` | bcftools 1.20 | left-align, split multi-allelics, set the VCF ID to `CHROM_POS_REF_ALT` | +| `VEP` | ensembl-vep 113.0 | consequence, gene, HGVS, ClinVar; optionally CADD and AlphaMissense | +| `LOAD_DB` | python + psycopg | parse VEP `--tab`, insert variants, mark the job succeeded | + +## Running it + +```bash +# dry run: annotate without touching the database +make pipeline VCF=data/example.vcf.gz + +# annotate for a real job created in the UI +make annotate JOB= VCF=data/example.vcf.gz + +# no 25 GB VEP cache? query Ensembl's public database instead +VEP_DATABASE=true make pipeline VCF=pipeline/tests/data/tiny.vcf +``` + +Database mode takes roughly 25–35 seconds per variant and returns no allele frequencies and no +plugin scores, so keep those runs to tens of variants. What that costs the ranking is described in +[Ranking](Ranking). + +## Variant identity + +**VEP's `Location` and `Allele` columns cannot be used as a key.** They trim indel alleles and +shift positions, so `CT>C` comes back as `->` at a different coordinate. `NORMALISE` therefore +writes identity into the VCF ID, which VEP passes through untouched: + +```bash +bcftools annotate --set-id '%CHROM\_%POS\_%REF\_%FIRST_ALT' +``` + +and `load_db.parse_variant_id` reads it back. It `rsplit`s on `_` from the right, because contig +names can contain underscores (`chrUn_KI270742v1`) while positions and alleles cannot. + +Verified on a real run: `22:42126611 CT>C` survives with its alleles intact. + +## Executors + +The workflow never names an executor. `pipeline/nextflow.config` maps profiles: + +- **local / docker** — a laptop, or the API's subprocess backend +- **gcp** — Google Batch, Spot VMs, needs `--project` and `--bucket` +- **Argo** — the same containers as workflow steps, see [Kubernetes and GitOps](Kubernetes-and-GitOps) + +## Stub runs + +Every process has a `stub:` block, so CI can run the whole workflow with no containers, no cache +and no database: + +```bash +cd pipeline && nextflow run main.nf -stub-run --vcf tests/data/tiny.vcf +``` + +This runs on every pull request. It checks wiring and channel shapes — a renamed output, a process +emitting the wrong cardinality — which is the class of breakage that otherwise only shows up in a +two-hour annotation run. + +## Secrets + +The database URL is passed by environment or as a Nextflow secret, **never on a command line**, so +it does not land in `.command.sh` or the workflow logs. `PIPELINE_DATABASE_URL` exists because the +loader runs inside a container where the API's own `localhost` would be the container itself. diff --git a/Ranking.md b/Ranking.md new file mode 100644 index 0000000..f297e1a --- /dev/null +++ b/Ranking.md @@ -0,0 +1,76 @@ +# Ranking + +`api/app/services/triage.py`. The rank is a **weighted mean of four components a reviewer can +audit**, not a black box. + +| Component | Weight | Scores when | +|---|---|---| +| phenotype fit | 0.35 | always | +| rarity in gnomAD | 0.25 | the run looked up frequencies | +| consequence severity | 0.20 | always | +| model P(pathogenic) | 0.20 | the run has CADD or AlphaMissense | + +**ClinVar is deliberately not an input.** It is displayed beside the result as independent +confirmation, so nothing ranks highly merely because ClinVar already called it pathogenic. + +Rarity and consequence **filter** — the usual first pass. Phenotype only **ranks**, because a real +diagnosis can sit in a gene nobody has annotated yet, and filtering on phenotype would hide exactly +that case. + +## Evidence that was never looked up abstains + +This is the most important rule in the file, and it replaced a real bug. + +Run without a VEP cache there are no allele frequencies. `rarity_score(None)` read that as *absent +from gnomAD, therefore maximally rare* and gave every variant a free 0.25. The model separately +returned the same 0.887 for every variant, from features it had never been given. Two of four +components were fiction and the total looked fully informed. + +Now `jobs.has_frequencies` and `jobs.has_effect_scores` record what the run actually produced, a +component with no evidence returns `None` rather than a number, and `combine()` renormalises the +weights over whatever is left. The score stays on 0–1 and means the same thing; the UI prints +"not looked up" instead of drawing a bar. + +```python +def combine(components): + weight = sum(WEIGHTS[n] for n, v in components.items() if v is not None) + return sum(WEIGHTS[n] * v for n, v in components.items() if v is not None) / weight +``` + +**If you add a component, decide what makes it abstain before you decide its weight.** + +## Phenotype fit + +Information-content-weighted recall: the share of the *total specificity* of the patient's terms +that this gene accounts for. + +``` +score = Σ IC(matched terms) / Σ IC(all the patient's terms) +``` + +`IC(term) = -ln(fraction of annotated genes carrying it)`, computed at load time. It matters: +"Dilated left subclavian artery" scores 7.88, "Global developmental delay" 0.93 — an 8.5× gap that +plain term counting threw away. + +A term HPO has never annotated to any gene gets `DEFAULT_IC`, keeping it in the denominator so it +depresses every gene equally. That is the neutral choice, not an oversight. + +## The model + +Trained on ClinVar, held out **by gene** rather than by variant. Allele frequency is deliberately +not a feature — it double-counted against `rarity_score` and was circular, since ACMG assigns +ClinVar's benign labels using frequency. Removing it dropped missense AUROC from 0.872 to **0.500, +exactly random**, which is why the model now abstains without CADD or AlphaMissense. +Full working in [Benchmarks](Benchmarks). + +## Worked example + +The published Loeys–Dietz case, run in VEP database mode: + +| | score | phenotype (0.64) | rarity | consequence (0.36) | model | +|---|---|---|---|---|---| +| *TGFBR2* 3:30672252 missense | **0.855** | 1.00 (30/30) | not looked up | 0.60 | not looked up | +| *OSBPL10* 3:31748090 missense | 0.218 | 0.00 | not looked up | 0.60 | not looked up | + +Both are rare missense variants, identical on every piece of evidence this run holds except one. +The phenotype is what separates a published diagnosis from an incidental variant. diff --git a/Roadmap.md b/Roadmap.md new file mode 100644 index 0000000..6cb2433 --- /dev/null +++ b/Roadmap.md @@ -0,0 +1,50 @@ +# Roadmap + +Known gaps, most consequential first. Nothing here is hidden in an issue tracker. + +## 1. 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** — `make training-set` builds from the ClinVar VCF, which carries neither. + +Today the model abstains, so nothing is broken. But installing the VEP plugins flips +`has_effect_scores` true, the model stops abstaining and takes 0.20 of every rank **while still +knowing nothing**. That is the same class of bug as the one the abstention rule fixed, one step +further out. + +*Fix:* log at training time which features actually had values, and refuse to score at serving +unless the model was trained on the features it is handed. + +## 2. AlphaMissense in the training set + +0.64 GB, tabix-indexed, tractable. It would close gap 1 and give the model something real to +contribute — missense is exactly where it currently scores 0.500. CADD's whole-genome file is +87.5 GB and is not worth it; its web API is an option for small volumes. + +## 3. gnomAD frequencies without the cache + +Proven to work (see [Data sources](Data-sources)) but not wired into the pipeline. It would make +the rarity component real in the default demonstration instead of abstaining. + +*Shape:* a pipeline step that streams the regions covering the case's variants into a small local +slice, then `vep --custom` against it alongside `--database`. Note VEP's Perl HTS cannot open the +remote index directly, and gnomAD's contigs need renaming to match. + +## 4. Decontaminating the benchmark + +HPO's `phenotype.hpoa` carries the source PMID for each annotation, so a **leave-one-publication- +out** rebuild is possible: score each case against a gene-to-phenotype table with that case's own +paper excluded. That would turn an upper bound into something much closer to a prospective number. + +## 5. A baseline to beat + +Scoring the same held-out rows with CADD and AlphaMissense would put the model's numbers next to +something meaningful rather than leaving them unanchored. + +## 6. Smaller things + +- ClinVar review status is enforced at training (2-star and above) but is not shown in the UI, + because VEP's `CLIN_SIG` does not carry `CLNREVSTAT`. A `--custom` annotation would fix it. +- The deployed demo has one shared credential and no per-user accounts. +- The funnel numbers come from tens of variants. A real exome has thousands of rare coding + variants; nothing here has been run at that scale. diff --git a/Testing-and-CI.md b/Testing-and-CI.md new file mode 100644 index 0000000..a824a81 --- /dev/null +++ b/Testing-and-CI.md @@ -0,0 +1,45 @@ +# Testing and CI + +```bash +make test # api (99 Python tests), ml, pipeline loader, web (49) +make lint # ruff + mypy for the API, svelte-check for the front end +``` + +## What is tested, and why there + +Tests are concentrated where a mistake would be silent rather than loud. + +- **`api/tests/test_triage.py`** — the ranking arithmetic, as pure logic. This is the scientific + claim the app makes, so it is tested without a database in the way. Includes the regression for + components scoring evidence nobody looked up. +- **`ml/tests/test_hpo.py`** — ontology propagation and information content, including a property + test against a reference transitive closure. That test exists because a subtly wrong graph walk + silently deleted 399 terms; see [Gotchas](Gotchas). +- **`pipeline/tests/test_load_db.py`** — the loader, against a real PostgreSQL. Idempotency is + asserted directly: loading twice must not duplicate variants. +- **`api/tests/test_security.py`** — the `vcf_uri` boundary. +- **`api/tests/test_migrations.py`** — Alembic upgrade *and* downgrade against a live database. + +Database tests use an embedded PostgreSQL (`pgserver`) locally and a service container in CI, +because the schema is part of the behaviour and SQLite would not exercise it. + +## CI + +Runs on pull requests and pushes to `main`: + +| Job | Does | +|---|---| +| api | ruff, mypy, pytest against a Postgres service | +| ml | pytest | +| web | `npm run check`, vitest, production build | +| pipeline | loader tests against Postgres, then `nextflow run main.nf -stub-run` | +| terraform | `fmt -check` and `validate` | +| images | builds and pushes to Artifact Registry, `main` only | + +CI authenticates to Google Cloud with Workload Identity Federation. **There is no service account +key in the repository or in repository secrets.** + +## A caution + +CI does not run on branch pushes, only on pull requests and `main`. Work on a long-lived branch is +therefore unverified by CI until the PR is opened — run `make test && make lint` locally. diff --git a/_Sidebar.md b/_Sidebar.md new file mode 100644 index 0000000..8d7d229 --- /dev/null +++ b/_Sidebar.md @@ -0,0 +1,23 @@ +### rarelens + +**Understanding it** +- [Home](Home) +- [Architecture](Architecture) +- [Ranking](Ranking) +- [Benchmarks](Benchmarks) +- [Data sources](Data-sources) + +**Working on it** +- [Local development](Local-development) +- [Pipeline](Pipeline) +- [Execution backends](Execution-backends) +- [Testing and CI](Testing-and-CI) + +**Running it** +- [Deployments](Deployments) +- [Kubernetes and GitOps](Kubernetes-and-GitOps) +- [Infrastructure](Infrastructure) + +**When it goes wrong** +- [Gotchas](Gotchas) +- [Roadmap](Roadmap)