From baac9af5c5685173998259874edadb72c7cac31f Mon Sep 17 00:00:00 2001 From: Kemal Yaylali Date: Tue, 29 Sep 2026 17:11:34 +0100 Subject: [PATCH] PGx CN: scan-date batch warning (PCN-006), CN-consistent BAF evidence (PCN-010, --bafs-in-cn), noise-outlier exclusion + robustness test (PCN-017/EDGE-026), platform from the CSV manifest (CLI-032) Also: star-allele calls are checked for one manifest/build per run (PGX-011), duplications at the CN cap are named xN (PGX-022), and tandem names attach supporting variants to their full copy. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01814GKVzZVL3FuF3DeiZXhv --- docs/messages.md | 1 + docs/validation/V2_pgx_cn.md | 13 ++++ src/simurg/cli/main.py | 10 ++- src/simurg/core/errors.py | 1 + src/simurg/pgx/cn.py | 32 ++++++++- src/simurg/pipeline/pgx_cn.py | 116 +++++++++++++++++++++++++++----- src/simurg/pipeline/pgx_star.py | 18 ++++- tests/unit/test_pgx_cn.py | 50 +++++++++++++- 8 files changed, 216 insertions(+), 25 deletions(-) diff --git a/docs/messages.md b/docs/messages.md index f2453f3..17635b5 100644 --- a/docs/messages.md +++ b/docs/messages.md @@ -46,6 +46,7 @@ Generated from `src/simurg/core/errors.py` — do not edit by hand. | `E_PGX_CN_FAILED` | error | PGx CN calling failed for the sample | PCN-002 | | `E_PGX_CN_MISSING` | error | A sample has no PGx CN VCF (give --no-cn to run without) | PGX-011 | | `E_PGX_GUIDELINE` | error | Guideline tables not present in the PGx DB | PGX-025 | +| `E_PGX_PLATFORM_UNKNOWN` | error | PGx CN platform (LCG/EX) unknown: give --csv or --platform | CLI-032 | | `E_PGX_SOURCE` | error | PGx source tables missing; run scripts/pgx/fetch_sources.sh | PGX-001 | | `E_PGX_TRAIN_TOO_SMALL` | error | Too few (passing) samples to train a PGx CN model | PCN-005 | | `E_PGX_VCF_UNINDEXED` | error | SNV VCFs for PGx must be bgzipped and indexed | PGX-011 | diff --git a/docs/validation/V2_pgx_cn.md b/docs/validation/V2_pgx_cn.md index 8b22b01..2b47c45 100644 --- a/docs/validation/V2_pgx_cn.md +++ b/docs/validation/V2_pgx_cn.md @@ -39,6 +39,19 @@ literature ≈ 0.4–0.5). Details: `V2_pgx_cn_demo36.json`. * PCN-017 target (≥ 95 % per region) is met for GSTM1, UGT2B17, SULT1A1, CYP2A6.exon.1; not (yet) shown for the CYP2D6 sub-regions. +## Additions + +* PCN-006: a warning when the batch's IDAT scan dates span more than 14 days (no sample-sheet input in `pgx + copy-number call`, so the lot/batch-column check does not apply). +* PCN-010: `pgx_cn_baf_evidence.tsv` — per region the log-likelihood of CN 0 … 5 from the informative SNP BAFs + (binomial weights from the training cohort's B-allele frequency); `--bafs-in-cn` adds it to the CN posterior + (experimental, off by default). On the demo it agrees poorly with the LRR calls (e.g. GSTM1 2/36): the SNPs in + these regions do not separate the CN states here, so it stays informational. +* PCN-009 / EDGE-026: training excludes samples whose genome-wide level **or** noise is an outlier (robust |z| > 5); + `test_outlier_robustness` injects 10 % noisy, shifted arrays and requires ≥ 99 % unchanged calls for the others. +* CLI-032: the platform (LCG/EX) comes from the CSV manifest's `Assay Format` (`--csv`), or `--platform`; the demo + model records `EX`. + ## Caveats (K-028) * Trained and evaluated on the same 30 individuals (the spec requires ≥ 96 training samples; the priors were diff --git a/src/simurg/cli/main.py b/src/simurg/cli/main.py index 9fd484b..685af70 100644 --- a/src/simurg/cli/main.py +++ b/src/simurg/cli/main.py @@ -935,6 +935,7 @@ def pgx_cn_train_cmd( 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, + csv: Annotated[Path | None, typer.Option("--csv", help="CSV manifest (Assay Format → platform).")] = 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, @@ -954,7 +955,8 @@ def pgx_cn_train_cmd( 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, + platform=platform, csv=csv, 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)) @@ -971,6 +973,9 @@ def pgx_cn_call_cmd( 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", + bafs_in_cn: Annotated[ + bool, typer.Option("--bafs-in-cn", help="Experimental: add CN-consistent BAF evidence to the CN posterior.") + ] = False, no_bgzip: Annotated[bool, typer.Option("--no-bgzip")] = False, force: ForceOpt = False, log_json: LogJson = False, @@ -987,7 +992,8 @@ def pgx_cn_call_cmd( 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 + cluster_mismatch=cluster_mismatch, bafs_in_cn=bafs_in_cn, bgzip=not no_bgzip, + force=force) # fmt: skip run_guarded(rep, lambda: run_call(o, rep)) diff --git a/src/simurg/core/errors.py b/src/simurg/core/errors.py index 3c9166c..dffe020 100644 --- a/src/simurg/core/errors.py +++ b/src/simurg/core/errors.py @@ -176,6 +176,7 @@ CATALOGUE: dict[str, MessageDef] = dict( _m("E_PGX_VCF_UNSUITABLE", "error", "SNV VCF unsuitable for PGx (squashing/locus filters/foreign)", "PGX-011"), _m("E_PGX_GUIDELINE", "error", "Guideline tables not present in the PGx DB", "PGX-025"), _m("W_PGX_NO_PROBE_MAP", "warning", "The PGx DB has no probe map for this manifest", "PGX-035"), + _m("E_PGX_PLATFORM_UNKNOWN", "error", "PGx CN platform (LCG/EX) unknown: give --csv or --platform", "CLI-032"), _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"), diff --git a/src/simurg/pgx/cn.py b/src/simurg/pgx/cn.py index a9fc740..2590edf 100644 --- a/src/simurg/pgx/cn.py +++ b/src/simurg/pgx/cn.py @@ -195,9 +195,37 @@ def posterior(x: float, fit: RegionFit) -> F64: return out -def call_cn(x: float, fit: RegionFit) -> tuple[int, int, F64]: - """PCN-009 steps 5–6: (CN, QUAL, posterior).""" +def baf_loglik(baf: F64, pfb: F64 | None = None, sd: float = 0.05) -> F64: + """PCN-010: per-CN-state log-likelihood of a region's SNP BAFs: for CN c the BAF is k/c with binomial weights + C(c, k) p^k (1 − p)^(c − k) from the probe's population B-allele frequency p (0.5 when unknown); CN 0 → uniform + [LIT: Wang et al. 2007 (PennCNV BAF emissions); OURS: region-level use].""" + ok = np.isfinite(baf) + b = baf[ok] + p = np.clip(np.nan_to_num(pfb[ok] if pfb is not None else np.full(len(b), 0.5), nan=0.5), 0.01, 0.99) + out = np.zeros(len(CNS)) + if len(b) == 0: + return out + for c in CNS: + if c == 0: + out[c] = 0.0 # uniform on [0, 1] + continue + k = np.arange(c + 1) + means = k / c + w = np.array([math.comb(c, int(j)) for j in k])[None, :] * p[:, None] ** k[None, :] * (1 - p[:, None]) ** ( + c - k[None, :]) # fmt: skip + dens = np.exp(-0.5 * ((b[:, None] - means[None, :]) / sd) ** 2) / (sd * math.sqrt(2 * math.pi)) + out[c] = float(np.log(0.98 * (w * dens).sum(axis=1) + 0.02).sum()) + return out + + +def call_cn(x: float, fit: RegionFit, baf_ll: F64 | None = None) -> tuple[int, int, F64]: + """PCN-009 steps 5–6: (CN, QUAL, posterior); ``baf_ll`` (PCN-010, experimental ``--bafs-in-cn``) is added to the + log posterior.""" p = posterior(x, fit) + if baf_ll is not None: + lp = np.log(np.maximum(p, 1e-300)) + baf_ll + p = np.exp(lp - lp.max()) + p = p / p.sum() 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)))))) diff --git a/src/simurg/pipeline/pgx_cn.py b/src/simurg/pipeline/pgx_cn.py index 0ca78e1..f329bcc 100644 --- a/src/simurg/pipeline/pgx_cn.py +++ b/src/simurg/pipeline/pgx_cn.py @@ -31,7 +31,7 @@ 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 +from simurg.pgx.cn import CN_MAX, Region, RegionFit, baf_loglik, call_cn, fit_mixture, parse_regions MODEL_SCHEMA = 1 REGIONS_FILE = "pgx_cn_regions.GRCh38.v1.bed" @@ -60,6 +60,7 @@ class PgxCnModel: 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) + probe_pfb: F32 | None = None # training B-allele frequency per probe (NaN: no called genotypes), PCN-010 def region_sel(self, name: str) -> npt.NDArray[np.int64]: return np.flatnonzero(np.asarray(self.probe_region) == name) @@ -111,8 +112,10 @@ def region_values(d: npt.NDArray[np.float64], model: PgxCnModel) -> dict[str, tu # ------------------------------------------------------------------------------------------------ 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 + cols = {"index": m.probe_idx, "region": m.probe_region, "ref": m.probe_ref, "name": m.header.pop("_names")} + if m.probe_pfb is not None: + cols["pfb"] = m.probe_pfb + table = pa.table(cols) 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, @@ -151,7 +154,9 @@ def read_model(path: Path) -> PgxCnModel: } 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 + table.column("ref").to_numpy().astype(np.float32), fits, + table.column("pfb").to_numpy().astype(np.float32) + if "pfb" in table.column_names else None) # fmt: skip # ------------------------------------------------------------------------------------------------ train @@ -165,6 +170,7 @@ class PgxTrainOptions: regions: Path | None = None summary: Path | None = None platform: str = "auto" + csv: Path | None = None # CSV manifest: source of the Assay Format heading (CLI-032) min_samples: int = 96 min_passing: int = 90 max_samples: int = 3000 @@ -219,6 +225,17 @@ def run_train(o: PgxTrainOptions, rep: Reporter) -> int: from simurg.pipeline.qc import write_tsv man = read_bpm(o.bpm) + csv_fmt = o.platform + if o.platform == "auto": + from simurg.formats.csv_manifest import assay_format, read_csv_manifest + + csv_fmt = assay_format(read_csv_manifest(o.csv)) if o.csv else "" + if csv_fmt not in ("LCG", "EX"): + raise SimurgError( + "E_PGX_PLATFORM_UNKNOWN", + "cannot tell the platform (LCG/EX)" + (f": Assay Format {csv_fmt!r}" if csv_fmt else ": no --csv"), + hint="give --csv or --platform LCG|EX", + ) 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) @@ -249,6 +266,8 @@ def run_train(o: PgxTrainOptions, rep: Reporter) -> int: 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) + alt_count = np.zeros(len(probe_idx)) + n_called = np.zeros(len(probe_idx)) exclusions: list[dict[str, str]] = [] passing = np.zeros(len(keys), dtype=np.bool_) cluster_files = set() @@ -259,6 +278,9 @@ def run_train(o: PgxTrainOptions, rep: Reporter) -> int: f"{man.num_loci}") # fmt: skip cluster_files.add(g.cluster_file) rows[i] = _centred(g.logr_ratios, sc_idx)[probe_idx] + gt = np.asarray(g.genotypes)[probe_idx] + n_called += gt > 0 + alt_count += np.where(gt == 2, 1, np.where(gt == 3, 2, 0)) reasons = [] if not g.logr_dev <= o.lrr_sd_max: reasons.append(f"logr_dev {g.logr_dev:.3f} > {o.lrr_sd_max}") @@ -280,22 +302,32 @@ def run_train(o: PgxTrainOptions, rep: Reporter) -> int: 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 + + def robust_z(v: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]: + med = float(np.nanmedian(v[passing])) + mad = float(np.nanmedian(np.abs(v[passing] - med))) * 1.4826 or 1e-9 + out: npt.NDArray[np.float64] = np.abs(v - med) / mad + return out + + # PCN-009 / EDGE-026: outlier samples out of the reference — genome-wide level or genome-wide noise + noise = 1.4826 * np.nanmedian(np.abs(d[:, sc_sel] - sc[:, None]), axis=1) + z, zn = robust_z(sc), robust_z(noise) + use = passing & (z <= 5) & (zn <= 5) 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)"}) + exclusions.append({"sample_id": keys[i], "reason": f"genome-wide NR outlier (level |z| = {z[i]:.1f}, " + f"noise |z| = {zn[i]:.1f}; limit 5)"}) # fmt: skip 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)) + with np.errstate(invalid="ignore", divide="ignore"): + pfb = np.where(n_called > 0, alt_count / (2 * n_called), np.nan).astype(np.float32) + model = PgxCnModel({}, regions, reg_text, probe_idx, probe_region, ref.astype(np.float32), probe_pfb=pfb) 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), @@ -331,6 +363,7 @@ class PgxCallOptions: recommended_passing: int = 96 lrr_sd_max: float = 0.20 cluster_mismatch: str = "error" + bafs_in_cn: bool = False # PCN-010 experimental bgzip: bool = True force: bool = False argv: list[str] = field(default_factory=list) @@ -359,8 +392,14 @@ def _vcf_header(sample: str, contigs: list[tuple[str, int]], meta: dict[str, str 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.""" +def call_sample_regions( + d: npt.NDArray[np.float64], + model: PgxCnModel, + passing: bool, + baf: npt.NDArray[np.float64] | None = None, + bafs_in_cn: bool = False, +) -> list[dict[str, Any]]: + """PCN-002/009/010: one result per region; ``baf`` = the sample's BAF at ``model.probe_idx``.""" vals = region_values(d, model) out = [] for r in model.regions: @@ -368,23 +407,44 @@ def call_sample_regions(d: npt.NDArray[np.float64], model: PgxCnModel, passing: n_all = len(model.region_sel(r.name)) fit = model.fits.get(r.name) filters: list[str] = [] + bll = None + n_snp = 0 + if baf is not None: + sel = model.region_sel(r.name) + rb = baf[sel] + pf = model.probe_pfb[sel] if model.probe_pfb is not None else None + if pf is not None: + rb = np.where(np.isfinite(pf) & (pf > 0.01) & (pf < 0.99), rb, np.nan) # informative SNPs only + n_snp = int(np.isfinite(rb).sum()) + bll = baf_loglik(rb, pf) 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) + cn, qual, _ = call_cn(x, fit, bll if bafs_in_cn and n_snp else None) # 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"]} - ) + out.append({"region": r, "cn": cn, "qual": qual, "nr": nr, "x": x, "probes": nf, "filters": filters or ["PASS"], + "n_snp": n_snp, + "baf_ll": None if bll is None or not n_snp else [round(float(v), 3) for v in bll]}) # fmt: skip return out +def _parse_date(text: str) -> Any: + from datetime import datetime as _dt + + for fmt in ("%m/%d/%Y %I:%M:%S %p", "%m/%d/%Y %H:%M:%S", "%Y-%m-%dT%H:%M:%S", "%Y-%m-%d %H:%M:%S"): + try: + return _dt.strptime(text.strip(), fmt) + except ValueError: + continue + return None + + 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 @@ -406,13 +466,20 @@ def run_call(o: PgxCallOptions, rep: Reporter) -> int: mcl = model.header.get("cluster_file_name") # pass 1: batch QC from GTC headers (PCN-004) heads = {} + scan_dates: list[Any] = [] 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)) + scan_dates.append(_parse_date(g.imaging_date)) 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 + dates = [x for x in scan_dates if x is not None] + if dates and (max(dates) - min(dates)).days > 14: # PCN-006 (sample-sheet lot/batch columns: no sheet input here) + span = (max(dates) - min(dates)).days + rep.warn("W_PGX_BATCH_COMPOSITION", f"scan dates span {span} days ({min(dates):%Y-%m-%d} … " + f"{max(dates):%Y-%m-%d}); batch statistics may mix processing batches") # 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(): @@ -445,6 +512,7 @@ def run_call(o: PgxCallOptions, rep: Reporter) -> int: } centre_idx = model.probe_idx[model.region_sel("")] # the scale probes (autosomal), as in training o.out.mkdir(parents=True, exist_ok=True) + baf_rows: list[dict[str, Any]] = [] 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 @@ -453,13 +521,25 @@ def run_call(o: PgxCallOptions, rep: Reporter) -> int: 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) + bafs = None + if len( + g.b_allele_freqs + ): # PCN-010: SNP probes only (called genotypes; intensity-only CNV probes never are) + bafs = np.asarray(g.b_allele_freqs, dtype=np.float64)[model.probe_idx] + bafs[np.asarray(g.genotypes)[model.probe_idx] == 0] = np.nan + res = call_sample_regions(d, model, passing, bafs, o.bafs_in_cn) 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 + for x in res: + if x["baf_ll"] is not None: + bl = x["baf_ll"] + baf_rows.append({"sample_id": k, "region": x["region"].name, "n_snp": x["n_snp"], + **{f"ll_cn{c}": bl[c] for c in range(len(bl))}, + "best_cn_by_baf": int(np.argmax(bl[1:]) + 1)}) # fmt: skip 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 @@ -502,6 +582,8 @@ def run_call(o: PgxCallOptions, rep: Reporter) -> int: 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) + write_tsv(o.out / "pgx_cn_baf_evidence.tsv", ["sample_id", "region", "n_snp"] + [f"ll_cn{c}" for c in range(6)] + + ["best_cn_by_baf"], baf_rows) # PCN-010: for the star-allele caller and review # fmt: skip 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") diff --git a/src/simurg/pipeline/pgx_star.py b/src/simurg/pipeline/pgx_star.py index e2dc116..b17ba0d 100644 --- a/src/simurg/pipeline/pgx_star.py +++ b/src/simurg/pipeline/pgx_star.py @@ -11,6 +11,7 @@ import contextlib import gzip import json import math +import re import sqlite3 from collections import defaultdict from dataclasses import dataclass, field @@ -205,7 +206,9 @@ def solution_detail(m: GeneModel, obs: Observations, recs: list[dict[str, Any] | hap_of = {m.haplotypes[i]: i for i in s.copies} reason = {NO_CALL: "no_call", MASKED: "masked", NOT_COVERED: "not_on_array"} for name in names: - core = name.split("x")[0] if "x" in name and name.split("x")[-1].isdigit() else name + unit = name.split("+")[-1] # "*68+*4" → the full copy "*4"; "*2x2"/"*1xN" → "*2"/"*1" + mx = re.fullmatch(r"(.+?)x(\d+|N)", unit) + core = mx.group(1) if mx else unit h = hap_of.get(core) sup, unres = [], [] if h is not None: @@ -300,6 +303,7 @@ def run_star_call(o: StarOptions, rep: Reporter) -> int: unres_rows: list[dict[str, Any]] = [] not_cov_rows: dict[tuple[str, str], dict[str, Any]] = {} calls_json = [] + first_ident: tuple[str, str] | None = None for sid, path in sorted(snv.items()): recs, meta = read_vcf_sites(path, wanted) opts = meta.get("simurg_options", "") @@ -313,6 +317,12 @@ def run_star_call(o: StarOptions, rep: Reporter) -> int: or kv.get("include_loci", "none").lower() != "none" ): # PGX-011 raise SimurgError("E_PGX_VCF_UNSUITABLE", f"{path}: made with unsquashed duplicates or a locus filter") + ident = (meta.get("simurg_manifest", "").split(";")[0], meta.get("simurg_genome_build", "")) + if first_ident is None: + first_ident = ident + elif ident != first_ident and not o.allow_foreign_vcf: # PGX-011 (4) + raise SimurgError("E_PGX_VCF_UNSUITABLE", f"{path}: manifest/build {ident} differs from {first_ident}", + hint="call one manifest and build per run") # fmt: skip cn = {} if o.no_cn or sid not in cnv else read_cn_vcf(cnv[sid]) masked = masks.get(meta.get("simurg_manifest", "").split(";")[0], set()) sample_name = sid @@ -375,6 +385,12 @@ def call_one_gene( feasible_ratio=o.feasible_ratio, structure=structure) # fmt: skip if not sols: return None, [], [], [], [] + if gcn is not None and gcn >= 5: # PGX-022: CN 5 means "5 or more" → the duplication count is unknown (xN) + import re as _re + + for s_ in sols: + s_.haplotypes = tuple(_re.sub(r"x\d+$", "xN", h) for h in s_.haplotypes) + s_.members = [tuple(_re.sub(r"x\d+$", "xN", h) for h in mm) for mm in s_.members] top = sols[0] details = [solution_detail(m, obs, vrec, s_, top, o.low_gs, MAX_EDIT_DISTANCE) for s_ in sols] flags = [] diff --git a/tests/unit/test_pgx_cn.py b/tests/unit/test_pgx_cn.py index 471327e..574a24c 100644 --- a/tests/unit/test_pgx_cn.py +++ b/tests/unit/test_pgx_cn.py @@ -57,7 +57,7 @@ def test_packaged_regions() -> None: # PCN-001 assert next(r for r in regs if r.gene == "GSTT1").chrom == "chr22_KI270879v1_alt" -def _dataset(tmp: Path, n_samples: int = 40) -> tuple[Path, Path, Path, Path]: +def _dataset(tmp: Path, n_samples: int = 40, n_outliers: int = 0) -> tuple[Path, Path, Path, Path]: rng = np.random.default_rng(5) loci = [] for c in ("1", "2"): @@ -94,6 +94,16 @@ def _dataset(tmp: Path, n_samples: int = 40) -> tuple[Path, Path, Path, Path]: cluster_file="c.egt") # fmt: skip (gtcs / sid).mkdir(parents=True) write_gtc(g, gtcs / sid / f"{sid}.gtc") + for j in range(n_outliers): # EDGE-026: noisy, shifted arrays that still pass the LRR SD gate + sid = f"2000000099{j:02d}_R01C01" + lrr = rng.normal(0.3, 0.19, n) + g = Gtc(num_snps=n, genotypes=np.ones(n, np.uint8), genotype_scores=np.full(n, 0.9, np.float32), + b_allele_freqs=np.zeros(n, np.float32), logr_ratios=lrr.astype(np.float32), + raw_x=np.full(n, 1000, np.uint16), raw_y=np.full(n, 1000, np.uint16), + base_calls=np.zeros((n, 2), np.uint8), gender=b"F", call_rate=0.995, logr_dev=0.19, + cluster_file="c.egt") # fmt: skip + (gtcs / sid).mkdir(parents=True) + write_gtc(g, gtcs / sid / f"{sid}.gtc") (tmp / "truth.txt").write_text("\n".join(f"{k}\t{a}\t{b}" for k, (a, b) in truth.items())) return bpm, fa, gtcs, regions @@ -102,10 +112,11 @@ def test_pgx_cn_train_and_call(tmp_path: Path) -> None: # CLI-031/032, PCN-002/ bpm, fa, gtcs, regions = _dataset(tmp_path) run = CliRunner() r = run.invoke(app, ["pgx", "copy-number", "train", "--gtcs", str(gtcs), "--bpm", str(bpm), "--fasta", str(fa), - "--out", str(tmp_path / "model"), "--regions", str(regions), "--quiet"]) # fmt: skip + "--out", str(tmp_path / "model"), "--regions", str(regions), "--platform", "EX", + "--quiet"]) # fmt: skip assert r.exit_code == 2 and "E_PGX_TRAIN_TOO_SMALL" in r.output # 40 < 96 (PCN-005) r = run.invoke(app, ["pgx", "copy-number", "train", "--gtcs", str(gtcs), "--bpm", str(bpm), "--fasta", str(fa), - "--out", str(tmp_path / "model"), "--regions", str(regions), "--min-samples", "30", + "--out", str(tmp_path / "model"), "--regions", str(regions), "--min-samples", "30", "--platform", "EX", "--min-passing", "30", "--quiet"]) # fmt: skip assert r.exit_code == 0, r.output model = tmp_path / "model" / "m.simurg-pgxcn" @@ -143,3 +154,36 @@ def test_pgx_cn_train_and_call(tmp_path: Path) -> None: # CLI-031/032, PCN-002/ head = (out / "pgx_cn_sample_summary.tsv").read_text().splitlines()[0].split("\t") assert head[:6] == ["sample_id", "sample_name", "sample_folder", "logr_dev", "passed_batch_qc", "used_in_reference"] and head[-1] == "flags" # fmt: skip + + +def _calls(out: Path) -> dict[str, str]: + rows = [ln.split("\t") for ln in (out / "pgx_cn_sample_summary.tsv").read_text().splitlines()] + head = rows[0] + return {r[0]: ":".join(r[head.index("cn_GENEA")] + r[head.index("cn_GENEB")]) for r in rows[1:]} + + +def test_outlier_robustness(tmp_path: Path) -> None: # PCN-017 / EDGE-026: 10 % outliers leave other calls alone + run = CliRunner() + res = {} + for tag, n_out in (("clean", 0), ("outliers", 4)): + d = tmp_path / tag + d.mkdir() + bpm, fa, gtcs, regions = _dataset(d, 40, n_out) + r = run.invoke(app, ["pgx", "copy-number", "train", "--gtcs", str(gtcs), "--bpm", str(bpm), "--fasta", str(fa), + "--out", str(d / "model"), "--regions", str(regions), "--min-samples", "30", + "--min-passing", "30", "--platform", "EX", "--quiet"]) # fmt: skip + assert r.exit_code == 0, r.output + r = run.invoke(app, ["pgx", "copy-number", "call", "--gtcs", str(gtcs), "--model", str(d / "model" / "m.simurg-pgxcn"), + "--fasta", str(fa), "--out", str(d / "call"), "--quiet"]) # fmt: skip + assert r.exit_code == 0, r.output + res[tag] = _calls(d / "call") + assert (d / "call" / "pgx_cn_baf_evidence.tsv").exists() # PCN-010 output + same = sum(res["clean"][k] == res["outliers"][k] for k in res["clean"]) + assert same >= 0.99 * len(res["clean"]) + + +def test_platform_required(tmp_path: Path) -> None: # CLI-032: --platform auto needs a CSV manifest + bpm, fa, gtcs, regions = _dataset(tmp_path, 30) + r = CliRunner().invoke(app, ["pgx", "copy-number", "train", "--gtcs", str(gtcs), "--bpm", str(bpm), "--fasta", + str(fa), "--out", str(tmp_path / "m"), "--regions", str(regions), "--quiet"]) # fmt: skip + assert r.exit_code == 2 and "E_PGX_PLATFORM_UNKNOWN" in r.output