PGx DB: pgx db build from CPIC (PGX-001..005, IN-053, CLI-035)
- scripts/pgx/fetch_sources.sh: CPIC API tables (CC0) into a dated folder with version and hashes. - simurg/pgx/db.py: CPIC HGVS locations and values -> left-normalised VCF variants (SNV, del, ins, dup, repeats, delins, chrM); IUPAC values are core-allele wildcards; alleles without definitions are flagged. - pipeline/pgx_db.py: deterministic SQLite (content-hash name), probe_map via the SNV VCF site mapper; GSA v4 PGx: 22 genes, 1354 alleles, 5002 probe links; two builds byte-identical.
This commit is contained in:
1 parent
c3e0681b6e
commit
0c5dd03ca8
8 files changed
+627
No files matched your search
@@ -810,6 +810,34 @@ pgx_cn_app = typer.Typer(help="Targeted PGx copy number (PCN-001).")
|
||||
pgx_app.add_typer(pgx_cn_app, name="copy-number")
|
||||
|
||||
|
||||
pgx_db_app = typer.Typer(help="PGx allele-definition database (PGX-001).")
|
||||
pgx_app.add_typer(pgx_db_app, name="db")
|
||||
|
||||
|
||||
@pgx_db_app.command("build")
|
||||
def pgx_db_build_cmd(
|
||||
sources: Annotated[Path, typer.Option("--sources", help="Folder from scripts/pgx/fetch_sources.sh (required).")],
|
||||
fasta: Annotated[Path, typer.Option("--fasta", help="GRCh38 FASTA for variant normalisation (required).")],
|
||||
out: Annotated[Path, typer.Option("--out", help="Output folder for pgx_db-<schema>-<sha12>.sqlite.")] = Path("."),
|
||||
manifests: Annotated[
|
||||
list[Path] | None, typer.Option("--manifests", help="BPM manifest(s) for the probe map (repeatable).")
|
||||
] = None,
|
||||
csv: Annotated[list[Path] | None, typer.Option("--csv", help="CSV manifest(s), same order (indels).")] = None,
|
||||
mask: Annotated[list[Path] | None, typer.Option("--mask", help="Probe mask TSV(s) (PGX-010).")] = None,
|
||||
force: ForceOpt = False,
|
||||
log_json: LogJson = False,
|
||||
debug: Debug = False,
|
||||
quiet: Quiet = False,
|
||||
) -> None:
|
||||
"""Build the PGx DB from public sources (CLI-035, PGX-002)."""
|
||||
from simurg.pipeline.pgx_db import PgxDbOptions, run_db_build
|
||||
|
||||
rep = _reporter("pgx db build", out, log_json, debug, quiet)
|
||||
o = PgxDbOptions(sources=sources, fasta=fasta, out=out, manifests=manifests or [], csvs=csv or [],
|
||||
masks=mask or [], force=force) # fmt: skip
|
||||
run_guarded(rep, lambda: run_db_build(o, rep))
|
||||
|
||||
|
||||
@pgx_cn_app.command("train")
|
||||
def pgx_cn_train_cmd(
|
||||
gtcs: Annotated[Path, typer.Option("--gtcs", help="Folder searched recursively for training GTCs (required).")],
|
||||
|
||||
@@ -171,6 +171,7 @@ CATALOGUE: dict[str, MessageDef] = dict(
|
||||
_m("W_PGX_NO_ALT_CONTIG", "warning", "FASTA lacks a contig carrying a PGx region; REF written as N", "IN-062"),
|
||||
_m("W_PGX_BATCH_COMPOSITION", "warning", "Batch spans several lots/batches or > 14 days of scans", "PCN-006"),
|
||||
_m("E_PGX_CN_FAILED", "error", "PGx CN calling failed for the sample", "PCN-002"),
|
||||
_m("E_PGX_SOURCE", "error", "PGx source tables missing; run scripts/pgx/fetch_sources.sh", "PGX-001"),
|
||||
# --- methylation control QC (MET-001 .. MET-018) -------------------------------------------------
|
||||
_m("E_METH_NO_ASSAY", "error", "Probe section not found", "MET-002"),
|
||||
_m("E_METH_NO_CONTROLS", "error", "Missing control probes in manifest", "MET-002"),
|
||||
|
||||
@@ -0,0 +1,223 @@
|
||||
"""PGx allele definitions from public sources → normalised VCF-style variants — spec PGX-001 … PGX-005.
|
||||
Pure functions: the reference sequence comes in as a ``fetch(contig, start0, end0)`` callable.
|
||||
|
||||
CPIC tables (CC0, https://api.cpicpgx.org) [DB: CPIC allele definition, allele function, diplotype→phenotype]:
|
||||
``sequence_location`` gives each defining position as HGVS on the chromosome (``g.``/``m.``); ``allele_location_value``
|
||||
gives each allele's value there. Values are literal sequences (``A``, ``AA`` for a dup, ``TA(7)`` for a repeat) or
|
||||
``del…``/``ins…`` forms; the reference allele definition carries the reference value. Every value becomes a
|
||||
replacement of a reference window, then a left-normalised, anchored VCF (POS, REF, ALT) — the representation the
|
||||
SNV VCFs use (IN-046), so definitions and array records compare directly.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import re
|
||||
from collections.abc import Callable
|
||||
from dataclasses import dataclass, field
|
||||
from typing import Any
|
||||
|
||||
Fetch = Callable[[str, int, int], str]
|
||||
HGVS_RE = re.compile(r"^([gm])\.(\d+)(?:_(\d+))?(.*)$")
|
||||
REPEAT_RE = re.compile(r"^([ACGT]+)\((\d+)\)$")
|
||||
IUPAC_WILDCARD = set("RYSWKMBDHVN") # CPIC core-allele definitions: the position does not constrain the allele
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class Window:
|
||||
"""Reference interval [start, end] (1-based, inclusive) replaced by allele values; ``insertion`` windows are the
|
||||
single anchor base S of ``g.S_EinsX`` (values are anchor + inserted sequence)."""
|
||||
|
||||
start: int
|
||||
end: int
|
||||
kind: str # snv | del | ins | dup | repeat | delins
|
||||
unit: str = "" # repeat unit
|
||||
|
||||
|
||||
def parse_location(hgvs: str) -> Window | None:
|
||||
"""First HGVS expression of a CPIC ``chromosomelocation`` (several are ``;``-joined for one site)."""
|
||||
first = hgvs.split(";")[0].strip()
|
||||
m = HGVS_RE.match(first)
|
||||
if not m:
|
||||
return None
|
||||
s = int(m.group(2))
|
||||
e = int(m.group(3)) if m.group(3) else s
|
||||
rest = m.group(4)
|
||||
if re.fullmatch(r"[ACGT]>[ACGT]", rest):
|
||||
return Window(s, s, "snv")
|
||||
if rest.startswith("delins"):
|
||||
return Window(s, e, "delins")
|
||||
if rest.startswith("del"):
|
||||
return Window(s, e, "del")
|
||||
if rest.startswith("ins"):
|
||||
return Window(s, s, "ins")
|
||||
if rest.startswith("dup"):
|
||||
return Window(s, e, "dup")
|
||||
rm = re.fullmatch(r"([ACGT]+)\[(\d+)\]", rest)
|
||||
if rm:
|
||||
return Window(s, s, "repeat", rm.group(1))
|
||||
return None
|
||||
|
||||
|
||||
def value_sequence(value: str, w: Window, fetch: Fetch, contig: str) -> tuple[int, int, str] | None:
|
||||
"""(start, end, sequence) that the value puts in place of the reference window; None = unusable/ambiguous."""
|
||||
v = value.strip()
|
||||
if " or " in v or ";" in v or not v:
|
||||
return None
|
||||
if w.kind == "repeat":
|
||||
m = REPEAT_RE.match(v)
|
||||
unit = w.unit
|
||||
n_ref = 0
|
||||
while fetch(contig, w.start - 1 + n_ref * len(unit), w.start - 1 + (n_ref + 1) * len(unit)) == unit:
|
||||
n_ref += 1
|
||||
end = w.start - 1 + n_ref * len(unit)
|
||||
if m:
|
||||
return w.start, end, m.group(1) * int(m.group(2))
|
||||
if re.fullmatch(r"[ACGT]+", v):
|
||||
return w.start, end, v
|
||||
return None
|
||||
if w.kind == "ins":
|
||||
anchor = fetch(contig, w.start - 1, w.start)
|
||||
if v == "del":
|
||||
return w.start, w.start, anchor
|
||||
if v.startswith("ins"):
|
||||
return w.start, w.start, anchor + v[3:]
|
||||
return None
|
||||
if w.kind == "delins":
|
||||
if v.startswith("del"):
|
||||
return w.start, w.end, fetch(contig, w.start - 1, w.end)
|
||||
if v.startswith("ins"):
|
||||
return w.start, w.end, v[3:]
|
||||
return None
|
||||
if v.startswith("del"):
|
||||
return w.start, w.end, ""
|
||||
if re.fullmatch(r"[ACGT]+", v):
|
||||
return w.start, w.end, v
|
||||
return None
|
||||
|
||||
|
||||
def normalize(fetch: Fetch, contig: str, start: int, end: int, alt: str) -> tuple[int, str, str]:
|
||||
"""Replacement of reference [start, end] (1-based; end = start − 1 for a pure insertion before start) by ``alt``
|
||||
→ left-normalised, anchored (POS, REF, ALT) [LIT: Tan et al. 2015 Bioinformatics 31:2202 (vt normalize)]."""
|
||||
ref = fetch(contig, start - 1, end)
|
||||
pos = start
|
||||
# trim the common suffix, then the common prefix
|
||||
while ref and alt and ref[-1] == alt[-1]:
|
||||
ref, alt = ref[:-1], alt[:-1]
|
||||
while len(ref) > 0 and len(alt) > 0 and ref[0] == alt[0]:
|
||||
ref, alt = ref[1:], alt[1:]
|
||||
pos += 1
|
||||
if ref and alt and len(ref) == len(alt):
|
||||
return pos, ref, alt # SNV / MNV
|
||||
# left shift the indel while the base before it equals its last base
|
||||
while pos > 1 and (ref or alt):
|
||||
prev = fetch(contig, pos - 2, pos - 1)
|
||||
body = ref or alt
|
||||
if prev != body[-1] or (ref and alt):
|
||||
break
|
||||
if ref:
|
||||
ref = prev + ref[:-1]
|
||||
else:
|
||||
alt = prev + alt[:-1]
|
||||
pos -= 1
|
||||
anchor = fetch(contig, pos - 2, pos - 1)
|
||||
return pos - 1, anchor + ref, anchor + alt
|
||||
|
||||
|
||||
@dataclass
|
||||
class VariantSite:
|
||||
location_id: int
|
||||
gene: str
|
||||
chrom: str # contig name as in the FASTA
|
||||
hgvs: str
|
||||
rsid: str | None
|
||||
ref_value: str | None
|
||||
alleles: dict[str, tuple[int, str, str]] = field(default_factory=dict) # value → normalised (pos, ref, alt)
|
||||
|
||||
|
||||
@dataclass
|
||||
class AlleleDef:
|
||||
gene: str
|
||||
name: str
|
||||
function: str | None
|
||||
activity: str | None
|
||||
reference: bool
|
||||
structural: bool
|
||||
variants: list[tuple[int, str]] = field(default_factory=list) # (location_id, value) where ≠ reference
|
||||
frequencies: dict[str, float] = field(default_factory=dict)
|
||||
wildcards: list[int] = field(default_factory=list) # location ids with an IUPAC value (either base)
|
||||
|
||||
@property
|
||||
def undefined(self) -> bool:
|
||||
"""A non-reference allele without any defining variant cannot be called (e.g. CPIC NUDT15 *2)."""
|
||||
return not self.reference and not self.variants
|
||||
|
||||
|
||||
def cpic_definitions(
|
||||
t: dict[str, list[dict[str, Any]]], fetch: Fetch, contig_of: Callable[[str], str | None]
|
||||
) -> tuple[dict[str, list[AlleleDef]], dict[int, VariantSite], list[str]]:
|
||||
"""CPIC tables → (gene → allele definitions, location id → site, problems)."""
|
||||
problems: list[str] = []
|
||||
gene_chr = {g["symbol"]: g.get("chr") for g in t["gene"]}
|
||||
locs = {loc["id"]: loc for loc in t["sequence_location"]}
|
||||
values: dict[int, list[tuple[int, str]]] = {}
|
||||
for v in t["allele_location_value"]:
|
||||
values.setdefault(v["alleledefinitionid"], []).append((v["locationid"], v["variantallele"]))
|
||||
func = {(a["genesymbol"], a["name"]): a for a in t["allele"]}
|
||||
defs: dict[str, list[AlleleDef]] = {}
|
||||
sites: dict[int, VariantSite] = {}
|
||||
ref_value: dict[int, str] = {}
|
||||
for d in sorted(t["allele_definition"], key=lambda d: (d["genesymbol"], d["name"])):
|
||||
if d.get("matchesreferencesequence"):
|
||||
for lid, val in values.get(d["id"], []):
|
||||
ref_value[lid] = val
|
||||
for d in sorted(t["allele_definition"], key=lambda d: (d["genesymbol"], d["name"])):
|
||||
g = d["genesymbol"]
|
||||
a = func.get((g, d["name"]), {})
|
||||
freq = {k: float(v) for k, v in (a.get("frequency") or {}).items() if isinstance(v, int | float)}
|
||||
ref_def, sv = bool(d.get("matchesreferencesequence")), bool(d.get("structuralvariation"))
|
||||
ad = AlleleDef(g, d["name"], a.get("clinicalfunctionalstatus"), a.get("activityvalue"), ref_def, sv, [], freq)
|
||||
for lid, val in sorted(values.get(d["id"], [])):
|
||||
loc = locs.get(lid)
|
||||
if loc is None:
|
||||
continue
|
||||
if lid not in sites:
|
||||
hg = loc.get("chromosomelocation") or ""
|
||||
chrom = "chrM" if hg.startswith("m.") else (gene_chr.get(g) or "")
|
||||
contig = contig_of(chrom) if chrom else None
|
||||
if contig is None:
|
||||
problems.append(f"{g} {loc['name']}: no contig for {chrom!r}")
|
||||
continue
|
||||
sites[lid] = VariantSite(lid, g, contig, hg, loc.get("dbsnpid"), ref_value.get(lid))
|
||||
site = sites[lid]
|
||||
if len(val) == 1 and val in IUPAC_WILDCARD:
|
||||
ad.wildcards.append(lid)
|
||||
elif val != site.ref_value:
|
||||
ad.variants.append((lid, val))
|
||||
defs.setdefault(g, []).append(ad)
|
||||
for lid, site in sites.items():
|
||||
w = parse_location(site.hgvs)
|
||||
if w is None:
|
||||
problems.append(f"{site.gene} {site.hgvs}: unparsed location")
|
||||
continue
|
||||
vals = {val for d in defs.get(site.gene, []) for (l2, val) in d.variants if l2 == lid}
|
||||
for val in sorted(vals):
|
||||
seq = value_sequence(val, w, fetch, site.chrom)
|
||||
if seq is None:
|
||||
problems.append(f"{site.gene} {site.hgvs}: value {val!r} not representable")
|
||||
continue
|
||||
s, e, alt = seq
|
||||
if w.kind == "ins":
|
||||
s, e, alt = s + 1, s, alt[1:]
|
||||
site.alleles[val] = normalize(fetch, site.chrom, s, e, alt)
|
||||
return defs, sites, problems
|
||||
|
||||
|
||||
def result_basis(gene: str, names: list[str]) -> str:
|
||||
"""PGX-002 ``genes.result_basis``."""
|
||||
if gene in ("HLA-A", "HLA-B"):
|
||||
return "tagging_variant"
|
||||
if gene in ("G6PD", "APOE"):
|
||||
return "haplotype_named"
|
||||
if names and all(n.startswith("*") for n in names):
|
||||
return "star_allele"
|
||||
return "variant_based"
|
||||
@@ -0,0 +1,223 @@
|
||||
"""``simurg pgx db build`` — the PGx allele-definition database (spec PGX-001 … PGX-005, PGX-010, IN-053, CLI-035).
|
||||
|
||||
I/O and orchestration; the parsing and normalisation live in ``simurg.pgx.db``. Sources are the dated folders written
|
||||
by ``scripts/pgx/fetch_sources.sh`` (``cpic-YYYYMMDD/`` with ``source.json``). The build is deterministic: rows are
|
||||
inserted in sorted order into a fresh SQLite file with fixed pragmas, and the file is named by its content hash.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import hashlib
|
||||
import json
|
||||
import sqlite3
|
||||
from dataclasses import dataclass, field
|
||||
from pathlib import Path
|
||||
from typing import Any
|
||||
|
||||
from simurg import __version__
|
||||
from simurg.core.errors import ExitCode, SimurgError
|
||||
from simurg.core.log import Reporter
|
||||
|
||||
DB_SCHEMA = 1
|
||||
CPIC_TABLES = ("gene", "allele", "allele_definition", "allele_location_value", "sequence_location", "gene_result",
|
||||
"gene_result_lookup") # fmt: skip
|
||||
GENE_PAD = 50_000 # probes within this distance of a gene's defining positions are resolved for probe_map
|
||||
|
||||
SCHEMA_SQL = """
|
||||
CREATE TABLE meta(key TEXT PRIMARY KEY, value TEXT NOT NULL);
|
||||
CREATE TABLE genes(gene TEXT PRIMARY KEY, result_basis TEXT NOT NULL, chrom TEXT);
|
||||
CREATE TABLE alleles(gene TEXT, allele TEXT, kind TEXT, function TEXT, activity TEXT, structural INTEGER,
|
||||
defined INTEGER, PRIMARY KEY(gene, allele));
|
||||
CREATE TABLE allele_variants(gene TEXT, allele TEXT, build TEXT, chrom TEXT, pos INTEGER, ref TEXT, alt TEXT,
|
||||
hgvs TEXT, rsid TEXT, location_id INTEGER, value TEXT);
|
||||
CREATE TABLE allele_wildcards(gene TEXT, allele TEXT, location_id INTEGER);
|
||||
CREATE TABLE sites(gene TEXT, location_id INTEGER, chrom TEXT, hgvs TEXT, rsid TEXT, ref_value TEXT,
|
||||
PRIMARY KEY(gene, location_id));
|
||||
CREATE TABLE site_values(gene TEXT, location_id INTEGER, value TEXT, pos INTEGER, ref TEXT, alt TEXT);
|
||||
CREATE TABLE probe_map(manifest TEXT, probe_name TEXT, chrom TEXT, pos INTEGER, ref TEXT, alt TEXT);
|
||||
CREATE TABLE masks(manifest TEXT, probe_name TEXT, reason TEXT);
|
||||
CREATE TABLE limitations(manifest TEXT, gene TEXT, allele TEXT, status TEXT, reported_as TEXT);
|
||||
CREATE TABLE frequencies(gene TEXT, allele TEXT, population TEXT, freq REAL);
|
||||
CREATE TABLE phenotypes(guideline TEXT, gene TEXT, lookup_key TEXT, activity_score TEXT, phenotype TEXT);
|
||||
CREATE INDEX av_pos ON allele_variants(chrom, pos);
|
||||
CREATE INDEX pm_pos ON probe_map(manifest, chrom, pos);
|
||||
"""
|
||||
|
||||
|
||||
@dataclass
|
||||
class PgxDbOptions:
|
||||
sources: Path
|
||||
fasta: Path
|
||||
out: Path
|
||||
manifests: list[Path] = field(default_factory=list)
|
||||
csvs: list[Path] = field(default_factory=list)
|
||||
masks: list[Path] = field(default_factory=list)
|
||||
force: bool = False
|
||||
|
||||
|
||||
def _latest(sources: Path, prefix: str) -> Path | None:
|
||||
cands = sorted(p for p in sources.glob(f"{prefix}-*") if (p / "source.json").exists())
|
||||
return cands[-1] if cands else None
|
||||
|
||||
|
||||
def _load_cpic(d: Path) -> dict[str, list[dict[str, Any]]]:
|
||||
t = {}
|
||||
for name in CPIC_TABLES:
|
||||
p = d / f"{name}.json"
|
||||
if not p.exists():
|
||||
raise SimurgError("E_PGX_SOURCE", f"{p} missing", hint="run scripts/pgx/fetch_sources.sh")
|
||||
t[name] = json.loads(p.read_text(encoding="utf-8"))
|
||||
return t
|
||||
|
||||
|
||||
def _probe_sites(
|
||||
bpm_path: Path,
|
||||
csv_path: Path | None,
|
||||
genome: Any,
|
||||
cmap: Any,
|
||||
windows: dict[str, list[tuple[int, int]]],
|
||||
rep: Reporter,
|
||||
) -> list[Any]:
|
||||
"""Plus-strand, normalised sites of the manifest probes inside the gene windows (the SNV VCF representation)."""
|
||||
from simurg.formats.bpm import read_bpm
|
||||
from simurg.harmonize.sites import ManifestView, map_locus
|
||||
|
||||
bpm = read_bpm(bpm_path)
|
||||
csvm = None
|
||||
if csv_path is not None:
|
||||
from simurg.formats.csv_manifest import read_csv_manifest
|
||||
|
||||
csvm = read_csv_manifest(csv_path)
|
||||
elif any("D" in s for s in bpm.snps):
|
||||
rep.warn("W_SKIP_INDEL_CONTEXT", f"{bpm_path.name}: no CSV manifest, indel probes are not mapped")
|
||||
src = csvm or bpm
|
||||
mv = ManifestView(bpm.names, bpm.snps, bpm.chroms, bpm.map_infos.tolist(), src.ilmn_strands, src.source_strands,
|
||||
src.ref_strands if any(src.ref_strands) else bpm.ref_strands,
|
||||
csvm.source_seqs if csvm else [""] * bpm.num_loci, bpm.assay_types.tolist(),
|
||||
bpm.intensity_only.tolist()) # fmt: skip
|
||||
out = []
|
||||
for i in range(bpm.num_loci):
|
||||
c = cmap.get(bpm.chroms[i])
|
||||
if c is None or c not in windows:
|
||||
continue
|
||||
p = int(bpm.map_infos[i])
|
||||
if not any(a - GENE_PAD <= p <= b + GENE_PAD for a, b in windows[c]):
|
||||
continue
|
||||
site, _ = map_locus(i, mv, genome.fetch, cmap.get, genome.build or "GRCh38")
|
||||
if site is not None:
|
||||
out.append(site)
|
||||
return out
|
||||
|
||||
|
||||
def run_db_build(o: PgxDbOptions, rep: Reporter) -> int:
|
||||
"""``simurg pgx db build`` (CLI-035)."""
|
||||
from simurg.formats.bpm import read_bpm
|
||||
from simurg.formats.fasta import ContigMap, Genome
|
||||
from simurg.pgx.db import cpic_definitions, result_basis
|
||||
|
||||
cdir = _latest(o.sources, "cpic")
|
||||
if cdir is None:
|
||||
raise SimurgError("E_PGX_SOURCE", f"no cpic-YYYYMMDD source folder in {o.sources}",
|
||||
hint="run scripts/pgx/fetch_sources.sh") # fmt: skip
|
||||
csrc = json.loads((cdir / "source.json").read_text())
|
||||
genome = Genome(o.fasta)
|
||||
if genome.build not in ("GRCh38", "unknown", ""):
|
||||
raise SimurgError("E_BUILD_MISMATCH", f"the PGx DB is GRCh38 only (FASTA is {genome.build})") # PCN-015
|
||||
cmap = ContigMap([e.name for e in genome.entries], "auto")
|
||||
t = _load_cpic(cdir)
|
||||
defs, sites, problems = cpic_definitions(t, genome.fetch, cmap.get)
|
||||
for p in problems:
|
||||
rep.info(f"definition note: {p}")
|
||||
rows: dict[str, list[tuple[Any, ...]]] = {k: [] for k in ("genes", "alleles", "allele_variants", "allele_wildcards",
|
||||
"sites", "site_values", "frequencies", "phenotypes",
|
||||
"probe_map", "masks")} # fmt: skip
|
||||
gene_chr = {g["symbol"]: g.get("chr") for g in t["gene"]}
|
||||
for gene, ds in sorted(defs.items()):
|
||||
rows["genes"].append((gene, result_basis(gene, [d.name for d in ds]), gene_chr.get(gene)))
|
||||
for d in ds:
|
||||
kind = "reference" if d.reference else ("star_allele" if d.name.startswith("*") else "named_haplotype")
|
||||
rows["alleles"].append(
|
||||
(gene, d.name, kind, d.function, d.activity, int(d.structural), int(not d.undefined))
|
||||
)
|
||||
for lid, val in d.variants:
|
||||
s = sites.get(lid)
|
||||
if s is None or val not in s.alleles:
|
||||
continue
|
||||
pos, ref, alt = s.alleles[val]
|
||||
rows["allele_variants"].append(
|
||||
(gene, d.name, "GRCh38", s.chrom, pos, ref, alt, s.hgvs, s.rsid, lid, val)
|
||||
)
|
||||
for lid in d.wildcards:
|
||||
rows["allele_wildcards"].append((gene, d.name, lid))
|
||||
for pop, f in sorted(d.frequencies.items()):
|
||||
rows["frequencies"].append((gene, d.name, pop, f))
|
||||
for lid, s in sorted(sites.items()):
|
||||
rows["sites"].append((s.gene, lid, s.chrom, s.hgvs, s.rsid, s.ref_value))
|
||||
for val, (pos, ref, alt) in sorted(s.alleles.items()):
|
||||
rows["site_values"].append((s.gene, lid, val, pos, ref, alt))
|
||||
results = {r["id"]: r for r in t["gene_result"]}
|
||||
for lk in t["gene_result_lookup"]:
|
||||
res = results.get(lk["phenotypeid"])
|
||||
if res is None:
|
||||
continue
|
||||
key = json.dumps(lk.get("lookupkey"), sort_keys=True, ensure_ascii=False)
|
||||
rows["phenotypes"].append(("cpic", res["genesymbol"], key, lk.get("totalactivityscore"), res["result"]))
|
||||
# probe map (PGX-003): manifest probes whose normalised site carries a DB variant
|
||||
windows: dict[str, list[tuple[int, int]]] = {}
|
||||
for av in rows["allele_variants"]:
|
||||
windows.setdefault(av[3], []).append((av[4], av[4] + len(av[5])))
|
||||
var_keys = {(av[3], av[4], av[5], av[6]) for av in rows["allele_variants"]}
|
||||
manifests = []
|
||||
csvs = dict(zip(o.manifests, o.csvs, strict=False))
|
||||
for bpm_path in o.manifests:
|
||||
name = read_bpm(bpm_path).manifest_name
|
||||
manifests.append(name)
|
||||
n = 0
|
||||
for ps in _probe_sites(bpm_path, csvs.get(bpm_path), genome, cmap, windows, rep):
|
||||
for alt in ps.alleles:
|
||||
if alt != ps.ref and (ps.contig, ps.pos, ps.ref, alt) in var_keys:
|
||||
rows["probe_map"].append((name, ps.name, ps.contig, ps.pos, ps.ref, alt))
|
||||
n += 1
|
||||
rep.info(f"{name}: {n} probe-to-variant link(s)")
|
||||
for mp in o.masks:
|
||||
for ln in mp.read_text(encoding="utf-8").splitlines():
|
||||
if ln.strip() and not ln.startswith("#"):
|
||||
cols = ln.split("\t")
|
||||
rows["masks"].append((cols[2] if len(cols) > 2 else mp.stem, cols[0], cols[1] if len(cols) > 1 else ""))
|
||||
meta = {"schema": "simurg.pgx_db", "schema_version": str(DB_SCHEMA), "builder": f"SimurgArray {__version__}",
|
||||
"genome_build": "GRCh38", "fasta_sha256": genome.sha256, "sqlite_version": sqlite3.sqlite_version,
|
||||
"allele_sources": json.dumps([{"name": "CPIC", "version": csrc["version"], "retrieved": csrc["retrieved"],
|
||||
"license": csrc.get("license")}]),
|
||||
"phenotype_sources": json.dumps([{"name": "CPIC", "version": csrc["version"]}]),
|
||||
"cpic_files_sha256": json.dumps(csrc.get("files", {}), sort_keys=True),
|
||||
"manifests": json.dumps(sorted(manifests)), "definition_notes": json.dumps(problems)} # fmt: skip
|
||||
# content hash over the sorted rows and meta (IN-053: the name and meta carry it; a file cannot hold its own hash)
|
||||
for rs in rows.values():
|
||||
rs.sort(key=lambda r: tuple(str(x) for x in r))
|
||||
sha = hashlib.sha256(json.dumps({"rows": rows, "meta": sorted(meta.items())}, sort_keys=True,
|
||||
default=str).encode()).hexdigest() # fmt: skip
|
||||
meta["content_sha256"] = sha
|
||||
o.out.mkdir(parents=True, exist_ok=True)
|
||||
tmp = o.out / ".pgx_db.part.sqlite"
|
||||
tmp.unlink(missing_ok=True)
|
||||
con = sqlite3.connect(tmp)
|
||||
con.execute("PRAGMA page_size=4096")
|
||||
con.execute("PRAGMA journal_mode=OFF")
|
||||
con.executescript(SCHEMA_SQL)
|
||||
for table, rs in rows.items():
|
||||
if rs:
|
||||
ph = ",".join("?" * len(rs[0]))
|
||||
con.executemany(f"INSERT INTO {table} VALUES ({ph})", rs)
|
||||
con.executemany("INSERT INTO meta VALUES (?, ?)", sorted(meta.items()))
|
||||
con.commit()
|
||||
con.execute("VACUUM")
|
||||
con.close()
|
||||
final = o.out / f"pgx_db-{DB_SCHEMA}-{sha[:12]}.sqlite"
|
||||
if final.exists() and not o.force:
|
||||
tmp.unlink()
|
||||
rep.info(f"identical PGx DB already present: {final}")
|
||||
return ExitCode.OK
|
||||
tmp.replace(final)
|
||||
rep.info(f"PGx DB: {len(defs)} genes, {len(rows['alleles'])} alleles, {len(rows['allele_variants'])} variant "
|
||||
f"rows, {len(rows['probe_map'])} probe links → {final} (content sha256 {sha})") # fmt: skip
|
||||
return ExitCode.OK
|
||||
Reference in new issue
Block a user