From 3eb51639cc38dde36184d9005da32b837f52d0ba Mon Sep 17 00:00:00 2001 From: Kemal Yaylali Date: Tue, 29 Sep 2026 00:13:22 +0100 Subject: [PATCH] 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). --- docs/messages.md | 11 +- src/simurg/cli/main.py | 66 ++++ src/simurg/core/errors.py | 32 +- src/simurg/formats/methyl_manifest.py | 157 +++++++++ src/simurg/formats/samplesheet.py | 10 +- src/simurg/pipeline/methyl.py | 443 ++++++++++++++++++++++++++ src/simurg/pipeline/qc.py | 27 +- src/simurg/pipeline/qc_report.py | 2 - src/simurg/qc/control_names.yaml | 14 + src/simurg/qc/controls.py | 22 +- src/simurg/signal/methyl.py | 288 +++++++++++++++++ tests/unit/test_methyl.py | 234 ++++++++++++++ tests/unit/test_qc.py | 7 +- 13 files changed, 1292 insertions(+), 21 deletions(-) create mode 100644 src/simurg/formats/methyl_manifest.py create mode 100644 src/simurg/pipeline/methyl.py create mode 100644 src/simurg/signal/methyl.py create mode 100644 tests/unit/test_methyl.py diff --git a/docs/messages.md b/docs/messages.md index 23d6255..b1d65c4 100644 --- a/docs/messages.md +++ b/docs/messages.md @@ -26,6 +26,12 @@ Generated from `src/simurg/core/errors.py` — do not edit by hand. | `E_IDAT_MANIFEST_MISMATCH` | error | IDATs and manifest belong to different products | IN-027 | | `E_IDAT_PAIR_MISMATCH` | error | Red and Grn IDATs of a sample do not match | IN-023 | | `E_INTERNAL` | error | Internal error (bug); rerun with --debug and report | CLI-014 | +| `E_METH_CONTROLS_FORMAT` | error | Controls not formatted correctly. Must have 4 columns (Address,Type,Color_Channel,Name) | MET-002 | +| `E_METH_MISSING_COLUMNS` | error | Missing required columns in the methylation manifest | MET-002 | +| `E_METH_NO_ASSAY` | error | Probe section not found | MET-002 | +| `E_METH_NO_CONTROLS` | error | Missing control probes in manifest | MET-002 | +| `E_METH_NO_IDATS` | error | No IDATs found | MET-004 | +| `E_METH_PREP_FAILED` | error | Prep failed (sample) | MET-010 | | `E_NO_INPUT` | error | Neither an input folder nor a sample sheet was given | CLI-006 | | `E_NO_INPUT_FOLDER` | error | Sample sheet has no RootFolder and no folder was given | CLI-006 | | `E_OUTPUT_EXISTS` | error | Outputs already exist; use --force or --resume | CLI-010 | @@ -37,7 +43,6 @@ Generated from `src/simurg/core/errors.py` — do not edit by hand. | `E_QC_FORMAT` | error | Invalid --output-format value | QC-063 | | `E_QC_LABEL_COUNT` | error | Number of labels does not match number of datasets | QC-042 | | `E_QC_LABEL_DUPLICATE` | error | Dataset labels are not unique | QC-042 | -| `E_QC_METHYLATION` | error | Methylation control metrics are not available yet (v2) | QC-006 | | `E_QC_MIXED_TYPES` | error | Genotyping and methylation datasets cannot be mixed | QC-043 | | `E_QC_NO_DATASETS` | error | No QC dataset folder found | QC-041 | | `E_QC_PROFILE` | error | Unknown QC profile | QC-024 | @@ -73,6 +78,10 @@ Generated from `src/simurg/core/errors.py` — do not edit by hand. | `W_MERGE_ALL_NC` | warning | Every probe at a merged site is a no-call | OUT-031 | | `W_MERGE_AMBIGUOUS` | warning | Probes at a site allow several genotypes; set missing | OUT-031 | | `W_MERGE_INCONSISTENT` | warning | Probes at a site disagree; genotype set missing | OUT-031 | +| `W_METH_MISSING_CONTROL_NAMES` | warning | Control names missing from the manifest; metric set to NA | MET-009 | +| `W_METH_MISSING_CONTROL_TYPES` | warning | Control types missing from the manifest; metric set to NA | MET-009 | +| `W_METH_MISSING_SAMPLE` | warning | Missing sample: IDATs absent; the sample fails, the run continues | MET-015 | +| `W_METH_NO_FEATURES` | warning | No features, skipping PCA plot | MET-013 | | `W_NO_BGZIP` | warning | Uncompressed VCF written; downstream commands need bgzip + index | OUT-010 | | `W_QC_LABEL_SANITIZED` | warning | Dataset label was sanitized | QC-042 | | `W_QC_NO_BAF` | warning | GTC has no B-allele frequencies; contamination not estimated | QC-064 | diff --git a/src/simurg/cli/main.py b/src/simurg/cli/main.py index e993de7..097a840 100644 --- a/src/simurg/cli/main.py +++ b/src/simurg/cli/main.py @@ -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, diff --git a/src/simurg/core/errors.py b/src/simurg/core/errors.py index 6606351..43b204a 100644 --- a/src/simurg/core/errors.py +++ b/src/simurg/core/errors.py @@ -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"), ] ) diff --git a/src/simurg/formats/methyl_manifest.py b/src/simurg/formats/methyl_manifest.py new file mode 100644 index 0000000..badf21f --- /dev/null +++ b/src/simurg/formats/methyl_manifest.py @@ -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=" 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 diff --git a/src/simurg/formats/samplesheet.py b/src/simurg/formats/samplesheet.py index 15cbe91..8b41d30 100644 --- a/src/simurg/formats/samplesheet.py +++ b/src/simurg/formats/samplesheet.py @@ -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])} diff --git a/src/simurg/pipeline/methyl.py b/src/simurg/pipeline/methyl.py new file mode 100644 index 0000000..a072e4e --- /dev/null +++ b/src/simurg/pipeline/methyl.py @@ -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("" + "".join(f"{html.escape(_cell(r.get(c)))}" for c in head) + "" for r in rows) + th_rows = "".join( + f"{m.column}{m.key}{_cell(th.get(m.key).value)}" for m in metrics + ) + generated = _dt.datetime.now(_dt.UTC).strftime("%Y-%m-%d %H:%M UTC") + page = f""" +Methylation QC + + +

Methylation control QC

+

{len(rows)} samples · manifest {html.escape(Path(man.path).name)} + · profile {html.escape(th.profile or "none")} · generated {generated} · SimurgArray

+
+

Samples

{"".join(f"" for c in head)}{trs}
{c}
+

Thresholds

{th_rows}
metrickeythreshold
+ +""" + 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 diff --git a/src/simurg/pipeline/qc.py b/src/simurg/pipeline/qc.py index 44ee319..979e41a 100644 --- a/src/simurg/pipeline/qc.py +++ b/src/simurg/pipeline/qc.py @@ -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)} diff --git a/src/simurg/pipeline/qc_report.py b/src/simurg/pipeline/qc_report.py index f1b6af9..2e4fc2e 100644 --- a/src/simurg/pipeline/qc_report.py +++ b/src/simurg/pipeline/qc_report.py @@ -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)") diff --git a/src/simurg/qc/control_names.yaml b/src/simurg/qc/control_names.yaml index 9967453..3d5e94b 100644 --- a/src/simurg/qc/control_names.yaml +++ b/src/simurg/qc/control_names.yaml @@ -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 diff --git a/src/simurg/qc/controls.py b/src/simurg/qc/controls.py index e879c15..fbd162c 100644 --- a/src/simurg/qc/controls.py +++ b/src/simurg/qc/controls.py @@ -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): diff --git a/src/simurg/signal/methyl.py b/src/simurg/signal/methyl.py new file mode 100644 index 0000000..0dae03e --- /dev/null +++ b/src/simurg/signal/methyl.py @@ -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) diff --git a/tests/unit/test_methyl.py b/tests/unit/test_methyl.py new file mode 100644 index 0000000..c78e65a --- /dev/null +++ b/tests/unit/test_methyl.py @@ -0,0 +1,234 @@ +"""Methylation control QC (spec MET-001 … MET-018).""" + +from __future__ import annotations + +import gzip +import math +from pathlib import Path + +import numpy as np +import pytest +from typer.testing import CliRunner + +from simurg.cli.main import app +from simurg.core.errors import FormatError +from simurg.formats.methyl_manifest import read_methyl_manifest, write_probe_tables +from simurg.qc.controls import ControlBead, canonical_id, compute_control_metrics +from simurg.signal import methyl as sig +from simurg.synth.files import SynthIdat, encode_idat + +HEAD = "Illumina, Inc.,,,\n[Heading],,,\nDescriptor File Name,SYNTH-METH_A1.csv,,\n[Assay],,,\n" +COLS = "IlmnID,Name,AddressA_ID,AlleleA_ProbeSeq,AddressB_ID,AlleleB_ProbeSeq,Infinium_Design_Type,Color_Channel\n" +CONTROLS = [ # (type, colour, name); one bead each + ("STAINING", "Red", "DNP (High)"), ("STAINING", "Purple", "DNP (Bkg)"), ("STAINING", "Green", "Biotin (High)"), + ("STAINING", "Blue", "Biotin (Bkg)"), ("EXTENSION", "Red", "Extension (A)"), ("EXTENSION", "Purple", + "Extension (T)"), ("EXTENSION", "Green", "Extension (C)"), ("EXTENSION", "Blue", "Extension (G)"), + ("HYBRIDIZATION", "Green", "Hyb (High)"), ("HYBRIDIZATION", "Blue", "Hyb (Medium)"), ("HYBRIDIZATION", "Black", + "Hyb (Low)"), ("TARGET REMOVAL", "Green", "Target Removal 1"), ("TARGET REMOVAL", "Blue", "Target Removal 2"), + ("RESTORATION", "Green", "Restore"), ("NON-POLYMORPHIC", "Red", "NP (A)"), ("NON-POLYMORPHIC", "Purple", "NP (T)"), + ("NON-POLYMORPHIC", "Green", "NP (C)"), ("NON-POLYMORPHIC", "Blue", "NP (G)"), + *[("SPECIFICITY I", "Green", f"GT Mismatch {i} (PM)") for i in range(1, 7)], + *[("SPECIFICITY I", "Blue", f"GT Mismatch {i} (MM)") for i in range(1, 7)], + *[("SPECIFICITY II", "Red", f"Specificity {i}") for i in range(1, 4)], + *[("BISULFITE CONVERSION I", "Green", f"BS Conversion I-C{i}") for i in range(1, 6)], + *[("BISULFITE CONVERSION I", "Blue", f"BS Conversion I-U{i}") for i in range(1, 6)], + *[("BISULFITE CONVERSION II", "Red", f"BS Conversion II-{i}") for i in range(1, 5)], + *[("NEGATIVE", "Red", f"Negative {i}") for i in range(1, 41)], + *[(f"NORM_{b}", "Red", f"Norm_{b}{i}") for b in "ACGT" for i in range(1, 6)], +] # fmt: skip +N_T1R, N_T1G, N_T2 = 60, 60, 180 + + +def _manifest(tmp: Path, *, controls: list[tuple[str, str, str]] | None = None, extra_field: bool = False) -> Path: + lines = [HEAD, COLS] + addr = 1000 + for i in range(N_T1R + N_T1G + N_T2): + if i < N_T1R + N_T1G: + ch = "Red" if i < N_T1R else "Grn" + lines.append(f"cg{i:08d},cg{i:08d},{addr},ACGT,{addr + 1},ACGT,I,{ch}\n") + addr += 2 + else: + lines.append(f"cg{i:08d},cg{i:08d},{addr},ACGT,,,II,\n") + addr += 1 + lines.append("rs0000001,rs0000001,9001,ACGT,,,II,\n") + lines.append("ch.1.1,ch.1.1,9002,ACGT,,,II,\n") + lines.append("[Controls],,,\n") + for k, (t, c, n) in enumerate(controls if controls is not None else CONTROLS): + tail = ",AVG" if extra_field and t == "NEGATIVE" else "" + lines.append(f"{20000 + k},{t},{c},{n}{tail},,,,\n") + p = tmp / "meth.csv" + p.write_text("".join(lines)) + return p + + +def _idats(folder: Path, man_path: Path, seed: int, barcode: str) -> None: + man = read_methyl_manifest(man_path) + rng = np.random.default_rng(seed) + addrs = np.unique(np.concatenate([man.u, man.m[man.m > 0], [c.address for c in man.controls]])).astype(np.int32) + n = len(addrs) + for ch in ("Grn", "Red"): + mean = rng.gamma(2.0, 2500.0, n).clip(50, 60000).astype(np.uint16) + folder.mkdir(parents=True, exist_ok=True) + (folder / f"{barcode}_R01C01_{ch}.idat").write_bytes( + encode_idat(SynthIdat(addrs, mean, n_beads=np.full(n, 8, np.uint8), barcode=barcode)) + ) + + +# ------------------------------------------------------------------------------------------------ manifest +def test_manifest_parse_and_tables(tmp_path: Path) -> None: # MET-002, MET-003 + man = read_methyl_manifest(_manifest(tmp_path, extra_field=True)) + assert man.num_probes == N_T1R + N_T1G + N_T2 + 2 + assert int(man.type_i.sum()) == N_T1R + N_T1G + assert (man.col[:N_T1R] == "R").all() and (man.col[N_T1R : N_T1R + N_T1G] == "G").all() + assert (man.col[N_T1R + N_T1G :] == "").all() + assert int(man.analytical().sum()) == N_T1R + N_T1G + N_T2 + 1 # cg + rs, not ch + assert len(man.controls) == len(CONTROLS) and man.controls[0].probe_id == "ctl_20000" + probes, ctl = write_probe_tables(man, tmp_path / "t") + lines = probes.read_text().splitlines() + assert lines[0] == "probe_id\tu\tm\tcol" + assert lines[1] == "cg00000000\t1000\t1001\tR" + assert lines[N_T1R + N_T1G + 1].endswith("\tNA\tNA") + assert ctl.read_text().splitlines()[1] == "20000\tSTAINING\tRed\tDNP (High)\tctl_20000" + + +@pytest.mark.parametrize( + ("text", "code"), + [ + ("[Heading],\nx,y\n[Controls]\n1,NEGATIVE,Red,N1\n", "E_METH_NO_ASSAY"), + ("[Assay]\n" + COLS + "cg1,cg1,1,A,,,II,\n", "E_METH_NO_CONTROLS"), + ("[Assay]\nIlmnID,Name,AddressA_ID\ncg1,cg1,1\n[Controls]\n1,NEGATIVE,Red,N1\n", "E_METH_MISSING_COLUMNS"), + ("[Assay]\n" + COLS + "cg1,cg1,1,A,,,II,\n[Controls]\n1,NEGATIVE,Red\n", "E_METH_CONTROLS_FORMAT"), + ], +) +def test_manifest_errors(tmp_path: Path, text: str, code: str) -> None: + p = tmp_path / "bad.csv" + p.write_text(text) + with pytest.raises(FormatError) as e: + read_methyl_manifest(p) + assert e.value.code == code + + +# ------------------------------------------------------------------------------------------------ numerics +def test_huber_matches_mass() -> None: + # MASS::huber(c(1, 2, 3, 4, 100)): mu is the fixed point of mu = (1 + 2 + 3 + 4 + (mu + 1.5 * s)) / 5 + mu, s = sig.huber(np.array([1.0, 2.0, 3.0, 4.0, 100.0])) + assert s == pytest.approx(1.4826) + assert mu == pytest.approx((10 + 1.5 * 1.4826) / 4, abs=1e-5) + with pytest.raises(ValueError): + sig.huber(np.ones(5)) + + +def test_mills_ratio() -> None: + assert sig._mills(np.array([0.0]))[0] == pytest.approx(math.sqrt(2 / math.pi)) + lo, hi = sig._mills(np.array([-30.0001, -29.9999])) + assert lo == pytest.approx(hi, rel=1e-5) + z = np.array([-5.0, -1.0, 2.0]) + phi = np.exp(-z * z / 2) / math.sqrt(2 * math.pi) + big = np.array([0.5 * math.erfc(-v / math.sqrt(2)) for v in z]) + assert np.allclose(sig._mills(z), phi / big, rtol=1e-12) + + +def test_normexp_signal() -> None: + x = np.array([1.0, 500.0, 5000.0, 60000.0]) + out = sig.normexp_signal(400.0, 150.0, 3000.0, x) + assert (out > 0).all() and (np.diff(out) > 0).all() + assert out[-1] == pytest.approx(60000.0 - 400.0 - 150.0**2 / 3000.0, rel=1e-9) + + +def test_detection_and_beta() -> None: # MET-010, MET-011 + s = sig.MethylSignals( + m=np.array([5000.0, 100.0, 3000.0, np.nan]), + u=np.array([100.0, 90.0, 4000.0, np.nan]), + oob_m=np.array([200.0, 150.0, np.nan, np.nan]), + oob_u=np.array([250.0, 300.0, np.nan, np.nan]), + masked=np.array([False, False, False, True]), + col=np.array(["R", "G", "", ""]), + ctrl_green=np.array([120.0, 130.0]), + ctrl_red=np.array([110.0, 140.0]), + ctrl_type=["NEGATIVE", "NEGATIVE"], + ) + p = sig.detection_p(s) + assert p[0] == 0.0 # far above the red background + assert p[1] > 0.5 # at the green background + assert p[3] == 1.0 # masked + b = sig.beta_values(s) + assert b[0] == pytest.approx(5000 / 5100) + assert sig.beta_values(s, 100.0)[0] == pytest.approx(5000 / 5200) + assert np.isnan(b[3]) + mv = sig.m_values(np.array([0.0, 0.5, 1.0])) + assert mv[1] == 0.0 and mv[0] == pytest.approx(math.log2(1e-6 / (1 - 1e-6))) + + +def test_methylation_metrics(tmp_path: Path) -> None: # MET-006 + man = read_methyl_manifest(_manifest(tmp_path)) + beads = [ControlBead(c.type, c.color_channel, c.name, c.address, canonical_id(c.name)) for c in man.controls] + assert all(b.canonical for b in beads if not b.category.startswith(("NEGATIVE", "NORM"))) + g = np.full(len(beads), 1000.0) + r = np.full(len(beads), 1000.0) + idx = {b.name: i for i, b in enumerate(beads)} + for nm, v in (("BS Conversion I-C3", 8000.0), ("BS Conversion I-U3", 400.0), ("BS Conversion I-U5", 500.0)): + r[idx[nm]] = v + r[idx["BS Conversion I-C4"]] = r[idx["BS Conversion I-C5"]] = 9000.0 + q = compute_control_metrics(beads, g, r, "methylation", 3000.0) + assert set(q.values) == {m for m in q.values} and len(q.values) == 21 + assert q.values["bs_conv_1_red"] == pytest.approx(8000 / 1000) # max(U3..U5) = 1000 (U4) + assert q.values["target_removal_2"] == pytest.approx((1000 + 3000) / 1000) + q2 = compute_control_metrics(beads, g, r, "methylation", 3000.0, "U4-U6") + assert q2.missing.get("bs_conv_1_red") == ["bs1_U6"] + + +# ------------------------------------------------------------------------------------------------ CLI +def test_methyl_qc_end_to_end(tmp_path: Path) -> None: # MET-013 + man = _manifest(tmp_path) + for k in range(3): + _idats(tmp_path / "idats" / f"20400000000{k}", man, k, f"20400000000{k}") + out = tmp_path / "out" + r = CliRunner().invoke(app, ["methyl", "qc", "--manifest", str(man), "--idats", str(tmp_path / "idats"), + "--out", str(out), "--threads", "1", "--quiet"]) # fmt: skip + assert r.exit_code == 0, r.output + for f in ("methyl_qc_summary.tsv", "methyl_qc_summary.xlsx", "methyl_qc_summary.html", "pcs.tsv.gz", + "probes.tsv", "controls.tsv", "controls_raw.tsv", "controls_qc.tsv", "run_manifest.json"): # fmt: skip + assert (out / f).exists(), f + rows = [ln.split("\t") for ln in (out / "methyl_qc_summary.tsv").read_text().splitlines()] + head = rows[0] + assert head[:3] == ["sentrix_id", "sentrix_position", "sample_id"] and len(rows) == 4 + for c in ("bs_conv_2_bg", "avg_green_norm", "n_probes", "prop_probes_passing", "pass_qc", "fail_codes"): + assert c in head + assert {r_[head.index("n_probes")] for r_ in rows[1:]} == {str(N_T1R + N_T1G + N_T2 + 1)} + sid = rows[1][2] + with gzip.open(out / "samples" / sid / f"{sid}_cgs.tsv.gz", "rt") as fh: + cgs = fh.read().splitlines() + assert cgs[0] == "probe_id\tbeta\tm_value\tdetection_p" and len(cgs) == N_T1R + N_T1G + N_T2 + 3 + with gzip.open(out / "pcs.tsv.gz", "rt") as fh: + assert len(fh.read().splitlines()) == 4 + # the folder is a valid qc report dataset (MET-014) + r = CliRunner().invoke(app, ["qc", "report", "--data", str(out), "--out", str(tmp_path / "rep"), "--quiet"]) + assert r.exit_code == 0, r.output + + +def test_methyl_qc_missing_types_and_no_idats(tmp_path: Path) -> None: # MET-009, MET-004 + man = _manifest(tmp_path, controls=[c for c in CONTROLS if c[0] != "RESTORATION"]) + _idats(tmp_path / "idats" / "204000000009", man, 9, "204000000009") + out = tmp_path / "out" + r = CliRunner().invoke(app, ["methyl", "qc", "--manifest", str(man), "--idats", str(tmp_path / "idats"), + "--out", str(out), "--threads", "1"]) # fmt: skip + assert r.exit_code == 0, r.output + assert "W_METH_MISSING_CONTROL_NAMES" in r.output + assert "restoration" in (out / "methyl_qc_summary.tsv").read_text().splitlines()[1] + (tmp_path / "empty").mkdir() + r = CliRunner().invoke(app, ["methyl", "qc", "--manifest", str(man), "--idats", str(tmp_path / "empty"), + "--out", str(tmp_path / "o2")]) # fmt: skip + assert r.exit_code == 2 and "E_METH_NO_IDATS" in r.output + + +def test_convert_manifest_and_qc_controls(tmp_path: Path) -> None: # MET-003, QC-012 + man = _manifest(tmp_path) + r = CliRunner().invoke(app, ["methyl", "convert-manifest", "--manifest", str(man), "--out", str(tmp_path / "c")]) + assert r.exit_code == 0, r.output + assert (tmp_path / "c" / "probes.tsv").exists() + _idats(tmp_path / "idats" / "204000000001", man, 1, "204000000001") + r = CliRunner().invoke(app, ["qc", "controls", "--array-type", "methylation", "--csv", str(man), "--idats", + str(tmp_path / "idats"), "--out", str(tmp_path / "q"), "--quiet"]) # fmt: skip + assert r.exit_code == 0, r.output + qc = (tmp_path / "q" / "controls_qc.tsv").read_text().splitlines() + assert "bs_conv_2_bg" in qc[0] and "stringency" not in qc[0] and qc[1].split("\t")[1] == "methylation" diff --git a/tests/unit/test_qc.py b/tests/unit/test_qc.py index aa45e00..4327395 100644 --- a/tests/unit/test_qc.py +++ b/tests/unit/test_qc.py @@ -167,10 +167,9 @@ def test_absent_addresses_are_skipped() -> None: # QC-008, K-021: two bead-pool assert q.missing["staining_green"] == ["biotin_bkg"] -def test_methylation_is_v2() -> None: - with pytest.raises(SimurgError) as e: - compute_control_metrics(_beads(), np.ones(21), np.ones(21), "methylation") - assert e.value.code == "E_QC_METHYLATION" +def test_methylation_metric_set() -> None: # QC-006: methylation metrics are computed (v2); missing ones are NaN + q = compute_control_metrics(_beads(), np.ones(21), np.ones(21), "methylation") + assert len(q.values) == 21 and "stringency" not in q.values and "bs_conv_2" in q.missing # ------------------------------------------------------------------------------------------------ thresholds