From 6aebb31991f830c6d5a3d0f1a788993cfb3e49e6 Mon Sep 17 00:00:00 2001 From: Kemal Yaylali Date: Tue, 29 Sep 2026 00:36:44 +0100 Subject: [PATCH] methylation: combine_neg option; oracle compares both background variants and logs sesame's negative controls --- scripts/v2/methyl_oracle_compare.py | 35 ++++++++++++++++------------- scripts/v2/sesame_oracle.R | 3 +++ src/simurg/signal/methyl.py | 16 +++++++------ 3 files changed, 32 insertions(+), 22 deletions(-) diff --git a/scripts/v2/methyl_oracle_compare.py b/scripts/v2/methyl_oracle_compare.py index 28d616d..9dddd6d 100644 --- a/scripts/v2/methyl_oracle_compare.py +++ b/scripts/v2/methyl_oracle_compare.py @@ -38,9 +38,6 @@ def main() -> int: for prefix in a.prefixes: # min_beads=0: the oracle runs without bead-count masking (see sesame_oracle.R) s = sig.assemble(read_idat(prefix + "_Grn.idat"), read_idat(prefix + "_Red.idat"), man, min_beads=0) - p = sig.detection_p(s) - n, fg, fr = sig.noob(s) - b = sig.beta_values(n) rows = list(csv.DictReader((a.sesame / (Path(prefix).name + ".tsv")).open(), delimiter="\t")) idx = np.array([pos[r["probe_id"]] for r in rows]) @@ -48,25 +45,33 @@ def main() -> int: return np.array([float(r[name]) if r[name] not in ("NA", "") else np.nan for r in rows]) mask_s = np.array([r["mask"] == "TRUE" for r in rows]) + variants = {} + for neg in (True, False): + p = sig.detection_p(s, combine_neg=neg) + n, fg, fr = sig.noob(s, combine_neg=neg) + variants[neg] = (p, n, fg, fr, sig.beta_values(n)) keep = ~mask_s & ~s.masked[idx] sp, sm, su, sb = col("pval"), col("m_noob"), col("u_noob"), col("beta_noob") - op, om, ou, ob = p[idx], n.m[idx], n.u[idx], b[idx] ana = man.analytical()[idx] & keep - res[Path(prefix).name] = { + out: dict[str, object] = { "probes": len(rows), "masked_sesame": int(mask_s.sum()), "masked_simurg": int(s.masked.sum()), - "mask_agreement": float((mask_s == s.masked[idx]).mean()), - "noob_fit": {"green": vars(fg), "red": vars(fr)}, - "pval_abs_diff": _q(np.abs(op - sp)[keep]), - "pval_identical_fraction": float((op[keep] == sp[keep]).mean()), - "detection_call_agreement_0.05": float(((op < 0.05) == (sp < 0.05))[keep].mean()), - "prop_passing": {"simurg": float((op[ana] < 0.05).mean()), "sesame": float((sp[ana] < 0.05).mean())}, - "noob_m_rel_diff": _q(np.abs(om - sm)[keep] / np.maximum(np.abs(sm[keep]), 1)), - "noob_u_rel_diff": _q(np.abs(ou - su)[keep] / np.maximum(np.abs(su[keep]), 1)), - "beta_abs_diff": _q(np.abs(ob - sb)[keep]), - "beta_corr": float(np.corrcoef(ob[keep & np.isfinite(sb)], sb[keep & np.isfinite(sb)])[0, 1]), } + for neg, (p, n, fg, fr, b) in variants.items(): + op, om, ou, ob = p[idx], n.m[idx], n.u[idx], b[idx] + out[f"combine_neg={neg}"] = { + "noob_fit": {"green": vars(fg), "red": vars(fr)}, + "pval_abs_diff": _q(np.abs(op - sp)[keep]), + "pval_identical_fraction": float((op[keep] == sp[keep]).mean()), + "detection_call_agreement_0.05": float(((op < 0.05) == (sp < 0.05))[keep].mean()), + "prop_passing": {"simurg": float((op[ana] < 0.05).mean()), "sesame": float((sp[ana] < 0.05).mean())}, + "noob_m_rel_diff": _q(np.abs(om - sm)[keep] / np.maximum(np.abs(sm[keep]), 1)), + "noob_u_rel_diff": _q(np.abs(ou - su)[keep] / np.maximum(np.abs(su[keep]), 1)), + "beta_abs_diff": _q(np.abs(ob - sb)[keep]), + "beta_corr": float(np.corrcoef(ob[keep & np.isfinite(sb)], sb[keep & np.isfinite(sb)])[0, 1]), + } + res[Path(prefix).name] = out a.out.write_text(json.dumps(res, indent=1)) print(json.dumps(res, indent=1)) return 0 diff --git a/scripts/v2/sesame_oracle.R b/scripts/v2/sesame_oracle.R index daec1af..08bc36c 100644 --- a/scripts/v2/sesame_oracle.R +++ b/scripts/v2/sesame_oracle.R @@ -32,6 +32,9 @@ for (prefix in args[-(1:3)]) { sesameDataCache("idatSignature") read() }) + neg <- tryCatch(negControls(sdf), error = function(e) NULL) + cat("negative controls found by sesame:", if (is.null(neg)) "error" else nrow(neg), "\n") + print(head(sdf[grepl("^ctl", sdf$Probe_ID), ], 3)) pv <- pOOBAH(sdf, return.pval = TRUE) sn <- noob(sdf) ib <- inband(sn) diff --git a/src/simurg/signal/methyl.py b/src/simurg/signal/methyl.py index 0dae03e..a28c195 100644 --- a/src/simurg/signal/methyl.py +++ b/src/simurg/signal/methyl.py @@ -88,9 +88,9 @@ def _finite(a: F64) -> F64: return a[np.isfinite(a)] -def background(s: MethylSignals) -> tuple[F64, F64]: - """Out-of-band background per channel (green, red): Infinium I OOB signal + negative controls.""" - neg = np.array([t.casefold() == "negative" for t in s.ctrl_type], dtype=np.bool_) +def background(s: MethylSignals, combine_neg: bool = True) -> tuple[F64, F64]: + """Out-of-band background per channel (green, red): Infinium I OOB signal (+ negative controls).""" + neg = np.array([combine_neg and t.casefold() == "negative" for t in s.ctrl_type], dtype=np.bool_) red_i, grn_i = s.col == "R", s.col == "G" oob_g = np.concatenate([s.oob_m[red_i], s.oob_u[red_i], s.ctrl_green[neg]]) oob_r = np.concatenate([s.oob_m[grn_i], s.oob_u[grn_i], s.ctrl_red[neg]]) @@ -118,9 +118,9 @@ def _ecdf(ref: F64): # type: ignore[no-untyped-def] return f -def detection_p(s: MethylSignals) -> F64: +def detection_p(s: MethylSignals, combine_neg: bool = True) -> F64: """pOOBAH detection p-value per probe; 1 where masked.""" - bg_g, bg_r = background(s) + bg_g, bg_r = background(s, combine_neg) fg, fr = _ecdf(bg_g), _ecdf(bg_r) red_i, t2 = s.col == "R", s.col == "" m, u = np.nan_to_num(s.m, nan=0.0), np.nan_to_num(s.u, nan=0.0) @@ -194,9 +194,11 @@ def noob_fit(ib: F64, oob: F64) -> NoobFit: return NoobFit(mu, s, max(huber(ib)[0] - mu, 10.0)) -def noob(s: MethylSignals, offset: float = NOOB_OFFSET) -> tuple[MethylSignals, NoobFit, NoobFit]: +def noob( + s: MethylSignals, offset: float = NOOB_OFFSET, combine_neg: bool = True +) -> tuple[MethylSignals, NoobFit, NoobFit]: """Background-corrected copy of ``s`` (probe and control signals); returns the green and red fits.""" - bg_g, bg_r = background(s) + bg_g, bg_r = background(s, combine_neg) ib_g, ib_r = in_band(s) fit_g, fit_r = noob_fit(ib_g, bg_g), noob_fit(ib_r, bg_r)