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'
  • {p.name}
  • ' for p in written if p.suffix == ".html") + atomic_write_text( + opts.out / "index.html", + f"{PRODUCT_NAME} QC reports\n", + ) + rep.info(f"QC report for {total} sample(s) in {len(datasets)} dataset(s): " + ", ".join(p.name for p in written)) + return ExitCode.OK diff --git a/src/simurg/qc/report/__init__.py b/src/simurg/qc/report/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/src/simurg/qc/report/report.css b/src/simurg/qc/report/report.css new file mode 100644 index 0000000..d8c6931 --- /dev/null +++ b/src/simurg/qc/report/report.css @@ -0,0 +1,47 @@ +:root { --blue:#3B75AF; --orange:#E8871E; --amber:#F2C14E; --grey:#9A9A9A; --bg:#fafafa; --line:#ddd; } +* { box-sizing: border-box; } +body { margin:0; font:14px/1.4 system-ui,-apple-system,"Segoe UI",Roboto,sans-serif; color:#222; background:var(--bg); } +header { padding:10px 16px; background:#fff; border-bottom:1px solid var(--line); } +header .title { font-size:20px; font-weight:600; } +.muted { color:#666; font-size:13px; } +.notice { margin-top:6px; padding:6px 10px; background:#FFF6DA; border:1px solid var(--amber); border-radius:4px; font-size:13px; } +#layout { display:flex; min-height:calc(100vh - 70px); } +nav { width:190px; flex:none; background:#fff; border-right:1px solid var(--line); padding:8px; } +nav button { display:block; width:100%; text-align:left; margin:2px 0; padding:8px 10px; border:0; background:none; border-radius:4px; cursor:pointer; font-size:14px; } +nav button.active { background:#E3EDF7; font-weight:600; } +nav button:hover { background:#eef3f8; } +main { flex:1; min-width:0; padding:12px 16px; } +.controls { display:flex; flex-wrap:wrap; gap:12px; align-items:end; margin-bottom:10px; } +.controls label { display:flex; flex-direction:column; font-size:12px; color:#555; } +select, input[type=text], input[type=number] { font-size:13px; padding:3px 5px; border:1px solid #bbb; border-radius:3px; } +button.btn { font-size:13px; padding:4px 10px; border:1px solid #999; background:#fff; border-radius:3px; cursor:pointer; } +button.btn:hover { background:#f0f0f0; } +.grid2 { display:grid; grid-template-columns:repeat(2, minmax(0,1fr)); gap:12px; } +@media (max-width: 1100px) { .grid2 { grid-template-columns:1fr; } } +.card { background:#fff; border:1px solid var(--line); border-radius:6px; padding:8px 10px; } +.card h3 { margin:0 0 2px; font-size:15px; } +.card .help { font-size:12px; color:#666; margin-bottom:4px; } +.plot { width:100%; height:300px; } +.plot.small { height:220px; } +.tabs { display:flex; gap:4px; flex-wrap:wrap; margin:6px 0; } +.tabs button { border:1px solid #bbb; background:#fff; padding:3px 8px; border-radius:3px; cursor:pointer; font-size:12px; } +.tabs button.on { background:#E3EDF7; border-color:var(--blue); } +table.dt { border-collapse:collapse; width:100%; font-size:12px; background:#fff; } +table.dt th, table.dt td { border:1px solid #e3e3e3; padding:3px 5px; white-space:nowrap; } +table.dt th { background:#f3f3f3; cursor:pointer; position:sticky; top:0; } +table.dt tfoot input { width:100%; min-width:50px; font-size:11px; } +table.dt tr.sel td { outline:1px solid var(--blue); background:#EAF2FB !important; } +.tablewrap { overflow:auto; max-height:520px; border:1px solid var(--line); } +td.PASS, td.CLEAR { background:#D6E6F5; } td.FAIL { background:#F9D3AC; } td.FLAGGED, td.code { background:#FFE7A3; } +td.fout { background:#F9D3AC; } td.cout { background:#FFF1C2; } +.pager { display:flex; gap:8px; align-items:center; margin:6px 0; font-size:12px; flex-wrap:wrap; } +.dropdown { position:relative; display:inline-block; } +.dropdown .menu { position:absolute; z-index:10; background:#fff; border:1px solid #bbb; padding:6px; max-height:300px; overflow:auto; min-width:180px; box-shadow:0 2px 6px rgba(0,0,0,.15); } +.modal { position:fixed; inset:0; background:rgba(0,0,0,.35); display:flex; align-items:flex-start; justify-content:center; padding-top:40px; z-index:50; } +.modal[hidden] { display:none; } +.modal-box { background:#fff; border-radius:6px; padding:14px 18px; width:min(1000px, 96vw); max-height:88vh; overflow:auto; } +.modal-box table { border-collapse:collapse; width:100%; font-size:13px; } +.modal-box td, .modal-box th { padding:3px 6px; border-bottom:1px solid #eee; text-align:left; } +.err { color:#B00020; font-size:12px; } +.legend-sw { display:inline-block; width:12px; height:12px; vertical-align:middle; margin:0 4px 0 10px; border:1px solid #888; } +.summarytable th { text-align:center; } diff --git a/src/simurg/qc/report/report.js b/src/simurg/qc/report/report.js new file mode 100644 index 0000000..1f9684c --- /dev/null +++ b/src/simurg/qc/report/report.js @@ -0,0 +1,645 @@ +/* SimurgArray QC report (spec QC-047 … QC-062). Plain JS + bundled Plotly.js; no network access. */ +(async function () { + "use strict"; + const COL = { PASS: "#3B75AF", CLEAR: "#3B75AF", FAIL: "#E8871E", FLAGGED: "#F2C14E", NA: "#9A9A9A" }; + const PALETTES = { + "Okabe-Ito": ["#0072B2", "#E69F00", "#009E73", "#CC79A7", "#56B4E9", "#D55E00", "#F0E442", "#000000"], + "Tol bright": ["#4477AA", "#EE6677", "#228833", "#CCBB44", "#66CCEE", "#AA3377", "#BBBBBB"], + "Tol muted": ["#332288", "#88CCEE", "#44AA99", "#117733", "#999933", "#DDCC77", "#CC6677", "#882255", "#AA4499"], + }; + const CATEGORIES = ["Staining", "Extension", "Target Removal", "Hybridization", "Restoration", "Stringency", + "Non-Specific Binding", "Non-Polymorphic"]; + const EXPECT = { + "Staining|Green": "Biotin (High) should be strong; DNP (High/Bkg) and Biotin (Bkg) should be near background.", + "Staining|Red": "DNP (High) should be strong; Biotin (High/Bkg) and DNP (Bkg) should be near background.", + "Extension|Green": "Extension (C) and (G) should be strong; (A) and (T) near background.", + "Extension|Red": "Extension (A) and (T) should be strong; (C) and (G) near background.", + "Target Removal|Green": "Should stay near background (target was washed off).", + "Target Removal|Red": "Not used; near background.", + "Hybridization|Green": "High > Medium > Low, all clearly above background.", + "Hybridization|Red": "Not used; near background.", + "Restoration|Green": "Above background for restored (FFPE) samples.", + "Restoration|Red": "Not used; near background.", + "Stringency|Green": "Not used; near background.", + "Stringency|Red": "PM should be clearly higher than MM.", + "Non-Specific Binding|Green": "All NSB probes should be low (background).", + "Non-Specific Binding|Red": "All NSB probes should be low (background).", + "Non-Polymorphic|Green": "NP (C) and (G) should be strong; (A) and (T) low.", + "Non-Polymorphic|Red": "NP (A) and (T) should be strong; (C) and (G) low.", + }; + const CONFIG = { displaylogo: false, responsive: true, + modeBarButtonsToAdd: ["select2d", "lasso2d"], toImageButtonOptions: { format: "png", scale: 2 } }; + + // ------------------------------------------------------------------------------------------ data + const b64 = document.getElementById("data").textContent.trim(); + const bin = Uint8Array.from(atob(b64), (c) => c.charCodeAt(0)); + const D = JSON.parse(await new Response(new Blob([bin]).stream().pipeThrough(new DecompressionStream("gzip"))).text()); + const TH = JSON.parse(JSON.stringify(D.thresholds)); + const MDEF = Object.fromEntries(D.metricDefs.map((m) => [m.column, m])); + let cur = 0; // current dataset index (shared by views A-C) + let selection = null; // Set of sample_ids or null (QC-051) + const nTotal = D.datasets.reduce((a, d) => a + d.samples.length, 0); + + document.getElementById("hdrinfo").textContent = + `${D.product} ${D.version} · generated ${D.generated} · ${D.datasets.length} dataset(s) · ${nTotal} sample(s)` + + (D.profile ? ` · profile ${D.profile}` : ""); + document.getElementById("metanotice").hidden = !D.metadataNotice; + if (nTotal > 20000) { + const n = document.getElementById("scalenotice"); + n.hidden = false; + n.textContent = `${nTotal} samples: scatter plots use WebGL and may offer a density mode (QC-057).`; + } + + // ------------------------------------------------------------------------------------------ helpers + const $ = (tag, attrs = {}, ...kids) => { + const e = document.createElement(tag); + for (const [k, v] of Object.entries(attrs)) { + if (k === "class") e.className = v; + else if (k.startsWith("on")) e.addEventListener(k.slice(2), v); + else if (v !== undefined && v !== null) e.setAttribute(k, v); + } + for (const k of kids.flat()) if (k !== null && k !== undefined) e.append(k instanceof Node ? k : String(k)); + return e; + }; + const select = (label, options, value, onchange) => { + const s = $("select", { onchange: (e) => onchange(e.target.value) }, + options.map(([v, t]) => { const o = $("option", { value: v }, t); if (String(v) === String(value)) o.selected = true; return o; })); + return $("label", {}, label, s); + }; + const fmt = (v, d = 4) => (v === null || v === undefined || Number.isNaN(v) ? "" : typeof v === "number" ? +v.toFixed(d) : v); + const ds = () => D.datasets[cur]; + const val = (s, col) => (col in s.m ? s.m[col] : s[col] !== undefined ? s[col] : s.meta[col]); + const trunc = (x) => { x = String(x ?? ""); return x.length > 200 ? x.slice(0, 200) + "…" : x; }; + const fails = (key, v) => { + const t = TH[key]; + if (!t || t.value === null || t.value === undefined || t.value === "") return false; + if (v === null || v === undefined || Number.isNaN(v)) return true; + return t.direction === "below" ? v < t.value : v > t.value; + }; + const enabled = (key) => TH[key] && TH[key].value !== null && TH[key].value !== undefined && TH[key].value !== ""; + function recompute() { // QC-030 / QC-031, client side (QC-055 Apply) + for (const d of D.datasets) for (const s of d.samples) { + if (!s.has_summary) s.functional_status = "NA"; + else if (s.status !== "OK") s.functional_status = "FAIL"; + else s.functional_status = fails("callRate", s.autosomal_call_rate) || fails("logRDev", s.logr_dev) ? "FAIL" : "PASS"; + const na = new Set((s.na_metrics || "").split(";").filter(Boolean)); + const codes = D.metricDefs.filter((m) => enabled(m.key) && m.column in s.m && !na.has(m.column) && fails(m.key, s.m[m.column])).map((m) => m.key); + s.control_status = codes.length ? "FLAGGED" : "CLEAR"; + s.control_flag_codes = codes.join(";"); + } + } + recompute(); + const eligible = (d) => d.metadata.filter((c) => { // QC-056 + const u = new Set(d.samples.map((s) => s.meta[c]).filter((x) => x !== "" && x !== undefined)); + return u.size >= 2 && u.size <= 20; + }); + const hover = (s, extra) => { + const d = ds(); + const lines = [`${s.sample_id}`, `dataset: ${d.label}`]; + if (s.barcode) lines.push(`barcode: ${s.barcode} · R${String(s.row).padStart(2, "0")}C${String(s.col).padStart(2, "0")}`); + if (s.autosomal_call_rate !== null && s.autosomal_call_rate !== undefined) lines.push(`autosomal call rate: ${fmt(s.autosomal_call_rate)}`); + if (extra) lines.push(extra); + d.metadata.slice(0, 12).forEach((c) => { if (s.meta[c]) lines.push(`${c}: ${trunc(s.meta[c])}`); }); + return lines.join("
    "); + }; + function sorted(samples, key) { + const arr = samples.slice(); + if (key === "barcode") arr.sort((a, b) => (a.barcode + a.sample_id).localeCompare(b.barcode + b.sample_id)); + else arr.sort((a, b) => { + const x = val(a, key), y = val(b, key); + if (typeof x === "number" && typeof y === "number") return x - y; + if (x === null || x === undefined) return 1; + if (y === null || y === undefined) return -1; + return String(x).localeCompare(String(y)); + }); + return arr; + } + const scatterType = (n) => (n > 12000 ? "scattergl" : "scatter"); + const sortOptions = (d) => [["barcode", "Sentrix barcode"], ...(d.hasSummary ? [["autosomal_call_rate", "Autosomal call rate"], ["logr_dev", "LRR SD"]] : []), + ...d.metrics.map((m) => [m, m]), ...d.metadata.map((c) => [c, c])]; + const datasetSelect = (onchange) => select("Dataset", D.datasets.map((d, i) => [i, `${d.label} (${d.samples.length})`]), cur, (v) => { cur = +v; selection = null; onchange(); }); + function download(name, text, mime) { + const a = $("a", { href: URL.createObjectURL(new Blob([text], { type: mime })), download: name }); + document.body.append(a); a.click(); a.remove(); + } + + // ------------------------------------------------------------------------------------------ linked selection + const linked = []; // {div, samples(order)} for every rendered scatter + function registerLinked(div, order) { + linked.push({ div, order }); + div.on("plotly_selected", (ev) => { + if (!ev || !ev.points) return; + selection = new Set(ev.points.map((p) => p.customdata).filter(Boolean)); + applySelection(div); + }); + div.on("plotly_deselect", () => { selection = null; applySelection(div); }); + div.on("plotly_doubleclick", () => { selection = null; applySelection(div); }); + } + function applySelection(origin) { + for (const { div, order } of linked) { + if (!document.body.contains(div) || div === origin || !div.data) continue; + const sp = div.data.map((tr) => (selection && tr.customdata ? tr.customdata.map((id, i) => (selection.has(id) ? i : -1)).filter((i) => i >= 0) : null)); + Plotly.restyle(div, { selectedpoints: sp }); + } + if (table) table.refresh(); + } + + // ------------------------------------------------------------------------------------------ box statistics (QC-054) + function quant(sortedV, p) { + if (!sortedV.length) return null; + const h = (sortedV.length - 1) * p, lo = Math.floor(h); + return sortedV[lo] + (h - lo) * ((sortedV[lo + 1] ?? sortedV[lo]) - sortedV[lo]); + } + function boxStats(values) { + const v = values.filter((x) => typeof x === "number" && !Number.isNaN(x)).sort((a, b) => a - b); + if (!v.length) return null; + const q1 = quant(v, 0.25), med = quant(v, 0.5), q3 = quant(v, 0.75), iqr = q3 - q1; + const lo = v.find((x) => x >= q1 - 1.5 * iqr), hi = [...v].reverse().find((x) => x <= q3 + 1.5 * iqr); + return { q1, med, q3, lo, hi, out: v.filter((x) => x < lo || x > hi), n: v.length, + mean: v.reduce((a, b) => a + b, 0) / v.length, min: v[0], max: v[v.length - 1], + sd: v.length > 1 ? Math.sqrt(v.reduce((a, b) => a + (b - v.reduce((c, e) => c + e, 0) / v.length) ** 2, 0) / (v.length - 1)) : 0 }; + } + function boxTraces(groups, horizontal = true) { // explicit Tukey rule: whiskers to extreme points within 1.5 IQR + const names = [], q1 = [], med = [], q3 = [], lo = [], hi = [], ox = [], oy = []; + for (const [name, values] of groups) { + const b = boxStats(values); if (!b) continue; + names.push(name); q1.push(b.q1); med.push(b.med); q3.push(b.q3); lo.push(b.lo); hi.push(b.hi); + b.out.forEach((x) => { ox.push(x); oy.push(name); }); + } + const box = { type: "box", q1, median: med, q3, lowerfence: lo, upperfence: hi, marker: { color: "#3B75AF" }, name: "", showlegend: false, hoverinfo: "x" }; + if (horizontal) box.y = names; else box.x = names; + box.orientation = horizontal ? "h" : "v"; + const pts = { type: "scatter", mode: "markers", marker: { color: "#3B75AF", size: 5 }, showlegend: false, name: "outliers" }; + if (horizontal) { pts.x = ox; pts.y = oy; } else { pts.x = oy; pts.y = ox; } + return [box, pts]; + } + const thLine = (key, horizontalAxis) => { + if (!enabled(key)) return { shapes: [], annotations: [] }; + const t = TH[key], txt = `${key === "callRate" || key === "logRDev" ? "Fail" : "Flag"} if ${t.direction} ${t.value}`; + const s = horizontalAxis + ? { type: "line", xref: "x", yref: "paper", x0: t.value, x1: t.value, y0: 0, y1: 1, line: { dash: "dash", color: "#C0392B" } } + : { type: "line", xref: "paper", yref: "y", x0: 0, x1: 1, y0: t.value, y1: t.value, line: { dash: "dash", color: "#C0392B" } }; + const a = horizontalAxis + ? { xref: "x", yref: "paper", x: t.value, y: 0.02, text: txt, showarrow: false, xanchor: "left", yanchor: "bottom", font: { size: 11, color: "#C0392B" } } + : { xref: "paper", yref: "y", x: 1, y: t.value, text: txt, showarrow: false, xanchor: "right", yanchor: "bottom", font: { size: 11, color: "#C0392B" } }; + return { shapes: [s], annotations: [a] }; + }; + + // ------------------------------------------------------------------------------------------ view A: dashboard + const dashState = { sort: "barcode", palette: "Okabe-Ito", facet: "", tab: null }; + function renderDash() { + const root = document.getElementById("v-dash"); root.replaceChildren(); + linked.length = 0; + const d = ds(); + root.append($("div", { class: "controls" }, datasetSelect(renderDash), + select("Sort plots by", sortOptions(d), dashState.sort, (v) => { dashState.sort = v; renderDash(); }), + select("Colour palette", Object.keys(PALETTES).map((k) => [k, k]), dashState.palette, (v) => { dashState.palette = v; renderDash(); }), + select("Facet by", [["", "None"], ...eligible(d).map((c) => [c, c])], dashState.facet, (v) => { dashState.facet = v; renderDash(); }))); + let order = sorted(d.samples, dashState.sort); + let facetShapes = [], facetNotes = []; + if (dashState.facet) { // QC-056: up to 12 facets, the rest grouped as "Other" + const counts = {}; + order.forEach((s) => { const v = s.meta[dashState.facet] || "(empty)"; counts[v] = (counts[v] || 0) + 1; }); + const top = Object.keys(counts).sort().slice(0, 11), lvl = (s) => { const v = s.meta[dashState.facet] || "(empty)"; return top.includes(v) || Object.keys(counts).length <= 12 ? v : "Other"; }; + order = order.slice().sort((a, b) => lvl(a).localeCompare(lvl(b))); + let start = 0; + for (let i = 1; i <= order.length; i++) if (i === order.length || lvl(order[i]) !== lvl(order[i - 1])) { + facetNotes.push({ x: order[Math.floor((start + i - 1) / 2)].sample_id, y: 1, xref: "x", yref: "paper", text: lvl(order[start]), showarrow: false, yanchor: "bottom", font: { size: 11 } }); + if (i < order.length) facetShapes.push({ type: "line", xref: "x", yref: "paper", x0: i - 0.5, x1: i - 0.5, y0: 0, y1: 1, line: { color: "#888", dash: "dot" } }); + start = i; + } + if (Object.keys(counts).length > 12) root.append($("div", { class: "muted" }, `${Object.keys(counts).length} levels: showing 11 and "Other".`)); + } + const ids = order.map((s) => s.sample_id); + if (d.hasSummary) { + const c = $("div", { class: "card" }, $("h3", {}, "Autosomal Call Rate"), $("div", { class: "plot small" })); + root.append(c); + Plotly.newPlot(c.lastChild, boxTraces([[d.label, d.samples.map((s) => s.autosomal_call_rate)]]), + { margin: { l: 90, r: 20, t: 20, b: 40 }, ...thLine("callRate", true), xaxis: { title: "autosomal call rate" } }, CONFIG); + } + const cats = CATEGORIES.filter((c) => d.probes.some((p) => p.category === c)).concat( + [...new Set(d.probes.map((p) => p.category))].filter((c) => !CATEGORIES.includes(c))); + const cards = []; + for (const cat of cats) for (const ch of ["Green", "Red"]) cards.push({ cat, ch }); + if (d.metrics.includes("assay_control_norm_r")) cards.push({ cat: "Assay control", ch: "" }); + const big = d.samples.length > 1000; // QC-057: tabs per category + if (big && !dashState.tab) dashState.tab = cards[0] && cards[0].cat; + if (big) root.append($("div", { class: "tabs" }, [...new Set(cards.map((c) => c.cat))].map((c) => + $("button", { class: c === dashState.tab ? "on" : "", onclick: () => { dashState.tab = c; renderDash(); } }, c)))); + const grid = $("div", { class: "grid2" }); root.append(grid); + const pal = PALETTES[dashState.palette]; + for (const { cat, ch } of cards) { + if (big && cat !== dashState.tab) continue; + const title = ch ? `${cat} – ${ch}` : cat; + const help = ch ? EXPECT[`${cat}|${ch}`] || "" : "Normalized R of the assay-control locus; should be ≥ the threshold."; + const card = $("div", { class: "card" }, $("h3", {}, title), $("div", { class: "help" }, help), $("div", { class: "plot" })); + grid.append(card); + let traces; + if (ch) { + traces = d.probes.filter((p) => p.category === cat && p.channel === ch).map((p, i) => ({ + type: scatterType(order.length), mode: "markers", name: p.name, x: ids, y: order.map((s) => s.raw[p.key] ?? null), + customdata: ids, text: order.map((s) => hover(s, `${p.name} ${ch}: ${s.raw[p.key] ?? ""}`)), hoverinfo: "text", + marker: { color: pal[i % pal.length], size: 6 } })); + } else { + traces = [{ type: scatterType(order.length), mode: "markers", name: "assay control", x: ids, y: order.map((s) => s.m.assay_control_norm_r), + customdata: ids, text: order.map((s) => hover(s, `norm R: ${fmt(s.m.assay_control_norm_r)}`)), hoverinfo: "text", marker: { color: pal[0], size: 6 } }]; + } + const tl = ch ? { shapes: [], annotations: [] } : thLine("assayControl", false); + Plotly.newPlot(card.lastChild, traces, { margin: { l: 60, r: 10, t: 24, b: 30 }, dragmode: "select", hovermode: "closest", + xaxis: { type: "category", showticklabels: false, title: "samples" }, yaxis: { rangemode: "tozero", title: ch ? "raw intensity" : "norm R" }, + showlegend: true, legend: { orientation: "h", y: -0.15 }, shapes: [...facetShapes, ...tl.shapes], annotations: [...facetNotes, ...tl.annotations] }, CONFIG); + registerLinked(card.lastChild, order); + } + } + + // ------------------------------------------------------------------------------------------ table (QC-050) + let table = null; + function makeTable(host, d) { + const cols = [ + ["dataset", "Dataset", (s) => d.label], ["sample_id", "Sample ID", (s) => s.sample_id], + ...(d.hasSummary ? [["functional_status", "Functional QC Status", (s) => s.functional_status]] : []), + ["control_status", "Control QC Status", (s) => s.control_status], + ["control_flag_codes", "Control QC Flag Codes", (s) => s.control_flag_codes], + ...(d.hasSummary ? [["autosomal_call_rate", "Autosomal Call Rate", (s) => s.autosomal_call_rate], ["call_rate", "Overall Call Rate", (s) => s.call_rate], + ["logr_dev", "LRR SD", (s) => s.logr_dev], ["sex_est", "Sex Estimate", (s) => ({ M: "Male", F: "Female" })[s.sex_est] || "Unknown"], + ["contamination", "Contamination", (s) => s.contamination], ["flags", "Flags", (s) => s.flags]] : []), + ...d.metrics.map((m) => [m, m, (s) => s.m[m]]), + ...d.metadata.map((c) => [c, c, (s) => s.meta[c]]), + ]; + const key = Object.fromEntries(D.metricDefs.map((m) => [m.column, m.key]).concat([["autosomal_call_rate", "callRate"], ["logr_dev", "logRDev"]])); + const st = { sort: null, asc: true, page: 0, per: 10, search: "", filt: {}, hidden: new Set(), rowsel: new Set() }; + const wrap = $("div"); + host.append(wrap); + function rows() { + let r = d.samples.filter((s) => !selection || selection.has(s.sample_id)); + if (st.search) { const q = st.search.toLowerCase(); r = r.filter((s) => cols.some(([, , f]) => String(f(s) ?? "").toLowerCase().includes(q))); } + for (const [c, f] of Object.entries(st.filt)) { + const col = cols.find((x) => x[0] === c); if (!col) continue; + r = r.filter((s) => { const v = col[2](s); + if (typeof f === "string") return !f || String(v ?? "").toLowerCase().includes(f.toLowerCase()); + if (f.min !== "" && f.min !== undefined && !(v >= +f.min)) return false; + if (f.max !== "" && f.max !== undefined && !(v <= +f.max)) return false; return true; }); + } + if (st.sort) { const col = cols.find((x) => x[0] === st.sort); r = r.slice().sort((a, b) => { + const x = col[2](a), y = col[2](b); const c = typeof x === "number" && typeof y === "number" ? x - y : String(x ?? "").localeCompare(String(y ?? "")); + return st.asc ? c : -c; }); } + return r; + } + function exportRows() { + const r = st.rowsel.size ? d.samples.filter((s) => st.rowsel.has(s.sample_id)) : rows(); + const vis = cols.filter((c) => !st.hidden.has(c[0])); + const esc = (v) => { v = v === null || v === undefined ? "" : String(v); return /[",\n]/.test(v) ? `"${v.replace(/"/g, '""')}"` : v; }; + download(`${d.label}_qc_table.csv`, [vis.map((c) => esc(c[1])).join(","), ...r.map((s) => vis.map((c) => esc(c[2](s))).join(","))].join("\n") + "\n", "text/csv"); + } + function refresh() { + const r = rows(), n = r.length, per = st.per === 0 ? n || 1 : st.per; + st.page = Math.min(st.page, Math.max(0, Math.ceil(n / per) - 1)); + const page = r.slice(st.page * per, st.page * per + per); + const vis = cols.filter((c) => !st.hidden.has(c[0])); + const colMenu = $("div", { class: "menu", hidden: "" }, cols.map((c) => $("label", { style: "display:block" }, + (() => { const cb = $("input", { type: "checkbox", onchange: (e) => { e.target.checked ? st.hidden.delete(c[0]) : st.hidden.add(c[0]); refresh(); } }); cb.checked = !st.hidden.has(c[0]); return cb; })(), " ", c[1]))); + const head = $("div", { class: "pager" }, + $("input", { type: "text", placeholder: "Search…", value: st.search, oninput: (e) => { st.search = e.target.value; st.page = 0; refreshSoon(); } }), + $("span", { class: "dropdown" }, $("button", { class: "btn", onclick: () => { colMenu.hidden = !colMenu.hidden; } }, "Columns ▾"), colMenu), + $("button", { class: "btn", onclick: exportRows }, "Export ▾ CSV"), + $("button", { class: "btn", onclick: () => { st.rowsel.clear(); selection = null; applySelection(null); refresh(); } }, "Select none"), + select("Show", [[10, 10], [25, 25], [50, 50], [100, 100], [0, "All"]], st.per, (v) => { st.per = +v; st.page = 0; refresh(); })); + const thead = $("tr", {}, $("th", {}, ""), vis.map((c) => $("th", { onclick: () => { st.asc = st.sort === c[0] ? !st.asc : true; st.sort = c[0]; refresh(); } }, + c[1] + (st.sort === c[0] ? (st.asc ? " ▲" : " ▼") : "")))); + const foot = $("tr", {}, $("td", {}, ""), vis.map((c) => { + const numeric = d.samples.some((s) => typeof c[2](s) === "number"); + if (!numeric) return $("td", {}, $("input", { type: "text", placeholder: "filter", value: st.filt[c[0]] || "", onchange: (e) => { st.filt[c[0]] = e.target.value; refresh(); } })); + const f = st.filt[c[0]] || {}; + return $("td", {}, $("input", { type: "number", step: "any", placeholder: "min", value: f.min ?? "", onchange: (e) => { st.filt[c[0]] = { ...(st.filt[c[0]] || {}), min: e.target.value }; refresh(); } }), + $("input", { type: "number", step: "any", placeholder: "max", value: f.max ?? "", onchange: (e) => { st.filt[c[0]] = { ...(st.filt[c[0]] || {}), max: e.target.value }; refresh(); } })); + })); + const body = page.map((s) => { + const cb = $("input", { type: "checkbox", onchange: (e) => { e.target.checked ? st.rowsel.add(s.sample_id) : st.rowsel.delete(s.sample_id); + selection = st.rowsel.size ? new Set(st.rowsel) : null; applySelection(null); } }); + cb.checked = st.rowsel.has(s.sample_id); + return $("tr", { class: st.rowsel.has(s.sample_id) ? "sel" : "" }, $("td", {}, cb), vis.map((c) => { + const v = c[2](s); let cls = ""; + if (c[0] === "functional_status" || c[0] === "control_status") cls = v; + else if (c[0] === "control_flag_codes" && v) cls = "code"; + else if (key[c[0]] && typeof v === "number" && enabled(key[c[0]]) && fails(key[c[0]], v)) cls = c[0] === "autosomal_call_rate" || c[0] === "logr_dev" ? "fout" : "cout"; + return $("td", { class: cls }, fmt(v)); + })); + }); + const pager = $("div", { class: "pager" }, `Showing ${n ? st.page * per + 1 : 0} to ${Math.min(n, st.page * per + per)} of ${n} entries`, + $("button", { class: "btn", onclick: () => { st.page = Math.max(0, st.page - 1); refresh(); } }, "Prev"), + $("button", { class: "btn", onclick: () => { st.page = Math.min(Math.ceil(n / per) - 1, st.page + 1); refresh(); } }, "Next")); + wrap.replaceChildren(head, $("div", { class: "tablewrap" }, $("table", { class: "dt" }, $("thead", {}, thead), $("tbody", {}, body), $("tfoot", {}, foot))), pager); + } + let timer = null; + const refreshSoon = () => { clearTimeout(timer); timer = setTimeout(() => { refresh(); const i = wrap.querySelector("input[type=text]"); if (i) { i.focus(); i.setSelectionRange(i.value.length, i.value.length); } }, 250); }; + refresh(); + return { refresh }; + } + + // ------------------------------------------------------------------------------------------ view B: automated QC + const autoState = { colour: "", sort: "barcode" }; + function renderAuto() { + const root = document.getElementById("v-auto"); root.replaceChildren(); + linked.length = 0; + const d = ds(); + root.append($("div", { class: "controls" }, datasetSelect(renderAuto), + select("Colour points by", [["", "Functional status"], ...eligible(d).map((c) => [c, c])], autoState.colour, (v) => { autoState.colour = v; renderAuto(); }), + select("Sort by", sortOptions(d), autoState.sort, (v) => { autoState.sort = v; renderAuto(); }))); + const hist = $("div", { class: "grid2" }); root.append(hist); + const bars = (title, levels, field) => { + const n = d.samples.length, counts = levels.map((l) => d.samples.filter((s) => s[field] === l).length); + const c = $("div", { class: "card" }, $("h3", {}, title), $("div", { class: "plot small" })); hist.append(c); + Plotly.newPlot(c.lastChild, [{ type: "bar", orientation: "h", y: levels, x: counts, marker: { color: levels.map((l) => COL[l]) }, + text: counts.map((k) => `${k} (${n ? Math.round((100 * k) / n) : 0}%)`), textposition: "auto", hoverinfo: "y+x" }], + { margin: { l: 80, r: 20, t: 10, b: 30 } }, CONFIG); + }; + if (d.hasSummary) bars("Functional QC", ["FAIL", "PASS"], "functional_status"); + bars("Control QC", ["FLAGGED", "CLEAR"], "control_status"); + const tcard = $("div", { class: "card", style: "margin:12px 0" }, $("h3", {}, "Samples")); root.append(tcard); + table = makeTable(tcard, d); + const order = sorted(d.samples, autoState.sort), ids = order.map((s) => s.sample_id); + const grid = $("div", { class: "grid2" }); root.append(grid); + const metrics = [...(d.hasSummary ? D.functionalKeys.map((f) => ({ column: f.column, key: f.key, label: f.label, formula: f.formula })) : []), + ...d.metrics.map((m) => ({ column: m, key: MDEF[m] ? MDEF[m].key : m, label: m, formula: MDEF[m] ? MDEF[m].formula : "" }))]; + const groupsOf = () => { + if (!autoState.colour) return ["PASS", "FAIL", "NA"].map((l) => [l, COL[l], (s) => (s.functional_status || "NA") === l]); + const lv = [...new Set(d.samples.map((s) => s.meta[autoState.colour] || "(empty)"))].sort(); + const pal = PALETTES["Okabe-Ito"]; + return lv.map((l, i) => [l, pal[i % pal.length], (s) => (s.meta[autoState.colour] || "(empty)") === l]); + }; + for (const m of metrics) { + const card = $("div", { class: "card" }, $("h3", {}, m.label), $("div", { class: "help" }, m.formula), $("div", { class: "plot" })); + grid.append(card); + const traces = groupsOf().map(([name, color, pred]) => { + const sub = order.filter(pred); + return { type: scatterType(order.length), mode: "markers", name, x: sub.map((s) => s.sample_id), y: sub.map((s) => val(s, m.column)), + customdata: sub.map((s) => s.sample_id), text: sub.map((s) => hover(s, `${m.label}: ${fmt(val(s, m.column))}`)), hoverinfo: "text", marker: { color, size: 6 } }; + }).filter((t) => t.x.length); + const tl = thLine(m.key, false); + Plotly.newPlot(card.lastChild, traces, { margin: { l: 60, r: 10, t: 24, b: 30 }, dragmode: "select", hovermode: "closest", + xaxis: { type: "category", categoryorder: "array", categoryarray: ids, showticklabels: false, title: "Sample ID" }, + yaxis: { title: m.label }, showlegend: true, legend: { orientation: "h", y: -0.15 }, shapes: tl.shapes, annotations: tl.annotations }, CONFIG); + registerLinked(card.lastChild, order); + } + } + + // ------------------------------------------------------------------------------------------ view C: heatmaps (QC-052) + const heatState = { scale: "Continuous", group: "Chip", per: 2, sort: "barcode", page: 0, plate: "" }; + function renderHeat() { + const root = document.getElementById("v-heat"); root.replaceChildren(); + const d = ds(); + if (!d.hasSummary) { root.append($("div", { class: "controls" }, datasetSelect(renderHeat)), $("p", {}, "Heatmaps need the sample summary (functional metrics).")); return; } + const plates = [...new Set(d.samples.map((s) => s.meta.sample_plate).filter(Boolean))].sort(); + root.append($("div", { class: "controls" }, datasetSelect(renderHeat), + select("Colour scale", [["Continuous", "Continuous"], ["Binary", "Binary"]], heatState.scale, (v) => { heatState.scale = v; renderHeat(); }), + select("Group by", [["Chip", "Chip"], ["Plate", "Plate"]], heatState.group, (v) => { heatState.group = v; heatState.page = 0; renderHeat(); }), + heatState.group === "Chip" ? select("Max chips per page", [[1, 1], [2, 2], [4, 4], [8, 8]], heatState.per, (v) => { heatState.per = +v; heatState.page = 0; renderHeat(); }) : null, + select("Sort by", [["barcode", heatState.group === "Chip" ? "Barcode" : "Plate"], ["mean", "Call-rate mean"]], heatState.sort, (v) => { heatState.sort = v; renderHeat(); }), + heatState.group === "Plate" ? select("Plate", [["", "All"], ...plates.map((p) => [p, p])], heatState.plate, (v) => { heatState.plate = v; heatState.page = 0; renderHeat(); }) : null)); + let groups; + if (heatState.group === "Chip") { + const by = {}; d.samples.filter((s) => s.barcode).forEach((s) => (by[s.barcode] = by[s.barcode] || []).push(s)); + groups = Object.entries(by); + } else { + const by = {}; d.samples.forEach((s) => { const p = s.meta.sample_plate; if (p && (!heatState.plate || heatState.plate === p)) (by[p] = by[p] || []).push(s); }); + groups = Object.entries(by); + const noPlate = d.samples.filter((s) => !s.meta.sample_plate); + if (noPlate.length) root.append($("div", { class: "muted" }, `No plate: ${noPlate.map((s) => s.sample_id).join(", ")}`)); + } + const mean = (ss) => ss.reduce((a, s) => a + (s.autosomal_call_rate || 0), 0) / ss.length; + groups.sort((a, b) => (heatState.sort === "mean" ? mean(a[1]) - mean(b[1]) : a[0].localeCompare(b[0]))); + const per = heatState.group === "Chip" ? heatState.per : 1, pages = Math.max(1, Math.ceil(groups.length / per)); + heatState.page = Math.min(heatState.page, pages - 1); + root.append($("div", { class: "pager" }, $("button", { class: "btn", onclick: () => { heatState.page = Math.max(0, heatState.page - 1); renderHeat(); } }, "Prev"), + `${heatState.page + 1} / ${pages}`, $("button", { class: "btn", onclick: () => { heatState.page = Math.min(pages - 1, heatState.page + 1); renderHeat(); } }, "Next"))); + if (heatState.scale === "Binary") root.append($("div", { class: "muted" }, + $("span", { class: "legend-sw", style: "background:#E8871E" }), "call rate below cut-off", $("span", { class: "legend-sw", style: "background:#1B7837" }), "call rate at/above cut-off", + $("span", { class: "legend-sw", style: "background:#2166AC" }), "LRR SD below cut-off", $("span", { class: "legend-sw", style: "background:#F2E14E" }), "LRR SD above cut-off")); + for (const [name, ss] of groups.slice(heatState.page * per, heatState.page * per + per)) { + let nr, nc, pos; + if (heatState.group === "Chip") { + const f = d.chipFormat; nr = f ? f[0] : Math.max(...ss.map((s) => s.row)); nc = f ? f[1] : Math.max(...ss.map((s) => s.col)); + pos = (s) => [s.row - 1, s.col - 1]; + } else { + nr = 8; nc = 12; + pos = (s) => { const m = /^([A-Ha-h])0?(\d{1,2})$/.exec(s.meta.sample_well || ""); return m ? [m[1].toUpperCase().charCodeAt(0) - 65, +m[2] - 1] : null; }; + } + const rowLab = heatState.group === "Chip" ? [...Array(nr)].map((_, i) => `R${String(i + 1).padStart(2, "0")}`) : "ABCDEFGH".split(""); + const colLab = heatState.group === "Chip" ? [...Array(nc)].map((_, i) => `C${String(i + 1).padStart(2, "0")}`) : [...Array(12)].map((_, i) => String(i + 1)); + const card = $("div", { class: "card", style: "margin-bottom:12px" }, $("h3", {}, name), $("div", { class: "grid2" })); + root.append(card); + for (const [col, key, title, dec, better] of [["autosomal_call_rate", "callRate", "Autosomal Call Rate (higher is better)", 3, "high"], ["logr_dev", "logRDev", "LRR SD (lower is better)", 2, "low"]]) { + const z = [...Array(nr)].map(() => Array(nc).fill(null)), txt = [...Array(nr)].map(() => Array(nc).fill("missing")), hov = [...Array(nr)].map(() => Array(nc).fill("missing")); + for (const s of ss) { const p = pos(s); if (!p || p[0] >= nr || p[1] >= nc) continue; const v = s[col]; + z[p[0]][p[1]] = heatState.scale === "Binary" ? (v === null ? null : fails(key, v) ? 1 : 0) : v; + txt[p[0]][p[1]] = v === null ? "" : v.toFixed(dec); hov[p[0]][p[1]] = `${s.sample_id}
    ${s.sample_name}
    ${rowLab[p[0]]} ${colLab[p[1]]}
    ${fmt(v)}`; } + const cs = heatState.scale === "Binary" + ? (better === "high" ? [[0, "#1B7837"], [1, "#E8871E"]] : [[0, "#2166AC"], [1, "#F2E14E"]]) + : (better === "high" ? [[0, "#E5F5E0"], [1, "#00441B"]] : [[0, "#08306B"], [1, "#DEEBF7"]]); + const div = $("div", { class: "plot" }); card.lastChild.append($("div", {}, $("div", { class: "help" }, `${title} · ${enabled(key) ? "Cut-off = " + TH[key].value : "Cut-off disabled"}`), div)); + Plotly.newPlot(div, [{ type: "heatmap", z, x: colLab, y: rowLab, text: txt, texttemplate: "%{text}", hovertext: hov, hoverinfo: "text", + colorscale: cs, zmin: heatState.scale === "Binary" ? 0 : undefined, zmax: heatState.scale === "Binary" ? 1 : undefined, showscale: heatState.scale !== "Binary", xgap: 2, ygap: 2 }], + { margin: { l: 50, r: 10, t: 10, b: 30 }, yaxis: { autorange: "reversed" }, plot_bgcolor: "#d9d9d9" }, CONFIG); + } + } + } + + // ------------------------------------------------------------------------------------------ view D: trend (QC-054, QC-059, QC-060) + const trendState = { group: "dataset", sort: "name" }; + function chipStats(d) { // QC-053 + const by = {}; d.samples.filter((s) => s.barcode).forEach((s) => (by[s.barcode] = by[s.barcode] || []).push(s)); + const f = d.chipFormat, chips = Object.values(by); + const exp = (ss) => (f ? f[0] * f[1] : Math.max(...ss.map((s) => s.row)) * Math.max(...ss.map((s) => s.col))); + return { n: chips.length, complete: chips.filter((ss) => ss.length >= exp(ss)).length, missing: chips.reduce((a, ss) => a + Math.max(0, exp(ss) - ss.length), 0) }; + } + function renderTrend() { + const root = document.getElementById("v-trend"); root.replaceChildren(); + root.append($("div", { class: "controls" }, + select("Group by", [["dataset", "Dataset"], ["chip", "Chip"]], trendState.group, (v) => { trendState.group = v; renderTrend(); }), + select("Sort by", [["name", "Name"], ["median", "Median call rate"]], trendState.sort, (v) => { trendState.sort = v; renderTrend(); }))); + const hdr1 = $("tr", {}, $("th", { rowspan: 2 }, "Dataset"), $("th", { colspan: 4 }, "Scan"), $("th", { colspan: 4 }, "Autosomal Call Rate"), + $("th", { colspan: 3 }, "LRR SD"), $("th", { colspan: 2 }, "Assay control"), $("th", { colspan: 3 }, "Sex estimate")); + const hdr2 = $("tr", {}, ["date range", "# samples", "# complete / # chips", "# missing positions", "min", "median", "mean (SD)", "# < threshold", + "mean (SD)", "median", "# > threshold", "mean (SD)", "# < threshold", "# M", "# F", "# U"].map((t) => $("th", {}, t))); + const body = D.datasets.map((d) => { + const dates = d.samples.map((s) => (s.meta.imaging_date || "").slice(0, 10)).filter(Boolean).sort(); + const cr = boxStats(d.samples.map((s) => s.autosomal_call_rate)), ld = boxStats(d.samples.map((s) => s.logr_dev)), ac = boxStats(d.samples.map((s) => s.m.assay_control_norm_r)); + const cs = chipStats(d), sex = (x) => d.samples.filter((s) => (s.sex_est || "U") === x).length; + return $("tr", {}, [d.label, dates.length ? (dates[0] === dates[dates.length - 1] ? dates[0] : `${dates[0]} to ${dates[dates.length - 1]}`) : "", d.samples.length, `${cs.complete} / ${cs.n}`, cs.missing, + cr ? fmt(cr.min) : "", cr ? fmt(cr.med) : "", cr ? `${fmt(cr.mean)} (${fmt(cr.sd)})` : "", d.samples.filter((s) => s.autosomal_call_rate !== null && fails("callRate", s.autosomal_call_rate)).length, + ld ? `${fmt(ld.mean)} (${fmt(ld.sd)})` : "", ld ? fmt(ld.med) : "", d.samples.filter((s) => s.logr_dev !== null && fails("logRDev", s.logr_dev)).length, + ac ? `${fmt(ac.mean)} (${fmt(ac.sd)})` : "", ac ? d.samples.filter((s) => s.m.assay_control_norm_r !== undefined && s.m.assay_control_norm_r !== null && fails("assayControl", s.m.assay_control_norm_r)).length : "", + sex("M"), sex("F"), sex("U")].map((v) => $("td", {}, v))); + }); + root.append($("div", { class: "card tablewrap" }, $("table", { class: "dt summarytable" }, $("thead", {}, hdr1, hdr2), $("tbody", {}, body)))); + const groups = []; + for (const d of D.datasets) { + if (trendState.group === "dataset") groups.push([d.label, d.samples]); + else { const by = {}; d.samples.forEach((s) => (by[`${d.label}:${s.barcode}`] = by[`${d.label}:${s.barcode}`] || []).push(s)); groups.push(...Object.entries(by)); } + } + const medOf = (ss) => (boxStats(ss.map((s) => s.autosomal_call_rate)) || { med: 0 }).med; + groups.sort((a, b) => (trendState.sort === "median" ? medOf(a[1]) - medOf(b[1]) : a[0].localeCompare(b[0]))); + const grid = $("div", { class: "grid2", style: "margin-top:12px" }); root.append(grid); + const card = (title, help) => { const c = $("div", { class: "card" }, $("h3", {}, title), $("div", { class: "help" }, help || ""), $("div", { class: "plot" })); grid.append(c); return c.lastChild; }; + const h = Math.max(300, 26 * groups.length + 80); + const boxCard = (title, col, key, get) => { + const div = card(title, "Box: Q1–Q3 with median; whiskers reach the most extreme points within 1.5 × IQR; points beyond are outliers."); + div.style.height = `${h}px`; + Plotly.newPlot(div, boxTraces(groups.map(([n, ss]) => [n, ss.map(get)])), { margin: { l: 150, r: 20, t: 20, b: 40 }, xaxis: { title }, ...thLine(key, true) }, CONFIG); + }; + if (D.datasets.some((d) => d.hasSummary)) { + boxCard("Autosomal call rate", "autosomal_call_rate", "callRate", (s) => s.autosomal_call_rate); + boxCard("LRR SD", "logr_dev", "logRDev", (s) => s.logr_dev); + } + if (D.datasets.some((d) => d.metrics.includes("assay_control_norm_r"))) boxCard("Assay control (norm R)", "assay_control_norm_r", "assayControl", (s) => s.m.assay_control_norm_r); + const cnt = card("Sample count"); cnt.style.height = `${h}px`; + Plotly.newPlot(cnt, [{ type: "bar", orientation: "h", y: groups.map((g) => g[0]), x: groups.map((g) => g[1].length), marker: { color: "#3B75AF" }, text: groups.map((g) => g[1].length), textposition: "auto" }], + { margin: { l: 150, r: 20, t: 20, b: 40 } }, CONFIG); + const sx = card("Sex estimates"); sx.style.height = `${h}px`; + Plotly.newPlot(sx, [["M", "Male", "#3B75AF"], ["F", "Female", "#E8871E"], ["U", "Unknown", "#9A9A9A"]].map(([k, n, c]) => ({ + type: "bar", orientation: "h", name: n, y: groups.map((g) => g[0]), x: groups.map((g) => g[1].filter((s) => (s.sex_est || "U") === k).length), marker: { color: c } })), + { barmode: "group", margin: { l: 150, r: 20, t: 20, b: 40 }, legend: { orientation: "h" } }, CONFIG); + // QC-059 sex view: X vs Y median R when the summary carries them + const all = D.datasets.flatMap((d) => d.samples); + const sexDiv = card("Sex check", "X vs Y median R per sample, coloured by the estimate; the symbol shows the check against the reported sex."); + if (all.some((s) => typeof s.meta.x_median_r === "string" && s.meta.x_median_r !== "")) { + Plotly.newPlot(sexDiv, ["M", "F", "U"].map((k) => { const ss = all.filter((s) => (s.sex_est || "U") === k); + return { type: "scatter", mode: "markers", name: k, x: ss.map((s) => +s.meta.x_median_r), y: ss.map((s) => +s.meta.y_median_r), text: ss.map((s) => s.sample_id), + marker: { color: { M: "#3B75AF", F: "#E8871E", U: "#9A9A9A" }[k], symbol: ss.map((s) => (s.sex_check === "mismatch" ? "x" : "circle")) } }; }), + { margin: { l: 60, r: 10, t: 10, b: 40 }, xaxis: { title: "X median R" }, yaxis: { title: "Y median R" } }, CONFIG); + } else { + const counts = ["match", "mismatch", "unknown"].map((k) => all.filter((s) => (s.sex_check || "unknown") === k).length); + Plotly.newPlot(sexDiv, [{ type: "bar", x: ["match", "mismatch", "unknown"], y: counts, marker: { color: ["#3B75AF", "#E8871E", "#9A9A9A"] }, text: counts, textposition: "auto" }], + { margin: { l: 50, r: 10, t: 10, b: 40 } }, CONFIG); + sexDiv.previousSibling.textContent = "Sex estimate checked against the reported sex (X/Y median R are not in this summary)."; + } + // QC-060 control-metric PCA of the current dataset + const d = ds(), pcaDiv = card(`Control-metric PCA (${d.label})`, "log-transformed, standardised derived control metrics"); + const mets = d.metrics.filter((m) => d.samples.every((s) => typeof s.m[m] === "number" && s.m[m] > 0)); + if (mets.length < 3 || d.samples.length < 10) { pcaDiv.style.height = "60px"; pcaDiv.textContent = "insufficient data"; } + else { + const X = d.samples.map((s) => mets.map((m) => Math.log(s.m[m]))); + const mu = mets.map((_, j) => X.reduce((a, r) => a + r[j], 0) / X.length); + const sd = mets.map((_, j) => Math.sqrt(X.reduce((a, r) => a + (r[j] - mu[j]) ** 2, 0) / (X.length - 1)) || 1); + const Z = X.map((r) => r.map((v, j) => (v - mu[j]) / sd[j])); + const C = mets.map((_, i) => mets.map((_, j) => Z.reduce((a, r) => a + r[i] * r[j], 0) / (Z.length - 1))); + const pcs = []; + for (let k = 0; k < 2; k++) { let v = mets.map((_, i) => 1 + i * 0.01); + for (let it = 0; it < 200; it++) { let w = C.map((row) => row.reduce((a, c, j) => a + c * v[j], 0)); + for (const p of pcs) { const dot = w.reduce((a, x, i) => a + x * p[i], 0); w = w.map((x, i) => x - dot * p[i]); } + const n = Math.hypot(...w) || 1; v = w.map((x) => x / n); } + pcs.push(v); } + const proj = Z.map((r) => pcs.map((p) => r.reduce((a, x, i) => a + x * p[i], 0))); + const field = d.metadata.find((c) => c.toLowerCase() === "sample_group") || eligible(d)[0]; + const lv = field ? [...new Set(d.samples.map((s) => s.meta[field] || "(empty)"))].sort() : ["all"]; + Plotly.newPlot(pcaDiv, lv.map((l, i) => { const idx = d.samples.map((s, j) => j).filter((j) => !field || (d.samples[j].meta[field] || "(empty)") === l); + return { type: "scatter", mode: "markers", name: l, x: idx.map((j) => proj[j][0]), y: idx.map((j) => proj[j][1]), text: idx.map((j) => d.samples[j].sample_id), + marker: { color: PALETTES["Okabe-Ito"][i % 8] } }; }), { margin: { l: 50, r: 10, t: 10, b: 40 }, xaxis: { title: "PC1" }, yaxis: { title: "PC2" }, + showlegend: true, legend: { orientation: "h", y: -0.2, title: { text: field || "" } } }, CONFIG); + } + } + + // ------------------------------------------------------------------------------------------ view E: thresholds (QC-055) + const RANGE_OK = (k, v) => (k === "callRate" ? v >= 0 && v <= 1 : v >= 0); + function renderModal() { + const box = document.getElementById("modal-box"); box.replaceChildren(); + const inputs = {}; + const row = (k, label, desc, verb) => { + const t = TH[k] || { value: null, direction: "below" }; + const inp = $("input", { type: "number", step: "any", value: t.value ?? "", style: "width:90px" }); + const dir = $("select", {}, ["below", "above"].map((x) => { const o = $("option", { value: x }, `${verb} if ${x}`); if (x === t.direction) o.selected = true; return o; })); + const err = $("span", { class: "err" }); + inputs[k] = { inp, dir, err }; + return $("tr", {}, $("td", {}, inp), $("td", {}, dir), $("td", {}, label, err), $("td", { class: "muted" }, desc)); + }; + box.append($("h2", {}, "QC Thresholds"), $("p", { class: "muted" }, "Blank disables a check. Changes apply to this browser session only."), + $("h3", {}, "Sample QC metrics"), $("table", {}, $("tr", {}, ["Threshold", "Fail if", "Metric", "Description"].map((x) => $("th", {}, x))), + D.functionalKeys.map((f) => row(f.key, f.label, f.formula, "Fail"))), + $("h3", {}, "Control QC metrics"), $("table", {}, $("tr", {}, ["Threshold", "Flag if", "Metric", "Description"].map((x) => $("th", {}, x))), + D.metricDefs.filter((m) => D.datasets.some((d) => d.metrics.includes(m.column))).map((m) => row(m.key, m.column, m.formula, "Flag")))); + const msg = $("div", { class: "err" }); + const fileIn = $("input", { type: "file", accept: ".yaml,.yml", hidden: "", onchange: async (e) => { + const f = e.target.files[0]; if (!f) return; msg.textContent = ""; + try { const parsed = parseYaml(await f.text()); + for (const [k, v] of Object.entries(parsed)) { if (!inputs[k]) continue; inputs[k].inp.value = v.value ?? ""; if (v.direction) inputs[k].dir.value = v.direction; } + } catch (err) { msg.textContent = String(err.message || err); } } }); + box.append($("div", { class: "pager", style: "margin-top:12px" }, + $("button", { class: "btn", onclick: () => { + let ok = true; + for (const [k, { inp, dir, err }] of Object.entries(inputs)) { err.textContent = ""; const raw = inp.value.trim(); + if (raw === "") continue; const v = Number(raw); + if (!Number.isFinite(v)) { err.textContent = " value must be numeric or blank"; ok = false; } + else if (!RANGE_OK(k, v)) { err.textContent = k === "callRate" ? " must be within [0, 1]" : " must be ≥ 0"; ok = false; } } + if (!ok) return; + for (const [k, { inp, dir }] of Object.entries(inputs)) TH[k] = { ...(TH[k] || {}), value: inp.value.trim() === "" ? null : Number(inp.value), direction: dir.value, source: "session" }; + recompute(); closeModal(); render(); } }, "Apply"), + $("button", { class: "btn", onclick: () => { + const lines = ["qc:", " report:"]; + for (const [k, { inp, dir }] of Object.entries(inputs)) { + const v = inp.value.trim() === "" ? "null" : inp.value.trim(), def = k === "logRDev" ? "above" : "below"; + lines.push(dir.value === def ? ` ${k}: ${v}` : ` ${k}: {threshold: ${v}, direction: ${dir.value}}`); } + if (TH.bgOffset && TH.bgOffset.value !== null) lines.push(` bgOffset: ${TH.bgOffset.value}`); + download("qc_thresholds.yaml", lines.join("\n") + "\n", "text/yaml"); } }, "Download YAML"), + $("button", { class: "btn", onclick: () => fileIn.click() }, "Upload YAML"), fileIn, + $("button", { class: "btn", onclick: closeModal }, "Close"), msg)); + } + function parseYaml(text) { // the QC-020 subset: qc/QC → report/Report → key: value | {threshold, direction} + const lines = text.split(/\r?\n/).map((l) => l.replace(/#.*$/, "")).filter((l) => l.trim()); + const ind = (l) => l.length - l.trimStart().length; + const top = lines.findIndex((l) => /^(qc|QC)\s*:\s*$/.test(l.trim()) && ind(l) === 0); + if (top < 0) throw new Error("missing top-level 'qc:' (or 'QC:') mapping"); + const rep = lines.findIndex((l, i) => i > top && /^(report|Report)\s*:\s*$/.test(l.trim())); + if (rep < 0) throw new Error("missing 'report:' (or 'Report:') section"); + const known = new Set(Object.keys(D.thresholds)), alias = { tgaControl: "assayControl", call_rate_min: "callRate", logr_dev_max: "logRDev" }; + const out = {}; + for (let i = rep + 1; i < lines.length && ind(lines[i]) > ind(lines[rep]); i++) { + const m = /^\s*([A-Za-z_][\w]*)\s*:\s*(.*)$/.exec(lines[i]); if (!m) throw new Error(`line not understood: ${lines[i].trim()}`); + const k = alias[m[1]] || m[1]; if (!known.has(k)) throw new Error(`unknown QC threshold key '${m[1]}'`); + let raw = m[2].trim(), direction; + const mm = /^\{\s*threshold\s*:\s*([^,}]*)\s*(?:,\s*direction\s*:\s*(\w+)\s*)?\}$/.exec(raw); + if (mm) { raw = mm[1].trim(); direction = mm[2]; if (direction && !["below", "above"].includes(direction)) throw new Error(`${k}: direction must be 'below' or 'above'`); } + if (["", "null", "~", "none"].includes(raw.toLowerCase())) { out[k] = { value: null, direction }; continue; } + const v = Number(raw); if (!Number.isFinite(v)) throw new Error(`${k}: value must be numeric or null`); + if (!RANGE_OK(k, v)) throw new Error(`${k}: value out of range`); + out[k] = { value: v, direction }; + } + return out; + } + const closeModal = () => { document.getElementById("modal").hidden = true; }; + document.getElementById("modal").addEventListener("click", (e) => { if (e.target.id === "modal") closeModal(); }); + + // ------------------------------------------------------------------------------------------ navigation + let view = ["dash", "auto", "heat", "trend"].includes(location.hash.slice(1)) ? location.hash.slice(1) : "dash"; + const renderers = { dash: renderDash, auto: renderAuto, heat: renderHeat, trend: renderTrend }; + function render() { + table = null; + for (const v of Object.keys(renderers)) document.getElementById(`v-${v}`).hidden = v !== view; + renderers[view](); + } + document.querySelectorAll("#nav button").forEach((b) => b.addEventListener("click", () => { + if (b.dataset.view === "thr") { renderModal(); document.getElementById("modal").hidden = false; return; } + document.querySelectorAll("#nav button").forEach((x) => x.classList.toggle("active", x === b)); + view = b.dataset.view; render(); + })); + document.querySelectorAll("#nav button").forEach((x) => x.classList.toggle("active", x.dataset.view === view)); + render(); + if (location.hash === "#selftest") { // exercised by the test suite through a headless browser + const res = {}; + try { parseYaml("qc:\n report:\n bogus: 1\n"); res.yamlUnknown = "accepted"; } catch (e) { res.yamlUnknown = "rejected"; } + const y = parseYaml("QC:\n Report:\n callRate: 0.9995\n stringency: {threshold: 2, direction: above}\n npRed: null\n"); + res.yaml = y.callRate.value === 0.9995 && y.stringency.direction === "above" && y.npRed.value === null; + renderModal(); + document.querySelector("#modal-box input").value = "0.9985"; + [...document.querySelectorAll("#modal-box button")].find((x) => x.textContent === "Apply").click(); + res.failAfterApply = D.datasets[0].samples.filter((s) => s.functional_status === "FAIL").length; + view = "auto"; render(); + selection = new Set(D.datasets[0].samples.slice(0, 3).map((s) => s.sample_id)); applySelection(null); + res.rowsWithSelection = document.querySelectorAll("#v-auto table.dt tbody tr").length; + res.linkedPlots = linked.length; + document.body.dataset.selftest = JSON.stringify(res); + } + document.body.dataset.ready = "1"; +})().catch((e) => { + document.body.prepend(Object.assign(document.createElement("pre"), { className: "err", textContent: "Report error: " + (e && e.stack || e) })); + throw e; +}); diff --git a/src/simurg/qc/report/template.html b/src/simurg/qc/report/template.html new file mode 100644 index 0000000..32e59e7 --- /dev/null +++ b/src/simurg/qc/report/template.html @@ -0,0 +1,37 @@ + + + + + +__TITLE__ + + + +
    +
    __TITLE__
    +
    loading…
    + + +
    +
    + +
    +
    + + + +
    +
    + + + + + + diff --git a/src/simurg/qc/report/vendor/plotly.LICENSE b/src/simurg/qc/report/vendor/plotly.LICENSE new file mode 100644 index 0000000..1196ce9 --- /dev/null +++ b/src/simurg/qc/report/vendor/plotly.LICENSE @@ -0,0 +1,21 @@ +The MIT License (MIT) + +Copyright (c) 2021 Plotly, Inc + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. diff --git a/src/simurg/qc/report/vendor/plotly.min.js.LICENSE.txt b/src/simurg/qc/report/vendor/plotly.min.js.LICENSE.txt new file mode 100644 index 0000000..970ecab --- /dev/null +++ b/src/simurg/qc/report/vendor/plotly.min.js.LICENSE.txt @@ -0,0 +1 @@ +Couldn't find the requested file /plotly.min.js.LICENSE.txt in plotly.js-dist-min. \ No newline at end of file diff --git a/src/simurg/qc/report/vendor/plotly.min.js.gz b/src/simurg/qc/report/vendor/plotly.min.js.gz new file mode 100644 index 0000000..8132396 Binary files /dev/null and b/src/simurg/qc/report/vendor/plotly.min.js.gz differ diff --git a/tests/unit/test_qc_report.py b/tests/unit/test_qc_report.py new file mode 100644 index 0000000..02f21d1 --- /dev/null +++ b/tests/unit/test_qc_report.py @@ -0,0 +1,167 @@ +"""QC report (spec QC-040 … QC-063): dataset discovery, labels, types, table, HTML determinism and error codes.""" + +from __future__ import annotations + +import datetime as dt +import json +import shutil +import subprocess +from pathlib import Path + +import pytest +from typer.testing import CliRunner + +from simurg.cli.main import app +from simurg.core.errors import SimurgError +from simurg.core.log import Reporter +from simurg.pipeline.qc_report import QcReportOptions, assign_labels, resolve_datasets, run_qc_report + +RAW_HDR = "sample_id\tcategory\tbead_address\tcontrol_name\tcolor_label\tgreen_intensity\tred_intensity\n" +PROBES = [("Staining", "1", "Biotin (High)", "Green"), ("Staining", "2", "Biotin (Bgnd)", "Blue"), + ("Non-Polymorphic", "3", "NP (A)", "Red"), ("NEGATIVE", "4", "neg", "Red")] # fmt: skip + + +def make_dataset(folder: Path, n: int = 4, array_type: str = "genotyping", summary: bool = True) -> Path: + folder.mkdir(parents=True, exist_ok=True) + ids = [f"2000000000{k}_R0{i + 1}C01" for k in (1,) for i in range(n)] + raw = RAW_HDR + "".join( + f"{s}\t{c}\t{a}\t{nm}\t{col}\t{1000 * (j + 1) + i}\t{500 + i}\n" + for i, s in enumerate(ids) + for j, (c, a, nm, col) in enumerate(PROBES) + ) + (folder / "controls_raw.tsv").write_text(raw) + qc = "sample_id\tarray_type\timaging_date\tsample_plate\tsample_well\tscanner_id\tscanner_version\tbg_offset\tstaining_green\tnp_red\tna_metrics\n" + qc += "".join( + f"{s}\t{array_type}\t2026-01-0{i + 1}T10:00:00\tP1\tA0{i + 1}\tN1\t1.0\t3000\t{6 - i}\t{10 + i}\t\n" + for i, s in enumerate(ids) + ) + (folder / "controls_qc.tsv").write_text(qc) + if summary: + cols = ["sample_id", "sample_name", "sample_folder", "status", "fail_reason", "autosomal_call_rate", + "call_rate", "logr_dev", "sex_est", "sex_reported", "sex_check", "het_rate", "contamination", + "functional_status", "control_status", "control_flag_codes", "flags", "Sample_Group"] # fmt: skip + lines = ["\t".join(cols)] + for i, s in enumerate(ids): + vals = {"sample_id": s, "sample_name": f"S{i}", "sample_folder": "x", "status": "OK", "fail_reason": "", + "autosomal_call_rate": f"{0.99 - 0.01 * i}", "call_rate": "0.99", "logr_dev": f"{0.1 + 0.05 * i}", + "sex_est": "MF"[i % 2], "sex_reported": "U", "sex_check": "unknown", "het_rate": "0.17", + "contamination": "0.001", "functional_status": "", "control_status": "", + "control_flag_codes": "", "flags": "", "Sample_Group": "g" + str(i % 2)} # fmt: skip + lines.append("\t".join(vals[c] for c in cols)) + (folder / "sample_summary.tsv").write_text("\n".join(lines) + "\n") + return folder + + +def run(tmp_path: Path, data: list[str], **kw: object) -> Path: + out = tmp_path / "out" + opts = QcReportOptions(data=data, out=out, report_name="r", now=dt.datetime(2026, 1, 2, tzinfo=dt.UTC), **kw) # type: ignore[arg-type] + assert run_qc_report(opts, Reporter(command="t", quiet=True)) == 0 + return out + + +def test_discovery_forms(tmp_path: Path) -> None: # QC-041 + a = make_dataset(tmp_path / "p" / "b") + b = make_dataset(tmp_path / "p" / "a") + assert resolve_datasets([str(tmp_path / "p")]) == [b, a] # sorted by path + assert resolve_datasets([f"{a},{b}"]) == [a, b] + with pytest.raises(SimurgError) as e: + resolve_datasets([str(tmp_path / "nope")]) + assert e.value.code == "E_QC_NO_DATASETS" + (tmp_path / "empty").mkdir() + with pytest.raises(SimurgError) as e: + resolve_datasets([str(tmp_path / "empty")]) + assert e.value.code == "E_QC_NO_DATASETS" + (b / "controls_raw.tsv").unlink() + with pytest.raises(SimurgError) as e: + resolve_datasets([f"{a},{b}"]) + assert e.value.code == "E_QC_DATASET_INCOMPLETE" + + +def test_labels(tmp_path: Path) -> None: # QC-042 + rep = Reporter(command="t", quiet=True) + f = [tmp_path / "x" / "run", tmp_path / "y" / "run"] + assert assign_labels(f, None, rep) == ["run", "run_2"] + assert assign_labels(f, ["My Label<1>", "b"], rep) == ["My_Label_1_", "b"] + assert rep.warning_counts().get("W_QC_LABEL_SANITIZED") == 1 + for labels, code in ((["a"], "E_QC_LABEL_COUNT"), (["a b", "a_b"], "E_QC_LABEL_DUPLICATE")): + with pytest.raises(SimurgError) as e: + assign_labels(f, labels, rep) + assert e.value.code == code + + +def test_report_outputs_and_determinism(tmp_path: Path) -> None: # QC-044, QC-045, QC-047 + d = make_dataset(tmp_path / "ds") + out = run(tmp_path, [str(d)], output_format="both") + html = (out / "r.html").read_text() + assert "SimurgArray QC Report" in html and "https://" not in html.split("plotly.js v")[0] + first = html + out2 = run(tmp_path, [str(d)], output_format="both") + assert (out2 / "r.html").read_text() == first # byte-deterministic with a fixed timestamp + rows = (out / "r_table.csv").read_text().splitlines() + head = rows[0].split(",") + assert head[:7] == ["dataset", "sample_id", "sample_name", "functional_status", "control_status", + "control_flag_codes", "na_metrics"] # fmt: skip + assert "neg_red" not in head and "Biotin (High)_green" in head # QC-010 exclusion + assert head[-1] == "sample_well" or "Sample_Group" in head + from openpyxl import load_workbook + + wb = load_workbook(out / "r_table.xlsx") + assert wb.sheetnames == ["Samples", "Thresholds", "Datasets", "About"] + ws = wb["Samples"] + status = [c.value for c in ws["D"]][1:] + assert status == ["PASS", "PASS", "FAIL", "FAIL"] # call rate 0.97 < 0.98; LRR SD 0.25 > 0.20 + + +def test_privacy_and_errors(tmp_path: Path) -> None: # QC-043, QC-046, QC-063 + d = make_dataset(tmp_path / "ds") + out = run(tmp_path, [str(d)], output_format="csv", no_metadata=True) + head = (out / "r_table.csv").read_text().splitlines()[0] + assert "Sample_Group" not in head and "sample_plate" not in head + runner = CliRunner() + res = runner.invoke(app, ["qc", "report", "--data", str(d), "--out", str(tmp_path / "o"), "--output-format", "pdf"]) + assert res.exit_code == 2 and "E_QC_FORMAT" in res.output + m = make_dataset(tmp_path / "meth", array_type="methylation") + res = runner.invoke(app, ["qc", "report", "--data", f"{d},{m}", "--out", str(tmp_path / "o")]) + assert res.exit_code == 2 and "E_QC_MIXED_TYPES" in res.output + bad = tmp_path / "bad.yaml" + bad.write_text("qc: {report: {callRate: 2}}") + res = runner.invoke(app, ["qc", "report", "--data", str(d), "--out", str(tmp_path / "o"), "--config", str(bad)]) + assert res.exit_code == 2 and "E_QC_THRESHOLD_RANGE" in res.output + + +def test_unmatched_samples_warn(tmp_path: Path) -> None: # QC-040 + d = make_dataset(tmp_path / "ds") + lines = (d / "sample_summary.tsv").read_text().splitlines() + (d / "sample_summary.tsv").write_text("\n".join(lines[:-1]) + "\n") + rep = Reporter(command="t", quiet=True) + run_qc_report(QcReportOptions(data=[str(d)], out=tmp_path / "o", output_format="csv"), rep) + assert rep.warning_counts().get("W_QC_UNMATCHED_SAMPLES") == 1 + + +@pytest.mark.slow +def test_report_selftest_in_browser(tmp_path: Path) -> None: # QC-051, QC-055 (client-side behaviour) + browser = shutil.which("chromium") or shutil.which("chromium-browser") or shutil.which("google-chrome") + if not browser: + pytest.skip("no Chromium for the headless report test") + d = make_dataset(tmp_path / "ds", n=6) + out = run(tmp_path, [str(d)], output_format="csv") + dom = subprocess.run( + [ + browser, + "--headless=new", + "--no-sandbox", + "--disable-gpu", + "--virtual-time-budget=20000", + "--dump-dom", + f"file://{out / 'r.html'}#selftest", + ], + capture_output=True, + text=True, + timeout=180, + ).stdout + marker = 'data-selftest="' + assert marker in dom, dom[-2000:] + res = json.loads(dom.split(marker, 1)[1].split('"', 1)[0].replace(""", '"')) + assert res["yamlUnknown"] == "rejected" and res["yaml"] is True + assert res["failAfterApply"] == 6 # callRate 0.9985 fails every synthetic sample (max 0.99) + assert res["rowsWithSelection"] == 3 and res["linkedPlots"] > 0