methylation: combine_neg option; oracle compares both background variants and logs sesame's negative controls

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

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)

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)