v2 methylation control QC (WIP): simurg methyl qc|convert-manifest (MET-001..018)
Methylation CSV manifest reader (MET-002/003), raw M/U assembly with bead masking (MET-005), pOOBAH detection p-values, noob background correction, linear/nonlinear dye bias, beta and M-values (MET-010..012), the 21 methylation control metrics (MET-006, bs-conv red-bg set option), per-sample ctrl/cgs tables, summary tsv/xlsx/html with the detection histogram and a beta PCA, pcs.tsv.gz, controls_raw/controls_qc.tsv; methylation datasets load in qc report (MET-014); qc controls --array-type methylation; minfi-style Sentrix_ID/Sentrix_Position sheet keys. Vendor EPIC manifests carry a 5th field on some control lines: >= 4 fields accepted. Runs on the public minfiDataEPIC set (3 x NA12878). Co-Authored-By: Claude Opus 5.5 <[email protected]> Claude-Session: https://claude.ai/code/session_01814GKVzZVL3FuF3DeiZXhv
This commit is contained in:
kemalandClaude Opus 5.5 committed
1 parent
7c0e348033
commit
c1c2c30c87
13 files changed
+1292
-21
No files matched your search
@@ -666,6 +666,72 @@ def qc_controls_cmd(
|
||||
run_guarded(rep, lambda: run_qc_controls(opts, rep))
|
||||
|
||||
|
||||
# --------------------------------------------------------------------------------------------------------
|
||||
methyl_app = typer.Typer(help="Methylation arrays: control QC, detection, beta values (MET-001).")
|
||||
app.add_typer(methyl_app, name="methyl")
|
||||
|
||||
|
||||
@methyl_app.command("qc")
|
||||
def methyl_qc_cmd(
|
||||
manifest: Annotated[Path, typer.Option("--manifest", help="Methylation CSV manifest (required).")],
|
||||
out: Annotated[Path, typer.Option("--out", help="Output folder (required).")],
|
||||
idats: Annotated[Path | None, typer.Option("--idats", help="Folder searched recursively for IDAT pairs.")] = None,
|
||||
sheet: Annotated[Path | None, typer.Option("--sheet", help="Sample sheet (CSV).")] = None,
|
||||
config: Annotated[Path | None, typer.Option("--config", help="QC threshold YAML.")] = None,
|
||||
profile: Annotated[str, typer.Option("--profile", help="QC profile name, 'auto' or 'none'.")] = "auto",
|
||||
threshold: Annotated[list[str] | None, typer.Option("--threshold", help="KEY=VALUE (repeatable).")] = None,
|
||||
bg_offset: Annotated[float | None, typer.Option("--bg-offset", help="Background offset x (3000).")] = None,
|
||||
bs_set: Annotated[str, typer.Option("--bs-conv-red-bg-set", help="U3-U5 (default) or U4-U6.")] = "U3-U5",
|
||||
pval: Annotated[float, typer.Option("--pval-threshold", help="Detection p threshold.")] = 0.05,
|
||||
beta_offset: Annotated[float, typer.Option("--beta-offset", help="alpha in M/(M+U+alpha).")] = 0.0,
|
||||
dye_bias: Annotated[str, typer.Option("--dye-bias", help="linear|nonlinear")] = "linear",
|
||||
min_prop: Annotated[
|
||||
float | None, typer.Option("--min-prop-passing", help="Fail samples below this detection rate.")
|
||||
] = None,
|
||||
pca_full: Annotated[int, typer.Option("--pca-max-samples-full", min=2)] = 100,
|
||||
pca_probes: Annotated[int, typer.Option("--pca-probes", min=2)] = 10_000,
|
||||
threads: Annotated[int, typer.Option("--threads", min=0, max=1024, help="Worker processes (0 = auto).")] = 0,
|
||||
force: Annotated[bool, typer.Option("--force", help="Replace existing outputs.")] = False,
|
||||
log_json: LogJson = False,
|
||||
debug: Debug = False,
|
||||
quiet: Quiet = False,
|
||||
) -> None:
|
||||
"""Methylation control QC: 21 control metrics, detection p-values, noob/dye-bias corrected β (MET-001)."""
|
||||
from simurg.pipeline.methyl import MethylQcOptions, run_methyl_qc
|
||||
|
||||
if dye_bias not in ("linear", "nonlinear"):
|
||||
raise typer.BadParameter("--dye-bias must be linear or nonlinear")
|
||||
if bs_set not in ("U3-U5", "U4-U6"):
|
||||
raise typer.BadParameter("--bs-conv-red-bg-set must be U3-U5 or U4-U6")
|
||||
if not 0 < pval <= 1:
|
||||
raise typer.BadParameter("--pval-threshold must be in (0, 1]")
|
||||
rep = _reporter("methyl qc", out, log_json, debug, quiet)
|
||||
opts = MethylQcOptions(
|
||||
out=out, manifest=manifest, idats=idats, sheet=sheet, config=config,
|
||||
profile=None if profile == "none" else profile, thresholds=threshold or [], bg_offset=bg_offset,
|
||||
bs_conv_red_bg_set=bs_set, pval_threshold=pval, beta_offset=beta_offset, dye_bias=dye_bias,
|
||||
min_prop_passing=min_prop, pca_max_samples_full=pca_full, pca_probes=pca_probes, threads=threads,
|
||||
force=force, argv=sys.argv[1:],
|
||||
) # fmt: skip
|
||||
run_guarded(rep, lambda: run_methyl_qc(opts, rep))
|
||||
|
||||
|
||||
@methyl_app.command("convert-manifest")
|
||||
def methyl_convert_cmd(
|
||||
manifest: Annotated[Path, typer.Option("--manifest", help="Methylation CSV manifest (required).")],
|
||||
out: Annotated[Path, typer.Option("--out", help="Output folder (required).")],
|
||||
force: Annotated[bool, typer.Option("--force", help="Replace existing outputs.")] = False,
|
||||
log_json: LogJson = False,
|
||||
debug: Debug = False,
|
||||
quiet: Quiet = False,
|
||||
) -> None:
|
||||
"""Methylation manifest -> probes.tsv and controls.tsv (MET-003)."""
|
||||
from simurg.pipeline.methyl import run_convert_manifest
|
||||
|
||||
rep = _reporter("methyl convert-manifest", out, log_json, debug, quiet)
|
||||
run_guarded(rep, lambda: run_convert_manifest(manifest, out, force, rep))
|
||||
|
||||
|
||||
def _qc_cohort(
|
||||
command: str,
|
||||
out: Path,
|
||||
|
||||
@@ -116,7 +116,6 @@ CATALOGUE: dict[str, MessageDef] = dict(
|
||||
_m("W_CTRL_MISSING", "warning", "Controls needed by a metric are absent; metric not evaluable", "QC-007"),
|
||||
_m("W_CTRL_ZERO_DIVISOR", "warning", "Control metric has a zero denominator (NaN)", "QC-008"),
|
||||
_m("W_CTRL_DEFS_DIFFER", "warning", "BPM and CSV control definitions differ; the BPM wins", "QC-001"),
|
||||
_m("E_QC_METHYLATION", "error", "Methylation control metrics are not available yet (v2)", "QC-006"),
|
||||
_m("E_QC_NO_DATASETS", "error", "No QC dataset folder found", "QC-041"),
|
||||
_m("E_QC_DATASET_INCOMPLETE", "error", "Dataset folder lacks a required file", "QC-041"),
|
||||
_m("E_QC_LABEL_COUNT", "error", "Number of labels does not match number of datasets", "QC-042"),
|
||||
@@ -136,6 +135,37 @@ CATALOGUE: dict[str, MessageDef] = dict(
|
||||
_m("W_QC_PEDIGREE", "warning", "Parent1/Parent2/Replicate value does not name a sample in the set", "IN-018"),
|
||||
_m("W_QC_NO_BAF", "warning", "GTC has no B-allele frequencies; contamination not estimated", "QC-064"),
|
||||
_m("W_QC_TOO_FEW", "warning", "Too few samples or sites for a cohort QC statistic", "QC-026"),
|
||||
# --- methylation control QC (MET-001 .. MET-018) -------------------------------------------------
|
||||
_m("E_METH_NO_ASSAY", "error", "Probe section not found", "MET-002"),
|
||||
_m("E_METH_NO_CONTROLS", "error", "Missing control probes in manifest", "MET-002"),
|
||||
_m("E_METH_MISSING_COLUMNS", "error", "Missing required columns in the methylation manifest", "MET-002"),
|
||||
_m(
|
||||
"E_METH_CONTROLS_FORMAT",
|
||||
"error",
|
||||
"Controls not formatted correctly. Must have 4 columns (Address,Type,Color_Channel,Name)",
|
||||
"MET-002",
|
||||
),
|
||||
_m("E_METH_NO_IDATS", "error", "No IDATs found", "MET-004"),
|
||||
_m(
|
||||
"W_METH_MISSING_SAMPLE",
|
||||
"warning",
|
||||
"Missing sample: IDATs absent; the sample fails, the run continues",
|
||||
"MET-015",
|
||||
),
|
||||
_m("E_METH_PREP_FAILED", "error", "Prep failed (sample)", "MET-010"),
|
||||
_m(
|
||||
"W_METH_MISSING_CONTROL_TYPES",
|
||||
"warning",
|
||||
"Control types missing from the manifest; metric set to NA",
|
||||
"MET-009",
|
||||
),
|
||||
_m(
|
||||
"W_METH_MISSING_CONTROL_NAMES",
|
||||
"warning",
|
||||
"Control names missing from the manifest; metric set to NA",
|
||||
"MET-009",
|
||||
),
|
||||
_m("W_METH_NO_FEATURES", "warning", "No features, skipping PCA plot", "MET-013"),
|
||||
]
|
||||
)
|
||||
|
||||
|
||||
@@ -0,0 +1,157 @@
|
||||
"""Methylation CSV manifest (EPIC v1/v2, MSA and custom designs) — spec MET-002, MET-003; plan §6.3.
|
||||
|
||||
Sections: ``[Heading]`` (key,value), ``[Assay]`` (header row + one row per probe; needs ``IlmnID``, ``AddressA_ID``,
|
||||
``AddressB_ID``, ``Color_Channel``) and ``[Controls]`` (``Address,Type,Color_Channel,Name``, no header row; the
|
||||
vendor files pad lines with empty fields). Infinium I probes have an ``AddressB_ID``; the colour channel is read
|
||||
from ``Color_Channel`` (``Red``/``Grn``), Infinium II probes leave both empty.
|
||||
[DOC: Illumina methylation manifest column descriptions; layout of the open-source toolkits' manifests (MET-003)]
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import csv
|
||||
from dataclasses import dataclass, field
|
||||
from pathlib import Path
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
from simurg.core.errors import FormatError
|
||||
from simurg.formats.csv_manifest import _open_text
|
||||
|
||||
REQUIRED = ("IlmnID", "AddressA_ID", "AddressB_ID", "Color_Channel")
|
||||
|
||||
|
||||
@dataclass
|
||||
class MethylControl:
|
||||
address: int
|
||||
type: str
|
||||
color_channel: str
|
||||
name: str
|
||||
|
||||
@property
|
||||
def probe_id(self) -> str:
|
||||
return f"ctl_{self.address}"
|
||||
|
||||
|
||||
@dataclass
|
||||
class MethylManifest:
|
||||
path: str
|
||||
heading: dict[str, str]
|
||||
probe_ids: list[str]
|
||||
u: npt.NDArray[np.int64] # AddressA_ID
|
||||
m: npt.NDArray[np.int64] # AddressB_ID; 0 for Infinium II
|
||||
col: npt.NDArray[np.str_] # "R" / "G" for Infinium I, "" for Infinium II
|
||||
controls: list[MethylControl] = field(default_factory=list)
|
||||
|
||||
@property
|
||||
def num_probes(self) -> int:
|
||||
return len(self.probe_ids)
|
||||
|
||||
@property
|
||||
def type_i(self) -> npt.NDArray[np.bool_]:
|
||||
return np.asarray(self.m != 0, dtype=np.bool_)
|
||||
|
||||
def analytical(self) -> npt.NDArray[np.bool_]:
|
||||
"""CG and SNP probes (MET-011): IDs starting with ``cg`` or ``rs``."""
|
||||
return np.array([p[:2] in ("cg", "rs") for p in self.probe_ids], dtype=np.bool_)
|
||||
|
||||
|
||||
def _address(s: str) -> int:
|
||||
s = s.strip()
|
||||
return int(s) if s else 0
|
||||
|
||||
|
||||
def _channel(s: str, *, infinium_i: bool) -> str:
|
||||
c = s.strip().casefold()
|
||||
if not infinium_i:
|
||||
return ""
|
||||
if c in ("red", "r"):
|
||||
return "R"
|
||||
if c in ("grn", "green", "g"):
|
||||
return "G"
|
||||
return "?"
|
||||
|
||||
|
||||
def read_methyl_manifest(path: str | Path) -> MethylManifest:
|
||||
"""Parse a methylation CSV manifest (plain or gzipped). Raises the MET-002 errors."""
|
||||
p = Path(path)
|
||||
spath = str(p)
|
||||
heading: dict[str, str] = {}
|
||||
header: list[str] | None = None
|
||||
idx: dict[str, int] = {}
|
||||
ids: list[str] = []
|
||||
u: list[int] = []
|
||||
m: list[int] = []
|
||||
col: list[str] = []
|
||||
controls: list[MethylControl] = []
|
||||
section = None
|
||||
saw_assay = saw_controls = False
|
||||
with _open_text(p) as fh:
|
||||
for row in csv.reader(fh):
|
||||
if not row or not any(c.strip() for c in row):
|
||||
continue
|
||||
first = row[0].strip()
|
||||
if first.startswith("[") and first.endswith("]"):
|
||||
section = first[1:-1].strip().casefold()
|
||||
saw_assay |= section == "assay"
|
||||
saw_controls |= section == "controls"
|
||||
continue
|
||||
if section == "heading":
|
||||
heading[first] = row[1].strip() if len(row) > 1 else ""
|
||||
elif section == "assay":
|
||||
if header is None:
|
||||
header = [c.strip() for c in row]
|
||||
idx = {c: i for i, c in enumerate(header)}
|
||||
miss = [c for c in REQUIRED if c not in idx]
|
||||
if miss:
|
||||
msg = f"{spath}: Missing required columns: {', '.join(miss)}"
|
||||
raise FormatError("E_METH_MISSING_COLUMNS", msg, path=spath)
|
||||
iu, im, ic, ii = idx["AddressA_ID"], idx["AddressB_ID"], idx["Color_Channel"], idx["IlmnID"]
|
||||
continue
|
||||
ids.append(row[ii].strip())
|
||||
a, b = _address(row[iu]), _address(row[im]) if im < len(row) else 0
|
||||
u.append(a)
|
||||
m.append(b)
|
||||
col.append(_channel(row[ic] if ic < len(row) else "", infinium_i=b != 0))
|
||||
elif section == "controls":
|
||||
fields = list(row)
|
||||
while fields and not fields[-1].strip():
|
||||
fields.pop()
|
||||
if len(fields) < 4: # vendor EPIC files add a 5th field ("AVG") to some NEGATIVE lines (K-026)
|
||||
raise FormatError(
|
||||
"E_METH_CONTROLS_FORMAT",
|
||||
f"{spath}: Controls not formatted correctly. Must have 4 columns "
|
||||
f"(Address,Type,Color_Channel,Name): {','.join(fields)[:80]!r}",
|
||||
path=spath,
|
||||
)
|
||||
addr, typ, colour, name = (f.strip() for f in fields[:4])
|
||||
controls.append(MethylControl(_address(addr), typ, colour, name))
|
||||
if not saw_assay or header is None:
|
||||
raise FormatError("E_METH_NO_ASSAY", f"{spath}: Probe section not found", path=spath)
|
||||
if not saw_controls or not controls:
|
||||
raise FormatError("E_METH_NO_CONTROLS", f"{spath}: Missing control probes in manifest", path=spath)
|
||||
return MethylManifest(
|
||||
path=spath,
|
||||
heading=heading,
|
||||
probe_ids=ids,
|
||||
u=np.asarray(u, dtype=np.int64),
|
||||
m=np.asarray(m, dtype=np.int64),
|
||||
col=np.asarray(col, dtype="<U1"),
|
||||
controls=controls,
|
||||
)
|
||||
|
||||
|
||||
def write_probe_tables(man: MethylManifest, out: Path) -> tuple[Path, Path]:
|
||||
"""MET-003: ``probes.tsv`` (probe_id, u, m, col; ``NA`` for Infinium II) and ``controls.tsv``."""
|
||||
out.mkdir(parents=True, exist_ok=True)
|
||||
probes, ctl = out / "probes.tsv", out / "controls.tsv"
|
||||
with probes.open("w", encoding="utf-8", newline="\n") as fh:
|
||||
fh.write("probe_id\tu\tm\tcol\n")
|
||||
for pid, a, b, c in zip(man.probe_ids, man.u.tolist(), man.m.tolist(), man.col.tolist(), strict=True):
|
||||
fh.write(f"{pid}\t{a}\t{b if b else 'NA'}\t{c or 'NA'}\n")
|
||||
with ctl.open("w", encoding="utf-8", newline="\n") as fh:
|
||||
fh.write("address\ttype\tcolor_channel\tname\tprobe_id\n")
|
||||
for c in man.controls:
|
||||
fh.write(f"{c.address}\t{c.type}\t{c.color_channel}\t{c.name}\t{c.probe_id}\n")
|
||||
return probes, ctl
|
||||
@@ -15,7 +15,12 @@ from pathlib import Path
|
||||
from simurg.core.errors import SimurgError
|
||||
|
||||
MAX_LEN = 500
|
||||
KEY_PAIRS = (("SentrixBarcode_A", "SentrixPosition_A"), ("beadChipName", "sampleSectionName"))
|
||||
# the third pair is the methylation toolkits' convention (minfi/sesame sheets), accepted for MET-004 (K-026)
|
||||
KEY_PAIRS = (
|
||||
("SentrixBarcode_A", "SentrixPosition_A"),
|
||||
("beadChipName", "sampleSectionName"),
|
||||
("Sentrix_ID", "Sentrix_Position"),
|
||||
)
|
||||
RECOGNISED = (
|
||||
"Sample_ID",
|
||||
"Sample_Name",
|
||||
@@ -246,7 +251,8 @@ def parse_samplesheet(path: str | Path, *, on_duplicate_sample: str = "error") -
|
||||
raise _err(
|
||||
"E_SHEET_NO_KEY_COLUMNS",
|
||||
spath,
|
||||
"missing SentrixBarcode_A/SentrixPosition_A (or beadChipName/sampleSectionName) columns",
|
||||
"missing SentrixBarcode_A/SentrixPosition_A (or beadChipName/sampleSectionName, "
|
||||
"Sentrix_ID/Sentrix_Position) columns",
|
||||
)
|
||||
|
||||
recognised = {fold(c) for c in RECOGNISED} | {fold(key[0]), fold(key[1])}
|
||||
|
||||
@@ -0,0 +1,443 @@
|
||||
"""``simurg methyl qc`` / ``simurg methyl convert-manifest`` — methylation control QC (spec MET-001 … MET-018).
|
||||
|
||||
Per sample (streamed, worker processes): control table, detection p-values, noob + dye-bias corrected β and M-values.
|
||||
Per batch: ``methyl_qc_summary.tsv/.xlsx/.html``, ``pcs.tsv.gz``, ``probes.tsv``/``controls.tsv`` and the shared
|
||||
``controls_raw.tsv``/``controls_qc.tsv`` (so ``simurg qc report`` accepts the folder).
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import datetime as _dt
|
||||
import gzip
|
||||
import html
|
||||
import json
|
||||
import math
|
||||
import multiprocessing as mp
|
||||
import os
|
||||
import pickle
|
||||
import tempfile
|
||||
from concurrent.futures import ProcessPoolExecutor
|
||||
from dataclasses import dataclass, field
|
||||
from pathlib import Path
|
||||
from typing import Any
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
from simurg.core.atomic import check_outputs_absent
|
||||
from simurg.core.errors import ExitCode, SimurgError
|
||||
from simurg.core.log import Reporter
|
||||
from simurg.core.provenance import RunManifest
|
||||
from simurg.formats.idat import check_pair, read_idat
|
||||
from simurg.formats.methyl_manifest import MethylManifest, read_methyl_manifest, write_probe_tables
|
||||
from simurg.pipeline.qc import ControlSample, idat_metadata, write_control_tables, write_tsv
|
||||
from simurg.qc.controls import ControlBead, canonical_id, compute_control_metrics, metrics_for
|
||||
from simurg.qc.thresholds import Thresholds, parse_cli_pairs, resolve_thresholds
|
||||
from simurg.signal import methyl as sig
|
||||
|
||||
SUMMARY = "methyl_qc_summary"
|
||||
CONTROL_TYPES = ("STAINING", "EXTENSION", "HYBRIDIZATION", "TARGET REMOVAL", "NON-POLYMORPHIC", "SPECIFICITY I",
|
||||
"SPECIFICITY II", "BISULFITE CONVERSION I", "BISULFITE CONVERSION II", "NEGATIVE", "NORM_A",
|
||||
"NORM_C", "NORM_G", "NORM_T") # fmt: skip
|
||||
AVG_COLUMNS = ("avg_green_raw", "avg_red_raw", "avg_green_norm", "avg_red_norm")
|
||||
|
||||
|
||||
@dataclass
|
||||
class MethylQcOptions:
|
||||
out: Path
|
||||
manifest: Path
|
||||
idats: Path | None = None
|
||||
sheet: Path | None = None
|
||||
config: Path | None = None
|
||||
profile: str | None = "auto"
|
||||
thresholds: list[str] = field(default_factory=list)
|
||||
bg_offset: float | None = None
|
||||
bs_conv_red_bg_set: str = "U3-U5"
|
||||
pval_threshold: float = 0.05
|
||||
beta_offset: float = 0.0
|
||||
dye_bias: str = "linear"
|
||||
min_prop_passing: float | None = None
|
||||
pca_max_samples_full: int = 100
|
||||
pca_probes: int = 10_000
|
||||
seed: int = 42
|
||||
threads: int = 0
|
||||
force: bool = False
|
||||
on_duplicate_idat: str = "error"
|
||||
argv: list[str] = field(default_factory=list)
|
||||
|
||||
|
||||
@dataclass
|
||||
class _Job:
|
||||
sample_id: str
|
||||
grn: str | None
|
||||
red: str | None
|
||||
sheet_values: dict[str, str]
|
||||
out_dir: str
|
||||
|
||||
|
||||
_W: dict[str, Any] = {}
|
||||
|
||||
|
||||
def _init_from_file(path: str) -> None:
|
||||
with open(path, "rb") as fh:
|
||||
_W.update(pickle.load(fh))
|
||||
|
||||
|
||||
def _write_gz_tsv(path: Path, header: str, lines: list[str]) -> None:
|
||||
tmp = path.with_suffix(path.suffix + ".part")
|
||||
with gzip.open(tmp, "wt", encoding="utf-8", newline="\n", compresslevel=6) as fh:
|
||||
fh.write(header + "\n")
|
||||
fh.write("\n".join(lines))
|
||||
if lines:
|
||||
fh.write("\n")
|
||||
tmp.replace(path)
|
||||
|
||||
|
||||
def _g6(v: float) -> str:
|
||||
return "NA" if not math.isfinite(v) else f"{v:.6g}"
|
||||
|
||||
|
||||
def _process(job: _Job) -> dict[str, Any]:
|
||||
man: MethylManifest = _W["man"]
|
||||
opts: MethylQcOptions = _W["opts"]
|
||||
beads: list[ControlBead] = _W["beads"]
|
||||
row: dict[str, Any] = {"sample_id": job.sample_id, "status": "OK", "fail_reason": ""}
|
||||
if not job.grn or not job.red:
|
||||
row.update(status="FAILED", fail_reason="W_METH_MISSING_SAMPLE: IDATs absent")
|
||||
return row
|
||||
step = "read"
|
||||
try:
|
||||
grn, red = read_idat(job.grn), read_idat(job.red)
|
||||
check_pair(grn, red)
|
||||
row["meta"] = idat_metadata(grn, red, job.sheet_values)
|
||||
row["sentrix_id"], row["sentrix_position"] = grn.barcode or red.barcode, grn.position or red.position
|
||||
row["scan_time"] = red.scan_time_iso or grn.scan_time_iso
|
||||
step = "signal assembly"
|
||||
raw = sig.assemble(grn, red, man)
|
||||
g, r = raw.ctrl_green, raw.ctrl_red
|
||||
qc = compute_control_metrics(beads, g, r, "methylation", _W["bg_offset"], opts.bs_conv_red_bg_set)
|
||||
row["control"] = ControlSample(job.sample_id, row["meta"], g, r, qc)
|
||||
step = "detection"
|
||||
pval = sig.detection_p(raw)
|
||||
step = "noob"
|
||||
corr, _, _ = sig.noob(raw)
|
||||
step = "dye bias"
|
||||
corr = sig.dye_bias_nonlinear(corr) if opts.dye_bias == "nonlinear" else sig.dye_bias_linear(corr)
|
||||
step = "beta"
|
||||
beta = sig.beta_values(corr, opts.beta_offset)
|
||||
mval = sig.m_values(beta)
|
||||
ana = _W["analytical"]
|
||||
row["n_probes"] = int(ana.sum())
|
||||
row["n_pass_detection"] = int((pval[ana] < opts.pval_threshold).sum())
|
||||
row["prop_probes_passing"] = row["n_pass_detection"] / row["n_probes"] if row["n_probes"] else math.nan
|
||||
row["avg_green_raw"], row["avg_red_raw"] = sig.average_signals(raw, ana)
|
||||
row["avg_green_norm"], row["avg_red_norm"] = sig.average_signals(corr, ana)
|
||||
row["pca"] = beta[_W["pca_idx"]].astype(np.float32)
|
||||
step = "write"
|
||||
d = Path(job.out_dir)
|
||||
d.mkdir(parents=True, exist_ok=True)
|
||||
_write_gz_tsv(
|
||||
d / f"{job.sample_id}_ctrl.tsv.gz",
|
||||
"address\ttype\tcolor_channel\tname\tprobe_id\tgreen\tred",
|
||||
[
|
||||
f"{c.address}\t{c.type}\t{c.color_channel}\t{c.name}\t{c.probe_id}\t{_g6(gv)}\t{_g6(rv)}"
|
||||
for c, gv, rv in zip(man.controls, g.tolist(), r.tolist(), strict=True)
|
||||
],
|
||||
)
|
||||
_write_gz_tsv(
|
||||
d / f"{job.sample_id}_cgs.tsv.gz",
|
||||
"probe_id\tbeta\tm_value\tdetection_p",
|
||||
[
|
||||
f"{pid}\t{_g6(b)}\t{_g6(mv)}\t{_g6(pv)}"
|
||||
for pid, b, mv, pv in zip(man.probe_ids, beta.tolist(), mval.tolist(), pval.tolist(), strict=True)
|
||||
],
|
||||
)
|
||||
except Exception as e: # MET-010: this sample fails, the others continue
|
||||
code = e.code if isinstance(e, SimurgError) else "E_METH_PREP_FAILED"
|
||||
row.update(status="FAILED", fail_reason=f"{code}: Prep failed ({step}): {e}")
|
||||
row.pop("control", None)
|
||||
return row
|
||||
|
||||
|
||||
def _pca(mat: npt.NDArray[np.float32]) -> npt.NDArray[np.float64] | None:
|
||||
"""First two principal components (samples × 2) of a samples × probes matrix; None when not computable."""
|
||||
ok = np.isfinite(mat).all(axis=0)
|
||||
x = mat[:, ok].astype(np.float64)
|
||||
if x.shape[0] < 2 or x.shape[1] < 2:
|
||||
return None
|
||||
x -= x.mean(axis=0)
|
||||
gram = x @ x.T
|
||||
w, v = np.linalg.eigh(gram)
|
||||
order = np.argsort(w)[::-1][:2]
|
||||
pcs = v[:, order] * np.sqrt(np.maximum(w[order], 0))
|
||||
for j in range(pcs.shape[1]): # deterministic sign: largest |loading| positive
|
||||
k = int(np.argmax(np.abs(pcs[:, j])))
|
||||
if pcs[k, j] < 0:
|
||||
pcs[:, j] = -pcs[:, j]
|
||||
if pcs.shape[1] < 2:
|
||||
pcs = np.column_stack([pcs, np.zeros(len(pcs))])
|
||||
return pcs
|
||||
|
||||
|
||||
def _thresholds(opts: MethylQcOptions, man: MethylManifest) -> Thresholds:
|
||||
cli: dict[str, Any] = parse_cli_pairs(opts.thresholds)
|
||||
if opts.bg_offset is not None:
|
||||
cli["bgOffset"] = opts.bg_offset
|
||||
return resolve_thresholds(
|
||||
config=opts.config,
|
||||
profile=opts.profile,
|
||||
cli=cli,
|
||||
array_type="methylation",
|
||||
manifest_name=man.heading.get("Descriptor File Name", Path(man.path).name),
|
||||
)
|
||||
|
||||
|
||||
def run_methyl_qc(opts: MethylQcOptions, rep: Reporter) -> int:
|
||||
from simurg.pipeline.inputs import check_paths, resolve_sample_set, sheet_values
|
||||
|
||||
rm = RunManifest(command="methyl qc", argv=opts.argv)
|
||||
rep.run_id = rm.run_id
|
||||
check_paths([(opts.manifest, "methylation manifest"), (opts.idats, "IDAT folder"), (opts.sheet, "sheet")])
|
||||
try:
|
||||
sset = resolve_sample_set(opts.idats, opts.sheet, opts.on_duplicate_idat, rep, rm)
|
||||
except SimurgError as e:
|
||||
if e.code == "E_NO_INPUT":
|
||||
raise SimurgError("E_METH_NO_IDATS", f"No IDATs found ({e.message})") from None
|
||||
raise
|
||||
outputs = [opts.out / f"{SUMMARY}.{x}" for x in ("tsv", "xlsx", "html")] + [opts.out / "pcs.tsv.gz"]
|
||||
check_outputs_absent(outputs, force=opts.force)
|
||||
rep.info(f"reading methylation manifest {opts.manifest}")
|
||||
man = read_methyl_manifest(opts.manifest)
|
||||
opts.out.mkdir(parents=True, exist_ok=True)
|
||||
for p in write_probe_tables(man, opts.out):
|
||||
rm.add_output(p, "simurg.methyl_manifest/1")
|
||||
th = _thresholds(opts, man)
|
||||
for rec in th.as_records():
|
||||
rm.set_option(f"qc.{rec['key']}", rec["value"], rec["source"])
|
||||
beads = [ControlBead(c.type, c.color_channel, c.name, c.address, canonical_id(c.name)) for c in man.controls]
|
||||
types = {c.type.upper() for c in man.controls}
|
||||
if miss := [t for t in CONTROL_TYPES if t not in types]:
|
||||
rep.warn("W_METH_MISSING_CONTROL_TYPES", f"manifest lacks control type(s): {', '.join(miss)}")
|
||||
analytical = man.analytical()
|
||||
n = len(sset.samples)
|
||||
cand = np.flatnonzero(analytical)
|
||||
if n > opts.pca_max_samples_full and len(cand) > opts.pca_probes:
|
||||
cand = np.sort(np.random.default_rng(opts.seed).choice(cand, opts.pca_probes, replace=False))
|
||||
jobs = []
|
||||
for s in sset.samples:
|
||||
grn = str(s.pair.grn) if s.pair is not None and s.pair.grn is not None else None
|
||||
red = str(s.pair.red) if s.pair is not None and s.pair.red is not None else None
|
||||
if not (grn and red):
|
||||
rep.warn("W_METH_MISSING_SAMPLE", f"Missing sample: {s.sample_id}", sample_id=s.sample_id)
|
||||
jobs.append(_Job(s.sample_id, grn, red, sheet_values(s, sset.sheet), str(opts.out / "samples" / s.sample_id)))
|
||||
state = {"man": man, "opts": opts, "beads": beads, "bg_offset": th.bg_offset, "analytical": analytical,
|
||||
"pca_idx": cand} # fmt: skip
|
||||
workers = max(1, min(opts.threads or (os.cpu_count() or 1), len(jobs) or 1))
|
||||
rep.info(f"methylation QC of {len(jobs)} sample(s) with {workers} worker(s)")
|
||||
if workers == 1:
|
||||
_W.update(state)
|
||||
try:
|
||||
rows = [_process(j) for j in jobs]
|
||||
finally:
|
||||
_W.clear()
|
||||
else:
|
||||
fd, path = tempfile.mkstemp(prefix=".simurg-workers-", suffix=".pkl", dir=opts.out)
|
||||
try:
|
||||
with os.fdopen(fd, "wb") as fh:
|
||||
pickle.dump(state, fh, protocol=5)
|
||||
ctx = mp.get_context("spawn")
|
||||
with ProcessPoolExecutor(workers, mp_context=ctx, initializer=_init_from_file, initargs=(path,)) as ex:
|
||||
rows = list(ex.map(_process, jobs))
|
||||
finally:
|
||||
Path(path).unlink(missing_ok=True)
|
||||
|
||||
metrics = metrics_for("methylation")
|
||||
col_key = {m.column: m.key for m in metrics}
|
||||
missing_names: dict[str, list[str]] = {}
|
||||
for r in rows:
|
||||
if r["status"] != "OK":
|
||||
rep.error(
|
||||
r["fail_reason"].split(":", 1)[0], f"{r['sample_id']}: {r['fail_reason']}", sample_id=r["sample_id"]
|
||||
)
|
||||
continue
|
||||
q = r["control"].qc
|
||||
missing_names.update(q.missing)
|
||||
if not r["scan_time"]:
|
||||
rep.warn("W_IDAT_NO_SCAN_TIME", f"{r['sample_id']}: No scan time available", sample_id=r["sample_id"])
|
||||
fails = [m.column for m in metrics if m.column not in q.missing and th.get(m.key).fails(q.values[m.column])]
|
||||
if opts.min_prop_passing is not None and not r["prop_probes_passing"] >= opts.min_prop_passing:
|
||||
fails.append("prop_probes_passing")
|
||||
r["fail_codes"] = ";".join(fails)
|
||||
r["na_metrics"] = ";".join(q.na_metrics)
|
||||
r["pass_qc"] = 0 if fails else 1
|
||||
r.update({m.column: q.values[m.column] for m in metrics})
|
||||
if missing_names:
|
||||
ids = sorted({c for v in missing_names.values() for c in v})
|
||||
rep.warn("W_METH_MISSING_CONTROL_NAMES", f"manifest lacks control(s) {', '.join(ids)}; "
|
||||
f"not evaluable: {', '.join(sorted(missing_names))}") # fmt: skip
|
||||
|
||||
meta_cols: list[str] = []
|
||||
for j in jobs:
|
||||
for k in j.sheet_values:
|
||||
if k not in meta_cols:
|
||||
meta_cols.append(k)
|
||||
cols = ["sentrix_id", "sentrix_position", "sample_id", *meta_cols, *(m.column for m in metrics), *AVG_COLUMNS,
|
||||
"scan_time", "n_probes", "n_pass_detection", "prop_probes_passing", "pass_qc", "fail_codes",
|
||||
"na_metrics", "status", "fail_reason"] # fmt: skip
|
||||
out_rows = []
|
||||
for j, r in zip(jobs, rows, strict=True):
|
||||
out_rows.append({**j.sheet_values, **{c: r.get(c, "") for c in cols if c not in j.sheet_values}})
|
||||
summary_tsv = opts.out / f"{SUMMARY}.tsv"
|
||||
write_tsv(summary_tsv, cols, out_rows)
|
||||
rm.add_output(summary_tsv, "simurg.methyl_qc_summary/1")
|
||||
_write_xlsx(opts.out / f"{SUMMARY}.xlsx", cols, out_rows, th, col_key)
|
||||
|
||||
ok_rows = [r for r in rows if r["status"] == "OK"]
|
||||
samples = [r["control"] for r in ok_rows]
|
||||
if samples:
|
||||
for p in write_control_tables(opts.out, samples, beads, th.bg_offset, rep, array_type="methylation"):
|
||||
rm.add_output(p, "simurg.controls/1")
|
||||
pcs = _pca(np.vstack([r["pca"] for r in ok_rows])) if len(ok_rows) >= 2 else None
|
||||
group_col = next((k for k in meta_cols if k.casefold() == "sample_group"), None)
|
||||
groups = {j.sample_id: (j.sheet_values.get(group_col, "") if group_col else "") for j in jobs}
|
||||
pc_rows = []
|
||||
if pcs is None:
|
||||
rep.warn("W_METH_NO_FEATURES", "No features, skipping PCA plot")
|
||||
else:
|
||||
for r, (p1, p2) in zip(ok_rows, pcs.tolist(), strict=True):
|
||||
pc_rows.append({"sample_id": r["sample_id"], "pc1": p1, "pc2": p2,
|
||||
"sample_group": groups[r["sample_id"]] or "NA"}) # fmt: skip
|
||||
with gzip.open(opts.out / "pcs.tsv.gz", "wt", encoding="utf-8", newline="\n") as fh:
|
||||
fh.write("sample_id\tpc1\tpc2\tsample_group\n")
|
||||
for pr in pc_rows:
|
||||
fh.write(f"{pr['sample_id']}\t{pr['pc1']:.6g}\t{pr['pc2']:.6g}\t{pr['sample_group']}\n")
|
||||
_write_html(opts.out / f"{SUMMARY}.html", out_rows, pc_rows, th, metrics, man, len(cand))
|
||||
n_ok = len(ok_rows)
|
||||
rm.summary = {"samples": len(rows), "prepared": n_ok, "pass_qc": sum(r.get("pass_qc") == 1 for r in ok_rows)}
|
||||
rm.status = "completed"
|
||||
rm.write(opts.out / "run_manifest.json")
|
||||
rep.info(f"done: {n_ok}/{len(rows)} samples prepared, {rm.summary['pass_qc']} pass QC; summary {summary_tsv}")
|
||||
if n_ok == len(rows):
|
||||
return ExitCode.OK
|
||||
return ExitCode.PARTIAL if n_ok else ExitCode.ALL_FAILED
|
||||
|
||||
|
||||
def _write_xlsx(
|
||||
path: Path, cols: list[str], rows: list[dict[str, Any]], th: Thresholds, col_key: dict[str, str]
|
||||
) -> None:
|
||||
from openpyxl import Workbook
|
||||
from openpyxl.styles import Font, PatternFill
|
||||
|
||||
bad = PatternFill("solid", fgColor="F5C28E")
|
||||
wb = Workbook()
|
||||
ws = wb.active
|
||||
ws.title = "Samples"
|
||||
ws.append(cols)
|
||||
for c in ws[1]:
|
||||
c.font = Font(bold=True)
|
||||
for r in rows:
|
||||
ws.append([_xl(r.get(c)) for c in cols])
|
||||
cur = ws.max_row
|
||||
for j, c in enumerate(cols, start=1):
|
||||
v = r.get(c)
|
||||
na = c in str(r.get("na_metrics", "")).split(";")
|
||||
metric_fails = c in col_key and isinstance(v, float) and not na and th.get(col_key[c]).fails(v)
|
||||
if metric_fails or (c == "pass_qc" and v == 0):
|
||||
ws.cell(cur, j).fill = bad
|
||||
ws.freeze_panes = "D2"
|
||||
ws.auto_filter.ref = ws.dimensions
|
||||
t2 = wb.create_sheet("Thresholds")
|
||||
t2.append(["key", "threshold", "direction", "source"])
|
||||
for k in col_key.values():
|
||||
t = th.get(k)
|
||||
t2.append([k, t.value, t.direction, t.source])
|
||||
tmp = path.with_suffix(".xlsx.part")
|
||||
wb.save(tmp)
|
||||
tmp.replace(path)
|
||||
|
||||
|
||||
def _xl(v: Any) -> Any:
|
||||
if isinstance(v, float) and not math.isfinite(v):
|
||||
return "NaN" if math.isnan(v) else ("Inf" if v > 0 else "-Inf")
|
||||
return "" if v is None else v
|
||||
|
||||
|
||||
def _write_html(
|
||||
path: Path,
|
||||
rows: list[dict[str, Any]],
|
||||
pcs: list[dict[str, Any]],
|
||||
th: Thresholds,
|
||||
metrics: list[Any],
|
||||
man: MethylManifest,
|
||||
n_pca_probes: int,
|
||||
) -> None:
|
||||
"""MET-013: self-contained page — summary table, detection histogram, PCA of β."""
|
||||
from simurg.pipeline.qc_report import _plotly
|
||||
|
||||
ok = [r for r in rows if r.get("status") == "OK"]
|
||||
fin = [r for r in ok if math.isfinite(r["prop_probes_passing"])]
|
||||
hist = {
|
||||
"pass": [r["prop_probes_passing"] for r in fin if r.get("pass_qc") == 1],
|
||||
"fail": [r["prop_probes_passing"] for r in fin if r.get("pass_qc") != 1],
|
||||
}
|
||||
data = {
|
||||
"hist": hist,
|
||||
"pcs": pcs,
|
||||
"n_pca_probes": n_pca_probes,
|
||||
}
|
||||
head = ["sample_id", "pass_qc", "prop_probes_passing", "fail_codes", "na_metrics", "scan_time"]
|
||||
trs = "".join("<tr>" + "".join(f"<td>{html.escape(_cell(r.get(c)))}</td>" for c in head) + "</tr>" for r in rows)
|
||||
th_rows = "".join(
|
||||
f"<tr><td>{m.column}</td><td>{m.key}</td><td>{_cell(th.get(m.key).value)}</td></tr>" for m in metrics
|
||||
)
|
||||
generated = _dt.datetime.now(_dt.UTC).strftime("%Y-%m-%d %H:%M UTC")
|
||||
page = f"""<!doctype html><html lang="en"><head><meta charset="utf-8">
|
||||
<meta name="viewport" content="width=device-width, initial-scale=1"><title>Methylation QC</title>
|
||||
<style>body{{font:14px system-ui,sans-serif;margin:16px;color:#222;background:#fff}}table{{border-collapse:collapse}}
|
||||
td,th{{border:1px solid #ccc;padding:3px 6px;text-align:left}}.plots{{display:flex;flex-wrap:wrap;gap:16px}}
|
||||
.plot{{width:min(560px,100%);height:380px}}.note{{color:#666}}</style>
|
||||
<script>{_plotly()}</script></head><body>
|
||||
<h1>Methylation control QC</h1>
|
||||
<p class="note">{len(rows)} samples · manifest {html.escape(Path(man.path).name)}
|
||||
· profile {html.escape(th.profile or "none")} · generated {generated} · SimurgArray</p>
|
||||
<div class="plots"><div id="hist" class="plot"></div><div id="pca" class="plot"></div></div>
|
||||
<h2>Samples</h2><table><tr>{"".join(f"<th>{c}</th>" for c in head)}</tr>{trs}</table>
|
||||
<h2>Thresholds</h2><table><tr><th>metric</th><th>key</th><th>threshold</th></tr>{th_rows}</table>
|
||||
<script>const D={json.dumps(data, allow_nan=False, default=_json_num)};
|
||||
Plotly.newPlot("hist",[{{x:D.hist.pass,type:"histogram",name:"pass",marker:{{color:"#4c78a8"}}}},
|
||||
{{x:D.hist.fail,type:"histogram",name:"fail",marker:{{color:"#e45756"}}}}],
|
||||
{{title:"Probes passing detection",barmode:"overlay",xaxis:{{title:"proportion"}},yaxis:{{title:"samples"}}}});
|
||||
if(D.pcs.length){{const g=[...new Set(D.pcs.map(p=>p.sample_group))];
|
||||
Plotly.newPlot("pca",g.map(k=>({{x:D.pcs.filter(p=>p.sample_group===k).map(p=>p.pc1),
|
||||
y:D.pcs.filter(p=>p.sample_group===k).map(p=>p.pc2),text:D.pcs.filter(p=>p.sample_group===k).map(p=>p.sample_id),
|
||||
mode:"markers",type:"scatter",name:k}})),{{title:"PCA of beta ("+D.n_pca_probes+" candidate probes)",
|
||||
xaxis:{{title:"PC1"}},yaxis:{{title:"PC2"}}}});}}
|
||||
else{{document.getElementById("pca").textContent="No features, PCA skipped";}}
|
||||
</script></body></html>
|
||||
"""
|
||||
tmp = path.with_suffix(".html.part")
|
||||
tmp.write_text(page, encoding="utf-8")
|
||||
tmp.replace(path)
|
||||
|
||||
|
||||
def _cell(v: Any) -> str:
|
||||
if v is None:
|
||||
return ""
|
||||
if isinstance(v, float):
|
||||
return "NaN" if math.isnan(v) else f"{v:.4g}"
|
||||
return str(v)
|
||||
|
||||
|
||||
def _json_num(v: Any) -> Any:
|
||||
if isinstance(v, np.floating | np.integer):
|
||||
return v.item()
|
||||
raise TypeError(type(v))
|
||||
|
||||
|
||||
def run_convert_manifest(manifest: Path, out: Path, force: bool, rep: Reporter) -> int:
|
||||
"""``simurg methyl convert-manifest`` (MET-003)."""
|
||||
check_outputs_absent([out / "probes.tsv", out / "controls.tsv"], force=force)
|
||||
man = read_methyl_manifest(manifest)
|
||||
probes, ctl = write_probe_tables(man, out)
|
||||
rep.info(f"{man.num_probes} probes -> {probes}; {len(man.controls)} controls -> {ctl}")
|
||||
return ExitCode.OK
|
||||
@@ -36,6 +36,7 @@ from simurg.qc.controls import (
|
||||
ControlBead,
|
||||
ControlQC,
|
||||
beads_from_definitions,
|
||||
canonical_id,
|
||||
compare_definitions,
|
||||
compute_control_metrics,
|
||||
control_intensities,
|
||||
@@ -123,9 +124,10 @@ def control_sample(
|
||||
beads: Sequence[ControlBead],
|
||||
bg_offset: float,
|
||||
sheet_values: Mapping[str, str] | None = None,
|
||||
array_type: str = "genotyping",
|
||||
) -> ControlSample:
|
||||
g, r = control_intensities(grn, red, beads)
|
||||
qc = compute_control_metrics(beads, g, r, "genotyping", bg_offset)
|
||||
qc = compute_control_metrics(beads, g, r, array_type, bg_offset)
|
||||
return ControlSample(sample_id, idat_metadata(grn, red, sheet_values), g, r, qc)
|
||||
|
||||
|
||||
@@ -250,8 +252,8 @@ def run_qc_controls(opts: QcControlsOptions, rep: Reporter) -> int:
|
||||
|
||||
rm = RunManifest(command="qc controls", argv=opts.argv)
|
||||
rep.run_id = rm.run_id
|
||||
if opts.array_type == "methylation":
|
||||
raise SimurgError("E_QC_METHYLATION", "methylation control metrics are planned for v2 (QC-006)")
|
||||
if opts.array_type == "methylation" and opts.csv is None:
|
||||
raise SimurgError("E_NO_INPUT", "methylation control QC needs the methylation CSV manifest (--csv)")
|
||||
if opts.bpm is None and opts.csv is None:
|
||||
raise SimurgError("E_NO_INPUT", "give --bpm or --csv for the control definitions (QC-001)")
|
||||
if opts.bg_offset < 0:
|
||||
@@ -261,8 +263,15 @@ def run_qc_controls(opts: QcControlsOptions, rep: Reporter) -> int:
|
||||
names = (COMPAT_RAW, COMPAT_QC) if opts.compat else (CONTROLS_RAW, CONTROLS_QC)
|
||||
targets = [opts.out / n for n in names]
|
||||
check_outputs_absent(targets, force=opts.force) # QC-014: other files in the folder are fine
|
||||
manifest = read_bpm(opts.bpm) if opts.bpm else None
|
||||
beads = load_control_beads(manifest, opts.csv, rep)
|
||||
if opts.array_type == "methylation":
|
||||
from simurg.formats.methyl_manifest import read_methyl_manifest
|
||||
|
||||
assert opts.csv is not None
|
||||
mman = read_methyl_manifest(opts.csv)
|
||||
beads = [ControlBead(c.type, c.color_channel, c.name, c.address, canonical_id(c.name)) for c in mman.controls]
|
||||
else:
|
||||
manifest = read_bpm(opts.bpm) if opts.bpm else None
|
||||
beads = load_control_beads(manifest, opts.csv, rep)
|
||||
if not beads:
|
||||
raise SimurgError("E_BPM_FORMAT", "the manifest has no control definitions")
|
||||
samples: list[ControlSample] = []
|
||||
@@ -273,7 +282,9 @@ def run_qc_controls(opts: QcControlsOptions, rep: Reporter) -> int:
|
||||
assert s.pair.grn is not None and s.pair.red is not None
|
||||
grn, red = read_idat(s.pair.grn), read_idat(s.pair.red)
|
||||
check_pair(grn, red)
|
||||
samples.append(control_sample(s.sample_id, grn, red, beads, opts.bg_offset, sheet_values(s, sset.sheet)))
|
||||
samples.append(
|
||||
control_sample(s.sample_id, grn, red, beads, opts.bg_offset, sheet_values(s, sset.sheet), opts.array_type)
|
||||
)
|
||||
# QC-014: sample IDs of an existing summary must match
|
||||
summary = opts.out / "sample_summary.tsv"
|
||||
if summary.exists():
|
||||
@@ -282,7 +293,9 @@ def run_qc_controls(opts: QcControlsOptions, rep: Reporter) -> int:
|
||||
if ids ^ mine:
|
||||
rep.warn("W_QC_UNMATCHED_SAMPLES", f"{len(ids ^ mine)} sample ID(s) differ from {summary}")
|
||||
opts.out.mkdir(parents=True, exist_ok=True)
|
||||
written = write_control_tables(opts.out, samples, beads, opts.bg_offset, rep, compat=opts.compat)
|
||||
written = write_control_tables(
|
||||
opts.out, samples, beads, opts.bg_offset, rep, array_type=opts.array_type, compat=opts.compat
|
||||
)
|
||||
for p in written:
|
||||
rm.add_output(p, "simurg.controls/1")
|
||||
rm.summary = {"samples": len(samples)}
|
||||
|
||||
@@ -593,8 +593,6 @@ def run_qc_report(opts: QcReportOptions, rep: Reporter) -> int:
|
||||
"Genotyping and methylation datasets cannot be mixed: "
|
||||
+ "; ".join(f"{t}: {', '.join(v)}" for t, v in by_type.items()),
|
||||
)
|
||||
if "methylation" in types:
|
||||
raise SimurgError("E_QC_METHYLATION", "methylation datasets in the report arrive with v2 (QC-062)")
|
||||
total = sum(len(ds.samples) for ds in datasets)
|
||||
if total > 20000:
|
||||
rep.warn("W_QC_TOO_FEW", f"{total} samples in total: the report may be slow to open (QC-057)")
|
||||
|
||||
@@ -51,6 +51,20 @@ canonical:
|
||||
bs conversion i u3: bs1_U3
|
||||
bs conversion i u4: bs1_U4
|
||||
bs conversion i u5: bs1_U5
|
||||
bs conversion i-c1: bs1_C1
|
||||
bs conversion i-c2: bs1_C2
|
||||
bs conversion i-c3: bs1_C3
|
||||
bs conversion i-c4: bs1_C4
|
||||
bs conversion i-c5: bs1_C5
|
||||
bs conversion i-c6: bs1_C6
|
||||
bs conversion i-u1: bs1_U1
|
||||
bs conversion i-u2: bs1_U2
|
||||
bs conversion i-u3: bs1_U3
|
||||
bs conversion i-u4: bs1_U4
|
||||
bs conversion i-u5: bs1_U5
|
||||
bs conversion i-u6: bs1_U6
|
||||
bs conversion i c6: bs1_C6
|
||||
bs conversion i u6: bs1_U6
|
||||
bs conversion ii-1: bs2_1
|
||||
bs conversion ii-2: bs2_2
|
||||
bs conversion ii-3: bs2_3
|
||||
|
||||
@@ -71,6 +71,10 @@ METRICS: tuple[Metric, ...] = (
|
||||
Metric("bs_conv_2", "bisulfiteConversionII", False, True),
|
||||
Metric("bs_conv_2_bg", "bisulfiteConversionIIBg", False, True),
|
||||
)
|
||||
GTMM_PM = [f"gtmm_{i}_pm" for i in range(1, 7)]
|
||||
GTMM_MM = [f"gtmm_{i}_mm" for i in range(1, 7)]
|
||||
SPEC2 = ["spec2_1", "spec2_2", "spec2_3"]
|
||||
BS2 = ["bs2_1", "bs2_2", "bs2_3", "bs2_4"]
|
||||
ASSAY_CONTROL_COLUMN = "assay_control_norm_r" # QC-005: copied from the sample summary
|
||||
COMPAT_METRIC_NAMES = {
|
||||
"restoration": "RestorationQC",
|
||||
@@ -223,14 +227,12 @@ def compute_control_metrics(
|
||||
red: npt.NDArray[np.float64],
|
||||
array_type: str = "genotyping",
|
||||
bg_offset: float = DEFAULT_BG_OFFSET,
|
||||
bs_conv_red_bg_set: str = "U3-U5",
|
||||
) -> ControlQC:
|
||||
"""QC-003/QC-004 metrics of one sample. Single-probe values are bead means; min/max run over all beads (QC-002)."""
|
||||
if array_type != "genotyping":
|
||||
from simurg.core.errors import SimurgError
|
||||
|
||||
raise SimurgError("E_QC_METHYLATION", "methylation control metrics are planned for v2 (QC-006)")
|
||||
p = _Probes(beads, green, red)
|
||||
x = bg_offset
|
||||
bs_red_u = ["bs1_U3", "bs1_U4", "bs1_U5"] if bs_conv_red_bg_set == "U3-U5" else ["bs1_U4", "bs1_U5", "bs1_U6"]
|
||||
formulas: dict[str, tuple[Any, Any]] = {
|
||||
# column: (numerator thunk, denominator thunk); each thunk records the probes it needs
|
||||
"restoration": (lambda: p.mean("restoration", "G"), lambda: p.max(["ext_A", "ext_T"], "G") + x),
|
||||
@@ -246,6 +248,18 @@ def compute_control_metrics(
|
||||
"nsb_red": (lambda: p.max(["ext_C", "ext_G"], "R") + x, lambda: p.max(["nsb_bkg"], "R")),
|
||||
"np_green": (lambda: p.min(["np_C", "np_G"], "G"), lambda: p.max(["np_A", "np_T"], "G")),
|
||||
"np_red": (lambda: p.min(["np_A", "np_T"], "R"), lambda: p.max(["np_C", "np_G"], "R")),
|
||||
# methylation (QC-004 rows 9, 15-24; MET-006)
|
||||
"target_removal_2": (lambda: p.max(["ext_A", "ext_T"], "G") + x, lambda: p.mean("target_removal_2", "G")),
|
||||
"specificity_1_green": (lambda: p.min(GTMM_PM[:3], "G"), lambda: p.max(GTMM_MM[:3], "G")),
|
||||
"specificity_1_red": (lambda: p.min(GTMM_PM[3:], "R"), lambda: p.max(GTMM_MM[3:], "R")),
|
||||
"specificity_2": (lambda: p.min(SPEC2, "R"), lambda: p.max(SPEC2, "G")),
|
||||
"specificity_2_bg": (lambda: p.max(["ext_A", "ext_T"], "G") + x, lambda: p.max(SPEC2, "G")),
|
||||
"bs_conv_1_green": (lambda: p.min(["bs1_C1", "bs1_C2"], "G"), lambda: p.max(["bs1_U1", "bs1_U2"], "G")),
|
||||
"bs_conv_1_green_bg": (lambda: p.max(["ext_A", "ext_T"], "G") + x, lambda: p.max(["bs1_U1", "bs1_U2"], "G")),
|
||||
"bs_conv_1_red": (lambda: p.min(["bs1_C3", "bs1_C4", "bs1_C5"], "R"), lambda: p.max(bs_red_u, "R")),
|
||||
"bs_conv_1_red_bg": (lambda: p.max(["ext_C", "ext_G"], "R") + x, lambda: p.max(bs_red_u, "R")),
|
||||
"bs_conv_2": (lambda: p.min(BS2, "R"), lambda: p.max(BS2, "G")),
|
||||
"bs_conv_2_bg": (lambda: p.max(["ext_A", "ext_T"], "G") + x, lambda: p.max(BS2, "G")),
|
||||
}
|
||||
res = ControlQC(values={})
|
||||
for m in metrics_for(array_type):
|
||||
|
||||
@@ -0,0 +1,288 @@
|
||||
"""Methylation signal processing — spec MET-005, MET-010 … MET-012 (plan §7, v2). Pure functions, no I/O.
|
||||
|
||||
* Raw signals (MET-005): Infinium I — M = channel[AddressB], U = channel[AddressA] in the probe's colour channel,
|
||||
out-of-band = the other channel; Infinium II — M = Grn[AddressA], U = Red[AddressA]. Addresses missing from the
|
||||
IDAT or read from < 3 beads are masked.
|
||||
* Detection p-value, pOOBAH [LIT: Zhou, Laird & Shen 2017 NAR 45:e22; OSS: sesame `pOOBAH`, MIT]: empirical CDF of
|
||||
the out-of-band Infinium I signal plus the negative controls per channel; p = 1 − max(F(M), F(U)) with each signal
|
||||
read in its own channel (Infinium II: M green, U red).
|
||||
* Background correction, noob [LIT: Triche et al. 2013 NAR 41:e90; OSS: sesame `noob`, MIT; limma `normexp.signal`,
|
||||
GPL — formula only]: normal-exponential convolution per channel, background = out-of-band + negative controls,
|
||||
µ and σ from Huber's M-estimator (MASS `huber`, k = 1.5, MAD scale), α = max(Huber(in-band) − µ, 10); offset 15.
|
||||
* Dye bias (MET-010 step 3) [OURS]: linear — scale red and green so the means of the normalization controls
|
||||
(NORM_A/T red, NORM_C/G green) meet at their average; nonlinear — quantile matching of the Infinium I red and green
|
||||
in-band signals to their average quantile function.
|
||||
* β = M / (M + U + α) with signals floored at 1; M-value = log2(β' / (1 − β')), β' clipped to [1e-6, 1 − 1e-6].
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import math
|
||||
from dataclasses import dataclass
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
from simurg.formats.idat import Idat
|
||||
from simurg.formats.methyl_manifest import MethylManifest
|
||||
from simurg.signal.raw_xy import AddressIndex
|
||||
|
||||
F64 = npt.NDArray[np.float64]
|
||||
Mask = npt.NDArray[np.bool_]
|
||||
MIN_BEADS = 3
|
||||
NOOB_OFFSET = 15.0
|
||||
|
||||
|
||||
@dataclass
|
||||
class MethylSignals:
|
||||
"""In-band M/U per probe (float64, NaN where masked), out-of-band signal of Infinium I probes, control
|
||||
intensities (manifest control order) and the probe layout needed to split them by channel."""
|
||||
|
||||
m: F64
|
||||
u: F64
|
||||
oob_m: F64 # Infinium I: the other channel at the M address; NaN for Infinium II
|
||||
oob_u: F64
|
||||
masked: Mask
|
||||
col: npt.NDArray[np.str_] # "R" / "G" / "" (Infinium II)
|
||||
ctrl_green: F64
|
||||
ctrl_red: F64
|
||||
ctrl_type: list[str]
|
||||
|
||||
@property
|
||||
def type_i(self) -> Mask:
|
||||
return np.asarray(self.col != "", dtype=np.bool_)
|
||||
|
||||
|
||||
def assemble(grn: Idat, red: Idat, man: MethylManifest, *, min_beads: int = MIN_BEADS) -> MethylSignals:
|
||||
"""MET-005 raw signal assembly."""
|
||||
gi, ri = AddressIndex(grn), AddressIndex(red)
|
||||
|
||||
def read(addrs: npt.NDArray[np.int64]) -> tuple[F64, F64, Mask]:
|
||||
rg, fg = gi.lookup(addrs)
|
||||
rr, fr = ri.lookup(addrs)
|
||||
g = grn.mean[rg].astype(np.float64)
|
||||
r = red.mean[rr].astype(np.float64)
|
||||
ok = fg & fr & (grn.n_beads[rg] >= min_beads) & (red.n_beads[rr] >= min_beads) & (addrs != 0)
|
||||
return g, r, ok
|
||||
|
||||
ga, ra, oka = read(man.u)
|
||||
gb, rb, okb = read(np.where(man.m != 0, man.m, man.u))
|
||||
t1 = man.m != 0
|
||||
is_r = man.col == "R"
|
||||
# Infinium II: M green, U red at address A; Infinium I: both addresses in the probe's channel
|
||||
m = np.where(t1, np.where(is_r, rb, gb), ga)
|
||||
u = np.where(t1, np.where(is_r, ra, ga), ra)
|
||||
oob_m = np.where(t1, np.where(is_r, gb, rb), np.nan)
|
||||
oob_u = np.where(t1, np.where(is_r, ga, ra), np.nan)
|
||||
masked = ~(oka & np.where(t1, okb, True)) | (t1 & ~np.isin(man.col, ("R", "G")))
|
||||
for a in (m, u, oob_m, oob_u):
|
||||
a[masked] = np.nan
|
||||
caddr = np.array([c.address for c in man.controls], dtype=np.int64)
|
||||
cg, cr, cok = read(caddr)
|
||||
cg[~cok] = np.nan
|
||||
cr[~cok] = np.nan
|
||||
return MethylSignals(m, u, oob_m, oob_u, masked, man.col, cg, cr, [c.type for c in man.controls])
|
||||
|
||||
|
||||
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_)
|
||||
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]])
|
||||
return _finite(oob_g), _finite(oob_r)
|
||||
|
||||
|
||||
def in_band(s: MethylSignals) -> tuple[F64, F64]:
|
||||
"""In-band foreground per channel (green, red): Infinium I M/U in their channel, Infinium II M green / U red."""
|
||||
red_i, grn_i, t2 = s.col == "R", s.col == "G", s.col == ""
|
||||
g = np.concatenate([s.m[grn_i], s.u[grn_i], s.m[t2]])
|
||||
r = np.concatenate([s.m[red_i], s.u[red_i], s.u[t2]])
|
||||
return _finite(g), _finite(r)
|
||||
|
||||
|
||||
# ------------------------------------------------------------------------------------------------ detection
|
||||
def _ecdf(ref: F64): # type: ignore[no-untyped-def]
|
||||
srt = np.sort(ref)
|
||||
n = len(srt)
|
||||
|
||||
def f(x: F64) -> F64:
|
||||
if n == 0:
|
||||
return np.zeros(len(x))
|
||||
return np.searchsorted(srt, x, side="right") / n # R ecdf: P(X <= x)
|
||||
|
||||
return f
|
||||
|
||||
|
||||
def detection_p(s: MethylSignals) -> F64:
|
||||
"""pOOBAH detection p-value per probe; 1 where masked."""
|
||||
bg_g, bg_r = background(s)
|
||||
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)
|
||||
fm = np.where(red_i, fr(m), fg(m)) # Infinium II M is green, Infinium I green M is green
|
||||
fu = np.where(red_i | t2, fr(u), fg(u)) # Infinium II U is red
|
||||
p = 1.0 - np.maximum(fm, fu)
|
||||
p[s.masked] = 1.0
|
||||
return np.asarray(p, dtype=np.float64)
|
||||
|
||||
|
||||
# ------------------------------------------------------------------------------------------------ noob
|
||||
def huber(y: F64, k: float = 1.5, tol: float = 1e-6) -> tuple[float, float]:
|
||||
"""MASS::huber: Huber M-estimator of location with MAD scale (returns the location before the last update)."""
|
||||
y = _finite(np.asarray(y, dtype=np.float64))
|
||||
n = len(y)
|
||||
if n == 0:
|
||||
raise ValueError("huber: no finite values")
|
||||
mu = float(np.median(y))
|
||||
s = float(1.4826 * np.median(np.abs(y - mu)))
|
||||
if s == 0:
|
||||
raise ValueError("huber: cannot estimate scale, MAD is zero")
|
||||
for _ in range(1000):
|
||||
yy = np.minimum(np.maximum(mu - k * s, y), mu + k * s)
|
||||
mu1 = float(yy.sum() / n)
|
||||
if abs(mu - mu1) < tol * s:
|
||||
break
|
||||
mu = mu1
|
||||
return mu, s
|
||||
|
||||
|
||||
_erfc = np.frompyfunc(math.erfc, 1, 1)
|
||||
|
||||
|
||||
def _mills(z: F64) -> F64:
|
||||
"""φ(z)/Φ(z) (inverse Mills ratio), stable far into the lower tail."""
|
||||
z = np.asarray(z, dtype=np.float64)
|
||||
out = np.empty_like(z)
|
||||
lo = z < -30.0
|
||||
zz = z[~lo]
|
||||
phi = np.exp(-0.5 * zz * zz) / math.sqrt(2.0 * math.pi)
|
||||
big_phi = 0.5 * np.asarray(_erfc(-zz / math.sqrt(2.0)), dtype=np.float64)
|
||||
out[~lo] = phi / big_phi
|
||||
zl = z[lo]
|
||||
r = 1.0 / (zl * zl)
|
||||
out[lo] = -zl / (1.0 - r + 3.0 * r * r - 15.0 * r**3 + 105.0 * r**4)
|
||||
return out
|
||||
|
||||
|
||||
def normexp_signal(mu: float, sigma: float, alpha: float, x: F64) -> F64:
|
||||
"""limma ``normexp.signal`` (formula): E[signal | observed x] under x = N(µ, σ²) + Exp(α)."""
|
||||
if alpha <= 0 or sigma <= 0:
|
||||
raise ValueError("normexp: alpha and sigma must be positive")
|
||||
mu_sf = x - mu - sigma * sigma / alpha
|
||||
sig = mu_sf + sigma * _mills(mu_sf / sigma)
|
||||
ok = np.isfinite(sig)
|
||||
if (sig[ok] < 0).any():
|
||||
sig[ok] = np.maximum(sig[ok], 1e-6)
|
||||
return sig
|
||||
|
||||
|
||||
@dataclass
|
||||
class NoobFit:
|
||||
mu: float
|
||||
sigma: float
|
||||
alpha: float
|
||||
|
||||
|
||||
def noob_fit(ib: F64, oob: F64) -> NoobFit:
|
||||
ib = np.where(ib < 1, 1.0, ib)
|
||||
mu, s = huber(oob)
|
||||
return NoobFit(mu, s, max(huber(ib)[0] - mu, 10.0))
|
||||
|
||||
|
||||
def noob(s: MethylSignals, offset: float = NOOB_OFFSET) -> 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)
|
||||
ib_g, ib_r = in_band(s)
|
||||
fit_g, fit_r = noob_fit(ib_g, bg_g), noob_fit(ib_r, bg_r)
|
||||
|
||||
def corr(x: F64, f: NoobFit) -> F64:
|
||||
y = np.where(np.isfinite(x), np.maximum(x, 1.0), np.nan)
|
||||
out = np.full_like(y, np.nan)
|
||||
ok = np.isfinite(y)
|
||||
out[ok] = normexp_signal(f.mu, f.sigma, f.alpha, y[ok]) + offset
|
||||
return out
|
||||
|
||||
red_i, t2 = s.col == "R", s.col == ""
|
||||
fm = np.where(red_i, 1, 0) # channel of M: red for Infinium I red, else green
|
||||
fu = np.where(red_i | t2, 1, 0)
|
||||
m = np.where(fm == 1, corr(s.m, fit_r), corr(s.m, fit_g))
|
||||
u = np.where(fu == 1, corr(s.u, fit_r), corr(s.u, fit_g))
|
||||
out = MethylSignals(
|
||||
m, u, s.oob_m, s.oob_u, s.masked, s.col, corr(s.ctrl_green, fit_g), corr(s.ctrl_red, fit_r), s.ctrl_type
|
||||
)
|
||||
return out, fit_g, fit_r
|
||||
|
||||
|
||||
# ------------------------------------------------------------------------------------------------ dye bias
|
||||
def _scale(s: MethylSignals, fg: float, fr: float) -> MethylSignals:
|
||||
red_i, t2 = s.col == "R", s.col == ""
|
||||
m = s.m * np.where(red_i, fr, fg)
|
||||
u = s.u * np.where(red_i | t2, fr, fg)
|
||||
return MethylSignals(m, u, s.oob_m, s.oob_u, s.masked, s.col, s.ctrl_green * fg, s.ctrl_red * fr, s.ctrl_type)
|
||||
|
||||
|
||||
def dye_bias_linear(s: MethylSignals) -> MethylSignals:
|
||||
"""Scale the channels so the normalization-control means (NORM_A/T red, NORM_C/G green) meet at their average."""
|
||||
t = np.array([x.upper() for x in s.ctrl_type])
|
||||
red = _finite(s.ctrl_red[np.isin(t, ("NORM_A", "NORM_T"))])
|
||||
grn = _finite(s.ctrl_green[np.isin(t, ("NORM_C", "NORM_G"))])
|
||||
if not len(red) or not len(grn):
|
||||
raise ValueError("dye bias: normalization controls missing")
|
||||
mr, mg = float(red.mean()), float(grn.mean())
|
||||
ref = (mr + mg) / 2.0
|
||||
return _scale(s, ref / mg, ref / mr)
|
||||
|
||||
|
||||
def dye_bias_nonlinear(s: MethylSignals) -> MethylSignals:
|
||||
"""Quantile-match the Infinium I red and green in-band signals to their average quantile function; every signal
|
||||
of a channel is mapped through that channel's curve (linear extrapolation beyond the Infinium I range)."""
|
||||
red_i, grn_i, t2 = s.col == "R", s.col == "G", s.col == ""
|
||||
r_ref = np.sort(_finite(np.concatenate([s.m[red_i], s.u[red_i]])))
|
||||
g_ref = np.sort(_finite(np.concatenate([s.m[grn_i], s.u[grn_i]])))
|
||||
if len(r_ref) < 2 or len(g_ref) < 2:
|
||||
raise ValueError("dye bias: too few Infinium I probes")
|
||||
q = np.linspace(0, 1, 1001)
|
||||
rq, gq = np.quantile(r_ref, q), np.quantile(g_ref, q)
|
||||
target = (rq + gq) / 2.0
|
||||
|
||||
def mapper(src_q: F64): # type: ignore[no-untyped-def]
|
||||
lo_slope = target[0] / src_q[0] if src_q[0] > 0 else 1.0
|
||||
hi_slope = target[-1] / src_q[-1] if src_q[-1] > 0 else 1.0
|
||||
|
||||
def f(x: F64) -> F64:
|
||||
y = np.interp(x, src_q, target)
|
||||
y = np.where(x < src_q[0], x * lo_slope, y)
|
||||
return np.where(x > src_q[-1], x * hi_slope, y)
|
||||
|
||||
return f
|
||||
|
||||
fr, fg = mapper(rq), mapper(gq)
|
||||
m = np.where(red_i, fr(s.m), fg(s.m))
|
||||
u = np.where(red_i | t2, fr(s.u), fg(s.u))
|
||||
return MethylSignals(m, u, s.oob_m, s.oob_u, s.masked, s.col, fg(s.ctrl_green), fr(s.ctrl_red), s.ctrl_type)
|
||||
|
||||
|
||||
# ------------------------------------------------------------------------------------------------ beta
|
||||
def beta_values(s: MethylSignals, offset: float = 0.0) -> F64:
|
||||
m = np.maximum(s.m, 1.0)
|
||||
u = np.maximum(s.u, 1.0)
|
||||
return m / (m + u + offset)
|
||||
|
||||
|
||||
def m_values(beta: F64) -> F64:
|
||||
b = np.clip(beta, 1e-6, 1 - 1e-6)
|
||||
return np.log2(b / (1 - b))
|
||||
|
||||
|
||||
def average_signals(s: MethylSignals, analytical: Mask) -> tuple[float, float]:
|
||||
"""MET-012: mean green and red in-band signal over analytical probes."""
|
||||
red_i, grn_i, t2 = s.col == "R", s.col == "G", s.col == ""
|
||||
a = analytical
|
||||
g = _finite(np.concatenate([s.m[grn_i & a], s.u[grn_i & a], s.m[t2 & a]]))
|
||||
r = _finite(np.concatenate([s.m[red_i & a], s.u[red_i & a], s.u[t2 & a]]))
|
||||
return (float(g.mean()) if len(g) else math.nan, float(r.mean()) if len(r) else math.nan)
|
||||
Reference in new issue
Block a user