CN: per-probe LRR noise from the batch reference (CYTO-042), --probe-lrr-sd/--no-probe-lrr-sd

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 |

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))

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

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)

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)

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)