P5a: QC — control metrics, thresholds, status, contamination, relatedness, locus summary
- qc/controls: canonical control names (QC-002), bkg/offset (QC-003), 13 genotyping metrics (QC-004/005), missing/zero-divisor handling (QC-007/008); controls_raw.tsv / controls_qc.tsv (QC-009/011) incl. compat CSVs - qc/thresholds + profiles: YAML schema, validation codes, precedence, shipped profiles, auto (QC-020..025) - qc/status: functional PASS/FAIL/NA, control CLEAR/FLAGGED, flags (QC-030..033) - qc/contamination: BAF regression; bit-compatible with gtc2vcf BAFregress (36/36, K-017), exact mode too - qc/relatedness: KING-robust, duplicates, Mendelian errors, pedigree from sheet (IN-018) - qc/locus_summary: OUT-072, bit-identical to BeadArrayFiles except ChiTest (K-018) - qc/suggest: suggest-thresholds (QC-026) - pipeline/qc: disk-backed cohort store; CLI qc controls|samples|contamination|relatedness|loci|suggest-thresholds - genotype: control tables (OUT-070), het_rate, QC columns, cohort QC - tests: unit + oracle (BAFregress plugin, BeadArrayFiles locus summary)
This commit is contained in:
kemal committed
1 parent
d11226ca31
commit
a8caa26aa4
28 files changed
+3355
-22
No files matched your search
@@ -300,6 +300,18 @@ def genotype_cmd(
|
||||
str | None, typer.Option("--autocall-date", help="Fixed value for GTC field 201 (reproducible output).")
|
||||
] = None,
|
||||
threads: Annotated[int, typer.Option("--threads", min=0, max=1024, help="Worker processes (0 = auto).")] = 0,
|
||||
controls: Annotated[
|
||||
bool, typer.Option("--controls/--no-controls", help="Write controls_raw.tsv and controls_qc.tsv.")
|
||||
] = True,
|
||||
cohort_qc: Annotated[
|
||||
bool, typer.Option("--cohort-qc/--no-cohort-qc", help="Contamination and relatedness after genotyping.")
|
||||
] = True,
|
||||
qc_config: Annotated[Path | None, typer.Option("--qc-config", help="QC threshold YAML (QC-020).")] = None,
|
||||
qc_profile: Annotated[str | None, typer.Option("--qc-profile", help="QC profile name or 'auto'.")] = None,
|
||||
threshold_call_rate: Annotated[float | None, typer.Option("--threshold-call-rate")] = None,
|
||||
threshold_lrr_sd: Annotated[float | None, typer.Option("--threshold-lrr-sd")] = None,
|
||||
threshold: Annotated[list[str] | None, typer.Option("--threshold", help="KEY=VALUE (repeatable).")] = None,
|
||||
bg_offset: Annotated[float | None, typer.Option("--bg-offset", help="Control background offset (3000).")] = None,
|
||||
force: Annotated[bool, typer.Option("--force", help="Replace existing outputs.")] = False,
|
||||
layout: Annotated[str, typer.Option("--layout", help="nested|flat")] = "nested",
|
||||
on_duplicate_idat: Annotated[str, typer.Option("--on-duplicate-idat", help="error|first|newest")] = "error",
|
||||
@@ -345,10 +357,276 @@ def genotype_cmd(
|
||||
sex_version=sex_version,
|
||||
gtc_version=gtc_version,
|
||||
autocall_date=autocall_date,
|
||||
controls=controls,
|
||||
cohort_qc=cohort_qc,
|
||||
qc_config=qc_config,
|
||||
qc_profile=qc_profile,
|
||||
threshold_call_rate=threshold_call_rate,
|
||||
threshold_lrr_sd=threshold_lrr_sd,
|
||||
thresholds=list(threshold or []),
|
||||
bg_offset=bg_offset,
|
||||
argv=sys.argv[1:],
|
||||
)
|
||||
run_guarded(rep, lambda: run_genotype(opts, rep))
|
||||
|
||||
|
||||
# --------------------------------------------------------------------------------------------------------
|
||||
qc_app = typer.Typer(help="Quality control: controls, samples, contamination, relatedness, loci, thresholds.")
|
||||
app.add_typer(qc_app, name="qc")
|
||||
|
||||
|
||||
def _reporter(command: str, out: Path | None, log_json: bool, debug: bool, quiet: bool) -> Reporter:
|
||||
return Reporter(
|
||||
command=command,
|
||||
level="debug" if debug else "info",
|
||||
json_stderr=log_json,
|
||||
debug=debug,
|
||||
quiet=quiet,
|
||||
log_file=(out / "logs" / f"{command.replace(' ', '_')}.log") if out is not None else None,
|
||||
)
|
||||
|
||||
|
||||
LogJson = Annotated[bool, typer.Option("--log-json", "--json-log", help="JSON Lines log on stderr.")]
|
||||
Debug = Annotated[bool, typer.Option("--debug", help="Debug logging with tracebacks.")]
|
||||
Quiet = Annotated[bool, typer.Option("--quiet", "--no-progress", help="Suppress info messages.")]
|
||||
|
||||
|
||||
@qc_app.command("controls")
|
||||
def qc_controls_cmd(
|
||||
out: Annotated[Path, typer.Option("--out", help="Output folder (required).")],
|
||||
array_type: Annotated[str, typer.Option("--array-type", help="genotyping|methylation (required).")],
|
||||
bpm: Annotated[Path | None, typer.Option("--bpm", help="BPM manifest (control definitions).")] = None,
|
||||
csv: Annotated[Path | None, typer.Option("--csv", help="CSV manifest (control definitions).")] = None,
|
||||
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,
|
||||
bg_offset: Annotated[float, typer.Option("--bg-offset", help="Background offset x (QC-003).")] = 3000.0,
|
||||
compat: Annotated[bool, typer.Option("--compat-column-names", help="Reference-style CSV names.")] = False,
|
||||
force: Annotated[bool, typer.Option("--force", help="Replace existing control tables.")] = False,
|
||||
log_json: LogJson = False,
|
||||
debug: Debug = False,
|
||||
quiet: Quiet = False,
|
||||
) -> None:
|
||||
"""Control-probe QC from IDATs without genotyping (QC-012)."""
|
||||
from simurg.pipeline.qc import QcControlsOptions, run_qc_controls
|
||||
|
||||
if array_type not in ("genotyping", "methylation"):
|
||||
raise typer.BadParameter("--array-type must be genotyping or methylation")
|
||||
rep = _reporter("qc controls", out, log_json, debug, quiet)
|
||||
opts = QcControlsOptions(
|
||||
out=out,
|
||||
bpm=bpm,
|
||||
csv=csv,
|
||||
idats=idats,
|
||||
sheet=sheet,
|
||||
array_type=array_type,
|
||||
bg_offset=bg_offset,
|
||||
compat=compat,
|
||||
force=force,
|
||||
argv=sys.argv[1:],
|
||||
)
|
||||
run_guarded(rep, lambda: run_qc_controls(opts, rep))
|
||||
|
||||
|
||||
def _qc_cohort(
|
||||
command: str,
|
||||
out: Path,
|
||||
bpm: Path,
|
||||
gtcs: Path | None,
|
||||
sheet: Path | None,
|
||||
egt: Path | None,
|
||||
threads: int,
|
||||
force: bool,
|
||||
log_json: bool,
|
||||
debug: bool,
|
||||
quiet: bool,
|
||||
*,
|
||||
all_pairs: bool = False,
|
||||
pedigree_from_sheet: bool = True,
|
||||
median: str = "reference",
|
||||
config: Path | None = None,
|
||||
profile: str | None = None,
|
||||
threshold_call_rate: float | None = None,
|
||||
threshold_lrr_sd: float | None = None,
|
||||
threshold: list[str] | None = None,
|
||||
) -> None:
|
||||
from simurg.pipeline.qc import QcCohortOptions, run_qc_cohort
|
||||
from simurg.qc.thresholds import parse_cli_pairs, resolve_thresholds
|
||||
|
||||
if median not in ("reference", "exact"):
|
||||
raise typer.BadParameter("--median must be reference or exact")
|
||||
rep = _reporter(f"qc {command}", out, log_json, debug, quiet)
|
||||
|
||||
def body() -> int:
|
||||
cli: dict[str, Any] = parse_cli_pairs(list(threshold or []))
|
||||
if threshold_call_rate is not None:
|
||||
cli["callRate"] = threshold_call_rate
|
||||
if threshold_lrr_sd is not None:
|
||||
cli["logRDev"] = threshold_lrr_sd
|
||||
th = resolve_thresholds(config=config, profile=profile, cli=cli)
|
||||
opts = QcCohortOptions(
|
||||
command=command,
|
||||
out=out,
|
||||
bpm=bpm,
|
||||
gtcs=gtcs,
|
||||
sheet=sheet,
|
||||
egt=egt,
|
||||
threads=threads,
|
||||
force=force,
|
||||
all_pairs=all_pairs,
|
||||
pedigree_from_sheet=pedigree_from_sheet,
|
||||
median=median,
|
||||
thresholds=th,
|
||||
argv=sys.argv[1:],
|
||||
)
|
||||
return run_qc_cohort(opts, rep)
|
||||
|
||||
run_guarded(rep, body)
|
||||
|
||||
|
||||
GtcsOpt = Annotated[Path | None, typer.Option("--gtcs", help="Folder searched recursively for GTC files.")]
|
||||
SheetOpt = Annotated[Path | None, typer.Option("--sheet", help="Sample sheet (CSV).")]
|
||||
BpmOpt = Annotated[Path, typer.Option("--bpm", help="BPM manifest (required).")]
|
||||
OutOpt = Annotated[Path, typer.Option("--out", help="Output folder (required).")]
|
||||
ThreadsOpt = Annotated[int, typer.Option("--threads", min=0, max=1024, help="Worker processes (0 = auto).")]
|
||||
ForceOpt = Annotated[bool, typer.Option("--force", help="Replace existing outputs.")]
|
||||
MedianOpt = Annotated[str, typer.Option("--median", help="reference (default, bit-compatible) or exact.")]
|
||||
|
||||
|
||||
@qc_app.command("samples")
|
||||
def qc_samples_cmd(
|
||||
out: OutOpt,
|
||||
bpm: BpmOpt,
|
||||
gtcs: GtcsOpt = None,
|
||||
sheet: SheetOpt = None,
|
||||
config: Annotated[Path | None, typer.Option("--config", help="QC threshold YAML.")] = None,
|
||||
profile: Annotated[str | None, typer.Option("--profile", help="QC profile name or 'auto'.")] = None,
|
||||
threshold_call_rate: Annotated[float | None, typer.Option("--threshold-call-rate")] = None,
|
||||
threshold_lrr_sd: Annotated[float | None, typer.Option("--threshold-lrr-sd")] = None,
|
||||
threshold: Annotated[list[str] | None, typer.Option("--threshold", help="KEY=VALUE (repeatable).")] = None,
|
||||
all_pairs: Annotated[bool, typer.Option("--all-pairs", help="Write every sample pair.")] = False,
|
||||
median: MedianOpt = "reference",
|
||||
threads: ThreadsOpt = 0,
|
||||
force: ForceOpt = False,
|
||||
log_json: LogJson = False,
|
||||
debug: Debug = False,
|
||||
quiet: Quiet = False,
|
||||
) -> None:
|
||||
"""Sample QC from GTCs: het rate, contamination, relatedness, flags (sample_qc.tsv)."""
|
||||
_qc_cohort(
|
||||
"samples",
|
||||
out,
|
||||
bpm,
|
||||
gtcs,
|
||||
sheet,
|
||||
None,
|
||||
threads,
|
||||
force,
|
||||
log_json,
|
||||
debug,
|
||||
quiet,
|
||||
all_pairs=all_pairs,
|
||||
median=median,
|
||||
config=config,
|
||||
profile=profile,
|
||||
threshold_call_rate=threshold_call_rate,
|
||||
threshold_lrr_sd=threshold_lrr_sd,
|
||||
threshold=threshold,
|
||||
)
|
||||
|
||||
|
||||
@qc_app.command("contamination")
|
||||
def qc_contamination_cmd(
|
||||
out: OutOpt,
|
||||
bpm: BpmOpt,
|
||||
gtcs: GtcsOpt = None,
|
||||
sheet: SheetOpt = None,
|
||||
median: MedianOpt = "reference",
|
||||
threads: ThreadsOpt = 0,
|
||||
force: ForceOpt = False,
|
||||
log_json: LogJson = False,
|
||||
debug: Debug = False,
|
||||
quiet: Quiet = False,
|
||||
) -> None:
|
||||
"""BAF-regression contamination estimate per sample (contamination.tsv)."""
|
||||
_qc_cohort("contamination", out, bpm, gtcs, sheet, None, threads, force, log_json, debug, quiet, median=median)
|
||||
|
||||
|
||||
@qc_app.command("relatedness")
|
||||
def qc_relatedness_cmd(
|
||||
out: OutOpt,
|
||||
bpm: BpmOpt,
|
||||
gtcs: GtcsOpt = None,
|
||||
sheet: SheetOpt = None,
|
||||
pedigree_from_sheet: Annotated[
|
||||
bool, typer.Option("--pedigree-from-sheet/--no-pedigree", help="Use Parent1/Parent2/Replicate.")
|
||||
] = True,
|
||||
all_pairs: Annotated[bool, typer.Option("--all-pairs", help="Write every sample pair.")] = False,
|
||||
threads: ThreadsOpt = 0,
|
||||
force: ForceOpt = False,
|
||||
log_json: LogJson = False,
|
||||
debug: Debug = False,
|
||||
quiet: Quiet = False,
|
||||
) -> None:
|
||||
"""KING-robust kinship, duplicates and Mendelian errors (relatedness.tsv, mendel.tsv)."""
|
||||
_qc_cohort(
|
||||
"relatedness",
|
||||
out,
|
||||
bpm,
|
||||
gtcs,
|
||||
sheet,
|
||||
None,
|
||||
threads,
|
||||
force,
|
||||
log_json,
|
||||
debug,
|
||||
quiet,
|
||||
all_pairs=all_pairs,
|
||||
pedigree_from_sheet=pedigree_from_sheet,
|
||||
)
|
||||
|
||||
|
||||
@qc_app.command("loci")
|
||||
def qc_loci_cmd(
|
||||
out: OutOpt,
|
||||
bpm: BpmOpt,
|
||||
egt: Annotated[Path, typer.Option("--egt", help="EGT cluster file (required).")],
|
||||
gtcs: GtcsOpt = None,
|
||||
sheet: SheetOpt = None,
|
||||
threads: ThreadsOpt = 0,
|
||||
force: ForceOpt = False,
|
||||
log_json: LogJson = False,
|
||||
debug: Debug = False,
|
||||
quiet: Quiet = False,
|
||||
) -> None:
|
||||
"""Per-locus summary statistics across samples (locus_summary.tsv, OUT-072)."""
|
||||
_qc_cohort("loci", out, bpm, gtcs, sheet, egt, threads, force, log_json, debug, quiet)
|
||||
|
||||
|
||||
@qc_app.command("suggest-thresholds")
|
||||
def qc_suggest_cmd(
|
||||
data: Annotated[list[Path], typer.Option("--data", help="Dataset folder(s) with controls_qc.tsv.")],
|
||||
out: Annotated[Path, typer.Option("--out", help="Output YAML file (required).")],
|
||||
k: Annotated[float, typer.Option("--k", min=0.0, help="Spread multiplier.")] = 4.0,
|
||||
min_samples: Annotated[int, typer.Option("--min-samples", min=1, help="Minimum samples per metric.")] = 30,
|
||||
config: Annotated[Path | None, typer.Option("--config", help="Current thresholds (YAML).")] = None,
|
||||
profile: Annotated[str | None, typer.Option("--profile", help="Current profile.")] = None,
|
||||
log_json: LogJson = False,
|
||||
debug: Debug = False,
|
||||
quiet: Quiet = False,
|
||||
) -> None:
|
||||
"""Data-driven QC threshold suggestion (QC-026)."""
|
||||
from simurg.pipeline.qc import run_suggest_thresholds
|
||||
from simurg.qc.thresholds import resolve_thresholds
|
||||
|
||||
rep = _reporter("qc suggest-thresholds", None, log_json, debug, quiet)
|
||||
folders = [Path(p) for d in data for p in str(d).split(",") if p]
|
||||
|
||||
def body() -> int:
|
||||
th = resolve_thresholds(config=config, profile=profile)
|
||||
return run_suggest_thresholds(folders, out, th, rep, k=k, min_samples=min_samples)
|
||||
|
||||
run_guarded(rep, body)
|
||||
|
||||
|
||||
if __name__ == "__main__": # pragma: no cover
|
||||
app()
|
||||
@@ -89,6 +89,30 @@ CATALOGUE: dict[str, MessageDef] = dict(
|
||||
_m("E_EGT_LOAD", "error", "Cluster (EGT) file cannot be loaded", "IN-044"),
|
||||
_m("W_EGT_NAME", "warning", "Cluster file name is not portable", "IN-045"),
|
||||
_m("W_MANIFEST_NO_REFSTRAND", "warning", "Manifest lacks RefStrand annotation", "IN-040"),
|
||||
# --- QC (QC-001 .. QC-064) -----------------------------------------------------------------
|
||||
_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"),
|
||||
_m("E_QC_LABEL_DUPLICATE", "error", "Dataset labels are not unique", "QC-042"),
|
||||
_m("W_QC_LABEL_SANITIZED", "warning", "Dataset label was sanitized", "QC-042"),
|
||||
_m("E_QC_MIXED_TYPES", "error", "Genotyping and methylation datasets cannot be mixed", "QC-043"),
|
||||
_m("W_QC_UNMATCHED_SAMPLES", "warning", "Sample IDs present in one QC file but not the other", "QC-040"),
|
||||
_m("E_QC_FORMAT", "error", "Invalid --output-format value", "QC-063"),
|
||||
_m("E_QC_CONFIG_MISSING", "error", "QC config file does not exist", "QC-021"),
|
||||
_m("E_QC_CONFIG_EXT", "error", "QC config file must be .yaml or .yml", "QC-021"),
|
||||
_m("E_QC_CONFIG_PARSE", "error", "QC config file is not valid YAML", "QC-021"),
|
||||
_m("E_QC_CONFIG_STRUCTURE", "error", "QC config lacks the qc/report (QC/Report) wrapper", "QC-021"),
|
||||
_m("E_QC_THRESHOLD_KEY", "error", "Unknown QC threshold key", "QC-021"),
|
||||
_m("E_QC_THRESHOLD_VALUE", "error", "QC threshold value is not numeric or null", "QC-021"),
|
||||
_m("E_QC_THRESHOLD_RANGE", "error", "QC threshold value out of range", "QC-021"),
|
||||
_m("E_QC_PROFILE", "error", "Unknown QC profile", "QC-024"),
|
||||
_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"),
|
||||
]
|
||||
)
|
||||
|
||||
|
||||
@@ -194,3 +194,19 @@ def resolve_samples(disc: Discovery, sheet: SampleSheet | None) -> tuple[list[Re
|
||||
if extra:
|
||||
notes.append(("I_UNLISTED_IDATS", f"{len(extra)} IDAT pair(s) not listed in the sample sheet are ignored"))
|
||||
return out, notes
|
||||
|
||||
|
||||
def discover_gtcs(root: str | Path) -> dict[str, Path]:
|
||||
"""``*.gtc`` under ``root`` keyed by file-name stem (IN-030); duplicate keys → ``E_DUP_GTC``."""
|
||||
root = Path(root)
|
||||
if not root.exists():
|
||||
raise SimurgError("E_FILE_NOT_FOUND", f"GTC folder not found: {root}", path=str(root))
|
||||
out: dict[str, Path] = {}
|
||||
for p in walk_files(root, (".gtc",)):
|
||||
key = p.name[: -len(".gtc")]
|
||||
if key in out:
|
||||
raise SimurgError("E_DUP_GTC", f"GTC files {out[key]} and {p} resolve to the same sample {key!r}")
|
||||
out[key] = p
|
||||
if not out:
|
||||
raise SimurgError("E_NO_INPUT", f"no GTC files found under {root}")
|
||||
return out
|
||||
+164
-11
@@ -22,6 +22,7 @@ from simurg.core.errors import ExitCode, SimurgError
|
||||
from simurg.core.log import Reporter
|
||||
from simurg.core.provenance import RunManifest
|
||||
from simurg.formats.bpm import read_bpm
|
||||
from simurg.formats.csv_manifest import assay_format
|
||||
from simurg.formats.egt import read_egt
|
||||
from simurg.formats.gtc import write_gtc
|
||||
from simurg.formats.idat import check_pair, read_idat
|
||||
@@ -37,6 +38,24 @@ from simurg.genotype.sample import (
|
||||
)
|
||||
from simurg.genotype.sex import SexParams
|
||||
from simurg.pipeline.inputs import check_paths, resolve_sample_set, sample_dir, sheet_values
|
||||
from simurg.pipeline.qc import (
|
||||
CONTROLS_QC,
|
||||
CONTROLS_RAW,
|
||||
MENDEL_COLUMNS,
|
||||
RELATEDNESS_COLUMNS,
|
||||
ControlSample,
|
||||
cohort_qc,
|
||||
control_sample,
|
||||
default_tmp,
|
||||
load_cohort,
|
||||
load_control_beads,
|
||||
pedigree_from_sheet,
|
||||
write_control_tables,
|
||||
write_tsv,
|
||||
)
|
||||
from simurg.qc.controls import ControlBead
|
||||
from simurg.qc.status import control_status, functional_status, het_zscores, sample_flags
|
||||
from simurg.qc.thresholds import Thresholds, parse_cli_pairs, resolve_thresholds
|
||||
from simurg.signal.normalization.apply import normalize, polar
|
||||
from simurg.signal.raw_xy import LocusMapCache, check_coverage
|
||||
|
||||
@@ -67,7 +86,7 @@ FIXED_COLUMNS = [
|
||||
"het_rate",
|
||||
"contamination",
|
||||
]
|
||||
QC_COLUMNS = ["functional_status", "control_status", "control_flag_codes", "flags"] # filled from phase P5
|
||||
QC_COLUMNS = ["functional_status", "control_status", "control_flag_codes", "flags"] # QC-030/031/033
|
||||
FLOAT_COLUMNS = {"autosomal_call_rate", "call_rate", "logr_dev", "gc10", "gc50", "het_rate", "contamination"}
|
||||
|
||||
|
||||
@@ -90,6 +109,14 @@ class GenotypeOptions:
|
||||
gtc_version: int | None = None
|
||||
autocall_date: str | None = None
|
||||
assay_control_probes: tuple[str, ...] = ("TGA_Ctrl_5716",)
|
||||
controls: bool = True # OUT-070
|
||||
cohort_qc: bool = True # contamination + relatedness after genotyping (QC-064)
|
||||
qc_config: Path | None = None
|
||||
qc_profile: str | None = None
|
||||
threshold_call_rate: float | None = None
|
||||
threshold_lrr_sd: float | None = None
|
||||
thresholds: list[str] = field(default_factory=list)
|
||||
bg_offset: float | None = None
|
||||
argv: list[str] = field(default_factory=list)
|
||||
|
||||
|
||||
@@ -100,6 +127,7 @@ class Job:
|
||||
grn: str
|
||||
red: str
|
||||
gtc_path: str
|
||||
sheet_values: dict[str, str] = field(default_factory=dict)
|
||||
|
||||
|
||||
_W: dict[str, Any] = {}
|
||||
@@ -111,8 +139,12 @@ def _init(
|
||||
auto_1_22: np.ndarray[Any, np.dtype[np.bool_]],
|
||||
probe_idx: dict[str, int],
|
||||
started: _dt.datetime,
|
||||
beads: list[ControlBead] | None = None,
|
||||
bg_offset: float = 3000.0,
|
||||
) -> None:
|
||||
_W.update(ctx=ctx, params=params, auto=auto_1_22, probes=probe_idx, started=started)
|
||||
_W.update(
|
||||
ctx=ctx, params=params, auto=auto_1_22, probes=probe_idx, started=started, beads=beads, bg_offset=bg_offset
|
||||
)
|
||||
|
||||
|
||||
def _process(job: Job) -> dict[str, Any]:
|
||||
@@ -134,6 +166,11 @@ def _process(job: Job) -> dict[str, Any]:
|
||||
write_gtc(g, Path(job.gtc_path))
|
||||
auto = _W["auto"] & ~ctx.manifest.intensity_only
|
||||
n_auto = int(auto.sum())
|
||||
g_auto = g.genotypes[auto]
|
||||
n_called = int((g_auto > 0).sum())
|
||||
row["het_rate"] = float((g_auto == 2).sum() / n_called) if n_called else math.nan
|
||||
if _W.get("beads"):
|
||||
row["control"] = control_sample(job.sample_id, grn, red, _W["beads"], _W["bg_offset"], job.sheet_values)
|
||||
row.update(
|
||||
autosomal_call_rate=float(np.float32(int((g.genotypes[auto] > 0).sum())) / np.float32(n_auto))
|
||||
if n_auto
|
||||
@@ -233,10 +270,22 @@ def run_genotype(opts: GenotypeOptions, rep: Reporter) -> int:
|
||||
continue
|
||||
stem = sset.stems[s.sample_id]
|
||||
gtc = sample_dir(opts.out, stem, opts.layout) / f"{stem}.gtc"
|
||||
jobs.append(Job(s.sample_id, s.sample_name, str(s.pair.grn), str(s.pair.red), str(gtc)))
|
||||
jobs.append(
|
||||
Job(
|
||||
s.sample_id,
|
||||
s.sample_name,
|
||||
str(s.pair.grn),
|
||||
str(s.pair.red),
|
||||
str(gtc),
|
||||
{k: str(v) for k, v in sheet_values(s, sset.sheet).items()},
|
||||
)
|
||||
)
|
||||
|
||||
summary_tsv, summary_json = opts.out / "sample_summary.tsv", opts.out / "sample_summary.json"
|
||||
check_outputs_absent([summary_tsv, summary_json] + [Path(j.gtc_path) for j in jobs], force=opts.force)
|
||||
extra = [opts.out / CONTROLS_RAW, opts.out / CONTROLS_QC] if opts.controls else []
|
||||
if opts.cohort_qc:
|
||||
extra.append(opts.out / "relatedness.tsv")
|
||||
check_outputs_absent([summary_tsv, summary_json, *extra] + [Path(j.gtc_path) for j in jobs], force=opts.force)
|
||||
|
||||
rep.info(f"reading manifest {opts.bpm} and cluster file {opts.egt}")
|
||||
manifest = read_bpm(opts.bpm)
|
||||
@@ -256,6 +305,10 @@ def run_genotype(opts: GenotypeOptions, rep: Reporter) -> int:
|
||||
"W_EGT_NAME", f"cluster file was built for {egt.manifest_name!r}, manifest is {manifest.manifest_name!r}"
|
||||
)
|
||||
clusters = align_clusters(egt, manifest)
|
||||
th = _thresholds(opts, manifest)
|
||||
for rec in th.as_records():
|
||||
rm.set_option(f"qc.{rec['key']}", rec["value"], rec["source"])
|
||||
beads = load_control_beads(manifest, opts.csv, rep) if opts.controls else []
|
||||
ctx = manifest_context(manifest, clusters, cluster_file_name=Path(opts.egt).name)
|
||||
lite = manifest.lite()
|
||||
ctx = replace(ctx, manifest=lite, cache=LocusMapCache(lite)) # workers get no string columns
|
||||
@@ -282,14 +335,14 @@ def run_genotype(opts: GenotypeOptions, rep: Reporter) -> int:
|
||||
rep.info(f"genotyping {len(jobs)} sample(s) with {workers} worker(s)")
|
||||
try:
|
||||
if workers == 1:
|
||||
_init(ctx, params, auto_1_22, probes, started)
|
||||
_init(ctx, params, auto_1_22, probes, started, beads, th.bg_offset)
|
||||
results = [_process(j) for j in jobs]
|
||||
else:
|
||||
with ProcessPoolExecutor(
|
||||
workers,
|
||||
mp_context=mp.get_context("spawn"),
|
||||
initializer=_init,
|
||||
initargs=(ctx, params, auto_1_22, probes, started),
|
||||
initargs=(ctx, params, auto_1_22, probes, started, beads, th.bg_offset),
|
||||
) as ex:
|
||||
results = list(ex.map(_process, jobs))
|
||||
except KeyboardInterrupt:
|
||||
@@ -315,14 +368,11 @@ def run_genotype(opts: GenotypeOptions, rep: Reporter) -> int:
|
||||
rows[r["sample_id"]].update(r)
|
||||
|
||||
probe_cols = [f"{p}_norm_r" for p in probes]
|
||||
_run_qc(opts, rep, rm, th, beads, sset, rows, results, jobs, manifest, probe_cols, params.gtc_version)
|
||||
sheet_cols = [f"sheet:{c}" for c in (sset.sheet.passthrough_columns() if sset.sheet else [])]
|
||||
out_rows = []
|
||||
for s in sset.samples:
|
||||
r = rows[s.sample_id]
|
||||
known = {"M", "F"}
|
||||
est, rep_sex = r.get("sex_est", "U"), r.get("sex_reported", "U")
|
||||
r["sex_check"] = ("match" if est == rep_sex else "mismatch") if est in known and rep_sex in known else "unknown"
|
||||
out_rows.append(r)
|
||||
out_rows.append(rows[s.sample_id])
|
||||
columns = FIXED_COLUMNS + probe_cols + QC_COLUMNS + sheet_cols
|
||||
header = [c.removeprefix("sheet:") for c in columns]
|
||||
lines = ["\t".join(header)] + ["\t".join(_fmt(c, r.get(c, "")) for c in columns) for r in out_rows]
|
||||
@@ -350,3 +400,106 @@ def run_genotype(opts: GenotypeOptions, rep: Reporter) -> int:
|
||||
if ok == len(out_rows):
|
||||
return ExitCode.OK
|
||||
return ExitCode.PARTIAL if ok else ExitCode.ALL_FAILED
|
||||
|
||||
|
||||
def _thresholds(opts: GenotypeOptions, manifest: Any) -> Thresholds:
|
||||
cli: dict[str, Any] = parse_cli_pairs(opts.thresholds)
|
||||
if opts.threshold_call_rate is not None:
|
||||
cli["callRate"] = opts.threshold_call_rate
|
||||
if opts.threshold_lrr_sd is not None:
|
||||
cli["logRDev"] = opts.threshold_lrr_sd
|
||||
if opts.bg_offset is not None:
|
||||
cli["bgOffset"] = opts.bg_offset
|
||||
return resolve_thresholds(
|
||||
config=opts.qc_config,
|
||||
profile=opts.qc_profile,
|
||||
cli=cli,
|
||||
array_type="genotyping",
|
||||
assay_format=assay_format(manifest),
|
||||
manifest_name=manifest.manifest_name,
|
||||
)
|
||||
|
||||
|
||||
def _run_qc(
|
||||
opts: GenotypeOptions,
|
||||
rep: Reporter,
|
||||
rm: RunManifest,
|
||||
th: Thresholds,
|
||||
beads: list[ControlBead],
|
||||
sset: Any,
|
||||
rows: dict[str, dict[str, Any]],
|
||||
results: list[dict[str, Any]],
|
||||
jobs: list[Job],
|
||||
manifest: Any,
|
||||
probe_cols: list[str],
|
||||
gtc_version: int,
|
||||
) -> None:
|
||||
"""Control tables (OUT-070), het rate, cohort QC and the QC-030/031/033 columns of the summary."""
|
||||
controls: dict[str, ControlSample] = {r["sample_id"]: r.pop("control") for r in results if "control" in r}
|
||||
for r in rows.values(): # results were merged into rows already; drop the heavy objects there too
|
||||
r.pop("control", None)
|
||||
if opts.controls and beads:
|
||||
assay = {
|
||||
sid: min((float(rows[sid].get(c, math.nan)) for c in probe_cols), default=math.nan) for sid in controls
|
||||
}
|
||||
ordered = [controls[s.sample_id] for s in sset.samples if s.sample_id in controls]
|
||||
for p in write_control_tables(
|
||||
opts.out, ordered, beads, th.bg_offset, rep, assay_control=assay if probe_cols else None
|
||||
):
|
||||
rm.add_output(p, "simurg.controls/1")
|
||||
ids = [s.sample_id for s in sset.samples]
|
||||
het_z = dict(zip(ids, het_zscores([rows[i].get("het_rate", math.nan) for i in ids]), strict=True))
|
||||
dups: dict[str, list[str]] = {}
|
||||
unexpected: set[str] = set()
|
||||
ok_jobs = [j for j in jobs if rows[j.sample_id].get("status") == "OK"]
|
||||
if opts.cohort_qc and len(ok_jobs) >= 2:
|
||||
tmp = default_tmp(opts.out)
|
||||
rep.info(f"cohort QC (contamination, relatedness) on {len(ok_jobs)} sample(s)")
|
||||
store, infos = load_cohort(
|
||||
[(j.sample_id, j.gtc_path) for j in ok_jobs],
|
||||
manifest.num_loci,
|
||||
tmp,
|
||||
baf=gtc_version >= 4,
|
||||
scores=False,
|
||||
workers=max(1, opts.threads or (os.cpu_count() or 1)),
|
||||
)
|
||||
try:
|
||||
ped = pedigree_from_sheet(sset.sheet, [j.sample_id for j in ok_jobs], rep)
|
||||
res = cohort_qc(
|
||||
store, infos, manifest, rep, pedigree=ped, workers=max(1, opts.threads or (os.cpu_count() or 1))
|
||||
)
|
||||
finally:
|
||||
store.close()
|
||||
if res.contamination is not None:
|
||||
for i, s in enumerate(infos):
|
||||
rows[s.sample_id]["contamination"] = float(res.contamination.estimate[i])
|
||||
dups, unexpected = res.duplicates, res.unexpected
|
||||
write_tsv(opts.out / "relatedness.tsv", RELATEDNESS_COLUMNS, res.pairs)
|
||||
rm.add_output(opts.out / "relatedness.tsv", "simurg.relatedness/1")
|
||||
if res.mendel:
|
||||
write_tsv(opts.out / "mendel.tsv", MENDEL_COLUMNS, res.mendel)
|
||||
rm.add_output(opts.out / "mendel.tsv", "simurg.mendel/1")
|
||||
for sid in ids:
|
||||
r = rows[sid]
|
||||
known = {"M", "F"}
|
||||
est, rep_sex = r.get("sex_est", "U"), r.get("sex_reported", "U")
|
||||
r["sex_check"] = ("match" if est == rep_sex else "mismatch") if est in known and rep_sex in known else "unknown"
|
||||
if r.get("status") == "MISSING_INPUT":
|
||||
continue
|
||||
r["functional_status"] = functional_status(r, th)
|
||||
if sid in controls:
|
||||
vals = dict(controls[sid].qc.values)
|
||||
if probe_cols:
|
||||
vals["assay_control_norm_r"] = min(float(r.get(c, math.nan)) for c in probe_cols)
|
||||
cs, codes = control_status(vals, controls[sid].qc.na_metrics, th)
|
||||
r["control_status"], r["control_flag_codes"] = cs, ";".join(codes)
|
||||
r["flags"] = ";".join(
|
||||
sample_flags(
|
||||
r,
|
||||
th,
|
||||
het_z=het_z.get(sid, math.nan),
|
||||
probe_columns=probe_cols,
|
||||
duplicates=dups.get(sid, []),
|
||||
unexpected_relatedness=sid in unexpected,
|
||||
)
|
||||
)
|
||||
@@ -0,0 +1,915 @@
|
||||
"""QC orchestration and I/O (spec QC-001 … QC-064, CLI-024, CLI-025, OUT-070, OUT-072; plan §9).
|
||||
|
||||
The science lives in ``simurg.qc``; this module reads IDATs/GTCs, streams cohort matrices through disk-backed
|
||||
arrays (so thousands of samples fit in a small memory budget) and writes the QC tables.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import datetime as _dt
|
||||
import math
|
||||
import multiprocessing as mp
|
||||
import os
|
||||
import shutil
|
||||
import tempfile
|
||||
from collections.abc import Iterator, Mapping, Sequence
|
||||
from concurrent.futures import ProcessPoolExecutor
|
||||
from dataclasses import dataclass, field
|
||||
from pathlib import Path
|
||||
from typing import IO, Any, cast
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
from simurg.core.atomic import atomic_write_text, 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.bpm import Manifest
|
||||
from simurg.formats.idat import Idat
|
||||
from simurg.formats.samplesheet import SampleSheet, fold
|
||||
from simurg.genotype.baf_lrr import chrom_core
|
||||
from simurg.qc.contamination import ContaminationResult, baf_regress
|
||||
from simurg.qc.controls import (
|
||||
ASSAY_CONTROL_COLUMN,
|
||||
COMPAT_METRIC_NAMES,
|
||||
ControlBead,
|
||||
ControlQC,
|
||||
beads_from_definitions,
|
||||
compare_definitions,
|
||||
compute_control_metrics,
|
||||
control_intensities,
|
||||
metrics_for,
|
||||
)
|
||||
from simurg.qc.relatedness import DEGREE_CUTS, DUPLICATE_CONCORDANCE, degree, mendel_errors, pair_stats, select_sites
|
||||
from simurg.qc.status import het_zscores
|
||||
from simurg.qc.thresholds import Thresholds
|
||||
|
||||
CONTROLS_RAW = "controls_raw.tsv"
|
||||
CONTROLS_QC = "controls_qc.tsv"
|
||||
COMPAT_RAW = "controls.raw_metrics.csv"
|
||||
COMPAT_QC = "controls.qc_metrics.csv"
|
||||
MIN_CONTAMINATION_SAMPLES = 10
|
||||
|
||||
|
||||
def fmt(v: Any, digits: int = 8) -> str:
|
||||
if v is None or v == "":
|
||||
return ""
|
||||
if isinstance(v, float | np.floating):
|
||||
f = float(v)
|
||||
return "NaN" if math.isnan(f) else f"{f:.{digits}f}"
|
||||
return str(v).replace("\t", " ").replace("\n", " ")
|
||||
|
||||
|
||||
def write_tsv(path: Path, columns: Sequence[str], rows: Sequence[Mapping[str, Any]], sep: str = "\t") -> None:
|
||||
lines = [sep.join(columns)] + [sep.join(fmt(r.get(c, "")) for c in columns) for r in rows]
|
||||
atomic_write_text(path, "\n".join(lines) + "\n")
|
||||
|
||||
|
||||
def autosome_mask(chroms: Sequence[str]) -> npt.NDArray[np.bool_]:
|
||||
"""Chromosomes 1–22 only (het rate, relatedness, contamination sites)."""
|
||||
auto = {str(i) for i in range(1, 23)}
|
||||
return np.array([chrom_core(c) in auto for c in chroms], dtype=np.bool_)
|
||||
|
||||
|
||||
# ================================================================================================ controls
|
||||
def load_control_beads(manifest: Manifest | None, csv_path: Path | None, rep: Reporter) -> list[ControlBead]:
|
||||
"""QC-001: BPM control config wins over the CSV ``[Controls]`` section; differences give one warning."""
|
||||
bpm_beads = beads_from_definitions(manifest.controls) if manifest is not None else []
|
||||
csv_beads: list[ControlBead] = []
|
||||
if csv_path is not None:
|
||||
from simurg.formats.csv_manifest import read_csv_manifest
|
||||
|
||||
csv_beads = beads_from_definitions(read_csv_manifest(csv_path, load_sequences=False).controls)
|
||||
if bpm_beads and csv_beads and (diffs := compare_definitions(bpm_beads, csv_beads)):
|
||||
rep.warn("W_CTRL_DEFS_DIFFER", "BPM and CSV control definitions differ (using the BPM): " + "; ".join(diffs))
|
||||
return bpm_beads or csv_beads
|
||||
|
||||
|
||||
@dataclass
|
||||
class ControlSample:
|
||||
sample_id: str
|
||||
meta: dict[str, str]
|
||||
green: npt.NDArray[np.float64]
|
||||
red: npt.NDArray[np.float64]
|
||||
qc: ControlQC
|
||||
|
||||
|
||||
def _param(pars: str, key: str) -> str:
|
||||
from simurg.genotype.sample import _param as p
|
||||
|
||||
return p(pars, key)
|
||||
|
||||
|
||||
def idat_metadata(grn: Idat, red: Idat, sheet_values: Mapping[str, str] | None = None) -> dict[str, str]:
|
||||
"""QC-011 metadata: IDAT values, plate and well overridden by the sheet."""
|
||||
scan = red.scan_entry() or grn.scan_entry()
|
||||
sv = {fold(k): v for k, v in (sheet_values or {}).items()}
|
||||
return {
|
||||
"imaging_date": red.scan_time_iso or grn.scan_time_iso,
|
||||
"sample_plate": sv.get(fold("Sample_Plate")) or red.plate or grn.plate,
|
||||
"sample_well": sv.get(fold("Sample_Well")) or red.well or grn.well,
|
||||
"scanner_id": _param(scan.block_pars, "sherlockID") if scan else "",
|
||||
"scanner_version": scan.code_version if scan else "",
|
||||
}
|
||||
|
||||
|
||||
def control_sample(
|
||||
sample_id: str,
|
||||
grn: Idat,
|
||||
red: Idat,
|
||||
beads: Sequence[ControlBead],
|
||||
bg_offset: float,
|
||||
sheet_values: Mapping[str, str] | None = None,
|
||||
) -> ControlSample:
|
||||
g, r = control_intensities(grn, red, beads)
|
||||
qc = compute_control_metrics(beads, g, r, "genotyping", bg_offset)
|
||||
return ControlSample(sample_id, idat_metadata(grn, red, sheet_values), g, r, qc)
|
||||
|
||||
|
||||
def _compat_date(iso: str) -> str:
|
||||
try:
|
||||
t = _dt.datetime.fromisoformat(iso)
|
||||
except ValueError:
|
||||
return iso
|
||||
from simurg.genotype.sample import invariant_datetime
|
||||
|
||||
return invariant_datetime(t)
|
||||
|
||||
|
||||
def write_control_tables(
|
||||
out: Path,
|
||||
samples: Sequence[ControlSample],
|
||||
beads: Sequence[ControlBead],
|
||||
bg_offset: float,
|
||||
rep: Reporter,
|
||||
*,
|
||||
array_type: str = "genotyping",
|
||||
assay_control: Mapping[str, float] | None = None,
|
||||
compat: bool = False,
|
||||
) -> list[Path]:
|
||||
"""QC-009 and QC-011 (or their compatibility CSVs). ``assay_control`` maps sample_id → assay-control norm R."""
|
||||
metrics = metrics_for(array_type)
|
||||
missing: dict[str, list[str]] = {}
|
||||
for s in samples:
|
||||
missing.update(s.qc.missing)
|
||||
for col in s.qc.zero_divisor:
|
||||
rep.warn("W_CTRL_ZERO_DIVISOR", f"{s.sample_id}: {col} has a zero denominator", sample_id=s.sample_id)
|
||||
if missing:
|
||||
ids = sorted({c for v in missing.values() for c in v})
|
||||
rep.warn(
|
||||
"W_CTRL_MISSING",
|
||||
f"manifest lacks control(s) {', '.join(ids)}; not evaluable: {', '.join(sorted(missing))}",
|
||||
)
|
||||
raw_rows = [
|
||||
{
|
||||
"sample_id": s.sample_id,
|
||||
"category": b.category,
|
||||
"bead_address": b.address,
|
||||
"control_name": b.name,
|
||||
"color_label": b.color,
|
||||
"green_intensity": "NaN" if math.isnan(gv) else int(gv),
|
||||
"red_intensity": "NaN" if math.isnan(rv) else int(rv),
|
||||
}
|
||||
for s in samples
|
||||
for b, gv, rv in zip(beads, s.green.tolist(), s.red.tolist(), strict=True)
|
||||
]
|
||||
qc_rows = []
|
||||
for s in samples:
|
||||
row: dict[str, Any] = {
|
||||
"sample_id": s.sample_id,
|
||||
"array_type": array_type,
|
||||
**s.meta,
|
||||
"bg_offset": float(bg_offset),
|
||||
**{m.column: s.qc.values.get(m.column, math.nan) for m in metrics},
|
||||
"na_metrics": ";".join(s.qc.na_metrics),
|
||||
}
|
||||
if assay_control is not None:
|
||||
row[ASSAY_CONTROL_COLUMN] = assay_control.get(s.sample_id, math.nan)
|
||||
qc_rows.append(row)
|
||||
meta_cols = ["imaging_date", "sample_plate", "sample_well", "scanner_id", "scanner_version"]
|
||||
qc_cols = ["sample_id", "array_type", *meta_cols, "bg_offset", *(m.column for m in metrics)]
|
||||
if assay_control is not None:
|
||||
qc_cols.append(ASSAY_CONTROL_COLUMN)
|
||||
qc_cols.append("na_metrics")
|
||||
if not compat:
|
||||
raw_cols = list(raw_rows[0]) if raw_rows else ["sample_id"]
|
||||
write_tsv(out / CONTROLS_RAW, raw_cols, raw_rows)
|
||||
write_tsv(out / CONTROLS_QC, qc_cols, qc_rows)
|
||||
return [out / CONTROLS_RAW, out / CONTROLS_QC]
|
||||
ren_raw = {
|
||||
"sample_id": "SampleId",
|
||||
"category": "Category",
|
||||
"bead_address": "Beadtype",
|
||||
"control_name": "Control",
|
||||
"color_label": "Color",
|
||||
"green_intensity": "GreenIntensity",
|
||||
"red_intensity": "RedIntensity",
|
||||
}
|
||||
ren_qc = {
|
||||
"sample_id": "SampleId",
|
||||
"imaging_date": "ImagingDate",
|
||||
"sample_plate": "SamplePlate",
|
||||
"sample_well": "SampleWell",
|
||||
"scanner_id": "ScannerId",
|
||||
"scanner_version": "ScannerVersion",
|
||||
**COMPAT_METRIC_NAMES,
|
||||
}
|
||||
raw_c = [{ren_raw[k]: v for k, v in r.items()} for r in raw_rows]
|
||||
qc_keep = [c for c in qc_cols if c in ren_qc]
|
||||
qc_c = [
|
||||
{ren_qc[c]: (_compat_date(r[c]) if c == "imaging_date" else r.get(c, "")) for c in qc_keep} for r in qc_rows
|
||||
]
|
||||
write_tsv(out / COMPAT_RAW, list(ren_raw.values()), raw_c, sep=",")
|
||||
write_tsv(out / COMPAT_QC, [ren_qc[c] for c in qc_keep], qc_c, sep=",")
|
||||
return [out / COMPAT_RAW, out / COMPAT_QC]
|
||||
|
||||
|
||||
@dataclass
|
||||
class QcControlsOptions:
|
||||
out: Path
|
||||
bpm: Path | None = None
|
||||
csv: Path | None = None
|
||||
idats: Path | None = None
|
||||
sheet: Path | None = None
|
||||
array_type: str = "genotyping"
|
||||
bg_offset: float = 3000.0
|
||||
compat: bool = False
|
||||
force: bool = False
|
||||
on_duplicate_idat: str = "error"
|
||||
argv: list[str] = field(default_factory=list)
|
||||
|
||||
|
||||
def run_qc_controls(opts: QcControlsOptions, rep: Reporter) -> int:
|
||||
"""``simurg qc controls`` — control tables from IDATs without genotyping (QC-012, QC-014)."""
|
||||
from simurg.formats.bpm import read_bpm
|
||||
from simurg.formats.idat import check_pair, read_idat
|
||||
from simurg.pipeline.inputs import check_paths, resolve_sample_set, sheet_values
|
||||
|
||||
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.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:
|
||||
raise SimurgError("E_QC_THRESHOLD_RANGE", f"--bg-offset must be >= 0, got {opts.bg_offset}")
|
||||
check_paths([(opts.bpm, "BPM"), (opts.csv, "CSV manifest"), (opts.idats, "IDAT folder"), (opts.sheet, "sheet")])
|
||||
sset = resolve_sample_set(opts.idats, opts.sheet, opts.on_duplicate_idat, rep, rm)
|
||||
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 not beads:
|
||||
raise SimurgError("E_BPM_FORMAT", "the manifest has no control definitions")
|
||||
samples: list[ControlSample] = []
|
||||
for s in sset.samples:
|
||||
if s.pair is None or not s.pair.complete:
|
||||
rep.warn("W_IDAT_UNPAIRED", f"{s.sample_id}: no complete IDAT pair; skipped", sample_id=s.sample_id)
|
||||
continue
|
||||
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)))
|
||||
# QC-014: sample IDs of an existing summary must match
|
||||
summary = opts.out / "sample_summary.tsv"
|
||||
if summary.exists():
|
||||
ids = {ln.split("\t", 1)[0] for ln in summary.read_text(encoding="utf-8").splitlines()[1:] if ln}
|
||||
mine = {s.sample_id for s in samples}
|
||||
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)
|
||||
for p in written:
|
||||
rm.add_output(p, "simurg.controls/1")
|
||||
rm.summary = {"samples": len(samples)}
|
||||
rm.status = "completed"
|
||||
rm.write(opts.out / "run_manifest.qc_controls.json")
|
||||
rep.info(f"control QC for {len(samples)} sample(s) written to {opts.out}")
|
||||
return ExitCode.OK
|
||||
|
||||
|
||||
# ================================================================================================ cohort store
|
||||
@dataclass
|
||||
class SampleInfo:
|
||||
sample_id: str
|
||||
path: str
|
||||
sample_name: str = ""
|
||||
call_rate: float = math.nan
|
||||
logr_dev: float = math.nan
|
||||
sex_est: str = "U"
|
||||
has_baf: bool = False
|
||||
error: str = ""
|
||||
|
||||
|
||||
class CohortStore:
|
||||
"""Disk-backed (samples × loci) matrices: genotypes (uint8), BAF and GenCall scores (float32, optional)."""
|
||||
|
||||
def __init__(self, root: Path, n_samples: int, n_loci: int, *, baf: bool, scores: bool) -> None:
|
||||
self.root, self.shape = root, (n_samples, n_loci)
|
||||
root.mkdir(parents=True, exist_ok=True)
|
||||
self.gt = np.lib.format.open_memmap(root / "gt.npy", "w+", np.uint8, self.shape)
|
||||
self.baf = np.lib.format.open_memmap(root / "baf.npy", "w+", np.float32, self.shape) if baf else None
|
||||
self.score = np.lib.format.open_memmap(root / "score.npy", "w+", np.float32, self.shape) if scores else None
|
||||
|
||||
def close(self) -> None:
|
||||
for name in ("gt", "baf", "score"):
|
||||
a = getattr(self, name)
|
||||
if a is not None:
|
||||
a.flush()
|
||||
self.gt = self.baf = self.score = None # type: ignore[assignment]
|
||||
shutil.rmtree(self.root, ignore_errors=True)
|
||||
|
||||
|
||||
def _load_one(task: tuple[int, str, str, str, bool, bool]) -> SampleInfo:
|
||||
from simurg.formats.gtc import read_gtc
|
||||
|
||||
i, sid, path, root, want_baf, want_scores = task
|
||||
info = SampleInfo(sid, path)
|
||||
try:
|
||||
g = read_gtc(path)
|
||||
info.sample_name, info.call_rate, info.logr_dev, info.sex_est = (
|
||||
g.sample_name,
|
||||
float(g.call_rate),
|
||||
float(g.logr_dev) if g.version >= 4 else math.nan,
|
||||
g.gender_char,
|
||||
)
|
||||
gt = np.load(Path(root) / "gt.npy", mmap_mode="r+")
|
||||
if len(g.genotypes) != gt.shape[1]:
|
||||
raise SimurgError("E_GTC_UNSUITABLE", f"{path}: {len(g.genotypes)} loci, manifest has {gt.shape[1]}")
|
||||
gt[i] = g.genotypes
|
||||
gt.flush()
|
||||
if want_baf:
|
||||
baf = np.load(Path(root) / "baf.npy", mmap_mode="r+")
|
||||
info.has_baf = len(g.b_allele_freqs) == gt.shape[1]
|
||||
baf[i] = g.b_allele_freqs if info.has_baf else np.nan
|
||||
baf.flush()
|
||||
if want_scores:
|
||||
sc = np.load(Path(root) / "score.npy", mmap_mode="r+")
|
||||
sc[i] = g.genotype_scores
|
||||
sc.flush()
|
||||
except SimurgError as e:
|
||||
info.error = f"{e.code}: {e.message}"
|
||||
except Exception as e:
|
||||
info.error = f"E_INTERNAL: {type(e).__name__}: {e}"
|
||||
return info
|
||||
|
||||
|
||||
def load_cohort(
|
||||
samples: Sequence[tuple[str, str]], n_loci: int, tmp: Path, *, baf: bool, scores: bool, workers: int
|
||||
) -> tuple[CohortStore, list[SampleInfo]]:
|
||||
"""Read GTCs (sample_id, path) into a :class:`CohortStore` with a process pool."""
|
||||
store = CohortStore(tmp, len(samples), n_loci, baf=baf, scores=scores)
|
||||
store.gt.flush()
|
||||
tasks = [(i, sid, p, str(tmp), baf, scores) for i, (sid, p) in enumerate(samples)]
|
||||
if workers <= 1 or len(tasks) <= 1:
|
||||
infos = [_load_one(t) for t in tasks]
|
||||
else:
|
||||
with ProcessPoolExecutor(workers, mp_context=mp.get_context("spawn")) as ex:
|
||||
infos = list(ex.map(_load_one, tasks, chunksize=4))
|
||||
return store, infos
|
||||
|
||||
|
||||
class SiteView:
|
||||
"""(sites × samples) row access onto a (samples × loci) matrix, optionally restricted to some samples."""
|
||||
|
||||
def __init__(self, m: npt.NDArray[Any], samples: npt.NDArray[np.intp] | None = None) -> None:
|
||||
self.m, self.samples = m, samples
|
||||
self.shape = (m.shape[1], m.shape[0] if samples is None else len(samples))
|
||||
|
||||
def __getitem__(self, rows: Any) -> npt.NDArray[Any]:
|
||||
block = np.asarray(self.m[:, rows])
|
||||
return (block if self.samples is None else block[self.samples]).T
|
||||
|
||||
|
||||
def _blocks(n: int, size: int) -> Iterator[slice]:
|
||||
for s in range(0, n, size):
|
||||
yield slice(s, min(n, s + size))
|
||||
|
||||
|
||||
def genotype_counts(
|
||||
gt: npt.NDArray[np.uint8], samples: npt.NDArray[np.intp] | None = None, block: int = 65536
|
||||
) -> npt.NDArray[np.int64]:
|
||||
"""(loci × 4) counts of no-call/AA/AB/BB over (a subset of) samples, streamed in locus blocks."""
|
||||
out = np.zeros((gt.shape[1], 4), np.int64)
|
||||
for sl in _blocks(gt.shape[1], block):
|
||||
g = np.asarray(gt[:, sl])
|
||||
if samples is not None:
|
||||
g = g[samples]
|
||||
for k in range(4):
|
||||
out[sl, k] = (g == k).sum(axis=0)
|
||||
return out
|
||||
|
||||
|
||||
def het_rates(gt: npt.NDArray[np.uint8], loci: npt.NDArray[np.intp], block: int = 65536) -> npt.NDArray[np.float64]:
|
||||
het = np.zeros(gt.shape[0], np.int64)
|
||||
called = np.zeros(gt.shape[0], np.int64)
|
||||
for s in range(0, len(loci), block):
|
||||
g = np.asarray(gt[:, loci[s : s + block]])
|
||||
het += (g == 2).sum(axis=1)
|
||||
called += (g > 0).sum(axis=1)
|
||||
with np.errstate(invalid="ignore", divide="ignore"):
|
||||
return np.where(called > 0, het / called, np.nan)
|
||||
|
||||
|
||||
# ================================================================================================ pedigree
|
||||
@dataclass
|
||||
class Pedigree:
|
||||
parents: dict[str, list[str]] # child sample_id -> genotyped parents (0–2)
|
||||
replicates: dict[str, str] # sample_id -> declared replicate sample_id
|
||||
|
||||
def expected(self, a: str, b: str) -> str:
|
||||
if b in self.parents.get(a, []) or a in self.parents.get(b, []):
|
||||
return "parent-offspring"
|
||||
if self.replicates.get(a) == b or self.replicates.get(b) == a:
|
||||
return "replicate"
|
||||
pa, pb = set(self.parents.get(a, [])), set(self.parents.get(b, []))
|
||||
if pa & pb:
|
||||
return "sibling"
|
||||
return ""
|
||||
|
||||
|
||||
def pedigree_from_sheet(sheet: SampleSheet | None, sample_ids: Sequence[str], rep: Reporter) -> Pedigree:
|
||||
"""IN-018: Parent1/Parent2/Replicate values name Sample_ID, Sample_Name or sample_id (in that order)."""
|
||||
ped = Pedigree({}, {})
|
||||
if sheet is None:
|
||||
return ped
|
||||
ids = set(sample_ids)
|
||||
by_col: dict[str, dict[str, str]] = {"sample_id_col": {}, "name": {}}
|
||||
rows = [r for r in sheet.rows if r.sample_id in ids]
|
||||
cols = {fold(c): c for c in sheet.columns}
|
||||
for r in rows:
|
||||
if (c := cols.get(fold("Sample_ID"))) and r.values.get(c):
|
||||
by_col["sample_id_col"].setdefault(r.values[c], r.sample_id)
|
||||
if r.sample_name:
|
||||
by_col["name"].setdefault(r.sample_name, r.sample_id)
|
||||
|
||||
def resolve(v: str) -> str | None:
|
||||
v = v.strip()
|
||||
if v in ("", "0", "NA", "na", "-"):
|
||||
return None
|
||||
for table in (by_col["sample_id_col"], by_col["name"]):
|
||||
if v in table:
|
||||
return table[v]
|
||||
if v in ids:
|
||||
return v
|
||||
rep.warn("W_QC_PEDIGREE", f"{v!r} does not name a sample in this set")
|
||||
return None
|
||||
|
||||
for r in rows:
|
||||
for key in ("Parent1", "Parent2"):
|
||||
if (c := cols.get(fold(key))) and (p := resolve(r.values.get(c, ""))) and p != r.sample_id:
|
||||
ped.parents.setdefault(r.sample_id, []).append(p)
|
||||
if (c := cols.get(fold("Replicate"))) and (p := resolve(r.values.get(c, ""))) and p != r.sample_id:
|
||||
ped.replicates[r.sample_id] = p
|
||||
return ped
|
||||
|
||||
|
||||
# ================================================================================================ cohort QC
|
||||
@dataclass
|
||||
class CohortQC:
|
||||
het_rate: npt.NDArray[np.float64]
|
||||
het_z: list[float]
|
||||
contamination: ContaminationResult | None = None
|
||||
pairs: list[dict[str, Any]] = field(default_factory=list)
|
||||
duplicates: dict[str, list[str]] = field(default_factory=dict)
|
||||
unexpected: set[str] = field(default_factory=set)
|
||||
mendel: list[dict[str, Any]] = field(default_factory=list)
|
||||
|
||||
|
||||
def cohort_qc(
|
||||
store: CohortStore,
|
||||
infos: Sequence[SampleInfo],
|
||||
manifest: Manifest,
|
||||
rep: Reporter,
|
||||
*,
|
||||
pedigree: Pedigree | None = None,
|
||||
contamination: bool = True,
|
||||
relatedness: bool = True,
|
||||
all_pairs: bool = False,
|
||||
median: str = "reference",
|
||||
workers: int = 1,
|
||||
) -> CohortQC:
|
||||
ids = [s.sample_id for s in infos]
|
||||
ok = np.array([not s.error for s in infos])
|
||||
auto = np.flatnonzero(autosome_mask(manifest.chroms) & ~manifest.intensity_only)
|
||||
het = het_rates(store.gt, auto)
|
||||
het[~ok] = np.nan
|
||||
res = CohortQC(het_rate=het, het_z=het_zscores(het.tolist()))
|
||||
if contamination and store.baf is not None:
|
||||
usable = ok & np.array([s.has_baf for s in infos])
|
||||
for s in infos:
|
||||
if not s.error and not s.has_baf:
|
||||
rep.warn("W_QC_NO_BAF", f"{s.sample_id}: GTC has no BAF; contamination not estimated")
|
||||
if usable.sum() < MIN_CONTAMINATION_SAMPLES:
|
||||
rep.warn(
|
||||
"W_QC_TOO_FEW",
|
||||
f"contamination needs >= {MIN_CONTAMINATION_SAMPLES} samples with BAF (have {int(usable.sum())})",
|
||||
)
|
||||
else:
|
||||
rows_idx = np.flatnonzero(usable)
|
||||
sel = None if usable.all() else rows_idx
|
||||
c = genotype_counts(store.gt, sel)[:, 1:].astype(np.float64)
|
||||
an = 2 * c.sum(axis=1)
|
||||
site_af = np.full(manifest.num_loci, np.nan)
|
||||
with np.errstate(invalid="ignore", divide="ignore"):
|
||||
af = np.where(an > 0, (c[:, 1] + 2 * c[:, 2]) / an, np.nan)
|
||||
site_af[auto] = af[auto] # contamination sites: autosomes only
|
||||
with _pool(workers) as ex:
|
||||
cr = baf_regress(
|
||||
SiteView(store.gt, sel), # type: ignore[arg-type]
|
||||
SiteView(store.baf, sel), # type: ignore[arg-type]
|
||||
site_af,
|
||||
median=median,
|
||||
map_fn=ex.map if ex is not None else map,
|
||||
)
|
||||
est = np.full(len(ids), np.nan)
|
||||
nh = np.zeros(len(ids), np.int64)
|
||||
est[rows_idx], nh[rows_idx] = cr.estimate, cr.n_hom
|
||||
res.contamination = ContaminationResult(est, nh, cr.n_sites)
|
||||
if relatedness and ok.sum() >= 2:
|
||||
usable_loci = np.zeros(manifest.num_loci, np.bool_)
|
||||
usable_loci[auto] = True
|
||||
g_ok = np.flatnonzero(ok)
|
||||
sites = select_sites(genotype_counts(store.gt, None if ok.all() else g_ok), usable_loci)
|
||||
G = SiteView(store.gt, None if ok.all() else g_ok)[sites]
|
||||
ps = pair_stats(G)
|
||||
ped = pedigree or Pedigree({}, {})
|
||||
sub = [ids[i] for i in g_ok]
|
||||
for a in range(len(sub)):
|
||||
for b in range(a + 1, len(sub)):
|
||||
phi, conc = float(ps.kinship[a, b]), float(ps.concordance[a, b])
|
||||
exp = ped.expected(sub[a], sub[b])
|
||||
dup = conc >= DUPLICATE_CONCORDANCE
|
||||
if dup:
|
||||
res.duplicates.setdefault(sub[a], []).append(sub[b])
|
||||
res.duplicates.setdefault(sub[b], []).append(sub[a])
|
||||
rel = math.isfinite(phi) and phi > DEGREE_CUTS[2][0] # second degree or closer
|
||||
if rel and not exp:
|
||||
res.unexpected.update((sub[a], sub[b]))
|
||||
if all_pairs or dup or exp or (math.isfinite(phi) and phi > DEGREE_CUTS[3][0]):
|
||||
res.pairs.append(
|
||||
{
|
||||
"sample_a": sub[a],
|
||||
"sample_b": sub[b],
|
||||
"n_sites": int(ps.n_sites[a, b]),
|
||||
"kinship": phi,
|
||||
"concordance": conc,
|
||||
"degree": degree(phi) if math.isfinite(phi) else "",
|
||||
"expected": exp,
|
||||
}
|
||||
)
|
||||
for child, parents in (ped.parents or {}).items():
|
||||
if child not in ids:
|
||||
continue
|
||||
ci = ids.index(child)
|
||||
prow = [np.asarray(store.gt[ids.index(p)])[auto] for p in parents[:2]]
|
||||
e, n = mendel_errors(
|
||||
np.asarray(store.gt[ci])[auto], prow[0] if prow else None, prow[1] if len(prow) > 1 else None
|
||||
)
|
||||
res.mendel.append(
|
||||
{
|
||||
"child": child,
|
||||
"parent1": parents[0] if parents else "",
|
||||
"parent2": parents[1] if len(parents) > 1 else "",
|
||||
"n_errors": e,
|
||||
"n_informative": n,
|
||||
"error_rate": e / n if n else math.nan,
|
||||
}
|
||||
)
|
||||
return res
|
||||
|
||||
|
||||
class _NullPool:
|
||||
def __enter__(self) -> None:
|
||||
return None
|
||||
|
||||
def __exit__(self, *a: object) -> None:
|
||||
return None
|
||||
|
||||
|
||||
def _pool(workers: int) -> Any:
|
||||
if workers <= 1:
|
||||
return _NullPool()
|
||||
return ProcessPoolExecutor(workers, mp_context=mp.get_context("spawn"))
|
||||
|
||||
|
||||
def default_tmp(out: Path) -> Path:
|
||||
base = Path(os.environ.get("SIMURG_TMP", out))
|
||||
base.mkdir(parents=True, exist_ok=True)
|
||||
return Path(tempfile.mkdtemp(prefix=".simurg_cohort_", dir=base))
|
||||
|
||||
|
||||
def contamination_rows(ids: Sequence[str], cr: ContaminationResult | None) -> dict[str, dict[str, Any]]:
|
||||
if cr is None:
|
||||
return {}
|
||||
return {sid: {"contamination": float(cr.estimate[i]), "n_hom_sites": int(cr.n_hom[i])} for i, sid in enumerate(ids)}
|
||||
|
||||
|
||||
def status_from_thresholds(th: Thresholds) -> list[dict[str, Any]]:
|
||||
return th.as_records()
|
||||
|
||||
|
||||
# ================================================================================================ locus summary
|
||||
def write_locus_summary(
|
||||
path: Path, store: CohortStore, manifest: Manifest, egt: Any, samples: npt.NDArray[np.intp] | None = None
|
||||
) -> None:
|
||||
"""OUT-072, streamed in locus blocks. ``egt`` is a parsed cluster file (records joined by locus name)."""
|
||||
from simurg.core.atomic import atomic_open
|
||||
from simurg.qc.locus_summary import COLUMNS, summarize_block
|
||||
|
||||
assert store.score is not None
|
||||
idx = egt.name_index()
|
||||
rec = np.array([idx.get(nm, -1) for nm in manifest.names], dtype=np.int64)
|
||||
has = rec >= 0
|
||||
safe = np.where(has, rec, 0)
|
||||
|
||||
def take(a: npt.NDArray[Any]) -> npt.NDArray[np.float64]:
|
||||
v = np.asarray(a)[safe].astype(np.float64)
|
||||
v[~has] = np.nan
|
||||
return v
|
||||
|
||||
cl = {
|
||||
"locus_id": np.where(has, np.asarray(egt.addresses)[safe], 0),
|
||||
"gentrain_score": take(egt.gentrain_score),
|
||||
"cluster_sep": take(egt.cluster_sep),
|
||||
}
|
||||
for j, gname in enumerate(("aa", "ab", "bb")):
|
||||
cl[f"{gname}_t_mean"] = take(egt.theta_mean[:, j])
|
||||
cl[f"{gname}_t_sd"] = take(egt.theta_dev[:, j])
|
||||
cl[f"{gname}_r_mean"] = take(egt.r_mean[:, j])
|
||||
cl[f"{gname}_r_sd"] = take(egt.r_dev[:, j])
|
||||
strands = [s if s in ("+", "-") else "U" for s in (manifest.ref_strands or [""] * manifest.num_loci)]
|
||||
|
||||
def f(v: Any) -> str:
|
||||
if isinstance(v, float | np.floating):
|
||||
return "NaN" if math.isnan(v) else f"{float(v):.8g}"
|
||||
return str(v)
|
||||
|
||||
with atomic_open(path, "w", encoding="utf-8") as raw:
|
||||
fh = cast(IO[str], raw)
|
||||
fh.write("\t".join(COLUMNS) + "\n")
|
||||
for sl in _blocks(manifest.num_loci, 65536):
|
||||
g = np.asarray(store.gt[:, sl])
|
||||
s = np.asarray(store.score[:, sl])
|
||||
if samples is not None:
|
||||
g, s = g[samples], s[samples]
|
||||
stats = summarize_block(g, s)
|
||||
for k in range(sl.stop - sl.start):
|
||||
i = sl.start + k
|
||||
row = {"locus_name": manifest.names[i], "plus_minus_strand": strands[i]}
|
||||
vals = [row[c] if c in row else (cl[c][i] if c in cl else stats[c][k]) for c in COLUMNS]
|
||||
fh.write("\t".join(f(v) for v in vals) + "\n")
|
||||
|
||||
|
||||
# ================================================================================================ GTC-based commands
|
||||
SAMPLE_QC_COLUMNS = [
|
||||
"sample_id",
|
||||
"sample_name",
|
||||
"status",
|
||||
"call_rate",
|
||||
"autosomal_call_rate",
|
||||
"logr_dev",
|
||||
"sex_est",
|
||||
"sex_reported",
|
||||
"sex_check",
|
||||
"het_rate",
|
||||
"het_z",
|
||||
"contamination",
|
||||
"n_hom_sites",
|
||||
"functional_status",
|
||||
"flags",
|
||||
]
|
||||
RELATEDNESS_COLUMNS = ["sample_a", "sample_b", "n_sites", "kinship", "concordance", "degree", "expected"]
|
||||
MENDEL_COLUMNS = ["child", "parent1", "parent2", "n_errors", "n_informative", "error_rate"]
|
||||
|
||||
|
||||
@dataclass
|
||||
class QcCohortOptions:
|
||||
"""Options shared by ``qc samples | contamination | relatedness | loci`` (CLI-025)."""
|
||||
|
||||
command: str
|
||||
out: Path
|
||||
bpm: Path
|
||||
gtcs: Path | None = None
|
||||
sheet: Path | None = None
|
||||
egt: Path | None = None
|
||||
threads: int = 0
|
||||
force: bool = False
|
||||
all_pairs: bool = False
|
||||
pedigree_from_sheet: bool = True
|
||||
median: str = "reference"
|
||||
thresholds: Thresholds | None = None
|
||||
argv: list[str] = field(default_factory=list)
|
||||
|
||||
|
||||
def _gtc_samples(opts: QcCohortOptions, rep: Reporter) -> tuple[list[tuple[str, str]], SampleSheet | None]:
|
||||
"""(sample_id, path) in sheet order (sheet given) or discovery order (IN-030, IN-016)."""
|
||||
from simurg.formats.discovery import discover_gtcs
|
||||
from simurg.formats.samplesheet import parse_samplesheet
|
||||
|
||||
sheet = parse_samplesheet(opts.sheet) if opts.sheet else None
|
||||
root = opts.gtcs or (sheet.root_folder if sheet else None)
|
||||
if root is None:
|
||||
raise SimurgError("E_NO_INPUT", "give --gtcs or a sample sheet with a RootFolder")
|
||||
found = discover_gtcs(root)
|
||||
if sheet is None:
|
||||
return [(k, str(v)) for k, v in found.items()], None
|
||||
out = []
|
||||
for r in sheet.rows:
|
||||
if r.sample_id in found:
|
||||
out.append((r.sample_id, str(found[r.sample_id])))
|
||||
else:
|
||||
rep.warn("W_QC_UNMATCHED_SAMPLES", f"{r.sample_id}: no GTC under {root}", sample_id=r.sample_id)
|
||||
if not out:
|
||||
raise SimurgError("E_NO_INPUT", f"none of the sheet's samples has a GTC under {root}")
|
||||
return out, sheet
|
||||
|
||||
|
||||
def run_qc_cohort(opts: QcCohortOptions, rep: Reporter) -> int:
|
||||
"""``simurg qc samples | contamination | relatedness | loci`` (CLI-025)."""
|
||||
from simurg.formats.bpm import read_bpm
|
||||
from simurg.pipeline.inputs import check_paths
|
||||
from simurg.qc.status import functional_status, sample_flags
|
||||
from simurg.qc.thresholds import resolve_thresholds
|
||||
|
||||
rm = RunManifest(command=f"qc {opts.command}", argv=opts.argv)
|
||||
rep.run_id = rm.run_id
|
||||
check_paths([(opts.bpm, "BPM"), (opts.gtcs, "GTC folder"), (opts.sheet, "sample sheet"), (opts.egt, "EGT")])
|
||||
if opts.command == "loci" and opts.egt is None:
|
||||
raise SimurgError("E_NO_INPUT", "qc loci needs --egt (cluster statistics, OUT-072)")
|
||||
names = {
|
||||
"samples": ["sample_qc.tsv", "contamination.tsv", "relatedness.tsv"],
|
||||
"contamination": ["contamination.tsv"],
|
||||
"relatedness": ["relatedness.tsv"],
|
||||
"loci": ["locus_summary.tsv"],
|
||||
}[opts.command]
|
||||
outputs = [opts.out / n for n in names]
|
||||
check_outputs_absent(outputs, force=opts.force)
|
||||
samples, sheet = _gtc_samples(opts, rep)
|
||||
rep.info(f"reading manifest {opts.bpm}")
|
||||
manifest = read_bpm(opts.bpm)
|
||||
workers = max(1, opts.threads or (os.cpu_count() or 1))
|
||||
opts.out.mkdir(parents=True, exist_ok=True)
|
||||
tmp = default_tmp(opts.out)
|
||||
want_baf = opts.command in ("samples", "contamination")
|
||||
rep.info(f"loading {len(samples)} GTC(s) with {workers} worker(s)")
|
||||
store, infos = load_cohort(
|
||||
samples, manifest.num_loci, tmp, baf=want_baf, scores=opts.command == "loci", workers=workers
|
||||
)
|
||||
try:
|
||||
for s in infos:
|
||||
if s.error:
|
||||
rep.error(s.error.split(":", 1)[0], f"{s.sample_id}: {s.error}", sample_id=s.sample_id)
|
||||
ids = [s.sample_id for s in infos]
|
||||
if opts.command == "loci":
|
||||
egt_path = opts.egt
|
||||
assert egt_path is not None
|
||||
from simurg.formats.egt import read_egt
|
||||
|
||||
ok = np.flatnonzero([not s.error for s in infos])
|
||||
write_locus_summary(opts.out / "locus_summary.tsv", store, manifest, read_egt(egt_path), ok)
|
||||
else:
|
||||
ped = pedigree_from_sheet(sheet, ids, rep) if opts.pedigree_from_sheet else None
|
||||
res = cohort_qc(
|
||||
store,
|
||||
infos,
|
||||
manifest,
|
||||
rep,
|
||||
pedigree=ped,
|
||||
contamination=opts.command in ("samples", "contamination"),
|
||||
relatedness=opts.command in ("samples", "relatedness"),
|
||||
all_pairs=opts.all_pairs,
|
||||
median=opts.median,
|
||||
workers=workers,
|
||||
)
|
||||
crow = contamination_rows(ids, res.contamination)
|
||||
if opts.command in ("samples", "contamination"):
|
||||
write_tsv(
|
||||
opts.out / "contamination.tsv",
|
||||
["sample_id", "contamination", "n_hom_sites"],
|
||||
[{"sample_id": sid, **crow.get(sid, {})} for sid in ids],
|
||||
)
|
||||
if opts.command in ("samples", "relatedness"):
|
||||
write_tsv(opts.out / "relatedness.tsv", RELATEDNESS_COLUMNS, res.pairs)
|
||||
if res.mendel:
|
||||
write_tsv(opts.out / "mendel.tsv", MENDEL_COLUMNS, res.mendel)
|
||||
outputs.append(opts.out / "mendel.tsv")
|
||||
if opts.command == "samples":
|
||||
th = opts.thresholds or resolve_thresholds()
|
||||
auto = autosome_mask(manifest.chroms) & ~manifest.intensity_only
|
||||
auto_idx = np.flatnonzero(auto)
|
||||
sex_rep = {r.sample_id: r.reported_sex for r in sheet.rows} if sheet else {}
|
||||
rows = []
|
||||
for i, s in enumerate(infos):
|
||||
g = np.asarray(store.gt[i])[auto_idx]
|
||||
row: dict[str, Any] = {
|
||||
"sample_id": s.sample_id,
|
||||
"sample_name": s.sample_name,
|
||||
"status": "FAILED" if s.error else "OK",
|
||||
"call_rate": s.call_rate,
|
||||
"autosomal_call_rate": float((g > 0).mean()) if len(g) and not s.error else math.nan,
|
||||
"logr_dev": s.logr_dev,
|
||||
"sex_est": s.sex_est,
|
||||
"sex_reported": sex_rep.get(s.sample_id, "U"),
|
||||
"het_rate": float(res.het_rate[i]),
|
||||
"het_z": res.het_z[i],
|
||||
**crow.get(s.sample_id, {}),
|
||||
}
|
||||
known = {"M", "F"}
|
||||
row["sex_check"] = (
|
||||
("match" if row["sex_est"] == row["sex_reported"] else "mismatch")
|
||||
if row["sex_est"] in known and row["sex_reported"] in known
|
||||
else "unknown"
|
||||
)
|
||||
row["functional_status"] = functional_status(row, th)
|
||||
row["flags"] = ";".join(
|
||||
sample_flags(
|
||||
row,
|
||||
th,
|
||||
het_z=res.het_z[i],
|
||||
duplicates=res.duplicates.get(s.sample_id, []),
|
||||
unexpected_relatedness=s.sample_id in res.unexpected,
|
||||
)
|
||||
)
|
||||
rows.append(row)
|
||||
write_tsv(opts.out / "sample_qc.tsv", SAMPLE_QC_COLUMNS, rows)
|
||||
finally:
|
||||
store.close()
|
||||
for p in outputs:
|
||||
if p.exists():
|
||||
rm.add_output(p, f"simurg.{p.stem}/1")
|
||||
rm.summary = {"samples": len(infos), "failed": sum(bool(s.error) for s in infos)}
|
||||
rm.status = "completed"
|
||||
rm.write(opts.out / f"run_manifest.qc_{opts.command}.json")
|
||||
failed = sum(bool(s.error) for s in infos)
|
||||
rep.info(f"done: qc {opts.command} on {len(infos) - failed}/{len(infos)} sample(s); outputs in {opts.out}")
|
||||
if failed == 0:
|
||||
return ExitCode.OK
|
||||
return ExitCode.PARTIAL if failed < len(infos) else ExitCode.ALL_FAILED
|
||||
|
||||
|
||||
# ================================================================================================ suggest-thresholds
|
||||
def _read_table(path: Path) -> list[dict[str, str]]:
|
||||
import csv
|
||||
|
||||
text = path.read_text(encoding="utf-8-sig")
|
||||
delim = "\t" if "\t" in text.split("\n", 1)[0] else ","
|
||||
return list(csv.DictReader(text.splitlines(), delimiter=delim))
|
||||
|
||||
|
||||
def dataset_values(folder: Path) -> dict[str, list[float]]:
|
||||
"""Per-key metric values of one QC dataset folder (QC-040 file contract)."""
|
||||
from simurg.qc.thresholds import KEY_TO_COLUMN
|
||||
|
||||
def num(v: str | None) -> float:
|
||||
try:
|
||||
return float(v) if v not in (None, "") else math.nan
|
||||
except ValueError:
|
||||
return math.nan
|
||||
|
||||
out: dict[str, list[float]] = {}
|
||||
qc = folder / CONTROLS_QC
|
||||
if not qc.exists():
|
||||
raise SimurgError(
|
||||
"E_QC_DATASET_INCOMPLETE", f"Dataset {folder.name} is missing required file(s): {CONTROLS_QC}"
|
||||
)
|
||||
rows = _read_table(qc)
|
||||
for key, col in KEY_TO_COLUMN.items():
|
||||
if rows and col in rows[0]:
|
||||
out[key] = [num(r.get(col)) for r in rows]
|
||||
summary = folder / "sample_summary.tsv"
|
||||
if summary.exists():
|
||||
srows = [r for r in _read_table(summary) if r.get("status", "OK") == "OK"]
|
||||
out["callRate"] = [num(r.get("autosomal_call_rate")) for r in srows]
|
||||
out["logRDev"] = [num(r.get("logr_dev")) for r in srows]
|
||||
return out
|
||||
|
||||
|
||||
def run_suggest_thresholds(
|
||||
data: Sequence[Path], out_yaml: Path, th: Thresholds, rep: Reporter, *, k: float = 4.0, min_samples: int = 30
|
||||
) -> int:
|
||||
from simurg.qc.suggest import render_yaml, suggest
|
||||
|
||||
if not data:
|
||||
raise SimurgError("E_QC_NO_DATASETS", "No valid dataset folders detected.")
|
||||
merged: dict[str, list[float]] = {}
|
||||
for d in data:
|
||||
if not d.exists():
|
||||
raise SimurgError("E_QC_NO_DATASETS", f"Input data folder does not exist: {d}", path=str(d))
|
||||
for key, vals in dataset_values(d).items():
|
||||
merged.setdefault(key, []).extend(vals)
|
||||
sugs = suggest(merged, th, k=k, min_samples=min_samples)
|
||||
for s in sugs:
|
||||
if s.note:
|
||||
rep.warn("W_QC_TOO_FEW", f"{s.key}: {s.note}; left unset")
|
||||
out_yaml.parent.mkdir(parents=True, exist_ok=True)
|
||||
atomic_write_text(out_yaml, render_yaml(sugs, k))
|
||||
write_tsv(
|
||||
out_yaml.with_suffix(".tsv"),
|
||||
["key", "n", "centre", "spread", "value", "pct_flagged_current", "pct_flagged_suggested", "note"],
|
||||
[s.__dict__ for s in sugs],
|
||||
)
|
||||
rep.info(f"suggested thresholds written to {out_yaml}")
|
||||
return ExitCode.OK
|
||||
@@ -0,0 +1,179 @@
|
||||
"""Contamination estimate by B-allele-frequency regression (spec QC-064, plan §9.3).
|
||||
|
||||
[LIT: Jun G. et al. 2012, AJHG 91:839 — BAF regression on homozygous calls]
|
||||
[OSS: gtc2vcf ``BAFregress.c`` (MIT) — reference numerics: per-site median adjustment of AA/BB BAF when a cluster
|
||||
has >= 5 calls, float32 accumulators updated with double-precision addends, OLS slope in double.]
|
||||
|
||||
For every site with cohort B-allele frequency p in [t, 1 - t] (t = 0.1f):
|
||||
AA call: x = p, y = BAF - median_AA(BAF)
|
||||
BB call: x = 1 - p, y = 1 - (BAF - (median_BB(BAF) - 1))
|
||||
The contamination fraction is the slope m = (n·Σxy − Σx·Σy) / (n·Σx² − (Σx)²).
|
||||
|
||||
QUIRK (K-017): for an even number of calls the reference averages the quickselect result with whatever element sits
|
||||
just below it after the partial partition (klib ``ks_ksmall``), which is not always the lower middle value.
|
||||
``median="reference"`` (default) reproduces it exactly; ``median="exact"`` uses the true median.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import math
|
||||
from collections.abc import Callable, Iterable
|
||||
from dataclasses import dataclass
|
||||
from typing import Any
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
F32 = np.float32
|
||||
CHUNK = 4096
|
||||
|
||||
|
||||
@dataclass
|
||||
class ContaminationResult:
|
||||
estimate: npt.NDArray[np.float64] # per sample; NaN when undefined
|
||||
n_hom: npt.NDArray[np.int64]
|
||||
n_sites: int
|
||||
|
||||
|
||||
def b_allele_frequency(genotypes: npt.NDArray[np.uint8]) -> npt.NDArray[np.float64]:
|
||||
"""Per-site B-allele frequency over called samples (rows = sites); NaN where nobody is called."""
|
||||
n_ab = (genotypes == 2).sum(axis=1)
|
||||
n_bb = (genotypes == 3).sum(axis=1)
|
||||
an = 2 * ((genotypes == 1).sum(axis=1) + n_ab + n_bb)
|
||||
with np.errstate(invalid="ignore", divide="ignore"):
|
||||
return np.where(an > 0, (n_ab + 2 * n_bb) / an.astype(np.float64), np.nan)
|
||||
|
||||
|
||||
def ks_ksmall(arr: list[float], kk: int) -> float:
|
||||
"""k-th smallest by in-place quickselect, partitioning exactly as klib ``ksort.h`` (MIT; htslib copy)
|
||||
so that the element left at ``arr[kk-1]`` matches the reference."""
|
||||
low, high, k = 0, len(arr) - 1, kk
|
||||
while True:
|
||||
if high <= low:
|
||||
return arr[k]
|
||||
if high == low + 1:
|
||||
if arr[high] < arr[low]:
|
||||
arr[low], arr[high] = arr[high], arr[low]
|
||||
return arr[k]
|
||||
mid = low + (high - low) // 2
|
||||
if arr[high] < arr[mid]:
|
||||
arr[mid], arr[high] = arr[high], arr[mid]
|
||||
if arr[high] < arr[low]:
|
||||
arr[low], arr[high] = arr[high], arr[low]
|
||||
if arr[low] < arr[mid]:
|
||||
arr[mid], arr[low] = arr[low], arr[mid]
|
||||
arr[mid], arr[low + 1] = arr[low + 1], arr[mid]
|
||||
ll, hh = low + 1, high
|
||||
pivot = arr[low]
|
||||
while True:
|
||||
ll += 1
|
||||
while arr[ll] < pivot:
|
||||
ll += 1
|
||||
hh -= 1
|
||||
while pivot < arr[hh]:
|
||||
hh -= 1
|
||||
if hh < ll:
|
||||
break
|
||||
arr[ll], arr[hh] = arr[hh], arr[ll]
|
||||
arr[low], arr[hh] = arr[hh], arr[low]
|
||||
if hh <= k:
|
||||
low = ll
|
||||
if hh >= k:
|
||||
high = hh - 1
|
||||
|
||||
|
||||
def _median_reference(values: npt.NDArray[np.float32]) -> np.float32:
|
||||
t = values.tolist() # float32 → exact Python floats; comparisons are unchanged
|
||||
j = len(t)
|
||||
m = F32(ks_ksmall(t, j // 2))
|
||||
if j % 2 == 0:
|
||||
m = F32((m + F32(t[j // 2 - 1])) * F32(0.5))
|
||||
return m
|
||||
|
||||
|
||||
def _median_exact(baf: npt.NDArray[np.float32], mask: npt.NDArray[np.bool_]) -> npt.NDArray[np.float32]:
|
||||
"""Row-wise true median of ``baf`` where ``mask``; float32, even counts average the middle pair as (a+b)*0.5f."""
|
||||
v = np.where(mask, baf, F32(np.inf)).astype(F32)
|
||||
v.sort(axis=1)
|
||||
k = mask.sum(axis=1)
|
||||
rows = np.arange(len(v))
|
||||
hi = v[rows, k // 2]
|
||||
lo = v[rows, np.maximum(k // 2 - 1, 0)]
|
||||
return np.where(k % 2 == 0, ((lo + hi) * F32(0.5)).astype(F32), hi).astype(F32)
|
||||
|
||||
|
||||
def cluster_adjustments(
|
||||
task: tuple[npt.NDArray[np.uint8], npt.NDArray[np.float32], int, str],
|
||||
) -> tuple[npt.NDArray[np.float32], npt.NDArray[np.float32]]:
|
||||
"""Per-site AA and BB BAF adjustments for one chunk of sites (picklable for process pools)."""
|
||||
g, b, adjust_min, median = task
|
||||
adj_aa = np.zeros(len(g), F32)
|
||||
adj_bb = np.zeros(len(g), F32)
|
||||
if adjust_min == -1:
|
||||
return adj_aa, adj_bb
|
||||
for code, out, shift in ((1, adj_aa, F32(0.0)), (3, adj_bb, F32(1.0))):
|
||||
mask = g == code
|
||||
use = np.flatnonzero(mask.sum(axis=1) >= adjust_min)
|
||||
if not len(use):
|
||||
continue
|
||||
if median == "exact":
|
||||
out[use] = _median_exact(b[use], mask[use]) - shift
|
||||
else:
|
||||
for i in use:
|
||||
out[i] = _median_reference(b[i][mask[i]]) - shift
|
||||
return adj_aa, adj_bb
|
||||
|
||||
|
||||
MapFn = Callable[[Callable[[Any], Any], Iterable[Any]], Iterable[Any]]
|
||||
|
||||
|
||||
def baf_regress(
|
||||
genotypes: npt.NDArray[np.uint8],
|
||||
baf: npt.NDArray[np.float32],
|
||||
af: npt.NDArray[np.float64] | None = None,
|
||||
*,
|
||||
threshold: float = 0.1,
|
||||
adjust_min: int = 5,
|
||||
median: str = "reference",
|
||||
map_fn: MapFn = map,
|
||||
) -> ContaminationResult:
|
||||
"""Contamination per sample. ``genotypes`` (0 no call, 1 AA, 2 AB, 3 BB) and ``baf`` are (sites × samples);
|
||||
``af`` defaults to the cohort B-allele frequency. Sites are used in row order (the reference's float32
|
||||
accumulation is order-dependent). ``map_fn`` (e.g. an executor's ``map``) parallelises the per-site medians."""
|
||||
if median not in ("reference", "exact"):
|
||||
raise ValueError(f"median must be 'reference' or 'exact', got {median!r}")
|
||||
n = genotypes.shape[1]
|
||||
if af is None:
|
||||
af = b_allele_frequency(genotypes)
|
||||
t = float(F32(threshold)) # the reference stores the threshold as float: 0.1f = 0.10000000149…
|
||||
keep = np.flatnonzero(np.isfinite(af) & (af >= t) & (af <= 1.0 - t))
|
||||
chunks = [keep[s : s + CHUNK] for s in range(0, len(keep), CHUNK)]
|
||||
tasks = ((np.asarray(genotypes[r]), np.asarray(baf[r], dtype=F32), adjust_min, median) for r in chunks)
|
||||
sx2, sxy, sx, sy = (np.zeros(n, F32) for _ in range(4))
|
||||
cnt = np.zeros(n, np.int64)
|
||||
for rows, (adj_aa, adj_bb) in zip(chunks, map_fn(cluster_adjustments, tasks), strict=True):
|
||||
g = np.asarray(genotypes[rows])
|
||||
b = np.asarray(baf[rows], dtype=F32)
|
||||
p = af[rows]
|
||||
is_aa, is_bb = g == 1, g == 3
|
||||
y_aa = (b - adj_aa[:, None]).astype(F32).astype(np.float64)
|
||||
y_bb = (b - adj_bb[:, None]).astype(F32).astype(np.float64)
|
||||
for i in range(len(rows)):
|
||||
a, m = is_aa[i], is_aa[i] | is_bb[i]
|
||||
if not m.any():
|
||||
continue
|
||||
pa, pb = p[i], 1.0 - p[i]
|
||||
x = np.where(a, pa, pb)[m]
|
||||
xy = np.where(a, pa * y_aa[i], pb * (1.0 - y_bb[i]))[m]
|
||||
y = np.where(a, y_aa[i], 1.0 - y_bb[i])[m]
|
||||
sx2[m] = (sx2[m].astype(np.float64) + x * x).astype(F32)
|
||||
sxy[m] = (sxy[m].astype(np.float64) + xy).astype(F32)
|
||||
sx[m] = (sx[m].astype(np.float64) + x).astype(F32)
|
||||
sy[m] = (sy[m].astype(np.float64) + y).astype(F32)
|
||||
cnt[m] += 1
|
||||
nf = cnt.astype(np.float64)
|
||||
fx2, fxy, fx, fy = (s.astype(np.float64) for s in (sx2, sxy, sx, sy))
|
||||
denom = nf * fx2 - fx * fx
|
||||
with np.errstate(invalid="ignore", divide="ignore"):
|
||||
est = np.where(denom != 0, (nf * fxy - fx * fy) / denom, math.nan)
|
||||
return ContaminationResult(estimate=est, n_hom=cnt, n_sites=len(keep))
|
||||
@@ -0,0 +1,84 @@
|
||||
# Canonical control-probe IDs (spec QC-002). Keys are *normalised* names: case-folded, whitespace collapsed,
|
||||
# no spaces next to or inside parentheses, "bgnd" written as "bkg". Several manifest spellings map to one ID.
|
||||
canonical:
|
||||
# staining
|
||||
biotin(high): biotin_high
|
||||
biotin(bkg): biotin_bkg
|
||||
biotin(5k): biotin_5k
|
||||
dnp(high): dnp_high
|
||||
dnp(bkg): dnp_bkg
|
||||
dnp(20k): dnp_20k
|
||||
# extension
|
||||
extension(a): ext_A
|
||||
extension(t): ext_T
|
||||
extension(c): ext_C
|
||||
extension(g): ext_G
|
||||
# hybridization
|
||||
hyb(high): hyb_high
|
||||
hyb(medium): hyb_medium
|
||||
hyb(low): hyb_low
|
||||
# target removal
|
||||
target removal: target_removal_1
|
||||
target removal(1): target_removal_1
|
||||
target removal 1: target_removal_1
|
||||
target removal(2): target_removal_2
|
||||
target removal 2: target_removal_2
|
||||
# stringency
|
||||
string(pm): string_pm
|
||||
string(mm): string_mm
|
||||
# non-specific binding (one bead type per colour label on genotyping chips)
|
||||
nsb(bkg): nsb_bkg
|
||||
nsb(bkg) red: nsb_bkg
|
||||
nsb(bkg) green: nsb_bkg
|
||||
nsb(bkg) blue: nsb_bkg
|
||||
nsb(bkg) purple: nsb_bkg
|
||||
# non-polymorphic
|
||||
np(a): np_A
|
||||
np(t): np_T
|
||||
np(c): np_C
|
||||
np(g): np_G
|
||||
# restoration (FFPE)
|
||||
restoration: restoration
|
||||
restore: restoration
|
||||
# methylation (metrics in v2; names mapped already so controls_raw.tsv is labelled consistently)
|
||||
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 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 ii-1: bs2_1
|
||||
bs conversion ii-2: bs2_2
|
||||
bs conversion ii-3: bs2_3
|
||||
bs conversion ii-4: bs2_4
|
||||
gt mismatch 1(pm): gtmm_1_pm
|
||||
gt mismatch 2(pm): gtmm_2_pm
|
||||
gt mismatch 3(pm): gtmm_3_pm
|
||||
gt mismatch 4(pm): gtmm_4_pm
|
||||
gt mismatch 5(pm): gtmm_5_pm
|
||||
gt mismatch 6(pm): gtmm_6_pm
|
||||
gt mismatch 1(mm): gtmm_1_mm
|
||||
gt mismatch 2(mm): gtmm_2_mm
|
||||
gt mismatch 3(mm): gtmm_3_mm
|
||||
gt mismatch 4(mm): gtmm_4_mm
|
||||
gt mismatch 5(mm): gtmm_5_mm
|
||||
gt mismatch 6(mm): gtmm_6_mm
|
||||
specificity 1: spec2_1
|
||||
specificity 2: spec2_2
|
||||
specificity 3: spec2_3
|
||||
|
||||
# Categories treated as methylation-only (dataset type detection, QC-043).
|
||||
methylation_categories:
|
||||
- bisulfite conversion i
|
||||
- bisulfite conversion ii
|
||||
- specificity i
|
||||
- specificity ii
|
||||
|
||||
# Report exclusion list (QC-010): raw-probe views leave these out; derived metrics still use them.
|
||||
report_exclude:
|
||||
categories: [negative, norm]
|
||||
names: ["np(g) 1", "np(g) 2", "np(g) 3", "np(g) 4", "np(g) 5", "biotin(5k)", "dnp(20k)"]
|
||||
@@ -0,0 +1,257 @@
|
||||
"""Control-probe QC metrics (spec QC-001 … QC-011; plan §9.1).
|
||||
|
||||
Controls are the non-polymorphic bead types listed in the manifest's control section. Each gives a raw green and
|
||||
red intensity (IDAT means, no normalization). Derived metrics are signal/background ratios per assay step.
|
||||
[DOC: Illumina BeadArray Controls Reporter software guide and Infinium controls training material — formulas only]
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import math
|
||||
import re
|
||||
from collections.abc import Sequence
|
||||
from dataclasses import dataclass, field
|
||||
from functools import cache
|
||||
from importlib import resources
|
||||
from typing import Any
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
import yaml
|
||||
|
||||
from simurg.formats.idat import Idat
|
||||
from simurg.signal.raw_xy import AddressIndex
|
||||
|
||||
DEFAULT_BG_OFFSET = 3000.0
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class ControlBead:
|
||||
"""One control bead type: a (manifest control line, distinct address) pair, in manifest order."""
|
||||
|
||||
category: str
|
||||
color: str
|
||||
name: str
|
||||
address: int
|
||||
canonical: str # canonical ID, "" when the name is not in control_names.yaml
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class Metric:
|
||||
column: str # controls_qc.tsv column
|
||||
key: str # threshold YAML key (QC-020)
|
||||
genotyping: bool
|
||||
methylation: bool
|
||||
|
||||
|
||||
# QC-004, in table order.
|
||||
METRICS: tuple[Metric, ...] = (
|
||||
Metric("restoration", "restoration", True, True),
|
||||
Metric("staining_green", "stainingGreen", True, True),
|
||||
Metric("staining_red", "stainingRed", True, True),
|
||||
Metric("extension_green", "extensionGreen", True, True),
|
||||
Metric("extension_red", "extensionRed", True, True),
|
||||
Metric("hyb_high_med", "hybGreenHighMed", True, True),
|
||||
Metric("hyb_med_low", "hybGreenMediumLow", True, True),
|
||||
Metric("target_removal_1", "targetRemovalI", True, True),
|
||||
Metric("target_removal_2", "targetRemovalII", False, True),
|
||||
Metric("stringency", "stringency", True, False),
|
||||
Metric("nsb_green", "nonSpecificBindingGreen", True, False),
|
||||
Metric("nsb_red", "nonSpecificBindingRed", True, False),
|
||||
Metric("np_green", "npGreen", True, True),
|
||||
Metric("np_red", "npRed", True, True),
|
||||
Metric("specificity_1_green", "specificityIGreen", False, True),
|
||||
Metric("specificity_1_red", "specificityIRed", False, True),
|
||||
Metric("specificity_2", "specificityII", False, True),
|
||||
Metric("specificity_2_bg", "specificityIIBg", False, True),
|
||||
Metric("bs_conv_1_green", "bisulfiteConversionIGreen", False, True),
|
||||
Metric("bs_conv_1_green_bg", "bisulfiteConversionIGreenBg", False, True),
|
||||
Metric("bs_conv_1_red", "bisulfiteConversionIRed", False, True),
|
||||
Metric("bs_conv_1_red_bg", "bisulfiteConversionIRedBg", False, True),
|
||||
Metric("bs_conv_2", "bisulfiteConversionII", False, True),
|
||||
Metric("bs_conv_2_bg", "bisulfiteConversionIIBg", False, True),
|
||||
)
|
||||
ASSAY_CONTROL_COLUMN = "assay_control_norm_r" # QC-005: copied from the sample summary
|
||||
COMPAT_METRIC_NAMES = {
|
||||
"restoration": "RestorationQC",
|
||||
"staining_green": "StainingGreenQC",
|
||||
"staining_red": "StainingRedQC",
|
||||
"extension_green": "ExtensionGreenQC",
|
||||
"extension_red": "ExtensionRedQC",
|
||||
"hyb_high_med": "HybridizationHighMediumQC",
|
||||
"hyb_med_low": "HybridizationMediumLowQC",
|
||||
"target_removal_1": "TargetRemovalQC",
|
||||
"target_removal_2": "TargetRemoval2QC",
|
||||
"stringency": "StringencyQC",
|
||||
"nsb_green": "NonSpecificBindingGreenQC",
|
||||
"nsb_red": "NonSpecificBindingRedQC",
|
||||
"np_green": "NonPolymorphicGreenQC",
|
||||
"np_red": "NonPolymorphicRedQC",
|
||||
ASSAY_CONTROL_COLUMN: "AssayControlQC",
|
||||
}
|
||||
|
||||
|
||||
def metrics_for(array_type: str) -> list[Metric]:
|
||||
"""QC-005 / QC-006 metric sets."""
|
||||
if array_type == "genotyping":
|
||||
return [m for m in METRICS if m.genotyping]
|
||||
return [m for m in METRICS if m.methylation]
|
||||
|
||||
|
||||
# ------------------------------------------------------------------------------------------------ names
|
||||
@cache
|
||||
def name_config() -> dict[str, Any]:
|
||||
text = resources.files("simurg.qc").joinpath("control_names.yaml").read_text(encoding="utf-8")
|
||||
cfg: dict[str, Any] = yaml.safe_load(text)
|
||||
return cfg
|
||||
|
||||
|
||||
def normalise_name(name: str) -> str:
|
||||
"""QC-002 normalisation: case-fold, collapse whitespace, no spaces before or inside parentheses, bgnd = bkg."""
|
||||
s = " ".join(name.casefold().split())
|
||||
s = re.sub(r"\s*\(\s*", "(", s)
|
||||
s = re.sub(r"\s*\)", ")", s)
|
||||
return re.sub(r"\bbgnd\b|\bbkgd\b", "bkg", s)
|
||||
|
||||
|
||||
def canonical_id(name: str) -> str:
|
||||
return str(name_config()["canonical"].get(normalise_name(name), ""))
|
||||
|
||||
|
||||
def excluded_from_report(bead: ControlBead) -> bool:
|
||||
"""QC-010."""
|
||||
ex = name_config()["report_exclude"]
|
||||
return bead.category.casefold() in ex["categories"] or normalise_name(bead.name) in ex["names"]
|
||||
|
||||
|
||||
def is_methylation_controls(beads: Sequence[ControlBead]) -> bool:
|
||||
cats = set(name_config()["methylation_categories"])
|
||||
return any(b.category.casefold() in cats for b in beads)
|
||||
|
||||
|
||||
def beads_from_definitions(defs: Sequence[dict[str, str]]) -> list[ControlBead]:
|
||||
"""Manifest control lines (``address`` may be ``a:b:…``) → bead types in manifest order (QC-001)."""
|
||||
out: list[ControlBead] = []
|
||||
for d in defs:
|
||||
seen: set[int] = set()
|
||||
for part in str(d["address"]).split(":"):
|
||||
part = part.strip()
|
||||
if not part:
|
||||
continue
|
||||
addr = int(part)
|
||||
if addr in seen:
|
||||
continue
|
||||
seen.add(addr)
|
||||
out.append(ControlBead(d["category"], d["color"], d["name"], addr, canonical_id(d["name"])))
|
||||
return out
|
||||
|
||||
|
||||
def compare_definitions(a: Sequence[ControlBead], b: Sequence[ControlBead]) -> list[str]:
|
||||
"""Differences between two control lists (BPM vs CSV, QC-001); empty when they agree."""
|
||||
ka = {(x.address, normalise_name(x.name)) for x in a}
|
||||
kb = {(x.address, normalise_name(x.name)) for x in b}
|
||||
diffs = [f"only in BPM: {addr} {nm}" for addr, nm in sorted(ka - kb)[:5]]
|
||||
diffs += [f"only in CSV: {addr} {nm}" for addr, nm in sorted(kb - ka)[:5]]
|
||||
return diffs
|
||||
|
||||
|
||||
# ------------------------------------------------------------------------------------------------ intensities
|
||||
def control_intensities(
|
||||
grn: Idat, red: Idat, beads: Sequence[ControlBead]
|
||||
) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]:
|
||||
"""Raw green/red means per bead type; NaN where the IDAT lacks the address."""
|
||||
addrs = np.array([b.address for b in beads], dtype=np.int64)
|
||||
out = []
|
||||
for idat in (grn, red):
|
||||
rows, found = AddressIndex(idat).lookup(addrs)
|
||||
v = idat.mean[rows].astype(np.float64)
|
||||
v[~found] = math.nan
|
||||
out.append(v)
|
||||
return out[0], out[1]
|
||||
|
||||
|
||||
# ------------------------------------------------------------------------------------------------ metrics
|
||||
@dataclass
|
||||
class ControlQC:
|
||||
values: dict[str, float] # column -> value (NaN when not evaluable or degenerate)
|
||||
missing: dict[str, list[str]] = field(default_factory=dict) # column -> canonical IDs absent from the manifest
|
||||
zero_divisor: list[str] = field(default_factory=list) # columns with a zero denominator (QC-008)
|
||||
|
||||
@property
|
||||
def na_metrics(self) -> list[str]:
|
||||
return list(self.missing)
|
||||
|
||||
|
||||
class _Probes:
|
||||
def __init__(self, beads: Sequence[ControlBead], green: npt.NDArray[np.float64], red: npt.NDArray[np.float64]):
|
||||
self.idx: dict[str, list[int]] = {}
|
||||
for i, b in enumerate(beads):
|
||||
if b.canonical:
|
||||
self.idx.setdefault(b.canonical, []).append(i)
|
||||
self.ch = {"G": green, "R": red}
|
||||
self.absent: set[str] = set()
|
||||
|
||||
def vals(self, ids: Sequence[str], ch: str) -> npt.NDArray[np.float64]:
|
||||
rows: list[int] = []
|
||||
for cid in ids:
|
||||
if cid not in self.idx:
|
||||
self.absent.add(cid)
|
||||
rows.extend(self.idx.get(cid, []))
|
||||
return self.ch[ch][rows]
|
||||
|
||||
def mean(self, cid: str, ch: str) -> float:
|
||||
v = self.vals([cid], ch)
|
||||
return float(v.mean()) if len(v) else math.nan
|
||||
|
||||
def min(self, ids: Sequence[str], ch: str) -> float:
|
||||
v = self.vals(ids, ch)
|
||||
return float(v.min()) if len(v) else math.nan
|
||||
|
||||
def max(self, ids: Sequence[str], ch: str) -> float:
|
||||
v = self.vals(ids, ch)
|
||||
return float(v.max()) if len(v) else math.nan
|
||||
|
||||
|
||||
def compute_control_metrics(
|
||||
beads: Sequence[ControlBead],
|
||||
green: npt.NDArray[np.float64],
|
||||
red: npt.NDArray[np.float64],
|
||||
array_type: str = "genotyping",
|
||||
bg_offset: float = DEFAULT_BG_OFFSET,
|
||||
) -> 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
|
||||
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),
|
||||
"staining_green": (lambda: p.mean("biotin_high", "G"), lambda: p.mean("biotin_bkg", "G")),
|
||||
"staining_red": (lambda: p.mean("dnp_high", "R"), lambda: p.mean("dnp_bkg", "R")),
|
||||
"extension_green": (lambda: p.min(["ext_C", "ext_G"], "G"), lambda: p.max(["ext_A", "ext_T"], "G")),
|
||||
"extension_red": (lambda: p.min(["ext_A", "ext_T"], "R"), lambda: p.max(["ext_C", "ext_G"], "R")),
|
||||
"hyb_high_med": (lambda: p.mean("hyb_high", "G"), lambda: p.mean("hyb_medium", "G")),
|
||||
"hyb_med_low": (lambda: p.mean("hyb_medium", "G"), lambda: p.mean("hyb_low", "G")),
|
||||
"target_removal_1": (lambda: p.max(["ext_A", "ext_T"], "G") + x, lambda: p.mean("target_removal_1", "G")),
|
||||
"stringency": (lambda: p.mean("string_pm", "R"), lambda: p.mean("string_mm", "R")),
|
||||
"nsb_green": (lambda: p.max(["ext_A", "ext_T"], "G") + x, lambda: p.max(["nsb_bkg"], "G")),
|
||||
"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")),
|
||||
}
|
||||
res = ControlQC(values={})
|
||||
for m in metrics_for(array_type):
|
||||
p.absent = set()
|
||||
num, den = formulas[m.column][0](), formulas[m.column][1]()
|
||||
if p.absent:
|
||||
res.missing[m.column] = sorted(p.absent)
|
||||
res.values[m.column] = math.nan
|
||||
elif den == 0:
|
||||
res.zero_divisor.append(m.column)
|
||||
res.values[m.column] = math.nan
|
||||
else:
|
||||
res.values[m.column] = num / den # NaN propagates when a bead's address is absent from the IDAT
|
||||
return res
|
||||
@@ -0,0 +1,135 @@
|
||||
"""Per-locus summary statistics across samples (spec OUT-072, CLI-025; plan §7.9).
|
||||
|
||||
[OSS: BeadArrayFiles ``examples/locus_summary.py`` (BSD-2) — formulas: call/genotype frequencies, minor frequency,
|
||||
k-regularised heterozygote excess with k = (0.013/q)^4, ChiTest_P100 = 2·Φ(±dh; 0, 0.1), GenomeStudio percentile of
|
||||
the called GenCall scores.] Vectorised over a block of loci; the caller streams blocks for large cohorts.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import math
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
COLUMNS = [
|
||||
"locus_name",
|
||||
"locus_id",
|
||||
"n_nocalls",
|
||||
"n_calls",
|
||||
"call_freq",
|
||||
"aa_freq",
|
||||
"ab_freq",
|
||||
"bb_freq",
|
||||
"minor_freq",
|
||||
"gentrain_score",
|
||||
"gc50",
|
||||
"gc10",
|
||||
"het_excess",
|
||||
"chi_test_p100",
|
||||
"cluster_sep",
|
||||
"aa_t_mean",
|
||||
"aa_t_sd",
|
||||
"ab_t_mean",
|
||||
"ab_t_sd",
|
||||
"bb_t_mean",
|
||||
"bb_t_sd",
|
||||
"aa_r_mean",
|
||||
"aa_r_sd",
|
||||
"ab_r_mean",
|
||||
"ab_r_sd",
|
||||
"bb_r_mean",
|
||||
"bb_r_sd",
|
||||
"plus_minus_strand",
|
||||
"het_excess_simple",
|
||||
]
|
||||
HET_EXCESS_EPS = 1e-9 # OQ-030: denominator guard of the simple form
|
||||
M_SQRT1_2 = 0.70710678118654752440 # C literal (1/sqrt(2) in Python is 1 ULP lower)
|
||||
|
||||
|
||||
def _norm_cdf_sd01(x: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]:
|
||||
"""Φ(x; 0, 0.1) with the branch structure of the cephes ``ndtr`` used by the reference."""
|
||||
z = (x / 0.1) * M_SQRT1_2
|
||||
out = np.empty_like(z)
|
||||
for i, v in enumerate(z.tolist()):
|
||||
a = abs(v)
|
||||
if a < M_SQRT1_2:
|
||||
y = 0.5 + 0.5 * math.erf(v)
|
||||
else:
|
||||
y = 0.5 * math.erfc(a)
|
||||
if v > 0:
|
||||
y = 1.0 - y
|
||||
out[i] = y
|
||||
return out
|
||||
|
||||
|
||||
def gs_percentile_columns(
|
||||
sorted_scores: npt.NDArray[np.float32], counts: npt.NDArray[np.int64], p: int
|
||||
) -> npt.NDArray[np.float64]:
|
||||
"""GenomeStudio percentile of the first ``counts[j]`` values of each column of ``sorted_scores`` (double math)."""
|
||||
n = counts.astype(np.int64)
|
||||
pos = n * p / 100.0
|
||||
idx = pos.astype(np.int64) # int(num*p/100) truncates like the reference
|
||||
idx = np.where(pos - idx < 0.5, idx - 1, idx)
|
||||
cols = np.arange(sorted_scores.shape[1])
|
||||
last = np.maximum(n - 1, 0)
|
||||
lo = np.clip(idx, 0, np.maximum(last - 1, 0))
|
||||
y1 = sorted_scores[lo, cols].astype(np.float64)
|
||||
y2 = sorted_scores[np.minimum(lo + 1, last), cols].astype(np.float64)
|
||||
with np.errstate(invalid="ignore", divide="ignore"):
|
||||
x1 = 100 * (lo + 0.5) / n
|
||||
x2 = 100 * (lo + 1 + 0.5) / n
|
||||
mid = y1 + (y2 - y1) / (x2 - x1) * (p - x1)
|
||||
first = sorted_scores[0, cols].astype(np.float64)
|
||||
lastv = sorted_scores[last, cols].astype(np.float64)
|
||||
out = np.where(idx < 0, first, np.where(idx >= n - 1, lastv, mid))
|
||||
return np.where(n == 0, np.nan, out)
|
||||
|
||||
|
||||
def summarize_block(genotypes: npt.NDArray[np.uint8], scores: npt.NDArray[np.float32]) -> dict[str, npt.NDArray]: # type: ignore[type-arg]
|
||||
"""Statistics for a block of loci; inputs are (samples × loci)."""
|
||||
g = np.asarray(genotypes)
|
||||
nc = (g == 0).sum(axis=0).astype(np.int64)
|
||||
aa = (g == 1).sum(axis=0).astype(np.int64)
|
||||
ab = (g == 2).sum(axis=0).astype(np.int64)
|
||||
bb = (g == 3).sum(axis=0).astype(np.int64)
|
||||
calls = aa + ab + bb
|
||||
with np.errstate(invalid="ignore", divide="ignore"):
|
||||
call_freq = np.where(calls + nc > 0, calls / (calls + nc), np.nan)
|
||||
cf = calls.astype(np.float64)
|
||||
aa_f = np.where(calls > 0, aa / cf, np.nan)
|
||||
ab_f = np.where(calls > 0, ab / cf, np.nan)
|
||||
bb_f = np.where(calls > 0, bb / cf, np.nan)
|
||||
a_freq = np.where(calls > 0, (aa * 2 + ab) / (2 * cf), np.nan)
|
||||
q = np.minimum(a_freq, 1.0 - a_freq)
|
||||
p = 1 - q
|
||||
temp = 0.013 / q
|
||||
k = temp * temp * temp * temp
|
||||
dh = ((ab / cf + k) / (2 * p * q + k)) - 1
|
||||
simple = (ab_f - 2 * p * q) / (2 * p * q + HET_EXCESS_EPS)
|
||||
mono = (aa + ab == 0) | (ab + bb == 0)
|
||||
regular = (calls > 0) & ~mono
|
||||
hw = np.zeros(len(calls))
|
||||
het = np.zeros(len(calls))
|
||||
if regular.any():
|
||||
d = dh[regular]
|
||||
cdf = _norm_cdf_sd01(d)
|
||||
hw[regular] = np.where(d < 0, 2 * cdf, 2 * (1 - cdf))
|
||||
het[regular] = d
|
||||
hw[(calls > 0) & mono] = 1.0
|
||||
called_scores = np.where(g != 0, scores.astype(np.float32), np.float32(np.inf))
|
||||
called_scores.sort(axis=0)
|
||||
return {
|
||||
"n_nocalls": nc,
|
||||
"n_calls": calls,
|
||||
"call_freq": call_freq,
|
||||
"aa_freq": aa_f,
|
||||
"ab_freq": ab_f,
|
||||
"bb_freq": bb_f,
|
||||
"minor_freq": q,
|
||||
"gc50": gs_percentile_columns(called_scores, calls, 50),
|
||||
"gc10": gs_percentile_columns(called_scores, calls, 10),
|
||||
"het_excess": het,
|
||||
"chi_test_p100": hw,
|
||||
"het_excess_simple": np.where(calls > 0, simple, np.nan),
|
||||
}
|
||||
@@ -0,0 +1,21 @@
|
||||
# Starting point for review, not a truth (spec QC-024): recalibrate on your own data with
|
||||
# `simurg qc suggest-thresholds`. For LCG/HTS genotyping chips.
|
||||
qc:
|
||||
report:
|
||||
callRate: 0.98
|
||||
logRDev: 0.20
|
||||
stainingGreen: 5
|
||||
stainingRed: 5
|
||||
extensionGreen: 5
|
||||
extensionRed: 5
|
||||
hybGreenHighMed: 1
|
||||
hybGreenMediumLow: 1
|
||||
targetRemovalI: 1
|
||||
restoration: null
|
||||
stringency: 2
|
||||
nonSpecificBindingGreen: 1
|
||||
nonSpecificBindingRed: 1
|
||||
npGreen: 5
|
||||
npRed: 5
|
||||
assayControl: 1
|
||||
bgOffset: 3000
|
||||
@@ -0,0 +1,21 @@
|
||||
# Starting point for review, not a truth (spec QC-024): recalibrate on your own data with
|
||||
# `simurg qc suggest-thresholds`. For EX-format genotyping chips (provisional, OQ-046).
|
||||
qc:
|
||||
report:
|
||||
callRate: 0.98
|
||||
logRDev: 0.20
|
||||
stainingGreen: 5
|
||||
stainingRed: 5
|
||||
extensionGreen: 5
|
||||
extensionRed: 5
|
||||
hybGreenHighMed: 1
|
||||
hybGreenMediumLow: 1
|
||||
targetRemovalI: 1
|
||||
restoration: null
|
||||
stringency: 2
|
||||
nonSpecificBindingGreen: 1
|
||||
nonSpecificBindingRed: 1
|
||||
npGreen: 5
|
||||
npRed: 5
|
||||
assayControl: 1
|
||||
bgOffset: 3000
|
||||
@@ -0,0 +1,20 @@
|
||||
# Alternative example values from the reference documentation (spec QC-024). Starting point for review only.
|
||||
qc:
|
||||
report:
|
||||
callRate: 0.98
|
||||
logRDev: 0.15
|
||||
stainingGreen: 5
|
||||
stainingRed: 5
|
||||
extensionGreen: 5
|
||||
extensionRed: 5
|
||||
hybGreenHighMed: 1.25
|
||||
hybGreenMediumLow: 1.25
|
||||
targetRemovalI: 5
|
||||
restoration: null
|
||||
stringency: 1.5
|
||||
nonSpecificBindingGreen: 5
|
||||
nonSpecificBindingRed: 5
|
||||
npGreen: 5
|
||||
npRed: 5
|
||||
assayControl: 1
|
||||
bgOffset: 3000
|
||||
@@ -0,0 +1,25 @@
|
||||
# Methylation EPIC profile (spec QC-024). Methylation metrics are computed from v2 on.
|
||||
qc:
|
||||
report:
|
||||
stainingGreen: 5
|
||||
stainingRed: 5
|
||||
extensionGreen: 5
|
||||
extensionRed: 5
|
||||
hybGreenHighMed: 1
|
||||
hybGreenMediumLow: 1
|
||||
targetRemovalI: 1
|
||||
targetRemovalII: 1
|
||||
restoration: 0
|
||||
npGreen: 5
|
||||
npRed: 5
|
||||
bisulfiteConversionIGreen: 1
|
||||
bisulfiteConversionIRed: 1
|
||||
bisulfiteConversionIGreenBg: 1
|
||||
bisulfiteConversionIRedBg: 1
|
||||
bisulfiteConversionII: 1
|
||||
bisulfiteConversionIIBg: 1
|
||||
specificityIGreen: 1
|
||||
specificityIRed: 1
|
||||
specificityII: 1
|
||||
specificityIIBg: 1
|
||||
bgOffset: 3000
|
||||
@@ -0,0 +1,25 @@
|
||||
# Methylation MSA profile (spec QC-024). Methylation metrics are computed from v2 on.
|
||||
qc:
|
||||
report:
|
||||
stainingGreen: 5
|
||||
stainingRed: 5
|
||||
extensionGreen: 5
|
||||
extensionRed: 5
|
||||
hybGreenHighMed: 1
|
||||
hybGreenMediumLow: 1
|
||||
targetRemovalI: 1
|
||||
targetRemovalII: 1
|
||||
restoration: 0
|
||||
npGreen: 2.5
|
||||
npRed: 3
|
||||
bisulfiteConversionIGreen: 1
|
||||
bisulfiteConversionIRed: 1
|
||||
bisulfiteConversionIGreenBg: 0.5
|
||||
bisulfiteConversionIRedBg: 0.5
|
||||
bisulfiteConversionII: 0.5
|
||||
bisulfiteConversionIIBg: 0.5
|
||||
specificityIGreen: 1
|
||||
specificityIRed: 1
|
||||
specificityII: 1
|
||||
specificityIIBg: 1
|
||||
bgOffset: 3000
|
||||
@@ -0,0 +1,113 @@
|
||||
"""Identity and relatedness (spec QC-033, QC-064; plan §9.4).
|
||||
|
||||
* KING-robust kinship between every pair of samples
|
||||
[LIT: Manichaikul A. et al. 2010, Bioinformatics 26:2867, eq. for the within-family (robust) estimator]:
|
||||
φ_ij = [N(Aa,Aa) − 2·N(AA,aa)] / [2·min(N_Aa(i), N_Aa(j))] + 1/2 − [N_Aa(i) + N_Aa(j)] / [4·min(N_Aa(i), N_Aa(j))]
|
||||
with every count taken over the sites called in both samples.
|
||||
* Duplicates: genotype concordance over jointly called sites >= 0.99 (plan §9.4) [OURS threshold].
|
||||
* Mendelian errors of child/parent trios from the sample sheet (Parent1/Parent2).
|
||||
|
||||
Counts are formed with matrix products of per-genotype indicator matrices, so a cohort is O(n²·m) arithmetic in a
|
||||
few BLAS calls. Relationship degrees use KING's published cut-offs (> 0.354 duplicate/MZ twin, > 0.177 first
|
||||
degree, > 0.0884 second, > 0.0442 third).
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
DEGREE_CUTS = ((0.354, "duplicate"), (0.177, "first"), (0.0884, "second"), (0.0442, "third"))
|
||||
DUPLICATE_CONCORDANCE = 0.99
|
||||
|
||||
|
||||
@dataclass
|
||||
class PairStats:
|
||||
kinship: npt.NDArray[np.float64] # (n, n); NaN on the diagonal or where undefined
|
||||
concordance: npt.NDArray[np.float64] # (n, n)
|
||||
n_sites: npt.NDArray[np.int64] # jointly called sites
|
||||
|
||||
|
||||
def degree(phi: float) -> str:
|
||||
for cut, name in DEGREE_CUTS:
|
||||
if phi > cut:
|
||||
return name
|
||||
return "unrelated"
|
||||
|
||||
|
||||
def select_sites(
|
||||
counts: npt.NDArray[np.int64],
|
||||
usable: npt.NDArray[np.bool_],
|
||||
*,
|
||||
min_call_rate: float = 0.95,
|
||||
min_maf: float = 0.01,
|
||||
max_sites: int = 100_000,
|
||||
) -> npt.NDArray[np.intp]:
|
||||
"""High-quality sites from per-locus (no-call, AA, AB, BB) ``counts``: ``usable`` (autosomal, not
|
||||
intensity-only), cohort call rate and MAF filters, thinned evenly to at most ``max_sites``."""
|
||||
called = counts[:, 1:].sum(axis=1)
|
||||
total = counts.sum(axis=1)
|
||||
with np.errstate(invalid="ignore", divide="ignore"):
|
||||
cr = np.where(total > 0, called / total, 0.0)
|
||||
af = np.where(called > 0, (counts[:, 2] + 2 * counts[:, 3]) / (2.0 * called), np.nan)
|
||||
maf = np.minimum(af, 1 - af)
|
||||
idx = np.flatnonzero(usable & (cr >= min_call_rate) & (maf >= min_maf))
|
||||
if len(idx) > max_sites:
|
||||
idx = idx[np.linspace(0, len(idx) - 1, max_sites).round().astype(np.intp)]
|
||||
return idx
|
||||
|
||||
|
||||
def pair_stats(genotypes: npt.NDArray[np.uint8]) -> PairStats:
|
||||
"""KING-robust kinship and concordance for every pair; ``genotypes`` is (sites × samples), 0 = no call."""
|
||||
g = np.asarray(genotypes)
|
||||
aa = (g == 1).astype(np.float32)
|
||||
ab = (g == 2).astype(np.float32)
|
||||
bb = (g == 3).astype(np.float32)
|
||||
c = (g > 0).astype(np.float32)
|
||||
# float32 products are exact while counts stay below 2**24 sites
|
||||
t = np.float64
|
||||
n_het_het = (ab.T @ ab).astype(t)
|
||||
n_opp = (aa.T @ bb).astype(t)
|
||||
n_opp = n_opp + n_opp.T
|
||||
het_i = (ab.T @ c).astype(t) # [i, j] = hets of i over sites called in j
|
||||
het_j = het_i.T
|
||||
joint = (c.T @ c).astype(t)
|
||||
same = (aa.T @ aa + ab.T @ ab + bb.T @ bb).astype(t)
|
||||
mn = np.minimum(het_i, het_j)
|
||||
with np.errstate(invalid="ignore", divide="ignore"):
|
||||
phi = (n_het_het - 2 * n_opp) / (2 * mn) + 0.5 - (het_i + het_j) / (4 * mn)
|
||||
conc = same / joint
|
||||
phi[mn == 0] = np.nan
|
||||
np.fill_diagonal(phi, np.nan)
|
||||
np.fill_diagonal(conc, np.nan)
|
||||
return PairStats(kinship=phi, concordance=conc, n_sites=joint.astype(np.int64))
|
||||
|
||||
|
||||
def mendel_errors(
|
||||
child: npt.NDArray[np.uint8], father: npt.NDArray[np.uint8] | None, mother: npt.NDArray[np.uint8] | None
|
||||
) -> tuple[int, int]:
|
||||
"""(errors, informative sites) for one child with one or two genotyped parents; autosomal sites only.
|
||||
|
||||
A site is informative when the child and every given parent are called. Errors: the child carries an allele
|
||||
that no parent can transmit (duo: opposite homozygotes; trio: child genotype outside the parental combinations).
|
||||
"""
|
||||
parents = [p for p in (father, mother) if p is not None]
|
||||
if not parents:
|
||||
return 0, 0
|
||||
ok = child > 0
|
||||
for p in parents:
|
||||
ok &= p > 0
|
||||
ch = child[ok].astype(np.int8) - 1 # B-allele count 0/1/2
|
||||
ps = [p[ok].astype(np.int8) - 1 for p in parents]
|
||||
if len(ps) == 1:
|
||||
par = ps[0]
|
||||
err = ((ch == 0) & (par == 2)) | ((ch == 2) & (par == 0))
|
||||
else:
|
||||
f, m = ps
|
||||
# B alleles a parent can transmit: hom-A → {0}, het → {0, 1}, hom-B → {1}
|
||||
lo = f // 2 + m // 2
|
||||
hi = (f + 1) // 2 + (m + 1) // 2
|
||||
err = (ch < lo) | (ch > hi)
|
||||
return int(err.sum()), int(ok.sum())
|
||||
@@ -0,0 +1,117 @@
|
||||
"""Two-axis sample QC status and informational flags (spec QC-030, QC-031, QC-033, QC-034).
|
||||
|
||||
* functional status (PASS/FAIL/NA) — call rate and LogR deviation only;
|
||||
* control status (CLEAR/FLAGGED) — control-probe metrics including the assay control;
|
||||
* flags — informational, never change either status.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import math
|
||||
from collections.abc import Mapping, Sequence
|
||||
from dataclasses import dataclass, field
|
||||
from typing import Any
|
||||
|
||||
import numpy as np
|
||||
|
||||
from simurg.qc.thresholds import CONTROL_KEYS, KEY_TO_COLUMN, Thresholds
|
||||
|
||||
CONTAMINATION_FLAG = 0.02 # plan §9.3 [OURS]
|
||||
HET_Z = 3.0
|
||||
ASSAY_CONTROL_MIN = 1.0 # OUT-062
|
||||
|
||||
|
||||
@dataclass
|
||||
class SampleStatus:
|
||||
functional_status: str = "NA"
|
||||
control_status: str = ""
|
||||
control_flag_codes: list[str] = field(default_factory=list)
|
||||
na_metrics: list[str] = field(default_factory=list)
|
||||
flags: list[str] = field(default_factory=list)
|
||||
|
||||
def columns(self) -> dict[str, str]:
|
||||
return {
|
||||
"functional_status": self.functional_status,
|
||||
"control_status": self.control_status,
|
||||
"control_flag_codes": ";".join(self.control_flag_codes),
|
||||
"flags": ";".join(self.flags),
|
||||
}
|
||||
|
||||
|
||||
def _num(v: Any) -> float:
|
||||
try:
|
||||
return float(v)
|
||||
except (TypeError, ValueError):
|
||||
return math.nan
|
||||
|
||||
|
||||
def functional_status(row: Mapping[str, Any] | None, th: Thresholds) -> str:
|
||||
"""QC-030. ``row`` is a sample-summary row; None means no genotyping summary (NA)."""
|
||||
if row is None:
|
||||
return "NA"
|
||||
if row.get("status", "OK") != "OK":
|
||||
return "FAIL"
|
||||
cr, ld = th.get("callRate"), th.get("logRDev")
|
||||
if cr.fails(_num(row.get("autosomal_call_rate"))) or ld.fails(_num(row.get("logr_dev"))):
|
||||
return "FAIL"
|
||||
return "PASS"
|
||||
|
||||
|
||||
def control_status(
|
||||
values: Mapping[str, float] | None, na_metrics: Sequence[str], th: Thresholds
|
||||
) -> tuple[str, list[str]]:
|
||||
"""QC-031: FLAGGED iff an enabled, evaluable control metric is out of range. Codes in QC-004 order."""
|
||||
if values is None:
|
||||
return "", []
|
||||
codes = []
|
||||
for key in CONTROL_KEYS:
|
||||
col = KEY_TO_COLUMN[key]
|
||||
t = th.get(key)
|
||||
if not t.enabled or col not in values or col in na_metrics:
|
||||
continue
|
||||
if t.fails(_num(values[col])):
|
||||
codes.append(key)
|
||||
return ("FLAGGED" if codes else "CLEAR"), codes
|
||||
|
||||
|
||||
def het_zscores(het: Sequence[float]) -> list[float]:
|
||||
"""z-scores of heterozygosity within the batch (NaN where undefined or with < 3 finite values)."""
|
||||
a = np.array([_num(h) for h in het], dtype=np.float64)
|
||||
ok = np.isfinite(a)
|
||||
if ok.sum() < 3:
|
||||
return [math.nan] * len(a)
|
||||
mu, sd = float(a[ok].mean()), float(a[ok].std(ddof=1))
|
||||
if sd == 0:
|
||||
return [0.0 if o else math.nan for o in ok]
|
||||
return [float((x - mu) / sd) if o else math.nan for x, o in zip(a, ok, strict=True)]
|
||||
|
||||
|
||||
def sample_flags(
|
||||
row: Mapping[str, Any],
|
||||
th: Thresholds,
|
||||
*,
|
||||
het_z: float = math.nan,
|
||||
probe_columns: Sequence[str] = (),
|
||||
duplicates: Sequence[str] = (),
|
||||
unexpected_relatedness: bool = False,
|
||||
) -> list[str]:
|
||||
"""QC-033 flags in a fixed order."""
|
||||
flags: list[str] = []
|
||||
ok = row.get("status", "OK") == "OK"
|
||||
if ok and th.get("callRate").fails(_num(row.get("autosomal_call_rate"))):
|
||||
flags.append("low_call_rate")
|
||||
if ok and th.get("logRDev").fails(_num(row.get("logr_dev"))):
|
||||
flags.append("high_lrr_sd")
|
||||
if row.get("sex_check") == "mismatch":
|
||||
flags.append("sex_mismatch")
|
||||
if math.isfinite(het_z) and abs(het_z) > HET_Z:
|
||||
flags.append("het_outlier")
|
||||
c = _num(row.get("contamination"))
|
||||
if math.isfinite(c) and c > CONTAMINATION_FLAG:
|
||||
flags.append("contaminated")
|
||||
if ok and any(_num(row.get(p)) < ASSAY_CONTROL_MIN for p in probe_columns):
|
||||
flags.append("assay_control_fail")
|
||||
flags.extend(f"duplicate_of:{d}" for d in duplicates)
|
||||
if unexpected_relatedness:
|
||||
flags.append("unexpected_relatedness")
|
||||
return flags
|
||||
@@ -0,0 +1,98 @@
|
||||
"""Data-driven threshold suggestion (spec QC-026, CLI-027) [OURS].
|
||||
|
||||
For each metric: robust centre m = median and spread s = 1.4826·MAD over samples, on the log scale for ratio
|
||||
metrics (every control metric and the assay-control norm R); suggestion exp(m − k·s) for ``below`` keys and
|
||||
exp(m + k·s) for ``above`` keys (linear scale for call rate and LogR deviation). At least ``min_samples`` finite
|
||||
values are required; otherwise the key is left unset with a note.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import math
|
||||
from collections.abc import Mapping, Sequence
|
||||
from dataclasses import dataclass
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
from simurg.qc.thresholds import ABOVE_KEYS, Thresholds
|
||||
|
||||
LINEAR_KEYS = {"callRate", "logRDev"}
|
||||
MAD_SCALE = 1.4826
|
||||
|
||||
|
||||
@dataclass
|
||||
class Suggestion:
|
||||
key: str
|
||||
n: int
|
||||
centre: float # on the analysis scale (log for ratio metrics)
|
||||
spread: float
|
||||
value: float | None
|
||||
pct_flagged_current: float
|
||||
pct_flagged_suggested: float
|
||||
note: str = ""
|
||||
|
||||
|
||||
def _flagged(vals: npt.NDArray[np.float64], thr: float | None, above: bool) -> float:
|
||||
if thr is None or len(vals) == 0:
|
||||
return math.nan
|
||||
bad = vals > thr if above else vals < thr
|
||||
return float(100.0 * bad.mean())
|
||||
|
||||
|
||||
def suggest(
|
||||
values: Mapping[str, Sequence[float]],
|
||||
current: Thresholds,
|
||||
*,
|
||||
k: float = 4.0,
|
||||
min_samples: int = 30,
|
||||
) -> list[Suggestion]:
|
||||
out = []
|
||||
for key, raw in values.items():
|
||||
v = np.array([float(x) for x in raw], dtype=np.float64)
|
||||
v = v[np.isfinite(v)]
|
||||
above = key in ABOVE_KEYS
|
||||
linear = key in LINEAR_KEYS
|
||||
cur = current.get(key).value
|
||||
if len(v) < min_samples:
|
||||
out.append(
|
||||
Suggestion(
|
||||
key,
|
||||
len(v),
|
||||
math.nan,
|
||||
math.nan,
|
||||
None,
|
||||
_flagged(v, cur, above),
|
||||
math.nan,
|
||||
f"fewer than {min_samples} samples",
|
||||
)
|
||||
)
|
||||
continue
|
||||
x = v if linear else np.log(np.maximum(v, 1e-12))
|
||||
m = float(np.median(x))
|
||||
s = MAD_SCALE * float(np.median(np.abs(x - m)))
|
||||
t = m + k * s if above else m - k * s
|
||||
val = t if linear else math.exp(t)
|
||||
if key == "callRate":
|
||||
val = min(max(val, 0.0), 1.0)
|
||||
val = max(val, 0.0)
|
||||
out.append(Suggestion(key, len(v), m, s, val, _flagged(v, cur, above), _flagged(v, val, above)))
|
||||
return out
|
||||
|
||||
|
||||
def render_yaml(sugs: Sequence[Suggestion], k: float) -> str:
|
||||
lines = [
|
||||
"# Suggested QC thresholds (simurg qc suggest-thresholds, spec QC-026).",
|
||||
f"# centre = median, spread = {MAD_SCALE}*MAD (log scale for ratio metrics), k = {k}.",
|
||||
"# Review before use: these are starting points computed from your own data.",
|
||||
"qc:",
|
||||
" report:",
|
||||
]
|
||||
for s in sugs:
|
||||
info = (
|
||||
f"n={s.n} m={s.centre:.4g} s={s.spread:.4g} flagged: current={s.pct_flagged_current:.2f}% "
|
||||
f"suggested={s.pct_flagged_suggested:.2f}%"
|
||||
)
|
||||
val = "null" if s.value is None else f"{s.value:.6g}"
|
||||
lines.append(f" {s.key}: {val} # {s.note or info}")
|
||||
return "\n".join(lines) + "\n"
|
||||
@@ -0,0 +1,190 @@
|
||||
"""QC thresholds: YAML schema, validation, profiles and precedence (spec QC-020 … QC-025).
|
||||
|
||||
Precedence, highest first: CLI > ``--config`` YAML > profile > built-in. Every effective value keeps its source.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import math
|
||||
from collections.abc import Mapping
|
||||
from dataclasses import dataclass, field
|
||||
from importlib import resources
|
||||
from pathlib import Path
|
||||
from typing import Any
|
||||
|
||||
import yaml
|
||||
|
||||
from simurg.core.errors import SimurgError
|
||||
from simurg.qc.controls import DEFAULT_BG_OFFSET, METRICS
|
||||
|
||||
FUNCTIONAL_KEYS = ("callRate", "logRDev")
|
||||
CONTROL_KEYS = (*(m.key for m in METRICS), "assayControl")
|
||||
KNOWN_KEYS = (*FUNCTIONAL_KEYS, *CONTROL_KEYS, "bgOffset")
|
||||
ALIASES = {"tgaControl": "assayControl", "call_rate_min": "callRate", "logr_dev_max": "logRDev"}
|
||||
ABOVE_KEYS = {"logRDev"} # flagged when value > threshold; every other key is flagged when value < threshold
|
||||
BUILTIN: dict[str, float | None] = {"callRate": 0.98, "logRDev": 0.20} # QC-023 [EMU]; every other key unset
|
||||
PROFILES = ("genotyping", "genotyping_ex", "genotyping_example", "methylation_epic", "methylation_msa")
|
||||
KEY_TO_COLUMN = {m.key: m.column for m in METRICS} | {"assayControl": "assay_control_norm_r"}
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class Threshold:
|
||||
value: float | None # None: check disabled (QC-022)
|
||||
direction: str # "below" | "above"
|
||||
source: str # builtin | profile:<name> | config | cli
|
||||
|
||||
@property
|
||||
def enabled(self) -> bool:
|
||||
return self.value is not None
|
||||
|
||||
def fails(self, v: float) -> bool:
|
||||
"""True when ``v`` is out of range. Equality passes (QC-020); NaN signal fails (QC-008)."""
|
||||
if self.value is None:
|
||||
return False
|
||||
if v is None or math.isnan(v):
|
||||
return True
|
||||
return v < self.value if self.direction == "below" else v > self.value
|
||||
|
||||
|
||||
@dataclass
|
||||
class Thresholds:
|
||||
items: dict[str, Threshold] = field(default_factory=dict)
|
||||
profile: str = ""
|
||||
|
||||
def get(self, key: str) -> Threshold:
|
||||
return self.items.get(key, Threshold(None, "above" if key in ABOVE_KEYS else "below", "unset"))
|
||||
|
||||
@property
|
||||
def bg_offset(self) -> float:
|
||||
t = self.items.get("bgOffset")
|
||||
return float(t.value) if t is not None and t.value is not None else DEFAULT_BG_OFFSET
|
||||
|
||||
def as_records(self) -> list[dict[str, Any]]:
|
||||
return [
|
||||
{"key": k, "value": t.value, "direction": t.direction, "source": t.source}
|
||||
for k, t in sorted(self.items.items(), key=lambda kv: KNOWN_KEYS.index(kv[0]))
|
||||
]
|
||||
|
||||
|
||||
def _canon_key(key: str) -> str:
|
||||
key = ALIASES.get(key, key)
|
||||
if key not in KNOWN_KEYS:
|
||||
raise SimurgError(
|
||||
"E_QC_THRESHOLD_KEY", f"unknown QC threshold key {key!r}; valid keys: {', '.join(KNOWN_KEYS)}"
|
||||
)
|
||||
return key
|
||||
|
||||
|
||||
def _parse_value(key: str, raw: Any, source: str) -> Threshold:
|
||||
direction = "above" if key in ABOVE_KEYS else "below"
|
||||
if isinstance(raw, Mapping):
|
||||
direction = str(raw.get("direction", direction))
|
||||
if direction not in ("below", "above"):
|
||||
raise SimurgError("E_QC_THRESHOLD_VALUE", f"{key}: direction must be 'below' or 'above', got {direction!r}")
|
||||
raw = raw.get("threshold")
|
||||
if raw is None or (isinstance(raw, str) and raw.strip().lower() in ("", "null", "none", "~")):
|
||||
return Threshold(None, direction, source)
|
||||
if isinstance(raw, bool):
|
||||
raise SimurgError("E_QC_THRESHOLD_VALUE", f"{key}: value must be numeric or null, got {raw!r}")
|
||||
try:
|
||||
v = float(raw)
|
||||
except (TypeError, ValueError):
|
||||
raise SimurgError("E_QC_THRESHOLD_VALUE", f"{key}: value must be numeric or null, got {raw!r}") from None
|
||||
if not math.isfinite(v):
|
||||
raise SimurgError("E_QC_THRESHOLD_VALUE", f"{key}: value must be finite, got {raw!r}")
|
||||
if key == "callRate" and not 0.0 <= v <= 1.0:
|
||||
raise SimurgError("E_QC_THRESHOLD_RANGE", f"callRate must be within [0, 1], got {v}")
|
||||
if v < 0:
|
||||
raise SimurgError("E_QC_THRESHOLD_RANGE", f"{key} must be >= 0, got {v}")
|
||||
return Threshold(v, direction, source)
|
||||
|
||||
|
||||
def parse_report_mapping(doc: Any, where: str) -> tuple[dict[str, Any], str | None]:
|
||||
"""``qc: {profile, report: {...}}`` (or ``QC:``/``Report:``) → (report mapping, profile)."""
|
||||
if not isinstance(doc, Mapping):
|
||||
raise SimurgError("E_QC_CONFIG_STRUCTURE", f"{where}: expected a 'qc:' mapping with a 'report:' section")
|
||||
top = doc.get("qc", doc.get("QC"))
|
||||
if not isinstance(top, Mapping):
|
||||
raise SimurgError("E_QC_CONFIG_STRUCTURE", f"{where}: missing top-level 'qc:' (or 'QC:') mapping")
|
||||
report = top.get("report", top.get("Report"))
|
||||
if report is None:
|
||||
report = {}
|
||||
if not isinstance(report, Mapping):
|
||||
raise SimurgError("E_QC_CONFIG_STRUCTURE", f"{where}: 'report:' (or 'Report:') must be a mapping")
|
||||
profile = top.get("profile")
|
||||
return dict(report), (str(profile) if profile else None)
|
||||
|
||||
|
||||
def load_config(path: Path) -> tuple[dict[str, Any], str | None]:
|
||||
"""Read and structurally validate a threshold YAML (QC-021)."""
|
||||
p = Path(path)
|
||||
if not p.exists():
|
||||
raise SimurgError("E_QC_CONFIG_MISSING", f"QC config file was provided but does not exist: {p}", path=str(p))
|
||||
if p.suffix.lower() not in (".yaml", ".yml"):
|
||||
raise SimurgError("E_QC_CONFIG_EXT", f"QC config file must have a .yaml or .yml extension: {p}", path=str(p))
|
||||
try:
|
||||
doc = yaml.safe_load(p.read_text(encoding="utf-8-sig"))
|
||||
except yaml.YAMLError as e:
|
||||
raise SimurgError("E_QC_CONFIG_PARSE", f"QC config file is not valid YAML: {p}: {e}", path=str(p)) from None
|
||||
return parse_report_mapping(doc, str(p))
|
||||
|
||||
|
||||
def profile_values(name: str) -> dict[str, Any]:
|
||||
"""Shipped profile (QC-024); ``+ffpe`` suffix enables restoration: 1."""
|
||||
base, _, suffix = name.partition("+")
|
||||
if base not in PROFILES or suffix not in ("", "ffpe"):
|
||||
raise SimurgError("E_QC_PROFILE", f"unknown QC profile {name!r}; choose from {', '.join(PROFILES)} or auto")
|
||||
text = resources.files("simurg.qc").joinpath("profiles", f"{base}.yaml").read_text(encoding="utf-8")
|
||||
report, _ = parse_report_mapping(yaml.safe_load(text), f"profile {base}")
|
||||
if suffix == "ffpe":
|
||||
report["restoration"] = 1
|
||||
return report
|
||||
|
||||
|
||||
def auto_profile(array_type: str, assay_format: str = "", manifest_name: str = "") -> str:
|
||||
"""QC-025."""
|
||||
if array_type == "methylation":
|
||||
return "methylation_msa" if "msa" in manifest_name.lower() else "methylation_epic"
|
||||
return "genotyping_ex" if assay_format.strip().upper() == "EX" else "genotyping"
|
||||
|
||||
|
||||
def parse_cli_pairs(pairs: list[str]) -> dict[str, Any]:
|
||||
out: dict[str, Any] = {}
|
||||
for item in pairs:
|
||||
key, sep, val = item.partition("=")
|
||||
if not sep:
|
||||
raise SimurgError("E_QC_THRESHOLD_VALUE", f"--threshold expects KEY=VALUE, got {item!r}")
|
||||
out[key.strip()] = val.strip()
|
||||
return out
|
||||
|
||||
|
||||
def resolve_thresholds(
|
||||
*,
|
||||
config: Path | None = None,
|
||||
profile: str | None = None,
|
||||
cli: Mapping[str, Any] | None = None,
|
||||
array_type: str = "genotyping",
|
||||
assay_format: str = "",
|
||||
manifest_name: str = "",
|
||||
) -> Thresholds:
|
||||
"""Merge built-in < profile < config < CLI (QC-023). A ``profile:`` in the YAML is used when no CLI profile."""
|
||||
layers: list[tuple[str, Mapping[str, Any]]] = [("builtin", BUILTIN)]
|
||||
cfg_report: dict[str, Any] = {}
|
||||
cfg_profile = None
|
||||
if config is not None:
|
||||
cfg_report, cfg_profile = load_config(config)
|
||||
chosen = profile or cfg_profile
|
||||
if chosen == "auto":
|
||||
chosen = auto_profile(array_type, assay_format, manifest_name)
|
||||
if chosen:
|
||||
layers.append((f"profile:{chosen}", profile_values(chosen)))
|
||||
if cfg_report:
|
||||
layers.append(("config", cfg_report))
|
||||
if cli:
|
||||
layers.append(("cli", cli))
|
||||
th = Thresholds(profile=chosen or "")
|
||||
for source, mapping in layers:
|
||||
for raw_key, raw in mapping.items():
|
||||
key = _canon_key(str(raw_key))
|
||||
th.items[key] = _parse_value(key, raw, source)
|
||||
return th
|
||||
Reference in new issue
Block a user