From 4c0ea5e04f25dee518790baf5b03943fc10a2282 Mon Sep 17 00:00:00 2001 From: Kemal Yaylali Date: Sun, 27 Sep 2026 20:54:06 +0100 Subject: [PATCH] P2 (WIP): raw X/Y assembly, canonical Parquet model, `simurg ingest` MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - signal/raw_xy.py: vectorised Infinium I/II raw X/Y and control assembly (plan §7.1) with a per-chip address-map cache; IN-027 coverage gate - core/schemas.py: Arrow schemas (loci, sample_intensity, controls), deterministic atomic Parquet writing - pipeline/ingest.py + CLI: discovery, sheet, IN-010 naming, overwrite check before any work, spawn-based workers, per-sample isolation, --resume by input hashes, run manifest, exit codes 0/3/4 - matches gtc2vcf raw X/Y and controls exactly on 18 GSA samples; outputs byte-identical across thread counts and layouts --- docs/messages.md | 3 + pyproject.toml | 1 + src/simurg/cli/main.py | 51 ++++- src/simurg/core/errors.py | 3 + src/simurg/core/schemas.py | 142 +++++++++++++ src/simurg/formats/bpm.py | 33 ++- src/simurg/pipeline/ingest.py | 343 +++++++++++++++++++++++++++++++ src/simurg/signal/raw_xy.py | 201 ++++++++++++++++++ tests/unit/test_raw_xy_ingest.py | 195 ++++++++++++++++++ 9 files changed, 970 insertions(+), 2 deletions(-) create mode 100644 src/simurg/core/schemas.py create mode 100644 src/simurg/pipeline/ingest.py create mode 100644 src/simurg/signal/raw_xy.py create mode 100644 tests/unit/test_raw_xy_ingest.py diff --git a/docs/messages.md b/docs/messages.md index 8685c4c..87c3ff2 100644 --- a/docs/messages.md +++ b/docs/messages.md @@ -18,6 +18,7 @@ Generated from `src/simurg/core/errors.py` — do not edit by hand. | `E_GTC_UNSUITABLE` | error | GTC lacks fields required by this command | IN-031 | | `E_IDAT_ENCRYPTED` | error | Encrypted (version 1) IDATs are not supported | IN-023 | | `E_IDAT_FORMAT` | error | IDAT file is malformed | IN-023 | +| `E_IDAT_MANIFEST_MISMATCH` | error | IDATs and manifest belong to different products | IN-027 | | `E_IDAT_PAIR_MISMATCH` | error | Red and Grn IDATs of a sample do not match | IN-023 | | `E_INTERNAL` | error | Internal error (bug); rerun with --debug and report | CLI-014 | | `E_NO_INPUT` | error | Neither an input folder nor a sample sheet was given | CLI-006 | @@ -37,6 +38,8 @@ Generated from `src/simurg/core/errors.py` — do not edit by hand. | `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 | +| `W_IDAT_LOW_ADDRESS_COVERAGE` | warning | Some manifest addresses are absent from the IDAT | IN-027 | +| `W_IDAT_MANIFEST_NAME` | warning | IDAT names a different manifest than the BPM | IN-027 | | `W_IDAT_NAME_MISMATCH` | warning | IDAT barcode/position disagrees with its file name | IN-023 | | `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 | diff --git a/pyproject.toml b/pyproject.toml index 5574e89..e6fdc56 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -15,6 +15,7 @@ dependencies = [ "typer>=0.12", # CLI "rich>=13", # progress / human-readable output "pyyaml>=6", # --config files + "pyarrow>=15", # canonical columnar model and Parquet outputs (plan §5) ] [project.optional-dependencies] diff --git a/src/simurg/cli/main.py b/src/simurg/cli/main.py index a814558..23fb9be 100644 --- a/src/simurg/cli/main.py +++ b/src/simurg/cli/main.py @@ -1,6 +1,6 @@ """``simurg`` command-line interface (spec §2, plan §14). -Phase P1 provides ``simurg version`` and ``simurg inspect``. Later phases add the other commands of CLI-001. +Available: ``version``, ``inspect`` (P1), ``ingest`` (P2). Later phases add the other commands of CLI-001. """ from __future__ import annotations @@ -203,5 +203,54 @@ def inspect_cmd( run_guarded(rep, body) +# -------------------------------------------------------------------------------------------------------- +@app.command("ingest") +def ingest_cmd( + bpm: Annotated[Path, typer.Option("--bpm", help="BPM manifest (required).")], + out: Annotated[Path, typer.Option("--out", help="Output folder (required).")], + 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, + csv: Annotated[Path | None, typer.Option("--csv", help="CSV manifest, checked against the BPM.")] = None, + threads: Annotated[int, typer.Option("--threads", min=0, max=1024, help="Worker processes (0 = auto).")] = 0, + force: Annotated[bool, typer.Option("--force", help="Replace existing outputs.")] = False, + resume: Annotated[bool, typer.Option("--resume", help="Keep outputs whose inputs are unchanged.")] = 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", + log_json: Annotated[bool, typer.Option("--log-json", "--json-log", help="JSON Lines log on stderr.")] = False, + debug: Annotated[bool, typer.Option("--debug", help="Debug logging with tracebacks.")] = False, + log_level: Annotated[str, typer.Option("--log-level", help="debug|info|warning|error")] = "info", + quiet: Annotated[bool, typer.Option("--quiet", "--no-progress", help="Suppress info messages.")] = False, +) -> None: + """IDAT pairs + BPM → canonical per-sample raw intensities (Parquet). Phase P2.""" + from simurg.pipeline.ingest import IngestOptions, run_ingest + + if layout not in ("nested", "flat"): + raise typer.BadParameter("--layout must be nested or flat") + if on_duplicate_idat not in ("error", "first", "newest"): + raise typer.BadParameter("--on-duplicate-idat must be error, first or newest") + rep = Reporter( + command="ingest", + level="debug" if debug else log_level, + json_stderr=log_json, + debug=debug, + quiet=quiet, + log_file=out / "logs" / "ingest.log", + ) + opts = IngestOptions( + out=out, + bpm=bpm, + idats=idats, + sheet=sheet, + csv=csv, + threads=threads, + force=force, + resume=resume, + layout=layout, + on_duplicate_idat=on_duplicate_idat, + argv=sys.argv[1:], + ) + run_guarded(rep, lambda: run_ingest(opts, rep)) + + if __name__ == "__main__": # pragma: no cover app() diff --git a/src/simurg/core/errors.py b/src/simurg/core/errors.py index 7f94df8..e21bfd3 100644 --- a/src/simurg/core/errors.py +++ b/src/simurg/core/errors.py @@ -72,6 +72,9 @@ CATALOGUE: dict[str, MessageDef] = dict( _m("E_IDAT_ENCRYPTED", "error", "Encrypted (version 1) IDATs are not supported", "IN-023"), _m("E_IDAT_PAIR_MISMATCH", "error", "Red and Grn IDATs of a sample do not match", "IN-023"), _m("W_IDAT_NAME_MISMATCH", "warning", "IDAT barcode/position disagrees with its file name", "IN-023"), + _m("E_IDAT_MANIFEST_MISMATCH", "error", "IDATs and manifest belong to different products", "IN-027"), + _m("W_IDAT_LOW_ADDRESS_COVERAGE", "warning", "Some manifest addresses are absent from the IDAT", "IN-027"), + _m("W_IDAT_MANIFEST_NAME", "warning", "IDAT names a different manifest than the BPM", "IN-027"), _m("W_IDAT_NO_SCAN_TIME", "warning", "IDAT has no scan entry in its run info", "IN-024"), _m("I_UNLISTED_IDATS", "info", "IDAT pairs not listed in the sample sheet are ignored", "IN-016"), _m("E_DUP_GTC", "error", "Several GTC files resolve to the same sample", "IN-030"), diff --git a/src/simurg/core/schemas.py b/src/simurg/core/schemas.py new file mode 100644 index 0000000..18caadb --- /dev/null +++ b/src/simurg/core/schemas.py @@ -0,0 +1,142 @@ +"""Canonical columnar data model (plan §5) as Apache Arrow schemas, stored as Parquet. + +Schemas carry a ``simurg.schema`` metadata key with a version, so readers can refuse incompatible files. +Parquet files are written without timestamps so identical inputs give byte-identical files (plan §15). +""" + +from __future__ import annotations + +import json +from pathlib import Path +from typing import Any + +import numpy as np +import pyarrow as pa +import pyarrow.parquet as pq + +from simurg.core.atomic import atomic_open +from simurg.formats.bpm import Manifest + +LOCI_SCHEMA_ID = "simurg.loci/1" +INTENSITY_SCHEMA_ID = "simurg.sample_intensity/1" +CONTROLS_SCHEMA_ID = "simurg.sample_controls/1" + +LOCI_SCHEMA = pa.schema( + [ + ("locus_idx", pa.int32()), + ("name", pa.string()), + ("ilmn_id", pa.string()), + ("chrom", pa.string()), + ("pos", pa.int64()), + ("snp", pa.string()), + ("ilmn_strand", pa.string()), + ("source_strand", pa.string()), + ("ref_strand", pa.string()), + ("address_a", pa.int32()), + ("address_b", pa.int32()), + ("assay_type", pa.uint8()), + ("norm_id", pa.uint8()), + ("norm_bin_id", pa.uint8()), + ("norm_lookup", pa.int32()), + ("exp_clusters", pa.int8()), + ("intensity_only", pa.bool_()), + ("ploidy", pa.string()), + ("genome_build", pa.string()), + ], + metadata={"simurg.schema": LOCI_SCHEMA_ID}, +) + +INTENSITY_SCHEMA = pa.schema( + [("raw_x", pa.uint16()), ("raw_y", pa.uint16()), ("missing", pa.bool_())], + metadata={"simurg.schema": INTENSITY_SCHEMA_ID}, +) + +CONTROLS_SCHEMA = pa.schema( + [("address", pa.int32()), ("red", pa.uint16()), ("green", pa.uint16())], + metadata={"simurg.schema": CONTROLS_SCHEMA_ID}, +) + +_WRITE_OPTS: dict[str, Any] = {"compression": "zstd", "use_dictionary": True, "write_statistics": True} + + +def _dict_cols(values: list[str]) -> pa.Array: + return pa.array(values, type=pa.string()) + + +def loci_table(m: Manifest) -> pa.Table: + n = m.num_loci + cols = { + "locus_idx": pa.array(np.arange(n, dtype=np.int32)), + "name": _dict_cols(m.names), + "ilmn_id": _dict_cols(m.ilmn_ids), + "chrom": _dict_cols(m.chroms), + "pos": pa.array(m.map_infos.astype(np.int64)), + "snp": _dict_cols(m.snps), + "ilmn_strand": _dict_cols(m.ilmn_strands), + "source_strand": _dict_cols(m.source_strands), + "ref_strand": _dict_cols(m.ref_strands), + "address_a": pa.array(m.address_a.astype(np.int32)), + "address_b": pa.array(m.address_b.astype(np.int32)), + "assay_type": pa.array(m.assay_types.astype(np.uint8)), + "norm_id": pa.array(m.norm_ids.astype(np.uint8)), + "norm_bin_id": pa.array(m.norm_bin_ids()), + "norm_lookup": pa.array(m.normalization_lookups().astype(np.int32)), + "exp_clusters": pa.array(m.exp_clusters.astype(np.int8)), + "intensity_only": pa.array(m.intensity_only.astype(np.bool_)), + "ploidy": _dict_cols(m.ploidies), + "genome_build": _dict_cols(m.genome_builds), + } + meta = {"simurg.schema": LOCI_SCHEMA_ID, "manifest_name": m.manifest_name, "manifest_format": m.format} + return pa.table(cols, schema=LOCI_SCHEMA.with_metadata(meta)) + + +def with_metadata(table: pa.Table, extra: dict[str, Any]) -> pa.Table: + meta = dict(table.schema.metadata or {}) + meta.update( + {k.encode(): (v if isinstance(v, str) else json.dumps(v, sort_keys=True)).encode() for k, v in extra.items()} + ) + return table.replace_schema_metadata(meta) + + +def intensity_table(raw_x: np.ndarray, raw_y: np.ndarray, missing: np.ndarray, meta: dict[str, Any]) -> pa.Table: + t = pa.table( + { + "raw_x": pa.array(raw_x.astype(np.uint16)), + "raw_y": pa.array(raw_y.astype(np.uint16)), + "missing": pa.array(missing.astype(np.bool_)), + }, + schema=INTENSITY_SCHEMA, + ) + return with_metadata(t, meta) + + +def controls_table(addresses: np.ndarray, red: np.ndarray, green: np.ndarray, meta: dict[str, Any]) -> pa.Table: + t = pa.table( + { + "address": pa.array(addresses.astype(np.int32)), + "red": pa.array(red.astype(np.uint16)), + "green": pa.array(green.astype(np.uint16)), + }, + schema=CONTROLS_SCHEMA, + ) + return with_metadata(t, meta) + + +def write_parquet(table: pa.Table, path: Path) -> None: + """Atomic, deterministic Parquet write.""" + with atomic_open(Path(path), "wb") as fh: + pq.write_table(table, fh, **_WRITE_OPTS) + + +def read_metadata(path: Path) -> dict[str, str]: + md = pq.read_schema(path).metadata or {} + return {k.decode(): v.decode() for k, v in md.items()} + + +def read_table(path: Path, expect_schema: str | None = None) -> pa.Table: + t = pq.read_table(path) + if expect_schema is not None: + got = (t.schema.metadata or {}).get(b"simurg.schema", b"").decode() + if got != expect_schema: + raise ValueError(f"{path}: schema {got!r}, expected {expect_schema!r}") + return t diff --git a/src/simurg/formats/bpm.py b/src/simurg/formats/bpm.py index a856e78..a445fd9 100644 --- a/src/simurg/formats/bpm.py +++ b/src/simurg/formats/bpm.py @@ -25,7 +25,7 @@ Normalization-bin rule (plan §6.2):: from __future__ import annotations import struct -from dataclasses import dataclass, field +from dataclasses import dataclass, field, replace from pathlib import Path from typing import Any @@ -103,6 +103,37 @@ class Manifest: def name_index(self) -> dict[str, int]: return {n: i for i, n in enumerate(self.names)} + def lite(self) -> Manifest: + """Copy keeping only what intensity assembly needs (cheap to pickle to worker processes).""" + n = self.num_loci + empty = [""] * 0 + return replace( + self, + names=list(self.names), + ilmn_ids=empty, + snps=empty, + chroms=empty, + ilmn_strands=empty, + source_strands=empty, + ref_strands=empty, + ploidies=empty, + species=empty, + genome_builds=empty, + sources=empty, + source_versions=empty, + customer_strands=empty, + source_seqs=empty, + top_genomic_seqs=empty, + allele_a_probe_seqs=empty, + allele_b_probe_seqs=empty, + map_infos=np.zeros(0, np.int64), + base_fractions=None, + header_lines=[], + heading={}, + exp_clusters=self.exp_clusters[:n], + intensity_only=self.intensity_only[:n], + ) + def summary(self) -> dict[str, Any]: at = self.assay_types return { diff --git a/src/simurg/pipeline/ingest.py b/src/simurg/pipeline/ingest.py new file mode 100644 index 0000000..6cd3b10 --- /dev/null +++ b/src/simurg/pipeline/ingest.py @@ -0,0 +1,343 @@ +"""``simurg ingest``: IDAT pairs + BPM → canonical per-sample raw intensities (plan §14, §7.1; phase P2). + +Outputs (``--layout nested``):: + + /loci.parquet manifest in canonical schema + /samples//.intensity.parquet raw_x, raw_y, missing (+ provenance metadata) + /samples//.controls.parquet control-probe red/green intensities + /ingest_summary.tsv one row per sample (OK / FAILED / MISSING_INPUT) + /run_manifest.json + +Per-sample failures are isolated (CLI-012); exit code 0/3/4 follows CLI-014. +""" + +from __future__ import annotations + +import multiprocessing as mp +import os +from concurrent.futures import ProcessPoolExecutor +from dataclasses import dataclass, field +from pathlib import Path +from typing import Any + +from simurg import __version__ +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, sha256_file +from simurg.core.schemas import ( + CONTROLS_SCHEMA_ID, + INTENSITY_SCHEMA_ID, + LOCI_SCHEMA_ID, + controls_table, + intensity_table, + loci_table, + read_metadata, + write_parquet, +) +from simurg.formats.bpm import Manifest, read_bpm +from simurg.formats.discovery import ResolvedSample, discover_idats, resolve_input_folder, resolve_samples +from simurg.formats.idat import check_pair, parse_idat_filename, read_idat +from simurg.formats.samplesheet import SampleSheet, output_stem, parse_samplesheet +from simurg.signal.raw_xy import LocusMapCache, assemble_raw_xy, check_coverage + +SUMMARY_COLUMNS = [ + "sample_id", + "sample_name", + "status", + "fail_reason", + "n_missing", + "chip_type", + "scan_time", + "address_coverage", + "grn_idat", + "red_idat", +] + + +@dataclass +class IngestOptions: + out: Path + bpm: Path + idats: Path | None = None + sheet: Path | None = None + csv: Path | None = None + threads: int = 0 + force: bool = False + resume: bool = False + layout: str = "nested" + on_duplicate_idat: str = "error" + argv: list[str] = field(default_factory=list) + + +@dataclass +class SampleJob: + sample_id: str + sample_name: str + stem: str + grn: str + red: str + intensity_path: str + controls_path: str + resume: bool + + +# Per-process state. Workers are started with "spawn": forking after pyarrow has started threads can deadlock. +_STATE: dict[str, Any] = {} + + +def _init_worker(manifest: Manifest, manifest_sha256: str) -> None: + _STATE.clear() + _STATE["manifest"] = manifest + _STATE["manifest_sha256"] = manifest_sha256 + + +def _sample_paths(out: Path, stem: str, layout: str) -> tuple[Path, Path]: + d = out / "samples" / stem if layout == "nested" else out + return d / f"{stem}.intensity.parquet", d / f"{stem}.controls.parquet" + + +def _process(job: SampleJob) -> dict[str, Any]: + manifest: Manifest = _STATE["manifest"] + cache: LocusMapCache = _STATE.setdefault("cache", LocusMapCache(manifest)) + row: dict[str, Any] = { + "sample_id": job.sample_id, + "sample_name": job.sample_name, + "status": "OK", + "fail_reason": "", + "grn_idat": job.grn, + "red_idat": job.red, + "notes": [], + } + try: + hashes = {"grn_sha256": sha256_file(job.grn), "red_sha256": sha256_file(job.red)} + if job.resume and Path(job.intensity_path).exists() and Path(job.controls_path).exists(): + md = read_metadata(Path(job.intensity_path)) + if ( + md.get("grn_sha256") == hashes["grn_sha256"] + and md.get("red_sha256") == hashes["red_sha256"] + and md.get("manifest_sha256") == _STATE["manifest_sha256"] + ): + row.update( + status="OK", + resumed=True, + n_missing=md.get("n_missing", ""), + chip_type=md.get("chip_type", ""), + scan_time=md.get("scan_time", ""), + address_coverage=md.get("address_coverage", ""), + ) + return row + grn, red = read_idat(job.grn), read_idat(job.red) + check_pair(grn, red) + for idat in (grn, red): + parsed = parse_idat_filename(idat.path) + if parsed and idat.barcode and (idat.barcode, idat.position.upper()) != parsed[:2]: + row["notes"].append( + ( + "W_IDAT_NAME_MISMATCH", + f"{idat.path}: internal barcode/position " + f"{idat.barcode}_{idat.position} differs from the file name", + ) + ) + if not idat.scan_time_raw: + row["notes"].append(("W_IDAT_NO_SCAN_TIME", f"{idat.path}: no scan entry in run info")) + cov, notes = check_coverage(grn, manifest) + row["notes"].extend(notes) + raw = assemble_raw_xy(grn, red, manifest, cache=cache) + meta = { + "sample_id": job.sample_id, + "sample_name": job.sample_name, + "barcode": grn.barcode, + "position": grn.position, + "chip_type": grn.chip_type, + "scan_time": grn.scan_time_iso, + "idat_sample_id": grn.sample_id, + "manifest_name": manifest.manifest_name, + "manifest_sha256": _STATE["manifest_sha256"], + "simurg_version": __version__, + "n_missing": str(raw.n_missing), + "address_coverage": f"{cov:.6f}", + **hashes, + } + write_parquet(intensity_table(raw.raw_x, raw.raw_y, raw.missing, meta), Path(job.intensity_path)) + write_parquet( + controls_table(raw.control_addresses, raw.controls_x, raw.controls_y, meta), Path(job.controls_path) + ) + row.update( + n_missing=raw.n_missing, chip_type=grn.chip_type, scan_time=grn.scan_time_iso, address_coverage=f"{cov:.6f}" + ) + except SimurgError as e: + row.update(status="FAILED", fail_reason=f"{e.code}: {e.message}", error_code=e.code) + except Exception as e: + row.update(status="FAILED", fail_reason=f"E_INTERNAL: {type(e).__name__}: {e}", error_code="E_INTERNAL") + return row + + +def _tsv(rows: list[dict[str, Any]]) -> str: + def cell(v: Any) -> str: + return str(v).replace("\t", " ").replace("\n", " ") if v is not None else "" + + lines = ["\t".join(SUMMARY_COLUMNS)] + lines += ["\t".join(cell(r.get(c, "")) for c in SUMMARY_COLUMNS) for r in rows] + return "\n".join(lines) + "\n" + + +def _missing_row(s: ResolvedSample) -> dict[str, Any]: + grn = str(s.pair.grn) if s.pair and s.pair.grn else "" + red = str(s.pair.red) if s.pair and s.pair.red else "" + return { + "sample_id": s.sample_id, + "sample_name": s.sample_name, + "status": "MISSING_INPUT", + "fail_reason": "no complete IDAT pair found", + "grn_idat": grn, + "red_idat": red, + } + + +def run_ingest(opts: IngestOptions, rep: Reporter) -> int: + rm = RunManifest(command="ingest", argv=opts.argv) + rep.run_id = rm.run_id + for k in ("bpm", "idats", "sheet", "csv", "threads", "force", "resume", "layout", "on_duplicate_idat"): + v = getattr(opts, k) + rm.set_option(k, str(v) if v is not None else None, "cli") + + for p, what in ( + (opts.bpm, "BPM"), + (opts.sheet, "sample sheet"), + (opts.csv, "CSV manifest"), + (opts.idats, "IDAT folder"), + ): + if p is not None and not Path(p).exists(): + raise SimurgError("E_FILE_NOT_FOUND", f"{what} not found: {p}", path=str(p)) + + # --- inputs (IN-001 … IN-022) -------------------------------------------------------------------- + sheet: SampleSheet | None = parse_samplesheet(opts.sheet) if opts.sheet else None + if sheet: + rm.sheet_header = sheet.header + for code, msg in sheet.warnings: + rep.warn(code, msg) + folder, notes = resolve_input_folder(opts.idats, sheet) + for code, msg in notes: + rep.warn(code, msg) if code.startswith("W_") else rep.info(msg, code) + disc = discover_idats(folder, on_duplicate=opts.on_duplicate_idat) + for code, msg in disc.warnings: + rep.warn(code, msg) + samples, notes = resolve_samples(disc, sheet) + for code, msg in notes: + rep.info(msg, code) + if not samples: + raise SimurgError("E_NO_INPUT", f"no IDAT pairs found under {folder}") + + # --- names and output paths (IN-010, IN-012) ----------------------------------------------------- + use_idat_names = sheet is None or not sheet.has_sample_name + stems: dict[str, str] = {} + jobs: list[SampleJob] = [] + rows: list[dict[str, Any]] = [] + for s in samples: + if use_idat_names and s.pair is not None and s.pair.grn is not None: + # sheet Sample_Name > IDAT sample name (field 406) > sample_id + try: + idat_name = read_idat(s.pair.grn, load_arrays=False).sample_id.strip() + except SimurgError: + idat_name = "" + if idat_name: + s.sample_name = idat_name + stem = output_stem(s.sample_name) + if stem in stems and stems[stem] != s.sample_id: + if not use_idat_names: # sheet names that sanitise to the same stem + raise SimurgError( + "E_SHEET_NAME_COLLISION", + f"samples {stems[stem]} and {s.sample_id} both map to output name {stem!r}", + ) + rep.info(f"{s.sample_id}: IDAT sample name {s.sample_name!r} is not unique; using the sample ID") + s.sample_name = s.sample_id + stem = output_stem(s.sample_id) + stems[stem] = s.sample_id + if s.status != "READY" or s.pair is None or not s.pair.complete: + rows.append(_missing_row(s)) + continue + ip, cp = _sample_paths(opts.out, stem, opts.layout) + jobs.append( + SampleJob(s.sample_id, s.sample_name, stem, str(s.pair.grn), str(s.pair.red), str(ip), str(cp), opts.resume) + ) + + # --- overwrite policy before any work (CLI-010) -------------------------------------------------- + loci_path = opts.out / "loci.parquet" + summary_path = opts.out / "ingest_summary.tsv" + if not opts.resume: + planned = [loci_path, summary_path] + [Path(p) for j in jobs for p in (j.intensity_path, j.controls_path)] + check_outputs_absent(planned, force=opts.force) + + # --- manifest (IN-040, IN-042) ------------------------------------------------------------------- + rep.info(f"reading manifest {opts.bpm}") + manifest = read_bpm(opts.bpm) + if opts.csv: + from simurg.formats.csv_manifest import compare_manifests, read_csv_manifest + + problems = compare_manifests(manifest, read_csv_manifest(opts.csv, load_sequences=False)) + if problems: + raise SimurgError("E_BPM_CSV_MISMATCH", "BPM and CSV manifests differ:\n " + "\n ".join(problems)) + rm.add_input("bpm", opts.bpm) + if opts.csv: + rm.add_input("csv_manifest", opts.csv) + if opts.sheet: + rm.add_input("sample_sheet", opts.sheet) + opts.out.mkdir(parents=True, exist_ok=True) + write_parquet(loci_table(manifest), loci_path) + + # --- per-sample work (CLI-012, CLI-013) ---------------------------------------------------------- + manifest_sha = sha256_file(opts.bpm) + workers = opts.threads or (os.cpu_count() or 1) + workers = max(1, min(workers, len(jobs) or 1)) + rep.info(f"ingesting {len(jobs)} sample(s) with {workers} worker(s)") + results: list[dict[str, Any]] = [] + lite = manifest.lite() + try: + if workers == 1: + _init_worker(lite, manifest_sha) + results = [_process(j) for j in jobs] + else: + ctx = mp.get_context("spawn") + with ProcessPoolExecutor( + workers, mp_context=ctx, initializer=_init_worker, initargs=(lite, manifest_sha) + ) as ex: + results = list(ex.map(_process, jobs)) + except KeyboardInterrupt: + rm.status = "interrupted" + rm.write(opts.out / "run_manifest.json") + raise + finally: + _STATE.clear() + + # --- reporting --------------------------------------------------------------------------------- + for r in results: + for code, msg in r.pop("notes", []): + rep.warn(code, msg, sample_id=r["sample_id"]) + if r["status"] == "FAILED": + rep.error( + r.get("error_code", "E_INTERNAL"), f"{r['sample_id']}: {r['fail_reason']}", sample_id=r["sample_id"] + ) + rm.failed_samples.append({"sample_id": r["sample_id"], "code": r.get("error_code", "")}) + order = {s.sample_id: i for i, s in enumerate(samples)} + all_rows = sorted(rows + results, key=lambda r: order[r["sample_id"]]) + atomic_write_text(summary_path, _tsv(all_rows)) + + ok = sum(r["status"] == "OK" for r in all_rows) + failed = len(all_rows) - ok + rm.summary = { + "samples": len(all_rows), + "succeeded": ok, + "failed_or_missing": failed, + "resumed": sum(bool(r.get("resumed")) for r in results), + "warnings": rep.warning_counts(), + } + rm.add_output(loci_path, LOCI_SCHEMA_ID) + rm.add_output(summary_path) + rm.schema_versions.update({"*.intensity.parquet": INTENSITY_SCHEMA_ID, "*.controls.parquet": CONTROLS_SCHEMA_ID}) + rm.status = "completed" + rm.write(opts.out / "run_manifest.json") + rep.info(f"done: {ok} succeeded, {failed} failed or missing; summary {summary_path}") + if ok == len(all_rows): + return ExitCode.OK + return ExitCode.PARTIAL if ok else ExitCode.ALL_FAILED diff --git a/src/simurg/signal/raw_xy.py b/src/simurg/signal/raw_xy.py new file mode 100644 index 0000000..73611d4 --- /dev/null +++ b/src/simurg/signal/raw_xy.py @@ -0,0 +1,201 @@ +"""Raw X/Y assembly from a Red/Grn IDAT pair — plan §7.1, spec IN-027. + +For every manifest locus (BPM names order, i.e. GTC order): + +* Infinium II (assay 0): X = Red.Mean[addrA], Y = Grn.Mean[addrA] +* Infinium I A/T (assay 1): X = Red.Mean[addrA], Y = Red.Mean[addrB] +* Infinium I C/G (assay 2): X = Grn.Mean[addrA], Y = Grn.Mean[addrB] + +A locus is (0, 0) — "missing" — unless every address it uses is present and has NBeads >= 2 in the channel(s) +read. Control probes (BPM control config, flattened in order) take X = Red, Y = Grn with no bead-count filter. +[OSS: gtc2vcf ``fill_array`` / ``fill_controls_array`` (MIT); LIT: Steemers 2006 for the chemistry] +""" + +from __future__ import annotations + +import hashlib +from dataclasses import dataclass + +import numpy as np +import numpy.typing as npt + +from simurg.core.errors import SimurgError +from simurg.formats.bpm import ASSAY_INF_I_AT, ASSAY_INF_I_GC, ASSAY_INF_II, Manifest +from simurg.formats.idat import Idat + +MIN_BEADS = 2 +COVERAGE_ERROR = 0.99 + + +@dataclass +class RawIntensities: + """Per-sample raw intensities in manifest order.""" + + raw_x: npt.NDArray[np.uint16] + raw_y: npt.NDArray[np.uint16] + missing: npt.NDArray[np.bool_] # address absent or NBeads < MIN_BEADS + controls_x: npt.NDArray[np.uint16] + controls_y: npt.NDArray[np.uint16] + control_addresses: npt.NDArray[np.int32] + + @property + def n_missing(self) -> int: + return int(self.missing.sum()) + + +class AddressIndex: + """Vectorised address → row lookup for one IDAT (sorted-array search, no per-address Python loop).""" + + def __init__(self, idat: Idat) -> None: + ids = np.asarray(idat.address_ids, dtype=np.int64) + self.order = np.argsort(ids, kind="stable") + self.sorted = ids[self.order] + + def lookup(self, addresses: npt.ArrayLike) -> tuple[npt.NDArray[np.intp], npt.NDArray[np.bool_]]: + """Row index of each address (0 where absent) and a presence mask.""" + a = np.asarray(addresses, dtype=np.int64) + if len(self.sorted) == 0: + return np.zeros(len(a), dtype=np.intp), np.zeros(len(a), dtype=np.bool_) + pos = np.minimum(np.searchsorted(self.sorted, a), len(self.sorted) - 1) + found = self.sorted[pos] == a + return np.where(found, self.order[pos], 0), found + + +def control_addresses(manifest: Manifest) -> npt.NDArray[np.int32]: + """Flattened control addresses from the BPM control config (``a:b:c:d,Category,Color,Name`` per line).""" + out: list[int] = [] + for line in manifest.control_config.splitlines(): + head = line.split(",", 1)[0].strip() + if not head: + continue + for part in head.split(":"): + part = part.strip() + if part: + out.append(int(part)) + return np.asarray(out, dtype=np.int32) + + +def manifest_addresses(manifest: Manifest) -> npt.NDArray[np.int32]: + """Every bead address the manifest's loci read (A for all loci, B for Infinium I).""" + b = manifest.address_b[manifest.address_b != 0] + return np.unique(np.concatenate([manifest.address_a, b])) + + +def address_coverage(idat: Idat, manifest: Manifest) -> float: + """Fraction of manifest bead addresses present in the IDAT (IN-027).""" + wanted = manifest_addresses(manifest) + if len(wanted) == 0: + return 1.0 + _, found = AddressIndex(idat).lookup(wanted) + return float(found.mean()) + + +def check_coverage(idat: Idat, manifest: Manifest) -> tuple[float, list[tuple[str, str]]]: + """IN-027: < 0.99 → E_IDAT_MANIFEST_MISMATCH (raised); < 1 → W_IDAT_LOW_ADDRESS_COVERAGE (returned).""" + cov = address_coverage(idat, manifest) + notes: list[tuple[str, str]] = [] + if cov < COVERAGE_ERROR: + raise SimurgError( + "E_IDAT_MANIFEST_MISMATCH", + f"{idat.path}: only {cov:.2%} of the addresses of manifest {manifest.manifest_name!r} are present " + f"(chip type {idat.chip_type!r}); the IDATs and the manifest belong to different products", + path=idat.path, + ) + if cov < 1.0: + notes.append( + ( + "W_IDAT_LOW_ADDRESS_COVERAGE", + f"{idat.path}: {1 - cov:.4%} of manifest addresses are absent from the IDAT", + ) + ) + if ( + idat.manifest + and manifest.manifest_name + and idat.manifest.split(".")[0].lower() not in manifest.manifest_name.lower() + ): + notes.append( + ( + "W_IDAT_MANIFEST_NAME", + f"{idat.path}: IDAT names manifest {idat.manifest!r}, BPM is {manifest.manifest_name!r}", + ) + ) + return cov, notes + + +@dataclass +class LocusMap: + """Row indices of every manifest address in one IDAT address layout (reusable across samples of a chip).""" + + a: npt.NDArray[np.intp] + a_ok: npt.NDArray[np.bool_] + b: npt.NDArray[np.intp] + b_ok: npt.NDArray[np.bool_] + c: npt.NDArray[np.intp] + c_ok: npt.NDArray[np.bool_] + + +def build_locus_map(idat: Idat, manifest: Manifest, caddr: npt.NDArray[np.int32]) -> LocusMap: + ix = AddressIndex(idat) + a, a_ok = ix.lookup(manifest.address_a) + b, b_ok = ix.lookup(manifest.address_b) + c, c_ok = ix.lookup(caddr) + return LocusMap(a, a_ok, b, b_ok, c, c_ok) + + +class LocusMapCache: + """Caches :class:`LocusMap` per (IDAT address layout, manifest). Samples of one chip type share a layout.""" + + def __init__(self, manifest: Manifest, max_entries: int = 8) -> None: + self.manifest = manifest + self.caddr = control_addresses(manifest) + self._maps: dict[bytes, LocusMap] = {} + self.max_entries = max_entries + + def get(self, idat: Idat) -> LocusMap: + key = hashlib.blake2b(np.ascontiguousarray(idat.address_ids).tobytes(), digest_size=16).digest() + lm = self._maps.get(key) + if lm is None: + if len(self._maps) >= self.max_entries: + self._maps.pop(next(iter(self._maps))) + lm = build_locus_map(idat, self.manifest, self.caddr) + self._maps[key] = lm + return lm + + +def assemble_raw_xy( + grn: Idat, red: Idat, manifest: Manifest, *, min_beads: int = MIN_BEADS, cache: LocusMapCache | None = None +) -> RawIntensities: + """Build per-locus raw X/Y and control intensities (plan §7.1).""" + n = manifest.num_loci + at = manifest.assay_types + if (at > ASSAY_INF_I_GC).any(): + bad = int(np.flatnonzero(at > ASSAY_INF_I_GC)[0]) + raise SimurgError("E_BPM_FORMAT", f"locus {manifest.names[bad]!r} has unknown assay type {int(at[bad])}") + cache = cache or LocusMapCache(manifest) + gm, rm = cache.get(grn), cache.get(red) + + red_mean, grn_mean = red.mean, grn.mean + red_nb, grn_nb = red.n_beads, grn.n_beads + x = np.zeros(n, dtype=np.uint16) + y = np.zeros(n, dtype=np.uint16) + ok = np.zeros(n, dtype=np.bool_) + + # (assay, X source, Y source): (mean, n_beads, row index, present) + plans = ( + (ASSAY_INF_II, (red_mean, red_nb, rm.a, rm.a_ok), (grn_mean, grn_nb, gm.a, gm.a_ok)), + (ASSAY_INF_I_AT, (red_mean, red_nb, rm.a, rm.a_ok), (red_mean, red_nb, rm.b, rm.b_ok)), + (ASSAY_INF_I_GC, (grn_mean, grn_nb, gm.a, gm.a_ok), (grn_mean, grn_nb, gm.b, gm.b_ok)), + ) + for assay, (xm, xnb, xi, xok), (ym, ynb, yi, yok) in plans: + sel = (at == assay) & xok & yok + sel[sel] &= (xnb[xi[sel]] >= min_beads) & (ynb[yi[sel]] >= min_beads) + x[sel] = xm[xi[sel]] + y[sel] = ym[yi[sel]] + ok |= sel + + both = rm.c_ok & gm.c_ok + cx = np.zeros(len(cache.caddr), dtype=np.uint16) + cy = np.zeros(len(cache.caddr), dtype=np.uint16) + cx[both] = red_mean[rm.c[both]] + cy[both] = grn_mean[gm.c[both]] + return RawIntensities(raw_x=x, raw_y=y, missing=~ok, controls_x=cx, controls_y=cy, control_addresses=cache.caddr) diff --git a/tests/unit/test_raw_xy_ingest.py b/tests/unit/test_raw_xy_ingest.py new file mode 100644 index 0000000..e219d14 --- /dev/null +++ b/tests/unit/test_raw_xy_ingest.py @@ -0,0 +1,195 @@ +"""Phase P2: raw X/Y assembly rules (plan §7.1, IN-027) and ``simurg ingest`` behaviour (CLI-010/012/014).""" + +from __future__ import annotations + +from pathlib import Path + +import numpy as np +import pyarrow.parquet as pq +import pytest +from typer.testing import CliRunner + +from simurg.cli.main import app +from simurg.core.errors import SimurgError +from simurg.formats.bpm import read_bpm +from simurg.formats.idat import read_idat +from simurg.signal.raw_xy import address_coverage, assemble_raw_xy, check_coverage, control_addresses +from simurg.synth.files import SynthIdat, SynthLocus, encode_bpm, encode_idat + +runner = CliRunner() +# addresses 1..8 exist; per-address means differ by channel so X/Y sources are identifiable +ADDR = np.arange(1, 9, dtype=np.int32) +RED = np.array([101, 102, 103, 104, 105, 106, 107, 108], np.uint16) +GRN = np.array([201, 202, 203, 204, 205, 206, 207, 208], np.uint16) +BEADS = np.array([10, 10, 10, 10, 1, 10, 10, 10], np.uint8) # address 5 has a single bead + + +def loci() -> list[SynthLocus]: + return [ + SynthLocus("inf2", address_a=1, assay_type=0), + SynthLocus("inf1_at", address_a=2, address_b=3, assay_type=1), + SynthLocus("inf1_gc", address_a=4, address_b=6, assay_type=2), + SynthLocus("low_beads", address_a=5, assay_type=0), + SynthLocus("absent", address_a=99, assay_type=0), + SynthLocus("inf1_absent_b", address_a=7, address_b=98, assay_type=1), + ] + + +def write_sample( + folder: Path, barcode: str, pos: str, *, addr=ADDR, red=RED, grn=GRN, beads=BEADS, sample_id: str = "" +) -> None: # type: ignore[no-untyped-def] + folder.mkdir(parents=True, exist_ok=True) + for ch, mean in (("Red", red), ("Grn", grn)): + s = SynthIdat(address_ids=addr, mean=mean, n_beads=beads, barcode=barcode, position=pos, sample_id=sample_id) + (folder / f"{barcode}_{pos}_{ch}.idat").write_bytes(encode_idat(s)) + + +def bpm(tmp_path: Path, control: str = "0000000008:0000000001,Staining,Red,DNP (High)\n") -> Path: + p = tmp_path / "m.bpm" + p.write_bytes(encode_bpm(loci(), control_config=control)) + return p + + +def test_assembly_rules(tmp_path) -> None: # type: ignore[no-untyped-def] + write_sample(tmp_path / "i", "200000000001", "R01C01") + m = read_bpm(bpm(tmp_path)) + g = read_idat(tmp_path / "i" / "200000000001_R01C01_Grn.idat") + r = read_idat(tmp_path / "i" / "200000000001_R01C01_Red.idat") + raw = assemble_raw_xy(g, r, m) + # inf2: X=Red[1], Y=Grn[1] | A/T: X=Red[2], Y=Red[3] | G/C: X=Grn[4], Y=Grn[6] + assert raw.raw_x.tolist() == [101, 102, 204, 0, 0, 0] + assert raw.raw_y.tolist() == [201, 103, 206, 0, 0, 0] + assert raw.missing.tolist() == [False, False, False, True, True, True] + # controls: flattened addresses in config order, X=Red, Y=Grn, no bead filter + assert control_addresses(m).tolist() == [8, 1] + assert raw.controls_x.tolist() == [108, 101] and raw.controls_y.tolist() == [208, 201] + + +def test_address_order_independent(tmp_path) -> None: # type: ignore[no-untyped-def] + perm = np.array([7, 2, 5, 0, 3, 6, 1, 4]) + write_sample(tmp_path / "a", "200000000001", "R01C01") + write_sample( + tmp_path / "b", "200000000001", "R02C01", addr=ADDR[perm], red=RED[perm], grn=GRN[perm], beads=BEADS[perm] + ) + m = read_bpm(bpm(tmp_path)) + out = [] + for d, pos in (("a", "R01C01"), ("b", "R02C01")): + g = read_idat(tmp_path / d / f"200000000001_{pos}_Grn.idat") + r = read_idat(tmp_path / d / f"200000000001_{pos}_Red.idat") + out.append(assemble_raw_xy(g, r, m)) + assert np.array_equal(out[0].raw_x, out[1].raw_x) and np.array_equal(out[0].raw_y, out[1].raw_y) + + +def test_coverage_gate(tmp_path) -> None: # type: ignore[no-untyped-def] # IN-027 + write_sample(tmp_path / "i", "200000000001", "R01C01") + g = read_idat(tmp_path / "i" / "200000000001_R01C01_Grn.idat") + m = read_bpm(bpm(tmp_path)) + assert address_coverage(g, m) == pytest.approx(7 / 9) # 99 and 98 absent + with pytest.raises(SimurgError) as e: + check_coverage(g, m) + assert e.value.code == "E_IDAT_MANIFEST_MISMATCH" + big = [SynthLocus(f"rs{i}", address_a=i + 1) for i in range(200)] + [SynthLocus("x", address_a=10_000)] + (tmp_path / "big.bpm").write_bytes(encode_bpm(big)) + addr = np.arange(1, 201, dtype=np.int32) + write_sample( + tmp_path / "j", + "200000000001", + "R02C01", + addr=addr, + red=np.ones(200, np.uint16), + grn=np.ones(200, np.uint16), + beads=np.full(200, 5, np.uint8), + ) + gj = read_idat(tmp_path / "j" / "200000000001_R02C01_Grn.idat") + cov, notes = check_coverage(gj, read_bpm(tmp_path / "big.bpm")) + assert cov == pytest.approx(200 / 201) and notes[0][0] == "W_IDAT_LOW_ADDRESS_COVERAGE" + + +def _full_bpm(tmp_path: Path) -> Path: + p = tmp_path / "full.bpm" + p.write_bytes( + encode_bpm( + [SynthLocus(f"rs{i}", address_a=int(a)) for i, a in enumerate(ADDR)], + control_config="0000000001,Staining,Red,DNP (High)\n", + ) + ) + return p + + +def run(*args: str): # type: ignore[no-untyped-def] + return runner.invoke(app, ["ingest", "--quiet", *args]) + + +def test_ingest_end_to_end(tmp_path) -> None: # type: ignore[no-untyped-def] + idats = tmp_path / "idats" + write_sample(idats / "c1", "200000000001", "R01C01", sample_id="NA12878") + write_sample(idats / "c1", "200000000001", "R02C01") + (idats / "c1" / "200000000001_R03C01_Grn.idat").write_bytes(b"IDAT garbage") # corrupt, unpaired + (idats / "c1" / "200000000001_R03C01_Red.idat").write_bytes(b"IDAT garbage") + out = tmp_path / "out" + r = run("--bpm", str(_full_bpm(tmp_path)), "--idats", str(idats), "--out", str(out), "--threads", "2") + assert r.exit_code == 3, r.output # partial failure + rows = [ln.split("\t") for ln in (out / "ingest_summary.tsv").read_text().splitlines()] + status = {row[0]: (row[1], row[2]) for row in rows[1:]} + assert status["200000000001_R01C01"] == ("NA12878", "OK") # IDAT name used (IN-010) + assert status["200000000001_R02C01"] == ("200000000001_R02C01", "OK") + assert status["200000000001_R03C01"][1] == "FAILED" + t = pq.read_table(out / "samples" / "NA12878" / "NA12878.intensity.parquet") + assert t.column("raw_x").to_pylist() == [*RED.tolist()[:4], 0, *RED.tolist()[5:]] + md = {k.decode(): v.decode() for k, v in t.schema.metadata.items()} + assert md["simurg.schema"] == "simurg.sample_intensity/1" and md["n_missing"] == "1" + assert (out / "run_manifest.json").exists() and (out / "loci.parquet").exists() + + # rerun without --force -> refused before any work + r = run("--bpm", str(_full_bpm(tmp_path)), "--idats", str(idats), "--out", str(out)) + assert r.exit_code == 2 and "E_OUTPUT_EXISTS" in r.output + + +def test_ingest_sheet_missing_and_all_failed(tmp_path) -> None: # type: ignore[no-untyped-def] + idats = tmp_path / "idats" + write_sample(idats, "200000000001", "R01C01") + sheet = tmp_path / "s.csv" + sheet.write_text( + "[Data]\nSentrixBarcode_A,SentrixPosition_A,Sample_Name\n200000000001,R01C01,alpha\n200000000001,R09C01,ghost\n" + ) + out = tmp_path / "o" + r = run("--bpm", str(_full_bpm(tmp_path)), "--idats", str(idats), "--sheet", str(sheet), "--out", str(out)) + assert r.exit_code == 3, r.output + txt = (out / "ingest_summary.tsv").read_text() + assert "alpha\tOK" in txt and "ghost\tMISSING_INPUT" in txt + assert (out / "samples" / "alpha" / "alpha.intensity.parquet").exists() + + sheet.write_text("[Data]\nSentrixBarcode_A,SentrixPosition_A\n200000000001,R09C01\n") + r = run( + "--bpm", str(_full_bpm(tmp_path)), "--idats", str(idats), "--sheet", str(sheet), "--out", str(tmp_path / "o2") + ) + assert r.exit_code == 4 + + +def test_ingest_resume_force_flat_and_determinism(tmp_path) -> None: # type: ignore[no-untyped-def] + idats = tmp_path / "idats" + for i in range(1, 5): + write_sample(idats, "200000000001", f"R0{i}C01") + b = str(_full_bpm(tmp_path)) + a_out, b_out = tmp_path / "a", tmp_path / "b" + assert run("--bpm", b, "--idats", str(idats), "--out", str(a_out), "--threads", "1").exit_code == 0 + assert ( + run("--bpm", b, "--idats", str(idats), "--out", str(b_out), "--threads", "3", "--layout", "flat").exit_code == 0 + ) + for f in sorted((a_out / "samples").rglob("*.parquet")): + assert f.read_bytes() == (b_out / f.name).read_bytes(), f.name # threads/layout never change bytes + victim = next((a_out / "samples").rglob("*.intensity.parquet")) + victim.unlink() + assert run("--bpm", b, "--idats", str(idats), "--out", str(a_out), "--resume").exit_code == 0 + import json + + assert json.loads((a_out / "run_manifest.json").read_text())["summary"]["resumed"] == 3 + assert victim.exists() + assert run("--bpm", b, "--idats", str(idats), "--out", str(a_out), "--force").exit_code == 0 + + +def test_ingest_usage_errors(tmp_path) -> None: # type: ignore[no-untyped-def] + r = run("--bpm", str(tmp_path / "nope.bpm"), "--idats", str(tmp_path), "--out", str(tmp_path / "o")) + assert r.exit_code == 2 and "E_FILE_NOT_FOUND" in r.output + r = run("--bpm", str(_full_bpm(tmp_path)), "--out", str(tmp_path / "o")) + assert r.exit_code == 2 and "E_NO_INPUT" in r.output