diff --git a/.gitignore b/.gitignore index 7fcc305..f87a07a 100644 --- a/.gitignore +++ b/.gitignore @@ -17,4 +17,5 @@ mlruns/ .terraform/ data/*.vcf* data/*.tsv +data/*.case.json !data/README.md diff --git a/Makefile b/Makefile index b85b5d0..3e149dc 100644 --- a/Makefile +++ b/Makefile @@ -1,4 +1,4 @@ -.PHONY: up down clean migrate test lint data hpo demo-case training-set train loader pipeline annotate images kind serverless-deploy serverless-destroy gcp-configure gcp-secrets +.PHONY: up down clean migrate test lint data hpo demo-case published-case training-set train loader pipeline annotate images kind serverless-deploy serverless-destroy gcp-configure gcp-secrets VCF ?= data/example.vcf.gz MLFLOW_URI ?= http://localhost:5001 @@ -37,6 +37,9 @@ hpo: ## load HPO gene-to-phenotype annotations, which the ranking matches again demo-case: ## build the simulated proband: GIAB background + one ClinVar pathogenic variant scripts/make-demo-case.sh +published-case: ## build a case from a published patient: a GA4GH phenopacket + GIAB background + scripts/make-published-case.py + training-set: ## build a ClinVar training table, shaped like VEP --tab output scripts/make-training-set.sh diff --git a/README.md b/README.md index d4851fd..bd5d2e3 100644 --- a/README.md +++ b/README.md @@ -33,6 +33,7 @@ make up # postgres + api + web + mlflow via docker-compose make migrate # alembic upgrade head make hpo # HPO gene-to-phenotype annotations: what the ranking matches against make demo-case # a simulated proband: GIAB background + one ClinVar pathogenic variant +make published-case # a real published patient: their reported phenotype and causal variant make test # api, ml, loader and web tests (no Docker needed for the DB tests) ``` @@ -41,6 +42,15 @@ give it the phenotype of the planted disease (for the default NF2 case: bilatera schwannoma, sensorineural hearing impairment, tinnitus, meningioma, cataract), and analyse it. The planted variant should come back ranked first. +`make published-case` is the same idea with nothing invented. It builds a case from a GA4GH +phenopacket curated from a peer-reviewed case report — by default the *TGFBR2* proband from +Loeys et al., *Nat Genet* 2005, [10.1038/ng1511](https://doi.org/10.1038/ng1511), the paper that +first described Loeys-Dietz syndrome. The patient's 30 reported HPO terms and their causal +variant come straight from the publication; the background variants come from GIAB HG002, because +the rest of that patient's genome is not public. It writes the phenotype list alongside the VCF, +so the case can be created exactly as reported. See [docs/data.md](docs/data.md) for the +provenance and for what this case does and does not demonstrate. + The docker-compose API has no Nextflow, so "Run VEP annotation" marks the job failed with the command to run instead. With Nextflow and Docker on the host, a VEP cache in `pipeline/cache/vep` and a VCF under `data/` (see [data/README.md](data/README.md)): diff --git a/data/README.md b/data/README.md index e3ed81a..16849ad 100644 --- a/data/README.md +++ b/data/README.md @@ -1,7 +1,9 @@ # Test data -`make data` fetches the demo slice and `make demo-case` builds the simulated proband. Provenance, licences, citations and the -evaluation plan live in [../docs/data.md](../docs/data.md). +`make data` fetches the demo slice, `make demo-case` builds the simulated proband, and +`make published-case` builds a case from a published patient (their reported phenotype and +causal variant, read from a GA4GH phenopacket, inside a GIAB background genome). Provenance, +licences, citations and the evaluation plan live in [../docs/data.md](../docs/data.md). No patient data. Use public sources only: diff --git a/docs/data.md b/docs/data.md index cc23999..e786c44 100644 --- a/docs/data.md +++ b/docs/data.md @@ -29,6 +29,7 @@ benchmarks variant callers against, so it is both realistic and unambiguously sh | **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 @@ -59,6 +60,66 @@ instead of the 25 GB cache. It returns no gnomAD frequencies, so every variant l gnomAD and the rarity term stops discriminating. 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, 21 "rare" (see the caveat below), **2** surviving the +consequence filter, **1** matching the phenotype. + +| | score | phenotype | rarity | consequence | model | +|---|---|---|---|---|---| +| *TGFBR2* 3:30672252 missense | **0.897** | 1.00 (30/30 terms) | 1.00 | 0.60 | 0.887 | +| *OSBPL10* 3:31748090 missense | 0.547 | 0.00 | 1.00 | 0.60 | 0.887 | + +This is the whole argument for phenotype-driven triage in one table. Both are rare missense +variants; the model scores them **identically**, to three decimal places, because nothing about +the variants themselves distinguishes them. 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. +- **Rarity is not doing any work without a VEP cache.** See the caveat above: in database mode + every variant looks absent from gnomAD, so the funnel's rarity step passes everything and every + variant scores a full 1.0 on rarity. `--af_gnomade` is rejected outright with `--database`, and + plain `--af` returns nothing even for common variants — checked against rs429358, roughly 15% + globally. Frequencies need the cache; there is no shortcut. +- **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 +``` + ## The model, and what its numbers mean `make training-set` builds a training table straight from ClinVar rather than running VEP over diff --git a/scripts/make-published-case.py b/scripts/make-published-case.py new file mode 100755 index 0000000..d1d7242 --- /dev/null +++ b/scripts/make-published-case.py @@ -0,0 +1,222 @@ +#!/usr/bin/env python3 +"""Build a demo case from a PUBLISHED patient described in a peer-reviewed case report. + +Nothing about the patient is invented here. Their phenotype terms and their causal variant are +the ones the authors reported, read from a GA4GH phenopacket in Monarch's phenopacket-store +(BSD-3-Clause), which curates published case reports into machine-readable phenopackets and keeps +the PMID on each one. + +What is *not* public is the rest of that patient's genome, and triage is only meaningful when the +causal variant has to be found among others. So a real public genome (GIAB HG002) supplies the +background variants, exactly as phenotype-driven triage tools are benchmarked. The resulting VCF +is therefore a published diagnosis inside a public background genome, not a real person's exome. + +The background is drawn from coding exons where it can be. Variants picked at random from a 3 Mb +window are almost all intronic, the consequence filter discards every one of them, and the causal +variant is then the only candidate left standing -- which demonstrates nothing. Real coding +variants from HG002 give the ranking something it has to rank *against*. + +Citations, licences and caveats: docs/data.md +""" +import argparse +import bisect +import json +import re +import subprocess +import sys +import urllib.request +from pathlib import Path + +STORE = "https://raw.githubusercontent.com/monarch-initiative/phenopacket-store/main/notebooks" +# Loeys et al., Nat Genet 2005 (PMID:15731757, doi:10.1038/ng1511): the paper that first defined +# Loeys-Dietz syndrome. Family 4 II-1 carries TGFBR2 c.1069G>T and has 30 reported HPO terms. +DEFAULT = f"{STORE}/TGFBR2/phenopackets/PMID_15731757_Family_4_II_1.json" + +GIAB = ( + "https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/release/AshkenazimTrio" + "/HG002_NA24385_son/NISTv4.2.1/GRCh38/HG002_GRCh38_1_22_v4.2.1_benchmark.vcf.gz" +) +BCFTOOLS = "quay.io/biocontainers/bcftools:1.20--h8b25389_0" +ENSEMBL = "https://rest.ensembl.org" # GRCh38 by default; the overlap endpoint caps a region at 5 Mb +# VEP's database mode takes roughly 25s per variant, so the case stays deliberately small. +BACKGROUND = 20 +WINDOW = 1_500_000 +SLUG = re.compile(r"[^a-z0-9]+") + + +def fetch(url: str) -> dict: + with urllib.request.urlopen(url) as response: + return json.load(response) + + +def causal_variant(packet: dict) -> dict: + """The one variant the authors called causal, with its GRCh38 coordinates and HGVS.""" + found = [] + for interpretation in packet.get("interpretations", []): + for genomic in interpretation.get("diagnosis", {}).get("genomicInterpretations", []): + descriptor = genomic.get("variantInterpretation", {}).get("variationDescriptor", {}) + if record := descriptor.get("vcfRecord"): + found.append((record, descriptor)) + if len(found) != 1: + sys.exit(f"expected exactly one variant with VCF coordinates, found {len(found)}") + record, descriptor = found[0] + if record.get("genomeAssembly") not in ("hg38", "GRCh38"): + sys.exit(f"{record.get('genomeAssembly')} is not GRCh38; the pipeline expects GRCh38") + return { + "chrom": record["chrom"].removeprefix("chr"), + "pos": int(record["pos"]), + "ref": record["ref"], + "alt": record["alt"], + "gene": descriptor.get("geneContext", {}).get("symbol"), + "hgvs": {e["syntax"]: e["value"] for e in descriptor.get("expressions", [])}, + } + + +def phenotype_terms(packet: dict) -> list[dict]: + """Observed terms only: an excluded term means the authors looked and did not find it.""" + return [ + {"hpo_id": f["type"]["id"], "label": f["type"]["label"]} + for f in packet.get("phenotypicFeatures", []) + if not f.get("excluded") + ] + + +def provenance(packet: dict) -> dict: + reference = (packet.get("metaData", {}).get("externalReferences") or [{}])[0] + disease = next( + ( + i.get("diagnosis", {}).get("disease", {}).get("label") + for i in packet.get("interpretations", []) + if i.get("diagnosis", {}).get("disease") + ), + None, + ) + return { + "phenopacket_id": packet.get("id"), + "pmid": reference.get("id"), + "title": reference.get("description"), + "disease": disease, + } + + +def coding_intervals(chrom: str, start: int, end: int) -> list[tuple[int, int]]: + """Merged coding-exon spans in the window, from Ensembl's public REST API.""" + url = f"{ENSEMBL}/overlap/region/human/{chrom}:{start}-{end}?feature=cds;content-type=application/json" + with urllib.request.urlopen(url) as response: + features = json.load(response) + merged: list[list[int]] = [] + for s, e in sorted((f["start"], f["end"]) for f in features): + if merged and s <= merged[-1][1] + 1: + merged[-1][1] = max(merged[-1][1], e) + else: + merged.append([s, e]) + return [(s, e) for s, e in merged] + + +def is_coding(pos: int, intervals: list[tuple[int, int]]) -> bool: + i = bisect.bisect_right(intervals, (pos, float("inf"))) - 1 + return i >= 0 and intervals[i][0] <= pos <= intervals[i][1] + + +def spread(rows: list, count: int) -> list: + """Take `count` items spaced across the list, not the first `count` in one gene.""" + if len(rows) <= count: + return rows + step = len(rows) / count + return [rows[int(i * step)] for i in range(count)] + + +def giab_variants(region: str) -> list[tuple[str, int, str, str]]: + """Biallelic short variants from GIAB HG002, streamed from the indexed public VCF.""" + print(f"==> background: GIAB HG002 {region}", file=sys.stderr) + out = subprocess.run( + ["docker", "run", "--rm", BCFTOOLS, "bcftools", "view", "-H", "-r", region, GIAB], + capture_output=True, text=True, check=True, + ).stdout + rows = [] + for line in out.splitlines(): + chrom, pos, _id, ref, alt, *_ = line.split("\t") + if "," in alt or len(ref) >= 20 or len(alt) >= 20: + continue + rows.append((chrom.removeprefix("chr"), int(pos), ref, alt)) + if not rows: + sys.exit(f"no GIAB variants in {region}") + return rows + + +def background(chrom: str, start: int, end: int, count: int) -> list[tuple[str, int, str, str]]: + """Prefer HG002's coding variants, then fill out the rest of the window. + + Coding ones are what make the run a triage rather than a formality: they survive the + consequence filter and have to be ranked below the causal variant on phenotype and rarity. + """ + rows = giab_variants(f"chr{chrom}:{start}-{end}") + intervals = coding_intervals(chrom, start, end) + coding = [r for r in rows if is_coding(r[1], intervals)] + print(f" {len(rows)} variants in the window, {len(coding)} in coding exons", file=sys.stderr) + chosen = spread(coding, count) + if len(chosen) < count: + rest = [r for r in rows if r not in set(chosen)] + chosen += spread(rest, count - len(chosen)) + return chosen + + +def write_vcf(path: Path, variants: list[tuple[str, int, str, str]], source: str) -> None: + body = "".join( + f"{c}\t{p}\t.\t{r}\t{a}\t.\tPASS\t.\n" for c, p, r, a in sorted(variants, key=lambda v: v[1]) + ) + contigs = "".join(f"##contig=\n" for c in sorted({v[0] for v in variants})) + header = f"##fileformat=VCFv4.2\n##source={source}\n{contigs}#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n" + raw = path.with_suffix("") + raw.write_text(header + body) + subprocess.run( + ["docker", "run", "--rm", "-v", f"{path.parent.resolve()}:/out", BCFTOOLS, + "bash", "-eu", "-c", f"bgzip -f /out/{raw.name} && tabix -f -p vcf /out/{path.name}"], + check=True, + ) + + +def main() -> None: + p = argparse.ArgumentParser() + p.add_argument("--phenopacket", default=DEFAULT, help="phenopacket-store JSON URL") + p.add_argument("--background", type=int, default=BACKGROUND) + p.add_argument("--window", type=int, default=WINDOW, help="bp either side of the variant") + p.add_argument("--out-dir", default="data") + a = p.parse_args() + + packet = fetch(a.phenopacket) + variant = causal_variant(packet) + terms = phenotype_terms(packet) + source = provenance(packet) + gene = variant["gene"] or "unknown" + + start, end = max(1, variant["pos"] - a.window), variant["pos"] + a.window + region = f"chr{variant['chrom']}:{start}-{end}" + variants = background(variant["chrom"], start, end, a.background) + causal = (variant["chrom"], variant["pos"], variant["ref"], variant["alt"]) + variants = [v for v in variants if v != causal] + [causal] + + slug = SLUG.sub("-", f"{gene} {source['pmid'] or ''}".lower()).strip("-") + out = Path(a.out_dir) + out.mkdir(parents=True, exist_ok=True) + vcf = out / f"proband-{slug}.vcf.gz" + write_vcf(vcf, variants, f"rarelens published case: {source['pmid']} {gene} + GIAB HG002 background") + + case = { + "name": f"{gene} {source['pmid']}", + "vcf_uri": str(vcf.resolve()), + "assembly": "GRCh38", + "phenotypes": terms, + "_provenance": source | {"causal_variant": variant, "background": {"source": "GIAB HG002 v4.2.1", "region": region, "n": len(variants) - 1}}, + } + (out / f"proband-{slug}.case.json").write_text(json.dumps(case, indent=2) + "\n") + + print(f"\nwrote {vcf} ({len(variants)} variants) and {out / f'proband-{slug}.case.json'}") + print(f" published case: {source['pmid']} - {source['disease']}") + print(f" reported diagnosis: {gene} {variant['hgvs'].get('hgvs.c', '')} {variant['hgvs'].get('hgvs.p', '')}") + print(f" at {variant['chrom']}:{variant['pos']} {variant['ref']}>{variant['alt']}, among {len(variants) - 1} GIAB background variants") + print(f" reported phenotype: {len(terms)} HPO terms") + + +if __name__ == "__main__": + main()