acmg-variant-curation

14 files shipped with this skill. Inline text is the catalog copy. The upstream link is the source repository.

• SKILL.md15 KBc2edb5a96193upstream ↗

  ---
  name: acmg-variant-curation
  description: Use when deriving ACMG codes for a variant from ClinVar, gnomAD, VEP and PubMed, or curating a VUS. Fetches all four sources live and prints the basis for every code.
  license: MIT
  compatibility: Requires python3 (stdlib only, no pip install) and outbound HTTPS to eutils.ncbi.nlm.nih.gov, gnomad.broadinstitute.org and rest.ensembl.org. A sandbox that blocks network access disables the whole skill.
  metadata:
    version: "2.1.0"
    author: Swapnil Padmane
    repository: https://github.com/darshak26/acmg-variant-curation
    spec: https://agentskills.io/specification
  ---
  
  # ACMG variant curation
  
  Derive an evidence-backed ACMG lean for a germline SNV/indel using four public sources
  and the ACMG/AMP 2015 rulebook. Built for the TP53 project but gene-agnostic.
  
  Scope, stated honestly: the seven codes below are the ones derivable from public
  databases. Curated criteria — PS4, PM6, PP1, PS2 — need case-level data and are **not**
  covered. This skill does not handle CNVs, structural variants or somatic curation; the
  PP3/BP4 machinery means nothing there. Out-of-scope input is rejected with exit code 2 and a
  pointer to the right standard, never half-curated.
  
  Five steps, in this order: **ClinVar → gnomAD → VEP → PubMed → ACMG**. Order matters —
  population frequency can end curation early (BA1), and functional evidence outranks every
  predictor.
  
  ## Hard rules
  
  1. **Never state a score you did not fetch in this session.** No recalled REVEL values, no
     remembered frequencies. Run the script, quote its output. An identifier is the same kind
     of value: resolve the ClinVar id and *derive* the gnomAD id from ClinVar's own SPDI — a
     coordinate written out by hand is a guess that looks like data.
  2. **NA is not zero.** A missing predictor is excluded from numerator *and* denominator. It
     never votes damaging, never votes tolerated.
  3. **PS3/BS3 require a paper you actually read.** A title match is not a functional assay.
     If the exact protein change returns zero hits, the honest answer is "no functional data".
  4. **Pin the versions.** gnomAD dataset (`gnomad_r4`), transcript accession, predictor score
     versions, and the date fetched. Unpinned frequencies are unreproducible.
  5. **Declare your scheme, in the sentence that reports the code.** PP3 has two defensible
     implementations (below) that disagree on real variants — D393V fires PP3 under `count3`
     and yields **no code at all** under `revel`. Say which scheme in Methods *and* name it
     beside the call wherever the code is quoted ("PP3 under count3, 4 of 4 available"). A
     bare "PP3 (4 of 4)" is refutable by anyone who runs the other scheme, and the person who
     refutes it will not believe the rest of the table. Never mix schemes between variants.
  6. **Check `consistency.ok` before you read any code.** The three identifiers must be
     proven to describe one variant. `fetch_evidence.py` now emits a `consistency` block;
     `ok: false` means the payload must not be curated. `acmg_codes.py` refuses PM2 on it.
  7. **A curation is a proposal, not a classification.** VUS work in a student project is
     evidence assembly; clinical use needs a ClinGen expert panel.
  
  ## Run it
  
  Paths below are relative to this skill's root, per the Agent Skills spec, so the block
  runs unchanged on any host. Hermes installs the skill at
  `~/.hermes/profiles/athena/skills/research/acmg-variant-curation/` — `cd` there first if
  your host does not make the skill root the working directory.
  
  ```bash
  # 0. confirm all four endpoints answer before trusting a batch
  python3 scripts/fetch_evidence.py --selftest
  
  # 1-3. fetch ClinVar + gnomAD + VEP (same variant, three identifiers)
  python3 scripts/fetch_evidence.py \
      --clinvar 3634421 \
      --gnomad 17-7669613-T-A \
      --hgvs "NM_000546.6:c.1178A>T" --transcript ENST00000269305 \
      --pubmed-terms "p.Asp393Val,D393V" \
      --out D393V.json
  
  # 4. read PubMed by hand (see the return payload's action_required)
  
  # 4b. trust nothing until this says ok: true
  python3 -c "import json;c=json.load(open('D393V.json'))['consistency'];print(c['ok'],c['problems'])"
  
  # 5. derive codes
  python3 scripts/acmg_codes.py D393V.json --scheme count3
  ```
  
  Identifiers must describe the same variant: ClinVar variation ID, GRCh38 gnomAD id
  (`chr-pos-ref-alt`), HGVS c. with accession, and the MANE Select ENST.
  
  ### Resolve identifiers, never compose them
  
  The guard catches bad identifiers; it cannot supply good ones. Get each id from a source of
  record, in this order — resolving takes ~15 s and removes the whole failure class.
  
  ```bash
  # 1. variation id from a tight gene + protein-change query (expect 1 hit, not 50)
  curl -s "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=clinvar&retmode=json&term=TP53%5Bgene%5D+AND+G245S"
  
  # 2. canonical GRCh38 SPDI for that id
  curl -s "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=clinvar&retmode=json&id=12365" \
   | python3 -c "import json,sys; r=json.load(sys.stdin)['result']['12365'];
  print(r['title']); print(r['variation_set'][0]['canonical_spdi'])"
  # NM_000546.6(TP53):c.733G>A (p.Gly245Ser)
  # NC_000017.11:7674229:C:T     <- 0-based position, ref, alt
  # gnomAD id = 17-7674230-C-T  <- position + 1, alleles unchanged
  ```
  
  - Read the ClinVar `title` and confirm the protein change is the one you meant *before*
    using its id. `esearch` on a protein change can return a neighbouring codon.
  - Derive the gnomAD id from the SPDI: **0-based position + 1**, alleles as printed. Never
    count base pairs from a remembered codon number.
  - Check the accession build: `NC_000017.11` is GRCh38 (gnomAD v4). A `.10` accession is
    GRCh37 and its coordinates will be shifted.
  - Take the `--hgvs` accession from the same ClinVar title, and the ENST from the transcript
    the record actually annotates — do not pair a fresh id with a cached transcript.
  
  ### What the consistency guard proves, and what it does not
  
  `check_identity()` compares three things: ClinVar's title HGVS against `--hgvs`, the
  gnomAD id's chromosome and position against VEP's resolved locus, and the gnomAD alleles
  against VEP's `allele_string` (complemented when `strand == -1`). That covers a wrong
  variation ID, a wrong coordinate, and a wrong alt allele.
  
  It cannot catch a wrong ClinVar id when `--hgvs` is omitted. Pass all four identifiers or
  the guard is blind on that axis.
  
  ### After editing either script, prove them this way
  
  A happy-path run exercises no guard — it only proves the endpoints answered. Run three
  payloads and check each one:
  
  | Case | Expect |
  |---|---|
  | all identifiers correct | `consistency.ok = true`; codes print normally |
  | wrong `--clinvar`, rest correct | `ok = false`, stderr `WRONG VARIATION ID` |
  | wrong `--gnomad` position | `ok = false`, stderr `DIFFERENT POSITION`, and PM2 withheld with an explanatory note |
  
  Then `python3 -m py_compile` both scripts and re-run `--selftest`. Allow 1–2 minutes per
  payload: the elapsed time is VEP retries, not a hang.
  
  ## The 7 derivable codes
  
  | Code | Direction | Strength | Source | Fires when |
  |---|---|---|---|---|
  | BA1 | Benign | Standalone | gnomAD | grpmax AF > 0.05 → curation ends immediately |
  | BS1 | Benign | Strong | gnomAD | AF > gene-specific prevalence (default 0.01) |
  | PM2 | Pathogenic | Supporting | gnomAD | Absent (AF = 0). "Ultra-rare" needs a declared threshold |
  | PP3 | Pathogenic | Supporting | VEP | ≥3 available predictors damaging, none tolerated |
  | BP4 | Benign | Supporting | VEP | ≥3 available predictors tolerated, none damaging |
  | PS3 | Pathogenic | Strong | literature | Functional assay shows loss of function (human-read) |
  | BS3 | Benign | Strong | literature | Functional assay shows no effect (human-read) |
  
  Predictor thresholds: REVEL ≥ 0.5 damaging / < 0.5 tolerated · CADD PHRED ≥ 20 / < 20 ·
  BayesDel addAF ≥ 0.13 damaging / ≤ −0.36 tolerated (between = indeterminate, votes neither) ·
  PolyPhen D or P damaging / B benign · SIFT deleterious / tolerated.
  
  Two PP3 schemes:
  - `count3` — project default, ≥3 of available tools damaging. Works without a calibrated
    meta-predictor; can fire when REVEL alone sits in the middle.
  - `revel` — single calibrated predictor, PP3 at REVEL ≥ 0.75, BP4 at ≤ 0.15, middle is
    indeterminate. ClinGen SVI 2022 preference; stricter.
  
  Decision logic — **points, never a code count** (ClinGen SVI / Tavtigian 2020 weights):
  Very Strong 8, Strong 4, Moderate 2, Supporting 1. BA1 → benign, done. Pathogenic at
  ≥ 10 points, benign at ≥ 7, otherwise **VUS**. Direction is reported separately and is
  not a classification.
  
  The seven public-data codes max out at 6 pathogenic points (PM2 1 + PP3 1 + PS3 4), so
  this skill **cannot** reach a pathogenic call alone. Do not "improve" the lean to fix
  that: the gap is closed by case data (PS4) or a read assay (PS3), not by re-weighting.
  When the codes and ClinVar point different ways, print `CONFLICT WITH CLINVAR` and do not
  resolve it in the tool's favour.
  
  ## Pitfalls (all verified against the live APIs)
  
  - **gnomAD v4 `joint` has no `af` field.** Compute `ac/an`. Use `dataset: gnomad_r4`;
    `popmax` lives under `faf95`, not on the block. A bare `variant(id:...)` query fails —
    the argument is `variantId` plus a required `dataset`.
  - **A differing AF is usually a different subset, not an older release.** The v4 payload
    carries `joint`, `exome` and `genome` blocks: D393V reads 2/1,613,494 = 1.24e-06 under
    `joint` and 2/1,461,238 = 1.37e-06 under `exome` — same release, same two carriers. Check
    a suspect number against each block before blaming version drift, and when it matches one
    exactly, the error was the subset. Curate from `joint` (the superset), and record the
    block beside the dataset and date; an unlabelled AF is unreproducible even when correct.
  - **VEP flags are camelCase and opt-in.** Lowercase flags are silently ignored. BayesDel is
    not a VEP plugin — it needs `dbNSFP=BayesDel_addAF_score`, and that param suppresses
    AlphaMissense on the same request, so fetch it in a second call and merge.
  - **VEP coverage is per transcript × per position.** PolyPhen returning nothing for one
    codon is normal, not a command error. Filter with `transcript_id=` to the MANE Select
    isoform; without it a variant spans 100+ transcripts.
  - **PubMed has a precision cliff.** `"TP53" AND (...) OR (...)` returned 18,103 hits of noise;
    tight `[tiab]` phrase queries over the exact protein change returned 0. Zero is a result:
    it means no one has published an assay. Never widen to gene-only and automate the reading.
  - **NCBI rate-limits at 3 requests/second** without an API key. Batch jobs need sleep or a key.
  - **VEP flaps.** Identical URLs return 200, 500 and 503 within seconds — measured at
    roughly one success in two. `_get()` skips retries on 4xx (a real answer) but retries
    5xx six times with jitter for VEP, so a run takes 1–2 min. A 500 is not a bad parameter.
  - **A typo'd gnomAD id is indistinguishable from a rare variant.** `variant(variantId:)`
    answers `Variant not found` for a wrong position *and* for a genuinely unobserved allele.
    PM2 is pathogenic evidence, so it may only fire on an id that `consistency.ok` passed.
    Worked example: `17-7676154-C-T` looks plausible, is 1.9 kb off, resolved to "absent",
    and manufactured a PM2 that the real variant (G245S, AC=5, AF 3.1e-06) does not deserve.
  - **VEP `allele_string` is transcript-oriented.** On a minus-strand gene it is the reverse
    complement of the gnomAD id: `c.1178A>T` is VEP `A/T` at 7669613 but gnomAD `T-A`.
    Complement when `strand == -1`, or every TP53 variant reads as a mismatch.
  - **`dbNSFP` must not ride on the AlphaMissense request.** Same URL = HTTP 500 and, worse,
    AlphaMissense silently disappears from the merged result. Two URLs, then merge.
  - **A swallowed dbNSFP failure reads as "no BayesDel coverage".** A 5xx on the flaky merge
    call leaves BayesDel `None`, which looks exactly like a position VEP has no dbNSFP data
    for. Since count3 counts only tools that returned data, this silently shrinks the
    denominator and can flip PP3/BP4. The payload now carries `vep.bayesdel_merge`
    (`ok` / `not returned` / `merge failed: ...`) — check it before quoting the "n/available"
    count, and re-run the fetch if it says the merge failed.
  - **A guard that compares nothing passes everything.** The identity check originally parsed
    only `c.` HGVS, so a ClinVar record filed as `NC_000017.11:g.7579817G>C` left
    `clinvar_hgvs: null` and returned `ok: true` without comparing a single value. Any "no
    problems found" result must be backed by a list of checks actually performed — the payload
    now carries `consistency.checks` and reports `no cross-check was possible` when empty.
  - **A wrong MANE transcript silently empties VEP.** Filtering `--transcript` to another
    gene's accession drops every consequence: `most_severe_consequence` becomes `"?"`, scores
    are `{}`, and downstream that reads as "this position has no predictors" instead of "you
    gave me the wrong transcript". The payload flags `vep.transcript_mismatch`. A one-character
    typo in an accession otherwise produces a confident wrong answer with no error.
  - **CNV and SV records live in ClinVar `obj_type`, not `variant_type`.** `variant_type`
    is null in esummary, so a check against it never fires. Match `copy number gain`,
    `copy number loss`, `inversion`, `translocation`, `complex` and abort with exit 2.
  - **Compare parsed HGVS coordinates, not raw strings.** ClinVar writes
    `NM_000546.6(TP53):c.733G>A` — accession outside the parens — so a string compare marks
    identical variants as mismatched and blocks a correct run.
  - **Recompute every code from the numbers printed beside it.** Hand-written evidence
    tables routinely pair a code with a value that contradicts its own threshold — a variant
    with 2 alleles in gnomAD is not "absent", so PM2 cannot fire under an AF = 0 threshold.
    Either drop the code or declare an ultra-rare cutoff; never carry both. When auditing an
    existing table, check each code against its threshold and its printed value *before*
    adding new evidence, and report the contradiction rather than silently adopting the
    table's lean.
  
  ## Deliverable
  
  An evidence matrix, one row per variant, columns: variant · ClinVar classification and
  review status · gnomAD AF, grpmax, homozygotes · predictor scores and how many called
  damaging · functional evidence · codes applied · lean. Every cell traceable to a fetched
  value with source and date.
  
  Methods sentence to carry: *"Variants were annotated against ClinVar, gnomAD (dataset
  gnomad_r4, accessed YYYY-MM-DD) and Ensembl VEP (transcript NM_…/ENST…). In-silico evidence
  was assessed on algorithms with annotation for that position; algorithms without data were
  excluded and recorded as NA. PP3/BP4 was applied on the {count3|revel} scheme."*
  
  ## Files
  
  - `scripts/fetch_evidence.py` — all four sources, stdlib only, NA never becomes 0
  - `scripts/acmg_codes.py` — code derivation, auditable basis line per code
  - `scripts/validate_skill.py` — Agent Skills spec conformance + host compatibility table
  - `tests/test_engine.py` — offline: does PM2 still refuse an unverified id? (`python3 tests/test_engine.py`)
  - `references/thresholds.md` — full tables, scheme comparison, two worked examples
    (D393V: scheme divergence; G245S: the identity trap that manufactures PM2)
  - `examples/g245s.json` — a real fetched payload, used as the test fixture
  
  Published as a standalone repo: https://github.com/darshak26/acmg-variant-curation
  (version in `metadata.version`; see CHANGELOG.md — threshold or guard changes are breaking).
  

• README.md17 KB5260e494d173upstream ↗

  # acmg-variant-curation
  
  **Derive ACMG/AMP 2015 codes for a germline variant from ClinVar, gnomAD, Ensembl VEP and
  PubMed — with the number that triggered every code printed beside it.**
  
  [![CI](https://github.com/darshak26/acmg-variant-curation/actions/workflows/validate.yml/badge.svg)](https://github.com/darshak26/acmg-variant-curation/actions/workflows/validate.yml)
  [![Release](https://img.shields.io/github/v/release/darshak26/acmg-variant-curation)](https://github.com/darshak26/acmg-variant-curation/releases)
  [![License: MIT](https://img.shields.io/badge/license-MIT-yellow.svg)](LICENSE)
  [![Python 3.10+](https://img.shields.io/badge/python-3.10%2B-blue.svg)](https://www.python.org/)
  [![Agent Skills](https://img.shields.io/badge/Agent%20Skills-spec%20compliant-6f42c1.svg)](https://agentskills.io/specification)
  
  An [Agent Skills](https://agentskills.io/specification) package: works as a skill in Claude Code,
  Codex CLI, Hermes and Antigravity. Stdlib-only Python, no install step, no API keys.
  
  > **This is evidence assembly, not a clinical classification.** It has not been validated against
  > a ClinGen expert-panel set. See [Limitations](#limitations) first.
  
  ---
  
  An LLM asked "what is the ACMG classification of this variant" answers from training data. The
  number sounds precise and cannot be checked. This tool converts that into a question you *can*
  check: **did every number come from a live query, and does it clear the threshold printed beside
  it?**
  
  Two properties make that enforceable rather than aspirational:
  
  - Every code prints its basis — `PP3 · 4/4 available tools damaging (REVEL, CADD, BayesDel, SIFT)`.
  - A **missing** predictor is `NA`, excluded from numerator *and* denominator. It never votes
    damaging, never votes tolerated. NA is not zero.
  
  ## Contents
  
  - [Installation](#installation)
  - [Quick start](#quick-start)
  - [CLI reference](#cli-reference)
  - [How it works](#how-it-works)
  - [ACMG codes](#acmg-codes)
  - [The lean](#the-lean)
  - [The identifier guard](#the-identifier-guard)
  - [Compatibility](#compatibility)
  - [Project structure](#project-structure)
  - [Development](#development)
  - [Contributing](#contributing)
  - [Limitations](#limitations)
  - [Versioning](#versioning)
  - [Citation](#citation)
  - [License](#license)
  
  ## Installation
  
  Requires **Python 3.10+** and outbound HTTPS. Nothing else — no `pip install`, no API key, no
  config file.
  
  The repository root *is* the skill. Copy the directory to your host's skill location:
  
  ```bash
  git clone https://github.com/darshak26/acmg-variant-curation.git
  cp -r acmg-variant-curation ~/.claude/skills/        # Claude Code
  ```
  
  | Host | Skill location | Invoke |
  |---|---|---|
  | Claude Code | `~/.claude/skills/acmg-variant-curation/` or `.claude/skills/…` in a repo | `/acmg-variant-curation` |
  | Codex CLI | `~/.agents/skills/acmg-variant-curation/` — **not** `~/.codex/skills` | `/skills` or `$acmg-variant-curation` |
  | Hermes | `~/.hermes/profiles/<profile>/skills/research/acmg-variant-curation/` | `/acmg-variant-curation` |
  | Antigravity | inside a plugin bundle: `agy plugin import claude` | `/skills` in-session |
  
  Not using an agent? Run the scripts directly — they are standalone CLI tools.
  
  > **Codex CLI note:** its default sandbox is `workspace-write` with network access **off**. Every
  > step here is a live HTTPS call, so the skill does nothing until you enable network for the
  > sandbox.
  
  ## Quick start
  
  ```bash
  # 0. confirm all four endpoints answer before trusting a batch
  python3 scripts/fetch_evidence.py --selftest
  
  # 1-3. ClinVar + gnomAD + VEP for one variant (three identifiers, same variant)
  python3 scripts/fetch_evidence.py \
      --clinvar 12365 \
      --gnomad 17-7674230-C-T \
      --hgvs "NM_000546.6:c.733G>A" --transcript ENST00000269305 \
      --pubmed-terms "p.Gly245Ser,G245S" \
      --out G245S.json
  
  # 4. read PubMed by hand — the payload's action_required explains why
  
  # 5. derive codes
  python3 scripts/acmg_codes.py G245S.json --scheme count3
  ```
  
  Real output for TP53 p.Gly245Ser, from [`examples/g245s.json`](examples/g245s.json):
  
  ```
  Variant: NM_000546.6(TP53):c.733G>A (p.Gly245Ser)
  ClinVar: Pathogenic (reviewed by expert panel)
  gnomAD gnomad_r4: AF=3.098419434278186e-06 homozygotes=0 faf95_popmax=6.8e-07 (nfe)
  
  Predictor calls:
    REVEL     0.947        damaging
    CADD      27.5         damaging
    BayesDel  0.550935     damaging
    SIFT      deleterious  damaging
  
  Codes:
    PP3  pathogenic supporting  4/4 available tools damaging (REVEL, CADD, BayesDel, SIFT)
  
  Points (ClinGen SVI): pathogenic 1 [PP3 1]  |  benign 0
  Lean: VUS - insufficient evidence (1 pathogenic / 0 benign points).
        Direction pathogenic, which is not a classification
  Verdict: VUS
  ```
  
  Read those last three lines first. **This variant is Pathogenic by ClinGen expert panel and the tool
  still says VUS** — deliberately. Richards 2015 has no combination rule under which one Supporting
  code classifies anything. Calling it "pathogenic-leaning" because one code outnumbered zero would
  be flattering, not accurate. The tool reports direction, states plainly that direction is not a
  classification, and names what would close the gap: PS4 case data or a read PS3 assay, neither of
  which is derivable from these databases.
  
  Two details also worth reading. **PM2 did not fire** — the variant is in gnomAD (AC=5), so under a
  strict "absent from gnomAD" reading it cannot, and the tool says so instead of quietly rounding five
  alleles down to zero. And the count is **4/4, not 3/4** — BayesDel comes from a second VEP call, and
  the payload records `vep.bayesdel_merge: "ok"`, because a failed merge there would shrink the
  denominator and could flip the call.
  
  ## CLI reference
  
  ### `scripts/fetch_evidence.py`
  
  | Flag | Required | Meaning |
  |---|---|---|
  | `--clinvar` | one of these | ClinVar variation ID, e.g. `12365` |
  | `--gnomad` | one of these | gnomAD v4 ID, GRCh38: `chr-pos-ref-alt`, e.g. `17-7674230-C-T` |
  | `--hgvs` | one of these | HGVS with accession — `NM_000546.6:c.733G>A` or `NC_000017.11:g.7674230C>T` |
  | `--transcript` | no | MANE Select ENST, filters VEP to one isoform |
  | `--pubmed-terms` | no | Comma-separated protein-change spellings, e.g. `"p.Gly245Ser,G245S"` |
  | `--out` | no | JSON output path (default: stdout) |
  | `--selftest` | no | Check all four endpoints and exit |
  
  Exit codes: `0` success, `2` out of scope (CNV/SV).
  
  ### `scripts/acmg_codes.py`
  
  | Flag | Default | Meaning |
  |---|---|---|
  | `jsonfile` | — | Payload from `fetch_evidence.py` |
  | `--scheme` | `count3` | `count3` or `revel` — pick one, declare it, never mix |
  | `--bs1` | `0.01` | BS1 frequency threshold; override with gene-specific prevalence |
  | `--ba1` | `0.05` | BA1 standalone-frequency threshold |
  | `--pm2-max-af` | `0.0` | PM2 "absent" threshold. Raising it is a methodological choice — declare it |
  | `--functional` | `none` | `damaging` / `benign`, from an assay **you read**. Never inferred |
  
  ## How it works
  
  Order matters — population frequency can end curation early (BA1), and functional evidence outranks
  every in-silico predictor.
  
  | Step | Source | Yields |
  |---|---|---|
  | 1 | ClinVar (NCBI E-utilities) | classification, review status, last evaluated, variant type |
  | 2 | gnomAD v4 (GraphQL) | AF, grpmax faf95, homozygotes, absence |
  | 3 | Ensembl VEP (REST) | REVEL, CADD, BayesDel, PolyPhen, SIFT, AlphaMissense |
  | 4 | PubMed (E-utilities) | whether a functional assay exists at all |
  | 5 | `acmg_codes.py` | the codes, each with its basis |
  
  **Binaries, not just variants.** `PS3`/`BS3` are never inferred. PubMed returns candidate papers and
  you read them; a title is not a functional assay.
  
  ## ACMG codes
  
  Seven codes are derivable from public data:
  
  | Code | Direction | Strength | Source | Fires when |
  |---|---|---|---|---|
  | BA1 | Benign | Standalone | gnomAD | grpmax AF > 0.05 — curation ends |
  | BS1 | Benign | Strong | gnomAD | AF > prevalence (default 0.01) |
  | PM2 | Pathogenic | Supporting | gnomAD | absent, **and the identifier verified** |
  | PP3 | Pathogenic | Supporting | VEP | ≥3 available predictors damaging |
  | BP4 | Benign | Supporting | VEP | ≥3 available predictors tolerated |
  | PS3 | Pathogenic | Strong | literature | assay shows loss of function (human-read) |
  | BS3 | Benign | Strong | literature | assay shows no effect (human-read) |
  
  Predictor thresholds: REVEL ≥ 0.5 · CADD PHRED ≥ 20 · BayesDel addAF ≥ 0.13 (≤ −0.36 tolerated,
  between = indeterminate, votes neither) · PolyPhen D/P · SIFT deleterious.
  
  **Not covered, deliberately:** PS4, PM6, PP1 and PS2 need case-level data this pipeline cannot see.
  CNVs, structural variants and somatic curation are out of scope — the PP3/BP4 machinery means
  nothing there. Clinical use needs a ClinGen expert panel.
  
  Out-of-scope input is **rejected, not half-curated** — exit code 2, the variant type and size, and a
  pointer to the right standard:
  
  ```
  [OUT OF SCOPE] Variant 150702 is a copy number gain, ~126 kb: GRCh38/hg38 17p13.1(chr17:7544902-7670581)x3
  acmg-variant-curation curates germline SNVs and small indels under Richards 2015.
  For CNVs use the ACMG/ClinGen 2020 technical standards (Riggs et al., Genet Med 22:245-257) ...
  ```
  
  ### Two in-silico schemes
  
  PP3 has two defensible implementations that disagree on real variants. Pick one, declare it in
  Methods, never mix.
  
  - **`count3`** (default) — PP3 when ≥3 of the *available* predictors call damaging. Works without a
    calibrated meta-predictor.
  - **`revel`** — single calibrated predictor: PP3 at REVEL ≥ 0.75, BP4 at ≤ 0.15, middle is
    indeterminate. The stricter ClinGen SVI 2022 reading.
  
  ## The lean
  
  Codes carry direction *and strength*. Strength converts to ClinGen SVI points (Tavtigian et al.
  2020): Very Strong 8, Strong 4, Moderate 2, Supporting 1. **Pathogenic at ≥ 10, benign at ≥ 7**,
  otherwise VUS. BA1 stays standalone.
  
  **The ceiling is a design property, not a bug.** The seven public-data codes max out at **6 points**
  (PM2 1 + PP3 1 + PS3 4), so this tool cannot reach a pathogenic classification on its own. A test
  pins that so it can't be "fixed" into over-calling. Expect `VUS` with a direction on almost
  everything — that is the honest answer for public data.
  
  This replaced a code-count lean that called any variant with more pathogenic codes than benign ones
  "pathogenic-leaning", which meant a variant merely absent from gnomAD and nothing else read as
  pathogenic — synonymous, intronic and UTR variants alike. When the codes and ClinVar disagree, the
  conflict is printed as `CONFLICT WITH CLINVAR` rather than resolved silently in the tool's favour.
  
  ## The identifier guard
  
  The failure this tool exists to prevent is not a wrong score — it is a **right score attached to the
  wrong variant**. ClinVar will happily return a pathogenic classification for variant A while gnomAD
  and VEP answer for variant B, and every code still fires.
  
  `fetch_evidence.py` cross-checks the identifiers and emits a `consistency` block:
  
  ```json
  {
    "ok": false,
    "problems": ["gnomAD id is at 17:7676154 but VEP resolved 17:7674230 -> DIFFERENT POSITION"],
    "checks": ["requested HGVS == ClinVar HGVS (c.733G>A)"],
    "identifiers": { "gnomad_id": "17-7676154-C-T", "vep_locus": "17:7674230 G/A" }
  }
  ```
  
  `acmg_codes.py` then **refuses PM2** on that payload. This matters because gnomAD answers
  `Variant not found` identically for a typo and for a genuinely ultra-rare allele — absence alone
  proves nothing about a variant nobody verified.
  
  A guard that compares nothing must not report success. If no cross-check is possible, `ok` is
  `false` and PM2 is withheld; `consistency.checks` lists the comparisons actually performed.
  
  ## Compatibility
  
  | Host | Loads | Notes |
  |---|---|---|
  | Claude Code | ✅ | Follows the Agent Skills spec; loads from personal, project or plugin level |
  | Codex CLI | ✅ | `~/.agents/skills/`; needs sandbox network access enabled |
  | Hermes | ✅ | Auto-registers `/<name>` as a slash command |
  | Antigravity | ✅ | Via a Claude-format plugin bundle |
  
  Frontmatter uses only the six spec fields (`name`, `description`, `license`, `compatibility`,
  `metadata`, `allowed-tools`). No host-specific extensions — those hard-error on some packaging
  paths. Paths inside `SKILL.md` are relative to the skill root.
  
  Run `python3 scripts/validate_skill.py` to check this yourself.
  
  ## Project structure
  
  ```
  acmg-variant-curation/
  ├── SKILL.md                  skill definition: rules, codes, pitfalls, deliverable
  ├── README.md                 this file
  ├── CHANGELOG.md              version history; threshold/guard changes are breaking
  ├── CITATION.cff              machine-readable citation
  ├── LICENSE · NOTICE          MIT + scope disclaimer
  ├── scripts/
  │   ├── fetch_evidence.py     ClinVar + gnomAD + VEP + PubMed → one payload
  │   ├── acmg_codes.py         payload → codes, points, lean
  │   └── validate_skill.py     spec conformance + host compatibility
  ├── references/
  │   └── thresholds.md         full tables, scheme comparison, worked examples
  ├── examples/
  │   └── g245s.json            real fetched payload, used as test fixture
  ├── tests/
  │   └── test_engine.py        offline assertions, no network
  └── .github/workflows/
      └── validate.yml          CI: Python 3.10 + 3.12
  ```
  
  ## Development
  
  ```bash
  python3 scripts/validate_skill.py              # spec conformance + compatibility table
  python3 tests/test_engine.py                   # 19 offline assertions, no network
  python3 scripts/fetch_evidence.py --selftest   # network: are all four endpoints up?
  ```
  
  CI runs the first two on Python 3.10 and 3.12. The live-endpoint check is manual
  (`workflow_dispatch`) because Ensembl VEP flaps roughly half the time and would make CI red for
  reasons unrelated to the code.
  
  **After editing either script, prove them with three payloads**, not one — a happy path exercises
  no guard:
  
  | Case | Expect |
  |---|---|
  | all identifiers correct | `consistency.ok = true`, codes print normally |
  | wrong `--clinvar` | `ok = false`, `WRONG VARIATION ID` on stderr |
  | wrong `--gnomad` position | `ok = false`, `DIFFERENT POSITION`, PM2 withheld |
  
  Allow 1–2 minutes per payload: the elapsed time is VEP retries, not a hang.
  
  ## Contributing
  
  Issues and pull requests are welcome. Please:
  
  1. **Add a test** for any change to the guard or the lean — those are the behaviours that matter.
  2. **Keep the scripts stdlib-only.** CI fails on a non-stdlib import; `pip install` defeats the
     point of a portable skill.
  3. **Never widen a threshold without a citation** in `references/thresholds.md`.
  4. Update `CHANGELOG.md` and bump `metadata.version`. A change to thresholds, schemes or guard
     behaviour is a **breaking** change.
  
  ## Limitations
  
  - **Not validated against expert calls.** Exercised on a handful of TP53 variants. No concordance
    study against a ClinGen panel-classified set has been run. Treat output as a draft for review, not
    a result.
  - **Cannot reach a pathogenic call.** The seven codes max out at 6 of the 10 points required. This
    is deliberate; see [The lean](#the-lean).
  - **PS3/BS3 are never inferred.** You read the papers.
  - **Ensembl VEP flaps.** Identical requests return 200, 500 and 503 within seconds — roughly one
    success in two, measured 2026-10-26. The client retries 5xx six times with jitter, so a run takes
    1–2 minutes. A 500 is not a bad parameter.
  - **Public data only.** No somatic, no CNV/SV, no phenotype, no segregation, no case-level data.
  - **Codes can change between versions.** The lean and guard are part of the contract; a v1.0.0 and
    a v2.0.0 run on the same variant may disagree.
  
  ## Versioning
  
  [Semantic versioning](https://semver.org/). The version lives in `metadata.version` in `SKILL.md`
  and in [CHANGELOG.md](CHANGELOG.md). Pin it in your write-up alongside the dataset and access date:
  
  > Variants were annotated against ClinVar, gnomAD (dataset `gnomad_r4`, accessed YYYY-MM-DD) and
  > Ensembl VEP (transcript NM_000546.6/ENST00000269305). In-silico evidence was assessed on
  > algorithms with annotation for that position; algorithms without data were excluded and recorded
  > as NA. PP3/BP4 was applied on the `count3` scheme. Lean on ClinGen SVI points (Tavtigian et al.
  > 2020). Tool: acmg-variant-curation v2.1.0.
  
  ## Citation
  
  > Padmane S. *acmg-variant-curation* (v2.1.0): ACMG code derivation from public variant databases
  > as an agent skill. https://github.com/darshak26/acmg-variant-curation
  
  Machine-readable metadata in [`CITATION.cff`](CITATION.cff).
  
  ## License
  
  [MIT](LICENSE) — see also [NOTICE](NOTICE) for the scope disclaimer.
  
  ## Acknowledgements
  
  ACMG/AMP 2015 criteria: Richards et al., *Genet Med* 17:405–424. Points weighting: Tavtigian et al.,
  *Genet Med* 22:1653–1657 (2020). CNV standards cited in out-of-scope messages: Riggs et al.,
  *Genet Med* 22:245–257 (2020).
  
  Data from [ClinVar](https://www.ncbi.nlm.nih.gov/clinvar/), [gnomAD](https://gnomad.broadinstitute.org/),
  [Ensembl VEP](https://www.ensembl.org/Tools/VEP) and [PubMed](https://pubmed.ncbi.nlm.nih.gov/).
  Not affiliated with or endorsed by any of them.
  

• scripts/acmg_codes.py13 KB360ad618a1aeupstream ↗

  #!/usr/bin/env python3
  """Apply the 7-code ACMG/AMP 2015 public-database subset to one fetched variant.
  
  Input is the JSON from fetch_evidence.py. Output is an auditable table: every code
  carries the number that triggered it, so a reader can recheck the call by hand.
  
  Two in-silico schemes, both defensible, and they can disagree:
    count3  - PP3 when >=3 of the AVAILABLE predictors call damaging (project default).
              NA is excluded from numerator AND denominator; 'indeterminate' scores
              vote neither way.
    revel   - PP3 at REVEL >= 0.75, BP4 at REVEL <= 0.15 (single calibrated predictor).
  Declare which scheme was used in Methods. Do not mix them.
  
  PS3/BS3 are never inferred by this script. Pass --functional from a paper you read.
  
  Usage:
    python3 acmg_codes.py D393V.json --scheme count3
    python3 acmg_codes.py D393V.json --special-bs1 0.001 --functional none
  """
  import argparse, json, sys
  
  DAMAGING = "pathogenic"
  BENIGN = "benign"
  
  # ClinGen SVI points weighting (Tavtigian et al. 2020). The vote-count lean this
  # replaced said "pathogenic-leaning" for one Supporting code, which no ACMG
  # combination rule supports: a single PP/PM code cannot classify anything.
  STRENGTH_POINTS = {"very_strong": 8, "strong": 4, "moderate": 2, "supporting": 1}
  PATHOGENIC_THRESHOLD = 10       # Tavtigian: >=10 points -> pathogenic
  BENIGN_THRESHOLD = 7            # >=7 benign points -> benign
  
  
  def _clinvar_direction(label):
      """'Pathogenic/Likely pathogenic' -> pathogenic, 'Benign' -> benign, else None."""
      s = (label or "").lower()
      if "conflict" in s or not s:
          return None
      has_p, has_b = "pathogenic" in s, "benign" in s
      if has_p == has_b:
          return None
      return "pathogenic" if has_p else "benign"
  
  
  def _polyphen_dir(label):
      if not label:
          return None
      l = label.lower()
      if any(k in l for k in ("damaging", "deleterious")):
          return "damaging"
      if "benign" in l or "tolerated" in l:
          return "tolerated"
      return "indeterminate"
  
  
  def _sift_dir(label):
      if not label:
          return None
      l = label.lower()
      if "deleter" in l:
          return "damaging"
      if "tolerat" in l:
          return "tolerated"
      return "indeterminate"
  
  
  def classify_predictors(scores):
      """Return per-tool direction and the three buckets. Missing tool -> None, not 'tolerated'."""
      rows, damaging, tolerated, indeterminate = [], [], [], []
  
      def add(tool, value, direction):
          rows.append({"tool": tool, "value": value, "call": direction})
          if direction == "damaging":
              damaging.append(tool)
          elif direction == "tolerated":
              tolerated.append(tool)
          elif direction == "indeterminate":
              indeterminate.append(tool)
  
      v = scores.get("revel")
      if v is not None:
          add("REVEL", v, "damaging" if v >= 0.5 else "tolerated")
      v = scores.get("cadd_phred")
      if v is not None:
          add("CADD", v, "damaging" if v >= 20 else "tolerated")
      v = scores.get("bayesdel_addaf")
      if v is not None:
          add("BayesDel", v, "damaging" if v >= 0.13 else ("tolerated" if v <= -0.36 else "indeterminate"))
      d = _polyphen_dir(scores.get("polyphen"))
      if d:
          add("PolyPhen", scores.get("polyphen"), d)
      d = _sift_dir(scores.get("sift"))
      if d:
          add("SIFT", scores.get("sift"), d)
      return rows, damaging, tolerated, indeterminate
  
  
  def main():
      ap = argparse.ArgumentParser()
      ap.add_argument("jsonfile")
      ap.add_argument("--scheme", choices=["count3", "revel"], default="count3")
      ap.add_argument("--bs1", type=float, default=0.01,
                      help="BS1 frequency threshold; override with gene-specific prevalence")
      ap.add_argument("--ba1", type=float, default=0.05)
      ap.add_argument("--functional", choices=["none", "damaging", "benign"], default="none",
                      help="from a functional assay you read, never inferred")
      ap.add_argument("--pm2-max-af", type=float, default=0.0,
                      help="PM2 'absent' threshold. 0.0 = strictly absent. Raising it is a "
                           "methodological choice and must be declared.")
      a = ap.parse_args()
      d = json.load(open(a.jsonfile, encoding="utf-8"))
  
      g = d.get("gnomad") or {}
      vep_obj = d.get("vep") or {}
      vep = vep_obj.get("scores") or {}
      af = g.get("af")
      faf95 = g.get("faf95_popmax")
      freq_for_ba1 = faf95 if faf95 is not None else af
  
      codes, notes = [], []
  
      # --- frequency codes
      if freq_for_ba1 is not None and freq_for_ba1 > a.ba1:
          codes.append({"code": "BA1", "direction": BENIGN, "strength": "standalone",
                        "basis": f"grpmax AF {freq_for_ba1:.3g} > {a.ba1}",
                        "dataset": g.get("dataset")})
      if af is not None and af > a.bs1:
          codes.append({"code": "BS1", "direction": BENIGN, "strength": "strong",
                        "basis": f"AF {af:.3g} > {a.bs1}", "dataset": g.get("dataset")})
      # gnomAD answers "Variant not found" identically for a typo and for a genuinely
      # ultra-rare allele. PM2 is a pathogenic code, so it may not rest on an
      # identifier nobody verified: require consistency.ok when the payload carries it.
      cons = d.get("consistency") or {}
      identity_ok = cons.get("ok", True)
      if (g.get("absent") or (af is not None and af <= a.pm2_max_af)) and not identity_ok:
          notes.append(
              "PM2 NOT applied: the gnomAD/ClinVar/VEP identifiers disagree "
              f"({'; '.join(cons.get('problems') or [])}). An unverified id returning "
              "'Variant not found' is indistinguishable from a typo, and a typo must "
              "never become pathogenic evidence. Re-run fetch with the correct id.")
      elif g.get("absent") or (af is not None and af <= a.pm2_max_af):
          basis = ("not found in " + str(g.get("dataset"))) if g.get("absent") else \
                  f"AF {af:.3g} <= {a.pm2_max_af}"
          codes.append({"code": "PM2", "direction": DAMAGING, "strength": "supporting",
                        "basis": basis, "dataset": g.get("dataset")})
      if af is not None and af > a.pm2_max_af and not g.get("absent"):
          notes.append(
              f"PM2 NOT applied: variant is present in {g.get('dataset')} "
              f"(ac={g.get('joint', {}).get('ac') if isinstance(g.get('joint'), dict) else '?'}, "
              f"AF {af:.3g}) and --pm2-max-af is {a.pm2_max_af}. Under the strict "
              f"'absent from gnomAD' reading PM2 cannot fire. If the project instead uses "
              f"'ultra-rare' as the PM2 condition, declare that threshold in Methods.")
  
      # --- in-silico codes
      rows, damaging, tolerated, indet = classify_predictors(vep)
      if a.scheme == "count3":
          available = damaging + tolerated + indet
          if len(available) >= 2 and len(damaging) >= 3 and not tolerated:
              codes.append({"code": "PP3", "direction": DAMAGING, "strength": "supporting",
                            "basis": f"{len(damaging)}/{len(available)} available tools damaging "
                                     f"({', '.join(damaging)})"})
          elif len(available) >= 2 and len(tolerated) >= 3 and not damaging:
              codes.append({"code": "BP4", "direction": BENIGN, "strength": "supporting",
                            "basis": f"{len(tolerated)}/{len(available)} available tools tolerated "
                                     f"({', '.join(tolerated)})"})
          else:
              notes.append(
                  f"No PP3/BP4 on the count3 scheme: {len(damaging)} damaging, "
                  f"{len(tolerated)} tolerated, {len(indet)} indeterminate, "
                  f"{len(rows) - len(available)} missing.")
          if len(available) < 2:
              notes.append("Fewer than 2 predictors had data: ClinGen SVI says do not assign "
                           "PP3/BP4 on one tool.")
      else:
          r = vep.get("revel")
          if r is None:
              notes.append("revel scheme chosen but REVEL has no value for this position.")
          elif r >= 0.75:
              codes.append({"code": "PP3", "direction": DAMAGING, "strength": "supporting",
                            "basis": f"REVEL {r} >= 0.75"})
          elif r <= 0.15:
              codes.append({"code": "BP4", "direction": BENIGN, "strength": "supporting",
                            "basis": f"REVEL {r} <= 0.15"})
          else:
              notes.append(f"REVEL {r} sits in the indeterminate band (0.15-0.75): "
                           f"no PP3/BP4 on the calibrated scheme.")
  
      # --- functional codes (human-entered only)
      if a.functional == "damaging":
          codes.append({"code": "PS3", "direction": DAMAGING, "strength": "strong",
                        "basis": "functional assay reported a damaging effect (supplied by curator)"})
      elif a.functional == "benign":
          codes.append({"code": "BS3", "direction": BENIGN, "strength": "strong",
                        "basis": "functional assay reported no effect (supplied by curator)"})
  
      # --- lean, on ClinGen SVI points (Tavtigian 2020), never a code count
      path = [c for c in codes if c["direction"] == DAMAGING]
      ben = [c for c in codes if c["direction"] == BENIGN]
      path_pts = sum(STRENGTH_POINTS.get(c["strength"], 0) for c in path)
      ben_pts = sum(STRENGTH_POINTS.get(c["strength"], 0) for c in ben)
  
      verdict = None
      if any(c["code"] == "BA1" for c in codes):
          verdict = "BENIGN"
          lean = "Benign (standalone BA1) - curation complete, no further evidence needed"
      elif path_pts >= PATHOGENIC_THRESHOLD:
          verdict = "PATHOGENIC"
          lean = f"Pathogenic - {path_pts} pathogenic points (threshold {PATHOGENIC_THRESHOLD})"
      elif ben_pts >= BENIGN_THRESHOLD:
          verdict = "BENIGN"
          lean = f"Benign - {ben_pts} benign points (threshold {BENIGN_THRESHOLD})"
      else:
          verdict = "VUS"
          direction = ("pathogenic" if path_pts > ben_pts else
                       "benign" if ben_pts > path_pts else "none")
          lean = (f"VUS - insufficient evidence ({path_pts} pathogenic / {ben_pts} benign points). "
                  f"Direction {direction}, which is not a classification")
  
      # A public-data-only pipeline rarely reaches 10 points. Say so rather than let a
      # reader mistake "pathogenic-leaning" for a call, and flag a ClinVar clash loudly.
      if verdict != "BENIGN" and path_pts < PATHOGENIC_THRESHOLD and path_pts:
          notes.append(
              f"{path_pts} pathogenic point(s) from public data only - below the "
              f"{PATHOGENIC_THRESHOLD}-point threshold. The codes support a direction, not a "
              "classification. Case-level evidence (PS4) or a read functional assay (PS3) is "
              "what closes the gap, and neither is derivable from these databases.")
      cv_dir = _clinvar_direction((d.get("clinvar") or {}).get("classification"))
      cv_label = (d.get("clinvar") or {}).get("classification")
      if cv_dir == "benign" and path_pts:
          notes.append(
              f"CONFLICT WITH CLINVAR: ClinVar calls this '{cv_label}' but the codes reach "
              f"{path_pts} pathogenic point(s). Supporting-strength frequency/in-silico codes "
              "do not outweigh a benign classification. Read the ClinVar record before "
              "overruling it.")
      elif cv_dir == "pathogenic" and ben_pts:
          notes.append(
              f"CONFLICT WITH CLINVAR: ClinVar calls this '{cv_label}' but the codes reach "
              f"{ben_pts} benign point(s). Read the ClinVar record before overruling it.")
      if vep_obj.get("transcript_mismatch"):
          notes.insert(0, "VEP: " + vep_obj["transcript_mismatch"])
  
      print(f"Variant: {(d.get('clinvar') or {}).get('title')}")
      print(f"ClinVar: {(d.get('clinvar') or {}).get('classification')} "
            f"({(d.get('clinvar') or {}).get('review_status')})")
      print(f"gnomAD {g.get('dataset')}: AF={af} homozygotes={g.get('homozygotes')} "
            f"faf95_popmax={faf95} ({g.get('popmax_population')})")
      print("\nPredictor calls:")
      for r in rows:
          print(f"  {r['tool']:<9} {str(r['value']):<12} {r['call']}")
      print("\nCodes:")
      for c in codes:
          print(f"  {c['code']:<4} {c['direction']:<9} {c['strength']:<11} {c['basis']}")
      if not codes:
          print("  (none derivable)")
      fmt_pts = ", ".join(f"{c['code']} {STRENGTH_POINTS.get(c['strength'], 0)}" for c in path)
      fmt_ben = ", ".join(f"{c['code']} {STRENGTH_POINTS.get(c['strength'], 0)}" for c in ben)
      print(f"\nPoints (ClinGen SVI): pathogenic {path_pts}"
            f"{' [' + fmt_pts + ']' if fmt_pts else ''}  |  benign {ben_pts}"
            f"{' [' + fmt_ben + ']' if fmt_ben else ''}")
      print(f"Lean: {lean}")
      if notes:
          print("\nNotes and caveats:")
          for n in notes:
              print(f"  - {n}")
      print(f"Verdict: {verdict}")
      print(f"\nScheme: {a.scheme} | BS1>{a.bs1} | BA1>{a.ba1} | PM2 max AF {a.pm2_max_af} | "
            f"functional={a.functional}")
      print("Thresholds and dataset versions must be reported in Methods.")
  
  
  if __name__ == "__main__":
      main()
  

• scripts/fetch_evidence.py21 KB2084cfdab699upstream ↗

  #!/usr/bin/env python3
  """Fetch public evidence for one ClinVar variant: ClinVar, gnomAD, Ensembl VEP, PubMed.
  
  Stdlib only. Every source degrades to None (never 0) when it has no data,
  because 0 is a real score and NA is not evidence either way.
  
  Usage (all identifiers point at the SAME variant):
    python3 fetch_evidence.py \
        --clinvar 3634421 \
        --gnomad 17-7669613-T-A \
        --hgvs "NM_000546.6:c.1178A>T" --transcript ENST00000269305 \
        --pubmed-terms "p.Asp393Val,D393V" \
        --out D393V.json
  
  Run `--selftest` to check all four endpoints answer before trusting a batch.
  """
  import argparse, json, random, re, sys, time, urllib.error, urllib.parse, urllib.request
  from datetime import datetime, timezone
  
  UA = {"User-Agent": "hermes-athena-vus/1.0 (variant curation research)"}
  EUTILS = "https://eutils.ncbi.nlm.nih.gov/entrez/eutils"
  GNOMAD_API = "https://gnomad.broadinstitute.org/api"
  VEP_REST = "https://rest.ensembl.org"
  VEP_TRIES = 6   # Ensembl VEP returns 500/503 on ~half of requests; retry hard
  TIMEOUT = 45
  
  
  def _get(url, accept=None, tries=3):
      hdrs = dict(UA)
      if accept:
          hdrs["Accept"] = accept
      err = None
      for i in range(tries):
          try:
              req = urllib.request.Request(url, headers=hdrs)
              with urllib.request.urlopen(req, timeout=TIMEOUT) as r:
                  return r.read().decode("utf-8", "replace")
          except urllib.error.HTTPError as e:
              err = e
              # 4xx other than 429 is a real answer (bad id, bad param) - do not burn retries
              if 400 <= e.code < 500 and e.code != 429:
                  raise
          except Exception as e:
              err = e
          # Ensembl VEP flaps: identical URLs return 200/500/503 within seconds.
          backoff = min(2.0 * (i + 1), 8.0) + random.uniform(0, 1.0)
          time.sleep(backoff)
      raise RuntimeError(f"GET {url} failed after {tries} tries: {err!r}")
  
  
  def _post_json(url, payload, tries=3):
      err = None
      for i in range(tries):
          try:
              req = urllib.request.Request(
                  url, data=json.dumps(payload).encode(),
                  headers={**UA, "Content-Type": "application/json", "Accept": "application/json"})
              with urllib.request.urlopen(req, timeout=TIMEOUT) as r:
                  return json.loads(r.read().decode())
          except Exception as e:
              err = e
              body = ""
              try:
                  body = e.read().decode()[:300]  # GraphQL errors explain schema drift
              except Exception:
                  pass
              if body:
                  err = RuntimeError(f"{e!r} {body}")
              time.sleep(1.5 * (i + 1))
      raise RuntimeError(f"POST {url} failed: {err!r}")
  
  
  class OutOfScope(Exception):
      """The variant is real but outside this tool's competence. Never guess around it."""
  
  
  # ClinVar files copy-number and structural variants under obj_type, NOT variant_type
  # (which comes back null in esummary). Checking variant_type would never fire.
  OUT_OF_SCOPE_TYPES = ("copy number gain", "copy number loss", "inversion",
                        "translocation", "complex")
  
  
  # ---------------------------------------------------------------- ClinVar
  def fetch_clinvar(variation_id):
      url = f"{EUTILS}/esummary.fcgi?db=clinvar&retmode=json&id={urllib.parse.quote(str(variation_id))}"
      rec = json.loads(_get(url))["result"].get(str(variation_id), {})
      if not rec:
          return {"variation_id": variation_id, "error": "not found", "classification": None}
      gc = rec.get("germline_classification", {}) or {}
      obj_type = (rec.get("obj_type") or "").strip()
      if any(t in obj_type.lower() for t in OUT_OF_SCOPE_TYPES):
          title = rec.get("title") or ""
          span = re.search(r"chr[^:]*:([0-9]+)-([0-9]+)", title)
          size = ""
          if span:
              kb = (int(span.group(2)) - int(span.group(1))) / 1000
              size = f", ~{kb:,.0f} kb" if kb < 10000 else f", ~{kb/1000:,.1f} Mb"
          raise OutOfScope(
              f"Variant {variation_id} is a {obj_type}{size}: {title}\n"
              "acmg-variant-curation curates germline SNVs and small indels under "
              "Richards 2015. This variant is out of scope and no codes will be derived.\n"
              "For CNVs use the ACMG/ClinGen 2020 technical standards (Riggs et al., "
              "Genet Med 22:245-257) with a purpose-built CNV tool.")
      return {
          "variation_id": str(variation_id),
          "obj_type": obj_type,
          "accession": rec.get("accession"),
          "title": rec.get("title"),
          "gene_sort": rec.get("gene_sort"),
          "classification": gc.get("description"),
          "review_status": gc.get("review_status"),
          "last_evaluated": gc.get("last_evaluated"),
          "source": url,
      }
  
  
  # ----------------------------------------------------------------- gnomAD
  GNOMAD_Q = '''query($id:String!,$ds:DatasetId!){
    variant(variantId:$id, dataset:$ds){ variant_id rsids
      exome{ac an homozygote_count faf95{popmax popmax_population}}
      genome{ac an homozygote_count faf95{popmax popmax_population}}
      joint{ac an homozygote_count faf95{popmax popmax_population}}
    }}'''
  
  
  def fetch_gnomad(variant_id, dataset="gnomad_r4"):
      if not variant_id:
          return {"variant_id": None, "note": "no gnomAD id supplied"}
      d = _post_json(GNOMAD_API, {"query": GNOMAD_Q, "variables": {"id": variant_id, "ds": dataset}})
      v = (d.get("data") or {}).get("variant")
      if not v:
          return {"variant_id": variant_id, "dataset": dataset,
                  "absent": True, "af": None, "source": GNOMAD_API,
                  "note": "variant not found in dataset -> candidate for PM2, but confirm the id is right"}
      out = {"variant_id": v["variant_id"], "rsids": v.get("rsids") or [],
             "dataset": dataset, "source": GNOMAD_API, "absent": False}
      for pop in ("joint", "exome", "genome"):
          blk = v.get(pop)
          if not blk:
              continue
          faf = blk.get("faf95") or {}
          ac, an = blk.get("ac"), blk.get("an")
          out[pop] = {"ac": ac, "an": an,
                      "af": (ac / an) if (ac is not None and an) else None,
                      "homozygotes": blk.get("homozygote_count"),
                      "faf95_popmax": faf.get("popmax"),
                      "popmax_population": faf.get("popmax_population")}
      primary = out.get("joint") or out.get("exome") or out.get("genome") or {}
      out["af"] = primary.get("af")
      out["homozygotes"] = primary.get("homozygotes")
      out["faf95_popmax"] = primary.get("faf95_popmax")
      out["popmax_population"] = primary.get("popmax_population")
      return out
  
  
  # -------------------------------------------------------------------- VEP
  # Flags are opt-in and camelCase. Lowercase flags are silently ignored by Ensembl.
  VEP_SCORE_FLAGS = "CADD=1&REVEL=1&PolyPhen=1&SIFT=1&AlphaMissense=1"
  VEP_DBNSFP = "dbNSFP=BayesDel_addAF_score"   # dbNSFP and AlphaMissense do not return together
  
  
  def fetch_vep(hgvs, transcript=None):
      if not hgvs:
          return {"hgvs": None, "note": "no HGVS supplied"}
      # NOTE: dbNSFP is deliberately NOT on this URL. It suppresses AlphaMissense,
      # so it gets its own call below and the two results are merged.
      if transcript:
          url = f"{VEP_REST}/vep/human/hgvs/{urllib.parse.quote(hgvs, safe='')}" \
                f"?transcript_id={transcript}&{VEP_SCORE_FLAGS}&content-type=application/json"
      else:
          url = f"{VEP_REST}/vep/human/hgvs/{urllib.parse.quote(hgvs, safe='')}" \
                f"?canonical=1&{VEP_SCORE_FLAGS}&content-type=application/json"
      def _one(u):
          data = json.loads(_get(u, accept="application/json", tries=VEP_TRIES))
          return data[0] if data else {}
  
      d = _one(url)
      if not d:
          return {"hgvs": hgvs, "error": "VEP returned nothing", "scores": {}}
      tcs = d.get("transcript_consequences") or []
      tc = dict(tcs[0]) if tcs else {}
      # A wrong MANE accession makes the filter drop every consequence, which reads
      # downstream as "this position has no predictors" and quietly weakens the call.
      transcript_mismatch = None
      if transcript and not any(t.get("transcript_id") == transcript for t in tcs):
          transcript_mismatch = (
              f"no consequence returned for transcript {transcript}. Either the accession "
              "is wrong or the variant is not in that gene - the empty predictor set is the "
              "filter's doing, not missing data. Refetch with the correct MANE transcript.")
      # BayesDel is not a VEP plugin: it needs the dbNSFP param, which suppresses
      # AlphaMissense on the same request. Second call, merged in.
      # A failed merge must never look like "this position has no dbNSFP data": count3
      # counts only tools that returned a value, so a silently missing BayesDel shrinks
      # the denominator and can flip PP3/BP4. Record the outcome instead of guessing.
      bayesdel_status = "not requested (no transcript)" if not transcript else "not returned"
      if transcript:
          u2 = (f"{VEP_REST}/vep/human/hgvs/{urllib.parse.quote(hgvs, safe='')}"
                f"?transcript_id={transcript}&{VEP_DBNSFP}&content-type=application/json")
          try:
              d2 = _one(u2)
              t2 = (d2.get("transcript_consequences") or [{}])[0]
              if t2.get("bayesdel_addaf_score") is not None:
                  tc["bayesdel_addaf_score"] = t2["bayesdel_addaf_score"]
                  bayesdel_status = "ok"
          except Exception as e:                       # endpoint flapped past the retries
              bayesdel_status = f"merge failed: {str(e)[:120]}"
      scores = {
          "revel": tc.get("revel"),
          "cadd_phred": tc.get("cadd_phred"),
          "bayesdel_addaf": tc.get("bayesdel_addaf_score"),
          "polyphen": tc.get("polyphen_prediction"),   # absent == no coverage at this position
          "sift": tc.get("sift_prediction"),
          "alphamissense": (tc.get("alphamissense") or {}).get("am_class")
          if isinstance(tc.get("alphamissense"), dict) else tc.get("alphamissense"),
      }
      return {"hgvs": hgvs, "transcript_used": tc.get("transcript_id"),
              "most_severe_consequence": d.get("most_severe_consequence"),
              # carried so check_identity() can prove gnomAD and VEP are the same locus
              "seq_region_name": d.get("seq_region_name"), "start": d.get("start"),
              "allele_string": d.get("allele_string"), "strand": d.get("strand"),
              "hgvsp": d.get("hgvsp"), "hgvsc": d.get("hgvsc"),
              "scores": scores, "bayesdel_merge": bayesdel_status,
              "transcript_mismatch": transcript_mismatch,
              "tools_with_data": sorted(k for k, v in scores.items() if v is not None),
              "tools_missing": sorted(k for k, v in scores.items() if v is None),
              "source": url}
  
  
  # ----------------------------------------------------------------- PubMed
  def fetch_pubmed(terms, retmax=20):
      """Tight phrase search first. Tight query returning 0 is a RESULT, not a failure:
      it means no publication uses that protein change, so PS3/BS3 cannot be assigned."""
      if not terms:
          return {"terms": [], "note": "no PubMed terms supplied"}
      termlist = [t.strip() for t in terms.split(",") if t.strip()] or [terms]
      q = " OR ".join(f'"{t}"[tiab]' for t in termlist)
      url = f"{EUTILS}/esearch.fcgi?db=pubmed&retmode=json&retmax={retmax}&term={urllib.parse.quote(q)}"
      r = json.loads(_get(url))["esearchresult"]
      out = {"terms": termlist, "query": q, "count": int(r.get("count", 0)),
             "pmids": r.get("idlist") or [], "source": url}
      if out["pmids"]:
          ids = ",".join(out["pmids"][:retmax])
          s = json.loads(_get(f"{EUTILS}/esummary.fcgi?db=pubmed&retmode=json&id={ids}"))["result"]
          out["hits"] = [{"pmid": i, "date": s.get(i, {}).get("pubdate"),
                          "title": s.get(i, {}).get("title")} for i in out["pmids"]]
      out["action_required"] = (
          "Read abstracts by hand. PS3/BS3 may ONLY be assigned from a functional assay "
          "reported in a paper you actually read. A title is not evidence."
          if out["count"] else
          "Zero hits for the exact protein change: no PS3/BS3. Try a looser gene+domain "
          "query and read by eye before concluding nothing exists.")
      return out
  
  
  # ------------------------------------------------- identifier consistency guard
  # ClinVar records carry whichever HGVS flavour the submitter used. matching only c.
  # made the guard return ok:true WITHOUT comparing anything for genomic records -
  # a silent bypass of the one check this tool exists to perform.
  _HGVS = re.compile(
      r"((?:[A-Z]{2}_[0-9]+\.[0-9]+:)?[cgn]\.[0-9_]+"
      r"(?:[ACGTacgt]*>[ACGTacgt]+|del[ACGT]*|ins[ACGT]+|dup[ACGT]*))")
  _P_HGVS = re.compile(r"p\.\(?([A-Z][a-z]{2})(\d+)([A-Z][a-z]{2}|Ter|=)")
  
  
  def _parse_hgvs(s):
      """-> {kind, raw, pos, ref, alt} for c./g./n. notation, or None if absent."""
      if not s:
          return None
      m = _HGVS.search(s)
      if not m:
          return None
      raw = m.group(1)
      kind = re.match(r"(?:[A-Z]{2}_[0-9]+\.[0-9]+:)?([cgn])\.", raw)
      body = raw.split(":", 1)[-1] if ":" in raw else raw
      out = {"kind": kind.group(1) if kind else None, "raw": raw, "pos": None,
             "ref": None, "alt": None}
      pm = re.search(r"(\d+)([ACGTacgt]*)>([ACGTacgt]+)", body)
      if pm:
          out["pos"], out["ref"], out["alt"] = int(pm.group(1)), pm.group(2).upper(), pm.group(3).upper()
      else:
          pm = re.search(r"(\d+)", body)
          if pm:
              out["pos"] = int(pm.group(1))
      return out
  
  
  def _vep_alleles_plus_strand(vep, req_kind):
      """VEP reports allele_string in TRANSCRIPT orientation for c./n. input, and in
      plus-strand orientation for g. input. Normalise to plus strand before comparing."""
      av = (vep or {}).get("allele_string") or ""
      if "/" not in av:
          return None, None
      r, _, a = av.partition("/")
      if req_kind in ("c", "n") and (vep or {}).get("strand") == -1:
          comp = str.maketrans("ACGT", "TGCA")
          r, a = r.translate(comp), a.translate(comp)
      return r, a
  
  
  def check_identity(clinvar, gnomad, vep, hgvs=None):
      """The three identifiers must describe ONE variant.
  
      Nothing downstream can detect a mismatch: ClinVar will happily report a
      pathogenic classification for variant A while gnomAD and VEP answer for
      variant B, and every code still fires. Refuse to blend them.
      """
      problems, seen, checks = [], {}, []
      req = _parse_hgvs(hgvs)
      seen["requested_hgvs"] = req["raw"] if req else hgvs
  
      cv = None
      if isinstance(clinvar, dict) and clinvar.get("title"):
          seen["clinvar_title"] = clinvar["title"]
          cv = _parse_hgvs(clinvar["title"])
          seen["clinvar_hgvs"] = cv["raw"] if cv else None
  
      # Check 1: the HGVS you asked for vs the HGVS ClinVar filed. Same notation only.
      if req and cv:
          if req["kind"] == cv["kind"]:
              # Compare PARSED coordinates, not raw strings: ClinVar writes
              # "NM_000546.6(TP53):c.733G>A" where the accession sits outside the
              # variant parens, so a raw compare flags identical variants as wrong.
              same = None
              if req["pos"] and cv["pos"]:
                  same = (req["pos"], req["ref"], req["alt"]) == (cv["pos"], cv["ref"], cv["alt"])
              if same is None:
                  same = (req["raw"].split(":")[-1].lower() == cv["raw"].split(":")[-1].lower())
              if not same:
                  problems.append(f"ClinVar record is {cv['raw']}, you asked for "
                                  f"{req['raw']} -> WRONG VARIATION ID")
              else:
                  checks.append(f"requested HGVS == ClinVar HGVS ({req['raw'].split(':')[-1]})")
          else:
              checks.append(f"ClinVar is {cv['kind']}. and you asked for {req['kind']}. : "
                            "cross-checked by genomic position instead of by string")
  
      vep = vep if isinstance(vep, dict) else {}
      vpos = vep.get("start")
      vchr = str(vep.get("seq_region_name") or "")
      if vpos:
          seen["vep_locus"] = f"{vchr}:{vpos} {vep.get('allele_string')}"
      vref, valt = _vep_alleles_plus_strand(vep, req["kind"] if req else None)
  
      # Check 2: ClinVar genomic record vs VEP's resolution of the same variant.
      if cv and cv["kind"] == "g" and vpos:
          if cv["pos"] != vpos:
              problems.append(f"ClinVar genomic position is {cv['pos']} but VEP resolved "
                              f"{vpos} -> DIFFERENT POSITION")
          elif cv["ref"] and (vref, valt) != (cv["ref"], cv["alt"]):
              problems.append(f"ClinVar alleles {cv['ref']}>{cv['alt']} vs VEP "
                              f"{vref}>{valt} at {vpos} -> DIFFERENT ALT")
          else:
              checks.append(f"ClinVar genomic locus == VEP locus ({vchr}:{vpos})")
  
      # Check 3: the gnomAD id vs VEP's resolution.
      if isinstance(gnomad, dict):
          gid = gnomad.get("variant_id") or ""
          seen["gnomad_id"] = gid
          gm = re.match(r"^(\w+)-(\d+)-([ACGT]+)-([ACGT]+)$", gid)
          if gm and vpos:
              gchr, gpos, gref, galt = gm.groups()
              if gchr != vchr or gpos != str(vpos):
                  problems.append(f"gnomAD id is at {gchr}:{gpos} but VEP resolved "
                                  f"{vchr}:{vpos} -> DIFFERENT POSITION")
              elif vref and (gref, galt) != (vref, valt):
                  problems.append(f"gnomAD alleles {gref}>{galt} vs VEP "
                                  f"{vref}>{valt} at the same position -> DIFFERENT ALT")
              else:
                  checks.append(f"gnomAD id == VEP locus ({gchr}:{gpos} {gref}>{galt})")
  
      # Nothing was actually compared -> the identifiers are unverified, not verified.
      # ok must not be true by default: that is how the genomic-HGVS bypass happened.
      if not problems and not checks:
          problems.append(
              "no cross-check was possible: nothing tied the supplied identifiers to one "
              "variant. Supply all three (--clinvar, --gnomad, --hgvs) so the guard can "
              "prove they agree. PM2 will be withheld until it can.")
  
      return {"ok": not problems, "problems": problems, "checks": checks,
              "identifiers": seen}
  
  
  def selftest():
      print("ClinVar: ", bool(fetch_clinvar(3634421).get("classification")))
      print("gnomAD:  ", fetch_gnomad("17-7669613-T-A").get("af"))
      print("VEP:     ", fetch_vep("NM_000546.6:c.1178A>T", "ENST00000269305")["scores"]["revel"])
      print("PubMed:  ", fetch_pubmed("p.Asp393Val,D393V")["count"], "hits")
  
  
  def main():
      ap = argparse.ArgumentParser()
      ap.add_argument("--clinvar", help="ClinVar variation id (e.g. 3634421)")
      ap.add_argument("--gnomad", help="gnomAD variant id, GRCh38, e.g. 17-7669613-T-A")
      ap.add_argument("--hgvs", help="HGVS with accession: NM_000546.6:c.1178A>T or "
                                     "NC_000017.11:g.7674230C>T (both are parsed)")
      ap.add_argument("--transcript", help="MANE Select ENST id (filters VEP output)")
      ap.add_argument("--pubmed-terms", help="comma-separated protein-change spellings")
      ap.add_argument("--out", help="JSON output path (default stdout)")
      ap.add_argument("--selftest", action="store_true")
      a = ap.parse_args()
      if a.selftest:
          return selftest()
  
      result = {"fetched_at": datetime.now(timezone.utc).isoformat(timespec="seconds")}
      try:
          _collect(a, result)
      except OutOfScope as e:
          print(f"[OUT OF SCOPE] {e}", file=sys.stderr)
          return 2
  
      result["consistency"] = check_identity(
          result.get("clinvar"), result.get("gnomad"), result.get("vep"), a.hgvs)
      if not result["consistency"]["ok"]:
          print("!! IDENTIFIER CHECK FAILED - DO NOT CURATE THIS PAYLOAD", file=sys.stderr)
          for pr in result["consistency"]["problems"]:
              print("!!   " + pr, file=sys.stderr)
      else:
          for ck in result["consistency"]["checks"]:
              print(f"ok   {ck}", file=sys.stderr)
  
      txt = json.dumps(result, indent=2, ensure_ascii=False)
      if a.out:
          with open(a.out, "w", encoding="utf-8") as f:
              f.write(txt + "\n")
          print(f"wrote {a.out}")
      else:
          print(txt)
      return 0
  
  
  def _collect(a, result):
      for name, fn in (("clinvar", lambda: fetch_clinvar(a.clinvar)),
                       ("gnomad", lambda: fetch_gnomad(a.gnomad)),
                       ("vep", lambda: fetch_vep(a.hgvs, a.transcript)),
                       ("pubmed", lambda: fetch_pubmed(a.pubmed_terms))):
          if (name == "clinvar" and not a.clinvar) or (name == "gnomad" and not a.gnomad) \
             or (name == "vep" and not a.hgvs) or (name == "pubmed" and not a.pubmed_terms):
              result[name] = {"skipped": "identifier not supplied"}
              continue
          try:
              result[name] = fn()
          except OutOfScope:                          # not an endpoint failure: abort
              raise
          except Exception as e:                      # a dead endpoint is NA, never a zero
              result[name] = {"error": str(e)[:400]}
  
  
  if __name__ == "__main__":
      sys.exit(main())
  

• scripts/validate_skill.py6.8 KB37a922f2d9ecupstream ↗

  #!/usr/bin/env python3
  """Validate this skill against the Agent Skills spec and report host compatibility.
  
  Stdlib only, no PyYAML: a skill that needs `pip install` before it can be checked
  is not portable, and portability is the whole claim being tested here.
  
  Checks:
    1. YAML frontmatter parses and carries the required `name` and `description`.
    2. `name` obeys the spec grammar and matches the parent directory.
    3. `description` is 1-1024 chars.
    4. Frontmatter uses ONLY the six spec fields (warn on host-specific extras, which
       some hosts hard-error on: "Unexpected key(s) in SKILL.md frontmatter").
    5. Every script imports from the standard library only.
    6. SKILL.md stays under the recommended 500 lines.
  
  Usage:
    python3 scripts/validate_skill.py            # validate + compatibility table
    python3 scripts/validate_skill.py --quiet    # exit code only
  Exit code 0 = valid, 1 = at least one error.
  """
  import argparse
  import ast
  import pathlib
  import re
  import sys
  
  SPEC_FIELDS = {"name", "description", "license", "compatibility", "metadata", "allowed-tools"}
  
  # Fields other hosts add. Legal in Claude Code, a hard error when packaging for
  # claude.ai / the Skills API, and meaningless elsewhere.
  HOST_ONLY_FIELDS = {
      "when_to_use", "argument-hint", "arguments", "disable-model-invocation",
      "user-invocable", "disallowed-tools", "model", "effort", "context", "agent",
      "background", "hooks", "paths", "shell",
  }
  
  NAME_RE = re.compile(r"^[a-z0-9]+(-[a-z0-9]+)*$")
  HOSTS = [
      # host, where it looks, how it resolves the skill
      ("Hermes", "~/.hermes/profiles/<profile>/skills/<category>/<name>/", "auto-registers /<name>"),
      ("Claude Code", "~/.claude/skills/<name>/  or  .claude/skills/<name>/", "/<name>"),
      ("Codex CLI", "~/.agents/skills/<name>/  or  <repo>/.agents/skills/<name>/", "/skills or $<name>"),
      ("Antigravity", "agy plugin import claude  (skill inside a plugin bundle)", "/skills in-session"),
  ]
  
  
  def parse_frontmatter(text):
      """Minimal YAML subset: flat `key: value` plus one level of nested map."""
      m = re.match(r"^---\n(.*?)\n---\n", text, re.S)
      if not m:
          return None, "no YAML frontmatter block at the top of SKILL.md"
      data, parent = {}, None
      for raw in m.group(1).splitlines():
          if not raw.strip() or raw.lstrip().startswith("#"):
              continue
          if raw.startswith(("  ", "\t")):          # nested key
              if parent is None:
                  return None, f"nested value without a parent key: {raw!r}"
              k, _, v = raw.strip().partition(":")
              data[parent][k.strip()] = v.strip().strip('"').strip("'")
              continue
          k, _, v = raw.partition(":")
          k, v = k.strip(), v.strip()
          if not v:                                  # opens a nested map
              parent = k
              data[k] = {}
          else:
              parent = None
              data[k] = v.strip('"').strip("'")
      return data, None
  
  
  def stdlib_imports(paths):
      """Return {script: [non-stdlib module names]}."""
      std = set(sys.stdlib_module_names)
      bad = {}
      for p in paths:
          mods = set()
          for n in ast.walk(ast.parse(p.read_text(encoding="utf-8"))):
              if isinstance(n, ast.Import):
                  mods |= {a.name.split(".")[0] for a in n.names}
              elif isinstance(n, ast.ImportFrom) and n.level == 0 and n.module:
                  mods.add(n.module.split(".")[0])
          extra = sorted(mods - std)
          if extra:
              bad[p.name] = extra
      return bad
  
  
  def main():
      ap = argparse.ArgumentParser()
      ap.add_argument("--quiet", action="store_true")
      a = ap.parse_args()
  
      root = pathlib.Path(__file__).resolve().parent.parent
      errors, warnings, notes = [], [], []
  
      skill = root / "SKILL.md"
      if not skill.exists():
          print(f"FAIL  SKILL.md missing in {root}", file=sys.stderr)
          return 1
      text = skill.read_text(encoding="utf-8")
  
      fm, err = parse_frontmatter(text)
      if err:
          errors.append(err)
          fm = {}
  
      # 1-3. required fields and their constraints
      name, desc = fm.get("name", ""), fm.get("description", "")
      if not name:
          errors.append("required field `name` is missing")
      else:
          if not NAME_RE.fullmatch(name):
              errors.append(f"name {name!r} must be lowercase alphanumerics and single hyphens")
          if len(name) > 64:
              errors.append(f"name is {len(name)} chars, spec max is 64")
          if name != root.name:
              errors.append(f"name {name!r} must match the directory name {root.name!r}")
      if not desc:
          errors.append("required field `description` is missing")
      elif not 1 <= len(desc) <= 1024:
          errors.append(f"description is {len(desc)} chars, spec allows 1-1024")
  
      # 4. field surface
      for k in fm:
          if k in HOST_ONLY_FIELDS:
              errors.append(
                  f"frontmatter field {k!r} is a host-specific extension, not in the spec; "
                  "portable skills use only the six spec fields")
          elif k not in SPEC_FIELDS:
              warnings.append(f"unknown frontmatter field {k!r}")
      if fm.get("license") in ("", None):
          notes.append("no `license` field - see LICENSE for the actual grant")
      if not fm.get("compatibility"):
          notes.append("no `compatibility` field: hosts cannot warn about network/python needs")
  
      # 5. dependency-free scripts
      py = sorted((root / "scripts").glob("*.py")) if (root / "scripts").is_dir() else []
      if not py:
          notes.append("no scripts/ directory: instruction-only skill")
      for script, mods in stdlib_imports(py).items():
          errors.append(f"{script} imports non-stdlib modules {mods}; the skill must run "
                        "without a pip install")
  
      # 6. size budget
      lines = text.count("\n") + 1
      if lines > 500:
          warnings.append(f"SKILL.md is {lines} lines; spec recommends under 500")
  
      if not a.quiet:
          n_scripts = len(py)
          print(f"skill      : {name or '(unnamed)'}")
          print(f"version    : {fm.get('metadata', {}).get('version', '(unset)')}")
          print(f"dir        : {root}")
          print(f"description: {len(desc)} chars   body: {lines} lines   scripts: {n_scripts} (stdlib only)")
          print(f"fields     : {', '.join(fm) or '(none)'}")
          print()
          print("Host compatibility")
          print("-" * 78)
          for host, where, how in HOSTS:
              verdict = "OK " if not any("host-specific extension" in e for e in errors) else "HOST-ONLY"
              print(f"{verdict:10} {host:13} {how:26} {where}")
          print()
          for n in notes:
              print(f"note   {n}")
          for w in warnings:
              print(f"warn   {w}")
          for e in errors:
              print(f"FAIL   {e}")
          print()
          print("VALID" if not errors else f"INVALID ({len(errors)} error(s))")
      return 0 if not errors else 1
  
  
  if __name__ == "__main__":
      sys.exit(main())
  

• references/thresholds.md6.0 KBc237802ed8a8upstream ↗

  # Thresholds, schemes, and a worked example
  
  Every number here was fetched live from this box; each block carries its own access date.
  Do not carry a number forward from this file into a new curation — re-fetch it.
  
  ## Frequency thresholds (gnomAD)
  
  | Code | Direction | Strength | Threshold | Note |
  |---|---|---|---|---|
  | BA1 | Benign | Standalone | grpmax AF > 0.05 | Fires alone; curation stops |
  | BS1 | Benign | Strong | AF > gene-specific prevalence | Default 0.01; TP53 panel threshold overrides |
  | PM2 | Pathogenic | Supporting | Absent (AF = 0) | "Ultra-rare" needs a declared cutoff |
  
  Use `faf95.popmax` for the grpmax comparison; use `ac/an` for the plain AF. gnomAD v4's
  `joint` block exposes `ac`, `an`, `homozygote_count`, `faf95{popmax, popmax_population}` —
  **not** `af`.
  
  **Same release, three blocks — the number depends on which one.** D393V in `gnomad_r4`:
  `joint` 2/1,613,494 = 1.24e-06 · `exome` 2/1,461,238 = 1.37e-06 · carriers identical. An AF
  that matches a different block exactly is not "a superseded release", it is the wrong
  subset. Quote `joint` and name the block; a subset error is invisible until someone
  re-fetches, and it looks like version drift when they do.
  
  ## Predictor thresholds
  
  | Tool | Damaging | Tolerated | Indeterminate | VEP field |
  |---|---|---|---|---|
  | REVEL | ≥ 0.5 | < 0.5 | — | `revel` (with `REVEL=1`) |
  | CADD PHRED | ≥ 20 | < 20 | — | `cadd_phred` (with `CADD=1`) |
  | BayesDel addAF | ≥ 0.13 | ≤ −0.36 | −0.36 … 0.13 | `bayesdel_addaf_score` (dbNSFP param) |
  | PolyPhen | D / P | B | — | `polyphen_prediction` (with `PolyPhen=1`) |
  | SIFT | deleterious | tolerated | — | `sift_prediction` (default, no flag) |
  | AlphaMissense | likely_pathogenic | likely_benign | ambiguous | `alphamissense.am_class` |
  
  BayesDel and AlphaMissense are not ACMG codes on their own — they are extra voices. The
  count3 scheme counts BayesDel; AlphaMissense is reported but not counted by either scheme.
  
  ## The two PP3/BP4 schemes
  
  | | count3 (project default) | revel (ClinGen SVI 2022 style) |
  |---|---|---|
  | Condition | ≥3 available tools damaging, none tolerated | REVEL ≥ 0.75 |
  | Denominator | only tools with data; NA excluded | REVEL only |
  | Minimum data | ≥2 tools | REVEL present |
  | Effect on D393V | PP3 fires (4/4 damaging) | PP3 does not fire (0.69 is mid-band) |
  
  That divergence is the point of declaring a scheme. Same variant, same scores, two honest
  answers. A paper that switches schemes between variants is not reproducible.
  
  Downstream of that: **any write-up quoting a PP3/BP4 call names the scheme beside the code**,
  not only in Methods. A bare "4 of 4 damaging" reads as a fact about the variant when it is a
  fact about the scheme; on D393V the same payload derives nothing at all under `revel`.
  
  ## Worked example — TP53 p.Asp393Val, NM_000546.6:c.1178A>T
  
  Fetched values (gnomad_r4, accessed 2026-10-24):
  
  | Source | Value |
  |---|---|
  | ClinVar 3634421 | VCV003634421 · Uncertain significance · criteria provided, single submitter |
  | gnomAD joint | ac 2 / an 1,613,494 → AF 1.24e-06 · homozygotes 0 · faf95 popmax 3.65e-06 (sas) |
  | VEP | REVEL 0.69 · CADD 23.4 · BayesDel addAF 0.312 · SIFT deleterious · AlphaMissense likely_pathogenic |
  | VEP | PolyPhen absent (no coverage at this position — verified, not a bug) |
  | PubMed | 0 hits for the exact protein change under `[tiab]` phrase search |
  
  Derived:
  
  - **PP3 fires** under count3 (4 of 4 available damaging). Does not fire under revel.
  - **PM2 does not fire.** Two alleles exist, so the variant is not absent. Under AF = 0 the
    condition is false. It would need a declared ultra-rare cutoff (e.g. AF < 1e-5) to apply —
    and then the Methods must say so.
  - PS3/BS3: no published assay for this change → no functional code.
  
  Lean: pathogenic-leaning, untested. Two versions of the write-up are defensible; what is
  not defensible is claiming PM2 while reporting two carriers in the same table.
  
  ## Worked example — TP53 p.Gly245Ser, NM_000546.6:c.733G>A
  
  The same variant run twice: once with every identifier correct, once with a plausible but
  wrong gnomAD id. The scores are identical because VEP was never the problem; only the
  frequency source changed, and that is enough to invent a pathogenic code.
  
  Correct identifiers (gnomad_r4, accessed 2026-10-26): ClinVar 12365, gnomAD
  `17-7674230-C-T`, ClinVar SPDI `NC_000017.11:7674229:C:T` (0-based → 1-based 7674230).
  
  | Source | Value |
  |---|---|
  | ClinVar 12365 | `c.733G>A (p.Gly245Ser)` · Pathogenic · reviewed by expert panel |
  | gnomAD r4 joint | ac 5 · AF 3.10e-06 · homozygotes 0 · faf95 popmax 6.8e-07 (nfe) |
  | VEP | REVEL 0.947 · CADD 27.5 · SIFT deleterious · AlphaMissense likely_pathogenic · PolyPhen absent |
  | PubMed | 54 hits for `p.Gly245Ser`/`G245S` — titles only, none read |
  
  Derived: **PP3 fires under both schemes** (REVEL 0.947 ≥ 0.75; 3–4 of available tools
  damaging). **PM2 does not fire** — five carriers exist, so "absent" is false, and no
  ultra-rare cutoff was declared. PS3/BS3: no assay read, so none assigned. Lean:
  pathogenic-leaning, untested. ClinVar's expert-panel call is context, not a code.
  
  ### The trap: same variant, wrong gnomAD id
  
  Substituting `17-7676154-C-T` — a real TP53 position, 1,924 bp away — produced
  `absent: true` and a PM2. Nothing else in the pipeline objected. gnomAD answers
  `Variant not found` identically for a wrong position and for a genuine absence, so the
  payload is internally consistent and reads as "ultra-rare pathogenic variant" when the
  variant is not even the one being curated. The consistency guard now blocks this by
  comparing the id against VEP's resolved locus; PM2 additionally requires `ok: true`.
  
  Lesson worth carrying past this gene: frequency databases treat "not found" as a normal
  answer, so an unverified identifier converts a lookup failure into evidence.
  
  ## Sources for the thresholds
  
  - ACMG/AMP 2015 standards (Richards et al.) — code definitions and strengths
  - ClinGen SVI recommendations — PP3/BP4 calibration, PM2 downgrade, NA handling
  - gnomAD v4 GraphQL, Ensembl VEP REST (REST endpoints verified live, this box)
  

• .github/workflows/validate.yml752 B052c0060ba53upstream ↗

  name: validate
  
  on:
    push:
      branches: [main]
    pull_request:
    workflow_dispatch:      # live API check is manual: Ensembl VEP flaps ~50% of the time
  
  jobs:
    offline:
      name: validate + tests (py${{ matrix.python }})
      runs-on: ubuntu-latest
      strategy:
        fail-fast: false
        matrix:
          python: ["3.10", "3.12"]
      steps:
        - uses: actions/checkout@v4
        - uses: actions/setup-python@v5
          with:
            python-version: ${{ matrix.python }}
  
        - name: Skill conforms to the Agent Skills spec
          run: python3 scripts/validate_skill.py
  
        - name: Scripts compile
          run: python3 -m py_compile scripts/*.py
  
        - name: Engine tests (offline, no network)
          run: python3 tests/test_engine.py
  

• .gitignore60 Bd93a8d65ededupstream ↗

  __pycache__/
  *.pyc
  .venv/
  *.json
  !examples/*.json
  .DS_Store
  

• CHANGELOG.md5.6 KB0d9a664c66edupstream ↗

  # Changelog
  
  Version lives in `metadata.version` in `SKILL.md`. Thresholds, schemes and guard
  behaviour are part of the contract: changing them changes the codes a given variant
  receives, so pin the tool version in Methods alongside the dataset and access date.
  
  ## [2.1.0] - 2026-10-26
  
  **Guidance only.** No change to thresholds, schemes, the guard, or the codes a variant receives —
  a 2.0.0 and 2.1.0 run on the same payload agree. Safe to upgrade.
  
  ### Corrected
  - **"gnomAD version drift" was a misdiagnosis.** A pitfall claimed the same variant reads
    1.37e-06 in one release and 1.24e-06 in `gnomad_r4`. Both numbers are `gnomad_r4`, same two
    carriers, different blocks: `joint` 2/1,613,494 = 1.2395e-06 and `exome` 2/1,461,238 =
    1.3687e-06 (re-verified against the live API 2026-10-26). Curate from `joint`, and record the
    block beside the dataset and date — an unlabelled AF is unreproducible even when correct.
  
  ### Changed
  - **Name the scheme beside the code, not only in Methods.** D393V fires PP3 under `count3` and
    derives nothing at all under `revel`; a bare "PP3, 4 of 4" is refutable by anyone who runs the
    other scheme. Quote the code as "PP3 under count3, 4 of 4 available".
  
  ## [2.0.0] - 2026-10-26
  
  **Breaking: the codes a variant receives can differ from v1.0.0.** The lean is now ClinGen SVI
  points (Tavtigian 2020) instead of a code count, and the identifier guard can no longer pass
  without comparing anything. Re-run any v1.0.0 curation you intend to quote.
  
  ### Changed
  - **Lean: points, not votes.** A single Supporting code used to read "pathogenic-leaning"
    because one pathogenic code outnumbered zero benign ones. Any variant that was merely absent
    from gnomAD and had no predictors therefore looked pathogenic-leaning — synonymous, intronic
    and UTR variants included. Reported case: CD68 `NC_000017.11:g.7579817G>C`, ClinVar Benign,
    absent from gnomAD, CADD 6.68 — produced "Pathogenic-leaning". Now: 1 point, VUS, with the
    ClinVar conflict printed.
  - **Direction is no longer called a classification.** The output states the point totals, the
    direction, and that direction is not a classification. `Verdict:` is printed as a stable line
    (`PATHOGENIC` / `BENIGN` / `VUS`).
  - **The 6-point ceiling is now documented and tested.** The seven public-data codes cannot reach
    the 10-point pathogenic threshold, by construction. A test pins this so it is not "fixed" into
    over-calling.
  
  ### Fixed
  - **The guard could be bypassed.** `check_identity` parsed only `c.` HGVS, so a ClinVar record
    filed as genomic (`NC_000017.11:g.7579817G>C`) left `clinvar_hgvs: null` and returned
    `ok: true` having compared nothing. It now parses `c.`/`g.`/`n.`, cross-checks ClinVar
    against VEP's resolved locus, and returns **not ok** when no cross-check was possible.
    `consistency.checks` lists the comparisons actually performed.
  - **Accession prefix caused false mismatches.** Comparing raw HGVS strings marked
    `NM_000546.6(TP53):c.733G>A` against `c.733G>A` as a wrong variant. Comparison now uses
    parsed coordinates.
  - **A wrong MANE transcript silently emptied VEP.** No consequence returned for the requested
    accession produced blank scores and read as "no predictors". Now flagged as
    `vep.transcript_mismatch`.
  - **CNVs and structural variants are rejected, not half-curated.** ClinVar files them under
    `obj_type` (`variant_type` is null, so a check against that never fires). The run stops with
    exit code 2, the variant type and size, and a pointer to Riggs 2020.
  
  ### Added
  - Genomic (`g.`) and non-coding (`n.`) HGVS support.
  - Offline guard tests: genomic cross-check, wrong locus, nothing-to-compare, accession prefix,
    and out-of-scope rejection with a canned ClinVar response.
  
  ## [1.0.0] - 2026-10-26
  
  First public release. Seven codes (BA1, BS1, PM2, PP3, BP4, PS3, BS3) derived from
  ClinVar, gnomAD v4, Ensembl VEP and PubMed; two declared PP3 schemes (`count3`,
  `revel`); identifier consistency guard; offline test suite; spec-conformance validator.
  
  ### Fixed during v1.0.0 testing
  
  These were found by running the tool against real variants, not by reading it.
  
  - **A typo'd gnomAD id silently became PM2.** gnomAD answers `Variant not found` for a
    wrong position *and* for a genuinely unobserved allele, so a mistyped id read as
    "absent" and manufactured pathogenic evidence. `acmg_codes.py` now refuses PM2 when
    `consistency.ok` is false. Demonstrated on TP53 p.Gly245Ser: the real variant has
    AC=5 (AF 3.1e-06) and does not deserve PM2.
  - **No cross-check between identifiers.** ClinVar could report a classification for
    variant A while gnomAD and VEP answered for variant B. Added `check_identity()`,
    which compares the ClinVar title's HGVS against the requested HGVS and VEP's locus
    and alleles against the gnomAD id.
  - **`dbNSFP` rode on the AlphaMissense request.** The same URL both returned HTTP 500
    and suppressed AlphaMissense, silently thinning the predictor set. Split into two
    calls and merged; BayesDel now appears correctly (4/4 tools on G245S).
  - **Retries too thin for a flapping endpoint.** Ensembl VEP returns 200/500/503 within
    seconds on identical requests. 5xx is now retried 6x with jitter; 4xx fails fast
    instead of burning retries.
  - **Minus-strand allele orientation.** VEP returns `allele_string` in transcript
    orientation, so on TP53 every variant compared as a mismatch until complemented
    against `strand == -1`.
  
  ### Added
  
  - `scripts/validate_skill.py` — spec conformance and host-compatibility check.
  - `tests/test_engine.py` — offline assertions that PM2 refuses unverified ids.
  - `examples/g245s.json` — real fetched payload used as a test fixture.
  - `compatibility` and `metadata.version` frontmatter fields.
  

• CITATION.cff1.0 KB96deaf6b96daupstream ↗

  cff-version: 1.2.0
  message: "If you use this tool, cite it as below."
  title: "acmg-variant-curation"
  abstract: >-
    An Agent Skills package that derives ACMG/AMP 2015 codes for germline SNVs/indels
    from ClinVar, gnomAD, Ensembl VEP and PubMed, printing the evidence basis for every
    code and refusing pathogenic codes on unverified variant identifiers.
  type: software
  version: 2.1.0
  date-released: 2026-10-26
  license: MIT
  repository-code: "https://github.com/darshak26/acmg-variant-curation"
  keywords:
    - acmg
    - variant-interpretation
    - vus
    - variant-of-uncertain-significance
    - clinvar
    - gnomad
    - bioinformatics
    - agent-skills
  authors:
    - family-names: Padmane
      given-names: Swapnil
      email: [email protected]
  preferred-citation:
    type: software
    title: "acmg-variant-curation (v1.0.0): ACMG code derivation from public variant databases as an agent skill"
    authors:
      - family-names: Padmane
        given-names: Swapnil
    year: 2026
    version: 2.1.0
    url: "https://github.com/darshak26/acmg-variant-curation"
  

• examples/g245s.json2.7 KBc210a93fd968upstream ↗

  {
    "fetched_at": "2026-10-26T04:24:56+00:00",
    "clinvar": {
      "variation_id": "12365",
      "obj_type": "single nucleotide variant",
      "accession": "VCV000012365",
      "title": "NM_000546.6(TP53):c.733G>A (p.Gly245Ser)",
      "gene_sort": "TP53",
      "classification": "Pathogenic",
      "review_status": "reviewed by expert panel",
      "last_evaluated": "2024/08/05 00:00",
      "source": "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=clinvar&retmode=json&id=12365"
    },
    "gnomad": {
      "variant_id": "17-7674230-C-T",
      "rsids": [
        "rs28934575"
      ],
      "dataset": "gnomad_r4",
      "source": "https://gnomad.broadinstitute.org/api",
      "absent": false,
      "joint": {
        "ac": 5,
        "an": 1613726,
        "af": 3.098419434278186e-06,
        "homozygotes": 0,
        "faf95_popmax": 6.800000000000001e-07,
        "popmax_population": "nfe"
      },
      "exome": {
        "ac": 4,
        "an": 1461784,
        "af": 2.7363823930211302e-06,
        "homozygotes": 0,
        "faf95_popmax": 7.200000000000001e-07,
        "popmax_population": "nfe"
      },
      "genome": {
        "ac": 1,
        "an": 151942,
        "af": 6.581458714509484e-06,
        "homozygotes": 0,
        "faf95_popmax": null,
        "popmax_population": null
      },
      "af": 3.098419434278186e-06,
      "homozygotes": 0,
      "faf95_popmax": 6.800000000000001e-07,
      "popmax_population": "nfe"
    },
    "vep": {
      "hgvs": "NM_000546.6:c.733G>A",
      "transcript_used": "ENST00000269305",
      "most_severe_consequence": "missense_variant",
      "seq_region_name": "17",
      "start": 7674230,
      "allele_string": "G/A",
      "strand": -1,
      "hgvsp": null,
      "hgvsc": null,
      "scores": {
        "revel": 0.947,
        "cadd_phred": 27.5,
        "bayesdel_addaf": 0.550935,
        "polyphen": null,
        "sift": "deleterious",
        "alphamissense": "likely_pathogenic"
      },
      "bayesdel_merge": "ok",
      "transcript_mismatch": null,
      "tools_with_data": [
        "alphamissense",
        "bayesdel_addaf",
        "cadd_phred",
        "revel",
        "sift"
      ],
      "tools_missing": [
        "polyphen"
      ],
      "source": "https://rest.ensembl.org/vep/human/hgvs/NM_000546.6%3Ac.733G%3EA?transcript_id=ENST00000269305&CADD=1&REVEL=1&PolyPhen=1&SIFT=1&AlphaMissense=1&content-type=application/json"
    },
    "pubmed": {
      "skipped": "identifier not supplied"
    },
    "consistency": {
      "ok": true,
      "problems": [],
      "checks": [
        "requested HGVS == ClinVar HGVS (c.733G>A)",
        "gnomAD id == VEP locus (17:7674230 C>T)"
      ],
      "identifiers": {
        "requested_hgvs": "NM_000546.6:c.733G>A",
        "clinvar_title": "NM_000546.6(TP53):c.733G>A (p.Gly245Ser)",
        "clinvar_hgvs": "c.733G>A",
        "vep_locus": "17:7674230 G/A",
        "gnomad_id": "17-7674230-C-T"
      }
    }
  }
  

• LICENSE1.0 KBcc791ba1e856upstream ↗

  MIT License
  
  Copyright (c) 2026 Swapnil Padmane
  
  Permission is hereby granted, free of charge, to any person obtaining a copy
  of this software and associated documentation files (the "Software"), to deal
  in the Software without restriction, including without limitation the rights
  to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  copies of the Software, and to permit persons to whom the Software is
  furnished to do so, subject to the following conditions:
  
  The above copyright notice and this permission notice shall be included in all
  copies or substantial portions of the Software.
  
  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  SOFTWARE.
  

• NOTICE539 B5776c15774eaupstream ↗

  # NOTICE
  
  This software implements parts of the ACMG/AMP 2015 variant interpretation framework
  (Richards et al., *Genet Med* 17:405-424).
  
  It retrieves data from ClinVar, gnomAD, Ensembl VEP and PubMed. It is **not affiliated
  with or endorsed by** NCBI, the Broad Institute, Ensembl, or any of those projects.
  Data remains subject to each source's own terms of use.
  
  **Output is research evidence assembly, not a clinical classification.** It must not be
  used for patient care without review by a qualified clinical genetics professional.
  

• tests/test_engine.py6.8 KB876e4994ec88upstream ↗

  #!/usr/bin/env python3
  """Offline tests for the ACMG code engine and the identifier guard. Stdlib only, no network.
  
  Asserts the behaviours that are easy to regress and expensive to get wrong:
    1. PM2 does NOT fire when the variant is present in gnomAD (AC > 0).
    2. PM2 does NOT fire when the identifiers were not cross-checked.
    3. The lean uses ClinGen SVI points, not a code count.
    4. A genomic (g.) ClinVar record is cross-checked, not silently skipped.
    5. A guard that could not compare anything reports NOT ok.
    6. CNVs and structural variants are rejected, not half-curated.
  
  The ceiling is tested deliberately: the seven public-data codes max out at 6
  pathogenic points (PM2 1 + PP3 1 + PS3 4), below the 10-point threshold. This
  tool therefore CANNOT reach a pathogenic classification on its own - by design,
  and the tests pin that so nobody "fixes" it into over-calling.
  
  Run: python3 tests/test_engine.py
  """
  import importlib.util
  import json
  import pathlib
  import subprocess
  import sys
  import tempfile
  
  ROOT = pathlib.Path(__file__).resolve().parent.parent
  ENGINE = ROOT / "scripts" / "acmg_codes.py"
  FETCH = ROOT / "scripts" / "fetch_evidence.py"
  FIXTURE = ROOT / "examples" / "g245s.json"
  
  failures = []
  
  
  def _load_fetch():
      spec = importlib.util.spec_from_file_location("fetch_evidence", FETCH)
      mod = importlib.util.module_from_spec(spec)
      spec.loader.exec_module(mod)          # module level is definitions only, no network
      return mod
  
  
  def run(payload, *extra) -> str:
      with tempfile.NamedTemporaryFile("w", suffix=".json", delete=False) as f:
          json.dump(payload, f)
          tmp = f.name
      r = subprocess.run([sys.executable, str(ENGINE), tmp, *extra],
                         capture_output=True, text=True)
      if r.returncode != 0:
          failures.append(f"engine exited {r.returncode}: {r.stderr[:300]}")
      return r.stdout
  
  
  def check(label, cond, detail=""):
      if cond:
          print(f"  ok   {label}")
      else:
          failures.append(f"{label} {detail}")
          print(f"  FAIL {label} {detail}")
  
  
  def payload(**over):
      """Synthetic payload: absent from gnomAD, no predictors, identifiers verified."""
      d = {"fetched_at": "synthetic",
           "clinvar": {"title": "NM_000546.6(TP53):c.733G>A (p.Gly245Ser)",
                       "classification": "Pathogenic", "review_status": "x"},
           "gnomad": {"variant_id": "17-7674230-C-T", "dataset": "gnomad_r4",
                      "absent": True, "af": None},
           "vep": {"scores": {}},
           "consistency": {"ok": True, "problems": [], "checks": ["synthetic"],
                           "identifiers": {}}}
      d.update(over)
      return d
  
  
  print("case 1: real payload - variant present in gnomAD")
  base = json.loads(FIXTURE.read_text())
  out = run(base)
  check("PP3 applied", "PP3" in out)
  check("PM2 withheld (AC>0)", "PM2 NOT applied" in out,
        f"gnomad ac={base['gnomad']['joint']['ac']}")
  
  print("case 2: gnomAD absent, identifiers disagree -> PM2 suppressed")
  bad = json.loads(json.dumps(base))
  bad["gnomad"] = {"variant_id": "17-9999999-C-T", "dataset": "gnomad_r4",
                   "absent": True, "af": None}
  bad["consistency"] = {"ok": False, "problems": ["synthetic mismatch"], "identifiers": {}}
  out = run(bad)
  check("PM2 suppressed on unverified id", "PM2 NOT applied" in out)
  check("suppression names the mismatch", "identifiers disagree" in out)
  
  print("case 3: absent AND verified -> PM2 fires")
  out = run(payload())
  check("PM2 applied", "PM2" in out and "PM2 NOT applied" not in out)
  
  print("case 4: lean is SVI points, not a vote count")
  out = run(payload())
  check("one Supporting code is VUS, not pathogenic-leaning",
        "Verdict: VUS" in out and "Pathogenic-leaning" not in out)
  check("points line printed", "Points (ClinGen SVI): pathogenic 1 [PM2 1]" in out)
  check("1 point explained as insufficient", "below the 10-point threshold" in out)
  
  ceiling = payload(vep={"scores": {"revel": 0.95}})
  out = run(ceiling, "--functional", "damaging", "--scheme", "revel")
  check("PM2+PP3+PS3 = 6 points still VUS (public-data ceiling)",
        "Verdict: VUS" in out and "pathogenic 6" in out, out.split("Points")[-1][:90])
  
  freq = payload(gnomad={"variant_id": "x", "dataset": "gnomad_r4", "absent": False,
                         "af": 0.2, "faf95_popmax": 0.2})
  out = run(freq)
  check("high frequency -> BA1 -> BENIGN", "Verdict: BENIGN" in out and "BA1" in out)
  
  print("case 5: ClinVar conflict is flagged, not hidden")
  out = run(payload(clinvar={"title": "x", "classification": "Benign",
                             "review_status": "criteria provided"}))
  check("conflict note printed", "CONFLICT WITH CLINVAR" in out)
  check("verdict stays VUS", "Verdict: VUS" in out)
  
  print("case 6: identifier guard (offline, no network)")
  fe = _load_fetch()
  cv_g = {"title": "NC_000017.11:g.7579817G>C", "classification": "Benign"}
  vep_g = {"start": 7579817, "seq_region_name": "17", "allele_string": "G/C", "strand": 1}
  r = fe.check_identity(cv_g, {"variant_id": "17-7579817-G-C"}, vep_g,
                        "NC_000017.11:g.7579817G>C")
  check("genomic HGVS is actually cross-checked", r["ok"] and len(r["checks"]) >= 2,
        str(r["checks"]))
  r = fe.check_identity(cv_g, {"variant_id": "17-9999999-G-C"}, vep_g,
                        "NC_000017.11:g.7579817G>C")
  check("genomic record at the wrong locus fails", not r["ok"])
  r = fe.check_identity(None, {"variant_id": "17-7579817-G-C"}, None, None)
  check("nothing to compare -> not ok",
        (not r["ok"]) and "no cross-check" in r["problems"][0])
  r = fe.check_identity({"title": "NM_000546.6(TP53):c.733G>A (p.Gly245Ser)"},
                        {"variant_id": "17-7674230-C-T"},
                        # live VEP returns G/A (transcript orientation) for this coding
                        # HGVS on minus-strand TP53; complement -> C/T matches gnomAD
                        {"start": 7674230, "seq_region_name": "17",
                         "allele_string": "G/A", "strand": -1},
                        "NM_000546.6:c.733G>A")
  check("accession prefix does not cause a false mismatch", r["ok"], str(r["problems"]))
  r = fe.check_identity({"title": "NM_000546.6(TP53):c.710T>A (p.Met237Lys)"},
                        None, None, "NM_000546.6:c.733G>A")
  check("wrong coding variant still fails", not r["ok"])
  
  print("case 7: out-of-scope variant type is rejected (canned ClinVar response)")
  fe._get = lambda url, **kw: json.dumps({"result": {"150702": {
      "title": "GRCh38/hg38 17p13.1(chr17:7544902-7670581)x3",
      "obj_type": "copy number gain",
      "germline_classification": {"description": "Pathogenic"}}}})
  try:
      fe.fetch_clinvar("150702")
      check("CNV raises OutOfScope", False, "no exception raised")
  except fe.OutOfScope as e:
      check("CNV raises OutOfScope", "out of scope" in str(e).lower())
      check("rejection cites Riggs 2020", "Riggs" in str(e))
  
  print()
  if failures:
      print(f"FAILED ({len(failures)}):")
      for f in failures:
          print(" -", f)
      sys.exit(1)
  print("all engine tests passed")
  

Plain-text copy, no login: /skill-bundles/darshak26-acmg-variant-curation.md