PGx CN: pgx copy-number train|call (CLI-031/032, PCN-001..017, IN-050), region BED v1, demo identities (K-028)
- simurg/pgx/cn.py: per-region CN mixture, mu_c = delta + gamma*log2((c+1)/3), free lower CN-0 component, Hardy-Weinberg-constrained weights over haplotype copy numbers, EM from a grid of label starts. - pipeline/pgx_cn.py: model file (.simurg-pgxcn), batch rules (PCN-004), per-sample CN VCF + tbi, batch and sample summaries. - The GSA v4 PGx demo arrays are 30 CEPH/1000G individuals (scripts/pgx/demo36_identity.tsv); against NYGC 30x SV RD_CN 224/240 region calls agree (GSTM1, UGT2B17 30/30), in-sample (docs/validation/V2_pgx_cn.md).
This commit is contained in:
kemal committed
1 parent
8c7414e894
commit
cc8c934167
14 files changed
+1306
-1
No files matched your search
@@ -803,6 +803,78 @@ def cn_annotate_cmd(
|
||||
run_guarded(rep, lambda: run_cn_annotate(opts, rep))
|
||||
|
||||
|
||||
# --------------------------------------------------------------------------------------------------------
|
||||
pgx_app = typer.Typer(help="Pharmacogenomics: targeted copy number, star alleles, phenotypes (v2).")
|
||||
app.add_typer(pgx_app, name="pgx")
|
||||
pgx_cn_app = typer.Typer(help="Targeted PGx copy number (PCN-001).")
|
||||
pgx_app.add_typer(pgx_cn_app, name="copy-number")
|
||||
|
||||
|
||||
@pgx_cn_app.command("train")
|
||||
def pgx_cn_train_cmd(
|
||||
gtcs: Annotated[Path, typer.Option("--gtcs", help="Folder searched recursively for training GTCs (required).")],
|
||||
bpm: Annotated[Path, typer.Option("--bpm", help="BPM manifest (required).")],
|
||||
fasta: Annotated[Path, typer.Option("--fasta", help="GRCh38 FASTA (required).")],
|
||||
out: Annotated[Path, typer.Option("--out", help="Output folder.")] = Path("."),
|
||||
mask: Annotated[Path | None, typer.Option("--mask", help="Probe mask TSV (IN-052).")] = None,
|
||||
regions: Annotated[Path | None, typer.Option("--regions", help="Region BED (default: built-in v1).")] = None,
|
||||
summary: Annotated[Path | None, typer.Option("--summary", help="sample_summary.tsv (assay control).")] = None,
|
||||
platform: Annotated[str, typer.Option("--platform", help="auto|LCG|EX")] = "auto",
|
||||
min_samples: Annotated[int, typer.Option("--min-samples", min=1)] = 96,
|
||||
min_passing: Annotated[int, typer.Option("--min-passing", min=1)] = 90,
|
||||
max_samples: Annotated[int, typer.Option("--max-samples", min=1)] = 3000,
|
||||
train_lrr_sd_max: Annotated[float, typer.Option("--train-lrr-sd-max")] = 0.20,
|
||||
train_call_rate_min: Annotated[float, typer.Option("--train-call-rate-min")] = 0.99,
|
||||
train_assay_control_min: Annotated[float, typer.Option("--train-assay-control-min")] = 1.0,
|
||||
force: ForceOpt = False,
|
||||
log_json: LogJson = False,
|
||||
debug: Debug = False,
|
||||
quiet: Quiet = False,
|
||||
) -> None:
|
||||
"""Train a PGx CN model (CLI-032, PCN-005)."""
|
||||
from simurg.pipeline.pgx_cn import PgxTrainOptions, run_train
|
||||
|
||||
if platform not in ("auto", "LCG", "EX"):
|
||||
raise typer.BadParameter("--platform auto|LCG|EX")
|
||||
rep = _reporter("pgx copy-number train", out, log_json, debug, quiet)
|
||||
o = PgxTrainOptions(gtcs=gtcs, bpm=bpm, fasta=fasta, out=out, mask=mask, regions=regions, summary=summary,
|
||||
platform=platform, min_samples=min_samples, min_passing=min_passing, max_samples=max_samples,
|
||||
lrr_sd_max=train_lrr_sd_max, call_rate_min=train_call_rate_min,
|
||||
assay_control_min=train_assay_control_min, force=force) # fmt: skip
|
||||
run_guarded(rep, lambda: run_train(o, rep))
|
||||
|
||||
|
||||
@pgx_cn_app.command("call")
|
||||
def pgx_cn_call_cmd(
|
||||
gtcs: Annotated[Path, typer.Option("--gtcs", help="Folder searched recursively for GTC files (required).")],
|
||||
model: Annotated[Path, typer.Option("--model", help="PGx CN model (.simurg-pgxcn, required).")],
|
||||
fasta: Annotated[Path, typer.Option("--fasta", help="GRCh38 FASTA with alt contigs (required).")],
|
||||
out: Annotated[Path, typer.Option("--out", help="Output folder (required).")],
|
||||
min_batch: Annotated[int, typer.Option("--min-batch", min=1)] = 24,
|
||||
min_passing: Annotated[int, typer.Option("--min-passing", min=1)] = 22,
|
||||
recommended_passing: Annotated[int, typer.Option("--recommended-passing", min=1)] = 96,
|
||||
lrr_sd_max: Annotated[float, typer.Option("--lrr-sd-max")] = 0.20,
|
||||
cluster_mismatch: Annotated[str, typer.Option("--cluster-mismatch", help="error|warn")] = "error",
|
||||
no_bgzip: Annotated[bool, typer.Option("--no-bgzip")] = False,
|
||||
force: ForceOpt = False,
|
||||
log_json: LogJson = False,
|
||||
debug: Debug = False,
|
||||
quiet: Quiet = False,
|
||||
) -> None:
|
||||
"""Call targeted PGx copy number per region (CLI-031, PCN-009)."""
|
||||
from simurg.pipeline.pgx_cn import PgxCallOptions, run_call
|
||||
|
||||
if cluster_mismatch not in ("error", "warn"):
|
||||
raise typer.BadParameter("--cluster-mismatch error|warn")
|
||||
if min_passing > min_batch:
|
||||
raise typer.BadParameter("--min-passing must be <= --min-batch")
|
||||
rep = _reporter("pgx copy-number call", out, log_json, debug, quiet)
|
||||
o = PgxCallOptions(gtcs=gtcs, model=model, fasta=fasta, out=out, min_batch=min_batch, min_passing=min_passing,
|
||||
recommended_passing=recommended_passing, lrr_sd_max=lrr_sd_max,
|
||||
cluster_mismatch=cluster_mismatch, bgzip=not no_bgzip, force=force) # fmt: skip
|
||||
run_guarded(rep, lambda: run_call(o, rep))
|
||||
|
||||
|
||||
# --------------------------------------------------------------------------------------------------------
|
||||
methyl_app = typer.Typer(help="Methylation arrays: control QC, detection, beta values (MET-001).")
|
||||
app.add_typer(methyl_app, name="methyl")
|
||||
|
||||
@@ -151,6 +151,26 @@ CATALOGUE: dict[str, MessageDef] = dict(
|
||||
_m("E_CN_SIDECAR_MISSING", "error", "cn_call_manifest.json not found in the cn call output", "CYTO-032"),
|
||||
_m("E_CN_VCF_UNINDEXED", "error", "CNV VCFs must be bgzipped and indexed for annotate", "CYTO-032"),
|
||||
_m("E_CN_FAILED", "error", "CN calling failed for the sample", "CYTO-040"),
|
||||
# --- targeted PGx copy number (PCN-001 .. PCN-017) ---------------------------------------------------
|
||||
_m("E_PGX_BATCH_TOO_SMALL", "error", "Too few (passing) samples in the batch for PGx CN calling", "PCN-004"),
|
||||
_m("W_PGX_BATCH_BELOW_RECOMMENDED", "warning", "Fewer passing samples than recommended", "PCN-004"),
|
||||
_m("E_PGX_TRAIN_TOO_SMALL", "error", "Too few (passing) samples to train a PGx CN model", "PCN-005"),
|
||||
_m(
|
||||
"E_MODEL_CLUSTER_MISMATCH",
|
||||
"error",
|
||||
"Input GTCs were called with another cluster file than the model",
|
||||
"IN-050",
|
||||
),
|
||||
_m(
|
||||
"W_MODEL_CLUSTER_MISMATCH",
|
||||
"warning",
|
||||
"Input GTCs were called with another cluster file than the model",
|
||||
"IN-050",
|
||||
),
|
||||
_m("W_PGX_REGION_FEW_PROBES", "warning", "PGx CN region has fewer than 10 unmasked probes", "PCN-008"),
|
||||
_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"),
|
||||
# --- 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,204 @@
|
||||
"""Targeted PGx copy number — spec PCN-001, PCN-009. Pure functions, no I/O.
|
||||
|
||||
Region value x = median over the region's unmasked probes of (sample-centred LRR − the probe's training reference),
|
||||
minus the sample's scale (median of the same quantity over non-target autosomal probes).
|
||||
|
||||
CN mixture [OURS; array intensity compression as in cnvPartition/PennCNV, LIT: Wang et al. 2007 Genome Res 17:1665]:
|
||||
component means µ_c = δ_r + γ_r · log2((c + 1) / 3) for CN c = 0 … 5 — a probe signal proportional to (c + β) with
|
||||
a background β = 1 copy-equivalent, compressed by γ_r. δ_r (the region's CN-2 location), γ_r, σ_r and the weights
|
||||
π_r are per region; γ_r has a Gaussian prior N(0.25, 0.08²). CN 0 may sit lower than the template (its own shift
|
||||
ε₀ ≤ 0 and SD): GSTT1 probes drop by ~2 at CN 0, GSTM1 probes by ~0.3. The weights are constrained to Hardy–Weinberg
|
||||
proportions of haplotype copy numbers (0 … 3 copies per haplotype, Dirichlet prior favouring one copy). Fitted by
|
||||
EM from a grid of starts (each observed cluster labelled CN 0 … 4, CN 0 on the template or on the lowest values).
|
||||
HWE and the template's spacing ratios identify the labels in multi-cluster regions (GSTM1 0/1/2 rather than 1/2/3,
|
||||
which would imply unobserved CN-0 samples). On GSA v4 PGx CNV probes γ_r ≈ 0.19–0.56 (K-028).
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import math
|
||||
from dataclasses import dataclass
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
F64 = npt.NDArray[np.float64]
|
||||
CN_MAX = 5 # "5 or more"
|
||||
CNS = np.arange(CN_MAX + 1)
|
||||
TEMPLATE = np.log2((CNS + 1) / 3.0)
|
||||
HAP_ALPHA = np.array([1.2, 3.0, 1.2, 1.05]) # Dirichlet prior on haplotype copy frequencies f_0..f_3 (f_1 favoured)
|
||||
SIGMA_FLOOR = 0.012
|
||||
GAMMA_PRIOR = (0.25, 0.08) # mean, SD
|
||||
GAMMA_RANGE = (0.15, 0.80)
|
||||
PI_FLOOR = 0.002
|
||||
CN0_PRIOR_N = 2.0
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class Region:
|
||||
chrom: str
|
||||
start0: int
|
||||
end: int
|
||||
name: str
|
||||
gene: str
|
||||
|
||||
|
||||
def parse_regions(text: str) -> list[Region]:
|
||||
"""PCN-001 BED: ``chrom start0 end region gene``; ``#`` comments."""
|
||||
out = []
|
||||
for ln in text.splitlines():
|
||||
if not ln.strip() or ln.startswith("#"):
|
||||
continue
|
||||
c, s, e, r, g = ln.split("\t")[:5]
|
||||
out.append(Region(c, int(s), int(e), r, g))
|
||||
return out
|
||||
|
||||
|
||||
@dataclass
|
||||
class RegionFit:
|
||||
delta: float
|
||||
gamma: float
|
||||
sigma: float
|
||||
weights: list[float]
|
||||
n: int
|
||||
cn0_shift: float = 0.0 # ≤ 0: extra drop of the CN-0 component below the template (probe-specific background)
|
||||
cn0_sigma: float = 0.0 # SD of the CN-0 component (≥ sigma; 0 = sigma)
|
||||
|
||||
def means(self) -> F64:
|
||||
m = self.delta + self.gamma * TEMPLATE
|
||||
m[0] += self.cn0_shift
|
||||
return m
|
||||
|
||||
def sds(self) -> F64:
|
||||
s = np.full(len(CNS), self.sigma)
|
||||
s[0] = max(self.sigma, self.cn0_sigma)
|
||||
return s
|
||||
|
||||
|
||||
def hwe_weights(f: F64) -> F64:
|
||||
"""CN weights under Hardy–Weinberg from haplotype copy-number frequencies f_k (k = 0 … 3 copies):
|
||||
π_c = Σ_{i+j=c} f_i f_j, CN ≥ 5 pooled into 5 [LIT: HWE; OURS: use as the CN mixture prior]."""
|
||||
pi = np.zeros(len(CNS))
|
||||
for i in range(len(f)):
|
||||
for j in range(len(f)):
|
||||
pi[min(i + j, CN_MAX)] += f[i] * f[j]
|
||||
return pi
|
||||
|
||||
|
||||
def _hap_update(r: F64, f: F64) -> F64:
|
||||
"""M-step for f: split each sample's CN responsibilities over the haplotype pairs (i, j) with i + j = c."""
|
||||
counts = HAP_ALPHA - 1.0
|
||||
counts = counts.astype(np.float64).copy()
|
||||
pi = hwe_weights(f)
|
||||
for i in range(len(f)):
|
||||
for j in range(len(f)):
|
||||
c = min(i + j, CN_MAX)
|
||||
share = f[i] * f[j] / max(pi[c], 1e-300)
|
||||
w = float(r[:, c].sum()) * share
|
||||
counts[i] += w
|
||||
counts[j] += w
|
||||
counts = np.maximum(counts, 1e-9)
|
||||
out: F64 = counts / counts.sum()
|
||||
return out
|
||||
|
||||
|
||||
def _norm_pdf(x: F64, mu: F64, sd: float | F64) -> F64:
|
||||
sd_ = np.broadcast_to(np.asarray(sd, dtype=np.float64), mu.shape)[None, :]
|
||||
return np.exp(-0.5 * ((x[:, None] - mu[None, :]) / sd_) ** 2) / (sd_ * math.sqrt(2 * math.pi))
|
||||
|
||||
|
||||
def _em_region(x: F64, gamma: float, delta: float, iters: int, eps0: float = 0.0) -> tuple[RegionFit, float]:
|
||||
"""EM for one region; returns the fit and its penalised log-likelihood (HWE weights from haplotype copy
|
||||
frequencies with a Dirichlet prior, Gaussian prior on γ).
|
||||
CN 0 has its own shift ε₀ ≤ 0 below the template and its own SD, both shrunk towards the template with a prior
|
||||
weight of ``CN0_PRIOR_N`` samples (on some probe sets CN 0 is far lower than the compressed template, GSTT1)."""
|
||||
n = len(x)
|
||||
f = np.array([0.1, 0.8, 0.07, 0.03])
|
||||
pi = hwe_weights(f)
|
||||
sd = max(SIGMA_FLOOR, float(np.std(x)) if n > 1 else 0.05)
|
||||
sd0 = sd
|
||||
if eps0 < 0: # a start with a far CN-0 cluster: the other components get the spread of the rest
|
||||
rest = x[x > delta + gamma * TEMPLATE[0] + eps0 / 2]
|
||||
sd = max(SIGMA_FLOOR, float(np.std(rest)) if len(rest) > 1 else 0.05)
|
||||
sd0 = sd
|
||||
g0, gs = GAMMA_PRIOR
|
||||
ll = -math.inf
|
||||
for _ in range(iters):
|
||||
mu = delta + gamma * TEMPLATE
|
||||
mu[0] += eps0
|
||||
sds = np.full(len(CNS), sd)
|
||||
sds[0] = sd0
|
||||
dens = pi[None, :] * _norm_pdf(x, mu, sds) + 1e-300
|
||||
tot = dens.sum(axis=1, keepdims=True)
|
||||
ll = float(np.log(tot).sum() + ((HAP_ALPHA - 1) * np.log(np.maximum(f, 1e-12))).sum()
|
||||
- 0.5 * ((gamma - g0) / gs) ** 2) # fmt: skip
|
||||
r = dens / tot
|
||||
f = _hap_update(r, f)
|
||||
pi = hwe_weights(f)
|
||||
# δ and γ: weighted least squares of x on (1, t_c) over CN ≥ 1, γ shrunk towards its prior
|
||||
rr = r.copy()
|
||||
rr[:, 0] = 0.0
|
||||
w = rr.sum(axis=0)
|
||||
sw, st, stt = float(w.sum()), float((w * TEMPLATE).sum()), float((w * TEMPLATE**2).sum())
|
||||
sx, stx = float((rr * x[:, None]).sum()), float((rr * x[:, None] * TEMPLATE[None, :]).sum())
|
||||
lam = sd**2 / gs**2
|
||||
det = sw * (stt + lam) - st * st
|
||||
if det > 0:
|
||||
delta = (sx * (stt + lam) - st * (stx + lam * g0)) / det
|
||||
gamma = float(np.clip((sw * (stx + lam * g0) - st * sx) / det, *GAMMA_RANGE))
|
||||
mu = delta + gamma * TEMPLATE
|
||||
n1 = float(w.sum())
|
||||
if n1 > 0:
|
||||
sd = max(SIGMA_FLOOR, math.sqrt(float((rr * (x[:, None] - mu[None, :]) ** 2).sum() / n1)))
|
||||
r0 = r[:, 0]
|
||||
n0 = float(r0.sum())
|
||||
eps0 = min(0.0, float((r0 * (x - mu[0])).sum()) / (n0 + CN0_PRIOR_N))
|
||||
v0 = float((r0 * (x - mu[0] - eps0) ** 2).sum())
|
||||
sd0 = max(sd, math.sqrt((v0 + CN0_PRIOR_N * sd**2) / (n0 + CN0_PRIOR_N)))
|
||||
return RegionFit(float(delta), float(gamma), sd, pi.tolist(), n, float(eps0), float(sd0)), ll
|
||||
|
||||
|
||||
def fit_region(x: F64, iters: int = 200) -> RegionFit:
|
||||
"""PCN-009 step 4 for one region: the best of EM runs started with each observed cluster labelled CN 0 … 4, each
|
||||
with CN 0 on the template or on the lowest values."""
|
||||
g0 = GAMMA_PRIOR[0]
|
||||
cands = sorted({round(float(q) - g0 * float(t), 4) for q in np.quantile(x, [0.1, 0.3, 0.5, 0.7, 0.9])
|
||||
for t in TEMPLATE[:5]}) # fmt: skip
|
||||
lo = float(np.quantile(x, 0.02))
|
||||
best: tuple[float, float, float] | None = None
|
||||
for d0 in cands:
|
||||
for e0 in sorted({0.0, round(min(0.0, lo - (d0 + g0 * float(TEMPLATE[0]))), 4)}):
|
||||
_, ll = _em_region(x, g0, d0, 80, e0)
|
||||
if best is None or ll > best[0] + 1e-9:
|
||||
best = (ll, d0, e0)
|
||||
assert best is not None
|
||||
return _em_region(x, g0, best[1], iters, best[2])[0]
|
||||
|
||||
|
||||
def fit_mixture(values: dict[str, F64], iters: int = 200) -> dict[str, RegionFit]:
|
||||
"""Per-region fits for every region with ≥ 3 finite values."""
|
||||
out = {}
|
||||
for name, x in values.items():
|
||||
x = x[np.isfinite(x)]
|
||||
if len(x) >= 3:
|
||||
out[name] = fit_region(x, iters)
|
||||
return out
|
||||
|
||||
|
||||
def posterior(x: float, fit: RegionFit) -> F64:
|
||||
"""P(CN = c | x) for c = 0 … 5, with the training weights (floored) as prior."""
|
||||
pi = np.maximum(np.asarray(fit.weights), PI_FLOOR)
|
||||
pi = pi / pi.sum()
|
||||
d = pi * _norm_pdf(np.array([x]), fit.means(), fit.sds())[0]
|
||||
s = float(d.sum())
|
||||
out: F64 = d / s if s > 0 else pi
|
||||
return out
|
||||
|
||||
|
||||
def call_cn(x: float, fit: RegionFit) -> tuple[int, int, F64]:
|
||||
"""PCN-009 steps 5–6: (CN, QUAL, posterior)."""
|
||||
p = posterior(x, fit)
|
||||
cn = int(np.argmax(p))
|
||||
pmax = float(p[cn])
|
||||
qual = 60 if pmax >= 1.0 else int(min(60, max(0, round(-10 * math.log10(max(1e-300, 1 - pmax))))))
|
||||
return cn, qual, p
|
||||
@@ -0,0 +1,11 @@
|
||||
# simurg PGx CN regions v1, GRCh38, from UCSC ncbiRefSeq (PCN-001)
|
||||
chr1 109687816 109693745 GSTM1 GSTM1
|
||||
chr4 68537172 68576322 UGT2B17 UGT2B17
|
||||
chr10 133527362 133539123 CYP2E1 CYP2E1
|
||||
chr16 28605257 28610044 SULT1A1 SULT1A1
|
||||
chr19 40844772 40845293 CYP2A6.intron.7 CYP2A6
|
||||
chr19 40849746 40850947 CYP2A6.exon.1 CYP2A6
|
||||
chr22 42126498 42126752 CYP2D6.exon.9 CYP2D6
|
||||
chr22 42129185 42129737 CYP2D6.intron.2 CYP2D6
|
||||
chr22 42130810 42132810 CYP2D6.p5 CYP2D6
|
||||
chr22_KI270879v1_alt 270313 278446 GSTT1 GSTT1
|
||||
@@ -0,0 +1,510 @@
|
||||
"""``simurg pgx copy-number train|call`` — targeted PGx copy number (spec PCN-001 … PCN-017, CLI-031/032, IN-050).
|
||||
|
||||
I/O and orchestration only; the numerics live in ``simurg.pgx.cn``.
|
||||
|
||||
Deviations (K-028): the per-probe reference is the training cohort's median, stored in the model (not the calling
|
||||
batch's median, which is not CN 2 in regions with common deletions such as GSTM1); no GC correction step (the region
|
||||
median over many probes plus the sample scale already removes most of it); GSTT1 needs alt-contig probe coordinates,
|
||||
which the GSA v4 PGx GRCh38 manifest carries (160 probes on `chr22_KI270879v1_alt`); manifests without them give
|
||||
`NoCall` for GSTT1.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import hashlib
|
||||
import io
|
||||
import json
|
||||
import math
|
||||
import zipfile
|
||||
from dataclasses import dataclass, field
|
||||
from datetime import UTC, datetime
|
||||
from importlib import resources
|
||||
from pathlib import Path
|
||||
from typing import Any
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
import pyarrow as pa
|
||||
import pyarrow.parquet as pq
|
||||
|
||||
from simurg import __version__
|
||||
from simurg.core.errors import ExitCode, SimurgError
|
||||
from simurg.core.log import Reporter
|
||||
from simurg.formats.bpm import read_bpm
|
||||
from simurg.pgx.cn import CN_MAX, Region, RegionFit, call_cn, fit_mixture, parse_regions
|
||||
|
||||
MODEL_SCHEMA = 1
|
||||
REGIONS_FILE = "pgx_cn_regions.GRCh38.v1.bed"
|
||||
N_SCALE_PROBES = 20_000
|
||||
MIN_REGION_PROBES = 10 # PCN-008
|
||||
MIN_FINITE_FRAC = 0.5 # PCN-002: fewer finite region values → NoCall [OURS]
|
||||
F32 = npt.NDArray[np.float32]
|
||||
|
||||
|
||||
def default_regions() -> tuple[str, str]:
|
||||
"""(file name, text) of the packaged region BED (PCN-001)."""
|
||||
return REGIONS_FILE, resources.files("simurg.pgx.data").joinpath(REGIONS_FILE).read_text(encoding="utf-8")
|
||||
|
||||
|
||||
def _core(chrom: str) -> str:
|
||||
c = chrom[3:] if chrom.lower().startswith("chr") else chrom
|
||||
return "X" if c == "XY" else c
|
||||
|
||||
|
||||
@dataclass
|
||||
class PgxCnModel:
|
||||
header: dict[str, Any]
|
||||
regions: list[Region]
|
||||
regions_text: str
|
||||
probe_idx: npt.NDArray[np.int64] # manifest indices of all probes used (region probes, then scale probes)
|
||||
probe_region: list[str] # region name, "" for scale probes
|
||||
probe_ref: F32 # training reference (median sample-centred LRR)
|
||||
fits: dict[str, RegionFit] = field(default_factory=dict)
|
||||
|
||||
def region_sel(self, name: str) -> npt.NDArray[np.int64]:
|
||||
return np.flatnonzero(np.asarray(self.probe_region) == name)
|
||||
|
||||
|
||||
def region_probes(
|
||||
regions: list[Region], chroms: list[str], pos: npt.NDArray[np.int64], masked: set[int]
|
||||
) -> dict[str, npt.NDArray[np.int64]]:
|
||||
"""Manifest indices inside each region (1-based position in (start0, end]), mask removed (PCN-008)."""
|
||||
core = np.array([_core(c) for c in chroms])
|
||||
out = {}
|
||||
for r in regions:
|
||||
sel = np.flatnonzero((core == _core(r.chrom)) & (pos > r.start0) & (pos <= r.end))
|
||||
out[r.name] = np.array([i for i in sel.tolist() if i not in masked], dtype=np.int64)
|
||||
return out
|
||||
|
||||
|
||||
def scale_probes(
|
||||
chroms: list[str], pos: npt.NDArray[np.int64], exclude: npt.NDArray[np.int64], n: int = N_SCALE_PROBES
|
||||
) -> npt.NDArray[np.int64]:
|
||||
"""Seeded random autosomal non-target probes for the per-sample scale (PCN-009 step 2)."""
|
||||
core = np.array([_core(c) for c in chroms])
|
||||
ok = np.isin(core, [str(i) for i in range(1, 23)]) & (pos > 0)
|
||||
ok[exclude] = False
|
||||
cand = np.flatnonzero(ok)
|
||||
rng = np.random.default_rng(20260929)
|
||||
return np.sort(rng.choice(cand, min(n, len(cand)), replace=False)).astype(np.int64)
|
||||
|
||||
|
||||
def _centred(lrr: npt.NDArray[np.float32], auto_idx: npt.NDArray[np.int64]) -> npt.NDArray[np.float64]:
|
||||
x = lrr.astype(np.float64)
|
||||
x[~np.isfinite(x)] = np.nan
|
||||
med = float(np.nanmedian(x[auto_idx])) if np.isfinite(x[auto_idx]).any() else 0.0
|
||||
return x - med
|
||||
|
||||
|
||||
def region_values(d: npt.NDArray[np.float64], model: PgxCnModel) -> dict[str, tuple[float, int]]:
|
||||
"""Per region: (x, n finite probes). ``d`` = centred LRR at ``model.probe_idx`` minus the reference."""
|
||||
scale_sel = model.region_sel("")
|
||||
sc = float(np.nanmedian(d[scale_sel])) if np.isfinite(d[scale_sel]).any() else 0.0
|
||||
out = {}
|
||||
for r in model.regions:
|
||||
sel = model.region_sel(r.name)
|
||||
v = d[sel]
|
||||
nf = int(np.isfinite(v).sum())
|
||||
out[r.name] = (float(np.nanmedian(v)) - sc if nf else math.nan, nf)
|
||||
return out
|
||||
|
||||
|
||||
# ------------------------------------------------------------------------------------------------ model file
|
||||
def write_model(path: Path, m: PgxCnModel) -> None:
|
||||
table = pa.table({"index": m.probe_idx, "region": m.probe_region, "ref": m.probe_ref,
|
||||
"name": m.header.pop("_names")}) # fmt: skip
|
||||
buf = io.BytesIO()
|
||||
pq.write_table(table, buf, compression="zstd")
|
||||
mixture = {name: {"delta": f.delta, "gamma": f.gamma, "sigma": f.sigma, "weights": f.weights, "n": f.n,
|
||||
"cn0_shift": f.cn0_shift, "cn0_sigma": f.cn0_sigma} for name, f in m.fits.items()} # fmt: skip
|
||||
path.parent.mkdir(parents=True, exist_ok=True)
|
||||
tmp = path.with_suffix(path.suffix + ".part")
|
||||
with zipfile.ZipFile(tmp, "w", compression=zipfile.ZIP_STORED) as z:
|
||||
z.writestr("header.json", json.dumps(m.header, indent=1))
|
||||
z.writestr("probes.parquet", buf.getvalue())
|
||||
z.writestr("regions.bed", m.regions_text)
|
||||
z.writestr("mixture.json", json.dumps(mixture, indent=1))
|
||||
tmp.replace(path)
|
||||
|
||||
|
||||
def read_model(path: Path) -> PgxCnModel:
|
||||
"""IN-050."""
|
||||
with zipfile.ZipFile(path) as z:
|
||||
header = json.loads(z.read("header.json"))
|
||||
table = pq.read_table(io.BytesIO(z.read("probes.parquet")))
|
||||
regions_text = z.read("regions.bed").decode()
|
||||
mixture = json.loads(z.read("mixture.json"))
|
||||
if int(header.get("schema_version", 0)) < MODEL_SCHEMA:
|
||||
raise SimurgError("E_MODEL_VERSION", f"{path}: PGx CN model schema {header.get('schema_version')} is too old",
|
||||
hint="retrain it with simurg pgx copy-number train") # fmt: skip
|
||||
fits = {
|
||||
k: RegionFit(
|
||||
v["delta"],
|
||||
v["gamma"],
|
||||
v["sigma"],
|
||||
list(v["weights"]),
|
||||
int(v["n"]),
|
||||
v.get("cn0_shift", 0.0),
|
||||
v.get("cn0_sigma", 0.0),
|
||||
)
|
||||
for k, v in mixture.items()
|
||||
}
|
||||
return PgxCnModel(header, parse_regions(regions_text), regions_text,
|
||||
table.column("index").to_numpy().astype(np.int64), table.column("region").to_pylist(),
|
||||
table.column("ref").to_numpy().astype(np.float32), fits) # fmt: skip
|
||||
|
||||
|
||||
# ------------------------------------------------------------------------------------------------ train
|
||||
@dataclass
|
||||
class PgxTrainOptions:
|
||||
gtcs: Path
|
||||
bpm: Path
|
||||
fasta: Path
|
||||
out: Path
|
||||
mask: Path | None = None
|
||||
regions: Path | None = None
|
||||
summary: Path | None = None
|
||||
platform: str = "auto"
|
||||
min_samples: int = 96
|
||||
min_passing: int = 90
|
||||
max_samples: int = 3000
|
||||
lrr_sd_max: float = 0.20
|
||||
call_rate_min: float = 0.99
|
||||
assay_control_min: float = 1.0
|
||||
force: bool = False
|
||||
argv: list[str] = field(default_factory=list)
|
||||
|
||||
|
||||
def read_mask_tsv(path: Path | None) -> tuple[set[str], str]:
|
||||
"""IN-052: ``probe_name<TAB>reason``; returns (names, sha256)."""
|
||||
if path is None:
|
||||
return set(), ""
|
||||
data = path.read_bytes()
|
||||
names = set()
|
||||
for ln in data.decode("utf-8").splitlines():
|
||||
if ln.strip() and not ln.startswith("#"):
|
||||
names.add(ln.split("\t")[0].strip())
|
||||
return names, hashlib.sha256(data).hexdigest()
|
||||
|
||||
|
||||
def _summary(path: Path | None) -> dict[str, dict[str, str]]:
|
||||
if path is None or not path.exists():
|
||||
return {}
|
||||
lines = path.read_text(encoding="utf-8").splitlines()
|
||||
head = lines[0].split("\t")
|
||||
return {r[0]: dict(zip(head, r, strict=False)) for r in (ln.split("\t") for ln in lines[1:] if ln)}
|
||||
|
||||
|
||||
def _assay_control(row: dict[str, str]) -> float:
|
||||
for k, v in row.items():
|
||||
if k.endswith("_norm_r") and k.startswith("TGA_Ctrl"):
|
||||
try:
|
||||
return float(v)
|
||||
except ValueError:
|
||||
return math.nan
|
||||
return math.nan
|
||||
|
||||
|
||||
def _probe_set_sha(names: list[str], chroms: list[str], pos: npt.NDArray[np.int64]) -> str:
|
||||
text = "\n".join(f"{n}\t{c}\t{p}" for n, c, p in zip(names, chroms, pos.tolist(), strict=True))
|
||||
return hashlib.sha256(text.encode()).hexdigest()
|
||||
|
||||
|
||||
def run_train(o: PgxTrainOptions, rep: Reporter) -> int:
|
||||
"""``simurg pgx copy-number train`` (CLI-032, PCN-005)."""
|
||||
from simurg.core.atomic import check_outputs_absent
|
||||
from simurg.formats.discovery import discover_gtcs
|
||||
from simurg.formats.fasta import Genome
|
||||
from simurg.formats.gtc import read_gtc
|
||||
from simurg.pipeline.qc import write_tsv
|
||||
|
||||
man = read_bpm(o.bpm)
|
||||
out_model = o.out / f"{o.bpm.stem}.simurg-pgxcn"
|
||||
check_outputs_absent([out_model, o.out / "train_exclusions.tsv"], force=o.force)
|
||||
genome = Genome(o.fasta)
|
||||
if genome.build not in ("GRCh38", "unknown", ""):
|
||||
raise SimurgError("E_BUILD_MISMATCH", f"PGx CN supports GRCh38 only (FASTA is {genome.build})") # PCN-015
|
||||
reg_name, reg_text = (o.regions.name, o.regions.read_text()) if o.regions else default_regions()
|
||||
regions = parse_regions(reg_text)
|
||||
mask_path = o.mask or (
|
||||
o.bpm.with_name(o.bpm.stem + ".mask.tsv") if (o.bpm.with_name(o.bpm.stem + ".mask.tsv")).exists() else None
|
||||
)
|
||||
mask_names, mask_sha = read_mask_tsv(mask_path)
|
||||
masked = {i for i, n in enumerate(man.names) if n in mask_names}
|
||||
pos = man.map_infos.astype(np.int64)
|
||||
rp = region_probes(regions, list(man.chroms), pos, masked)
|
||||
for r in regions:
|
||||
if len(rp[r.name]) < MIN_REGION_PROBES:
|
||||
rep.warn("W_PGX_REGION_FEW_PROBES", f"region {r.name}: {len(rp[r.name])} unmasked probe(s) on "
|
||||
f"{man.manifest_name}; it will be NoCall") # fmt: skip
|
||||
reg_idx = np.concatenate([rp[r.name] for r in regions]) if regions else np.zeros(0, np.int64)
|
||||
sc_idx = scale_probes(list(man.chroms), pos, reg_idx)
|
||||
probe_idx = np.concatenate([reg_idx, sc_idx])
|
||||
probe_region = [r.name for r in regions for _ in range(len(rp[r.name]))] + [""] * len(sc_idx)
|
||||
gtcs = discover_gtcs(o.gtcs)
|
||||
keys = sorted(gtcs)
|
||||
if len(keys) > o.max_samples:
|
||||
rng = np.random.default_rng(42)
|
||||
keys = sorted(rng.choice(keys, o.max_samples, replace=False).tolist())
|
||||
rep.info(f"training on a seeded random subset of {o.max_samples} of {len(gtcs)} samples (PCN-005)")
|
||||
summ = _summary(o.summary or (o.gtcs / "sample_summary.tsv"))
|
||||
rows = np.full((len(keys), len(probe_idx)), np.nan, dtype=np.float32)
|
||||
exclusions: list[dict[str, str]] = []
|
||||
passing = np.zeros(len(keys), dtype=np.bool_)
|
||||
cluster_files = set()
|
||||
for i, k in enumerate(keys):
|
||||
g = read_gtc(gtcs[k])
|
||||
if len(g.logr_ratios) != man.num_loci:
|
||||
raise SimurgError("E_MODEL_MANIFEST_MISMATCH", f"{gtcs[k]} has {len(g.logr_ratios)} loci, the manifest "
|
||||
f"{man.num_loci}") # fmt: skip
|
||||
cluster_files.add(g.cluster_file)
|
||||
rows[i] = _centred(g.logr_ratios, sc_idx)[probe_idx]
|
||||
reasons = []
|
||||
if not g.logr_dev <= o.lrr_sd_max:
|
||||
reasons.append(f"logr_dev {g.logr_dev:.3f} > {o.lrr_sd_max}")
|
||||
if not g.call_rate >= o.call_rate_min:
|
||||
reasons.append(f"call_rate {g.call_rate:.4f} < {o.call_rate_min}")
|
||||
ac = _assay_control(summ.get(k, {}))
|
||||
if math.isfinite(ac) and ac < o.assay_control_min:
|
||||
reasons.append(f"assay_control_norm_r {ac:.3f} < {o.assay_control_min}")
|
||||
if reasons:
|
||||
exclusions.append({"sample_id": k, "reason": "; ".join(reasons)})
|
||||
else:
|
||||
passing[i] = True
|
||||
n_pass = int(passing.sum())
|
||||
if len(keys) < o.min_samples or n_pass < o.min_passing:
|
||||
raise SimurgError("E_PGX_TRAIN_TOO_SMALL", f"{n_pass} of {len(keys)} sample(s) passed QC; training requires "
|
||||
f"at least {o.min_samples} samples of which {o.min_passing} pass",
|
||||
hint="lower --min-samples/--min-passing only for research use") # fmt: skip
|
||||
ref = np.nanmedian(rows[passing], axis=0)
|
||||
d = rows - ref
|
||||
sc_sel = np.flatnonzero(np.asarray(probe_region) == "")
|
||||
sc = np.nanmedian(d[:, sc_sel], axis=1)
|
||||
mad = float(np.nanmedian(np.abs(sc[passing] - np.nanmedian(sc[passing])))) * 1.4826 or 1e-9
|
||||
z = np.abs(sc - np.nanmedian(sc[passing])) / mad
|
||||
use = passing & (z <= 5) # PCN-009: outlier samples out of the reference
|
||||
for i in np.flatnonzero(passing & ~use).tolist():
|
||||
exclusions.append({"sample_id": keys[i], "reason": f"genome-wide NR outlier (|z| = {z[i]:.1f} > 5)"})
|
||||
if (~use & passing).any():
|
||||
ref = np.nanmedian(rows[use], axis=0)
|
||||
d = rows - ref
|
||||
model = PgxCnModel({}, regions, reg_text, probe_idx, probe_region, ref.astype(np.float32))
|
||||
values = {}
|
||||
for r in regions:
|
||||
if len(rp[r.name]) < MIN_REGION_PROBES:
|
||||
continue
|
||||
values[r.name] = np.array([region_values(d[i], model)[r.name][0] for i in np.flatnonzero(use)])
|
||||
model.fits = fit_mixture(values)
|
||||
csv_fmt = "unknown" if o.platform == "auto" else o.platform
|
||||
model.header = {
|
||||
"schema": "simurg.pgx_cn_model", "schema_version": MODEL_SCHEMA, "software_version": __version__,
|
||||
"manifest_name": man.manifest_name, "manifest_revision": getattr(man, "version", None),
|
||||
"probe_set_sha256": _probe_set_sha(man.names, list(man.chroms), pos),
|
||||
"n_probes": man.num_loci,
|
||||
"cluster_file_name": sorted(cluster_files)[0] if len(cluster_files) == 1 else sorted(cluster_files),
|
||||
"genome_build": "GRCh38", "fasta_contigs_sha256": genome.sha256, "platform": csv_fmt,
|
||||
"mask_sha256": mask_sha, "regions_bed_version": reg_name, "n_training_samples": int(use.sum()),
|
||||
"training_exclusions": len(exclusions), "created": datetime.now(UTC).strftime("%Y-%m-%dT%H:%M:%SZ"),
|
||||
"regions_few_probes": [r.name for r in regions if len(rp[r.name]) < MIN_REGION_PROBES],
|
||||
"model": "mu_c = delta_r + gamma_r * log2((c + 1) / 3), c = 0..5 (simurg.pgx.cn)",
|
||||
"_names": [man.names[i] for i in probe_idx.tolist()],
|
||||
} # fmt: skip
|
||||
o.out.mkdir(parents=True, exist_ok=True)
|
||||
write_model(out_model, model)
|
||||
write_tsv(o.out / "train_exclusions.tsv", ["sample_id", "reason"], exclusions)
|
||||
for name, f in model.fits.items():
|
||||
rep.info(f"{name}: CN2 at {f.delta:+.3f}, gamma {f.gamma:.3f}, sd {f.sigma:.3f}, "
|
||||
f"weights {' '.join(f'{w:.2f}' for w in f.weights)}") # fmt: skip
|
||||
rep.info(f"PGx CN model from {int(use.sum())} sample(s) written to {out_model}")
|
||||
return ExitCode.OK
|
||||
|
||||
|
||||
# ------------------------------------------------------------------------------------------------ call
|
||||
@dataclass
|
||||
class PgxCallOptions:
|
||||
gtcs: Path
|
||||
model: Path
|
||||
fasta: Path
|
||||
out: Path
|
||||
min_batch: int = 24
|
||||
min_passing: int = 22
|
||||
recommended_passing: int = 96
|
||||
lrr_sd_max: float = 0.20
|
||||
cluster_mismatch: str = "error"
|
||||
bgzip: bool = True
|
||||
force: bool = False
|
||||
argv: list[str] = field(default_factory=list)
|
||||
|
||||
|
||||
def _vcf_header(sample: str, contigs: list[tuple[str, int]], meta: dict[str, str]) -> str:
|
||||
h = ["##fileformat=VCFv4.3", f"##source=SimurgArray {__version__} pgx-cn", f"##reference=file://{meta['fasta']}",
|
||||
"##simurg_genome_build=GRCh38", f"##simurg_manifest={meta['manifest']}",
|
||||
f"##simurg_pgx_cn_model={meta['model']}", f"##simurg_batch={meta['batch']}",
|
||||
f"##OverallPloidy={meta['ploidy']}", "##GCCorrect=false"] # fmt: skip
|
||||
h += [f"##contig=<ID={c},length={n}>" for c, n in contigs]
|
||||
h += ['##ALT=<ID=CNV,Description="Copy number variable region (targeted genotype)">',
|
||||
'##FILTER=<ID=PASS,Description="All filters passed">',
|
||||
'##FILTER=<ID=Q7,Description="Quality below 7">',
|
||||
'##FILTER=<ID=SampleQuality,Description="Sample LRR SD above 0.2">',
|
||||
'##FILTER=<ID=NoCall,Description="Region could not be assessed">',
|
||||
'##INFO=<ID=END,Number=1,Type=Integer,Description="Last base of the region (1-based, inclusive)">',
|
||||
'##INFO=<ID=SVTYPE,Number=1,Type=String,Description="CNV">',
|
||||
'##INFO=<ID=CNVLEN,Number=1,Type=Integer,Description="Region length in bp">',
|
||||
'##INFO=<ID=PROBE,Number=1,Type=Integer,Description="Number of probes used for the region">',
|
||||
'##INFO=<ID=REGION,Number=1,Type=String,Description="Target region name">',
|
||||
'##INFO=<ID=GENE,Number=1,Type=String,Description="Gene symbol">',
|
||||
'##FORMAT=<ID=CN,Number=1,Type=Integer,Description="Copy number; 5 means 5 or more">',
|
||||
'##FORMAT=<ID=NR,Number=1,Type=Float,Description="Aggregated normalized intensity (about CN/2)">',
|
||||
"\t".join(["#CHROM", "POS", "ID", "REF", "ALT", "QUAL", "FILTER", "INFO", "FORMAT", sample])] # fmt: skip
|
||||
return "\n".join(h) + "\n"
|
||||
|
||||
|
||||
def call_sample_regions(d: npt.NDArray[np.float64], model: PgxCnModel, passing: bool) -> list[dict[str, Any]]:
|
||||
"""PCN-002/009: one result per region."""
|
||||
vals = region_values(d, model)
|
||||
out = []
|
||||
for r in model.regions:
|
||||
x, nf = vals[r.name]
|
||||
n_all = len(model.region_sel(r.name))
|
||||
fit = model.fits.get(r.name)
|
||||
filters: list[str] = []
|
||||
if fit is None or n_all < MIN_REGION_PROBES or nf < MIN_FINITE_FRAC * n_all or not math.isfinite(x):
|
||||
cn, qual, nr = None, 0, math.nan
|
||||
filters.append("NoCall")
|
||||
else:
|
||||
cn, qual, _ = call_cn(x, fit)
|
||||
# NR ≈ CN / 2: the continuous CN that the region's mixture means imply for x [OURS]
|
||||
nr = max(0.0, 3.0 * 2.0 ** ((x - fit.delta) / fit.gamma) - 1.0) / 2.0 if cn else 0.0
|
||||
if qual < 7:
|
||||
filters.append("Q7")
|
||||
if not passing:
|
||||
filters.append("SampleQuality")
|
||||
out.append(
|
||||
{"region": r, "cn": cn, "qual": qual, "nr": nr, "x": x, "probes": nf, "filters": filters or ["PASS"]}
|
||||
)
|
||||
return out
|
||||
|
||||
|
||||
def run_call(o: PgxCallOptions, rep: Reporter) -> int:
|
||||
"""``simurg pgx copy-number call`` (CLI-031)."""
|
||||
from simurg.core.atomic import atomic_write_text, check_outputs_absent
|
||||
from simurg.formats.discovery import discover_gtcs
|
||||
from simurg.formats.fasta import ContigMap, Genome
|
||||
from simurg.formats.gtc import read_gtc
|
||||
from simurg.pipeline.cn import write_vcf
|
||||
from simurg.pipeline.qc import write_tsv
|
||||
|
||||
model = read_model(o.model)
|
||||
genome = Genome(o.fasta)
|
||||
if genome.build not in ("GRCh38", "unknown", ""):
|
||||
raise SimurgError("E_BUILD_MISMATCH", f"PGx CN supports GRCh38 only (FASTA is {genome.build})")
|
||||
cmap = ContigMap([e.name for e in genome.entries], "auto")
|
||||
gtcs = discover_gtcs(o.gtcs)
|
||||
ext = ".pgx_cn.vcf.gz" if o.bgzip else ".pgx_cn.vcf"
|
||||
check_outputs_absent([o.out / f"{k}{ext}" for k in gtcs] + [o.out / "pgx_cn_batch_summary.tsv"], force=o.force)
|
||||
n_probes = int(model.header["n_probes"])
|
||||
mcl = model.header.get("cluster_file_name")
|
||||
# pass 1: batch QC from GTC headers (PCN-004)
|
||||
heads = {}
|
||||
for k, p in gtcs.items():
|
||||
g = read_gtc(p)
|
||||
heads[k] = (g.logr_dev, g.sample_name or k, g.cluster_file, len(g.logr_ratios))
|
||||
n_pass = sum(1 for v in heads.values() if v[0] < o.lrr_sd_max)
|
||||
if len(heads) < o.min_batch or n_pass < o.min_passing:
|
||||
raise SimurgError("E_PGX_BATCH_TOO_SMALL", f"{n_pass} sample(s) passed QC. Requires at least "
|
||||
f"{max(o.min_passing, 1)} samples to proceed.") # fmt: skip
|
||||
if n_pass < o.recommended_passing:
|
||||
rep.warn("W_PGX_BATCH_BELOW_RECOMMENDED", f"{n_pass} passing sample(s); {o.recommended_passing} recommended")
|
||||
for k, (_, _, cf, nl) in heads.items():
|
||||
if nl != n_probes:
|
||||
raise SimurgError("E_MODEL_MANIFEST_MISMATCH", f"{k} has {nl} loci, the model's manifest {n_probes}")
|
||||
if isinstance(mcl, str) and cf and cf != mcl:
|
||||
msg = f"Cluster file ({cf}) differs from the CN model's cluster file ({mcl})"
|
||||
if o.cluster_mismatch == "error":
|
||||
raise SimurgError("E_MODEL_CLUSTER_MISMATCH", msg, hint="retrain the model (PCN-005)")
|
||||
rep.warn("W_MODEL_CLUSTER_MISMATCH", msg)
|
||||
break
|
||||
contigs: list[tuple[str, int]] = []
|
||||
contig_name: dict[str, str | None] = {}
|
||||
for r in model.regions:
|
||||
if r.chrom in contig_name:
|
||||
continue
|
||||
c = cmap.get(r.chrom)
|
||||
contig_name[r.chrom] = c
|
||||
if c is None:
|
||||
rep.warn("W_PGX_NO_ALT_CONTIG", f"{r.chrom} is not in {o.fasta.name}; REF written as N (IN-062)")
|
||||
contigs.append((r.chrom, genome.contig_length(c) if c else r.end))
|
||||
order = {c: i for i, (c, _) in enumerate(contigs)}
|
||||
model_sha = hashlib.sha256(o.model.read_bytes()).hexdigest()
|
||||
mh = model.header
|
||||
meta = {
|
||||
"fasta": str(o.fasta.resolve()),
|
||||
"manifest": f"{mh['manifest_name']};sha256={mh['probe_set_sha256']}",
|
||||
"model": f"{o.model.name};sha256={model_sha};cluster_file={mcl};regions={mh['regions_bed_version']}",
|
||||
"batch": f"n_samples={len(heads)};n_passing={n_pass};lrr_sd_max={o.lrr_sd_max}",
|
||||
}
|
||||
centre_idx = model.probe_idx[model.region_sel("")] # the scale probes (autosomal), as in training
|
||||
o.out.mkdir(parents=True, exist_ok=True)
|
||||
srows, per_region = [], {r.name: {"n_gain": 0, "n_loss": 0, "n_neutral": 0, "n_nocall": 0, "n_filtered": 0}
|
||||
for r in model.regions} # fmt: skip
|
||||
failed = 0
|
||||
for k, p in gtcs.items():
|
||||
try:
|
||||
g = read_gtc(p)
|
||||
d = _centred(g.logr_ratios, centre_idx)[model.probe_idx] - model.probe_ref
|
||||
passing = heads[k][0] < o.lrr_sd_max
|
||||
res = call_sample_regions(d, model, passing)
|
||||
except SimurgError:
|
||||
raise
|
||||
except Exception as e: # PCN: one sample never stops the batch
|
||||
failed += 1
|
||||
rep.error("E_PGX_CN_FAILED", f"{k}: {e}", sample_id=k)
|
||||
continue
|
||||
called = [(x["region"], x["cn"]) for x in res if x["cn"] is not None]
|
||||
tot = sum(r.end - r.start0 for r, _ in called)
|
||||
ploidy = sum((r.end - r.start0) * min(cn, CN_MAX) for r, cn in called) / tot if tot else math.nan
|
||||
meta_s = {**meta, "ploidy": f"{ploidy:.3f}" if math.isfinite(ploidy) else "."}
|
||||
recs = []
|
||||
for x in sorted(res, key=lambda y: (order[y["region"].chrom], y["region"].start0)):
|
||||
r = x["region"]
|
||||
c = contig_name.get(r.chrom)
|
||||
ref = (genome.fetch(c, r.start0 - 1, r.start0) or "N").upper() if c and r.start0 >= 1 else "N"
|
||||
s, e_ = r.start0 + 1, r.end
|
||||
cn = "." if x["cn"] is None else str(x["cn"])
|
||||
nr = "." if not math.isfinite(x["nr"]) else f"{x['nr']:.6f}"
|
||||
line = "\t".join([r.chrom, str(s - 1), f"CNV:{r.name}:{r.chrom}:{s}-{e_}", ref, "<CNV>", str(x["qual"]),
|
||||
";".join(x["filters"]), f"END={e_};SVTYPE=CNV;CNVLEN={e_ - s + 1};PROBE={x['probes']};"
|
||||
f"REGION={r.name};GENE={r.gene}", "CN:NR", f"{cn}:{nr}"]) # fmt: skip
|
||||
recs.append((r.chrom, s - 2, e_, line))
|
||||
pr = per_region[r.name]
|
||||
if x["cn"] is None:
|
||||
pr["n_nocall"] += 1
|
||||
elif x["cn"] > 2:
|
||||
pr["n_gain"] += 1
|
||||
elif x["cn"] < 2:
|
||||
pr["n_loss"] += 1
|
||||
else:
|
||||
pr["n_neutral"] += 1
|
||||
pr["n_filtered"] += x["filters"] != ["PASS"]
|
||||
write_vcf(o.out / f"{k}{ext}", _vcf_header(heads[k][1], contigs, meta_s), recs, o.bgzip)
|
||||
srows.append({"sample_id": k, "sample_name": heads[k][1], "sample_folder": str(Path(p).parent),
|
||||
"logr_dev": f"{heads[k][0]:.6f}", "passed_batch_qc": str(passing).lower(),
|
||||
"used_in_reference": "false",
|
||||
**{f"cn_{x['region'].name}": ("." if x["cn"] is None else x["cn"]) for x in res},
|
||||
"flags": ""}) # fmt: skip
|
||||
cols = ["region", "gene", "n_gain", "n_loss", "n_neutral", "n_nocall", "n_filtered"]
|
||||
rows: list[dict[str, Any]] = [
|
||||
{"region": r.name, "gene": r.gene, **per_region[r.name]} for r in sorted(model.regions, key=lambda r: r.name)
|
||||
]
|
||||
totals = {c: sum(per_region[r.name][c] for r in model.regions) for c in cols[2:]}
|
||||
rows.append({"region": "ALL", "gene": "", **totals})
|
||||
write_tsv(o.out / "pgx_cn_batch_summary.tsv", cols, rows)
|
||||
scols = ["sample_id", "sample_name", "sample_folder", "logr_dev", "passed_batch_qc", "used_in_reference"]
|
||||
scols += [f"cn_{r.name}" for r in model.regions] + ["flags"]
|
||||
write_tsv(o.out / "pgx_cn_sample_summary.tsv", scols, srows)
|
||||
run = {"software": f"SimurgArray {__version__}", **meta, "samples": len(gtcs), "failed": failed}
|
||||
atomic_write_text(o.out / "pgx_cn_run.json", json.dumps(run, indent=1))
|
||||
rep.info(f"PGx CN: {len(srows)}/{len(gtcs)} sample(s) called")
|
||||
if failed == 0:
|
||||
return ExitCode.OK
|
||||
return ExitCode.PARTIAL if srows else ExitCode.ALL_FAILED
|
||||
Reference in new issue
Block a user