diff --git a/SIMURGARRAY_FUNCTIONAL_SPEC.md b/SIMURGARRAY_FUNCTIONAL_SPEC.md index 4e47c04..4a6e3ab 100644 --- a/SIMURGARRAY_FUNCTIONAL_SPEC.md +++ b/SIMURGARRAY_FUNCTIONAL_SPEC.md @@ -2447,6 +2447,11 @@ about **algorithms or numerics**, the plan wins. Each such conflict found during (CYTO-020) to a constitutional DEL with `HIGHFRACTION`. This catches constitutional deletions whose LRR is compressed towards the mosaic states. +- **CYTO-042** `§10.2 · v2` **Per-probe LRR noise [OURS].** A model built with `--reference-gtcs` stores each + probe's robust LRR SD across the reference (1.4826·MAD after removing each sample's median and the probe's + reference median). `cn call` scales the sample's LRR emission SD per probe by that spread relative to the autosomal + median spread, clipped to [0.5, 3]. Without the column, or with `--no-probe-lrr-sd`, every probe uses the sample SD. + --- ## 7. Targeted PGx copy number (v2) @@ -3723,6 +3728,7 @@ Generated from the item tags (`build.py`). Phases: P0–P10 as in plan §18; `v2 | CYTO-039 | Batch scale | §10, §19 | v2 | D-9, E-35 | | CYTO-040 | Degenerate samples | §10 | v2 | E-35 | | CYTO-041 | MF refinements | §10.3 | v2 | K-027 | +| CYTO-042 | Per-probe LRR noise | §10.2 | v2 | K-027 | | PCN-001 | Region catalogue | §11 | v2 | B-23, C-G25, D-10 | | PCN-002 | Genotype every region [EMU] | §11 | v2 | B-23 | | PCN-003 | Inputs | §11 | v2 | A-G37 | diff --git a/src/simurg/cli/main.py b/src/simurg/cli/main.py index 09e1955..38f9d6b 100644 --- a/src/simurg/cli/main.py +++ b/src/simurg/cli/main.py @@ -739,6 +739,9 @@ def cn_call_cmd( switch_distance: Annotated[ float, typer.Option("--switch-distance", help="Advanced: HMM switch distance D (bp).") ] = 5e8, + probe_lrr_sd: Annotated[ + bool, typer.Option("--probe-lrr-sd/--no-probe-lrr-sd", help="Advanced: per-probe LRR noise from the model.") + ] = True, threads: ThreadsOpt = 0, force: ForceOpt = False, log_json: LogJson = False, @@ -768,7 +771,7 @@ def cn_call_cmd( out=out, bpm=bpm, model=model, gtcs=gtcs, summary=summary, fasta=fasta, thresholds=th, smoothing=smoothing, cn_cap=cn_cap, gc_correct=gc_correct, unknown_sex_policy=unknown_sex_policy, gt_mode=gt_mode, svlen_mode=svlen_mode, mask_bed=mask_bed, bgzip=not no_bgzip, mosaic_prior=mosaic_prior, - switch_d=switch_distance, threads=threads, force=force, + switch_d=switch_distance, probe_lrr_sd=probe_lrr_sd, threads=threads, force=force, argv=sys.argv[1:], ) # fmt: skip run_guarded(rep, lambda: run_cn_call(opts, rep)) diff --git a/src/simurg/cn/hmm.py b/src/simurg/cn/hmm.py index 3ee8fb1..f685d61 100644 --- a/src/simurg/cn/hmm.py +++ b/src/simurg/cn/hmm.py @@ -8,7 +8,8 @@ weights, switch model] * Diploid regions: CN0, CN1, CN2, CN2-LOH, CN3, CN3-LOH, CN4 and mosaic loss / gain / CN-LOH at cell fractions 0.2, 0.4, 0.6, 0.8. Haploid regions (male X outside PAR, male Y): CN0, CN1, CN2 and mosaic gain/loss relative to the region's median LRR. -* LRR ~ (1−ε)·N(µ_state, σ) + ε·U(−5, 3); σ = 1.4826·MAD of the sample's autosomal LRR (floor 0.08); CN0 σ = 2. +* LRR ~ (1−ε)·N(µ_state, σ) + ε·U(−5, 3); σ = 1.4826·MAD of the sample's autosomal LRR (floor 0.08), scaled per + probe by the reference's LRR spread (CYTO-042); CN0 σ = 2. Mosaic means interpolate linearly between the constitutional means by the cell fraction. * BAF ~ (1−ε)·Σ_g w_g N(b_g, σ_b) + ε·U(0, 1); genotype weights from the probe's PFB p (binomial in the copy number; LOH: {0: 1−p, 1: p}); mosaic heterozygotes split into the two bands 0.5 ± d(f). @@ -106,7 +107,7 @@ def _norm(x: F64, mu: float, sd: float | F64) -> F64: @dataclass class Params: - lrr_sd: float + lrr_sd: float | F64 # scalar or per probe (CYTO-042) baf_sd: float lrr_offset: float = 0.0 # added to every state mean (haploid regions: the region's median LRR) switch_d: float = SWITCH_D diff --git a/src/simurg/cn/model.py b/src/simurg/cn/model.py index cc205a5..aaf8c3d 100644 --- a/src/simurg/cn/model.py +++ b/src/simurg/cn/model.py @@ -86,3 +86,18 @@ def reference_het_spread(baf: npt.NDArray[np.float32]) -> npt.NDArray[np.float32 mad = np.nanmedian(np.abs(het - med), axis=0) n = np.isfinite(het).sum(axis=0) return np.where(n >= 5, 1.4826 * mad, np.nan).astype(np.float32) + + +def reference_lrr_spread(lrr: npt.NDArray[np.float32], ref_lrr: npt.NDArray[np.float32]) -> npt.NDArray[np.float32]: + """Per-probe robust SD (1.4826·MAD) of LRR across reference samples after removing each sample's median and the + probe's reference median [OURS; per-probe noise as in PennCNV's per-marker SD, LIT: Wang et al. 2007]; NaN with + fewer than 5 finite values.""" + import warnings + + x = np.where(np.isfinite(lrr), lrr, np.nan).astype(np.float32) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", RuntimeWarning) + x = x - np.nanmedian(x, axis=1, keepdims=True) - ref_lrr[None, :] + mad = np.nanmedian(np.abs(x - np.nanmedian(x, axis=0)), axis=0) + n = np.isfinite(x).sum(axis=0) + return np.where(n >= 5, 1.4826 * mad, np.nan).astype(np.float32) diff --git a/src/simurg/pipeline/cn.py b/src/simurg/pipeline/cn.py index c4a7ea4..64ca663 100644 --- a/src/simurg/pipeline/cn.py +++ b/src/simurg/pipeline/cn.py @@ -47,6 +47,7 @@ class CytoModel: ref_lrr_m: npt.NDArray[np.float32] | None = None # same, X/Y from male reference samples ref_het: npt.NDArray[np.float32] | None = None # per-probe median heterozygous BAF ref_het_sd: npt.NDArray[np.float32] | None = None # per-probe robust SD of heterozygous BAF + ref_lrr_sd: npt.NDArray[np.float32] | None = None # per-probe robust SD of LRR across the reference def run_build_model( @@ -120,10 +121,11 @@ def run_build_model( ref_by_sex[sex] = col n_f, n_m = sexes.count("F"), sexes.count("M") rep.info(f"batch reference sexes: {n_f} F, {n_m} M, {len(sexes) - n_f - n_m} U") - from simurg.cn.model import reference_het_spread + from simurg.cn.model import reference_het_spread, reference_lrr_spread ref_cols = {"ref_lrr": ref_lrr, "ref_lrr_f": ref_by_sex["F"], "ref_lrr_m": ref_by_sex["M"], - "ref_het_baf": ref_het, "ref_het_sd": reference_het_spread(baf), "ref_n_het": n_het} # fmt: skip + "ref_het_baf": ref_het, "ref_het_sd": reference_het_spread(baf), "ref_n_het": n_het, + "ref_lrr_sd": reference_lrr_spread(lrr, ref_lrr)} # fmt: skip n_ref = len(paths) rep.info(f"batch reference from {n_ref} GTC(s)") header = { @@ -177,7 +179,7 @@ def read_model(path: Path, manifest: Manifest) -> CytoModel: assert pfb is not None return CytoModel(header, table.column("name").to_pylist(), gc.astype(np.float32), pfb, ref_lrr=col("ref_lrr"), ref_lrr_f=col("ref_lrr_f"), ref_lrr_m=col("ref_lrr_m"), ref_het=col("ref_het_baf"), - ref_het_sd=col("ref_het_sd")) # fmt: skip + ref_het_sd=col("ref_het_sd"), ref_lrr_sd=col("ref_lrr_sd")) # fmt: skip # ================================================================================================ cn call @@ -208,6 +210,7 @@ class CnCallOptions: bgzip: bool = True mosaic_prior: float = 0.004 # [OURS] tuning knobs (hmm.MOSAIC_PRIOR, hmm.SWITCH_D) switch_d: float = 5e8 + probe_lrr_sd: bool = True # CYTO-042: per-probe LRR noise from the reference threads: int = 0 force: bool = False argv: list[str] = field(default_factory=list) @@ -288,6 +291,22 @@ def _resolve_sex(si: SampleInput, chroms: npt.NDArray[np.str_], lrr: F64, geno: return "U", "unresolved", ["sex_unresolved"] +def _lrr_noise( + seg_sd: float, auto: npt.NDArray[np.bool_], ref_lrr_sd: npt.NDArray[np.float32] | None +) -> npt.NDArray[np.float64]: + """CYTO-042: per-probe emission SD = the sample's LRR SD × the probe's reference LRR spread relative to the + autosomal median spread, clipped to [0.5, 3]; the sample SD without a reference spread.""" + n = len(auto) + if ref_lrr_sd is None: + return np.full(n, seg_sd) + h = np.asarray(ref_lrr_sd, dtype=np.float64) + ok = np.isfinite(h) & (h > 0) + med = float(np.median(h[ok & auto])) if (ok & auto).any() else 0.0 + if med <= 0: + return np.full(n, seg_sd) + return seg_sd * np.where(ok, np.clip(h / med, 0.5, 3.0), 1.0) + + def _baf_noise( baf: npt.NDArray[np.float64], geno: npt.NDArray[np.uint8], auto: npt.NDArray[np.bool_], ref_het_sd: npt.NDArray[np.float32] | None, baf_sd: float, @@ -364,6 +383,7 @@ def call_sample(si: SampleInput, ctx: dict[str, Any], opts: CnCallOptions) -> di lrr_s[sel] = moving_median(lrr[sel], opts.smoothing) # emission SD from the unsmoothed LRR: smoothed values are autocorrelated, their own SD would overstate evidence seg_sd = _robust_sd(lrr[auto], 0.05) + lrr_sd_probe = _lrr_noise(seg_sd, auto, model.ref_lrr_sd if opts.probe_lrr_sd else None) # CYTO-042 b_het = baf[auto & (baf > 0.3) & (baf < 0.7)] baf_sd = _robust_sd(b_het, 0.02) sigma = _baf_noise(baf, geno, auto, model.ref_het_sd, baf_sd) # CYTO-041 @@ -405,7 +425,7 @@ def call_sample(si: SampleInput, ctx: dict[str, Any], opts: CnCallOptions) -> di offset = float(np.nanmedian(x_raw)) if haploid else 0.0 if not math.isfinite(offset): offset = 0.0 - prm = hmm.Params(seg_sd, baf_sd, lrr_offset=offset, switch_d=opts.switch_d) + prm = hmm.Params(lrr_sd_probe[idx], baf_sd, lrr_offset=offset, switch_d=opts.switch_d) het_sd = model.ref_het_sd[idx].astype(np.float64) if model.ref_het_sd is not None else None e = hmm.emissions(states, x_seg, np.full(len(b), np.nan) if haploid else b, pfb[idx], prm, het_sd) sw = hmm.switch_probs(p, opts.switch_d) diff --git a/tests/unit/test_cn.py b/tests/unit/test_cn.py index 7e0494e..e60cd3d 100644 --- a/tests/unit/test_cn.py +++ b/tests/unit/test_cn.py @@ -249,3 +249,18 @@ def test_mf_refinements() -> None: # CYTO-018 / CYTO-041 mixed = hom.copy() mixed[:8] = 2 assert not complete_loss(mixed, np.full(30, 0.3)) + + +def test_probe_lrr_noise() -> None: # CYTO-042 + from simurg.cn.model import reference_lrr_spread + from simurg.pipeline.cn import _lrr_noise + + rng = np.random.default_rng(5) + sd = np.array([0.1] * 50 + [0.3] * 50) + lrr = (rng.normal(0, 1, (40, 100)) * sd + rng.normal(0, 0.2, (40, 1))).astype(np.float32) # + sample offsets + spread = reference_lrr_spread(lrr, np.zeros(100, dtype=np.float32)) + assert abs(float(np.median(spread[:50])) - 0.1) < 0.02 and abs(float(np.median(spread[50:])) - 0.3) < 0.05 + auto = np.ones(100, dtype=np.bool_) + s = _lrr_noise(0.2, auto, spread) + assert np.median(s[50:]) / np.median(s[:50]) > 2 + assert np.all(_lrr_noise(0.2, auto, None) == 0.2)