P5a: QC — control metrics, thresholds, status, contamination, relatedness, locus summary

- qc/controls: canonical control names (QC-002), bkg/offset (QC-003), 13 genotyping metrics (QC-004/005),
  missing/zero-divisor handling (QC-007/008); controls_raw.tsv / controls_qc.tsv (QC-009/011) incl. compat CSVs
- qc/thresholds + profiles: YAML schema, validation codes, precedence, shipped profiles, auto (QC-020..025)
- qc/status: functional PASS/FAIL/NA, control CLEAR/FLAGGED, flags (QC-030..033)
- qc/contamination: BAF regression; bit-compatible with gtc2vcf BAFregress (36/36, K-017), exact mode too
- qc/relatedness: KING-robust, duplicates, Mendelian errors, pedigree from sheet (IN-018)
- qc/locus_summary: OUT-072, bit-identical to BeadArrayFiles except ChiTest (K-018)
- qc/suggest: suggest-thresholds (QC-026)
- pipeline/qc: disk-backed cohort store; CLI qc controls|samples|contamination|relatedness|loci|suggest-thresholds
- genotype: control tables (OUT-070), het_rate, QC columns, cohort QC
- tests: unit + oracle (BAFregress plugin, BeadArrayFiles locus summary)

README.md

@@ -12,7 +12,11 @@ Personal research toolkit for Illumina Infinium genotyping BeadChip data
 Phases **P0 (setup)**, **P1 (file parsers)**, **P2 (canonical model, raw X/Y)**, **P3 (normalization)** and **P4 (genotype
 calling, GTC output)** are complete — see `docs/validation/`. Each is validated on 2141 vendor GTCs of the 1000 Genomes Omni2.5 set
 (normalization bit-identical; genotypes/scores identical in 2138/2141 samples).
-Available commands: `simurg version`, `simurg inspect <files…> [--json]`, `simurg ingest`, `simurg normalize`, `simurg genotype`.
+Phase **P5 (QC)** is in progress: control-probe metrics, thresholds/profiles, two-axis sample status, BAF-regression
+contamination, KING-robust relatedness, Mendelian errors and the locus summary are implemented and match their open-source
+references on the GSA demo set (`docs/validation/P5.md`); the cohort-scale validation and the HTML report are pending.
+Available commands: `simurg version`, `simurg inspect <files…> [--json]`, `simurg ingest`, `simurg normalize`, `simurg genotype`,
+`simurg qc controls|samples|contamination|relatedness|loci|suggest-thresholds`.
 
 ## Quick start
 
@@ -22,6 +26,7 @@ uv run pytest tests/unit                  # no downloads needed
 bash scripts/bootstrap_refs.sh            # reference open-source repos → refs/ (oracles)
 bash scripts/fetch_data.sh p0             # ~1.4 GB of public test data → data/
 bash scripts/make_oracle_gtcs.sh          # builds gtc2vcf, writes oracle GTCs
+bash scripts/make_oracle_qc.sh            # locus-summary oracle (BeadArrayFiles, separate venv)
 uv run pytest tests/oracle                # real-data acceptance tests (add -m "not slow" to skip the long ones)
 uv run simurg inspect data/p0/gsa_v4_pgx/idats_Demo_36/idats_Demo_36/207895980019/207895980019_R03C01_Grn.idat
 ```

THIRD_PARTY_NOTICES.md

@@ -5,8 +5,9 @@ used only as external subprocess oracles and are not copied or linked.
 
 | Source | Licence | How it is used |
 |---|---|---|
-| [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`; test oracle |
+| [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`); 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 |
+| [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 |
 | [illuminaio](https://github.com/HenrikBengtsson/illuminaio) | GPL-2 | Documentation of IDAT field codes only; oracle via subprocess |
@@ -67,3 +68,31 @@ 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.
 ```
+
+## MIT (klib `ksort.h` in htslib, licence header verbatim)
+
+```
+The MIT License
+
+Copyright (c) 2008, 2012-2013, 2017-2019 Genome Research Ltd (GRL).
+
+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.
+```

docs/known_unknowns.md

@@ -4,6 +4,9 @@ Plan §20 lists the open questions; this file records what the implementation ha
 
 | ID | Date | Topic | Finding | Status |
 |---|---|---|---|---|
+| 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 |
 | K-016 | 2026-09-28 | Sex rule of AutoCall 1.6.2.2 | Rule 1 (X heterozygosity) matches all 2141 vendor GTCs; rule 2 (X/Y intensity) differs in 12. | resolved |
 | K-015 | 2026-09-28 | Legacy GTC v3 statistics | AutoCall 1.6.2.2 writes GC10 = GC50 = 0 and #calls = #no-calls = 0 in GTC v3 (it does fill LogR dev and #intensity-only). gtc2vcf's AutoConvert 1.6.3.1 preset fills them, so this is specific to 1.6.2.x. Emulated by `--preset autocall-1.6.2`; the sample summary still reports the real statistics. | resolved; cloud re-run pending |
 | K-014 | 2026-09-28 | NaN sign | NaNs produced by libm (e.g. log2 of a negative extrapolated R) take the platform's default sign: + on aarch64, − on x86-64. SimurgArray always writes 0xFFC00000 (the reference's `-NAN`, x86 behaviour). The aarch64-built oracle differs in 1 locus of 9 GSA samples; oracle tests compare with NaNs canonicalised. | resolved |

docs/messages.md

@@ -25,6 +25,21 @@ Generated from `src/simurg/core/errors.py` — do not edit by hand.
 | `E_NO_INPUT` | error | Neither an input folder nor a sample sheet was given | CLI-006 |
 | `E_NO_INPUT_FOLDER` | error | Sample sheet has no RootFolder and no folder was given | CLI-006 |
 | `E_OUTPUT_EXISTS` | error | Outputs already exist; use --force or --resume | CLI-010 |
+| `E_QC_CONFIG_EXT` | error | QC config file must be .yaml or .yml | QC-021 |
+| `E_QC_CONFIG_MISSING` | error | QC config file does not exist | QC-021 |
+| `E_QC_CONFIG_PARSE` | error | QC config file is not valid YAML | QC-021 |
+| `E_QC_CONFIG_STRUCTURE` | error | QC config lacks the qc/report (QC/Report) wrapper | QC-021 |
+| `E_QC_DATASET_INCOMPLETE` | error | Dataset folder lacks a required file | QC-041 |
+| `E_QC_FORMAT` | error | Invalid --output-format value | QC-063 |
+| `E_QC_LABEL_COUNT` | error | Number of labels does not match number of datasets | QC-042 |
+| `E_QC_LABEL_DUPLICATE` | error | Dataset labels are not unique | QC-042 |
+| `E_QC_METHYLATION` | error | Methylation control metrics are not available yet (v2) | QC-006 |
+| `E_QC_MIXED_TYPES` | error | Genotyping and methylation datasets cannot be mixed | QC-043 |
+| `E_QC_NO_DATASETS` | error | No QC dataset folder found | QC-041 |
+| `E_QC_PROFILE` | error | Unknown QC profile | QC-024 |
+| `E_QC_THRESHOLD_KEY` | error | Unknown QC threshold key | QC-021 |
+| `E_QC_THRESHOLD_RANGE` | error | QC threshold value out of range | QC-021 |
+| `E_QC_THRESHOLD_VALUE` | error | QC threshold value is not numeric or null | QC-021 |
 | `E_SHEET_BARCODE_FORMAT` | error | Invalid SentrixBarcode_A / SentrixPosition_A value | IN-007 |
 | `E_SHEET_DELIMITER` | error | Sample sheet is not comma-delimited CSV | IN-001 |
 | `E_SHEET_DUP_COLUMN` | error | Duplicate column (column names are case-insensitive) | IN-003 |
@@ -36,6 +51,9 @@ Generated from `src/simurg/core/errors.py` — do not edit by hand.
 | `E_SHEET_RESERVED` | error | Sample sheet column collides with an output column | IN-005 |
 | `E_SHEET_TOO_LONG` | error | Sample sheet name or value longer than 500 characters | IN-006 |
 | `I_UNLISTED_IDATS` | info | IDAT pairs not listed in the sample sheet are ignored | IN-016 |
+| `W_CTRL_DEFS_DIFFER` | warning | BPM and CSV control definitions differ; the BPM wins | QC-001 |
+| `W_CTRL_MISSING` | warning | Controls needed by a metric are absent; metric not evaluable | QC-007 |
+| `W_CTRL_ZERO_DIVISOR` | warning | Control metric has a zero denominator (NaN) | QC-008 |
 | `W_DUP_IDAT_RESOLVED` | warning | Duplicate IDAT pair resolved by --on-duplicate-idat | IN-022 |
 | `W_EGT_NAME` | warning | Cluster file name is not portable | IN-045 |
 | `W_GTC_FOREIGN` | warning | GTC was produced by another tool | IN-031 |
@@ -45,6 +63,11 @@ Generated from `src/simurg/core/errors.py` — do not edit by hand.
 | `W_IDAT_NO_SCAN_TIME` | warning | IDAT has no scan entry in its run info | IN-024 |
 | `W_IDAT_UNPAIRED` | warning | IDAT channel file without its partner | IN-021 |
 | `W_MANIFEST_NO_REFSTRAND` | warning | Manifest lacks RefStrand annotation | IN-040 |
+| `W_QC_LABEL_SANITIZED` | warning | Dataset label was sanitized | QC-042 |
+| `W_QC_NO_BAF` | warning | GTC has no B-allele frequencies; contamination not estimated | QC-064 |
+| `W_QC_PEDIGREE` | warning | Parent1/Parent2/Replicate value does not name a sample in the set | IN-018 |
+| `W_QC_TOO_FEW` | warning | Too few samples or sites for a cohort QC statistic | QC-026 |
+| `W_QC_UNMATCHED_SAMPLES` | warning | Sample IDs present in one QC file but not the other | QC-040 |
 | `W_SHEET_BARCODE_REPAIRED` | warning | Barcode written as a decimal was repaired | IN-007 |
 | `W_SHEET_EMPTY_COLUMN` | warning | Columns with empty headers are ignored | IN-004 |
 | `W_SHEET_ROOTFOLDER_OVERRIDDEN` | warning | CLI folder overrides the sheet RootFolder | IN-008 |

docs/validation/P5.md

@@ -0,0 +1,22 @@
+# Phase P5 validation — QC
+
+Status: **part 1 (local, GSA v4 PGx demo, 36 samples) passed**; part 2 (1000 Genomes Omni2.5, 2141 samples, cloud) pending.
+
+## Part 1 — local, against open-source references (`tests/oracle/test_p5_qc.py`)
+
+| Check | Reference | Result |
+|---|---|---|
+| BAF-regression contamination, 36 samples, all 683 054 loci as sites | gtc2vcf `bcftools +BAFregress` (MIT), built by `scripts/make_oracle_gtcs.sh` | 36/36 identical at the plugin's printed precision (`%.4f`), homozygous-site counts identical (K-017: float threshold 0.1f and klib quickselect median reproduced) |
+| Locus summary, 683 054 loci | BeadArrayFiles `examples/locus_summary.py` (BSD-2), `scripts/make_oracle_qc.sh` | every column bit-identical except ChiTest_P100 (≤ 3.6e-15 relative, K-018) |
+| Control metrics | QC-004 formulas | unit tests with hand-computed values; GSA demo values plausible (staining, extension, NP > 5; hyb > 1) |
+| Thresholds / profiles / status / flags | QC-020 … QC-033 | unit tests for every QC-021 error code, precedence, aliases, profiles, both status axes |
+| KING-robust kinship, Mendelian errors | Manichaikul 2010 | simulated trio: duplicate φ = 0.50, parent–child φ ≈ 0.25, unrelated ≈ 0, 0 Mendelian errors in the true trio |
+| End to end | `simurg qc controls`, `simurg genotype` with QC | control tables identical between the two paths; QC columns filled |
+
+GSA demo observations: all 36 samples PASS (autosomal call rate ≥ 0.997, LogR dev 0.13–0.17), control status CLEAR;
+contamination 0.0012–0.0040; replicate pairs (φ ≈ 0.50) and first-degree relatives (φ ≈ 0.25) are present (K-019).
+
+## Part 2 — cohort scale (pending)
+
+1000 Genomes Omni2.5, 2141 samples: relatedness and Mendelian errors against the 1000G pedigree (trios), sex against the
+pedigree, contamination and control-metric distributions, runtime/memory of the disk-backed cohort store.

scripts/make_oracle_gtcs.sh

@@ -9,14 +9,14 @@ cd "$(dirname "$0")/.."
 V=1.22
 T=refs/tools
 B=$T/bcftools-$V
-if [[ ! -x $B/bcftools || ! -f $B/plugins/idat2gtc.so ]]; then
+if [[ ! -x $B/bcftools || ! -f $B/plugins/idat2gtc.so || ! -f $B/plugins/BAFregress.so ]]; then
   [[ -d refs/gtc2vcf ]] || bash scripts/bootstrap_refs.sh
   mkdir -p $T
   [[ -f $T/bcftools.tar.bz2 ]] || curl -sL -o $T/bcftools.tar.bz2 \
     "https://github.com/samtools/bcftools/releases/download/$V/bcftools-$V.tar.bz2"
   tar -xjf $T/bcftools.tar.bz2 -C $T
-  rm -f $B/plugins/{idat2gtc,gtc2vcf,affy2vcf}.{c,h}
+  rm -f $B/plugins/{idat2gtc,gtc2vcf,affy2vcf,BAFregress}.{c,h}
-  cp refs/gtc2vcf/{idat2gtc,gtc2vcf,affy2vcf}.c refs/gtc2vcf/*.h $B/plugins/
+  cp refs/gtc2vcf/{idat2gtc,gtc2vcf,affy2vcf,BAFregress}.c refs/gtc2vcf/*.h $B/plugins/
   (cd $B/htslib-$V && ./configure --disable-libcurl --disable-bz2 --disable-lzma --disable-ref-cache >/dev/null)
   (cd $B && ./configure --without-gsl >/dev/null && make -j"$(nproc)" CFLAGS="-g -O2 -ffp-contract=off" >/dev/null)
 fi

scripts/make_oracle_qc.sh

@@ -0,0 +1,19 @@
+#!/usr/bin/env bash
+# Oracle outputs for phase P5 (QC) on the P0 GSA demo data:
+#  * BeadArrayFiles (BSD-2) examples/locus_summary.py over the 36 oracle GTCs, run in a separate venv with
+#    scipy/pandas (oracle-only dependencies; never needed by SimurgArray itself);
+#  * the gtc2vcf BAFregress plugin (MIT) is built by scripts/make_oracle_gtcs.sh and run by the tests directly.
+set -euo pipefail
+cd "$(dirname "$0")/.."
+V=refs/tools/oracle-venv
+[[ -x $V/bin/python ]] || { uv venv -q $V && uv pip install -q --python $V/bin/python numpy scipy pandas; }
+SHIM=refs/tools/shim
+mkdir -p $SHIM && ln -sfn "$PWD/refs/BeadArrayFiles/module" $SHIM/IlluminaBeadArrayFiles
+D=data/p0/gsa_v4_pgx
+OUT=data/p0/oracle/locus_summary
+mkdir -p $OUT
+PYTHONPATH=$SHIM $V/bin/python refs/BeadArrayFiles/examples/locus_summary.py data/p0/oracle/gsa_gt3_nofma \
+  $D/GSA-PGx-48v4-0_20079540_E2_bpm/GSA-PGx-48v4-0_20079540_E2.bpm \
+  $D/GSA-PGx-48v4-0_20079540_E2_ClusterFile/GSA-PGx-48v4-0_20079540_E2_ClusterFile.egt \
+  $OUT/gsa36_beadarrayfiles.csv
+echo "wrote $OUT/gsa36_beadarrayfiles.csv"

src/simurg/cli/main.py

@@ -300,6 +300,18 @@ def genotype_cmd(
         str | None, typer.Option("--autocall-date", help="Fixed value for GTC field 201 (reproducible output).")
     ] = None,
     threads: Annotated[int, typer.Option("--threads", min=0, max=1024, help="Worker processes (0 = auto).")] = 0,
+    controls: Annotated[
+        bool, typer.Option("--controls/--no-controls", help="Write controls_raw.tsv and controls_qc.tsv.")
+    ] = True,
+    cohort_qc: Annotated[
+        bool, typer.Option("--cohort-qc/--no-cohort-qc", help="Contamination and relatedness after genotyping.")
+    ] = True,
+    qc_config: Annotated[Path | None, typer.Option("--qc-config", help="QC threshold YAML (QC-020).")] = None,
+    qc_profile: Annotated[str | None, typer.Option("--qc-profile", help="QC profile name or 'auto'.")] = None,
+    threshold_call_rate: Annotated[float | None, typer.Option("--threshold-call-rate")] = None,
+    threshold_lrr_sd: Annotated[float | None, typer.Option("--threshold-lrr-sd")] = None,
+    threshold: Annotated[list[str] | None, typer.Option("--threshold", help="KEY=VALUE (repeatable).")] = None,
+    bg_offset: Annotated[float | None, typer.Option("--bg-offset", help="Control background offset (3000).")] = None,
     force: Annotated[bool, typer.Option("--force", help="Replace existing outputs.")] = False,
     layout: Annotated[str, typer.Option("--layout", help="nested|flat")] = "nested",
     on_duplicate_idat: Annotated[str, typer.Option("--on-duplicate-idat", help="error|first|newest")] = "error",
@@ -345,10 +357,276 @@ def genotype_cmd(
         sex_version=sex_version,
         gtc_version=gtc_version,
         autocall_date=autocall_date,
+        controls=controls,
+        cohort_qc=cohort_qc,
+        qc_config=qc_config,
+        qc_profile=qc_profile,
+        threshold_call_rate=threshold_call_rate,
+        threshold_lrr_sd=threshold_lrr_sd,
+        thresholds=list(threshold or []),
+        bg_offset=bg_offset,
         argv=sys.argv[1:],
     )
     run_guarded(rep, lambda: run_genotype(opts, rep))
 
 
+# --------------------------------------------------------------------------------------------------------
+qc_app = typer.Typer(help="Quality control: controls, samples, contamination, relatedness, loci, thresholds.")
+app.add_typer(qc_app, name="qc")
+
+
+def _reporter(command: str, out: Path | None, log_json: bool, debug: bool, quiet: bool) -> Reporter:
+    return Reporter(
+        command=command,
+        level="debug" if debug else "info",
+        json_stderr=log_json,
+        debug=debug,
+        quiet=quiet,
+        log_file=(out / "logs" / f"{command.replace(' ', '_')}.log") if out is not None else None,
+    )
+
+
+LogJson = Annotated[bool, typer.Option("--log-json", "--json-log", help="JSON Lines log on stderr.")]
+Debug = Annotated[bool, typer.Option("--debug", help="Debug logging with tracebacks.")]
+Quiet = Annotated[bool, typer.Option("--quiet", "--no-progress", help="Suppress info messages.")]
+
+
+@qc_app.command("controls")
+def qc_controls_cmd(
+    out: Annotated[Path, typer.Option("--out", help="Output folder (required).")],
+    array_type: Annotated[str, typer.Option("--array-type", help="genotyping|methylation (required).")],
+    bpm: Annotated[Path | None, typer.Option("--bpm", help="BPM manifest (control definitions).")] = None,
+    csv: Annotated[Path | None, typer.Option("--csv", help="CSV manifest (control definitions).")] = None,
+    idats: Annotated[Path | None, typer.Option("--idats", help="Folder searched recursively for IDAT pairs.")] = None,
+    sheet: Annotated[Path | None, typer.Option("--sheet", help="Sample sheet (CSV).")] = None,
+    bg_offset: Annotated[float, typer.Option("--bg-offset", help="Background offset x (QC-003).")] = 3000.0,
+    compat: Annotated[bool, typer.Option("--compat-column-names", help="Reference-style CSV names.")] = False,
+    force: Annotated[bool, typer.Option("--force", help="Replace existing control tables.")] = False,
+    log_json: LogJson = False,
+    debug: Debug = False,
+    quiet: Quiet = False,
+) -> None:
+    """Control-probe QC from IDATs without genotyping (QC-012)."""
+    from simurg.pipeline.qc import QcControlsOptions, run_qc_controls
+
+    if array_type not in ("genotyping", "methylation"):
+        raise typer.BadParameter("--array-type must be genotyping or methylation")
+    rep = _reporter("qc controls", out, log_json, debug, quiet)
+    opts = QcControlsOptions(
+        out=out,
+        bpm=bpm,
+        csv=csv,
+        idats=idats,
+        sheet=sheet,
+        array_type=array_type,
+        bg_offset=bg_offset,
+        compat=compat,
+        force=force,
+        argv=sys.argv[1:],
+    )
+    run_guarded(rep, lambda: run_qc_controls(opts, rep))
+
+
+def _qc_cohort(
+    command: str,
+    out: Path,
+    bpm: Path,
+    gtcs: Path | None,
+    sheet: Path | None,
+    egt: Path | None,
+    threads: int,
+    force: bool,
+    log_json: bool,
+    debug: bool,
+    quiet: bool,
+    *,
+    all_pairs: bool = False,
+    pedigree_from_sheet: bool = True,
+    median: str = "reference",
+    config: Path | None = None,
+    profile: str | None = None,
+    threshold_call_rate: float | None = None,
+    threshold_lrr_sd: float | None = None,
+    threshold: list[str] | None = None,
+) -> None:
+    from simurg.pipeline.qc import QcCohortOptions, run_qc_cohort
+    from simurg.qc.thresholds import parse_cli_pairs, resolve_thresholds
+
+    if median not in ("reference", "exact"):
+        raise typer.BadParameter("--median must be reference or exact")
+    rep = _reporter(f"qc {command}", out, log_json, debug, quiet)
+
+    def body() -> int:
+        cli: dict[str, Any] = parse_cli_pairs(list(threshold or []))
+        if threshold_call_rate is not None:
+            cli["callRate"] = threshold_call_rate
+        if threshold_lrr_sd is not None:
+            cli["logRDev"] = threshold_lrr_sd
+        th = resolve_thresholds(config=config, profile=profile, cli=cli)
+        opts = QcCohortOptions(
+            command=command,
+            out=out,
+            bpm=bpm,
+            gtcs=gtcs,
+            sheet=sheet,
+            egt=egt,
+            threads=threads,
+            force=force,
+            all_pairs=all_pairs,
+            pedigree_from_sheet=pedigree_from_sheet,
+            median=median,
+            thresholds=th,
+            argv=sys.argv[1:],
+        )
+        return run_qc_cohort(opts, rep)
+
+    run_guarded(rep, body)
+
+
+GtcsOpt = Annotated[Path | None, typer.Option("--gtcs", help="Folder searched recursively for GTC files.")]
+SheetOpt = Annotated[Path | None, typer.Option("--sheet", help="Sample sheet (CSV).")]
+BpmOpt = Annotated[Path, typer.Option("--bpm", help="BPM manifest (required).")]
+OutOpt = Annotated[Path, typer.Option("--out", help="Output folder (required).")]
+ThreadsOpt = Annotated[int, typer.Option("--threads", min=0, max=1024, help="Worker processes (0 = auto).")]
+ForceOpt = Annotated[bool, typer.Option("--force", help="Replace existing outputs.")]
+MedianOpt = Annotated[str, typer.Option("--median", help="reference (default, bit-compatible) or exact.")]
+
+
+@qc_app.command("samples")
+def qc_samples_cmd(
+    out: OutOpt,
+    bpm: BpmOpt,
+    gtcs: GtcsOpt = None,
+    sheet: SheetOpt = None,
+    config: Annotated[Path | None, typer.Option("--config", help="QC threshold YAML.")] = None,
+    profile: Annotated[str | None, typer.Option("--profile", help="QC profile name or 'auto'.")] = None,
+    threshold_call_rate: Annotated[float | None, typer.Option("--threshold-call-rate")] = None,
+    threshold_lrr_sd: Annotated[float | None, typer.Option("--threshold-lrr-sd")] = None,
+    threshold: Annotated[list[str] | None, typer.Option("--threshold", help="KEY=VALUE (repeatable).")] = None,
+    all_pairs: Annotated[bool, typer.Option("--all-pairs", help="Write every sample pair.")] = False,
+    median: MedianOpt = "reference",
+    threads: ThreadsOpt = 0,
+    force: ForceOpt = False,
+    log_json: LogJson = False,
+    debug: Debug = False,
+    quiet: Quiet = False,
+) -> None:
+    """Sample QC from GTCs: het rate, contamination, relatedness, flags (sample_qc.tsv)."""
+    _qc_cohort(
+        "samples",
+        out,
+        bpm,
+        gtcs,
+        sheet,
+        None,
+        threads,
+        force,
+        log_json,
+        debug,
+        quiet,
+        all_pairs=all_pairs,
+        median=median,
+        config=config,
+        profile=profile,
+        threshold_call_rate=threshold_call_rate,
+        threshold_lrr_sd=threshold_lrr_sd,
+        threshold=threshold,
+    )
+
+
+@qc_app.command("contamination")
+def qc_contamination_cmd(
+    out: OutOpt,
+    bpm: BpmOpt,
+    gtcs: GtcsOpt = None,
+    sheet: SheetOpt = None,
+    median: MedianOpt = "reference",
+    threads: ThreadsOpt = 0,
+    force: ForceOpt = False,
+    log_json: LogJson = False,
+    debug: Debug = False,
+    quiet: Quiet = False,
+) -> None:
+    """BAF-regression contamination estimate per sample (contamination.tsv)."""
+    _qc_cohort("contamination", out, bpm, gtcs, sheet, None, threads, force, log_json, debug, quiet, median=median)
+
+
+@qc_app.command("relatedness")
+def qc_relatedness_cmd(
+    out: OutOpt,
+    bpm: BpmOpt,
+    gtcs: GtcsOpt = None,
+    sheet: SheetOpt = None,
+    pedigree_from_sheet: Annotated[
+        bool, typer.Option("--pedigree-from-sheet/--no-pedigree", help="Use Parent1/Parent2/Replicate.")
+    ] = True,
+    all_pairs: Annotated[bool, typer.Option("--all-pairs", help="Write every sample pair.")] = False,
+    threads: ThreadsOpt = 0,
+    force: ForceOpt = False,
+    log_json: LogJson = False,
+    debug: Debug = False,
+    quiet: Quiet = False,
+) -> None:
+    """KING-robust kinship, duplicates and Mendelian errors (relatedness.tsv, mendel.tsv)."""
+    _qc_cohort(
+        "relatedness",
+        out,
+        bpm,
+        gtcs,
+        sheet,
+        None,
+        threads,
+        force,
+        log_json,
+        debug,
+        quiet,
+        all_pairs=all_pairs,
+        pedigree_from_sheet=pedigree_from_sheet,
+    )
+
+
+@qc_app.command("loci")
+def qc_loci_cmd(
+    out: OutOpt,
+    bpm: BpmOpt,
+    egt: Annotated[Path, typer.Option("--egt", help="EGT cluster file (required).")],
+    gtcs: GtcsOpt = None,
+    sheet: SheetOpt = None,
+    threads: ThreadsOpt = 0,
+    force: ForceOpt = False,
+    log_json: LogJson = False,
+    debug: Debug = False,
+    quiet: Quiet = False,
+) -> None:
+    """Per-locus summary statistics across samples (locus_summary.tsv, OUT-072)."""
+    _qc_cohort("loci", out, bpm, gtcs, sheet, egt, threads, force, log_json, debug, quiet)
+
+
+@qc_app.command("suggest-thresholds")
+def qc_suggest_cmd(
+    data: Annotated[list[Path], typer.Option("--data", help="Dataset folder(s) with controls_qc.tsv.")],
+    out: Annotated[Path, typer.Option("--out", help="Output YAML file (required).")],
+    k: Annotated[float, typer.Option("--k", min=0.0, help="Spread multiplier.")] = 4.0,
+    min_samples: Annotated[int, typer.Option("--min-samples", min=1, help="Minimum samples per metric.")] = 30,
+    config: Annotated[Path | None, typer.Option("--config", help="Current thresholds (YAML).")] = None,
+    profile: Annotated[str | None, typer.Option("--profile", help="Current profile.")] = None,
+    log_json: LogJson = False,
+    debug: Debug = False,
+    quiet: Quiet = False,
+) -> None:
+    """Data-driven QC threshold suggestion (QC-026)."""
+    from simurg.pipeline.qc import run_suggest_thresholds
+    from simurg.qc.thresholds import resolve_thresholds
+
+    rep = _reporter("qc suggest-thresholds", None, log_json, debug, quiet)
+    folders = [Path(p) for d in data for p in str(d).split(",") if p]
+
+    def body() -> int:
+        th = resolve_thresholds(config=config, profile=profile)
+        return run_suggest_thresholds(folders, out, th, rep, k=k, min_samples=min_samples)
+
+    run_guarded(rep, body)
+
+
 if __name__ == "__main__":  # pragma: no cover
     app()

src/simurg/core/errors.py

@@ -89,6 +89,30 @@ CATALOGUE: dict[str, MessageDef] = dict(
         _m("E_EGT_LOAD", "error", "Cluster (EGT) file cannot be loaded", "IN-044"),
         _m("W_EGT_NAME", "warning", "Cluster file name is not portable", "IN-045"),
         _m("W_MANIFEST_NO_REFSTRAND", "warning", "Manifest lacks RefStrand annotation", "IN-040"),
+        # --- QC (QC-001 .. QC-064) -----------------------------------------------------------------
+        _m("W_CTRL_MISSING", "warning", "Controls needed by a metric are absent; metric not evaluable", "QC-007"),
+        _m("W_CTRL_ZERO_DIVISOR", "warning", "Control metric has a zero denominator (NaN)", "QC-008"),
+        _m("W_CTRL_DEFS_DIFFER", "warning", "BPM and CSV control definitions differ; the BPM wins", "QC-001"),
+        _m("E_QC_METHYLATION", "error", "Methylation control metrics are not available yet (v2)", "QC-006"),
+        _m("E_QC_NO_DATASETS", "error", "No QC dataset folder found", "QC-041"),
+        _m("E_QC_DATASET_INCOMPLETE", "error", "Dataset folder lacks a required file", "QC-041"),
+        _m("E_QC_LABEL_COUNT", "error", "Number of labels does not match number of datasets", "QC-042"),
+        _m("E_QC_LABEL_DUPLICATE", "error", "Dataset labels are not unique", "QC-042"),
+        _m("W_QC_LABEL_SANITIZED", "warning", "Dataset label was sanitized", "QC-042"),
+        _m("E_QC_MIXED_TYPES", "error", "Genotyping and methylation datasets cannot be mixed", "QC-043"),
+        _m("W_QC_UNMATCHED_SAMPLES", "warning", "Sample IDs present in one QC file but not the other", "QC-040"),
+        _m("E_QC_FORMAT", "error", "Invalid --output-format value", "QC-063"),
+        _m("E_QC_CONFIG_MISSING", "error", "QC config file does not exist", "QC-021"),
+        _m("E_QC_CONFIG_EXT", "error", "QC config file must be .yaml or .yml", "QC-021"),
+        _m("E_QC_CONFIG_PARSE", "error", "QC config file is not valid YAML", "QC-021"),
+        _m("E_QC_CONFIG_STRUCTURE", "error", "QC config lacks the qc/report (QC/Report) wrapper", "QC-021"),
+        _m("E_QC_THRESHOLD_KEY", "error", "Unknown QC threshold key", "QC-021"),
+        _m("E_QC_THRESHOLD_VALUE", "error", "QC threshold value is not numeric or null", "QC-021"),
+        _m("E_QC_THRESHOLD_RANGE", "error", "QC threshold value out of range", "QC-021"),
+        _m("E_QC_PROFILE", "error", "Unknown QC profile", "QC-024"),
+        _m("W_QC_PEDIGREE", "warning", "Parent1/Parent2/Replicate value does not name a sample in the set", "IN-018"),
+        _m("W_QC_NO_BAF", "warning", "GTC has no B-allele frequencies; contamination not estimated", "QC-064"),
+        _m("W_QC_TOO_FEW", "warning", "Too few samples or sites for a cohort QC statistic", "QC-026"),
     ]
 )
 

src/simurg/formats/discovery.py

@@ -194,3 +194,19 @@ def resolve_samples(disc: Discovery, sheet: SampleSheet | None) -> tuple[list[Re
     if extra:
         notes.append(("I_UNLISTED_IDATS", f"{len(extra)} IDAT pair(s) not listed in the sample sheet are ignored"))
     return out, notes
+
+
+def discover_gtcs(root: str | Path) -> dict[str, Path]:
+    """``*.gtc`` under ``root`` keyed by file-name stem (IN-030); duplicate keys → ``E_DUP_GTC``."""
+    root = Path(root)
+    if not root.exists():
+        raise SimurgError("E_FILE_NOT_FOUND", f"GTC folder not found: {root}", path=str(root))
+    out: dict[str, Path] = {}
+    for p in walk_files(root, (".gtc",)):
+        key = p.name[: -len(".gtc")]
+        if key in out:
+            raise SimurgError("E_DUP_GTC", f"GTC files {out[key]} and {p} resolve to the same sample {key!r}")
+        out[key] = p
+    if not out:
+        raise SimurgError("E_NO_INPUT", f"no GTC files found under {root}")
+    return out

src/simurg/pipeline/genotype.py

@@ -22,6 +22,7 @@ from simurg.core.errors import ExitCode, SimurgError
 from simurg.core.log import Reporter
 from simurg.core.provenance import RunManifest
 from simurg.formats.bpm import read_bpm
+from simurg.formats.csv_manifest import assay_format
 from simurg.formats.egt import read_egt
 from simurg.formats.gtc import write_gtc
 from simurg.formats.idat import check_pair, read_idat
@@ -37,6 +38,24 @@ from simurg.genotype.sample import (
 )
 from simurg.genotype.sex import SexParams
 from simurg.pipeline.inputs import check_paths, resolve_sample_set, sample_dir, sheet_values
+from simurg.pipeline.qc import (
+    CONTROLS_QC,
+    CONTROLS_RAW,
+    MENDEL_COLUMNS,
+    RELATEDNESS_COLUMNS,
+    ControlSample,
+    cohort_qc,
+    control_sample,
+    default_tmp,
+    load_cohort,
+    load_control_beads,
+    pedigree_from_sheet,
+    write_control_tables,
+    write_tsv,
+)
+from simurg.qc.controls import ControlBead
+from simurg.qc.status import control_status, functional_status, het_zscores, sample_flags
+from simurg.qc.thresholds import Thresholds, parse_cli_pairs, resolve_thresholds
 from simurg.signal.normalization.apply import normalize, polar
 from simurg.signal.raw_xy import LocusMapCache, check_coverage
 
@@ -67,7 +86,7 @@ FIXED_COLUMNS = [
     "het_rate",
     "contamination",
 ]
-QC_COLUMNS = ["functional_status", "control_status", "control_flag_codes", "flags"]  # filled from phase P5
+QC_COLUMNS = ["functional_status", "control_status", "control_flag_codes", "flags"]  # QC-030/031/033
 FLOAT_COLUMNS = {"autosomal_call_rate", "call_rate", "logr_dev", "gc10", "gc50", "het_rate", "contamination"}
 
 
@@ -90,6 +109,14 @@ class GenotypeOptions:
     gtc_version: int | None = None
     autocall_date: str | None = None
     assay_control_probes: tuple[str, ...] = ("TGA_Ctrl_5716",)
+    controls: bool = True  # OUT-070
+    cohort_qc: bool = True  # contamination + relatedness after genotyping (QC-064)
+    qc_config: Path | None = None
+    qc_profile: str | None = None
+    threshold_call_rate: float | None = None
+    threshold_lrr_sd: float | None = None
+    thresholds: list[str] = field(default_factory=list)
+    bg_offset: float | None = None
     argv: list[str] = field(default_factory=list)
 
 
@@ -100,6 +127,7 @@ class Job:
     grn: str
     red: str
     gtc_path: str
+    sheet_values: dict[str, str] = field(default_factory=dict)
 
 
 _W: dict[str, Any] = {}
@@ -111,8 +139,12 @@ def _init(
     auto_1_22: np.ndarray[Any, np.dtype[np.bool_]],
     probe_idx: dict[str, int],
     started: _dt.datetime,
+    beads: list[ControlBead] | None = None,
+    bg_offset: float = 3000.0,
 ) -> None:
-    _W.update(ctx=ctx, params=params, auto=auto_1_22, probes=probe_idx, started=started)
+    _W.update(
+        ctx=ctx, params=params, auto=auto_1_22, probes=probe_idx, started=started, beads=beads, bg_offset=bg_offset
+    )
 
 
 def _process(job: Job) -> dict[str, Any]:
@@ -134,6 +166,11 @@ def _process(job: Job) -> dict[str, Any]:
         write_gtc(g, Path(job.gtc_path))
         auto = _W["auto"] & ~ctx.manifest.intensity_only
         n_auto = int(auto.sum())
+        g_auto = g.genotypes[auto]
+        n_called = int((g_auto > 0).sum())
+        row["het_rate"] = float((g_auto == 2).sum() / n_called) if n_called else math.nan
+        if _W.get("beads"):
+            row["control"] = control_sample(job.sample_id, grn, red, _W["beads"], _W["bg_offset"], job.sheet_values)
         row.update(
             autosomal_call_rate=float(np.float32(int((g.genotypes[auto] > 0).sum())) / np.float32(n_auto))
             if n_auto
@@ -233,10 +270,22 @@ def run_genotype(opts: GenotypeOptions, rep: Reporter) -> int:
             continue
         stem = sset.stems[s.sample_id]
         gtc = sample_dir(opts.out, stem, opts.layout) / f"{stem}.gtc"
-        jobs.append(Job(s.sample_id, s.sample_name, str(s.pair.grn), str(s.pair.red), str(gtc)))
+        jobs.append(
+            Job(
+                s.sample_id,
+                s.sample_name,
+                str(s.pair.grn),
+                str(s.pair.red),
+                str(gtc),
+                {k: str(v) for k, v in sheet_values(s, sset.sheet).items()},
+            )
+        )
 
     summary_tsv, summary_json = opts.out / "sample_summary.tsv", opts.out / "sample_summary.json"
-    check_outputs_absent([summary_tsv, summary_json] + [Path(j.gtc_path) for j in jobs], force=opts.force)
+    extra = [opts.out / CONTROLS_RAW, opts.out / CONTROLS_QC] if opts.controls else []
+    if opts.cohort_qc:
+        extra.append(opts.out / "relatedness.tsv")
+    check_outputs_absent([summary_tsv, summary_json, *extra] + [Path(j.gtc_path) for j in jobs], force=opts.force)
 
     rep.info(f"reading manifest {opts.bpm} and cluster file {opts.egt}")
     manifest = read_bpm(opts.bpm)
@@ -256,6 +305,10 @@ def run_genotype(opts: GenotypeOptions, rep: Reporter) -> int:
             "W_EGT_NAME", f"cluster file was built for {egt.manifest_name!r}, manifest is {manifest.manifest_name!r}"
         )
     clusters = align_clusters(egt, manifest)
+    th = _thresholds(opts, manifest)
+    for rec in th.as_records():
+        rm.set_option(f"qc.{rec['key']}", rec["value"], rec["source"])
+    beads = load_control_beads(manifest, opts.csv, rep) if opts.controls else []
     ctx = manifest_context(manifest, clusters, cluster_file_name=Path(opts.egt).name)
     lite = manifest.lite()
     ctx = replace(ctx, manifest=lite, cache=LocusMapCache(lite))  # workers get no string columns
@@ -282,14 +335,14 @@ def run_genotype(opts: GenotypeOptions, rep: Reporter) -> int:
     rep.info(f"genotyping {len(jobs)} sample(s) with {workers} worker(s)")
     try:
         if workers == 1:
-            _init(ctx, params, auto_1_22, probes, started)
+            _init(ctx, params, auto_1_22, probes, started, beads, th.bg_offset)
             results = [_process(j) for j in jobs]
         else:
             with ProcessPoolExecutor(
                 workers,
                 mp_context=mp.get_context("spawn"),
                 initializer=_init,
-                initargs=(ctx, params, auto_1_22, probes, started),
+                initargs=(ctx, params, auto_1_22, probes, started, beads, th.bg_offset),
             ) as ex:
                 results = list(ex.map(_process, jobs))
     except KeyboardInterrupt:
@@ -315,14 +368,11 @@ def run_genotype(opts: GenotypeOptions, rep: Reporter) -> int:
         rows[r["sample_id"]].update(r)
 
     probe_cols = [f"{p}_norm_r" for p in probes]
+    _run_qc(opts, rep, rm, th, beads, sset, rows, results, jobs, manifest, probe_cols, params.gtc_version)
     sheet_cols = [f"sheet:{c}" for c in (sset.sheet.passthrough_columns() if sset.sheet else [])]
     out_rows = []
     for s in sset.samples:
-        r = rows[s.sample_id]
+        out_rows.append(rows[s.sample_id])
-        known = {"M", "F"}
-        est, rep_sex = r.get("sex_est", "U"), r.get("sex_reported", "U")
-        r["sex_check"] = ("match" if est == rep_sex else "mismatch") if est in known and rep_sex in known else "unknown"
-        out_rows.append(r)
     columns = FIXED_COLUMNS + probe_cols + QC_COLUMNS + sheet_cols
     header = [c.removeprefix("sheet:") for c in columns]
     lines = ["\t".join(header)] + ["\t".join(_fmt(c, r.get(c, "")) for c in columns) for r in out_rows]
@@ -350,3 +400,106 @@ def run_genotype(opts: GenotypeOptions, rep: Reporter) -> int:
     if ok == len(out_rows):
         return ExitCode.OK
     return ExitCode.PARTIAL if ok else ExitCode.ALL_FAILED
+
+
+def _thresholds(opts: GenotypeOptions, manifest: Any) -> Thresholds:
+    cli: dict[str, Any] = parse_cli_pairs(opts.thresholds)
+    if opts.threshold_call_rate is not None:
+        cli["callRate"] = opts.threshold_call_rate
+    if opts.threshold_lrr_sd is not None:
+        cli["logRDev"] = opts.threshold_lrr_sd
+    if opts.bg_offset is not None:
+        cli["bgOffset"] = opts.bg_offset
+    return resolve_thresholds(
+        config=opts.qc_config,
+        profile=opts.qc_profile,
+        cli=cli,
+        array_type="genotyping",
+        assay_format=assay_format(manifest),
+        manifest_name=manifest.manifest_name,
+    )
+
+
+def _run_qc(
+    opts: GenotypeOptions,
+    rep: Reporter,
+    rm: RunManifest,
+    th: Thresholds,
+    beads: list[ControlBead],
+    sset: Any,
+    rows: dict[str, dict[str, Any]],
+    results: list[dict[str, Any]],
+    jobs: list[Job],
+    manifest: Any,
+    probe_cols: list[str],
+    gtc_version: int,
+) -> None:
+    """Control tables (OUT-070), het rate, cohort QC and the QC-030/031/033 columns of the summary."""
+    controls: dict[str, ControlSample] = {r["sample_id"]: r.pop("control") for r in results if "control" in r}
+    for r in rows.values():  # results were merged into rows already; drop the heavy objects there too
+        r.pop("control", None)
+    if opts.controls and beads:
+        assay = {
+            sid: min((float(rows[sid].get(c, math.nan)) for c in probe_cols), default=math.nan) for sid in controls
+        }
+        ordered = [controls[s.sample_id] for s in sset.samples if s.sample_id in controls]
+        for p in write_control_tables(
+            opts.out, ordered, beads, th.bg_offset, rep, assay_control=assay if probe_cols else None
+        ):
+            rm.add_output(p, "simurg.controls/1")
+    ids = [s.sample_id for s in sset.samples]
+    het_z = dict(zip(ids, het_zscores([rows[i].get("het_rate", math.nan) for i in ids]), strict=True))
+    dups: dict[str, list[str]] = {}
+    unexpected: set[str] = set()
+    ok_jobs = [j for j in jobs if rows[j.sample_id].get("status") == "OK"]
+    if opts.cohort_qc and len(ok_jobs) >= 2:
+        tmp = default_tmp(opts.out)
+        rep.info(f"cohort QC (contamination, relatedness) on {len(ok_jobs)} sample(s)")
+        store, infos = load_cohort(
+            [(j.sample_id, j.gtc_path) for j in ok_jobs],
+            manifest.num_loci,
+            tmp,
+            baf=gtc_version >= 4,
+            scores=False,
+            workers=max(1, opts.threads or (os.cpu_count() or 1)),
+        )
+        try:
+            ped = pedigree_from_sheet(sset.sheet, [j.sample_id for j in ok_jobs], rep)
+            res = cohort_qc(
+                store, infos, manifest, rep, pedigree=ped, workers=max(1, opts.threads or (os.cpu_count() or 1))
+            )
+        finally:
+            store.close()
+        if res.contamination is not None:
+            for i, s in enumerate(infos):
+                rows[s.sample_id]["contamination"] = float(res.contamination.estimate[i])
+        dups, unexpected = res.duplicates, res.unexpected
+        write_tsv(opts.out / "relatedness.tsv", RELATEDNESS_COLUMNS, res.pairs)
+        rm.add_output(opts.out / "relatedness.tsv", "simurg.relatedness/1")
+        if res.mendel:
+            write_tsv(opts.out / "mendel.tsv", MENDEL_COLUMNS, res.mendel)
+            rm.add_output(opts.out / "mendel.tsv", "simurg.mendel/1")
+    for sid in ids:
+        r = rows[sid]
+        known = {"M", "F"}
+        est, rep_sex = r.get("sex_est", "U"), r.get("sex_reported", "U")
+        r["sex_check"] = ("match" if est == rep_sex else "mismatch") if est in known and rep_sex in known else "unknown"
+        if r.get("status") == "MISSING_INPUT":
+            continue
+        r["functional_status"] = functional_status(r, th)
+        if sid in controls:
+            vals = dict(controls[sid].qc.values)
+            if probe_cols:
+                vals["assay_control_norm_r"] = min(float(r.get(c, math.nan)) for c in probe_cols)
+            cs, codes = control_status(vals, controls[sid].qc.na_metrics, th)
+            r["control_status"], r["control_flag_codes"] = cs, ";".join(codes)
+        r["flags"] = ";".join(
+            sample_flags(
+                r,
+                th,
+                het_z=het_z.get(sid, math.nan),
+                probe_columns=probe_cols,
+                duplicates=dups.get(sid, []),
+                unexpected_relatedness=sid in unexpected,
+            )
+        )

src/simurg/pipeline/qc.py

@@ -0,0 +1,915 @@
+"""QC orchestration and I/O (spec QC-001 … QC-064, CLI-024, CLI-025, OUT-070, OUT-072; plan §9).
+
+The science lives in ``simurg.qc``; this module reads IDATs/GTCs, streams cohort matrices through disk-backed
+arrays (so thousands of samples fit in a small memory budget) and writes the QC tables.
+"""
+
+from __future__ import annotations
+
+import datetime as _dt
+import math
+import multiprocessing as mp
+import os
+import shutil
+import tempfile
+from collections.abc import Iterator, Mapping, Sequence
+from concurrent.futures import ProcessPoolExecutor
+from dataclasses import dataclass, field
+from pathlib import Path
+from typing import IO, Any, cast
+
+import numpy as np
+import numpy.typing as npt
+
+from simurg.core.atomic import atomic_write_text, check_outputs_absent
+from simurg.core.errors import ExitCode, SimurgError
+from simurg.core.log import Reporter
+from simurg.core.provenance import RunManifest
+from simurg.formats.bpm import Manifest
+from simurg.formats.idat import Idat
+from simurg.formats.samplesheet import SampleSheet, fold
+from simurg.genotype.baf_lrr import chrom_core
+from simurg.qc.contamination import ContaminationResult, baf_regress
+from simurg.qc.controls import (
+    ASSAY_CONTROL_COLUMN,
+    COMPAT_METRIC_NAMES,
+    ControlBead,
+    ControlQC,
+    beads_from_definitions,
+    compare_definitions,
+    compute_control_metrics,
+    control_intensities,
+    metrics_for,
+)
+from simurg.qc.relatedness import DEGREE_CUTS, DUPLICATE_CONCORDANCE, degree, mendel_errors, pair_stats, select_sites
+from simurg.qc.status import het_zscores
+from simurg.qc.thresholds import Thresholds
+
+CONTROLS_RAW = "controls_raw.tsv"
+CONTROLS_QC = "controls_qc.tsv"
+COMPAT_RAW = "controls.raw_metrics.csv"
+COMPAT_QC = "controls.qc_metrics.csv"
+MIN_CONTAMINATION_SAMPLES = 10
+
+
+def fmt(v: Any, digits: int = 8) -> str:
+    if v is None or v == "":
+        return ""
+    if isinstance(v, float | np.floating):
+        f = float(v)
+        return "NaN" if math.isnan(f) else f"{f:.{digits}f}"
+    return str(v).replace("\t", " ").replace("\n", " ")
+
+
+def write_tsv(path: Path, columns: Sequence[str], rows: Sequence[Mapping[str, Any]], sep: str = "\t") -> None:
+    lines = [sep.join(columns)] + [sep.join(fmt(r.get(c, "")) for c in columns) for r in rows]
+    atomic_write_text(path, "\n".join(lines) + "\n")
+
+
+def autosome_mask(chroms: Sequence[str]) -> npt.NDArray[np.bool_]:
+    """Chromosomes 1–22 only (het rate, relatedness, contamination sites)."""
+    auto = {str(i) for i in range(1, 23)}
+    return np.array([chrom_core(c) in auto for c in chroms], dtype=np.bool_)
+
+
+# ================================================================================================ controls
+def load_control_beads(manifest: Manifest | None, csv_path: Path | None, rep: Reporter) -> list[ControlBead]:
+    """QC-001: BPM control config wins over the CSV ``[Controls]`` section; differences give one warning."""
+    bpm_beads = beads_from_definitions(manifest.controls) if manifest is not None else []
+    csv_beads: list[ControlBead] = []
+    if csv_path is not None:
+        from simurg.formats.csv_manifest import read_csv_manifest
+
+        csv_beads = beads_from_definitions(read_csv_manifest(csv_path, load_sequences=False).controls)
+    if bpm_beads and csv_beads and (diffs := compare_definitions(bpm_beads, csv_beads)):
+        rep.warn("W_CTRL_DEFS_DIFFER", "BPM and CSV control definitions differ (using the BPM): " + "; ".join(diffs))
+    return bpm_beads or csv_beads
+
+
+@dataclass
+class ControlSample:
+    sample_id: str
+    meta: dict[str, str]
+    green: npt.NDArray[np.float64]
+    red: npt.NDArray[np.float64]
+    qc: ControlQC
+
+
+def _param(pars: str, key: str) -> str:
+    from simurg.genotype.sample import _param as p
+
+    return p(pars, key)
+
+
+def idat_metadata(grn: Idat, red: Idat, sheet_values: Mapping[str, str] | None = None) -> dict[str, str]:
+    """QC-011 metadata: IDAT values, plate and well overridden by the sheet."""
+    scan = red.scan_entry() or grn.scan_entry()
+    sv = {fold(k): v for k, v in (sheet_values or {}).items()}
+    return {
+        "imaging_date": red.scan_time_iso or grn.scan_time_iso,
+        "sample_plate": sv.get(fold("Sample_Plate")) or red.plate or grn.plate,
+        "sample_well": sv.get(fold("Sample_Well")) or red.well or grn.well,
+        "scanner_id": _param(scan.block_pars, "sherlockID") if scan else "",
+        "scanner_version": scan.code_version if scan else "",
+    }
+
+
+def control_sample(
+    sample_id: str,
+    grn: Idat,
+    red: Idat,
+    beads: Sequence[ControlBead],
+    bg_offset: float,
+    sheet_values: Mapping[str, str] | None = None,
+) -> ControlSample:
+    g, r = control_intensities(grn, red, beads)
+    qc = compute_control_metrics(beads, g, r, "genotyping", bg_offset)
+    return ControlSample(sample_id, idat_metadata(grn, red, sheet_values), g, r, qc)
+
+
+def _compat_date(iso: str) -> str:
+    try:
+        t = _dt.datetime.fromisoformat(iso)
+    except ValueError:
+        return iso
+    from simurg.genotype.sample import invariant_datetime
+
+    return invariant_datetime(t)
+
+
+def write_control_tables(
+    out: Path,
+    samples: Sequence[ControlSample],
+    beads: Sequence[ControlBead],
+    bg_offset: float,
+    rep: Reporter,
+    *,
+    array_type: str = "genotyping",
+    assay_control: Mapping[str, float] | None = None,
+    compat: bool = False,
+) -> list[Path]:
+    """QC-009 and QC-011 (or their compatibility CSVs). ``assay_control`` maps sample_id → assay-control norm R."""
+    metrics = metrics_for(array_type)
+    missing: dict[str, list[str]] = {}
+    for s in samples:
+        missing.update(s.qc.missing)
+        for col in s.qc.zero_divisor:
+            rep.warn("W_CTRL_ZERO_DIVISOR", f"{s.sample_id}: {col} has a zero denominator", sample_id=s.sample_id)
+    if missing:
+        ids = sorted({c for v in missing.values() for c in v})
+        rep.warn(
+            "W_CTRL_MISSING",
+            f"manifest lacks control(s) {', '.join(ids)}; not evaluable: {', '.join(sorted(missing))}",
+        )
+    raw_rows = [
+        {
+            "sample_id": s.sample_id,
+            "category": b.category,
+            "bead_address": b.address,
+            "control_name": b.name,
+            "color_label": b.color,
+            "green_intensity": "NaN" if math.isnan(gv) else int(gv),
+            "red_intensity": "NaN" if math.isnan(rv) else int(rv),
+        }
+        for s in samples
+        for b, gv, rv in zip(beads, s.green.tolist(), s.red.tolist(), strict=True)
+    ]
+    qc_rows = []
+    for s in samples:
+        row: dict[str, Any] = {
+            "sample_id": s.sample_id,
+            "array_type": array_type,
+            **s.meta,
+            "bg_offset": float(bg_offset),
+            **{m.column: s.qc.values.get(m.column, math.nan) for m in metrics},
+            "na_metrics": ";".join(s.qc.na_metrics),
+        }
+        if assay_control is not None:
+            row[ASSAY_CONTROL_COLUMN] = assay_control.get(s.sample_id, math.nan)
+        qc_rows.append(row)
+    meta_cols = ["imaging_date", "sample_plate", "sample_well", "scanner_id", "scanner_version"]
+    qc_cols = ["sample_id", "array_type", *meta_cols, "bg_offset", *(m.column for m in metrics)]
+    if assay_control is not None:
+        qc_cols.append(ASSAY_CONTROL_COLUMN)
+    qc_cols.append("na_metrics")
+    if not compat:
+        raw_cols = list(raw_rows[0]) if raw_rows else ["sample_id"]
+        write_tsv(out / CONTROLS_RAW, raw_cols, raw_rows)
+        write_tsv(out / CONTROLS_QC, qc_cols, qc_rows)
+        return [out / CONTROLS_RAW, out / CONTROLS_QC]
+    ren_raw = {
+        "sample_id": "SampleId",
+        "category": "Category",
+        "bead_address": "Beadtype",
+        "control_name": "Control",
+        "color_label": "Color",
+        "green_intensity": "GreenIntensity",
+        "red_intensity": "RedIntensity",
+    }
+    ren_qc = {
+        "sample_id": "SampleId",
+        "imaging_date": "ImagingDate",
+        "sample_plate": "SamplePlate",
+        "sample_well": "SampleWell",
+        "scanner_id": "ScannerId",
+        "scanner_version": "ScannerVersion",
+        **COMPAT_METRIC_NAMES,
+    }
+    raw_c = [{ren_raw[k]: v for k, v in r.items()} for r in raw_rows]
+    qc_keep = [c for c in qc_cols if c in ren_qc]
+    qc_c = [
+        {ren_qc[c]: (_compat_date(r[c]) if c == "imaging_date" else r.get(c, "")) for c in qc_keep} for r in qc_rows
+    ]
+    write_tsv(out / COMPAT_RAW, list(ren_raw.values()), raw_c, sep=",")
+    write_tsv(out / COMPAT_QC, [ren_qc[c] for c in qc_keep], qc_c, sep=",")
+    return [out / COMPAT_RAW, out / COMPAT_QC]
+
+
+@dataclass
+class QcControlsOptions:
+    out: Path
+    bpm: Path | None = None
+    csv: Path | None = None
+    idats: Path | None = None
+    sheet: Path | None = None
+    array_type: str = "genotyping"
+    bg_offset: float = 3000.0
+    compat: bool = False
+    force: bool = False
+    on_duplicate_idat: str = "error"
+    argv: list[str] = field(default_factory=list)
+
+
+def run_qc_controls(opts: QcControlsOptions, rep: Reporter) -> int:
+    """``simurg qc controls`` — control tables from IDATs without genotyping (QC-012, QC-014)."""
+    from simurg.formats.bpm import read_bpm
+    from simurg.formats.idat import check_pair, read_idat
+    from simurg.pipeline.inputs import check_paths, resolve_sample_set, sheet_values
+
+    rm = RunManifest(command="qc controls", argv=opts.argv)
+    rep.run_id = rm.run_id
+    if opts.array_type == "methylation":
+        raise SimurgError("E_QC_METHYLATION", "methylation control metrics are planned for v2 (QC-006)")
+    if opts.bpm is None and opts.csv is None:
+        raise SimurgError("E_NO_INPUT", "give --bpm or --csv for the control definitions (QC-001)")
+    if opts.bg_offset < 0:
+        raise SimurgError("E_QC_THRESHOLD_RANGE", f"--bg-offset must be >= 0, got {opts.bg_offset}")
+    check_paths([(opts.bpm, "BPM"), (opts.csv, "CSV manifest"), (opts.idats, "IDAT folder"), (opts.sheet, "sheet")])
+    sset = resolve_sample_set(opts.idats, opts.sheet, opts.on_duplicate_idat, rep, rm)
+    names = (COMPAT_RAW, COMPAT_QC) if opts.compat else (CONTROLS_RAW, CONTROLS_QC)
+    targets = [opts.out / n for n in names]
+    check_outputs_absent(targets, force=opts.force)  # QC-014: other files in the folder are fine
+    manifest = read_bpm(opts.bpm) if opts.bpm else None
+    beads = load_control_beads(manifest, opts.csv, rep)
+    if not beads:
+        raise SimurgError("E_BPM_FORMAT", "the manifest has no control definitions")
+    samples: list[ControlSample] = []
+    for s in sset.samples:
+        if s.pair is None or not s.pair.complete:
+            rep.warn("W_IDAT_UNPAIRED", f"{s.sample_id}: no complete IDAT pair; skipped", sample_id=s.sample_id)
+            continue
+        assert s.pair.grn is not None and s.pair.red is not None
+        grn, red = read_idat(s.pair.grn), read_idat(s.pair.red)
+        check_pair(grn, red)
+        samples.append(control_sample(s.sample_id, grn, red, beads, opts.bg_offset, sheet_values(s, sset.sheet)))
+    # QC-014: sample IDs of an existing summary must match
+    summary = opts.out / "sample_summary.tsv"
+    if summary.exists():
+        ids = {ln.split("\t", 1)[0] for ln in summary.read_text(encoding="utf-8").splitlines()[1:] if ln}
+        mine = {s.sample_id for s in samples}
+        if ids ^ mine:
+            rep.warn("W_QC_UNMATCHED_SAMPLES", f"{len(ids ^ mine)} sample ID(s) differ from {summary}")
+    opts.out.mkdir(parents=True, exist_ok=True)
+    written = write_control_tables(opts.out, samples, beads, opts.bg_offset, rep, compat=opts.compat)
+    for p in written:
+        rm.add_output(p, "simurg.controls/1")
+    rm.summary = {"samples": len(samples)}
+    rm.status = "completed"
+    rm.write(opts.out / "run_manifest.qc_controls.json")
+    rep.info(f"control QC for {len(samples)} sample(s) written to {opts.out}")
+    return ExitCode.OK
+
+
+# ================================================================================================ cohort store
+@dataclass
+class SampleInfo:
+    sample_id: str
+    path: str
+    sample_name: str = ""
+    call_rate: float = math.nan
+    logr_dev: float = math.nan
+    sex_est: str = "U"
+    has_baf: bool = False
+    error: str = ""
+
+
+class CohortStore:
+    """Disk-backed (samples × loci) matrices: genotypes (uint8), BAF and GenCall scores (float32, optional)."""
+
+    def __init__(self, root: Path, n_samples: int, n_loci: int, *, baf: bool, scores: bool) -> None:
+        self.root, self.shape = root, (n_samples, n_loci)
+        root.mkdir(parents=True, exist_ok=True)
+        self.gt = np.lib.format.open_memmap(root / "gt.npy", "w+", np.uint8, self.shape)
+        self.baf = np.lib.format.open_memmap(root / "baf.npy", "w+", np.float32, self.shape) if baf else None
+        self.score = np.lib.format.open_memmap(root / "score.npy", "w+", np.float32, self.shape) if scores else None
+
+    def close(self) -> None:
+        for name in ("gt", "baf", "score"):
+            a = getattr(self, name)
+            if a is not None:
+                a.flush()
+        self.gt = self.baf = self.score = None  # type: ignore[assignment]
+        shutil.rmtree(self.root, ignore_errors=True)
+
+
+def _load_one(task: tuple[int, str, str, str, bool, bool]) -> SampleInfo:
+    from simurg.formats.gtc import read_gtc
+
+    i, sid, path, root, want_baf, want_scores = task
+    info = SampleInfo(sid, path)
+    try:
+        g = read_gtc(path)
+        info.sample_name, info.call_rate, info.logr_dev, info.sex_est = (
+            g.sample_name,
+            float(g.call_rate),
+            float(g.logr_dev) if g.version >= 4 else math.nan,
+            g.gender_char,
+        )
+        gt = np.load(Path(root) / "gt.npy", mmap_mode="r+")
+        if len(g.genotypes) != gt.shape[1]:
+            raise SimurgError("E_GTC_UNSUITABLE", f"{path}: {len(g.genotypes)} loci, manifest has {gt.shape[1]}")
+        gt[i] = g.genotypes
+        gt.flush()
+        if want_baf:
+            baf = np.load(Path(root) / "baf.npy", mmap_mode="r+")
+            info.has_baf = len(g.b_allele_freqs) == gt.shape[1]
+            baf[i] = g.b_allele_freqs if info.has_baf else np.nan
+            baf.flush()
+        if want_scores:
+            sc = np.load(Path(root) / "score.npy", mmap_mode="r+")
+            sc[i] = g.genotype_scores
+            sc.flush()
+    except SimurgError as e:
+        info.error = f"{e.code}: {e.message}"
+    except Exception as e:
+        info.error = f"E_INTERNAL: {type(e).__name__}: {e}"
+    return info
+
+
+def load_cohort(
+    samples: Sequence[tuple[str, str]], n_loci: int, tmp: Path, *, baf: bool, scores: bool, workers: int
+) -> tuple[CohortStore, list[SampleInfo]]:
+    """Read GTCs (sample_id, path) into a :class:`CohortStore` with a process pool."""
+    store = CohortStore(tmp, len(samples), n_loci, baf=baf, scores=scores)
+    store.gt.flush()
+    tasks = [(i, sid, p, str(tmp), baf, scores) for i, (sid, p) in enumerate(samples)]
+    if workers <= 1 or len(tasks) <= 1:
+        infos = [_load_one(t) for t in tasks]
+    else:
+        with ProcessPoolExecutor(workers, mp_context=mp.get_context("spawn")) as ex:
+            infos = list(ex.map(_load_one, tasks, chunksize=4))
+    return store, infos
+
+
+class SiteView:
+    """(sites × samples) row access onto a (samples × loci) matrix, optionally restricted to some samples."""
+
+    def __init__(self, m: npt.NDArray[Any], samples: npt.NDArray[np.intp] | None = None) -> None:
+        self.m, self.samples = m, samples
+        self.shape = (m.shape[1], m.shape[0] if samples is None else len(samples))
+
+    def __getitem__(self, rows: Any) -> npt.NDArray[Any]:
+        block = np.asarray(self.m[:, rows])
+        return (block if self.samples is None else block[self.samples]).T
+
+
+def _blocks(n: int, size: int) -> Iterator[slice]:
+    for s in range(0, n, size):
+        yield slice(s, min(n, s + size))
+
+
+def genotype_counts(
+    gt: npt.NDArray[np.uint8], samples: npt.NDArray[np.intp] | None = None, block: int = 65536
+) -> npt.NDArray[np.int64]:
+    """(loci × 4) counts of no-call/AA/AB/BB over (a subset of) samples, streamed in locus blocks."""
+    out = np.zeros((gt.shape[1], 4), np.int64)
+    for sl in _blocks(gt.shape[1], block):
+        g = np.asarray(gt[:, sl])
+        if samples is not None:
+            g = g[samples]
+        for k in range(4):
+            out[sl, k] = (g == k).sum(axis=0)
+    return out
+
+
+def het_rates(gt: npt.NDArray[np.uint8], loci: npt.NDArray[np.intp], block: int = 65536) -> npt.NDArray[np.float64]:
+    het = np.zeros(gt.shape[0], np.int64)
+    called = np.zeros(gt.shape[0], np.int64)
+    for s in range(0, len(loci), block):
+        g = np.asarray(gt[:, loci[s : s + block]])
+        het += (g == 2).sum(axis=1)
+        called += (g > 0).sum(axis=1)
+    with np.errstate(invalid="ignore", divide="ignore"):
+        return np.where(called > 0, het / called, np.nan)
+
+
+# ================================================================================================ pedigree
+@dataclass
+class Pedigree:
+    parents: dict[str, list[str]]  # child sample_id -> genotyped parents (0–2)
+    replicates: dict[str, str]  # sample_id -> declared replicate sample_id
+
+    def expected(self, a: str, b: str) -> str:
+        if b in self.parents.get(a, []) or a in self.parents.get(b, []):
+            return "parent-offspring"
+        if self.replicates.get(a) == b or self.replicates.get(b) == a:
+            return "replicate"
+        pa, pb = set(self.parents.get(a, [])), set(self.parents.get(b, []))
+        if pa & pb:
+            return "sibling"
+        return ""
+
+
+def pedigree_from_sheet(sheet: SampleSheet | None, sample_ids: Sequence[str], rep: Reporter) -> Pedigree:
+    """IN-018: Parent1/Parent2/Replicate values name Sample_ID, Sample_Name or sample_id (in that order)."""
+    ped = Pedigree({}, {})
+    if sheet is None:
+        return ped
+    ids = set(sample_ids)
+    by_col: dict[str, dict[str, str]] = {"sample_id_col": {}, "name": {}}
+    rows = [r for r in sheet.rows if r.sample_id in ids]
+    cols = {fold(c): c for c in sheet.columns}
+    for r in rows:
+        if (c := cols.get(fold("Sample_ID"))) and r.values.get(c):
+            by_col["sample_id_col"].setdefault(r.values[c], r.sample_id)
+        if r.sample_name:
+            by_col["name"].setdefault(r.sample_name, r.sample_id)
+
+    def resolve(v: str) -> str | None:
+        v = v.strip()
+        if v in ("", "0", "NA", "na", "-"):
+            return None
+        for table in (by_col["sample_id_col"], by_col["name"]):
+            if v in table:
+                return table[v]
+        if v in ids:
+            return v
+        rep.warn("W_QC_PEDIGREE", f"{v!r} does not name a sample in this set")
+        return None
+
+    for r in rows:
+        for key in ("Parent1", "Parent2"):
+            if (c := cols.get(fold(key))) and (p := resolve(r.values.get(c, ""))) and p != r.sample_id:
+                ped.parents.setdefault(r.sample_id, []).append(p)
+        if (c := cols.get(fold("Replicate"))) and (p := resolve(r.values.get(c, ""))) and p != r.sample_id:
+            ped.replicates[r.sample_id] = p
+    return ped
+
+
+# ================================================================================================ cohort QC
+@dataclass
+class CohortQC:
+    het_rate: npt.NDArray[np.float64]
+    het_z: list[float]
+    contamination: ContaminationResult | None = None
+    pairs: list[dict[str, Any]] = field(default_factory=list)
+    duplicates: dict[str, list[str]] = field(default_factory=dict)
+    unexpected: set[str] = field(default_factory=set)
+    mendel: list[dict[str, Any]] = field(default_factory=list)
+
+
+def cohort_qc(
+    store: CohortStore,
+    infos: Sequence[SampleInfo],
+    manifest: Manifest,
+    rep: Reporter,
+    *,
+    pedigree: Pedigree | None = None,
+    contamination: bool = True,
+    relatedness: bool = True,
+    all_pairs: bool = False,
+    median: str = "reference",
+    workers: int = 1,
+) -> CohortQC:
+    ids = [s.sample_id for s in infos]
+    ok = np.array([not s.error for s in infos])
+    auto = np.flatnonzero(autosome_mask(manifest.chroms) & ~manifest.intensity_only)
+    het = het_rates(store.gt, auto)
+    het[~ok] = np.nan
+    res = CohortQC(het_rate=het, het_z=het_zscores(het.tolist()))
+    if contamination and store.baf is not None:
+        usable = ok & np.array([s.has_baf for s in infos])
+        for s in infos:
+            if not s.error and not s.has_baf:
+                rep.warn("W_QC_NO_BAF", f"{s.sample_id}: GTC has no BAF; contamination not estimated")
+        if usable.sum() < MIN_CONTAMINATION_SAMPLES:
+            rep.warn(
+                "W_QC_TOO_FEW",
+                f"contamination needs >= {MIN_CONTAMINATION_SAMPLES} samples with BAF (have {int(usable.sum())})",
+            )
+        else:
+            rows_idx = np.flatnonzero(usable)
+            sel = None if usable.all() else rows_idx
+            c = genotype_counts(store.gt, sel)[:, 1:].astype(np.float64)
+            an = 2 * c.sum(axis=1)
+            site_af = np.full(manifest.num_loci, np.nan)
+            with np.errstate(invalid="ignore", divide="ignore"):
+                af = np.where(an > 0, (c[:, 1] + 2 * c[:, 2]) / an, np.nan)
+            site_af[auto] = af[auto]  # contamination sites: autosomes only
+            with _pool(workers) as ex:
+                cr = baf_regress(
+                    SiteView(store.gt, sel),  # type: ignore[arg-type]
+                    SiteView(store.baf, sel),  # type: ignore[arg-type]
+                    site_af,
+                    median=median,
+                    map_fn=ex.map if ex is not None else map,
+                )
+            est = np.full(len(ids), np.nan)
+            nh = np.zeros(len(ids), np.int64)
+            est[rows_idx], nh[rows_idx] = cr.estimate, cr.n_hom
+            res.contamination = ContaminationResult(est, nh, cr.n_sites)
+    if relatedness and ok.sum() >= 2:
+        usable_loci = np.zeros(manifest.num_loci, np.bool_)
+        usable_loci[auto] = True
+        g_ok = np.flatnonzero(ok)
+        sites = select_sites(genotype_counts(store.gt, None if ok.all() else g_ok), usable_loci)
+        G = SiteView(store.gt, None if ok.all() else g_ok)[sites]
+        ps = pair_stats(G)
+        ped = pedigree or Pedigree({}, {})
+        sub = [ids[i] for i in g_ok]
+        for a in range(len(sub)):
+            for b in range(a + 1, len(sub)):
+                phi, conc = float(ps.kinship[a, b]), float(ps.concordance[a, b])
+                exp = ped.expected(sub[a], sub[b])
+                dup = conc >= DUPLICATE_CONCORDANCE
+                if dup:
+                    res.duplicates.setdefault(sub[a], []).append(sub[b])
+                    res.duplicates.setdefault(sub[b], []).append(sub[a])
+                rel = math.isfinite(phi) and phi > DEGREE_CUTS[2][0]  # second degree or closer
+                if rel and not exp:
+                    res.unexpected.update((sub[a], sub[b]))
+                if all_pairs or dup or exp or (math.isfinite(phi) and phi > DEGREE_CUTS[3][0]):
+                    res.pairs.append(
+                        {
+                            "sample_a": sub[a],
+                            "sample_b": sub[b],
+                            "n_sites": int(ps.n_sites[a, b]),
+                            "kinship": phi,
+                            "concordance": conc,
+                            "degree": degree(phi) if math.isfinite(phi) else "",
+                            "expected": exp,
+                        }
+                    )
+        for child, parents in (ped.parents or {}).items():
+            if child not in ids:
+                continue
+            ci = ids.index(child)
+            prow = [np.asarray(store.gt[ids.index(p)])[auto] for p in parents[:2]]
+            e, n = mendel_errors(
+                np.asarray(store.gt[ci])[auto], prow[0] if prow else None, prow[1] if len(prow) > 1 else None
+            )
+            res.mendel.append(
+                {
+                    "child": child,
+                    "parent1": parents[0] if parents else "",
+                    "parent2": parents[1] if len(parents) > 1 else "",
+                    "n_errors": e,
+                    "n_informative": n,
+                    "error_rate": e / n if n else math.nan,
+                }
+            )
+    return res
+
+
+class _NullPool:
+    def __enter__(self) -> None:
+        return None
+
+    def __exit__(self, *a: object) -> None:
+        return None
+
+
+def _pool(workers: int) -> Any:
+    if workers <= 1:
+        return _NullPool()
+    return ProcessPoolExecutor(workers, mp_context=mp.get_context("spawn"))
+
+
+def default_tmp(out: Path) -> Path:
+    base = Path(os.environ.get("SIMURG_TMP", out))
+    base.mkdir(parents=True, exist_ok=True)
+    return Path(tempfile.mkdtemp(prefix=".simurg_cohort_", dir=base))
+
+
+def contamination_rows(ids: Sequence[str], cr: ContaminationResult | None) -> dict[str, dict[str, Any]]:
+    if cr is None:
+        return {}
+    return {sid: {"contamination": float(cr.estimate[i]), "n_hom_sites": int(cr.n_hom[i])} for i, sid in enumerate(ids)}
+
+
+def status_from_thresholds(th: Thresholds) -> list[dict[str, Any]]:
+    return th.as_records()
+
+
+# ================================================================================================ locus summary
+def write_locus_summary(
+    path: Path, store: CohortStore, manifest: Manifest, egt: Any, samples: npt.NDArray[np.intp] | None = None
+) -> None:
+    """OUT-072, streamed in locus blocks. ``egt`` is a parsed cluster file (records joined by locus name)."""
+    from simurg.core.atomic import atomic_open
+    from simurg.qc.locus_summary import COLUMNS, summarize_block
+
+    assert store.score is not None
+    idx = egt.name_index()
+    rec = np.array([idx.get(nm, -1) for nm in manifest.names], dtype=np.int64)
+    has = rec >= 0
+    safe = np.where(has, rec, 0)
+
+    def take(a: npt.NDArray[Any]) -> npt.NDArray[np.float64]:
+        v = np.asarray(a)[safe].astype(np.float64)
+        v[~has] = np.nan
+        return v
+
+    cl = {
+        "locus_id": np.where(has, np.asarray(egt.addresses)[safe], 0),
+        "gentrain_score": take(egt.gentrain_score),
+        "cluster_sep": take(egt.cluster_sep),
+    }
+    for j, gname in enumerate(("aa", "ab", "bb")):
+        cl[f"{gname}_t_mean"] = take(egt.theta_mean[:, j])
+        cl[f"{gname}_t_sd"] = take(egt.theta_dev[:, j])
+        cl[f"{gname}_r_mean"] = take(egt.r_mean[:, j])
+        cl[f"{gname}_r_sd"] = take(egt.r_dev[:, j])
+    strands = [s if s in ("+", "-") else "U" for s in (manifest.ref_strands or [""] * manifest.num_loci)]
+
+    def f(v: Any) -> str:
+        if isinstance(v, float | np.floating):
+            return "NaN" if math.isnan(v) else f"{float(v):.8g}"
+        return str(v)
+
+    with atomic_open(path, "w", encoding="utf-8") as raw:
+        fh = cast(IO[str], raw)
+        fh.write("\t".join(COLUMNS) + "\n")
+        for sl in _blocks(manifest.num_loci, 65536):
+            g = np.asarray(store.gt[:, sl])
+            s = np.asarray(store.score[:, sl])
+            if samples is not None:
+                g, s = g[samples], s[samples]
+            stats = summarize_block(g, s)
+            for k in range(sl.stop - sl.start):
+                i = sl.start + k
+                row = {"locus_name": manifest.names[i], "plus_minus_strand": strands[i]}
+                vals = [row[c] if c in row else (cl[c][i] if c in cl else stats[c][k]) for c in COLUMNS]
+                fh.write("\t".join(f(v) for v in vals) + "\n")
+
+
+# ================================================================================================ GTC-based commands
+SAMPLE_QC_COLUMNS = [
+    "sample_id",
+    "sample_name",
+    "status",
+    "call_rate",
+    "autosomal_call_rate",
+    "logr_dev",
+    "sex_est",
+    "sex_reported",
+    "sex_check",
+    "het_rate",
+    "het_z",
+    "contamination",
+    "n_hom_sites",
+    "functional_status",
+    "flags",
+]
+RELATEDNESS_COLUMNS = ["sample_a", "sample_b", "n_sites", "kinship", "concordance", "degree", "expected"]
+MENDEL_COLUMNS = ["child", "parent1", "parent2", "n_errors", "n_informative", "error_rate"]
+
+
+@dataclass
+class QcCohortOptions:
+    """Options shared by ``qc samples | contamination | relatedness | loci`` (CLI-025)."""
+
+    command: str
+    out: Path
+    bpm: Path
+    gtcs: Path | None = None
+    sheet: Path | None = None
+    egt: Path | None = None
+    threads: int = 0
+    force: bool = False
+    all_pairs: bool = False
+    pedigree_from_sheet: bool = True
+    median: str = "reference"
+    thresholds: Thresholds | None = None
+    argv: list[str] = field(default_factory=list)
+
+
+def _gtc_samples(opts: QcCohortOptions, rep: Reporter) -> tuple[list[tuple[str, str]], SampleSheet | None]:
+    """(sample_id, path) in sheet order (sheet given) or discovery order (IN-030, IN-016)."""
+    from simurg.formats.discovery import discover_gtcs
+    from simurg.formats.samplesheet import parse_samplesheet
+
+    sheet = parse_samplesheet(opts.sheet) if opts.sheet else None
+    root = opts.gtcs or (sheet.root_folder if sheet else None)
+    if root is None:
+        raise SimurgError("E_NO_INPUT", "give --gtcs or a sample sheet with a RootFolder")
+    found = discover_gtcs(root)
+    if sheet is None:
+        return [(k, str(v)) for k, v in found.items()], None
+    out = []
+    for r in sheet.rows:
+        if r.sample_id in found:
+            out.append((r.sample_id, str(found[r.sample_id])))
+        else:
+            rep.warn("W_QC_UNMATCHED_SAMPLES", f"{r.sample_id}: no GTC under {root}", sample_id=r.sample_id)
+    if not out:
+        raise SimurgError("E_NO_INPUT", f"none of the sheet's samples has a GTC under {root}")
+    return out, sheet
+
+
+def run_qc_cohort(opts: QcCohortOptions, rep: Reporter) -> int:
+    """``simurg qc samples | contamination | relatedness | loci`` (CLI-025)."""
+    from simurg.formats.bpm import read_bpm
+    from simurg.pipeline.inputs import check_paths
+    from simurg.qc.status import functional_status, sample_flags
+    from simurg.qc.thresholds import resolve_thresholds
+
+    rm = RunManifest(command=f"qc {opts.command}", argv=opts.argv)
+    rep.run_id = rm.run_id
+    check_paths([(opts.bpm, "BPM"), (opts.gtcs, "GTC folder"), (opts.sheet, "sample sheet"), (opts.egt, "EGT")])
+    if opts.command == "loci" and opts.egt is None:
+        raise SimurgError("E_NO_INPUT", "qc loci needs --egt (cluster statistics, OUT-072)")
+    names = {
+        "samples": ["sample_qc.tsv", "contamination.tsv", "relatedness.tsv"],
+        "contamination": ["contamination.tsv"],
+        "relatedness": ["relatedness.tsv"],
+        "loci": ["locus_summary.tsv"],
+    }[opts.command]
+    outputs = [opts.out / n for n in names]
+    check_outputs_absent(outputs, force=opts.force)
+    samples, sheet = _gtc_samples(opts, rep)
+    rep.info(f"reading manifest {opts.bpm}")
+    manifest = read_bpm(opts.bpm)
+    workers = max(1, opts.threads or (os.cpu_count() or 1))
+    opts.out.mkdir(parents=True, exist_ok=True)
+    tmp = default_tmp(opts.out)
+    want_baf = opts.command in ("samples", "contamination")
+    rep.info(f"loading {len(samples)} GTC(s) with {workers} worker(s)")
+    store, infos = load_cohort(
+        samples, manifest.num_loci, tmp, baf=want_baf, scores=opts.command == "loci", workers=workers
+    )
+    try:
+        for s in infos:
+            if s.error:
+                rep.error(s.error.split(":", 1)[0], f"{s.sample_id}: {s.error}", sample_id=s.sample_id)
+        ids = [s.sample_id for s in infos]
+        if opts.command == "loci":
+            egt_path = opts.egt
+            assert egt_path is not None
+            from simurg.formats.egt import read_egt
+
+            ok = np.flatnonzero([not s.error for s in infos])
+            write_locus_summary(opts.out / "locus_summary.tsv", store, manifest, read_egt(egt_path), ok)
+        else:
+            ped = pedigree_from_sheet(sheet, ids, rep) if opts.pedigree_from_sheet else None
+            res = cohort_qc(
+                store,
+                infos,
+                manifest,
+                rep,
+                pedigree=ped,
+                contamination=opts.command in ("samples", "contamination"),
+                relatedness=opts.command in ("samples", "relatedness"),
+                all_pairs=opts.all_pairs,
+                median=opts.median,
+                workers=workers,
+            )
+            crow = contamination_rows(ids, res.contamination)
+            if opts.command in ("samples", "contamination"):
+                write_tsv(
+                    opts.out / "contamination.tsv",
+                    ["sample_id", "contamination", "n_hom_sites"],
+                    [{"sample_id": sid, **crow.get(sid, {})} for sid in ids],
+                )
+            if opts.command in ("samples", "relatedness"):
+                write_tsv(opts.out / "relatedness.tsv", RELATEDNESS_COLUMNS, res.pairs)
+                if res.mendel:
+                    write_tsv(opts.out / "mendel.tsv", MENDEL_COLUMNS, res.mendel)
+                    outputs.append(opts.out / "mendel.tsv")
+            if opts.command == "samples":
+                th = opts.thresholds or resolve_thresholds()
+                auto = autosome_mask(manifest.chroms) & ~manifest.intensity_only
+                auto_idx = np.flatnonzero(auto)
+                sex_rep = {r.sample_id: r.reported_sex for r in sheet.rows} if sheet else {}
+                rows = []
+                for i, s in enumerate(infos):
+                    g = np.asarray(store.gt[i])[auto_idx]
+                    row: dict[str, Any] = {
+                        "sample_id": s.sample_id,
+                        "sample_name": s.sample_name,
+                        "status": "FAILED" if s.error else "OK",
+                        "call_rate": s.call_rate,
+                        "autosomal_call_rate": float((g > 0).mean()) if len(g) and not s.error else math.nan,
+                        "logr_dev": s.logr_dev,
+                        "sex_est": s.sex_est,
+                        "sex_reported": sex_rep.get(s.sample_id, "U"),
+                        "het_rate": float(res.het_rate[i]),
+                        "het_z": res.het_z[i],
+                        **crow.get(s.sample_id, {}),
+                    }
+                    known = {"M", "F"}
+                    row["sex_check"] = (
+                        ("match" if row["sex_est"] == row["sex_reported"] else "mismatch")
+                        if row["sex_est"] in known and row["sex_reported"] in known
+                        else "unknown"
+                    )
+                    row["functional_status"] = functional_status(row, th)
+                    row["flags"] = ";".join(
+                        sample_flags(
+                            row,
+                            th,
+                            het_z=res.het_z[i],
+                            duplicates=res.duplicates.get(s.sample_id, []),
+                            unexpected_relatedness=s.sample_id in res.unexpected,
+                        )
+                    )
+                    rows.append(row)
+                write_tsv(opts.out / "sample_qc.tsv", SAMPLE_QC_COLUMNS, rows)
+    finally:
+        store.close()
+    for p in outputs:
+        if p.exists():
+            rm.add_output(p, f"simurg.{p.stem}/1")
+    rm.summary = {"samples": len(infos), "failed": sum(bool(s.error) for s in infos)}
+    rm.status = "completed"
+    rm.write(opts.out / f"run_manifest.qc_{opts.command}.json")
+    failed = sum(bool(s.error) for s in infos)
+    rep.info(f"done: qc {opts.command} on {len(infos) - failed}/{len(infos)} sample(s); outputs in {opts.out}")
+    if failed == 0:
+        return ExitCode.OK
+    return ExitCode.PARTIAL if failed < len(infos) else ExitCode.ALL_FAILED
+
+
+# ================================================================================================ suggest-thresholds
+def _read_table(path: Path) -> list[dict[str, str]]:
+    import csv
+
+    text = path.read_text(encoding="utf-8-sig")
+    delim = "\t" if "\t" in text.split("\n", 1)[0] else ","
+    return list(csv.DictReader(text.splitlines(), delimiter=delim))
+
+
+def dataset_values(folder: Path) -> dict[str, list[float]]:
+    """Per-key metric values of one QC dataset folder (QC-040 file contract)."""
+    from simurg.qc.thresholds import KEY_TO_COLUMN
+
+    def num(v: str | None) -> float:
+        try:
+            return float(v) if v not in (None, "") else math.nan
+        except ValueError:
+            return math.nan
+
+    out: dict[str, list[float]] = {}
+    qc = folder / CONTROLS_QC
+    if not qc.exists():
+        raise SimurgError(
+            "E_QC_DATASET_INCOMPLETE", f"Dataset {folder.name} is missing required file(s): {CONTROLS_QC}"
+        )
+    rows = _read_table(qc)
+    for key, col in KEY_TO_COLUMN.items():
+        if rows and col in rows[0]:
+            out[key] = [num(r.get(col)) for r in rows]
+    summary = folder / "sample_summary.tsv"
+    if summary.exists():
+        srows = [r for r in _read_table(summary) if r.get("status", "OK") == "OK"]
+        out["callRate"] = [num(r.get("autosomal_call_rate")) for r in srows]
+        out["logRDev"] = [num(r.get("logr_dev")) for r in srows]
+    return out
+
+
+def run_suggest_thresholds(
+    data: Sequence[Path], out_yaml: Path, th: Thresholds, rep: Reporter, *, k: float = 4.0, min_samples: int = 30
+) -> int:
+    from simurg.qc.suggest import render_yaml, suggest
+
+    if not data:
+        raise SimurgError("E_QC_NO_DATASETS", "No valid dataset folders detected.")
+    merged: dict[str, list[float]] = {}
+    for d in data:
+        if not d.exists():
+            raise SimurgError("E_QC_NO_DATASETS", f"Input data folder does not exist: {d}", path=str(d))
+        for key, vals in dataset_values(d).items():
+            merged.setdefault(key, []).extend(vals)
+    sugs = suggest(merged, th, k=k, min_samples=min_samples)
+    for s in sugs:
+        if s.note:
+            rep.warn("W_QC_TOO_FEW", f"{s.key}: {s.note}; left unset")
+    out_yaml.parent.mkdir(parents=True, exist_ok=True)
+    atomic_write_text(out_yaml, render_yaml(sugs, k))
+    write_tsv(
+        out_yaml.with_suffix(".tsv"),
+        ["key", "n", "centre", "spread", "value", "pct_flagged_current", "pct_flagged_suggested", "note"],
+        [s.__dict__ for s in sugs],
+    )
+    rep.info(f"suggested thresholds written to {out_yaml}")
+    return ExitCode.OK

src/simurg/qc/contamination.py

@@ -0,0 +1,179 @@
+"""Contamination estimate by B-allele-frequency regression (spec QC-064, plan §9.3).
+
+[LIT: Jun G. et al. 2012, AJHG 91:839 — BAF regression on homozygous calls]
+[OSS: gtc2vcf ``BAFregress.c`` (MIT) — reference numerics: per-site median adjustment of AA/BB BAF when a cluster
+has >= 5 calls, float32 accumulators updated with double-precision addends, OLS slope in double.]
+
+For every site with cohort B-allele frequency p in [t, 1 - t] (t = 0.1f):
+  AA call: x = p,     y = BAF - median_AA(BAF)
+  BB call: x = 1 - p, y = 1 - (BAF - (median_BB(BAF) - 1))
+The contamination fraction is the slope m = (n·Σxy − Σx·Σy) / (n·Σx² − (Σx)²).
+
+QUIRK (K-017): for an even number of calls the reference averages the quickselect result with whatever element sits
+just below it after the partial partition (klib ``ks_ksmall``), which is not always the lower middle value.
+``median="reference"`` (default) reproduces it exactly; ``median="exact"`` uses the true median.
+"""
+
+from __future__ import annotations
+
+import math
+from collections.abc import Callable, Iterable
+from dataclasses import dataclass
+from typing import Any
+
+import numpy as np
+import numpy.typing as npt
+
+F32 = np.float32
+CHUNK = 4096
+
+
+@dataclass
+class ContaminationResult:
+    estimate: npt.NDArray[np.float64]  # per sample; NaN when undefined
+    n_hom: npt.NDArray[np.int64]
+    n_sites: int
+
+
+def b_allele_frequency(genotypes: npt.NDArray[np.uint8]) -> npt.NDArray[np.float64]:
+    """Per-site B-allele frequency over called samples (rows = sites); NaN where nobody is called."""
+    n_ab = (genotypes == 2).sum(axis=1)
+    n_bb = (genotypes == 3).sum(axis=1)
+    an = 2 * ((genotypes == 1).sum(axis=1) + n_ab + n_bb)
+    with np.errstate(invalid="ignore", divide="ignore"):
+        return np.where(an > 0, (n_ab + 2 * n_bb) / an.astype(np.float64), np.nan)
+
+
+def ks_ksmall(arr: list[float], kk: int) -> float:
+    """k-th smallest by in-place quickselect, partitioning exactly as klib ``ksort.h`` (MIT; htslib copy)
+    so that the element left at ``arr[kk-1]`` matches the reference."""
+    low, high, k = 0, len(arr) - 1, kk
+    while True:
+        if high <= low:
+            return arr[k]
+        if high == low + 1:
+            if arr[high] < arr[low]:
+                arr[low], arr[high] = arr[high], arr[low]
+            return arr[k]
+        mid = low + (high - low) // 2
+        if arr[high] < arr[mid]:
+            arr[mid], arr[high] = arr[high], arr[mid]
+        if arr[high] < arr[low]:
+            arr[low], arr[high] = arr[high], arr[low]
+        if arr[low] < arr[mid]:
+            arr[mid], arr[low] = arr[low], arr[mid]
+        arr[mid], arr[low + 1] = arr[low + 1], arr[mid]
+        ll, hh = low + 1, high
+        pivot = arr[low]
+        while True:
+            ll += 1
+            while arr[ll] < pivot:
+                ll += 1
+            hh -= 1
+            while pivot < arr[hh]:
+                hh -= 1
+            if hh < ll:
+                break
+            arr[ll], arr[hh] = arr[hh], arr[ll]
+        arr[low], arr[hh] = arr[hh], arr[low]
+        if hh <= k:
+            low = ll
+        if hh >= k:
+            high = hh - 1
+
+
+def _median_reference(values: npt.NDArray[np.float32]) -> np.float32:
+    t = values.tolist()  # float32 → exact Python floats; comparisons are unchanged
+    j = len(t)
+    m = F32(ks_ksmall(t, j // 2))
+    if j % 2 == 0:
+        m = F32((m + F32(t[j // 2 - 1])) * F32(0.5))
+    return m
+
+
+def _median_exact(baf: npt.NDArray[np.float32], mask: npt.NDArray[np.bool_]) -> npt.NDArray[np.float32]:
+    """Row-wise true median of ``baf`` where ``mask``; float32, even counts average the middle pair as (a+b)*0.5f."""
+    v = np.where(mask, baf, F32(np.inf)).astype(F32)
+    v.sort(axis=1)
+    k = mask.sum(axis=1)
+    rows = np.arange(len(v))
+    hi = v[rows, k // 2]
+    lo = v[rows, np.maximum(k // 2 - 1, 0)]
+    return np.where(k % 2 == 0, ((lo + hi) * F32(0.5)).astype(F32), hi).astype(F32)
+
+
+def cluster_adjustments(
+    task: tuple[npt.NDArray[np.uint8], npt.NDArray[np.float32], int, str],
+) -> tuple[npt.NDArray[np.float32], npt.NDArray[np.float32]]:
+    """Per-site AA and BB BAF adjustments for one chunk of sites (picklable for process pools)."""
+    g, b, adjust_min, median = task
+    adj_aa = np.zeros(len(g), F32)
+    adj_bb = np.zeros(len(g), F32)
+    if adjust_min == -1:
+        return adj_aa, adj_bb
+    for code, out, shift in ((1, adj_aa, F32(0.0)), (3, adj_bb, F32(1.0))):
+        mask = g == code
+        use = np.flatnonzero(mask.sum(axis=1) >= adjust_min)
+        if not len(use):
+            continue
+        if median == "exact":
+            out[use] = _median_exact(b[use], mask[use]) - shift
+        else:
+            for i in use:
+                out[i] = _median_reference(b[i][mask[i]]) - shift
+    return adj_aa, adj_bb
+
+
+MapFn = Callable[[Callable[[Any], Any], Iterable[Any]], Iterable[Any]]
+
+
+def baf_regress(
+    genotypes: npt.NDArray[np.uint8],
+    baf: npt.NDArray[np.float32],
+    af: npt.NDArray[np.float64] | None = None,
+    *,
+    threshold: float = 0.1,
+    adjust_min: int = 5,
+    median: str = "reference",
+    map_fn: MapFn = map,
+) -> ContaminationResult:
+    """Contamination per sample. ``genotypes`` (0 no call, 1 AA, 2 AB, 3 BB) and ``baf`` are (sites × samples);
+    ``af`` defaults to the cohort B-allele frequency. Sites are used in row order (the reference's float32
+    accumulation is order-dependent). ``map_fn`` (e.g. an executor's ``map``) parallelises the per-site medians."""
+    if median not in ("reference", "exact"):
+        raise ValueError(f"median must be 'reference' or 'exact', got {median!r}")
+    n = genotypes.shape[1]
+    if af is None:
+        af = b_allele_frequency(genotypes)
+    t = float(F32(threshold))  # the reference stores the threshold as float: 0.1f = 0.10000000149…
+    keep = np.flatnonzero(np.isfinite(af) & (af >= t) & (af <= 1.0 - t))
+    chunks = [keep[s : s + CHUNK] for s in range(0, len(keep), CHUNK)]
+    tasks = ((np.asarray(genotypes[r]), np.asarray(baf[r], dtype=F32), adjust_min, median) for r in chunks)
+    sx2, sxy, sx, sy = (np.zeros(n, F32) for _ in range(4))
+    cnt = np.zeros(n, np.int64)
+    for rows, (adj_aa, adj_bb) in zip(chunks, map_fn(cluster_adjustments, tasks), strict=True):
+        g = np.asarray(genotypes[rows])
+        b = np.asarray(baf[rows], dtype=F32)
+        p = af[rows]
+        is_aa, is_bb = g == 1, g == 3
+        y_aa = (b - adj_aa[:, None]).astype(F32).astype(np.float64)
+        y_bb = (b - adj_bb[:, None]).astype(F32).astype(np.float64)
+        for i in range(len(rows)):
+            a, m = is_aa[i], is_aa[i] | is_bb[i]
+            if not m.any():
+                continue
+            pa, pb = p[i], 1.0 - p[i]
+            x = np.where(a, pa, pb)[m]
+            xy = np.where(a, pa * y_aa[i], pb * (1.0 - y_bb[i]))[m]
+            y = np.where(a, y_aa[i], 1.0 - y_bb[i])[m]
+            sx2[m] = (sx2[m].astype(np.float64) + x * x).astype(F32)
+            sxy[m] = (sxy[m].astype(np.float64) + xy).astype(F32)
+            sx[m] = (sx[m].astype(np.float64) + x).astype(F32)
+            sy[m] = (sy[m].astype(np.float64) + y).astype(F32)
+            cnt[m] += 1
+    nf = cnt.astype(np.float64)
+    fx2, fxy, fx, fy = (s.astype(np.float64) for s in (sx2, sxy, sx, sy))
+    denom = nf * fx2 - fx * fx
+    with np.errstate(invalid="ignore", divide="ignore"):
+        est = np.where(denom != 0, (nf * fxy - fx * fy) / denom, math.nan)
+    return ContaminationResult(estimate=est, n_hom=cnt, n_sites=len(keep))

src/simurg/qc/control_names.yaml

@@ -0,0 +1,84 @@
+# Canonical control-probe IDs (spec QC-002). Keys are *normalised* names: case-folded, whitespace collapsed,
+# no spaces next to or inside parentheses, "bgnd" written as "bkg". Several manifest spellings map to one ID.
+canonical:
+  # staining
+  biotin(high): biotin_high
+  biotin(bkg): biotin_bkg
+  biotin(5k): biotin_5k
+  dnp(high): dnp_high
+  dnp(bkg): dnp_bkg
+  dnp(20k): dnp_20k
+  # extension
+  extension(a): ext_A
+  extension(t): ext_T
+  extension(c): ext_C
+  extension(g): ext_G
+  # hybridization
+  hyb(high): hyb_high
+  hyb(medium): hyb_medium
+  hyb(low): hyb_low
+  # target removal
+  target removal: target_removal_1
+  target removal(1): target_removal_1
+  target removal 1: target_removal_1
+  target removal(2): target_removal_2
+  target removal 2: target_removal_2
+  # stringency
+  string(pm): string_pm
+  string(mm): string_mm
+  # non-specific binding (one bead type per colour label on genotyping chips)
+  nsb(bkg): nsb_bkg
+  nsb(bkg) red: nsb_bkg
+  nsb(bkg) green: nsb_bkg
+  nsb(bkg) blue: nsb_bkg
+  nsb(bkg) purple: nsb_bkg
+  # non-polymorphic
+  np(a): np_A
+  np(t): np_T
+  np(c): np_C
+  np(g): np_G
+  # restoration (FFPE)
+  restoration: restoration
+  restore: restoration
+  # methylation (metrics in v2; names mapped already so controls_raw.tsv is labelled consistently)
+  bs conversion i c1: bs1_C1
+  bs conversion i c2: bs1_C2
+  bs conversion i c3: bs1_C3
+  bs conversion i c4: bs1_C4
+  bs conversion i c5: bs1_C5
+  bs conversion i u1: bs1_U1
+  bs conversion i u2: bs1_U2
+  bs conversion i u3: bs1_U3
+  bs conversion i u4: bs1_U4
+  bs conversion i u5: bs1_U5
+  bs conversion ii-1: bs2_1
+  bs conversion ii-2: bs2_2
+  bs conversion ii-3: bs2_3
+  bs conversion ii-4: bs2_4
+  gt mismatch 1(pm): gtmm_1_pm
+  gt mismatch 2(pm): gtmm_2_pm
+  gt mismatch 3(pm): gtmm_3_pm
+  gt mismatch 4(pm): gtmm_4_pm
+  gt mismatch 5(pm): gtmm_5_pm
+  gt mismatch 6(pm): gtmm_6_pm
+  gt mismatch 1(mm): gtmm_1_mm
+  gt mismatch 2(mm): gtmm_2_mm
+  gt mismatch 3(mm): gtmm_3_mm
+  gt mismatch 4(mm): gtmm_4_mm
+  gt mismatch 5(mm): gtmm_5_mm
+  gt mismatch 6(mm): gtmm_6_mm
+  specificity 1: spec2_1
+  specificity 2: spec2_2
+  specificity 3: spec2_3
+
+# Categories treated as methylation-only (dataset type detection, QC-043).
+methylation_categories:
+  - bisulfite conversion i
+  - bisulfite conversion ii
+  - specificity i
+  - specificity ii
+
+# Report exclusion list (QC-010): raw-probe views leave these out; derived metrics still use them.
+report_exclude:
+  categories: [negative, norm]
+  names: ["np(g) 1", "np(g) 2", "np(g) 3", "np(g) 4", "np(g) 5", "biotin(5k)", "dnp(20k)"]

src/simurg/qc/controls.py

@@ -0,0 +1,257 @@
+"""Control-probe QC metrics (spec QC-001 … QC-011; plan §9.1).
+
+Controls are the non-polymorphic bead types listed in the manifest's control section. Each gives a raw green and
+red intensity (IDAT means, no normalization). Derived metrics are signal/background ratios per assay step.
+[DOC: Illumina BeadArray Controls Reporter software guide and Infinium controls training material — formulas only]
+"""
+
+from __future__ import annotations
+
+import math
+import re
+from collections.abc import Sequence
+from dataclasses import dataclass, field
+from functools import cache
+from importlib import resources
+from typing import Any
+
+import numpy as np
+import numpy.typing as npt
+import yaml
+
+from simurg.formats.idat import Idat
+from simurg.signal.raw_xy import AddressIndex
+
+DEFAULT_BG_OFFSET = 3000.0
+
+
+@dataclass(frozen=True)
+class ControlBead:
+    """One control bead type: a (manifest control line, distinct address) pair, in manifest order."""
+
+    category: str
+    color: str
+    name: str
+    address: int
+    canonical: str  # canonical ID, "" when the name is not in control_names.yaml
+
+
+@dataclass(frozen=True)
+class Metric:
+    column: str  # controls_qc.tsv column
+    key: str  # threshold YAML key (QC-020)
+    genotyping: bool
+    methylation: bool
+
+
+# QC-004, in table order.
+METRICS: tuple[Metric, ...] = (
+    Metric("restoration", "restoration", True, True),
+    Metric("staining_green", "stainingGreen", True, True),
+    Metric("staining_red", "stainingRed", True, True),
+    Metric("extension_green", "extensionGreen", True, True),
+    Metric("extension_red", "extensionRed", True, True),
+    Metric("hyb_high_med", "hybGreenHighMed", True, True),
+    Metric("hyb_med_low", "hybGreenMediumLow", True, True),
+    Metric("target_removal_1", "targetRemovalI", True, True),
+    Metric("target_removal_2", "targetRemovalII", False, True),
+    Metric("stringency", "stringency", True, False),
+    Metric("nsb_green", "nonSpecificBindingGreen", True, False),
+    Metric("nsb_red", "nonSpecificBindingRed", True, False),
+    Metric("np_green", "npGreen", True, True),
+    Metric("np_red", "npRed", True, True),
+    Metric("specificity_1_green", "specificityIGreen", False, True),
+    Metric("specificity_1_red", "specificityIRed", False, True),
+    Metric("specificity_2", "specificityII", False, True),
+    Metric("specificity_2_bg", "specificityIIBg", False, True),
+    Metric("bs_conv_1_green", "bisulfiteConversionIGreen", False, True),
+    Metric("bs_conv_1_green_bg", "bisulfiteConversionIGreenBg", False, True),
+    Metric("bs_conv_1_red", "bisulfiteConversionIRed", False, True),
+    Metric("bs_conv_1_red_bg", "bisulfiteConversionIRedBg", False, True),
+    Metric("bs_conv_2", "bisulfiteConversionII", False, True),
+    Metric("bs_conv_2_bg", "bisulfiteConversionIIBg", False, True),
+)
+ASSAY_CONTROL_COLUMN = "assay_control_norm_r"  # QC-005: copied from the sample summary
+COMPAT_METRIC_NAMES = {
+    "restoration": "RestorationQC",
+    "staining_green": "StainingGreenQC",
+    "staining_red": "StainingRedQC",
+    "extension_green": "ExtensionGreenQC",
+    "extension_red": "ExtensionRedQC",
+    "hyb_high_med": "HybridizationHighMediumQC",
+    "hyb_med_low": "HybridizationMediumLowQC",
+    "target_removal_1": "TargetRemovalQC",
+    "target_removal_2": "TargetRemoval2QC",
+    "stringency": "StringencyQC",
+    "nsb_green": "NonSpecificBindingGreenQC",
+    "nsb_red": "NonSpecificBindingRedQC",
+    "np_green": "NonPolymorphicGreenQC",
+    "np_red": "NonPolymorphicRedQC",
+    ASSAY_CONTROL_COLUMN: "AssayControlQC",
+}
+
+
+def metrics_for(array_type: str) -> list[Metric]:
+    """QC-005 / QC-006 metric sets."""
+    if array_type == "genotyping":
+        return [m for m in METRICS if m.genotyping]
+    return [m for m in METRICS if m.methylation]
+
+
+# ------------------------------------------------------------------------------------------------ names
+@cache
+def name_config() -> dict[str, Any]:
+    text = resources.files("simurg.qc").joinpath("control_names.yaml").read_text(encoding="utf-8")
+    cfg: dict[str, Any] = yaml.safe_load(text)
+    return cfg
+
+
+def normalise_name(name: str) -> str:
+    """QC-002 normalisation: case-fold, collapse whitespace, no spaces before or inside parentheses, bgnd = bkg."""
+    s = " ".join(name.casefold().split())
+    s = re.sub(r"\s*\(\s*", "(", s)
+    s = re.sub(r"\s*\)", ")", s)
+    return re.sub(r"\bbgnd\b|\bbkgd\b", "bkg", s)
+
+
+def canonical_id(name: str) -> str:
+    return str(name_config()["canonical"].get(normalise_name(name), ""))
+
+
+def excluded_from_report(bead: ControlBead) -> bool:
+    """QC-010."""
+    ex = name_config()["report_exclude"]
+    return bead.category.casefold() in ex["categories"] or normalise_name(bead.name) in ex["names"]
+
+
+def is_methylation_controls(beads: Sequence[ControlBead]) -> bool:
+    cats = set(name_config()["methylation_categories"])
+    return any(b.category.casefold() in cats for b in beads)
+
+
+def beads_from_definitions(defs: Sequence[dict[str, str]]) -> list[ControlBead]:
+    """Manifest control lines (``address`` may be ``a:b:…``) → bead types in manifest order (QC-001)."""
+    out: list[ControlBead] = []
+    for d in defs:
+        seen: set[int] = set()
+        for part in str(d["address"]).split(":"):
+            part = part.strip()
+            if not part:
+                continue
+            addr = int(part)
+            if addr in seen:
+                continue
+            seen.add(addr)
+            out.append(ControlBead(d["category"], d["color"], d["name"], addr, canonical_id(d["name"])))
+    return out
+
+
+def compare_definitions(a: Sequence[ControlBead], b: Sequence[ControlBead]) -> list[str]:
+    """Differences between two control lists (BPM vs CSV, QC-001); empty when they agree."""
+    ka = {(x.address, normalise_name(x.name)) for x in a}
+    kb = {(x.address, normalise_name(x.name)) for x in b}
+    diffs = [f"only in BPM: {addr} {nm}" for addr, nm in sorted(ka - kb)[:5]]
+    diffs += [f"only in CSV: {addr} {nm}" for addr, nm in sorted(kb - ka)[:5]]
+    return diffs
+
+
+# ------------------------------------------------------------------------------------------------ intensities
+def control_intensities(
+    grn: Idat, red: Idat, beads: Sequence[ControlBead]
+) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]:
+    """Raw green/red means per bead type; NaN where the IDAT lacks the address."""
+    addrs = np.array([b.address for b in beads], dtype=np.int64)
+    out = []
+    for idat in (grn, red):
+        rows, found = AddressIndex(idat).lookup(addrs)
+        v = idat.mean[rows].astype(np.float64)
+        v[~found] = math.nan
+        out.append(v)
+    return out[0], out[1]
+
+
+# ------------------------------------------------------------------------------------------------ metrics
+@dataclass
+class ControlQC:
+    values: dict[str, float]  # column -> value (NaN when not evaluable or degenerate)
+    missing: dict[str, list[str]] = field(default_factory=dict)  # column -> canonical IDs absent from the manifest
+    zero_divisor: list[str] = field(default_factory=list)  # columns with a zero denominator (QC-008)
+
+    @property
+    def na_metrics(self) -> list[str]:
+        return list(self.missing)
+
+
+class _Probes:
+    def __init__(self, beads: Sequence[ControlBead], green: npt.NDArray[np.float64], red: npt.NDArray[np.float64]):
+        self.idx: dict[str, list[int]] = {}
+        for i, b in enumerate(beads):
+            if b.canonical:
+                self.idx.setdefault(b.canonical, []).append(i)
+        self.ch = {"G": green, "R": red}
+        self.absent: set[str] = set()
+
+    def vals(self, ids: Sequence[str], ch: str) -> npt.NDArray[np.float64]:
+        rows: list[int] = []
+        for cid in ids:
+            if cid not in self.idx:
+                self.absent.add(cid)
+            rows.extend(self.idx.get(cid, []))
+        return self.ch[ch][rows]
+
+    def mean(self, cid: str, ch: str) -> float:
+        v = self.vals([cid], ch)
+        return float(v.mean()) if len(v) else math.nan
+
+    def min(self, ids: Sequence[str], ch: str) -> float:
+        v = self.vals(ids, ch)
+        return float(v.min()) if len(v) else math.nan
+
+    def max(self, ids: Sequence[str], ch: str) -> float:
+        v = self.vals(ids, ch)
+        return float(v.max()) if len(v) else math.nan
+
+
+def compute_control_metrics(
+    beads: Sequence[ControlBead],
+    green: npt.NDArray[np.float64],
+    red: npt.NDArray[np.float64],
+    array_type: str = "genotyping",
+    bg_offset: float = DEFAULT_BG_OFFSET,
+) -> ControlQC:
+    """QC-003/QC-004 metrics of one sample. Single-probe values are bead means; min/max run over all beads (QC-002)."""
+    if array_type != "genotyping":
+        from simurg.core.errors import SimurgError
+
+        raise SimurgError("E_QC_METHYLATION", "methylation control metrics are planned for v2 (QC-006)")
+    p = _Probes(beads, green, red)
+    x = bg_offset
+    formulas: dict[str, tuple[Any, Any]] = {
+        # column: (numerator thunk, denominator thunk); each thunk records the probes it needs
+        "restoration": (lambda: p.mean("restoration", "G"), lambda: p.max(["ext_A", "ext_T"], "G") + x),
+        "staining_green": (lambda: p.mean("biotin_high", "G"), lambda: p.mean("biotin_bkg", "G")),
+        "staining_red": (lambda: p.mean("dnp_high", "R"), lambda: p.mean("dnp_bkg", "R")),
+        "extension_green": (lambda: p.min(["ext_C", "ext_G"], "G"), lambda: p.max(["ext_A", "ext_T"], "G")),
+        "extension_red": (lambda: p.min(["ext_A", "ext_T"], "R"), lambda: p.max(["ext_C", "ext_G"], "R")),
+        "hyb_high_med": (lambda: p.mean("hyb_high", "G"), lambda: p.mean("hyb_medium", "G")),
+        "hyb_med_low": (lambda: p.mean("hyb_medium", "G"), lambda: p.mean("hyb_low", "G")),
+        "target_removal_1": (lambda: p.max(["ext_A", "ext_T"], "G") + x, lambda: p.mean("target_removal_1", "G")),
+        "stringency": (lambda: p.mean("string_pm", "R"), lambda: p.mean("string_mm", "R")),
+        "nsb_green": (lambda: p.max(["ext_A", "ext_T"], "G") + x, lambda: p.max(["nsb_bkg"], "G")),
+        "nsb_red": (lambda: p.max(["ext_C", "ext_G"], "R") + x, lambda: p.max(["nsb_bkg"], "R")),
+        "np_green": (lambda: p.min(["np_C", "np_G"], "G"), lambda: p.max(["np_A", "np_T"], "G")),
+        "np_red": (lambda: p.min(["np_A", "np_T"], "R"), lambda: p.max(["np_C", "np_G"], "R")),
+    }
+    res = ControlQC(values={})
+    for m in metrics_for(array_type):
+        p.absent = set()
+        num, den = formulas[m.column][0](), formulas[m.column][1]()
+        if p.absent:
+            res.missing[m.column] = sorted(p.absent)
+            res.values[m.column] = math.nan
+        elif den == 0:
+            res.zero_divisor.append(m.column)
+            res.values[m.column] = math.nan
+        else:
+            res.values[m.column] = num / den  # NaN propagates when a bead's address is absent from the IDAT
+    return res

src/simurg/qc/locus_summary.py

@@ -0,0 +1,135 @@
+"""Per-locus summary statistics across samples (spec OUT-072, CLI-025; plan §7.9).
+
+[OSS: BeadArrayFiles ``examples/locus_summary.py`` (BSD-2) — formulas: call/genotype frequencies, minor frequency,
+k-regularised heterozygote excess with k = (0.013/q)^4, ChiTest_P100 = 2·Φ(±dh; 0, 0.1), GenomeStudio percentile of
+the called GenCall scores.] Vectorised over a block of loci; the caller streams blocks for large cohorts.
+"""
+
+from __future__ import annotations
+
+import math
+
+import numpy as np
+import numpy.typing as npt
+
+COLUMNS = [
+    "locus_name",
+    "locus_id",
+    "n_nocalls",
+    "n_calls",
+    "call_freq",
+    "aa_freq",
+    "ab_freq",
+    "bb_freq",
+    "minor_freq",
+    "gentrain_score",
+    "gc50",
+    "gc10",
+    "het_excess",
+    "chi_test_p100",
+    "cluster_sep",
+    "aa_t_mean",
+    "aa_t_sd",
+    "ab_t_mean",
+    "ab_t_sd",
+    "bb_t_mean",
+    "bb_t_sd",
+    "aa_r_mean",
+    "aa_r_sd",
+    "ab_r_mean",
+    "ab_r_sd",
+    "bb_r_mean",
+    "bb_r_sd",
+    "plus_minus_strand",
+    "het_excess_simple",
+]
+HET_EXCESS_EPS = 1e-9  # OQ-030: denominator guard of the simple form
+M_SQRT1_2 = 0.70710678118654752440  # C literal (1/sqrt(2) in Python is 1 ULP lower)
+
+
+def _norm_cdf_sd01(x: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]:
+    """Φ(x; 0, 0.1) with the branch structure of the cephes ``ndtr`` used by the reference."""
+    z = (x / 0.1) * M_SQRT1_2
+    out = np.empty_like(z)
+    for i, v in enumerate(z.tolist()):
+        a = abs(v)
+        if a < M_SQRT1_2:
+            y = 0.5 + 0.5 * math.erf(v)
+        else:
+            y = 0.5 * math.erfc(a)
+            if v > 0:
+                y = 1.0 - y
+        out[i] = y
+    return out
+
+
+def gs_percentile_columns(
+    sorted_scores: npt.NDArray[np.float32], counts: npt.NDArray[np.int64], p: int
+) -> npt.NDArray[np.float64]:
+    """GenomeStudio percentile of the first ``counts[j]`` values of each column of ``sorted_scores`` (double math)."""
+    n = counts.astype(np.int64)
+    pos = n * p / 100.0
+    idx = pos.astype(np.int64)  # int(num*p/100) truncates like the reference
+    idx = np.where(pos - idx < 0.5, idx - 1, idx)
+    cols = np.arange(sorted_scores.shape[1])
+    last = np.maximum(n - 1, 0)
+    lo = np.clip(idx, 0, np.maximum(last - 1, 0))
+    y1 = sorted_scores[lo, cols].astype(np.float64)
+    y2 = sorted_scores[np.minimum(lo + 1, last), cols].astype(np.float64)
+    with np.errstate(invalid="ignore", divide="ignore"):
+        x1 = 100 * (lo + 0.5) / n
+        x2 = 100 * (lo + 1 + 0.5) / n
+        mid = y1 + (y2 - y1) / (x2 - x1) * (p - x1)
+    first = sorted_scores[0, cols].astype(np.float64)
+    lastv = sorted_scores[last, cols].astype(np.float64)
+    out = np.where(idx < 0, first, np.where(idx >= n - 1, lastv, mid))
+    return np.where(n == 0, np.nan, out)
+
+
+def summarize_block(genotypes: npt.NDArray[np.uint8], scores: npt.NDArray[np.float32]) -> dict[str, npt.NDArray]:  # type: ignore[type-arg]
+    """Statistics for a block of loci; inputs are (samples × loci)."""
+    g = np.asarray(genotypes)
+    nc = (g == 0).sum(axis=0).astype(np.int64)
+    aa = (g == 1).sum(axis=0).astype(np.int64)
+    ab = (g == 2).sum(axis=0).astype(np.int64)
+    bb = (g == 3).sum(axis=0).astype(np.int64)
+    calls = aa + ab + bb
+    with np.errstate(invalid="ignore", divide="ignore"):
+        call_freq = np.where(calls + nc > 0, calls / (calls + nc), np.nan)
+        cf = calls.astype(np.float64)
+        aa_f = np.where(calls > 0, aa / cf, np.nan)
+        ab_f = np.where(calls > 0, ab / cf, np.nan)
+        bb_f = np.where(calls > 0, bb / cf, np.nan)
+        a_freq = np.where(calls > 0, (aa * 2 + ab) / (2 * cf), np.nan)
+        q = np.minimum(a_freq, 1.0 - a_freq)
+        p = 1 - q
+        temp = 0.013 / q
+        k = temp * temp * temp * temp
+        dh = ((ab / cf + k) / (2 * p * q + k)) - 1
+        simple = (ab_f - 2 * p * q) / (2 * p * q + HET_EXCESS_EPS)
+    mono = (aa + ab == 0) | (ab + bb == 0)
+    regular = (calls > 0) & ~mono
+    hw = np.zeros(len(calls))
+    het = np.zeros(len(calls))
+    if regular.any():
+        d = dh[regular]
+        cdf = _norm_cdf_sd01(d)
+        hw[regular] = np.where(d < 0, 2 * cdf, 2 * (1 - cdf))
+        het[regular] = d
+    hw[(calls > 0) & mono] = 1.0
+    called_scores = np.where(g != 0, scores.astype(np.float32), np.float32(np.inf))
+    called_scores.sort(axis=0)
+    return {
+        "n_nocalls": nc,
+        "n_calls": calls,
+        "call_freq": call_freq,
+        "aa_freq": aa_f,
+        "ab_freq": ab_f,
+        "bb_freq": bb_f,
+        "minor_freq": q,
+        "gc50": gs_percentile_columns(called_scores, calls, 50),
+        "gc10": gs_percentile_columns(called_scores, calls, 10),
+        "het_excess": het,
+        "chi_test_p100": hw,
+        "het_excess_simple": np.where(calls > 0, simple, np.nan),
+    }

src/simurg/qc/profiles/genotyping.yaml

@@ -0,0 +1,21 @@
+# Starting point for review, not a truth (spec QC-024): recalibrate on your own data with
+# `simurg qc suggest-thresholds`. For LCG/HTS genotyping chips.
+qc:
+  report:
+    callRate: 0.98
+    logRDev: 0.20
+    stainingGreen: 5
+    stainingRed: 5
+    extensionGreen: 5
+    extensionRed: 5
+    hybGreenHighMed: 1
+    hybGreenMediumLow: 1
+    targetRemovalI: 1
+    restoration: null
+    stringency: 2
+    nonSpecificBindingGreen: 1
+    nonSpecificBindingRed: 1
+    npGreen: 5
+    npRed: 5
+    assayControl: 1
+    bgOffset: 3000

src/simurg/qc/profiles/genotyping_ex.yaml

@@ -0,0 +1,21 @@
+# Starting point for review, not a truth (spec QC-024): recalibrate on your own data with
+# `simurg qc suggest-thresholds`. For EX-format genotyping chips (provisional, OQ-046).
+qc:
+  report:
+    callRate: 0.98
+    logRDev: 0.20
+    stainingGreen: 5
+    stainingRed: 5
+    extensionGreen: 5
+    extensionRed: 5
+    hybGreenHighMed: 1
+    hybGreenMediumLow: 1
+    targetRemovalI: 1
+    restoration: null
+    stringency: 2
+    nonSpecificBindingGreen: 1
+    nonSpecificBindingRed: 1
+    npGreen: 5
+    npRed: 5
+    assayControl: 1
+    bgOffset: 3000

src/simurg/qc/profiles/genotyping_example.yaml

@@ -0,0 +1,20 @@
+# Alternative example values from the reference documentation (spec QC-024). Starting point for review only.
+qc:
+  report:
+    callRate: 0.98
+    logRDev: 0.15
+    stainingGreen: 5
+    stainingRed: 5
+    extensionGreen: 5
+    extensionRed: 5
+    hybGreenHighMed: 1.25
+    hybGreenMediumLow: 1.25
+    targetRemovalI: 5
+    restoration: null
+    stringency: 1.5
+    nonSpecificBindingGreen: 5
+    nonSpecificBindingRed: 5
+    npGreen: 5
+    npRed: 5
+    assayControl: 1
+    bgOffset: 3000

src/simurg/qc/profiles/methylation_epic.yaml

@@ -0,0 +1,25 @@
+# Methylation EPIC profile (spec QC-024). Methylation metrics are computed from v2 on.
+qc:
+  report:
+    stainingGreen: 5
+    stainingRed: 5
+    extensionGreen: 5
+    extensionRed: 5
+    hybGreenHighMed: 1
+    hybGreenMediumLow: 1
+    targetRemovalI: 1
+    targetRemovalII: 1
+    restoration: 0
+    npGreen: 5
+    npRed: 5
+    bisulfiteConversionIGreen: 1
+    bisulfiteConversionIRed: 1
+    bisulfiteConversionIGreenBg: 1
+    bisulfiteConversionIRedBg: 1
+    bisulfiteConversionII: 1
+    bisulfiteConversionIIBg: 1
+    specificityIGreen: 1
+    specificityIRed: 1
+    specificityII: 1
+    specificityIIBg: 1
+    bgOffset: 3000

src/simurg/qc/profiles/methylation_msa.yaml

@@ -0,0 +1,25 @@
+# Methylation MSA profile (spec QC-024). Methylation metrics are computed from v2 on.
+qc:
+  report:
+    stainingGreen: 5
+    stainingRed: 5
+    extensionGreen: 5
+    extensionRed: 5
+    hybGreenHighMed: 1
+    hybGreenMediumLow: 1
+    targetRemovalI: 1
+    targetRemovalII: 1
+    restoration: 0
+    npGreen: 2.5
+    npRed: 3
+    bisulfiteConversionIGreen: 1
+    bisulfiteConversionIRed: 1
+    bisulfiteConversionIGreenBg: 0.5
+    bisulfiteConversionIRedBg: 0.5
+    bisulfiteConversionII: 0.5
+    bisulfiteConversionIIBg: 0.5
+    specificityIGreen: 1
+    specificityIRed: 1
+    specificityII: 1
+    specificityIIBg: 1
+    bgOffset: 3000

src/simurg/qc/relatedness.py

@@ -0,0 +1,113 @@
+"""Identity and relatedness (spec QC-033, QC-064; plan §9.4).
+
+* KING-robust kinship between every pair of samples
+  [LIT: Manichaikul A. et al. 2010, Bioinformatics 26:2867, eq. for the within-family (robust) estimator]:
+      φ_ij = [N(Aa,Aa) − 2·N(AA,aa)] / [2·min(N_Aa(i), N_Aa(j))] + 1/2 − [N_Aa(i) + N_Aa(j)] / [4·min(N_Aa(i), N_Aa(j))]
+  with every count taken over the sites called in both samples.
+* Duplicates: genotype concordance over jointly called sites >= 0.99 (plan §9.4) [OURS threshold].
+* Mendelian errors of child/parent trios from the sample sheet (Parent1/Parent2).
+
+Counts are formed with matrix products of per-genotype indicator matrices, so a cohort is O(n²·m) arithmetic in a
+few BLAS calls. Relationship degrees use KING's published cut-offs (> 0.354 duplicate/MZ twin, > 0.177 first
+degree, > 0.0884 second, > 0.0442 third).
+"""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+
+import numpy as np
+import numpy.typing as npt
+
+DEGREE_CUTS = ((0.354, "duplicate"), (0.177, "first"), (0.0884, "second"), (0.0442, "third"))
+DUPLICATE_CONCORDANCE = 0.99
+
+
+@dataclass
+class PairStats:
+    kinship: npt.NDArray[np.float64]  # (n, n); NaN on the diagonal or where undefined
+    concordance: npt.NDArray[np.float64]  # (n, n)
+    n_sites: npt.NDArray[np.int64]  # jointly called sites
+
+
+def degree(phi: float) -> str:
+    for cut, name in DEGREE_CUTS:
+        if phi > cut:
+            return name
+    return "unrelated"
+
+
+def select_sites(
+    counts: npt.NDArray[np.int64],
+    usable: npt.NDArray[np.bool_],
+    *,
+    min_call_rate: float = 0.95,
+    min_maf: float = 0.01,
+    max_sites: int = 100_000,
+) -> npt.NDArray[np.intp]:
+    """High-quality sites from per-locus (no-call, AA, AB, BB) ``counts``: ``usable`` (autosomal, not
+    intensity-only), cohort call rate and MAF filters, thinned evenly to at most ``max_sites``."""
+    called = counts[:, 1:].sum(axis=1)
+    total = counts.sum(axis=1)
+    with np.errstate(invalid="ignore", divide="ignore"):
+        cr = np.where(total > 0, called / total, 0.0)
+        af = np.where(called > 0, (counts[:, 2] + 2 * counts[:, 3]) / (2.0 * called), np.nan)
+    maf = np.minimum(af, 1 - af)
+    idx = np.flatnonzero(usable & (cr >= min_call_rate) & (maf >= min_maf))
+    if len(idx) > max_sites:
+        idx = idx[np.linspace(0, len(idx) - 1, max_sites).round().astype(np.intp)]
+    return idx
+
+
+def pair_stats(genotypes: npt.NDArray[np.uint8]) -> PairStats:
+    """KING-robust kinship and concordance for every pair; ``genotypes`` is (sites × samples), 0 = no call."""
+    g = np.asarray(genotypes)
+    aa = (g == 1).astype(np.float32)
+    ab = (g == 2).astype(np.float32)
+    bb = (g == 3).astype(np.float32)
+    c = (g > 0).astype(np.float32)
+    # float32 products are exact while counts stay below 2**24 sites
+    t = np.float64
+    n_het_het = (ab.T @ ab).astype(t)
+    n_opp = (aa.T @ bb).astype(t)
+    n_opp = n_opp + n_opp.T
+    het_i = (ab.T @ c).astype(t)  # [i, j] = hets of i over sites called in j
+    het_j = het_i.T
+    joint = (c.T @ c).astype(t)
+    same = (aa.T @ aa + ab.T @ ab + bb.T @ bb).astype(t)
+    mn = np.minimum(het_i, het_j)
+    with np.errstate(invalid="ignore", divide="ignore"):
+        phi = (n_het_het - 2 * n_opp) / (2 * mn) + 0.5 - (het_i + het_j) / (4 * mn)
+        conc = same / joint
+    phi[mn == 0] = np.nan
+    np.fill_diagonal(phi, np.nan)
+    np.fill_diagonal(conc, np.nan)
+    return PairStats(kinship=phi, concordance=conc, n_sites=joint.astype(np.int64))
+
+
+def mendel_errors(
+    child: npt.NDArray[np.uint8], father: npt.NDArray[np.uint8] | None, mother: npt.NDArray[np.uint8] | None
+) -> tuple[int, int]:
+    """(errors, informative sites) for one child with one or two genotyped parents; autosomal sites only.
+
+    A site is informative when the child and every given parent are called. Errors: the child carries an allele
+    that no parent can transmit (duo: opposite homozygotes; trio: child genotype outside the parental combinations).
+    """
+    parents = [p for p in (father, mother) if p is not None]
+    if not parents:
+        return 0, 0
+    ok = child > 0
+    for p in parents:
+        ok &= p > 0
+    ch = child[ok].astype(np.int8) - 1  # B-allele count 0/1/2
+    ps = [p[ok].astype(np.int8) - 1 for p in parents]
+    if len(ps) == 1:
+        par = ps[0]
+        err = ((ch == 0) & (par == 2)) | ((ch == 2) & (par == 0))
+    else:
+        f, m = ps
+        # B alleles a parent can transmit: hom-A → {0}, het → {0, 1}, hom-B → {1}
+        lo = f // 2 + m // 2
+        hi = (f + 1) // 2 + (m + 1) // 2
+        err = (ch < lo) | (ch > hi)
+    return int(err.sum()), int(ok.sum())

src/simurg/qc/status.py

@@ -0,0 +1,117 @@
+"""Two-axis sample QC status and informational flags (spec QC-030, QC-031, QC-033, QC-034).
+
+* functional status (PASS/FAIL/NA) — call rate and LogR deviation only;
+* control status (CLEAR/FLAGGED) — control-probe metrics including the assay control;
+* flags — informational, never change either status.
+"""
+
+from __future__ import annotations
+
+import math
+from collections.abc import Mapping, Sequence
+from dataclasses import dataclass, field
+from typing import Any
+
+import numpy as np
+
+from simurg.qc.thresholds import CONTROL_KEYS, KEY_TO_COLUMN, Thresholds
+
+CONTAMINATION_FLAG = 0.02  # plan §9.3 [OURS]
+HET_Z = 3.0
+ASSAY_CONTROL_MIN = 1.0  # OUT-062
+
+
+@dataclass
+class SampleStatus:
+    functional_status: str = "NA"
+    control_status: str = ""
+    control_flag_codes: list[str] = field(default_factory=list)
+    na_metrics: list[str] = field(default_factory=list)
+    flags: list[str] = field(default_factory=list)
+
+    def columns(self) -> dict[str, str]:
+        return {
+            "functional_status": self.functional_status,
+            "control_status": self.control_status,
+            "control_flag_codes": ";".join(self.control_flag_codes),
+            "flags": ";".join(self.flags),
+        }
+
+
+def _num(v: Any) -> float:
+    try:
+        return float(v)
+    except (TypeError, ValueError):
+        return math.nan
+
+
+def functional_status(row: Mapping[str, Any] | None, th: Thresholds) -> str:
+    """QC-030. ``row`` is a sample-summary row; None means no genotyping summary (NA)."""
+    if row is None:
+        return "NA"
+    if row.get("status", "OK") != "OK":
+        return "FAIL"
+    cr, ld = th.get("callRate"), th.get("logRDev")
+    if cr.fails(_num(row.get("autosomal_call_rate"))) or ld.fails(_num(row.get("logr_dev"))):
+        return "FAIL"
+    return "PASS"
+
+
+def control_status(
+    values: Mapping[str, float] | None, na_metrics: Sequence[str], th: Thresholds
+) -> tuple[str, list[str]]:
+    """QC-031: FLAGGED iff an enabled, evaluable control metric is out of range. Codes in QC-004 order."""
+    if values is None:
+        return "", []
+    codes = []
+    for key in CONTROL_KEYS:
+        col = KEY_TO_COLUMN[key]
+        t = th.get(key)
+        if not t.enabled or col not in values or col in na_metrics:
+            continue
+        if t.fails(_num(values[col])):
+            codes.append(key)
+    return ("FLAGGED" if codes else "CLEAR"), codes
+
+
+def het_zscores(het: Sequence[float]) -> list[float]:
+    """z-scores of heterozygosity within the batch (NaN where undefined or with < 3 finite values)."""
+    a = np.array([_num(h) for h in het], dtype=np.float64)
+    ok = np.isfinite(a)
+    if ok.sum() < 3:
+        return [math.nan] * len(a)
+    mu, sd = float(a[ok].mean()), float(a[ok].std(ddof=1))
+    if sd == 0:
+        return [0.0 if o else math.nan for o in ok]
+    return [float((x - mu) / sd) if o else math.nan for x, o in zip(a, ok, strict=True)]
+
+
+def sample_flags(
+    row: Mapping[str, Any],
+    th: Thresholds,
+    *,
+    het_z: float = math.nan,
+    probe_columns: Sequence[str] = (),
+    duplicates: Sequence[str] = (),
+    unexpected_relatedness: bool = False,
+) -> list[str]:
+    """QC-033 flags in a fixed order."""
+    flags: list[str] = []
+    ok = row.get("status", "OK") == "OK"
+    if ok and th.get("callRate").fails(_num(row.get("autosomal_call_rate"))):
+        flags.append("low_call_rate")
+    if ok and th.get("logRDev").fails(_num(row.get("logr_dev"))):
+        flags.append("high_lrr_sd")
+    if row.get("sex_check") == "mismatch":
+        flags.append("sex_mismatch")
+    if math.isfinite(het_z) and abs(het_z) > HET_Z:
+        flags.append("het_outlier")
+    c = _num(row.get("contamination"))
+    if math.isfinite(c) and c > CONTAMINATION_FLAG:
+        flags.append("contaminated")
+    if ok and any(_num(row.get(p)) < ASSAY_CONTROL_MIN for p in probe_columns):
+        flags.append("assay_control_fail")
+    flags.extend(f"duplicate_of:{d}" for d in duplicates)
+    if unexpected_relatedness:
+        flags.append("unexpected_relatedness")
+    return flags

src/simurg/qc/suggest.py

@@ -0,0 +1,98 @@
+"""Data-driven threshold suggestion (spec QC-026, CLI-027) [OURS].
+
+For each metric: robust centre m = median and spread s = 1.4826·MAD over samples, on the log scale for ratio
+metrics (every control metric and the assay-control norm R); suggestion exp(m − k·s) for ``below`` keys and
+exp(m + k·s) for ``above`` keys (linear scale for call rate and LogR deviation). At least ``min_samples`` finite
+values are required; otherwise the key is left unset with a note.
+"""
+
+from __future__ import annotations
+
+import math
+from collections.abc import Mapping, Sequence
+from dataclasses import dataclass
+
+import numpy as np
+import numpy.typing as npt
+
+from simurg.qc.thresholds import ABOVE_KEYS, Thresholds
+
+LINEAR_KEYS = {"callRate", "logRDev"}
+MAD_SCALE = 1.4826
+
+
+@dataclass
+class Suggestion:
+    key: str
+    n: int
+    centre: float  # on the analysis scale (log for ratio metrics)
+    spread: float
+    value: float | None
+    pct_flagged_current: float
+    pct_flagged_suggested: float
+    note: str = ""
+
+
+def _flagged(vals: npt.NDArray[np.float64], thr: float | None, above: bool) -> float:
+    if thr is None or len(vals) == 0:
+        return math.nan
+    bad = vals > thr if above else vals < thr
+    return float(100.0 * bad.mean())
+
+
+def suggest(
+    values: Mapping[str, Sequence[float]],
+    current: Thresholds,
+    *,
+    k: float = 4.0,
+    min_samples: int = 30,
+) -> list[Suggestion]:
+    out = []
+    for key, raw in values.items():
+        v = np.array([float(x) for x in raw], dtype=np.float64)
+        v = v[np.isfinite(v)]
+        above = key in ABOVE_KEYS
+        linear = key in LINEAR_KEYS
+        cur = current.get(key).value
+        if len(v) < min_samples:
+            out.append(
+                Suggestion(
+                    key,
+                    len(v),
+                    math.nan,
+                    math.nan,
+                    None,
+                    _flagged(v, cur, above),
+                    math.nan,
+                    f"fewer than {min_samples} samples",
+                )
+            )
+            continue
+        x = v if linear else np.log(np.maximum(v, 1e-12))
+        m = float(np.median(x))
+        s = MAD_SCALE * float(np.median(np.abs(x - m)))
+        t = m + k * s if above else m - k * s
+        val = t if linear else math.exp(t)
+        if key == "callRate":
+            val = min(max(val, 0.0), 1.0)
+        val = max(val, 0.0)
+        out.append(Suggestion(key, len(v), m, s, val, _flagged(v, cur, above), _flagged(v, val, above)))
+    return out
+
+
+def render_yaml(sugs: Sequence[Suggestion], k: float) -> str:
+    lines = [
+        "# Suggested QC thresholds (simurg qc suggest-thresholds, spec QC-026).",
+        f"# centre = median, spread = {MAD_SCALE}*MAD (log scale for ratio metrics), k = {k}.",
+        "# Review before use: these are starting points computed from your own data.",
+        "qc:",
+        "  report:",
+    ]
+    for s in sugs:
+        info = (
+            f"n={s.n} m={s.centre:.4g} s={s.spread:.4g} flagged: current={s.pct_flagged_current:.2f}% "
+            f"suggested={s.pct_flagged_suggested:.2f}%"
+        )
+        val = "null" if s.value is None else f"{s.value:.6g}"
+        lines.append(f"    {s.key}: {val}  # {s.note or info}")
+    return "\n".join(lines) + "\n"

src/simurg/qc/thresholds.py

@@ -0,0 +1,190 @@
+"""QC thresholds: YAML schema, validation, profiles and precedence (spec QC-020 … QC-025).
+
+Precedence, highest first: CLI > ``--config`` YAML > profile > built-in. Every effective value keeps its source.
+"""
+
+from __future__ import annotations
+
+import math
+from collections.abc import Mapping
+from dataclasses import dataclass, field
+from importlib import resources
+from pathlib import Path
+from typing import Any
+
+import yaml
+
+from simurg.core.errors import SimurgError
+from simurg.qc.controls import DEFAULT_BG_OFFSET, METRICS
+
+FUNCTIONAL_KEYS = ("callRate", "logRDev")
+CONTROL_KEYS = (*(m.key for m in METRICS), "assayControl")
+KNOWN_KEYS = (*FUNCTIONAL_KEYS, *CONTROL_KEYS, "bgOffset")
+ALIASES = {"tgaControl": "assayControl", "call_rate_min": "callRate", "logr_dev_max": "logRDev"}
+ABOVE_KEYS = {"logRDev"}  # flagged when value > threshold; every other key is flagged when value < threshold
+BUILTIN: dict[str, float | None] = {"callRate": 0.98, "logRDev": 0.20}  # QC-023 [EMU]; every other key unset
+PROFILES = ("genotyping", "genotyping_ex", "genotyping_example", "methylation_epic", "methylation_msa")
+KEY_TO_COLUMN = {m.key: m.column for m in METRICS} | {"assayControl": "assay_control_norm_r"}
+
+
+@dataclass(frozen=True)
+class Threshold:
+    value: float | None  # None: check disabled (QC-022)
+    direction: str  # "below" | "above"
+    source: str  # builtin | profile:<name> | config | cli
+
+    @property
+    def enabled(self) -> bool:
+        return self.value is not None
+
+    def fails(self, v: float) -> bool:
+        """True when ``v`` is out of range. Equality passes (QC-020); NaN signal fails (QC-008)."""
+        if self.value is None:
+            return False
+        if v is None or math.isnan(v):
+            return True
+        return v < self.value if self.direction == "below" else v > self.value
+
+
+@dataclass
+class Thresholds:
+    items: dict[str, Threshold] = field(default_factory=dict)
+    profile: str = ""
+
+    def get(self, key: str) -> Threshold:
+        return self.items.get(key, Threshold(None, "above" if key in ABOVE_KEYS else "below", "unset"))
+
+    @property
+    def bg_offset(self) -> float:
+        t = self.items.get("bgOffset")
+        return float(t.value) if t is not None and t.value is not None else DEFAULT_BG_OFFSET
+
+    def as_records(self) -> list[dict[str, Any]]:
+        return [
+            {"key": k, "value": t.value, "direction": t.direction, "source": t.source}
+            for k, t in sorted(self.items.items(), key=lambda kv: KNOWN_KEYS.index(kv[0]))
+        ]
+
+
+def _canon_key(key: str) -> str:
+    key = ALIASES.get(key, key)
+    if key not in KNOWN_KEYS:
+        raise SimurgError(
+            "E_QC_THRESHOLD_KEY", f"unknown QC threshold key {key!r}; valid keys: {', '.join(KNOWN_KEYS)}"
+        )
+    return key
+
+
+def _parse_value(key: str, raw: Any, source: str) -> Threshold:
+    direction = "above" if key in ABOVE_KEYS else "below"
+    if isinstance(raw, Mapping):
+        direction = str(raw.get("direction", direction))
+        if direction not in ("below", "above"):
+            raise SimurgError("E_QC_THRESHOLD_VALUE", f"{key}: direction must be 'below' or 'above', got {direction!r}")
+        raw = raw.get("threshold")
+    if raw is None or (isinstance(raw, str) and raw.strip().lower() in ("", "null", "none", "~")):
+        return Threshold(None, direction, source)
+    if isinstance(raw, bool):
+        raise SimurgError("E_QC_THRESHOLD_VALUE", f"{key}: value must be numeric or null, got {raw!r}")
+    try:
+        v = float(raw)
+    except (TypeError, ValueError):
+        raise SimurgError("E_QC_THRESHOLD_VALUE", f"{key}: value must be numeric or null, got {raw!r}") from None
+    if not math.isfinite(v):
+        raise SimurgError("E_QC_THRESHOLD_VALUE", f"{key}: value must be finite, got {raw!r}")
+    if key == "callRate" and not 0.0 <= v <= 1.0:
+        raise SimurgError("E_QC_THRESHOLD_RANGE", f"callRate must be within [0, 1], got {v}")
+    if v < 0:
+        raise SimurgError("E_QC_THRESHOLD_RANGE", f"{key} must be >= 0, got {v}")
+    return Threshold(v, direction, source)
+
+
+def parse_report_mapping(doc: Any, where: str) -> tuple[dict[str, Any], str | None]:
+    """``qc: {profile, report: {...}}`` (or ``QC:``/``Report:``) → (report mapping, profile)."""
+    if not isinstance(doc, Mapping):
+        raise SimurgError("E_QC_CONFIG_STRUCTURE", f"{where}: expected a 'qc:' mapping with a 'report:' section")
+    top = doc.get("qc", doc.get("QC"))
+    if not isinstance(top, Mapping):
+        raise SimurgError("E_QC_CONFIG_STRUCTURE", f"{where}: missing top-level 'qc:' (or 'QC:') mapping")
+    report = top.get("report", top.get("Report"))
+    if report is None:
+        report = {}
+    if not isinstance(report, Mapping):
+        raise SimurgError("E_QC_CONFIG_STRUCTURE", f"{where}: 'report:' (or 'Report:') must be a mapping")
+    profile = top.get("profile")
+    return dict(report), (str(profile) if profile else None)
+
+
+def load_config(path: Path) -> tuple[dict[str, Any], str | None]:
+    """Read and structurally validate a threshold YAML (QC-021)."""
+    p = Path(path)
+    if not p.exists():
+        raise SimurgError("E_QC_CONFIG_MISSING", f"QC config file was provided but does not exist: {p}", path=str(p))
+    if p.suffix.lower() not in (".yaml", ".yml"):
+        raise SimurgError("E_QC_CONFIG_EXT", f"QC config file must have a .yaml or .yml extension: {p}", path=str(p))
+    try:
+        doc = yaml.safe_load(p.read_text(encoding="utf-8-sig"))
+    except yaml.YAMLError as e:
+        raise SimurgError("E_QC_CONFIG_PARSE", f"QC config file is not valid YAML: {p}: {e}", path=str(p)) from None
+    return parse_report_mapping(doc, str(p))
+
+
+def profile_values(name: str) -> dict[str, Any]:
+    """Shipped profile (QC-024); ``+ffpe`` suffix enables restoration: 1."""
+    base, _, suffix = name.partition("+")
+    if base not in PROFILES or suffix not in ("", "ffpe"):
+        raise SimurgError("E_QC_PROFILE", f"unknown QC profile {name!r}; choose from {', '.join(PROFILES)} or auto")
+    text = resources.files("simurg.qc").joinpath("profiles", f"{base}.yaml").read_text(encoding="utf-8")
+    report, _ = parse_report_mapping(yaml.safe_load(text), f"profile {base}")
+    if suffix == "ffpe":
+        report["restoration"] = 1
+    return report
+
+
+def auto_profile(array_type: str, assay_format: str = "", manifest_name: str = "") -> str:
+    """QC-025."""
+    if array_type == "methylation":
+        return "methylation_msa" if "msa" in manifest_name.lower() else "methylation_epic"
+    return "genotyping_ex" if assay_format.strip().upper() == "EX" else "genotyping"
+
+
+def parse_cli_pairs(pairs: list[str]) -> dict[str, Any]:
+    out: dict[str, Any] = {}
+    for item in pairs:
+        key, sep, val = item.partition("=")
+        if not sep:
+            raise SimurgError("E_QC_THRESHOLD_VALUE", f"--threshold expects KEY=VALUE, got {item!r}")
+        out[key.strip()] = val.strip()
+    return out
+
+
+def resolve_thresholds(
+    *,
+    config: Path | None = None,
+    profile: str | None = None,
+    cli: Mapping[str, Any] | None = None,
+    array_type: str = "genotyping",
+    assay_format: str = "",
+    manifest_name: str = "",
+) -> Thresholds:
+    """Merge built-in < profile < config < CLI (QC-023). A ``profile:`` in the YAML is used when no CLI profile."""
+    layers: list[tuple[str, Mapping[str, Any]]] = [("builtin", BUILTIN)]
+    cfg_report: dict[str, Any] = {}
+    cfg_profile = None
+    if config is not None:
+        cfg_report, cfg_profile = load_config(config)
+    chosen = profile or cfg_profile
+    if chosen == "auto":
+        chosen = auto_profile(array_type, assay_format, manifest_name)
+    if chosen:
+        layers.append((f"profile:{chosen}", profile_values(chosen)))
+    if cfg_report:
+        layers.append(("config", cfg_report))
+    if cli:
+        layers.append(("cli", cli))
+    th = Thresholds(profile=chosen or "")
+    for source, mapping in layers:
+        for raw_key, raw in mapping.items():
+            key = _canon_key(str(raw_key))
+            th.items[key] = _parse_value(key, raw, source)
+    return th

tests/oracle/test_p5_qc.py

@@ -0,0 +1,152 @@
+"""Phase P5 acceptance on P0 data (plan §18 P5).
+
+* BAF-regression contamination equals gtc2vcf ``+BAFregress`` (MIT, built by scripts/make_oracle_gtcs.sh) at its
+  printed precision for all 36 GSA oracle samples, with identical homozygous-site counts.
+* Locus summary equals BeadArrayFiles ``examples/locus_summary.py`` (BSD-2; scripts/make_oracle_qc.sh) exactly for
+  every column except ChiTest_P100 (≤ 1e-13 relative; scipy's normal CDF vs libm, known_unknowns K-018).
+* ``qc controls`` and ``genotype`` write the QC-009/QC-011 tables and the QC columns of the summary.
+"""
+
+from __future__ import annotations
+
+import csv
+import math
+import os
+import subprocess
+from pathlib import Path
+
+import numpy as np
+import pytest
+from tests.conftest import DATA, GSA_BPM, GSA_EGT, GSA_IDATS, ORACLE_GTCS, REFS, need
+
+from simurg.formats.gtc import read_gtc
+from simurg.qc.contamination import baf_regress
+from simurg.qc.locus_summary import summarize_block
+
+pytestmark = pytest.mark.data
+BCFTOOLS = REFS / "tools" / "bcftools-1.22"
+LOCUS_ORACLE = DATA / "oracle" / "locus_summary" / "gsa36_beadarrayfiles.csv"
+
+
+@pytest.fixture(scope="module")
+def cohort():  # type: ignore[no-untyped-def]
+    need(ORACLE_GTCS)
+    files = sorted(ORACLE_GTCS.glob("*.gtc"))
+    gs = [read_gtc(f) for f in files]
+    names = [f.stem for f in files]
+    gt = np.stack([g.genotypes for g in gs], axis=1)  # loci × samples
+    baf = np.stack([g.b_allele_freqs for g in gs], axis=1).astype(np.float32)
+    scores = np.stack([g.genotype_scores for g in gs], axis=0)  # samples × loci
+    return names, gt, baf, scores
+
+
+@pytest.mark.oracle
+def test_baf_regress_matches_gtc2vcf(cohort, tmp_path: Path) -> None:  # type: ignore[no-untyped-def]
+    if not (BCFTOOLS / "plugins" / "BAFregress.so").exists():
+        pytest.skip("BAFregress plugin not built (run scripts/make_oracle_gtcs.sh)")
+    names, gt, baf, _ = cohort
+    vcf = tmp_path / "t.vcf"
+    gtmap = {0: "./.", 1: "0/0", 2: "0/1", 3: "1/1"}
+    with vcf.open("w") as fh:
+        fh.write(
+            "##fileformat=VCFv4.2\n##contig=<ID=1>\n"
+            '##INFO=<ID=ALLELE_A,Number=1,Type=Integer,Description="A">\n'
+            '##INFO=<ID=ALLELE_B,Number=1,Type=Integer,Description="B">\n'
+            '##FORMAT=<ID=GT,Number=1,Type=String,Description="GT">\n'
+            '##FORMAT=<ID=BAF,Number=1,Type=Float,Description="BAF">\n'
+            "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\t" + "\t".join(names) + "\n"
+        )
+        for i in range(gt.shape[0]):
+            cells = "\t".join(
+                f"{gtmap[int(g)]}:{f'{b:.9g}' if np.isfinite(b) else '.'}" for g, b in zip(gt[i], baf[i], strict=True)
+            )
+            fh.write(f"1\t{i + 1}\t.\tA\tC\t.\t.\tALLELE_A=0;ALLELE_B=1\tGT:BAF\t{cells}\n")
+    env = dict(os.environ, BCFTOOLS_PLUGINS=str(BCFTOOLS / "plugins"))
+    out = subprocess.run(
+        [str(BCFTOOLS / "bcftools"), "+BAFregress", str(vcf)], capture_output=True, text=True, env=env, check=True
+    ).stdout
+    ref = {ln.split("\t")[0]: ln.split("\t")[1:] for ln in out.strip().splitlines()[1:]}
+    res = baf_regress(gt, baf)
+    for j, nm in enumerate(names):
+        assert f"{res.estimate[j]:.4f}" == ref[nm][0], nm
+        assert int(res.n_hom[j]) == int(ref[nm][1]), nm
+
+
+@pytest.mark.oracle
+def test_locus_summary_matches_beadarrayfiles(cohort) -> None:  # type: ignore[no-untyped-def]
+    if not LOCUS_ORACLE.exists():
+        pytest.skip("locus-summary oracle missing (run scripts/make_oracle_qc.sh)")
+    _, gt, _, scores = cohort
+    st = summarize_block(gt.T, scores)
+    exact = {
+        "#No_Calls": "n_nocalls",
+        "#Calls": "n_calls",
+        "Call_Freq": "call_freq",
+        "A/A_Freq": "aa_freq",
+        "A/B_Freq": "ab_freq",
+        "B/B_Freq": "bb_freq",
+        "Minor_Freq": "minor_freq",
+        "50%_GC_Score": "gc50",
+        "10%_GC_Score": "gc10",
+        "Het_Excess_Freq": "het_excess",
+    }
+    worst = 0.0
+    with LOCUS_ORACLE.open() as fh:
+        r = csv.reader(fh)
+        next(r)
+        next(r)
+        hdr = next(r)
+        for i, row in enumerate(r):
+            d = dict(zip(hdr, row, strict=True))
+            for col, key in exact.items():
+                ref, ours = float(d[col]), float(st[key][i])
+                assert ref == ours or (math.isnan(ref) and math.isnan(ours)), (i, col, d[col], ours)
+            ref, ours = float(d["ChiTest_P100"]), float(st["chi_test_p100"][i])
+            if ref != ours:
+                worst = max(worst, abs(ours - ref) / abs(ref))
+    assert worst < 1e-13
+
+
+def test_qc_controls_and_genotype_qc_columns(tmp_path: Path) -> None:
+    need(GSA_BPM, GSA_EGT, GSA_IDATS)
+    from typer.testing import CliRunner
+
+    from simurg.cli.main import app
+
+    sub = tmp_path / "idats"
+    chip = sorted(p for p in GSA_IDATS.iterdir() if p.is_dir())[0]
+    (sub / chip.name).mkdir(parents=True)
+    for f in sorted(chip.glob("*.idat"))[:6]:  # three samples
+        (sub / chip.name / f.name).symlink_to(f)
+    runner = CliRunner()
+    out1 = tmp_path / "ctl"
+    res = runner.invoke(
+        app,
+        ["qc", "controls", "--out", str(out1), "--array-type", "genotyping", "--bpm", str(GSA_BPM),
+         "--idats", str(sub), "--quiet"],
+    )  # fmt: skip
+    assert res.exit_code == 0, res.output
+    raw = (out1 / "controls_raw.tsv").read_text().splitlines()
+    assert raw[0].split("\t") == [
+        "sample_id", "category", "bead_address", "control_name", "color_label", "green_intensity", "red_intensity"
+    ]  # fmt: skip
+    assert len(raw) == 1 + 3 * 23
+    qc = list(csv.DictReader((out1 / "controls_qc.tsv").open(), delimiter="\t"))
+    assert len(qc) == 3 and float(qc[0]["staining_green"]) > 5 and qc[0]["na_metrics"] == ""
+
+    out2 = tmp_path / "gt"
+    res = runner.invoke(
+        app,
+        ["genotype", "--bpm", str(GSA_BPM), "--egt", str(GSA_EGT), "--idats", str(sub), "--out", str(out2),
+         "--threads", "1", "--qc-profile", "genotyping", "--quiet"],
+    )  # fmt: skip
+    assert res.exit_code == 0, res.output
+    rows = list(csv.DictReader((out2 / "sample_summary.tsv").open(), delimiter="\t"))
+    for r in rows:
+        assert r["functional_status"] == "PASS" and r["control_status"] in ("CLEAR", "FLAGGED")
+        assert 0.1 < float(r["het_rate"]) < 0.3
+        assert r["contamination"] == ""  # < 10 samples: not estimated (W_QC_TOO_FEW)
+    gqc = list(csv.DictReader((out2 / "controls_qc.tsv").open(), delimiter="\t"))
+    assert [r["staining_green"] for r in gqc] == [r["staining_green"] for r in qc]
+    assert "assay_control_norm_r" in gqc[0]
+    assert (out2 / "relatedness.tsv").exists()

tests/unit/test_genotype.py

@@ -240,9 +240,25 @@ def test_legacy_preset(tmp_path: Path) -> None:
 
     f = _synthetic_run(tmp_path)
     out = tmp_path / "legacy"
-    r = CliRunner().invoke(app, ["genotype", "--quiet", "--bpm", str(f["bpm"]), "--egt", str(f["egt"]),
+    r = CliRunner().invoke(
-                                 "--idats", str(f["idats"]), "--out", str(out), "--threads", "1",
+        app,
-                                 "--preset", "autocall-1.6.2"])
+        [
+            "genotype",
+            "--quiet",
+            "--bpm",
+            str(f["bpm"]),
+            "--egt",
+            str(f["egt"]),
+            "--idats",
+            str(f["idats"]),
+            "--out",
+            str(out),
+            "--threads",
+            "1",
+            "--preset",
+            "autocall-1.6.2",
+        ],
+    )
     assert r.exit_code == 0, r.output
     gtc = next(out.rglob("*.gtc"))
     g = read_gtc(gtc)
@@ -252,6 +268,20 @@ def test_legacy_preset(tmp_path: Path) -> None:
     assert s0["n_calls"] > 0 and s0["gc10"] > 0  # the summary still reports the real statistics
     run = json.loads((out / "run_manifest.json").read_text())
     assert run["effective_options"]["effective_norm_version"]["value"] == "1.1.2"
-    bad = CliRunner().invoke(app, ["genotype", "--bpm", str(f["bpm"]), "--egt", str(f["egt"]), "--idats",
+    bad = CliRunner().invoke(
-                                   str(f["idats"]), "--out", str(tmp_path / "x"), "--preset", "nope"])
+        app,
+        [
+            "genotype",
+            "--bpm",
+            str(f["bpm"]),
+            "--egt",
+            str(f["egt"]),
+            "--idats",
+            str(f["idats"]),
+            "--out",
+            str(tmp_path / "x"),
+            "--preset",
+            "nope",
+        ],
+    )
     assert bad.exit_code == 2

tests/unit/test_qc.py

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