Working in Python¶

python3 -m venv .venv
. .venv/bin/activate
pip install -q "ga4gh.vrs==2.4.0a6" "ga4gh.cat_vrs==0.8.0a6" "ga4gh.va_spec==0.5.0a7" "ga4gh.gkm==0.4.1"
pip install -q duckdb ijson    # for the Parquet and streaming examples

from importlib.metadata import version

from ga4gh.cat_vrs import CATVRS_VERSION
from ga4gh.core import CORE_VERSION
from ga4gh.va_spec import VASPEC_VERSION
from ga4gh.vrs import VRS_VERSION

for package in ("ga4gh.vrs", "ga4gh.cat_vrs", "ga4gh.va_spec", "ga4gh.gkm"):
    print(f"{package:14} {version(package)}")
print(f"VRS {VRS_VERSION}, Cat-VRS {CATVRS_VERSION}, VA-Spec {VASPEC_VERSION}, GKM-Core {CORE_VERSION}")

ga4gh.vrs      2.4.0a6
ga4gh.cat_vrs  0.8.0a6
ga4gh.va_spec  0.5.0a7
ga4gh.gkm      0.4.1
VRS 2.1.1, Cat-VRS 1.1.1, VA-Spec 1.1.0, GKM-Core 1.3.0

curl -sLO https://github.com/clingen-data-model/clinvar-gkm/raw/main/nch-igm-bundle.json.gz
curl -sO https://data.clinvar-gkm.org/v1/1.1/deltas/clinvar-gkm-delta_2026-1018.json.gz

import gzip
import json

with gzip.open("nch-igm-bundle.json.gz", "rt") as f:
    bundle = json.load(f)

print(bundle["schema_version"])
print({section: len(records) for section, records in bundle.items() if isinstance(records, dict)})


def resolve(pointer):
    """Return the record that a '#/section/id' pointer names."""
    section, key = pointer.removeprefix("#/").split("/", 1)
    return bundle[section][key]

1.1.0
{'sequenceReference': 24, 'location': 1913, 'allele': 2027, 'copyNumberChange': 1, 'gene': 460, 'variation': 2032, 'condition': 390, 'conditionSet': 33, 'submitter': 1, 'propositionPathogenicity': 943, 'propositionSomaticClinicalImpact': 2994, 'propositionClinvarCustom': 3, 'evidenceLine': 1497, 'scv': 2443}

scv = bundle["scv"]["clinvar.submission:SCV000864262.1"]
proposition = resolve(scv["proposition"])
variation = resolve(proposition["subject"])
allele = resolve(variation["constraints"][0]["allele"])
location = resolve(allele["location"])
sequence = resolve(location["sequenceReference"])

print(scv["classification"]["name"], "|", resolve(scv["contributions"][0]["contributor"])["name"])
print(proposition["type"], proposition["predicate"], "->", resolve(proposition["object"])["name"])
print(variation["id"], variation["name"])
print(allele["id"], allele["name"])
print(location["id"], location["start"], location["end"])
print(sequence["refgetAccession"], sequence["name"])

Benign | Institute for Genomic Medicine (IGM) Clinical Laboratory, Nationwide Children's Hospital
VariantPathogenicityProposition isCausalFor -> not specified
clinvar:599425 NM_018486.3(HDAC8):c.20C>T (p.Pro7Leu)
ga4gh:VA.yXEFCuYP7t_V3ALpYIMvWTy4PhfvT8ts NC_000023.11:g.72572742G>A
ga4gh:SL.sk5GccwF1vh_M4zGnjDYkoYNGrRKc6uW 72572741 72572742
SQ.w0WZEvgJF0zf_P4yyTzjjv9oW1z61HHP NC_000023.11

from ga4gh.core import ga4gh_identify
from ga4gh.vrs import models as vrs

vrs_allele = vrs.Allele(**allele)
print(type(vrs_allele.location).__name__, vrs_allele.location.root)


def content(record):
    """A record without its stored identifiers."""
    return {k: v for k, v in record.items() if k not in ("id", "digest")}


# Recompute the VRS identifier from content alone: inline the location and its sequence reference
fresh = vrs.Allele(**content(allele) | {"location": content(location) | {"sequenceReference": sequence}})
print(ga4gh_identify(fresh), ga4gh_identify(fresh) == allele["id"])

iriReference #/location/ga4gh:SL.sk5GccwF1vh_M4zGnjDYkoYNGrRKc6uW
ga4gh:VA.yXEFCuYP7t_V3ALpYIMvWTy4PhfvT8ts True

from ga4gh.cat_vrs import models as catvrs

categorical = catvrs.CategoricalVariant(**variation)
constraint = categorical.constraints[0].root
print(type(constraint).__name__, constraint.allele.root)
print([relation.primaryCoding.code.root for relation in constraint.relations])
print({e.name: e.value for e in categorical.extensions if isinstance(e.value, str)})

DefiningAlleleConstraint #/allele/ga4gh:VA.yXEFCuYP7t_V3ALpYIMvWTy4PhfvT8ts
['liftover_to', 'transcribed_to']
{'categoricalVariationType': 'CanonicalAllele', 'definingVrsVariationType': 'Allele', 'clinvarVariationType': 'single nucleotide variant', 'clinvarSubclassType': 'SimpleAllele', 'clinvarCytogeneticLocation': 'Xq13.1'}

from ga4gh.va_spec import base as va

pathogenicity = va.VariantPathogenicityProposition(**proposition)
statement = va.Statement(**scv)
print(pathogenicity.predicate, pathogenicity.subject.root, pathogenicity.object.root)
print(statement.direction, statement.strength.name, statement.classification.primaryCoding.code.root)
print(statement.proposition.root)
print([e.name for e in statement.extensions])

isCausalFor #/variation/clinvar:599425 #/condition/clinvar.trait:16789
disputes Definitive benign
#/propositionPathogenicity/SCV000864262-PATH
['clinvarScvId', 'clinvarScvVersion', 'submittedCondition', 'clinvarScvReviewStatus', 'submittedScvLocalKey', 'submissionLevel']

git clone -q --depth 1 --branch 1.1.0 --recurse-submodules --shallow-submodules \
    https://github.com/clingen-data-model/clinvar-gkm.git clinvar-gkm-1.1.0

import os
from collections import Counter
from pathlib import Path

from jsonschema import Draft202012Validator
from referencing import Registry, Resource
from ga4gh.gkm.bundles import BundleValidationError
from ga4gh.gkm.bundles.schema_validation import validate_bundle_schema

# Every schema in the 1.1.0 tree, keyed by $id: ClinVar-GKM's own and the GA4GH schemas they reference
schemas = {}
for directory, _, files in os.walk("clinvar-gkm-1.1.0/schema", followlinks=True):
    for name in files:
        try:
            document = json.loads(Path(directory, name).read_text())
        except (ValueError, UnicodeDecodeError):
            continue  # Makefiles and YAML sources
        if isinstance(document, dict) and "$id" in document:
            schemas[document["$id"]] = Resource.from_contents(document)
print(sorted(Counter("/".join(uri.split("/")[5:7]) for uri in schemas).items()))

bundle_schema = json.loads(Path("clinvar-gkm-1.1.0/schema/clinvar-gkm/clinvar-gkm-bundle.schema.json").read_text())
validator = Draft202012Validator(bundle_schema, registry=Registry().with_resources(schemas.items()))

validate_bundle_schema(bundle, bundle_schema, validator=validator)  # raises BundleValidationError on a violation
print("valid:", sum(len(records) for records in bundle.values() if isinstance(records, dict)), "records")

broken = json.loads(json.dumps({"schema_version": "1.1.0", "scv": {scv["id"]: scv}}))
broken["scv"][scv["id"]]["direction"] = "sideways"  # not a VA-Spec direction
try:
    validate_bundle_schema(broken, bundle_schema, validator=validator)
except BundleValidationError as error:
    print(error)

[('cat-vrs/1.1.1', 16), ('clinvar-gkm/1.1.0', 28), ('gkm-core/1.3.0', 11), ('va-spec/1.1.0', 40), ('vrs/2.1.1', 26)]
valid: 14761 records
Bundle does not match its JSON Schema at scv/clinvar.submission:SCV000864262.1/direction: 'sideways' is not one of ['supports', 'neutral', 'disputes']

import functools

import requests


@functools.cache
def fetch(uri):
    return Resource.from_contents(requests.get(uri, timeout=30).json())


fetching = Draft202012Validator(bundle_schema, registry=Registry(retrieve=fetch))
part = {"schema_version": bundle["schema_version"]} | {
    section: dict(list(records.items())[:25]) for section, records in bundle.items() if isinstance(records, dict)
}
validate_bundle_schema(part, bundle_schema, validator=fetching)
print("valid:", sum(len(records) for records in part.values() if isinstance(records, dict)), "records;",
      fetch.cache_info().currsize, "schemas fetched")

valid: 279 records; 74 schemas fetched

import duckdb

P = "https://data.clinvar-gkm.org/v1/1.1/parquet"

rows = duckdb.sql(f"SELECT id, location_id, data FROM '{P}/allele.parquet' LIMIT 3").fetchall()
for allele_id, location_id, data in rows:
    print(allele_id, vrs.Allele(**json.loads(data)).name)

ga4gh:VA.gaVM-oJ-UMO3cZMxe8-U0ktIKu7YvQK4 NC_000015.10:g.78513242G>A
ga4gh:VA.plnkuDHIFWLbLgn729ZU4dZlFhtM1Jco NC_000016.10:g.31224570C>A
ga4gh:VA.0Om9Ephu36fbTCIM_H4ushaqsJGXPOes NC_000017.11:g.41715616G>A

def by_id(table, ids, columns="id, data"):
    """Fetch the rows of one table whose id is in ids, keyed by id."""
    marks = ", ".join("?" * len(ids))
    query = f"SELECT {columns} FROM '{P}/{table}.parquet' WHERE id IN ({marks})"
    return {row[0]: row[1:] for row in duckdb.execute(query, list(ids)).fetchall()}


locations = by_id("location", [location_id for _, location_id, _ in rows], "id, sequence_reference_id, data")
sequences = by_id("sequenceReference", [sequence_id for sequence_id, _ in locations.values()])

for allele_id, location_id, data in rows:
    sequence_id, location_data = locations[location_id]
    inline = content(json.loads(location_data)) | {"sequenceReference": json.loads(sequences[sequence_id][0])}
    recomputed = ga4gh_identify(vrs.Allele(**content(json.loads(data)) | {"location": inline}))
    print(allele_id, recomputed == allele_id)

ga4gh:VA.gaVM-oJ-UMO3cZMxe8-U0ktIKu7YvQK4 True
ga4gh:VA.plnkuDHIFWLbLgn729ZU4dZlFhtM1Jco True
ga4gh:VA.0Om9Ephu36fbTCIM_H4ushaqsJGXPOes True

(prop_id, prop_data), = duckdb.sql(
    f"SELECT id, data FROM '{P}/propositionClinvarCustom.parquet' "
    "WHERE type = 'ClinvarRiskFactorProposition' AND gene_context_id IS NOT NULL LIMIT 1"
).fetchall()
(scv_data,), = duckdb.execute(
    f"SELECT data FROM '{P}/scv.parquet' WHERE proposition_id = ?", [prop_id]).fetchall()
risk = json.loads(prop_data)
variant_id = risk["subject"].removeprefix("#/variation/")
(variant_data,) = by_id("variation", [variant_id])[variant_id]

print(prop_id, risk["predicate"], risk["geneContextQualifier"]["name"])
print(va.Statement(**json.loads(scv_data)).classification.name)
print(catvrs.CategoricalVariant(**json.loads(variant_data)).name)

SCV000821694-RF isRiskFactorFor APC
risk factor
NM_000038.6(APC):c.3920T>A (p.Ile1307Lys)

from collections import Counter

from ga4gh.core.models import MappableConcept

MODELS = {
    "sequenceReference": vrs.SequenceReference, "location": vrs.SequenceLocation, "allele": vrs.Allele,
    "copyNumberCount": vrs.CopyNumberCount, "copyNumberChange": vrs.CopyNumberChange,
    "gene": MappableConcept, "condition": MappableConcept, "therapy": MappableConcept,
    "conditionSet": va.ConditionSet, "therapyGroup": va.TherapyGroup, "submitter": va.Agent,
    "variation": catvrs.CategoricalVariant, "evidenceLine": va.EvidenceLine,
    "scv": va.Statement, "vcv": va.Statement, "rcv": va.Statement,
}


def model_for(section, record):
    if section.startswith("proposition"):
        return getattr(va, record["type"])
    return MODELS[section]


tally = Counter()
for section, records in bundle.items():
    if isinstance(records, dict) and section != "propositionClinvarCustom":  # plain JSON, no library class
        for record in records.values():
            dumped = model_for(section, record)(**record).model_dump(mode="json", exclude_none=True)
            tally[section, "identical" if dumped == record else "differs"] += 1
for (section, result), n in tally.items():
    print(f"{section:33} {result:9} {n:5}")

sequenceReference                 identical    24
location                          identical  1913
allele                            identical  2027
copyNumberChange                  identical     1
gene                              identical   460
variation                         identical  2032
condition                         identical   390
conditionSet                      identical    33
submitter                         identical     1
propositionPathogenicity          identical   943
propositionSomaticClinicalImpact  identical  2994
evidenceLine                      identical  1497
scv                               differs    2443

import ijson

classifications = Counter()
with gzip.open("clinvar-gkm-delta_2026-1018.json.gz", "rb") as f:
    for scv_id, record in ijson.kvitems(f, "scv", use_float=True):
        classifications[va.Statement(**record).classification.name] += 1
print(sum(classifications.values()), classifications.most_common(4))

18138 [('Uncertain significance', 9724), ('Likely benign', 2766), ('Likely pathogenic', 2482), ('Pathogenic', 2337)]