diff --git a/THIRD_PARTY_NOTICES.md b/THIRD_PARTY_NOTICES.md index 4bce67e..6612fa5 100644 --- a/THIRD_PARTY_NOTICES.md +++ b/THIRD_PARTY_NOTICES.md @@ -7,6 +7,8 @@ used only as external subprocess oracles and are not copied or linked. |---|---|---| | [Illumina BeadArrayFiles](https://github.com/Illumina/BeadArrayFiles) — Copyright (c) 2020, Illumina | BSD-2-Clause | File-format layouts (GTC, BPM, EGT), normalization-bin rule, `GTC_File_Format_v5.pdf`; locus-summary formulas (`examples/locus_summary.py`, ported to `simurg/qc/locus_summary.py`); test oracle | | [gtc2vcf / idat2gtc](https://github.com/freeseek/gtc2vcf) — Copyright (C) 2018-2025 Giulio Genovese | MIT | IDAT field codes, BPM locus-entry order, EGT data-block variants, GTC write order, assay-type inference (`infer_assay_type`); BAF-regression numerics (`BAFregress.c`, `simurg/qc/contamination.py`); test oracle | +| [plotly.js](https://github.com/plotly/plotly.js) 2.35.2 — Copyright (c) 2021 Plotly, Inc | MIT | Bundled verbatim in `simurg/qc/report/vendor/plotly.min.js` and inlined into QC reports (licence files alongside) | +| [openpyxl](https://foss.heptapod.net/openpyxl/openpyxl) | MIT | Runtime dependency: `.xlsx` QC table (not vendored) | | [klib `ksort.h`](https://github.com/attractivechaos/klib) as shipped in htslib — Copyright (c) 2008, 2012-2013, 2017-2019 Genome Research Ltd | MIT | Quickselect `ks_ksmall` partition order, ported to `simurg/qc/contamination.py` (`ks_ksmall`) to reproduce BAFregress medians | | [Illumina GTCtoVCF](https://github.com/Illumina/GTCtoVCF) | Apache-2.0 | VCF conversion rules (phase P6) | | [Picard `picard.arrays`](https://github.com/broadinstitute/picard) | MIT | Second format reference | diff --git a/docs/known_unknowns.md b/docs/known_unknowns.md index 2c91553..75fe30f 100644 --- a/docs/known_unknowns.md +++ b/docs/known_unknowns.md @@ -4,6 +4,7 @@ Plan §20 lists the open questions; this file records what the implementation ha | ID | Date | Topic | Finding | Status | |---|---|---|---|---| +| K-020 | 2026-09-28 | QC report scope (P5b) | Implemented: views A–E, linked selection, thresholds modal (apply/download/upload YAML), chip/plate heatmaps, trend table and Tukey box plots, PCA card, xlsx/csv table, privacy options, `--split-by dataset`. Not yet: client-side xlsx export (CSV only), density mode for > 12 000 points (WebGL scatter is used), `--split-by plate|run`, X/Y median R for the sex card (the summary does not carry them; the card shows the sex check instead), methylation datasets (v2). | open (deviations listed) | | K-019 | 2026-09-28 | GSA demo relatedness | The 36-sample GSA v4 PGx demo set contains duplicate/replicate pairs (KING φ ≈ 0.50, concordance ≥ 0.9999) and first-degree relatives (φ ≈ 0.25); without a pedigree they are flagged `duplicate_of:`/`unexpected_relatedness` as designed. Useful as a small relatedness fixture. | documented | | K-018 | 2026-09-28 | Locus-summary ChiTest_P100 | Against BeadArrayFiles `locus_summary.py` all columns are bit-identical on 36 GSA samples (683 054 loci) except ChiTest_P100, which differs by ≤ 3.6e-15 relative in ~15 % of loci: scipy 1.18's `ndtr` matches neither libm nor a cephes erf/erfc port bit-for-bit. Accepted (p-value). | accepted | | K-017 | 2026-09-28 | BAF-regression median quirk | gtc2vcf `BAFregress` (1) compares AF with the threshold as a float (0.1f), and (2) for an even number of AA/BB calls averages the quickselect result with whatever element klib `ks_ksmall` leaves just below it, not the lower middle value. Both reproduced (`--median reference`, default): 36/36 GSA samples identical to the plugin at its printed precision; `--median exact` gives the true median. SimurgArray uses autosomal sites only in pipelines; the plugin uses every VCF record. | resolved | diff --git a/pyproject.toml b/pyproject.toml index e6fdc56..174d9fb 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -16,6 +16,7 @@ dependencies = [ "rich>=13", # progress / human-readable output "pyyaml>=6", # --config files "pyarrow>=15", # canonical columnar model and Parquet outputs (plan §5) + "openpyxl>=3.1", # QC table .xlsx (QC-045) ] [project.optional-dependencies] diff --git a/src/simurg/cli/main.py b/src/simurg/cli/main.py index 7c4ccd0..6fec45e 100644 --- a/src/simurg/cli/main.py +++ b/src/simurg/cli/main.py @@ -602,6 +602,62 @@ def qc_loci_cmd( _qc_cohort("loci", out, bpm, gtcs, sheet, egt, threads, force, log_json, debug, quiet) +@qc_app.command("report") +def qc_report_cmd( + data: Annotated[list[str], typer.Option("--data", help="Dataset folder(s): one, a comma list, or a parent.")], + out: Annotated[Path, typer.Option("--out", help="Output folder.")] = Path("."), + label: Annotated[ + list[str] | None, typer.Option("--label", help="Dataset labels (list or discovery order).") + ] = None, + config: Annotated[Path | None, typer.Option("--config", help="QC threshold YAML (QC-020).")] = 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, + output_format: Annotated[str, typer.Option("--output-format", help="xlsx|csv|both")] = "xlsx", + report_name: Annotated[str | None, typer.Option("--report-name", help="Fixed output file stem.")] = None, + drop_metadata: Annotated[ + list[str] | None, typer.Option("--drop-metadata", help="Metadata columns to drop.") + ] = None, + keep_metadata: Annotated[list[str] | None, typer.Option("--keep-metadata", help="Metadata allowlist.")] = None, + no_metadata: Annotated[bool, typer.Option("--no-metadata", help="Embed no metadata columns.")] = False, + redact_paths: Annotated[bool, typer.Option("--redact-paths", help="Blank sample_folder.")] = False, + split_by: Annotated[str, typer.Option("--split-by", help="none|dataset")] = "none", + log_json: LogJson = False, + debug: Debug = False, + quiet: Quiet = False, +) -> None: + """Self-contained HTML QC report and QC table from dataset folders (QC-040 … QC-063).""" + from simurg.pipeline.qc_report import QcReportOptions, run_qc_report + from simurg.qc.thresholds import parse_cli_pairs, resolve_thresholds + + rep = _reporter("qc report", None, 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 = QcReportOptions( + data=list(data), + out=out, + labels=list(label) if label else None, + output_format=output_format, + report_name=report_name, + drop_metadata=list(drop_metadata or []), + keep_metadata=list(keep_metadata) if keep_metadata else None, + no_metadata=no_metadata, + redact_paths=redact_paths, + split_by=split_by, + thresholds=th, + ) + return run_qc_report(opts, rep) + + run_guarded(rep, body) + + @qc_app.command("suggest-thresholds") def qc_suggest_cmd( data: Annotated[list[Path], typer.Option("--data", help="Dataset folder(s) with controls_qc.tsv.")], diff --git a/src/simurg/pipeline/qc_report.py b/src/simurg/pipeline/qc_report.py new file mode 100644 index 0000000..913ff82 --- /dev/null +++ b/src/simurg/pipeline/qc_report.py @@ -0,0 +1,623 @@ +"""``simurg qc report`` — QC table (xlsx/csv) and a self-contained HTML report (spec QC-040 … QC-063, CLI-026). + +The report reads only dataset folders (control tables, optionally the sample summary), never raw data. Data is +embedded as gzip-compressed JSON (decoded in the page with the browser's DecompressionStream) next to the bundled +Plotly.js (MIT), so the file works offline. +""" + +from __future__ import annotations + +import base64 +import csv +import datetime as _dt +import gzip +import hashlib +import io +import json +import math +import re +from collections import Counter +from collections.abc import Mapping, Sequence +from dataclasses import dataclass, field +from importlib import resources +from pathlib import Path +from typing import Any + +from simurg import PRODUCT_NAME, __version__ +from simurg.core.atomic import atomic_write_bytes, atomic_write_text +from simurg.core.errors import ExitCode, SimurgError +from simurg.core.log import Reporter +from simurg.qc.controls import COMPAT_METRIC_NAMES, METRICS, beads_from_definitions, excluded_from_report, metrics_for +from simurg.qc.status import control_status, functional_status +from simurg.qc.thresholds import ABOVE_KEYS, FUNCTIONAL_KEYS, KEY_TO_COLUMN, KNOWN_KEYS, Thresholds, resolve_thresholds + +REQUIRED = (("controls_raw.tsv", "controls.raw_metrics.csv"), ("controls_qc.tsv", "controls.qc_metrics.csv")) +SUMMARY_NAMES = ("sample_summary.tsv", "sample_summary.csv", "gt_sample_summary.csv") +QC_META = ("imaging_date", "sample_plate", "sample_well", "scanner_id", "scanner_version") +SUMMARY_CORE = { + "sample_id", "sample_name", "sample_folder", "status", "fail_reason", "autosomal_call_rate", "call_rate", + "logr_dev", "sex_est", "sex_reported", "sex_check", "gc10", "gc50", "n_calls", "n_nocalls", "n_intensity_only", + "p05_x", "p50_x", "p95_x", "p05_y", "p50_y", "p95_y", "het_rate", "contamination", "functional_status", + "control_status", "control_flag_codes", "flags", +} # fmt: skip +COMPAT_QC = { + "SampleId": "sample_id", "ImagingDate": "imaging_date", "SamplePlate": "sample_plate", + "SampleWell": "sample_well", "ScannerId": "scanner_id", "ScannerVersion": "scanner_version", +} # fmt: skip +COMPAT_RAW = { + "SampleId": "sample_id", "Category": "category", "Beadtype": "bead_address", "Control": "control_name", + "Color": "color_label", "GreenIntensity": "green_intensity", "RedIntensity": "red_intensity", +} # fmt: skip +COMPAT_SUMMARY = { + "Sample ID": "sample_id", "Sample Name": "sample_name", "Sample Folder": "sample_folder", + "Autosomal Call Rate": "autosomal_call_rate", "Call Rate": "call_rate", "Log R Ratio Std Dev": "logr_dev", + "Sex Estimate": "sex_est", +} # fmt: skip +FORMULAS = { + "restoration": "Restoration Green / (Bkg Green + offset)", + "staining_green": "Biotin (High) Green / Biotin (Bkg) Green", + "staining_red": "DNP (High) Red / DNP (Bkg) Red", + "extension_green": "lowest Extension (C, G) Green / highest Extension (A, T) Green", + "extension_red": "lowest Extension (A, T) Red / highest Extension (C, G) Red", + "hyb_high_med": "Hyb (High) Green / Hyb (Medium) Green", + "hyb_med_low": "Hyb (Medium) Green / Hyb (Low) Green", + "target_removal_1": "(Bkg Green + offset) / Target Removal Green", + "target_removal_2": "(Bkg Green + offset) / Target Removal 2 Green", + "stringency": "String (PM) Red / String (MM) Red", + "nsb_green": "(Bkg Green + offset) / highest NSB Green", + "nsb_red": "(Bkg Red + offset) / highest NSB Red", + "np_green": "lowest NP (C, G) Green / highest NP (A, T) Green", + "np_red": "lowest NP (A, T) Red / highest NP (C, G) Red", + "assay_control_norm_r": "normalized R of the assay-control locus (lowest over listed probes)", + "autosomal_call_rate": "calls / loci on chromosomes 1-22", + "logr_dev": "standard deviation of autosomal Log R ratio", +} +# QC-052 chip formats: samples per chip -> (rows, columns) +CHIP_FORMATS = {4: (4, 1), 8: (8, 1), 12: (12, 1), 16: (8, 2), 24: (12, 2), 48: (12, 4), 96: (12, 8)} +POS_RE = re.compile(r"^(\d+)_R(\d{2})C(\d{2})$") +LABEL_RE = re.compile(r"[^A-Za-z0-9_.-]") + + +# ------------------------------------------------------------------------------------------------ reading +def read_table(path: Path) -> list[dict[str, str]]: + """QC-061: UTF-8 (BOM allowed), delimiter by extension, RFC-4180 quoting; errors name file and line.""" + delim = "\t" if path.suffix.lower() == ".tsv" else "," + try: + text = path.read_text(encoding="utf-8-sig") + except UnicodeDecodeError as e: + raise SimurgError("E_QC_DATASET_INCOMPLETE", f"{path}: not UTF-8 ({e})", path=str(path)) from None + reader = csv.reader(io.StringIO(text), delimiter=delim) + try: + header = next(reader) + except StopIteration: + return [] + rows = [] + for i, rec in enumerate(reader, start=2): + if not rec: + continue + if len(rec) != len(header): + raise SimurgError( + "E_QC_DATASET_INCOMPLETE", + f"{path}: line {i} has {len(rec)} columns, header has {len(header)}", + path=str(path), + ) + rows.append(dict(zip(header, rec, strict=True))) + return rows + + +def _rename(rows: list[dict[str, str]], mapping: Mapping[str, str]) -> list[dict[str, str]]: + return [{mapping.get(k, k): v for k, v in r.items()} for r in rows] + + +def _find(folder: Path, names: Sequence[str]) -> Path | None: + for n in names: + if (folder / n).exists(): + return folder / n + return None + + +def _is_dataset(folder: Path) -> bool: + return all(_find(folder, pair) for pair in REQUIRED) + + +def resolve_datasets(data: Sequence[str]) -> list[Path]: + """QC-041: one folder, a comma list, or one parent folder searched recursively (sorted by path).""" + items = [Path(p.strip()) for d in data for p in str(d).split(",") if p.strip()] + if not items: + raise SimurgError("E_QC_NO_DATASETS", "No valid dataset folders detected.") + for p in items: + if not p.exists(): + raise SimurgError("E_QC_NO_DATASETS", f"Input data folder does not exist: {p}", path=str(p)) + if len(items) == 1 and not _is_dataset(items[0]): + found = sorted({f.parent for f in items[0].rglob("controls_qc.tsv")} | { + f.parent for f in items[0].rglob("controls.qc_metrics.csv")}) # fmt: skip + found = [f for f in found if _is_dataset(f)] + if not found: + raise SimurgError( + "E_QC_NO_DATASETS", + "No valid dataset folders detected. Each folder must include at least: controls_raw.tsv, " + "controls_qc.tsv", + ) + return found + for p in items: + missing = [pair[0] for pair in REQUIRED if not _find(p, pair)] + if missing: + raise SimurgError( + "E_QC_DATASET_INCOMPLETE", f"Dataset {p.name} is missing required file(s): {', '.join(missing)}" + ) + return items + + +def assign_labels(folders: Sequence[Path], labels: Sequence[str] | None, rep: Reporter) -> list[str]: + """QC-042.""" + if labels: + labels = [x for lab in labels for x in str(lab).split(",") if x != ""] + if len(labels) != len(folders): + raise SimurgError( + "E_QC_LABEL_COUNT", + f"Number of labels ({len(labels)}) does not match number of datasets ({len(folders)}).", + ) + raw = list(labels) + else: + seen: Counter[str] = Counter() + raw = [] + for f in folders: + seen[f.name] += 1 + raw.append(f.name if seen[f.name] == 1 else f"{f.name}_{seen[f.name]}") + out = [] + for lab in raw: + clean = LABEL_RE.sub("_", lab) + if clean != lab: + rep.warn("W_QC_LABEL_SANITIZED", f"Label sanitized: {lab} → {clean}") + out.append(clean) + dups = [k for k, c in Counter(out).items() if c > 1] + if dups: + raise SimurgError( + "E_QC_LABEL_DUPLICATE", f"Dataset labels are not unique after sanitisation: {', '.join(dups)}" + ) + return out + + +def _num(v: Any) -> float | None: + if v is None or v == "": + return None + try: + f = float(v) + except (TypeError, ValueError): + return None + return None if math.isnan(f) or math.isinf(f) else f + + +@dataclass +class Dataset: + label: str + path: Path + type: str + samples: list[dict[str, Any]] + probes: list[dict[str, str]] # raw probe columns: key, name, category, channel, color + metrics: list[str] # derived metric columns present + metadata: list[str] + has_summary: bool + chip_format: tuple[int, int] | None + profile: str = "" + files: dict[str, str] = field(default_factory=dict) + + +def _chip_format(folder: Path) -> tuple[int, int] | None: + for mf in sorted(folder.glob("run_manifest*.json")): + try: + doc = json.loads(mf.read_text(encoding="utf-8")) + except (OSError, ValueError): + continue + for inp in doc.get("inputs", []): + m = re.search(r"-(\d+)v\d", str(inp.get("name", ""))) + if inp.get("role") in ("bpm", "csv_manifest") and m and int(m.group(1)) in CHIP_FORMATS: + return CHIP_FORMATS[int(m.group(1))] + return None + + +def load_dataset(folder: Path, label: str, rep: Reporter) -> Dataset: + raw_p = _find(folder, REQUIRED[0]) + qc_p = _find(folder, REQUIRED[1]) + assert raw_p is not None and qc_p is not None + raw = _rename(read_table(raw_p), COMPAT_RAW) + qc = _rename(read_table(qc_p), COMPAT_QC) + sum_p = _find(folder, SUMMARY_NAMES) + summary = _rename(read_table(sum_p), COMPAT_SUMMARY) if sum_p else [] + for rows, name in ((raw, raw_p), (qc, qc_p)): + if rows and "sample_id" not in rows[0]: + raise SimurgError("E_QC_DATASET_INCOMPLETE", f"{name}: no sample_id column", path=str(name)) + cats = {r.get("category", "").casefold() for r in raw} + array_type = (qc[0].get("array_type") if qc else "") or ( + "methylation" if cats & {"bisulfite conversion i", "bisulfite conversion ii"} else "genotyping" + ) + compat_metric = {v: k for k, v in COMPAT_METRIC_NAMES.items()} + qc = [{compat_metric.get(k, k): v for k, v in r.items()} for r in qc] + metrics = [m.column for m in metrics_for(array_type) if qc and m.column in qc[0]] + if qc and "assay_control_norm_r" in qc[0]: + metrics.append("assay_control_norm_r") + # raw probes, wide + defs: dict[str, dict[str, str]] = {} + name_count = ( + Counter((r["control_name"], r["color_label"]) for r in raw if r["sample_id"] == raw[0]["sample_id"]) + if raw + else Counter() + ) + names_dup = Counter(n for n, _ in name_count) + wide: dict[str, dict[str, float | None]] = {} + for r in raw: + bead = beads_from_definitions( + [ + { + "address": r.get("bead_address") or "0", + "category": r["category"], + "color": r["color_label"], + "name": r["control_name"], + } + ] + )[0] + if excluded_from_report(bead): + continue + base = r["control_name"] if names_dup[r["control_name"]] <= 1 else f"{r['control_name']} {r['color_label']}" + for ch, col in (("Green", "green_intensity"), ("Red", "red_intensity")): + key = f"{base}_{ch.lower()}" + defs.setdefault(key, {"key": key, "name": base, "category": r["category"], "channel": ch, + "color": r["color_label"]}) # fmt: skip + wide.setdefault(r["sample_id"], {})[key] = _num(r[col]) + by_sum = {r["sample_id"]: r for r in summary} + by_qc = {r["sample_id"]: r for r in qc} + ids = list(by_qc) + unmatched = sorted(set(by_sum) ^ set(by_qc)) if summary else [] + if unmatched: + rep.warn( + "W_QC_UNMATCHED_SAMPLES", + f"{label}: {len(unmatched)} sample ID(s) in only one of the QC files, e.g. {', '.join(unmatched[:20])}", + ) + ids += [s for s in by_sum if s not in by_qc] + meta_cols: list[str] = [] + if summary: + cols = list(summary[0]) + tail = cols[cols.index("flags") + 1 :] if "flags" in cols else [c for c in cols if c not in SUMMARY_CORE] + meta_cols = [c for c in tail if not c.endswith("_norm_r")] + meta_cols.append("sample_folder") + meta_cols += [c for c in QC_META] + probe_cols = [c for c in (summary[0] if summary else {}) if c.endswith("_norm_r")] + samples = [] + for sid in ids: + s_row, q_row = by_sum.get(sid, {}), by_qc.get(sid, {}) + m = POS_RE.match(sid) + rec: dict[str, Any] = { + "sample_id": sid, + "sample_name": s_row.get("sample_name", sid), + "barcode": m.group(1) if m else "", + "row": int(m.group(2)) if m else None, + "col": int(m.group(3)) if m else None, + "status": s_row.get("status", "OK" if summary else ""), + "autosomal_call_rate": _num(s_row.get("autosomal_call_rate")), + "call_rate": _num(s_row.get("call_rate")), + "logr_dev": _num(s_row.get("logr_dev")), + "sex_est": s_row.get("sex_est", ""), + "sex_check": s_row.get("sex_check", ""), + "het_rate": _num(s_row.get("het_rate")), + "contamination": _num(s_row.get("contamination")), + "flags": s_row.get("flags", ""), + "na_metrics": q_row.get("na_metrics", ""), + "has_summary": bool(s_row), + "m": {c: _num(q_row.get(c)) for c in metrics}, + "raw": wide.get(sid, {}), + "meta": {c: (s_row.get(c) if c in s_row else q_row.get(c, "")) or "" for c in meta_cols}, + } + if "assay_control_norm_r" not in q_row and probe_cols: + vals = [v for v in (_num(s_row.get(c)) for c in probe_cols) if v is not None] + rec["m"]["assay_control_norm_r"] = min(vals) if vals else None + if "assay_control_norm_r" not in metrics: + metrics.append("assay_control_norm_r") + samples.append(rec) + files = {p.name: hashlib.sha256(p.read_bytes()).hexdigest() for p in (raw_p, qc_p, sum_p) if p is not None} + return Dataset( + label=label, + path=folder, + type=array_type, + samples=samples, + probes=list(defs.values()), + metrics=metrics, + metadata=meta_cols, + has_summary=bool(summary), + chip_format=_chip_format(folder), + files=files, + ) + + +def apply_privacy( + datasets: Sequence[Dataset], drop: Sequence[str], keep: Sequence[str] | None, no_metadata: bool, redact_paths: bool +) -> bool: + """QC-046; returns True when any metadata remains.""" + drop_f = {d.casefold() for x in drop for d in x.split(",") if d} + keep_f = {k.casefold() for x in keep for k in x.split(",") if k} if keep else None + any_left = False + for ds in datasets: + cols = [] + for c in ds.metadata: + if no_metadata or c.casefold() in drop_f or (keep_f is not None and c.casefold() not in keep_f): + continue + cols.append(c) + if redact_paths: + for s in ds.samples: + s["meta"]["sample_folder"] = "" + for s in ds.samples: + s["meta"] = {c: s["meta"].get(c, "") for c in cols} + ds.metadata = sorted(cols, key=str.casefold) + any_left |= any(any(s["meta"].values()) for s in ds.samples) + return any_left + + +# ------------------------------------------------------------------------------------------------ status +def evaluate(ds: Dataset, th: Thresholds) -> None: + for s in ds.samples: + row = {"status": s["status"], "autosomal_call_rate": s["autosomal_call_rate"], "logr_dev": s["logr_dev"]} + row = {k: (math.nan if v is None else v) for k, v in row.items()} + s["functional_status"] = functional_status(row if s["has_summary"] else None, th) + vals = {k: (math.nan if v is None else v) for k, v in s["m"].items()} + cs, codes = control_status(vals, [c for c in s["na_metrics"].split(";") if c], th) + s["control_status"], s["control_flag_codes"] = cs, ";".join(codes) + + +# ------------------------------------------------------------------------------------------------ table +def table_rows(datasets: Sequence[Dataset]) -> tuple[list[str], list[dict[str, Any]]]: + """QC-045 Samples sheet: fixed columns, functional metrics, derived metrics, raw intensities, metadata.""" + metrics: list[str] = [] + probes: list[str] = [] + meta: list[str] = [] + for ds in datasets: + metrics += [m for m in ds.metrics if m not in metrics and m != "assay_control_norm_r"] + probes += [p["key"] for p in ds.probes if p["key"] not in probes] + meta += [c for c in ds.metadata if c not in meta] + meta = sorted(meta, key=str.casefold) + cols = ["dataset", "sample_id", "sample_name", "functional_status", "control_status", "control_flag_codes", + "na_metrics", "autosomal_call_rate", "call_rate", "logr_dev", "sex_est", "assay_control_norm_r", + *metrics, *probes, *meta] # fmt: skip + rows = [] + for ds in datasets: + for s in ds.samples: + r: dict[str, Any] = {k: s.get(k) for k in cols[:11]} + r["dataset"] = ds.label + r["assay_control_norm_r"] = s["m"].get("assay_control_norm_r") + r.update({m: s["m"].get(m) for m in metrics}) + r.update({p: s["raw"].get(p) for p in probes}) + r.update({c: s["meta"].get(c, "") for c in meta}) + rows.append(r) + return cols, rows + + +def write_csv_table(path: Path, cols: Sequence[str], rows: Sequence[Mapping[str, Any]]) -> None: + buf = io.StringIO() + w = csv.writer(buf, lineterminator="\n") + w.writerow(cols) + for r in rows: + w.writerow(["" if r.get(c) is None else r.get(c) for c in cols]) + atomic_write_text(path, buf.getvalue()) + + +def write_xlsx_table( + path: Path, + cols: Sequence[str], + rows: Sequence[Mapping[str, Any]], + th: Thresholds, + datasets: Sequence[Dataset], + generated: str, +) -> None: + from openpyxl import Workbook + from openpyxl.styles import Font, PatternFill + + fills = { + "PASS": PatternFill("solid", fgColor="BBD7F0"), + "FAIL": PatternFill("solid", fgColor="F5C28E"), + "FLAGGED": PatternFill("solid", fgColor="FFE08A"), + "CLEAR": PatternFill("solid", fgColor="BBD7F0"), + "func_out": PatternFill("solid", fgColor="F5C28E"), + "ctrl_out": PatternFill("solid", fgColor="FFF1C2"), + "code": PatternFill("solid", fgColor="FFE08A"), + } + col_key = {v: k for k, v in KEY_TO_COLUMN.items()} | {"autosomal_call_rate": "callRate", "logr_dev": "logRDev"} + wb = Workbook() + ws = wb.active + ws.title = "Samples" + ws.append(list(cols)) + for c in ws[1]: + c.font = Font(bold=True) + for r in rows: + ws.append(["" if r.get(c) is None else r.get(c) for c in cols]) + cur = ws.max_row + for j, c in enumerate(cols, start=1): + v = r.get(c) + cell = ws.cell(cur, j) + if c in ("functional_status", "control_status") and v in fills: + cell.fill = fills[str(v)] + elif c == "control_flag_codes" and v: + cell.fill = fills["code"] + elif c in col_key and isinstance(v, int | float): + t = th.get(col_key[c]) + if t.enabled and t.fails(float(v)): + cell.fill = fills["func_out" if col_key[c] in FUNCTIONAL_KEYS else "ctrl_out"] + ws.freeze_panes = "A2" + ws.auto_filter.ref = ws.dimensions + t2 = wb.create_sheet("Thresholds") + t2.append(["key", "description", "threshold", "direction", "enabled", "source"]) + for k in KNOWN_KEYS: + t = th.get(k) + desc = FORMULAS.get( + KEY_TO_COLUMN.get(k, {"callRate": "autosomal_call_rate", "logRDev": "logr_dev"}.get(k, "")), "" + ) + t2.append( + [k, desc, t.value, t.direction, t.enabled, "default" if t.source in ("builtin", "unset") else t.source] + ) + t3 = wb.create_sheet("Datasets") + t3.append(["label", "path", "type", "n_samples", "profile"]) + for ds in datasets: + t3.append([ds.label, str(ds.path), ds.type, len(ds.samples), th.profile]) + t4 = wb.create_sheet("About") + t4.append(["software", f"{PRODUCT_NAME} {__version__}"]) + t4.append(["generated", generated]) + for ds in datasets: + for name, h in sorted(ds.files.items()): + t4.append([f"{ds.label}/{name}", h]) + buf = io.BytesIO() + wb.save(buf) + atomic_write_bytes(path, buf.getvalue()) + + +# ------------------------------------------------------------------------------------------------ html +def _asset(name: str) -> str: + return resources.files("simurg.qc.report").joinpath(name).read_text(encoding="utf-8") + + +def payload(datasets: Sequence[Dataset], th: Thresholds, generated: str, metadata_notice: bool) -> dict[str, Any]: + metric_defs = [ + { + "column": m.column, + "key": m.key, + "formula": FORMULAS.get(m.column, ""), + "direction": "above" if m.key in ABOVE_KEYS else "below", + } + for m in METRICS + ] + [ + { + "column": "assay_control_norm_r", + "key": "assayControl", + "formula": FORMULAS["assay_control_norm_r"], + "direction": "below", + } + ] + return { + "product": PRODUCT_NAME, + "version": __version__, + "generated": generated, + "metadataNotice": metadata_notice, + "profile": th.profile, + "thresholds": { + k: {"value": th.get(k).value, "direction": th.get(k).direction, "source": th.get(k).source} + for k in KNOWN_KEYS + }, + "functionalKeys": [ + { + "key": "callRate", + "column": "autosomal_call_rate", + "label": "Autosomal Call Rate", + "formula": FORMULAS["autosomal_call_rate"], + }, + {"key": "logRDev", "column": "logr_dev", "label": "LRR SD", "formula": FORMULAS["logr_dev"]}, + ], + "metricDefs": metric_defs, + "datasets": [ + { + "label": ds.label, + "type": ds.type, + "hasSummary": ds.has_summary, + "chipFormat": list(ds.chip_format) if ds.chip_format else None, + "metrics": ds.metrics, + "probes": ds.probes, + "metadata": ds.metadata, + "samples": ds.samples, + } + for ds in datasets + ], + } + + +def render_html(data: Mapping[str, Any]) -> str: + blob = json.dumps(data, sort_keys=True, ensure_ascii=False, separators=(",", ":"), allow_nan=False) + packed = base64.b64encode(gzip.compress(blob.encode("utf-8"), compresslevel=9, mtime=0)).decode("ascii") + page = _asset("template.html").replace("__TITLE__", f"{PRODUCT_NAME} QC Report").replace("__DATA__", packed) + for marker, text in ( + ("/*__STYLE__*/", _asset("report.css")), + ("/*__APP__*/", _asset("report.js")), + ("/*__PLOTLY__*/", _asset("vendor/plotly.min.js").replace(" Any: + if isinstance(o, float): + return None if math.isnan(o) or math.isinf(o) else o + if isinstance(o, dict): + return {k: _sanitize_json(v) for k, v in o.items()} + if isinstance(o, list): + return [_sanitize_json(v) for v in o] + return o + + +def run_qc_report(opts: QcReportOptions, rep: Reporter) -> int: + if opts.output_format not in ("xlsx", "csv", "both"): + raise SimurgError("E_QC_FORMAT", "Invalid --output-format value. Only csv, xlsx or both are accepted.") + if opts.split_by not in ("none", "dataset"): + raise SimurgError("E_QC_FORMAT", "--split-by supports none and dataset in this version") + th = opts.thresholds or resolve_thresholds() + folders = resolve_datasets(opts.data) + labels = assign_labels(folders, opts.labels, rep) + datasets = [load_dataset(f, lab, rep) for f, lab in zip(folders, labels, strict=True)] + types = {ds.type for ds in datasets} + if len(types) > 1: + by_type = {t: [ds.label for ds in datasets if ds.type == t] for t in sorted(types)} + raise SimurgError( + "E_QC_MIXED_TYPES", + "Genotyping and methylation datasets cannot be mixed: " + + "; ".join(f"{t}: {', '.join(v)}" for t, v in by_type.items()), + ) + if "methylation" in types: + raise SimurgError("E_QC_METHYLATION", "methylation datasets in the report arrive with v2 (QC-062)") + total = sum(len(ds.samples) for ds in datasets) + if total > 20000: + rep.warn("W_QC_TOO_FEW", f"{total} samples in total: the report may be slow to open (QC-057)") + notice = apply_privacy(datasets, opts.drop_metadata, opts.keep_metadata, opts.no_metadata, opts.redact_paths) + for ds in datasets: + evaluate(ds, th) + now = opts.now or _dt.datetime.now(_dt.UTC) + generated = now.strftime("%Y-%m-%dT%H:%M:%SZ") + day = now.strftime("%Y_%m_%d") + opts.out.mkdir(parents=True, exist_ok=True) + groups = [(None, datasets)] if opts.split_by == "none" else [(ds.label, [ds]) for ds in datasets] + written = [] + for suffix, group in groups: + stem_r = opts.report_name or f"qc_report_{day}" + stem_t = opts.report_name or f"qc_table_{day}" + if opts.report_name: + stem_r, stem_t = f"{opts.report_name}", f"{opts.report_name}_table" + if suffix: + stem_r, stem_t = f"{stem_r}_{suffix}", f"{stem_t}_{suffix}" + cols, rows = table_rows(group) + if opts.output_format in ("xlsx", "both"): + write_xlsx_table(opts.out / f"{stem_t}.xlsx", cols, rows, th, group, generated) + written.append(opts.out / f"{stem_t}.xlsx") + if opts.output_format in ("csv", "both"): + write_csv_table(opts.out / f"{stem_t}.csv", cols, rows) + written.append(opts.out / f"{stem_t}.csv") + html = render_html(_sanitize_json(payload(group, th, generated, notice))) + atomic_write_text(opts.out / f"{stem_r}.html", html) + written.append(opts.out / f"{stem_r}.html") + if opts.split_by != "none": + links = "".join(f'
--no-metadata, --drop-metadata or --keep-metadata to remove columns).0?[0]:[]);if(o.enter().append("g").classed(h.containerClassName,!0).style("cursor","pointer"),o.exit().each((function(){n.select(this).selectAll("g."+h.headerGroupClassName).each(a)})).remove(),0!==r.length){var l=o.selectAll("g."+h.headerGroupClassName).data(r,p);l.enter().append("g").classed(h.headerGroupClassName,!0);for(var c=s.ensureSingle(o,"g",h.dropdownButtonGroupClassName,(function(t){t.style("pointer-events","all")})),u=0;u 90&&i.log("Long binary search..."),f-1},e.sorterAsc=function(t,e){return t-e},e.sorterDes=function(t,e){return e-t},e.distinctVals=function(t){var r,n=t.slice();for(n.sort(e.sorterAsc),r=n.length-1;r>-1&&n[r]===o;r--);for(var i,a=n[r]-n[0]||1,s=a/(r||1)/1e4,l=[],c=0;c<=r;c++){var u=n[c],h=u-i;void 0===i?(l.push(u),i=u):h>s&&(a=Math.min(a,h),l.push(u),i=u)}return{vals:l,minDiff:a}},e.roundUp=function(t,e,r){for(var n,i=0,a=e.length-1,o=0,s=r?0:1,l=r?1:0,c=r?Math.ceil:Math.floor;i0&&(n=1),r&&n)return t.sort(e)}return n?t:t.reverse()},e.findIndexOfMin=function(t,e){e=e||a;for(var r,n=1/0,i=0;il?r.y-l:0;return Math.sqrt(c*c+h*h)}for(var p=f(c);p;){if((c+=p+r)>h)return;p=f(c)}for(p=f(h);p;){if(c>(h-=p+r))return;p=f(h)}return{min:c,max:h,len:h-c,total:u,isClosed:0===c&&h===u&&Math.abs(n.x-i.x)<.1&&Math.abs(n.y-i.y)<.1}},e.findPointOnPath=function(t,e,r,n){for(var i,a,o,s=(n=n||{}).pathLength||t.getTotalLength(),l=n.tolerance||.001,c=n.iterationLimit||30,u=t.getPointAtLength(0)[r]>t.getPointAtLength(s)[r]?-1:1,h=0,f=0,p=s;h
/i;e.BR_TAG_ALL=/
/gi;var b=/(^|[\s"'])style\s*=\s*("([^"]*);?"|'([^']*);?')/i,w=/(^|[\s"'])href\s*=\s*("([^"]*)"|'([^']*)')/i,T=/(^|[\s"'])target\s*=\s*("([^"\s]*)"|'([^'\s]*)')/i,k=/(^|[\s"'])popup\s*=\s*("([\w=,]*)"|'([\w=,]*)')/i;function A(t,e){if(!t)return null;var r=t.match(e),n=r&&(r[3]||r[4]);return n&&C(n)}var M=/(^|;)\s*color:/;e.plainText=function(t,e){for(var r=void 0!==(e=e||{}).len&&-1!==e.len?e.len:1/0,n=void 0!==e.allowedTags?e.allowedTags:["br"],i=t.split(v),a=[],o="",s=0,l=0;l
"+l;e.text=c}(t,o,r,c):"log"===u?function(t,e,r,n,a){var o=t.dtick,l=e.x,c=t.tickformat,u="string"==typeof o&&o.charAt(0);if("never"===a&&(a=""),n&&"L"!==u&&(o="L3",u="L"),c||"L"===u)e.text=wt(Math.pow(10,l),t,a,n);else if(i(o)||"D"===u&&s.mod(l+.01,1)<.1){var h=Math.round(l),f=Math.abs(h),p=t.exponentformat;"power"===p||_t(p)&&bt(h)?(e.text=0===h?1:1===h?"10":"10"+(h>1?"":z)+f+"",e.fontSize*=1.25):("e"===p||"E"===p)&&f>2?e.text="1"+p+(h>0?"+":z)+f:(e.text=wt(Math.pow(10,l),t,"","fakehover"),"D1"===o&&"y"===t._id.charAt(0)&&(e.dy-=e.fontSize/6))}else{if("D"!==u)throw"unrecognized dtick "+String(o);e.text=String(Math.round(Math.pow(10,s.mod(l,1)))),e.fontSize*=.75}if("D1"===t.dtick){var d=String(e.text).charAt(0);"0"!==d&&"1"!==d||("y"===t._id.charAt(0)?e.dx-=e.fontSize/4:(e.dy+=e.fontSize/2,e.dx+=(t.range[1]>t.range[0]?1:-1)*e.fontSize*(l<0?.5:.25)))}}(t,o,0,c,g):"category"===u?function(t,e){var r=t._categories[Math.round(e.x)];void 0===r&&(r=""),e.text=String(r)}(t,o):"multicategory"===u?function(t,e,r){var n=Math.round(e.x),i=t._categories[n]||[],a=void 0===i[1]?"":String(i[1]),o=void 0===i[0]?"":String(i[0]);r?e.text=o+" - "+a:(e.text=a,e.text2=o)}(t,o,r):Rt(t)?function(t,e,r,n,i){if("radians"!==t.thetaunit||r)e.text=wt(e.x,t,i,n);else{var a=e.x/180;if(0===a)e.text="0";else{var o=function(t){function e(t,e){return Math.abs(t-e)<=1e-6}var r=function(t){for(var r=1;!e(Math.round(t*r)/r,t);)r*=10;return r}(t),n=t*r,i=Math.abs(function t(r,n){return e(n,0)?r:t(n,r%n)}(n,r));return[Math.round(n/i),Math.round(r/i)]}(a);if(o[1]>=100)e.text=wt(s.deg2rad(e.x),t,i,n);else{var l=e.x<0;1===o[1]?1===o[0]?e.text="π":e.text=o[0]+"π":e.text=["",o[0],"","⁄","",o[1],"","π"].join(""),l&&(e.text=z+e.text)}}}}(t,o,r,c,g):function(t,e,r,n,i){"never"===i?i="":"all"===t.showexponent&&Math.abs(e.x/t.dtick)<1e-6&&(i="hide"),e.text=wt(e.x,t,i,n)}(t,o,0,c,g),n||(t.tickprefix&&!m(t.showtickprefix)&&(o.text=t.tickprefix+o.text),t.ticksuffix&&!m(t.showticksuffix)&&(o.text+=t.ticksuffix)),t.labelalias&&t.labelalias.hasOwnProperty(o.text)){var y=t.labelalias[o.text];"string"==typeof y&&(o.text=y)}return("boundaries"===t.tickson||t.showdividers)&&(o.xbnd=[f(o.x-.5),f(o.x+t.dtick-.5)]),o},Z.hoverLabelText=function(t,e,r){r&&(t=s.extendFlat({},t,{hoverformat:r}));var n=s.isArrayOrTypedArray(e)?e[0]:e,i=s.isArrayOrTypedArray(e)?e[1]:void 0;if(void 0!==i&&i!==n)return Z.hoverLabelText(t,n,r)+" - "+Z.hoverLabelText(t,i,r);var a="log"===t.type&&n<=0,o=Z.tickText(t,t.c2l(a?-n:n),"hover").text;return a?0===n?"0":z+o:o};var xt=["f","p","n","μ","m","","k","M","G","T"];function _t(t){return"SI"===t||"B"===t}function bt(t){return t>14||t<-15}function wt(t,e,r,n){var a=t<0,o=e._tickround,l=r||e.exponentformat||"B",c=e._tickexponent,u=Z.getTickFormat(e),h=e.separatethousands;if(n){var f={exponentformat:l,minexponent:e.minexponent,dtick:"none"===e.showexponent?e.dtick:i(t)&&Math.abs(t)||1,range:"none"===e.showexponent?e.range.map(e.r2d):[0,t||1]};yt(f),o=(Number(f._tickround)||0)+4,c=f._tickexponent,e.hoverformat&&(u=e.hoverformat)}if(u)return e._numFormat(u)(t).replace(/-/g,z);var p,d=Math.pow(10,-o)/2;if("none"===l&&(c=0),(t=Math.abs(t))
")):x=f.textLabel;var C={x:f.traceCoordinate[0],y:f.traceCoordinate[1],z:f.traceCoordinate[2],data:b._input,fullData:b,curveNumber:b.index,pointNumber:T};d.appendArrayPointValue(C,b,T),t._module.eventData&&(C=b._module.eventData(C,f,b,{},T));var L={points:[C]};if(e.fullSceneLayout.hovermode){var I=[];d.loneHover({trace:b,x:(.5+.5*v[0]/v[3])*s,y:(.5-.5*v[1]/v[3])*l,xLabel:k.xLabel,yLabel:k.yLabel,zLabel:k.zLabel,text:x,name:u.name,color:d.castHoverOption(b,T,"bgcolor")||u.color,borderColor:d.castHoverOption(b,T,"bordercolor"),fontFamily:d.castHoverOption(b,T,"font.family"),fontSize:d.castHoverOption(b,T,"font.size"),fontColor:d.castHoverOption(b,T,"font.color"),nameLength:d.castHoverOption(b,T,"namelength"),textAlign:d.castHoverOption(b,T,"align"),hovertemplate:h.castOption(b,T,"hovertemplate"),hovertemplateLabels:h.extendFlat({},C,k),eventData:[C]},{container:n,gd:r,inOut_bbox:I}),C.bbox=I[0]}f.distance<5&&(f.buttons||w)?r.emit("plotly_click",L):r.emit("plotly_hover",L),this.oldEventData=L}else d.loneUnhover(n),this.oldEventData&&r.emit("plotly_unhover",this.oldEventData),this.oldEventData=void 0;e.drawAnnotations(e)},k.recoverContext=function(){var t=this;t.glplot.dispose();var e=function(){t.glplot.gl.isContextLost()?requestAnimationFrame(e):t.initializeGLPlot()?t.plot.apply(t,t.plotArgs):h.error("Catastrophic and unrecoverable WebGL error. Context lost.")};requestAnimationFrame(e)};var M=["xaxis","yaxis","zaxis"];function S(t,e,r){for(var n=t.fullSceneLayout,i=0;i<3;i++){var a=M[i],o=a.charAt(0),s=n[a],l=e[o],c=e[o+"calendar"],u=e["_"+o+"length"];if(h.isArrayOrTypedArray(l))for(var f,p=0;p<(u||l.length);p++)if(h.isArrayOrTypedArray(l[p]))for(var d=0;d
");y.text(_).attr("data-unformatted",_).call(u.convertToTspans,t),v=c.bBox(y.node())}y.attr("transform",i(-3,8-v.height)),g.insert("rect",".static-attribution").attr({x:-v.width-6,y:-v.height-3,width:v.width+6,height:v.height+3,fill:"rgba(255, 255, 255, 0.75)"});var b=1;v.width+6>x&&(b=x/(v.width+6));var w=[n.l+n.w*p.x[1],n.t+n.h*(1-p.y[0])];g.attr("transform",i(w[0],w[1])+a(b))}},e.updateFx=function(t){for(var e=t._fullLayout,r=e._subplots[f],n=0;n
");_.text(T).attr("data-unformatted",T).call(h.convertToTspans,t),b=u.bBox(_.node())}_.attr("transform",a(-3,8-b.height)),x.insert("rect",".static-attribution").attr({x:-b.width-6,y:-b.height-3,width:b.width+6,height:b.height+3,fill:"rgba(255, 255, 255, 0.75)"});var k=1;b.width+6>w&&(k=w/(b.width+6));var A=[n.l+n.w*f.x[1],n.t+n.h*(1-f.y[0])];x.attr("transform",a(A[0],A[1])+o(k))}},e.updateFx=function(t){for(var e=t._fullLayout,r=e._subplots[p],n=0;n