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).
This commit is contained in:
kemal committed
1 parent
7f3cd96fc4
commit
1fc27a866a
13 files changed
+1292
-21
No files matched your search
@@ -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"
|
||||
Reference in new issue
Block a user