rarelens/scripts/make-published-case.py

#!/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=<ID={c}>\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()