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).
396 lines
17 KiB
Python
396 lines
17 KiB
Python
"""Phase P5 QC unit tests: control metrics (QC-001…QC-008), thresholds (QC-020…QC-025), status (QC-030…QC-033),
|
||
contamination (plan §9.3), relatedness (plan §9.4), locus summary (OUT-072), suggestions (QC-026)."""
|
||
|
||
from __future__ import annotations
|
||
|
||
import math
|
||
from pathlib import Path
|
||
|
||
import numpy as np
|
||
import pytest
|
||
|
||
from simurg.core.errors import SimurgError
|
||
from simurg.qc.contamination import _median_exact, _median_reference, b_allele_frequency, baf_regress, ks_ksmall
|
||
from simurg.qc.controls import (
|
||
beads_from_definitions,
|
||
canonical_id,
|
||
compute_control_metrics,
|
||
excluded_from_report,
|
||
normalise_name,
|
||
)
|
||
from simurg.qc.locus_summary import gs_percentile_columns, summarize_block
|
||
from simurg.qc.relatedness import degree, mendel_errors, pair_stats, select_sites
|
||
from simurg.qc.status import control_status, functional_status, het_zscores, sample_flags
|
||
from simurg.qc.suggest import suggest
|
||
from simurg.qc.thresholds import auto_profile, parse_cli_pairs, profile_values, resolve_thresholds
|
||
|
||
GSA_CONTROLS = [
|
||
("1", "Staining", "Red", "DNP (High)"),
|
||
("2", "Staining", "Purple", "DNP (Bgnd)"),
|
||
("3", "Staining", "Green", "Biotin (High)"),
|
||
("4", "Staining", "Blue", "Biotin (Bgnd)"),
|
||
("5", "Extension", "Red", "Extension (A)"),
|
||
("6", "Extension", "Purple", "Extension (T)"),
|
||
("7", "Extension", "Green", "Extension (C)"),
|
||
("8", "Extension", "Blue", "Extension (G)"),
|
||
("9", "Target Removal", "Green", "Target Removal"),
|
||
("10", "Hybridization", "Green", "Hyb (High)"),
|
||
("11", "Hybridization", "Blue", "Hyb (Medium)"),
|
||
("12", "Hybridization", "Black", "Hyb (Low)"),
|
||
("13", "Stringency", "Red", "String (PM)"),
|
||
("14", "Stringency", "Purple", "String (MM)"),
|
||
("15:15:15", "Non-Specific Binding", "Red", "NSB (Bgnd)"),
|
||
("16", "Non-Specific Binding", "Purple", "NSB (Bgnd)"),
|
||
("17", "Non-Polymorphic", "Red", "NP(A)"),
|
||
("18", "Non-Polymorphic", "Purple", "NP(T)"),
|
||
("19", "Non-Polymorphic", "Green", "NP(C)"),
|
||
("20", "Non-Polymorphic", "Blue", "NP(G)"),
|
||
("21", "Restoration", "Green", "Restore"),
|
||
]
|
||
|
||
|
||
def _beads(): # type: ignore[no-untyped-def]
|
||
return beads_from_definitions(
|
||
[{"address": a, "category": c, "color": col, "name": n} for a, c, col, n in GSA_CONTROLS]
|
||
)
|
||
|
||
|
||
# ------------------------------------------------------------------------------------------------ controls
|
||
@pytest.mark.parametrize(
|
||
("name", "cid"),
|
||
[
|
||
("Extension (A)", "ext_A"),
|
||
("Biotin (High)", "biotin_high"),
|
||
("DNP (Bgnd)", "dnp_bkg"),
|
||
("NP (A)", "np_A"),
|
||
("NP(A)", "np_A"),
|
||
("String (PM)", "string_pm"),
|
||
("NSB (Bgnd) Red", "nsb_bkg"),
|
||
("Hyb ( High )", "hyb_high"),
|
||
("Target Removal", "target_removal_1"),
|
||
("Restore", "restoration"),
|
||
("GT Mismatch 1 (PM)", "gtmm_1_pm"),
|
||
("Something else", ""),
|
||
],
|
||
)
|
||
def test_canonical_names(name: str, cid: str) -> None: # QC-002
|
||
assert canonical_id(name) == cid
|
||
|
||
|
||
def test_normalise_and_exclusion() -> None: # QC-002, QC-010
|
||
assert normalise_name(" NP (G) 1 ") == "np(g) 1"
|
||
b = beads_from_definitions([{"address": "1", "category": "NEGATIVE", "color": "Red", "name": "x"}])[0]
|
||
assert excluded_from_report(b)
|
||
|
||
|
||
def test_bead_expansion_dedups_addresses() -> None: # QC-001
|
||
beads = _beads()
|
||
assert len(beads) == len(GSA_CONTROLS) # "15:15:15" is one bead type
|
||
assert beads[14].address == 15
|
||
|
||
|
||
def test_control_metrics_formulas() -> None: # QC-003, QC-004
|
||
beads = _beads()
|
||
g = np.zeros(len(beads))
|
||
r = np.zeros(len(beads))
|
||
val = {b.name + b.color: i for i, b in enumerate(beads)}
|
||
|
||
def put(name: str, color: str, green: float, red: float) -> None:
|
||
g[val[name + color]], r[val[name + color]] = green, red
|
||
|
||
put("DNP (High)", "Red", 100, 8000)
|
||
put("DNP (Bgnd)", "Purple", 90, 800)
|
||
put("Biotin (High)", "Green", 9000, 50)
|
||
put("Biotin (Bgnd)", "Blue", 900, 50)
|
||
put("Extension (A)", "Red", 300, 9000)
|
||
put("Extension (T)", "Purple", 400, 8000)
|
||
put("Extension (C)", "Green", 7000, 200)
|
||
put("Extension (G)", "Blue", 6000, 250)
|
||
put("Target Removal", "Green", 1700, 10)
|
||
put("Hyb (High)", "Green", 6000, 10)
|
||
put("Hyb (Medium)", "Blue", 3000, 10)
|
||
put("Hyb (Low)", "Black", 1000, 10)
|
||
put("String (PM)", "Red", 10, 3000)
|
||
put("String (MM)", "Purple", 10, 1000)
|
||
put("NSB (Bgnd)", "Red", 500, 400)
|
||
put("NSB (Bgnd)", "Purple", 700, 300)
|
||
put("NP(A)", "Red", 100, 5000)
|
||
put("NP(T)", "Purple", 120, 6000)
|
||
put("NP(C)", "Green", 4000, 300)
|
||
put("NP(G)", "Blue", 5000, 200)
|
||
put("Restore", "Green", 1360, 0)
|
||
q = compute_control_metrics(beads, g, r, bg_offset=3000)
|
||
v = q.values
|
||
bkg_grn, bkg_red = 400.0, 250.0
|
||
assert v["restoration"] == pytest.approx(1360 / (bkg_grn + 3000))
|
||
assert v["staining_green"] == pytest.approx(10.0)
|
||
assert v["staining_red"] == pytest.approx(10.0)
|
||
assert v["extension_green"] == pytest.approx(6000 / 400)
|
||
assert v["extension_red"] == pytest.approx(8000 / 250)
|
||
assert v["hyb_high_med"] == pytest.approx(2.0)
|
||
assert v["hyb_med_low"] == pytest.approx(3.0)
|
||
assert v["target_removal_1"] == pytest.approx(3400 / 1700)
|
||
assert v["stringency"] == pytest.approx(3.0)
|
||
assert v["nsb_green"] == pytest.approx((bkg_grn + 3000) / 700) # max over all NSB beads
|
||
assert v["nsb_red"] == pytest.approx((bkg_red + 3000) / 400) # NSB red uses bkg_red
|
||
assert v["np_green"] == pytest.approx(4000 / 120)
|
||
assert v["np_red"] == pytest.approx(5000 / 300)
|
||
assert not q.missing and not q.zero_divisor
|
||
assert "target_removal_2" not in v # QC-007: not a genotyping column
|
||
|
||
|
||
def test_control_missing_and_zero_divisor() -> None: # QC-007, QC-008
|
||
beads = [b for b in _beads() if b.canonical != "restoration"]
|
||
g = np.ones(len(beads))
|
||
r = np.ones(len(beads))
|
||
i = next(k for k, b in enumerate(beads) if b.canonical == "string_mm")
|
||
r[i] = 0
|
||
q = compute_control_metrics(beads, g, r)
|
||
assert math.isnan(q.values["restoration"]) and q.missing["restoration"] == ["restoration"]
|
||
assert q.zero_divisor == ["stringency"] and q.values["stringency"] == math.inf # QC-008: signal / 0 passes
|
||
r[i - 1] = 0
|
||
q = compute_control_metrics(beads, g, r)
|
||
assert math.isnan(q.values["stringency"]) # 0 / 0 is NaN (flagged)
|
||
|
||
|
||
def test_absent_addresses_are_skipped() -> None: # QC-008, K-021: two bead-pool addresses, one on the chip
|
||
beads = beads_from_definitions(
|
||
[
|
||
{"address": "5:105", "category": "Staining", "color": "Green", "name": "Biotin (High)"},
|
||
{"address": "6:106", "category": "Staining", "color": "Blue", "name": "Biotin (Bgnd)"},
|
||
]
|
||
)
|
||
g = np.array([np.nan, 900.0, 100.0, np.nan])
|
||
q = compute_control_metrics(beads, g, np.zeros(4))
|
||
assert q.values["staining_green"] == pytest.approx(9.0)
|
||
q = compute_control_metrics(beads, np.array([np.nan, 900.0, np.nan, np.nan]), np.zeros(4))
|
||
assert q.missing["staining_green"] == ["biotin_bkg"]
|
||
|
||
|
||
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
|
||
def test_builtin_and_precedence(tmp_path: Path) -> None: # QC-023
|
||
th = resolve_thresholds()
|
||
assert th.get("callRate").value == 0.98 and th.get("logRDev").value == 0.20
|
||
assert not th.get("stainingGreen").enabled
|
||
cfg = tmp_path / "t.yaml"
|
||
cfg.write_text("qc:\n profile: genotyping\n report:\n callRate: 0.90\n npRed: null\n")
|
||
th = resolve_thresholds(config=cfg, cli={"callRate": "0.95"})
|
||
assert th.get("callRate").value == 0.95 and th.get("callRate").source == "cli"
|
||
assert th.get("npRed").value is None and th.get("npRed").source == "config" # QC-022 disabled
|
||
assert th.get("stringency").value == 2 and th.get("stringency").source == "profile:genotyping"
|
||
|
||
|
||
def test_yaml_spellings_and_aliases(tmp_path: Path) -> None: # QC-020
|
||
cfg = tmp_path / "t.yml"
|
||
cfg.write_text("QC:\n Report:\n tgaControl: 1\n stringency: {threshold: 2, direction: above}\n")
|
||
th = resolve_thresholds(config=cfg)
|
||
assert th.get("assayControl").value == 1
|
||
assert th.get("stringency").direction == "above"
|
||
|
||
|
||
@pytest.mark.parametrize(
|
||
("content", "suffix", "code"),
|
||
[
|
||
(None, ".yaml", "E_QC_CONFIG_MISSING"),
|
||
("qc: {report: {}}", ".txt", "E_QC_CONFIG_EXT"),
|
||
("qc: [unclosed", ".yaml", "E_QC_CONFIG_PARSE"),
|
||
("report: {callRate: 1}", ".yaml", "E_QC_CONFIG_STRUCTURE"),
|
||
("qc: {report: {bogus: 1}}", ".yaml", "E_QC_THRESHOLD_KEY"),
|
||
("qc: {report: {callRate: high}}", ".yaml", "E_QC_THRESHOLD_VALUE"),
|
||
("qc: {report: {callRate: 1.5}}", ".yaml", "E_QC_THRESHOLD_RANGE"),
|
||
("qc: {report: {npRed: -1}}", ".yaml", "E_QC_THRESHOLD_RANGE"),
|
||
("qc: {profile: nope, report: {}}", ".yaml", "E_QC_PROFILE"),
|
||
],
|
||
)
|
||
def test_config_validation(tmp_path: Path, content: str | None, suffix: str, code: str) -> None: # QC-021
|
||
p = tmp_path / f"c{suffix}"
|
||
if content is not None:
|
||
p.write_text(content)
|
||
with pytest.raises(SimurgError) as e:
|
||
resolve_thresholds(config=p)
|
||
assert e.value.code == code
|
||
|
||
|
||
def test_profiles() -> None: # QC-024, QC-025
|
||
assert profile_values("genotyping_example")["hybGreenHighMed"] == 1.25
|
||
assert profile_values("methylation_epic+ffpe")["restoration"] == 1
|
||
assert profile_values("methylation_msa")["bisulfiteConversionIIBg"] == 0.5
|
||
assert auto_profile("genotyping", "EX") == "genotyping_ex"
|
||
assert auto_profile("genotyping", "LCG") == "genotyping"
|
||
assert auto_profile("methylation", manifest_name="MSA-48v1") == "methylation_msa"
|
||
with pytest.raises(SimurgError):
|
||
parse_cli_pairs(["callRate"])
|
||
|
||
|
||
# ------------------------------------------------------------------------------------------------ status
|
||
def test_status_axes() -> None: # QC-030, QC-031
|
||
th = resolve_thresholds(profile="genotyping")
|
||
ok = {"status": "OK", "autosomal_call_rate": 0.98, "logr_dev": 0.2} # equality passes
|
||
assert functional_status(ok, th) == "PASS"
|
||
assert functional_status({**ok, "logr_dev": 0.21}, th) == "FAIL"
|
||
assert functional_status({**ok, "logr_dev": float("nan")}, th) == "FAIL"
|
||
assert functional_status({**ok, "status": "FAILED"}, th) == "FAIL"
|
||
assert functional_status(None, th) == "NA"
|
||
vals = {"np_red": 4.0, "stringency": 1.0, "staining_green": 9.0, "restoration": 0.1}
|
||
st, codes = control_status(vals, [], th)
|
||
assert st == "FLAGGED" and codes == ["stringency", "npRed"] # QC-004 order; restoration disabled
|
||
st, codes = control_status({**vals, "np_red": 6, "stringency": 3}, ["np_red"], th)
|
||
assert st == "CLEAR"
|
||
|
||
|
||
def test_flags() -> None: # QC-033
|
||
th = resolve_thresholds()
|
||
row = {
|
||
"status": "OK",
|
||
"autosomal_call_rate": 0.9,
|
||
"logr_dev": 0.3,
|
||
"sex_check": "mismatch",
|
||
"contamination": 0.05,
|
||
"TGA_norm_r": 0.5,
|
||
}
|
||
f = sample_flags(row, th, het_z=4.0, probe_columns=["TGA_norm_r"], duplicates=["x"], unexpected_relatedness=True)
|
||
assert f == [
|
||
"low_call_rate",
|
||
"high_lrr_sd",
|
||
"sex_mismatch",
|
||
"het_outlier",
|
||
"contaminated",
|
||
"assay_control_fail",
|
||
"duplicate_of:x",
|
||
"unexpected_relatedness",
|
||
]
|
||
z = het_zscores([0.1, 0.2, 0.3, float("nan")])
|
||
assert z[1] == pytest.approx(0.0, abs=1e-12) and math.isnan(z[3])
|
||
|
||
|
||
# ------------------------------------------------------------------------------------------------ contamination
|
||
def test_ksmall_matches_selection() -> None:
|
||
rng = np.random.default_rng(0)
|
||
for n in range(1, 40):
|
||
a = rng.random(n).tolist()
|
||
assert ks_ksmall(list(a), n // 2) == sorted(a)[n // 2]
|
||
|
||
|
||
def test_median_modes() -> None:
|
||
b = np.array([[0.1, 0.3, 0.2, 0.4, 0.9]], dtype=np.float32)
|
||
m = np.array([[True, True, True, True, False]])
|
||
assert _median_exact(b, m)[0] == np.float32((np.float32(0.2) + np.float32(0.3)) * np.float32(0.5))
|
||
assert np.isfinite(_median_reference(b[0][m[0]]))
|
||
|
||
|
||
def _simulate(alpha: float, n: int = 60, m: int = 20000, seed: int = 1): # type: ignore[no-untyped-def]
|
||
"""Genotypes and BAF of n samples; sample 0 is contaminated with fraction alpha."""
|
||
rng = np.random.default_rng(seed)
|
||
p = rng.uniform(0.1, 0.9, m)
|
||
g = rng.binomial(2, p[:, None], (m, n))
|
||
baf = g / 2 + rng.normal(0, 0.02, (m, n))
|
||
other = rng.binomial(2, p) / 2
|
||
baf[:, 0] = (1 - alpha) * g[:, 0] / 2 + alpha * other + rng.normal(0, 0.02, m)
|
||
return (g + 1).astype(np.uint8), baf.astype(np.float32)
|
||
|
||
|
||
@pytest.mark.parametrize("median", ["reference", "exact"])
|
||
def test_baf_regress_recovers_alpha(median: str) -> None: # plan §9.3
|
||
gt, baf = _simulate(0.1)
|
||
res = baf_regress(gt, baf, median=median)
|
||
assert res.estimate[0] == pytest.approx(0.1, abs=0.015)
|
||
assert np.all(np.abs(res.estimate[1:]) < 0.01)
|
||
assert res.n_sites > 0 and res.n_hom[0] > 0
|
||
|
||
|
||
def test_af_threshold_is_float32() -> None:
|
||
gt = np.array([[1] * 9 + [3]], dtype=np.uint8) # B frequency 0.1 exactly → excluded (0.1f > 0.1)
|
||
assert b_allele_frequency(gt)[0] == 0.1
|
||
assert baf_regress(gt, np.zeros_like(gt, dtype=np.float32)).n_sites == 0
|
||
|
||
|
||
# ------------------------------------------------------------------------------------------------ relatedness
|
||
def test_king_duplicate_parent_unrelated() -> None:
|
||
rng = np.random.default_rng(3)
|
||
m = 20000
|
||
p = rng.uniform(0.05, 0.95, m)
|
||
hap = lambda: (rng.random(m) < p).astype(np.int8) # noqa: E731
|
||
f1, f2, m1, m2 = hap(), hap(), hap(), hap()
|
||
father, mother = f1 + f2, m1 + m2
|
||
child = np.where(rng.random(m) < 0.5, f1, f2) + np.where(rng.random(m) < 0.5, m1, m2)
|
||
unrel = hap() + hap()
|
||
G = np.stack([father, mother, child, child.copy(), unrel], axis=1).astype(np.uint8) + 1
|
||
ps = pair_stats(G)
|
||
assert ps.kinship[2, 3] == pytest.approx(0.5, abs=0.01)
|
||
assert ps.kinship[0, 2] == pytest.approx(0.25, abs=0.03)
|
||
assert abs(ps.kinship[0, 4]) < 0.03 and abs(ps.kinship[0, 1]) < 0.03
|
||
assert ps.concordance[2, 3] == 1.0
|
||
assert degree(0.5) == "duplicate" and degree(0.25) == "first" and degree(0.0) == "unrelated"
|
||
e, n = mendel_errors(G[:, 2], G[:, 0], G[:, 1])
|
||
assert e == 0 and n == m
|
||
e, _ = mendel_errors(G[:, 4], G[:, 0], G[:, 1])
|
||
assert e > 1000
|
||
|
||
|
||
def test_mendel_rules() -> None:
|
||
c = np.array([1, 2, 3, 2, 1, 0], np.uint8)
|
||
f = np.array([1, 1, 1, 3, 3, 1], np.uint8)
|
||
mo = np.array([1, 1, 3, 3, 2, 1], np.uint8)
|
||
assert mendel_errors(c, f, mo) == (4, 5) # AA×AA→AB, AA×BB→BB, BB×BB→AB, BB×AB→AA; no-call skipped
|
||
assert mendel_errors(c, f, None) == (2, 5) # duo: opposite homozygotes only (BB/AA, AA/BB)
|
||
|
||
|
||
def test_select_sites() -> None:
|
||
counts = np.array([[0, 50, 0, 50], [10, 45, 0, 45], [0, 100, 0, 0], [0, 99, 1, 0]], np.int64)
|
||
idx = select_sites(counts, np.array([True, True, True, False]))
|
||
assert idx.tolist() == [0] # 1: call rate 0.9; 2: monomorphic; 3: not usable
|
||
|
||
|
||
# ------------------------------------------------------------------------------------------------ locus summary
|
||
def test_gs_percentile_matches_scalar_reference() -> None:
|
||
def ref(scores: list[float], p: int) -> float: # BeadArrayFiles ScoreStatistics.percentile
|
||
n = len(scores)
|
||
if n == 0:
|
||
return math.nan
|
||
idx = int(n * p / 100)
|
||
if n * p / 100.0 - idx < 0.5:
|
||
idx -= 1
|
||
if idx < 0:
|
||
return scores[0]
|
||
if idx >= n - 1:
|
||
return scores[-1]
|
||
x1, x2 = 100 * (idx + 0.5) / float(n), 100 * (idx + 1 + 0.5) / float(n)
|
||
y1, y2 = float(scores[idx]), float(scores[idx + 1])
|
||
return y1 + (y2 - y1) / (x2 - x1) * (p - x1)
|
||
|
||
rng = np.random.default_rng(5)
|
||
cols = [np.sort(rng.random(k).astype(np.float32)) for k in range(0, 25)]
|
||
mat = np.full((25, len(cols)), np.inf, np.float32)
|
||
for j, c in enumerate(cols):
|
||
mat[: len(c), j] = c
|
||
counts = np.array([len(c) for c in cols])
|
||
for p in (10, 50):
|
||
got = gs_percentile_columns(mat, counts, p)
|
||
want = [ref(c.tolist(), p) for c in cols]
|
||
np.testing.assert_array_equal(got, np.array(want))
|
||
|
||
|
||
def test_summarize_block_edge_cases() -> None:
|
||
g = np.array([[0, 1, 1, 2], [0, 1, 3, 2], [0, 1, 2, 2]], np.uint8) # samples × loci
|
||
s = np.full(g.shape, 0.8, np.float32)
|
||
st = summarize_block(g, s)
|
||
assert st["chi_test_p100"][0] == 0 and st["het_excess"][0] == 0 # no calls
|
||
assert st["chi_test_p100"][1] == 1 and st["het_excess"][1] == 0 # monomorphic
|
||
assert st["het_excess"][3] > 0.9 and st["chi_test_p100"][3] < 1e-10 # all AB: not monomorphic in the reference
|
||
assert np.isnan(st["gc50"][0]) and st["gc50"][1] == pytest.approx(0.8)
|
||
assert st["minor_freq"][2] == pytest.approx(0.5)
|
||
|
||
|
||
# ------------------------------------------------------------------------------------------------ suggest
|
||
def test_suggest_thresholds() -> None: # QC-026
|
||
rng = np.random.default_rng(2)
|
||
vals = {"npRed": list(np.exp(rng.normal(np.log(10), 0.1, 200))), "logRDev": [0.1] * 10}
|
||
sugs = {s.key: s for s in suggest(vals, resolve_thresholds(), k=4)}
|
||
assert sugs["npRed"].value == pytest.approx(10 * math.exp(-0.4), rel=0.1)
|
||
assert sugs["logRDev"].value is None and "fewer" in sugs["logRDev"].note
|