kemal
1fc27a866a
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).
235 lines
12 KiB
Python
235 lines
12 KiB
Python
"""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"
|