diff --git a/docs/known_unknowns.md b/docs/known_unknowns.md new file mode 100644 index 0000000..ec4725d --- /dev/null +++ b/docs/known_unknowns.md @@ -0,0 +1,13 @@ +# Known unknowns and findings log + +Plan §20 lists the open questions; this file records what the implementation has learned. Newest first. + +| ID | Date | Topic | Finding | Status | +|---|---|---|---|---| +| K-007 | 2026-09-27 | Demo data | The "GDA v1 Demo 36" download is a GenomeStudio project (`.bsc` + `Data/*.bin`), not IDATs. The GSA v4 PGx demo (`idats_Demo_36.zip`) is used as the P0 IDAT fixture instead. Bioconductor `hapmap370k` is no longer in the current release. | resolved (plan §17 updated) | +| K-006 | 2026-09-27 | Sample sheet | Recognised input columns (`Sample_ID`, `Sample_Name`, …) collide case-insensitively with output columns `sample_id`/`sample_name` (IN-005). Recognised columns are exempt from the reserved-name check and are passed through verbatim. | decision, revisit in P4 (OUT-060 pass-through) | +| K-005 | 2026-09-27 | Sample sheet | Illumina's GSA v4 template sheet has an empty `Sample_Name` column → `E_SHEET_NAME_EMPTY` (IN-010 behaviour as specified). | as specified | +| K-004 | 2026-09-27 | GTC scanner data | PMT gains are 4 bytes read as int32 by BeadArrayFiles and written as float by gtc2vcf. We keep the raw 4 bytes (byte-exact round-trip) and expose the int32 view. | resolved | +| K-003 | 2026-09-27 | Assay type from sequences | gtc2vcf's rule for inferring Infinium I colour from probe/source sequences (needed only for BPM locus entries < v6) agrees with the explicit BPM value for 99.71 % of GSA v4 Infinium I loci. The BPM is authoritative whenever it carries the assay type. | documented; test threshold 99.5 % | +| K-002 | 2026-09-27 | EGT data block (plan U3) | GSA v4 and GDA v1 EGTs are data-block v9 with 15 extra floats per record; the first extra float is the intensity threshold (non-zero for 16 015 GDA loci). The reader tries candidate record sizes and accepts the one whose trailing cluster counts match. dbv 7/8 files still untested. | partly resolved | +| K-001 | 2026-09-27 | IDAT MidBlock (plan U2) | In GSA v4 (BeadChip 48x4), HumanCNV370 (2x12) and 450k (12x8) IDATs the MidBlock count equals N and the section is 4 bytes per entry (int32, as in illuminaio). gtc2vcf reads uint8 but never uses the values. We infer the element size from the section length. | resolved for these chips | diff --git a/docs/messages.md b/docs/messages.md new file mode 100644 index 0000000..8685c4c --- /dev/null +++ b/docs/messages.md @@ -0,0 +1,47 @@ +# SimurgArray message catalogue + +Generated from `src/simurg/core/errors.py` — do not edit by hand. + +| Code | Level | Meaning | Spec | +|---|---|---|---| +| `E_BPM_CSV_MISMATCH` | error | BPM and CSV manifests describe different loci | IN-042 | +| `E_BPM_FORMAT` | error | BPM manifest is malformed or unsupported | IN-040 | +| `E_CONFIG_KEY` | error | Unknown key in --config YAML | CLI-019 | +| `E_CONFIG_PARSE` | error | --config file is not valid YAML | CLI-019 | +| `E_CSV_MANIFEST_FORMAT` | error | CSV manifest is malformed | IN-041 | +| `E_DUP_GTC` | error | Several GTC files resolve to the same sample | IN-030 | +| `E_DUP_IDAT` | error | Several IDAT pairs resolve to the same sample | IN-022 | +| `E_EGT_LOAD` | error | Cluster (EGT) file cannot be loaded | IN-044 | +| `E_FILE_NOT_FOUND` | error | Input file or folder does not exist | CLI-016 | +| `E_FILE_UNREADABLE` | error | Input file cannot be read | CLI-016 | +| `E_GTC_CORRUPT` | warning | GTC file failed validation; sample skipped | IN-032 | +| `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_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 | +| `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_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 | +| `E_SHEET_DUP_SAMPLE` | error | The same barcode/position appears in several rows | IN-009 | +| `E_SHEET_EMPTY` | error | Sample sheet has no data rows | IN-001 | +| `E_SHEET_NAME_COLLISION` | error | Sample names are not unique or clash with sample IDs | IN-010 | +| `E_SHEET_NAME_EMPTY` | error | Sample_Name is empty in at least one row | IN-010 | +| `E_SHEET_NO_KEY_COLUMNS` | error | Sample sheet lacks SentrixBarcode_A/SentrixPosition_A | IN-002 | +| `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_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_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 | +| `W_MANIFEST_NO_REFSTRAND` | warning | Manifest lacks RefStrand annotation | IN-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 | +| `W_SHEET_SEX_CONFLICT` | warning | Reported-sex columns disagree | IN-014 | diff --git a/docs/validation/P1.md b/docs/validation/P1.md new file mode 100644 index 0000000..81198be --- /dev/null +++ b/docs/validation/P1.md @@ -0,0 +1,58 @@ +# Phase P0/P1 validation report + +Date: 2026-09-27 · Machine: Raspberry Pi 5 (aarch64, 4 cores, 8 GB) · Python 3.13, NumPy 2.2. +Commands: `pytest tests/unit` (every CI run), `pytest tests/oracle -m "not slow"` and `-m slow` (needs +`scripts/fetch_data.sh p0`, `scripts/bootstrap_refs.sh`, `scripts/make_oracle_gtcs.sh`). + +## Scope delivered + +| Plan §18 P1 item | Module | Status | +|---|---|---| +| Little-endian binary + .NET string codec | `formats/binutil.py` | done | +| IDAT v3 reader (gzip, metadata, RunInfo → ISO date, encrypted detection) | `formats/idat.py` | done | +| BPM reader (v3–5, locus entries v4/6/7/8, norm-bin rule incl. mod-256 wrap) | `formats/bpm.py` | done | +| CSV manifest reader + BPM↔CSV consistency (IN-042) | `formats/csv_manifest.py` | done | +| EGT reader (data-block variants, count cross-check) | `formats/egt.py` | done | +| GTC v3–v5 reader and writer (layout-preserving) | `formats/gtc.py` | done | +| Sample sheet (IN-001 … IN-014) | `formats/samplesheet.py` | done | +| IDAT discovery, pairing, duplicates, input-set rule (IN-020 … IN-022, CLI-006) | `formats/discovery.py` | done | +| Message catalogue, exit codes, JSON logging, rate limiting | `core/errors.py`, `core/log.py` | done | +| Atomic writes, overwrite policy, run manifest | `core/atomic.py`, `core/provenance.py` | done (used from P4 on) | +| `simurg version`, `simurg inspect [--json]` | `cli/main.py` | done | +| Synthetic IDAT/BPM/EGT writers for tests | `synth/files.py` | done | + +## Acceptance results (plan §18 P1 gate) + +| Criterion | Data | Result | +|---|---|---| +| GTC fields equal to BeadArrayFiles (23 fields incl. raw X/Y, genotypes, scores, BAF, LRR, base calls, sample stats, percentiles, 48 normalization transforms) | 18 GSA v4 GTCs from gtc2vcf `idat2gtc` (preset: Array Analysis CLI) | **EXACT, 18/18** | +| GTC read → write byte-identical | same 18 files; synthetic v3/v5 incl. non-canonical layouts; 25 random property cases | **EXACT** | +| BPM equal to BeadArrayFiles (names, SNP, Chr, MapInfo, AddressA, assay type, normalization lookups, RefStrand, manifest name, control config) | GSA v4 PGx BPM (683 054 loci) | **EXACT** | +| EGT equal to BeadArrayFiles (N, θ/R mean/SD ×3, intensity threshold, GenTrain score, cluster separation, address) | GSA v4 (683 054 records) and GDA v1 (1 904 599 records), ~5 000 evenly spaced records each | **EXACT** | +| BPM ↔ CSV consistency (IN-042: names, counts, Chr, MapInfo, SNP, AddressA/B) | GSA v4 | **0 differences** | +| All IDATs parse and pair; Red/Grn address sets equal; same address set on every chip | 72 GSA IDATs (36 samples, 2 chips) | **pass** | +| Manifest addresses present in IDATs | GSA v4 | **100 %** | +| Bioconductor fixtures: 381 079 / 622 399 addresses (plain + gzip); 2 encrypted files → `E_IDAT_ENCRYPTED` | IlluminaDataTestFiles 1.50.0 | **pass** | +| Malformed inputs never raise uncaught exceptions (truncation, bad magic, bad TOC offsets, bad string lengths, duplicates) | synthetic matrix (EDGE-002, EDGE-004) | **pass** | +| Sample-sheet matrix EDGE-007 (27 cases) incl. 10 000 rows < 2 s (PERF-007) | synthetic | **pass** | +| Locale independence of `inspect --json` (C, tr_TR, de_DE) | synthetic | **pass** | + +## Performance on the Pi (single thread, pure Python/NumPy) + +| Operation | Time | Peak RSS | +|---|---|---| +| IDAT read (734 656 addresses) | 3 ms | – | +| GTC read (683 054 loci) | 34 ms | – | +| EGT read GSA / GDA | 2.2 s / 6.1 s | – | +| CSV manifest read GSA | 12.5 s | 713 MB | +| BPM read GSA | 17.8 s | 665 MB | + +BPM/CSV parsing is the only slow step; it is a one-off per run and moves to Rust in P9 (or is cached). + +## Not yet covered (needs the larger machine / later phases) + +- GEO sweep of 254 genotyping IDAT series and the 1000 Genomes Omni2.5 set (plan §18 P1 wording) — P1 is + validated on P0 data only; rerun `tests/oracle` after `scripts/fetch_data.sh p1`. +- GTCs written by vendor software (only gtc2vcf-written GTCs are available locally); the 1000 Genomes Broad set + contains vendor GTCs. +- IN-027 (IDAT ↔ manifest coverage gate) is P2; the coverage computation is already exercised by a test. diff --git a/src/simurg/__init__.py b/src/simurg/__init__.py new file mode 100644 index 0000000..e9f6f70 --- /dev/null +++ b/src/simurg/__init__.py @@ -0,0 +1,8 @@ +"""SimurgArray — open research toolkit for Illumina Infinium genotyping BeadChip data. + +Specification: SIMURGARRAY_MASTER_PLAN.md (algorithms, roadmap) and +SIMURGARRAY_FUNCTIONAL_SPEC.md (externally visible behaviour, item IDs such as CLI-021). +""" + +__version__ = "0.1.0.dev0" +PRODUCT_NAME = "SimurgArray" diff --git a/src/simurg/cli/__init__.py b/src/simurg/cli/__init__.py new file mode 100644 index 0000000..975b5f5 --- /dev/null +++ b/src/simurg/cli/__init__.py @@ -0,0 +1 @@ +"""SimurgArray cli package.""" diff --git a/src/simurg/cli/main.py b/src/simurg/cli/main.py new file mode 100644 index 0000000..a814558 --- /dev/null +++ b/src/simurg/cli/main.py @@ -0,0 +1,207 @@ +"""``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. +""" + +from __future__ import annotations + +import json +import sys +from collections.abc import Callable +from pathlib import Path +from typing import Annotated, Any, TypeVar + +import numpy as np +import typer + +from simurg import PRODUCT_NAME, __version__ +from simurg.core.errors import ExitCode, SimurgError +from simurg.core.log import Reporter +from simurg.core.provenance import git_revision + +app = typer.Typer( + name="simurg", + help=f"{PRODUCT_NAME}: analysis toolkit for Illumina Infinium genotyping BeadChip data.", + no_args_is_help=True, + add_completion=False, + pretty_exceptions_enable=False, + context_settings={"help_option_names": ["-h", "--help"]}, +) + +F = TypeVar("F", bound=Callable[..., Any]) + + +def version_string() -> str: + rev = git_revision() + return f"{PRODUCT_NAME} {__version__}" + (f"+{rev}" if rev else "") + + +def _version_cb(value: bool) -> None: + if value: + typer.echo(version_string()) + raise typer.Exit(0) + + +@app.callback() +def main( + version: Annotated[ + bool, typer.Option("--version", callback=_version_cb, is_eager=True, help="Print the version and exit.") + ] = False, +) -> None: + """Top-level options.""" + + +def run_guarded(reporter: Reporter, fn: Callable[[], int]) -> None: + """Run a command body mapping errors to catalogue codes and exit codes (CLI-014).""" + reporter.open() + code = ExitCode.OK + try: + code = ExitCode(fn()) + except SimurgError as e: + reporter.error(e.code, e.message + (f"\n hint: {e.hint}" if e.hint else ""), exc=e, file=e.path) + code = e.exit_code + except KeyboardInterrupt: + reporter.error("E_INTERNAL", "interrupted") + code = ExitCode.SIGINT + except Exception as e: + reporter.error( + "E_INTERNAL", f"internal error: {type(e).__name__}: {e}. Rerun with --debug and report it.", exc=e + ) + code = ExitCode.INTERNAL + finally: + reporter.close() + raise typer.Exit(int(code)) + + +# -------------------------------------------------------------------------------------------------------- +@app.command("version") +def version_cmd() -> None: + """Print the version.""" + typer.echo(version_string()) + + +# -------------------------------------------------------------------------------------------------------- +def _jsonable(o: Any) -> Any: + if isinstance(o, np.generic): + return o.item() + if isinstance(o, np.ndarray): + return o.tolist() + if isinstance(o, Path): + return str(o) + if isinstance(o, bytes): + return o.hex() + if isinstance(o, float) and o != o: # NaN + return None + return str(o) + + +def inspect_file(path: Path) -> dict[str, Any]: + from simurg.formats.detect import detect_kind + + kind = detect_kind(path) + info: dict[str, Any] = {"file": str(path), "kind": kind} + if kind == "idat": + from simurg.formats.idat import read_idat + + md = read_idat(path).metadata() + md.pop("path", None) + info.update(md) + elif kind == "bpm": + from simurg.formats.bpm import read_bpm + + info.update(read_bpm(path).summary()) + elif kind == "csv_manifest": + from simurg.formats.csv_manifest import assay_format, read_csv_manifest + + m = read_csv_manifest(path, load_sequences=False) + s = m.summary() + s.pop("normalization_bins", None) + s.pop("infinium_i_at", None) + s.pop("infinium_i_gc", None) + s["heading"] = m.heading + s["assay_format"] = assay_format(m) + info.update(s) + elif kind == "egt": + from simurg.formats.egt import read_egt + + info.update(read_egt(path).summary()) + elif kind == "gtc": + from simurg.formats.gtc import read_gtc, validate_gtc + + g = read_gtc(path, check_write_complete=False) + info.update(g.summary()) + info["problems"] = validate_gtc(g) + elif kind == "samplesheet": + from simurg.formats.samplesheet import parse_samplesheet + + sh = parse_samplesheet(path) + info.update( + { + "rows": len(sh.rows), + "columns": sh.columns, + "key_columns": [sh.barcode_column, sh.position_column], + "root_folder": str(sh.root_folder) if sh.root_folder else None, + "header": sh.header, + "has_sample_name": sh.has_sample_name, + "warnings": [f"{c}: {m}" for c, m in sh.warnings], + } + ) + else: + raise SimurgError("E_FILE_UNREADABLE", f"{path}: unrecognised file type", path=str(path)) + info.pop("path", None) + return info + + +def _print_human(info: dict[str, Any]) -> None: + typer.echo(f"== {info['file']} ({info['kind']})") + for k, v in info.items(): + if k in ("file", "kind"): + continue + if isinstance(v, list) and len(v) > 12: + v = [*v[:12], f"… ({len(v)} total)"] + if isinstance(v, (dict, list)): + v = json.dumps(v, default=_jsonable, ensure_ascii=False) + typer.echo(f" {k}: {v}") + + +@app.command("inspect") +def inspect_cmd( + files: Annotated[list[Path], typer.Argument(help="IDAT, BPM, CSV manifest, EGT, GTC or sample-sheet files.")], + as_json: Annotated[bool, typer.Option("--json", help="Machine-readable JSON output.")] = False, + 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: + """Print the header / table of contents of array files (spec CLI-020).""" + rep = Reporter( + command="inspect", level="debug" if debug else log_level, json_stderr=log_json, debug=debug, quiet=quiet + ) + + def body() -> int: + results = [] + failures = 0 + for f in files: + if not f.exists(): + raise SimurgError("E_FILE_NOT_FOUND", f"{f}: no such file", path=str(f)) + try: + info = inspect_file(f) + except SimurgError as e: + rep.error(e.code, e.message, file=str(f)) + failures += 1 + continue + results.append(info) + if not as_json: + _print_human(info) + if as_json: + out = results[0] if len(results) == 1 else results + sys.stdout.write(json.dumps(out, default=_jsonable, ensure_ascii=False, indent=2) + "\n") + if failures and not results: + return ExitCode.USAGE + return ExitCode.PARTIAL if failures else ExitCode.OK + + run_guarded(rep, body) + + +if __name__ == "__main__": # pragma: no cover + app() diff --git a/src/simurg/cn/__init__.py b/src/simurg/cn/__init__.py new file mode 100644 index 0000000..422b33c --- /dev/null +++ b/src/simurg/cn/__init__.py @@ -0,0 +1 @@ +"""SimurgArray cn package.""" diff --git a/src/simurg/core/__init__.py b/src/simurg/core/__init__.py new file mode 100644 index 0000000..94d5931 --- /dev/null +++ b/src/simurg/core/__init__.py @@ -0,0 +1 @@ +"""SimurgArray core package.""" diff --git a/src/simurg/core/atomic.py b/src/simurg/core/atomic.py new file mode 100644 index 0000000..559e359 --- /dev/null +++ b/src/simurg/core/atomic.py @@ -0,0 +1,86 @@ +"""Atomic output writing and overwrite policy (spec CLI-010, CLI-011). + +Every output is written to ``.part`` in the same directory, fsynced, then renamed, so an interrupted run +never leaves a truncated file under its final name. +""" + +from __future__ import annotations + +import contextlib +import os +from collections.abc import Iterable, Iterator +from pathlib import Path +from typing import IO + +from simurg.core.errors import SimurgError + +PART_SUFFIX = ".part" + + +def _fsync_dir(d: Path) -> None: + with contextlib.suppress(OSError): + fd = os.open(d, os.O_RDONLY) + try: + os.fsync(fd) + finally: + os.close(fd) + + +@contextlib.contextmanager +def atomic_open(path: Path, mode: str = "wb", encoding: str | None = None) -> Iterator[IO[bytes] | IO[str]]: + """Open ``path.part`` for writing; on success fsync and rename to ``path``; on error remove the part file.""" + path = Path(path) + path.parent.mkdir(parents=True, exist_ok=True) + tmp = path.with_name(path.name + PART_SUFFIX) + fh: IO[bytes] | IO[str] = ( + open(tmp, mode) # noqa: SIM115 - closed below + if "b" in mode + else open(tmp, mode, encoding=encoding or "utf-8", newline="") # noqa: SIM115 + ) + try: + yield fh + fh.flush() + os.fsync(fh.fileno()) + fh.close() + os.replace(tmp, path) + _fsync_dir(path.parent) + except BaseException: + fh.close() + with contextlib.suppress(FileNotFoundError): + tmp.unlink() + raise + + +def atomic_write_bytes(path: Path, data: bytes) -> None: + with atomic_open(path, "wb") as fh: + fh.write(data) # type: ignore[arg-type] + + +def atomic_write_text(path: Path, text: str) -> None: + with atomic_open(path, "w", encoding="utf-8") as fh: + fh.write(text) # type: ignore[call-overload] + + +def check_outputs_absent(paths: Iterable[Path], *, force: bool) -> None: + """Pre-flight overwrite check (CLI-010): raise E_OUTPUT_EXISTS listing up to 10 existing paths.""" + if force: + return + existing = [str(p) for p in paths if Path(p).exists()] + if existing: + shown = "\n ".join(existing[:10]) + more = f"\n … and {len(existing) - 10} more" if len(existing) > 10 else "" + raise SimurgError( + "E_OUTPUT_EXISTS", + f"{len(existing)} output file(s) already exist:\n {shown}{more}", + hint="use --force to replace them or --resume to keep completed samples", + ) + + +def remove_stale_parts(directory: Path) -> list[Path]: + """Delete leftover ``*.part`` files (e.g. after a crash); returns what was removed.""" + removed = [] + for p in sorted(Path(directory).rglob(f"*{PART_SUFFIX}")): + with contextlib.suppress(OSError): + p.unlink() + removed.append(p) + return removed diff --git a/src/simurg/core/errors.py b/src/simurg/core/errors.py new file mode 100644 index 0000000..7f94df8 --- /dev/null +++ b/src/simurg/core/errors.py @@ -0,0 +1,134 @@ +"""Stable message catalogue and exception types (spec CLI-014, CLI-016). + +Every error or warning SimurgArray emits has a stable ``code``. The human-readable text may change between +versions; the code may not. ``docs/messages.md`` is generated from :data:`CATALOGUE` +(``python -m simurg.core.errors``). +""" + +from __future__ import annotations + +from dataclasses import dataclass +from enum import IntEnum +from typing import Literal + +Level = Literal["error", "warning", "info"] + + +class ExitCode(IntEnum): + """Process exit codes (spec CLI-014).""" + + OK = 0 + INTERNAL = 1 + USAGE = 2 # usage or pre-flight validation error; no per-sample work started + PARTIAL = 3 # >= 1 sample failed and >= 1 succeeded + ALL_FAILED = 4 # run completed but no sample succeeded + SIGINT = 130 + SIGTERM = 143 + + +@dataclass(frozen=True) +class MessageDef: + code: str + level: Level + summary: str + spec: str + + +def _m(code: str, level: Level, summary: str, spec: str) -> tuple[str, MessageDef]: + return code, MessageDef(code, level, summary, spec) + + +CATALOGUE: dict[str, MessageDef] = dict( + [ + # --- CLI / run control ------------------------------------------------------------------- + _m("E_NO_INPUT", "error", "Neither an input folder nor a sample sheet was given", "CLI-006"), + _m("E_NO_INPUT_FOLDER", "error", "Sample sheet has no RootFolder and no folder was given", "CLI-006"), + _m("E_OUTPUT_EXISTS", "error", "Outputs already exist; use --force or --resume", "CLI-010"), + _m("E_CONFIG_KEY", "error", "Unknown key in --config YAML", "CLI-019"), + _m("E_CONFIG_PARSE", "error", "--config file is not valid YAML", "CLI-019"), + _m("E_FILE_NOT_FOUND", "error", "Input file or folder does not exist", "CLI-016"), + _m("E_FILE_UNREADABLE", "error", "Input file cannot be read", "CLI-016"), + _m("E_INTERNAL", "error", "Internal error (bug); rerun with --debug and report", "CLI-014"), + # --- sample sheet (IN-001 .. IN-017) ------------------------------------------------------- + _m("E_SHEET_DELIMITER", "error", "Sample sheet is not comma-delimited CSV", "IN-001"), + _m("E_SHEET_EMPTY", "error", "Sample sheet has no data rows", "IN-001"), + _m("E_SHEET_NO_KEY_COLUMNS", "error", "Sample sheet lacks SentrixBarcode_A/SentrixPosition_A", "IN-002"), + _m("E_SHEET_DUP_COLUMN", "error", "Duplicate column (column names are case-insensitive)", "IN-003"), + _m("W_SHEET_EMPTY_COLUMN", "warning", "Columns with empty headers are ignored", "IN-004"), + _m("E_SHEET_RESERVED", "error", "Sample sheet column collides with an output column", "IN-005"), + _m("E_SHEET_TOO_LONG", "error", "Sample sheet name or value longer than 500 characters", "IN-006"), + _m("E_SHEET_BARCODE_FORMAT", "error", "Invalid SentrixBarcode_A / SentrixPosition_A value", "IN-007"), + _m("W_SHEET_BARCODE_REPAIRED", "warning", "Barcode written as a decimal was repaired", "IN-007"), + _m("W_SHEET_ROOTFOLDER_OVERRIDDEN", "warning", "CLI folder overrides the sheet RootFolder", "IN-008"), + _m("E_SHEET_DUP_SAMPLE", "error", "The same barcode/position appears in several rows", "IN-009"), + _m("E_SHEET_NAME_EMPTY", "error", "Sample_Name is empty in at least one row", "IN-010"), + _m("E_SHEET_NAME_COLLISION", "error", "Sample names are not unique or clash with sample IDs", "IN-010"), + _m("W_SHEET_SEX_CONFLICT", "warning", "Reported-sex columns disagree", "IN-014"), + # --- IDAT / GTC discovery & parsing (IN-020 .. IN-033) -------------------------------------- + _m("E_DUP_IDAT", "error", "Several IDAT pairs resolve to the same sample", "IN-022"), + _m("W_DUP_IDAT_RESOLVED", "warning", "Duplicate IDAT pair resolved by --on-duplicate-idat", "IN-022"), + _m("W_IDAT_UNPAIRED", "warning", "IDAT channel file without its partner", "IN-021"), + _m("E_IDAT_FORMAT", "error", "IDAT file is malformed", "IN-023"), + _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("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"), + _m("E_GTC_CORRUPT", "warning", "GTC file failed validation; sample skipped", "IN-032"), + _m("W_GTC_FOREIGN", "warning", "GTC was produced by another tool", "IN-031"), + _m("E_GTC_UNSUITABLE", "error", "GTC lacks fields required by this command", "IN-031"), + # --- manifests & cluster files (IN-040 .. IN-047) ------------------------------------------- + _m("E_BPM_FORMAT", "error", "BPM manifest is malformed or unsupported", "IN-040"), + _m("E_CSV_MANIFEST_FORMAT", "error", "CSV manifest is malformed", "IN-041"), + _m("E_BPM_CSV_MISMATCH", "error", "BPM and CSV manifests describe different loci", "IN-042"), + _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"), + ] +) + + +class SimurgError(Exception): + """An error with a stable catalogue code and an exit code.""" + + def __init__( + self, + code: str, + message: str, + *, + exit_code: ExitCode = ExitCode.USAGE, + hint: str | None = None, + path: str | None = None, + ) -> None: + if code not in CATALOGUE: + raise ValueError(f"unknown message code {code!r}") + self.code = code + self.message = message + self.exit_code = exit_code + self.hint = hint + self.path = path + super().__init__(f"[{code}] {message}" + (f" (hint: {hint})" if hint else "")) + + +class FormatError(SimurgError): + """Malformed binary/text input. ``code`` names the file type (E_IDAT_FORMAT, E_GTC_CORRUPT, ...).""" + + +def render_markdown() -> str: + """Render the catalogue as ``docs/messages.md``.""" + lines = [ + "# SimurgArray message catalogue", + "", + "Generated from `src/simurg/core/errors.py` — do not edit by hand.", + "", + "| Code | Level | Meaning | Spec |", + "|---|---|---|---|", + ] + for d in sorted(CATALOGUE.values(), key=lambda m: m.code): + lines.append(f"| `{d.code}` | {d.level} | {d.summary} | {d.spec} |") + return "\n".join(lines) + "\n" + + +if __name__ == "__main__": # pragma: no cover + print(render_markdown(), end="") diff --git a/src/simurg/core/log.py b/src/simurg/core/log.py new file mode 100644 index 0000000..fffa80a --- /dev/null +++ b/src/simurg/core/log.py @@ -0,0 +1,128 @@ +"""Structured message reporting — spec CLI-007 (logging), CLI-008 (rate limiting), CLI-009 (run summary). + +Messages carry a stable catalogue ``code``. Human-readable lines go to stderr; with ``--log-json`` each record is +one JSON object per line (schema ``simurg.log/1``). A log file, when configured, always receives JSON Lines. +""" + +from __future__ import annotations + +import datetime as _dt +import json +import sys +import traceback +from collections import Counter +from dataclasses import dataclass, field +from pathlib import Path +from typing import IO, Any + +from simurg.core.errors import CATALOGUE + +LOG_SCHEMA = "simurg.log/1" +LEVELS = {"debug": 10, "info": 20, "warning": 30, "error": 40} +MAX_PER_SAMPLE = 20 + + +def _now() -> str: + return _dt.datetime.now(_dt.UTC).strftime("%Y-%m-%dT%H:%M:%S.%f")[:-3] + "Z" + + +@dataclass +class Reporter: + command: str = "" + run_id: str = "" + level: str = "info" + json_stderr: bool = False + debug: bool = False + quiet: bool = False + log_file: Path | None = None + stream: IO[str] = field(default_factory=lambda: sys.stderr) + counts: Counter[str] = field(default_factory=Counter) + _per_sample: Counter[tuple[str, str]] = field(default_factory=Counter) + _suppressed: Counter[tuple[str, str]] = field(default_factory=Counter) + _fh: IO[str] | None = None + + def open(self) -> None: + if self.log_file is not None: + self.log_file.parent.mkdir(parents=True, exist_ok=True) + self._fh = self.log_file.open("a", encoding="utf-8") + + def close(self) -> None: + self._flush_suppressed() + if self._fh is not None: + self._fh.close() + self._fh = None + + # -- core --------------------------------------------------------------------------------------------- + def emit(self, level: str, code: str | None, msg: str, *, exc: BaseException | None = None, **ctx: Any) -> None: + if code: + self.counts[code] += 1 + sample = ctx.get("sample_id") or "" + if code and level == "warning": + key = (code, sample) + self._per_sample[key] += 1 + if self._per_sample[key] > MAX_PER_SAMPLE: + self._suppressed[key] += 1 + self._write_file(level, code, msg, exc, ctx) + return + record_exc = "".join(traceback.format_exception(exc)) if (exc is not None and self.debug) else None + self._write_file(level, code, msg, exc, ctx) + if LEVELS.get(level, 20) < LEVELS.get(self.level, 20): + return + if self.quiet and level in ("info", "debug"): + return + if self.json_stderr: + self.stream.write(json.dumps(self._record(level, code, msg, record_exc, ctx), ensure_ascii=False) + "\n") + else: + prefix = {"error": "ERROR", "warning": "WARNING", "info": "INFO", "debug": "DEBUG"}.get(level, level) + code_s = f" [{code}]" if code else "" + self.stream.write(f"{prefix}{code_s}: {msg}\n") + if record_exc: + self.stream.write(record_exc) + self.stream.flush() + + def _record(self, level: str, code: str | None, msg: str, exc: str | None, ctx: dict[str, Any]) -> dict[str, Any]: + rec: dict[str, Any] = { + "schema": LOG_SCHEMA, + "ts": _now(), + "level": level, + "code": code, + "msg": msg, + "command": self.command, + "run_id": self.run_id, + "step": ctx.pop("step", None), + "sample_id": ctx.pop("sample_id", None), + "locus": ctx.pop("locus", None), + "file": ctx.pop("file", None), + "count": ctx.pop("count", 1), + "elapsed_ms": ctx.pop("elapsed_ms", None), + "exc": exc, + } + rec.update(ctx) + return rec + + def _write_file( + self, level: str, code: str | None, msg: str, exc: BaseException | None, ctx: dict[str, Any] + ) -> None: + if self._fh is None: + return + ex = "".join(traceback.format_exception(exc)) if exc is not None else None + self._fh.write(json.dumps(self._record(level, code, msg, ex, dict(ctx)), ensure_ascii=False) + "\n") + + def _flush_suppressed(self) -> None: + for (code, sample), n in sorted(self._suppressed.items()): + where = f" for {sample}" if sample else "" + self.emit("info", None, f"{code} repeated {n} more times{where}") + self._suppressed.clear() + + # -- helpers ------------------------------------------------------------------------------------------ + def warn(self, code: str, msg: str, **ctx: Any) -> None: + self.emit(CATALOGUE[code].level if code in CATALOGUE else "warning", code, msg, **ctx) + + def info(self, msg: str, code: str | None = None, **ctx: Any) -> None: + self.emit("info", code, msg, **ctx) + + def error(self, code: str, msg: str, *, exc: BaseException | None = None, **ctx: Any) -> None: + self.emit("error", code, msg, exc=exc, **ctx) + + def warning_counts(self) -> dict[str, int]: + return dict(sorted(self.counts.items())) diff --git a/src/simurg/core/provenance.py b/src/simurg/core/provenance.py new file mode 100644 index 0000000..70ceeab --- /dev/null +++ b/src/simurg/core/provenance.py @@ -0,0 +1,127 @@ +"""Provenance and the run manifest — plan §15, spec IN-047, OUT-073, CLI-005. + +Timestamps and the run ID live only in ``run_manifest.json`` so that data outputs stay byte-identical across runs. +""" + +from __future__ import annotations + +import datetime as _dt +import hashlib +import json +import os +import platform +import secrets +import sys +from dataclasses import dataclass, field +from pathlib import Path +from typing import Any + +import numpy as np + +from simurg import PRODUCT_NAME, __version__ +from simurg.core.atomic import atomic_write_text + +MANIFEST_SCHEMA = "simurg.run_manifest/1" + + +def sha256_file(path: str | Path, chunk: int = 1 << 20) -> str: + h = hashlib.sha256() + with Path(path).open("rb") as fh: + while b := fh.read(chunk): + h.update(b) + return h.hexdigest() + + +def new_run_id() -> str: + """Sortable, unique run ID: UTC timestamp + 8 random hex digits.""" + return _dt.datetime.now(_dt.UTC).strftime("%Y%m%dT%H%M%SZ") + "-" + secrets.token_hex(4) + + +def git_revision() -> str | None: + here = Path(__file__).resolve() + for parent in here.parents: + head = parent / ".git" / "HEAD" + if head.is_file(): + ref = head.read_text().strip() + if ref.startswith("ref: "): + ref_file = parent / ".git" / ref[5:] + return ref_file.read_text().strip()[:12] if ref_file.is_file() else None + return ref[:12] + return None + + +def environment() -> dict[str, Any]: + return { + "python": sys.version.split()[0], + "numpy": np.__version__, + "platform": platform.platform(), + "machine": platform.machine(), + "cpu_count": os.cpu_count(), + } + + +@dataclass +class RunManifest: + command: str + argv: list[str] + run_id: str = field(default_factory=new_run_id) + started: str = field(default_factory=lambda: _dt.datetime.now(_dt.UTC).isoformat()) + effective_options: dict[str, dict[str, Any]] = field(default_factory=dict) + inputs: list[dict[str, Any]] = field(default_factory=list) + outputs: list[dict[str, Any]] = field(default_factory=list) + sheet_header: dict[str, str] = field(default_factory=dict) + failed_samples: list[dict[str, str]] = field(default_factory=list) + summary: dict[str, Any] = field(default_factory=dict) + schema_versions: dict[str, str] = field(default_factory=dict) + status: str = "running" + seed: int = 42 + + def set_option(self, name: str, value: Any, source: str) -> None: + """Record an effective option with its source ``cli|yaml:|profile:|default`` (CLI-005).""" + self.effective_options[name] = {"value": value, "source": source} + + def add_input(self, role: str, path: str | Path, *, hash_file: bool = True) -> None: + p = Path(path) + entry: dict[str, Any] = {"role": role, "path": str(p.resolve()), "name": p.name} + if p.is_file(): + entry["size"] = p.stat().st_size + if hash_file: + entry["sha256"] = sha256_file(p) + self.inputs.append(entry) + + def add_output(self, path: str | Path, schema: str | None = None) -> None: + p = Path(path) + entry: dict[str, Any] = {"path": str(p.resolve())} + if p.is_file(): + entry["size"] = p.stat().st_size + entry["sha256"] = sha256_file(p) + if schema: + entry["schema"] = schema + self.schema_versions[p.name] = schema + self.outputs.append(entry) + + def to_dict(self) -> dict[str, Any]: + return { + "schema": MANIFEST_SCHEMA, + "product": PRODUCT_NAME, + "version": __version__, + "git_revision": git_revision(), + "run_id": self.run_id, + "command": self.command, + "argv": self.argv, + "status": self.status, + "started": self.started, + "finished": _dt.datetime.now(_dt.UTC).isoformat(), + "seed": self.seed, + "environment": environment(), + "effective_options": self.effective_options, + "inputs": self.inputs, + "outputs": self.outputs, + "sheet_header": self.sheet_header, + "failed_samples": self.failed_samples, + "summary": self.summary, + "schema_versions": self.schema_versions, + } + + def write(self, path: str | Path) -> None: + atomic_write_text(Path(path), json.dumps(self.to_dict(), indent=2, ensure_ascii=False) + "\n") diff --git a/src/simurg/formats/__init__.py b/src/simurg/formats/__init__.py new file mode 100644 index 0000000..3c5d750 --- /dev/null +++ b/src/simurg/formats/__init__.py @@ -0,0 +1 @@ +"""SimurgArray formats package.""" diff --git a/src/simurg/formats/binutil.py b/src/simurg/formats/binutil.py new file mode 100644 index 0000000..2a40ed3 --- /dev/null +++ b/src/simurg/formats/binutil.py @@ -0,0 +1,208 @@ +"""Little-endian binary reading/writing helpers shared by the IDAT, BPM, EGT and GTC codecs. + +All Illumina binary formats are little-endian. Strings are .NET ``BinaryWriter`` strings: a 7-bit varint +(LEB128) length prefix followed by UTF-8 bytes (plan §6; BeadArrayFiles ``read_string``, BSD-2). +Decoding is tolerant (spec IN-033): invalid UTF-8 becomes U+FFFD, while a length running past the end of the +buffer is a :class:`~simurg.core.errors.FormatError`, never an uncaught exception. +""" + +from __future__ import annotations + +import gzip +import io +import mmap +import struct +from pathlib import Path +from typing import Any + +import numpy as np +import numpy.typing as npt + +from simurg.core.errors import FormatError + +_I8 = struct.Struct(" bytes | mmap.mmap: + """Return the file content: memory-mapped when plain, decompressed when gzip.""" + p = Path(path) + with p.open("rb") as fh: + head = fh.read(2) + if head == GZIP_MAGIC: + fh.seek(0) + with gzip.GzipFile(fileobj=fh) as gz: + return gz.read() + fh.seek(0, io.SEEK_END) + if fh.tell() == 0: + return b"" + return mmap.mmap(fh.fileno(), 0, access=mmap.ACCESS_READ) + + +class BinaryReader: + """Bounds-checked sequential reader over an immutable buffer. + + ``err_code`` is the catalogue code raised on malformed input (e.g. ``E_IDAT_FORMAT``). + """ + + __slots__ = ("_buf", "_view", "err_code", "path", "pos", "size") + + def __init__(self, buf: bytes | mmap.mmap | memoryview, err_code: str, path: str | None = None) -> None: + self._buf = buf + self._view = memoryview(buf) + self.size = len(self._view) + self.pos = 0 + self.err_code = err_code + self.path = path + + # -- errors ------------------------------------------------------------------------------------ + def fail(self, message: str) -> FormatError: + where = f"{self.path}: " if self.path else "" + return FormatError(self.err_code, f"{where}{message}", path=self.path) + + def _need(self, n: int, what: str) -> int: + if n < 0 or self.pos + n > self.size: + raise self.fail(f"unexpected end of file reading {what} at offset {self.pos} (need {n} bytes)") + start = self.pos + self.pos += n + return start + + # -- positioning ------------------------------------------------------------------------------- + def seek(self, offset: int, what: str = "section") -> None: + if offset < 0 or offset > self.size: + raise self.fail(f"{what} offset {offset} outside file (size {self.size})") + self.pos = offset + + def remaining(self) -> int: + return self.size - self.pos + + # -- scalars ----------------------------------------------------------------------------------- + def _unpack(self, s: struct.Struct, what: str) -> Any: + start = self._need(s.size, what) + return s.unpack_from(self._view, start)[0] + + def i8(self, what: str = "int8") -> int: + return int(self._unpack(_I8, what)) + + def u8(self, what: str = "uint8") -> int: + return int(self._unpack(_U8, what)) + + def i16(self, what: str = "int16") -> int: + return int(self._unpack(_I16, what)) + + def u16(self, what: str = "uint16") -> int: + return int(self._unpack(_U16, what)) + + def i32(self, what: str = "int32") -> int: + return int(self._unpack(_I32, what)) + + def u32(self, what: str = "uint32") -> int: + return int(self._unpack(_U32, what)) + + def i64(self, what: str = "int64") -> int: + return int(self._unpack(_I64, what)) + + def f32(self, what: str = "float32") -> float: + return float(self._unpack(_F32, what)) + + def raw(self, n: int, what: str = "bytes") -> bytes: + start = self._need(n, what) + return bytes(self._view[start : start + n]) + + # -- strings ----------------------------------------------------------------------------------- + def varint_length(self, what: str = "string length") -> int: + """.NET 7-bit encoded int (at most 5 bytes).""" + result = 0 + shift = 0 + for _ in range(5): + b = self.u8(what) + result |= (b & 0x7F) << shift + if not b & 0x80: + return result + shift += 7 + raise self.fail(f"invalid 7-bit length prefix for {what} at offset {self.pos}") + + def string(self, what: str = "string") -> str: + n = self.varint_length(what) + start = self._need(n, what) + return bytes(self._view[start : start + n]).decode("utf-8", errors="replace") + + def skip_string(self, what: str = "string") -> None: + n = self.varint_length(what) + self._need(n, what) + + # -- arrays ------------------------------------------------------------------------------------ + def array(self, dtype: npt.DTypeLike, count: int, what: str = "array") -> npt.NDArray[Any]: + """Zero-copy little-endian array view (read-only).""" + dt = np.dtype(dtype).newbyteorder("<") + if count < 0: + raise self.fail(f"negative element count {count} for {what}") + start = self._need(count * dt.itemsize, what) + return np.frombuffer(self._view, dtype=dt, count=count, offset=start) + + def counted_array(self, dtype: npt.DTypeLike, what: str = "array") -> npt.NDArray[Any]: + """int32 element count followed by the elements.""" + return self.array(dtype, self.i32(f"{what} count"), what) + + +class BinaryWriter: + """Little-endian writer into an in-memory buffer.""" + + def __init__(self) -> None: + self._buf = io.BytesIO() + + def tell(self) -> int: + return self._buf.tell() + + def getvalue(self) -> bytes: + return self._buf.getvalue() + + def write(self, b: bytes) -> None: + self._buf.write(b) + + def u8(self, v: int) -> None: + self._buf.write(_U8.pack(v)) + + def i16(self, v: int) -> None: + self._buf.write(_I16.pack(v)) + + def u16(self, v: int) -> None: + self._buf.write(_U16.pack(v)) + + def i32(self, v: int) -> None: + self._buf.write(_I32.pack(v)) + + def u32(self, v: int) -> None: + self._buf.write(_U32.pack(v)) + + def f32(self, v: float) -> None: + self._buf.write(_F32.pack(v)) + + def string(self, s: str | bytes) -> None: + data = s.encode("utf-8") if isinstance(s, str) else s + n = len(data) + while True: + b = n & 0x7F + n >>= 7 + if n: + self._buf.write(bytes([b | 0x80])) + else: + self._buf.write(bytes([b])) + break + self._buf.write(data) + + def array(self, arr: npt.ArrayLike, dtype: npt.DTypeLike) -> None: + self._buf.write(np.ascontiguousarray(arr, dtype=np.dtype(dtype).newbyteorder("<")).tobytes()) + + def counted_array(self, arr: npt.ArrayLike, dtype: npt.DTypeLike) -> None: + a = np.asarray(arr) + self.i32(len(a)) + self.array(a, dtype) diff --git a/src/simurg/formats/bpm.py b/src/simurg/formats/bpm.py new file mode 100644 index 0000000..a856e78 --- /dev/null +++ b/src/simurg/formats/bpm.py @@ -0,0 +1,449 @@ +"""BPM (binary bead-pool manifest) reader — plan §6.2, spec IN-040. + +Layout (little-endian, .NET strings):: + + char[3] "BPM"; uint8 1 + int32 version (flag 0x1000 is XOR-ed off; 3..5 supported) + string manifest_name + string control_config (version > 1) + int32 num_loci + int32[num_loci] indexes + string[num_loci] names + uint8[num_loci] norm_ids (< 100) + num_loci × LocusEntry + [int32 m; m × string] trailer header lines (optional) + +LocusEntry field order follows gtc2vcf ``locusentry_read`` (MIT) and BeadArrayFiles ``LocusEntry`` (BSD-2). +Loci are returned in **names order**, which is the order of every per-locus array in a GTC file. + +Normalization-bin rule (plan §6.2):: + + norm_id' = (norm_id + 100 * assay_type) mod 256 # mimics the vendor's byte wrap (Omni5) + lookup = rank of norm_id' among the sorted distinct values +""" + +from __future__ import annotations + +import struct +from dataclasses import dataclass, field +from pathlib import Path +from typing import Any + +import numpy as np +import numpy.typing as npt + +from simurg.core.errors import FormatError +from simurg.formats.binutil import BinaryReader, load_bytes + +ASSAY_INF_II = 0 +ASSAY_INF_I_AT = 1 # red channel +ASSAY_INF_I_GC = 2 # green channel + +_I32 = struct.Struct("= 7) + header_lines: list[str] = field(default_factory=list) + controls: list[dict[str, str]] = field(default_factory=list) + heading: dict[str, str] = field(default_factory=dict) + + @property + def num_loci(self) -> int: + return len(self.names) + + @property + def has_ref_strand(self) -> bool: + return any(s for s in self.ref_strands) + + def norm_bin_ids(self) -> npt.NDArray[np.uint8]: + """``(norm_id + 100 * assay_type) mod 256`` per locus.""" + return ((self.norm_ids.astype(np.int64) + 100 * self.assay_types.astype(np.int64)) % 256).astype(np.uint8) + + def normalization_lookups(self) -> npt.NDArray[np.int32]: + """Index into the GTC's normalization-transform list for every locus.""" + ids = self.norm_bin_ids() + uniq = np.unique(ids) # sorted + return np.searchsorted(uniq, ids).astype(np.int32) + + def name_index(self) -> dict[str, int]: + return {n: i for i, n in enumerate(self.names)} + + def summary(self) -> dict[str, Any]: + at = self.assay_types + return { + "path": self.path, + "format": self.format, + "version": self.version, + "manifest_name": self.manifest_name, + "num_loci": self.num_loci, + "entry_versions": sorted({int(v) for v in np.unique(self.entry_versions)}), + "infinium_ii": int((at == ASSAY_INF_II).sum()), + "infinium_i_at": int((at == ASSAY_INF_I_AT).sum()), + "infinium_i_gc": int((at == ASSAY_INF_I_GC).sum()), + "intensity_only": int(self.intensity_only.sum()), + "normalization_bins": len(np.unique(self.norm_bin_ids())), + "has_ref_strand": self.has_ref_strand, + "n_controls": len(self.controls), + "genome_builds": sorted(set(self.genome_builds) - {""}), + "chromosomes": sorted(set(self.chroms), key=_chrom_key), + } + + +def _chrom_key(c: str) -> tuple[int, str]: + return (int(c), "") if c.isdigit() else (1000, c) + + +def parse_control_config(text: str) -> list[dict[str, str]]: + """Control-config lines ``addr[:addr..],Category,Color,Name`` → dicts.""" + out: list[dict[str, str]] = [] + for line in text.splitlines(): + parts = [p.strip() for p in line.split(",")] + if len(parts) < 4 or not parts[0]: + continue + out.append({"address": parts[0], "category": parts[1], "color": parts[2], "name": ",".join(parts[3:])}) + return out + + +def read_bpm(path: str | Path) -> Manifest: + """Parse a BPM file. Raises FormatError(E_BPM_FORMAT).""" + spath = str(path) + buf = load_bytes(path) + r = BinaryReader(buf, "E_BPM_FORMAT", spath) + if r.size < 8 or r.raw(3, "magic") != b"BPM": + raise r.fail("bad magic (not a BPM file)") + if (fmt := r.u8("format byte")) != 1: + raise r.fail(f"unknown BPM format byte {fmt}") + version = r.i32("version") + if version & 0x1000: + version ^= 0x1000 + if not 3 <= version <= 5: + raise r.fail(f"unsupported BPM version {version}") + manifest_name = r.string("manifest name") + control_config = r.string("control config") if version > 1 else "" + n = r.i32("number of loci") + if n <= 0 or n > 50_000_000: + raise r.fail(f"implausible number of loci {n}") + r.array("= 100).any(): + raise r.fail("invalid normalization ID (>= 100)") + name_to_idx: dict[str, int] = {} + for i, nm in enumerate(names): + if nm in name_to_idx: + raise r.fail(f"locus name {nm!r} present more than once") + name_to_idx[nm] = i + + cols = _parse_locus_entries(r, n, name_to_idx) + + header_lines: list[str] = [] + if r.remaining() >= 4: + m = r.i32("trailer count") + if 0 <= m <= 1_000_000: + header_lines = [r.string("trailer line") for _ in range(m)] + + return Manifest( + path=spath, + format="bpm", + version=version, + manifest_name=manifest_name, + control_config=control_config, + names=names, + norm_ids=norm_ids, + header_lines=header_lines, + controls=parse_control_config(control_config), + **cols, + ) + + +def _parse_locus_entries(r: BinaryReader, n: int, name_to_idx: dict[str, int]) -> dict[str, Any]: + """Hot loop: parse ``n`` LocusEntry records (optimised pure Python; Rust in phase P9).""" + view = memoryview(r._buf) + size = r.size + pos = r.pos + intern: dict[bytes, str] = {} + unpack_i32 = _I32.unpack_from + unpack_4f = _F32X4.unpack_from + + def fail(msg: str) -> FormatError: + r.pos = pos + return r.fail(msg) + + def s_at(p: int) -> tuple[str, int]: + # .NET 7-bit length-prefixed UTF-8 string, interned + if p >= size: + raise fail("unexpected end of file in locus entry") + b = view[p] + p += 1 + ln = b & 0x7F + shift = 7 + while b & 0x80: + if p >= size or shift > 28: + raise fail("bad string length in locus entry") + b = view[p] + p += 1 + ln |= (b & 0x7F) << shift + shift += 7 + end = p + ln + if end > size: + raise fail("string runs past end of file in locus entry") + raw = bytes(view[p:end]) + s = intern.get(raw) + if s is None: + s = raw.decode("utf-8", errors="replace") + if ln <= 64: + intern[raw] = s + return s, end + + str_cols = [ + "ilmn_ids", + "snps", + "chroms", + "ilmn_strands", + "source_strands", + "ref_strands", + "ploidies", + "species", + "genome_builds", + "sources", + "source_versions", + "customer_strands", + "source_seqs", + "top_genomic_seqs", + "allele_a_probe_seqs", + "allele_b_probe_seqs", + ] + out: dict[str, list[str]] = {c: [""] * n for c in str_cols} + map_infos = np.zeros(n, dtype=np.int64) + address_a = np.zeros(n, dtype=np.int32) + address_b = np.zeros(n, dtype=np.int32) + assay = np.zeros(n, dtype=np.uint8) + exp_clusters = np.zeros(n, dtype=np.int8) + intensity_only = np.zeros(n, dtype=np.bool_) + versions = np.zeros(n, dtype=np.int8) + fracs = np.zeros((n, 4), dtype=np.float32) + seen = np.zeros(n, dtype=np.bool_) + any_v7 = False + + o_ilmn, o_snp, o_chr = out["ilmn_ids"], out["snps"], out["chroms"] + o_istr, o_sstr, o_rstr = out["ilmn_strands"], out["source_strands"], out["ref_strands"] + o_ploidy, o_species, o_build = out["ploidies"], out["species"], out["genome_builds"] + o_source, o_sver, o_cstr = out["sources"], out["source_versions"], out["customer_strands"] + o_sseq, o_top, o_pa, o_pb = ( + out["source_seqs"], + out["top_genomic_seqs"], + out["allele_a_probe_seqs"], + out["allele_b_probe_seqs"], + ) + + for _ in range(n): + if pos + 4 > size: + raise fail("unexpected end of file at locus entry") + ver = unpack_i32(view, pos)[0] + pos += 4 + if ver not in (4, 6, 7, 8): + raise fail(f"unsupported locus entry version {ver}") + ilmn_id, pos = s_at(pos) + name, pos = s_at(pos) + _, pos = s_at(pos) + _, pos = s_at(pos) + _, pos = s_at(pos) + if pos + 4 > size: + raise fail("unexpected end of file at locus index") + pos += 4 # index (1-based); BeadArrayFiles maps entries by name, as we do + _, pos = s_at(pos) + ilmn_strand, pos = s_at(pos) + snp, pos = s_at(pos) + chrom, pos = s_at(pos) + ploidy, pos = s_at(pos) + species, pos = s_at(pos) + map_info, pos = s_at(pos) + top_seq, pos = s_at(pos) + cust_strand, pos = s_at(pos) + if pos + 8 > size: + raise fail("unexpected end of file at addresses") + a_addr = unpack_i32(view, pos)[0] + b_addr = unpack_i32(view, pos + 4)[0] + pos += 8 + probe_a, pos = s_at(pos) + probe_b, pos = s_at(pos) + build, pos = s_at(pos) + source, pos = s_at(pos) + source_ver, pos = s_at(pos) + source_strand, pos = s_at(pos) + source_seq, pos = s_at(pos) + at = -1 + exp_cl = 0 + int_only = 0 + if ver >= 6: + if pos + 4 > size: + raise fail("unexpected end of file at assay type") + exp_cl = int.from_bytes(view[pos + 1 : pos + 2], "little", signed=True) + int_only = view[pos + 2] + at = view[pos + 3] + pos += 4 + if at > 2: + raise fail(f"invalid assay type {at}") + if (b_addr == 0) != (at == 0): + raise fail("assay type is inconsistent with address B") + f4 = None + if ver >= 7: + if pos + 16 > size: + raise fail("unexpected end of file at base fractions") + f4 = unpack_4f(view, pos) + pos += 16 + ref_strand = "" + if ver >= 8: + ref_strand, pos = s_at(pos) + + idx = name_to_idx.get(name) + if idx is None: + raise fail(f"locus entry {name!r} not in the names table") + if seen[idx]: + raise fail(f"locus entry {name!r} appears twice") + seen[idx] = True + if at < 0: + at = infer_assay_type(probe_a, probe_b, source_seq) + o_ilmn[idx] = ilmn_id + o_snp[idx] = snp + o_chr[idx] = chrom + o_istr[idx] = ilmn_strand + o_sstr[idx] = source_strand + o_rstr[idx] = ref_strand + o_ploidy[idx] = ploidy + o_species[idx] = species + o_build[idx] = build + o_source[idx] = source + o_sver[idx] = source_ver + o_cstr[idx] = cust_strand + o_sseq[idx] = source_seq + o_top[idx] = top_seq + o_pa[idx] = probe_a + o_pb[idx] = probe_b + try: + map_infos[idx] = int(map_info) if map_info else 0 + except ValueError: + raise fail(f"non-integer MapInfo {map_info!r} for {name!r}") from None + address_a[idx] = a_addr + address_b[idx] = b_addr + assay[idx] = at + exp_clusters[idx] = exp_cl + intensity_only[idx] = bool(int_only) + versions[idx] = ver + if f4 is not None: + fracs[idx] = f4 + any_v7 = True + + if not seen.all(): + raise fail(f"{int((~seen).sum())} loci in the names table have no locus entry") + r.pos = pos + return { + **out, + "map_infos": map_infos, + "address_a": address_a, + "address_b": address_b, + "assay_types": assay, + "exp_clusters": exp_clusters, + "intensity_only": intensity_only, + "entry_versions": versions, + "base_fractions": fracs if any_v7 else None, + } + + +_IUPAC = { + "A": 1, + "C": 2, + "G": 4, + "T": 8, + "R": 5, + "Y": 10, + "S": 6, + "W": 9, + "K": 12, + "M": 3, + "B": 14, + "D": 13, + "H": 11, + "V": 7, + "N": 15, +} + + +def infer_assay_type(probe_a: str, probe_b: str, source_seq: str) -> int: + """Assay type for locus entries older than v6 (gtc2vcf ``get_assay_type``, MIT). + + Returns 0 for Infinium II, 1 for Infinium I red (A/T), 2 for Infinium I green (C/G), 255 if unknown. + """ + if not probe_a or not source_seq: + return 255 + if not probe_b: + return 0 + left = source_seq.find("[") + right = source_seq.find("]") + if left <= 0 or right < 0 or right + 1 >= len(source_seq): + return 255 + tl = source_seq[left - 1].upper() + tr = source_seq[right + 1].upper() + if tl in "AT" and tr in "AT": + return 1 + if tl in "CG" and tr in "CG": + return 2 + i = 2 + while i <= min(len(probe_a), len(probe_b)) and not ( + _IUPAC.get(probe_a[-i].upper(), 0) & _IUPAC.get(probe_b[-i].upper(), 0) + ): + i += 1 + if i > min(len(probe_a), len(probe_b)): + return 255 + t = probe_a[-i].upper() + if t in "CGS": + return 1 + if t in "ATW": + return 2 + rules = { + ("Y", "G"): 1, + ("Y", "T"): 1, + ("Y", "A"): 2, + ("K", "C"): 1, + ("K", "A"): 2, + ("M", "G"): 1, + ("M", "T"): 2, + ("R", "C"): 1, + ("R", "T"): 2, + } + return rules.get((t, tr), 255) diff --git a/src/simurg/formats/csv_manifest.py b/src/simurg/formats/csv_manifest.py new file mode 100644 index 0000000..d5cdcb8 --- /dev/null +++ b/src/simurg/formats/csv_manifest.py @@ -0,0 +1,207 @@ +"""CSV manifest reader — plan §6.3, spec IN-041/IN-042. + +Sections: ``[Heading]`` (key,value), ``[Assay]`` (header row + one row per locus), ``[Controls]`` +(``addr[:addr..],Category,Color,Name``). The CSV carries ``SourceSeq``/``TopGenomicSeq`` (needed for indels) +but **no normalization IDs** — genotyping always needs the BPM. +""" + +from __future__ import annotations + +import csv +import io +from pathlib import Path +from typing import Any + +import numpy as np + +from simurg.core.errors import FormatError +from simurg.formats.bpm import Manifest, infer_assay_type + +REQUIRED_COLUMNS = ("IlmnID", "Name", "SNP", "AddressA_ID", "Chr", "MapInfo") +_SEQUENCE_COLUMNS = ("SourceSeq", "TopGenomicSeq", "AlleleA_ProbeSeq", "AlleleB_ProbeSeq") + + +def _open_text(path: Path) -> io.TextIOBase: + raw = path.open("rb") + head = raw.read(2) + raw.seek(0) + if head == b"\x1f\x8b": + import gzip + + return io.TextIOWrapper(gzip.GzipFile(fileobj=raw), encoding="utf-8-sig", errors="replace", newline="") + return io.TextIOWrapper(raw, encoding="utf-8-sig", errors="replace", newline="") + + +def _int_or_zero(s: str) -> int: + s = s.strip() + if not s: + return 0 + try: + return int(s) + except ValueError: + return int(float(s)) + + +def read_csv_manifest(path: str | Path, *, load_sequences: bool = True) -> Manifest: + """Parse a CSV manifest into the same columnar :class:`Manifest` as :func:`read_bpm`. + + ``load_sequences=False`` skips probe/source sequences to save memory. Raises E_CSV_MANIFEST_FORMAT. + """ + p = Path(path) + spath = str(p) + + def fail(msg: str) -> FormatError: + return FormatError("E_CSV_MANIFEST_FORMAT", f"{spath}: {msg}", path=spath) + + heading: dict[str, str] = {} + controls: list[dict[str, str]] = [] + cols: dict[str, list[str]] = {} + header: list[str] | None = None + section = None + intern: dict[str, str] = {} + + def it(s: str) -> str: + v = intern.get(s) + if v is None: + intern[s] = s + v = s + return v + + with _open_text(p) as fh: + for row in csv.reader(fh): + if not row or all(not c.strip() for c in row): + continue + first = row[0].strip() + if first.startswith("[") and first.endswith("]"): + section = first[1:-1].strip().lower() + continue + if section == "heading" and len(row) >= 2: + heading[first.rstrip()] = row[1].strip() + elif section == "heading": + continue + elif section == "assay": + if header is None: + header = [c.strip() for c in row] + missing = [c for c in REQUIRED_COLUMNS if c not in header] + if missing: + raise fail(f"[Assay] section lacks required columns: {', '.join(missing)}") + keep = [c for c in header if load_sequences or c not in _SEQUENCE_COLUMNS] + cols = {c: [] for c in keep} + idx = [(c, header.index(c)) for c in keep] + continue + if len(row) < len(header): + row = row + [""] * (len(header) - len(row)) + for c, i in idx: + v = row[i] + cols[c].append( + v + if c in ("IlmnID", "Name", "SourceSeq", "TopGenomicSeq", "AlleleA_ProbeSeq", "AlleleB_ProbeSeq") + else it(v) + ) + elif section == "controls" and len(row) >= 4: + controls.append( + { + "address": first, + "category": row[1].strip(), + "color": row[2].strip(), + "name": ",".join(row[3:]).strip(), + } + ) + if header is None: + raise fail("no [Assay] section") + n = len(cols["Name"]) + if n == 0: + raise fail("[Assay] section has no rows") + + def col(name: str) -> list[str]: + return cols.get(name) or [""] * n + + try: + addr_a = np.array([_int_or_zero(v) for v in cols["AddressA_ID"]], dtype=np.int32) + addr_b = np.array([_int_or_zero(v) for v in col("AddressB_ID")], dtype=np.int32) + map_info = np.array([_int_or_zero(v) for v in cols["MapInfo"]], dtype=np.int64) + except ValueError as e: + raise fail(f"non-numeric address or MapInfo value ({e})") from None + probe_a, probe_b, src = col("AlleleA_ProbeSeq"), col("AlleleB_ProbeSeq"), col("SourceSeq") + assay = np.zeros(n, dtype=np.uint8) + for k in np.flatnonzero(addr_b).tolist(): + assay[k] = infer_assay_type(probe_a[k], probe_b[k], src[k]) if load_sequences else 255 + exp = np.array([_int_or_zero(v) for v in col("Exp_Clusters")], dtype=np.int8) + int_only = np.array([_int_or_zero(v) != 0 for v in col("Intensity_Only")], dtype=np.bool_) + + return Manifest( + path=spath, + format="csv", + version=0, + manifest_name=heading.get("Descriptor File Name", p.name), + control_config="", + names=cols["Name"], + ilmn_ids=cols["IlmnID"], + snps=cols["SNP"], + chroms=cols["Chr"], + map_infos=map_info, + ilmn_strands=col("IlmnStrand"), + source_strands=col("SourceStrand"), + ref_strands=col("RefStrand"), + ploidies=col("Ploidy"), + species=col("Species"), + genome_builds=col("GenomeBuild"), + sources=col("Source"), + source_versions=col("SourceVersion"), + customer_strands=[""] * n, + source_seqs=src, + top_genomic_seqs=col("TopGenomicSeq"), + allele_a_probe_seqs=probe_a, + allele_b_probe_seqs=probe_b, + address_a=addr_a, + address_b=addr_b, + assay_types=assay, + exp_clusters=exp, + intensity_only=int_only, + norm_ids=np.zeros(n, dtype=np.uint8), + entry_versions=np.zeros(n, dtype=np.int8), + controls=controls, + heading=heading, + ) + + +def assay_format(manifest: Manifest) -> str: + """Platform code from the ``Assay Format`` heading (e.g. ``Infinium EX`` → ``EX``; ``Infinium LCG`` → ``LCG``).""" + v = manifest.heading.get("Assay Format", "") + return v.split()[-1] if v else "" + + +def compare_manifests(bpm: Manifest, csvm: Manifest, max_examples: int = 10) -> list[str]: + """BPM ↔ CSV consistency (spec IN-042). Returns human-readable differences (empty list = consistent).""" + problems: list[str] = [] + if bpm.num_loci != csvm.num_loci: + problems.append(f"locus count differs: BPM {bpm.num_loci}, CSV {csvm.num_loci}") + bidx = bpm.name_index() + cidx = csvm.name_index() + only_b = [n for n in bpm.names if n not in cidx] + only_c = [n for n in csvm.names if n not in bidx] + if only_b: + problems.append(f"{len(only_b)} loci only in BPM, e.g. {', '.join(only_b[:max_examples])}") + if only_c: + problems.append(f"{len(only_c)} loci only in CSV, e.g. {', '.join(only_c[:max_examples])}") + fields: list[tuple[str, Any, Any]] = [ + ("Chr", bpm.chroms, csvm.chroms), + ("MapInfo", bpm.map_infos, csvm.map_infos), + ("SNP", bpm.snps, csvm.snps), + ("AddressA", bpm.address_a, csvm.address_a), + ("AddressB", bpm.address_b, csvm.address_b), + ] + for label, bcol, ccol in fields: + bad: list[str] = [] + count = 0 + for name, i in bidx.items(): + j = cidx.get(name) + if j is None: + continue + if bcol[i] != ccol[j]: + count += 1 + if len(bad) < max_examples: + bad.append(f"{name} ({bcol[i]} vs {ccol[j]})") + if count: + problems.append(f"{label} differs for {count} loci, e.g. {'; '.join(bad)}") + return problems diff --git a/src/simurg/formats/detect.py b/src/simurg/formats/detect.py new file mode 100644 index 0000000..eec47f1 --- /dev/null +++ b/src/simurg/formats/detect.py @@ -0,0 +1,42 @@ +"""File-type sniffing for ``simurg inspect`` (spec CLI-020).""" + +from __future__ import annotations + +import gzip +from pathlib import Path + +KINDS = ("idat", "bpm", "egt", "gtc", "csv_manifest", "samplesheet", "vcf", "unknown") + + +def _head(path: Path, n: int = 4096) -> bytes: + with path.open("rb") as fh: + b = fh.read(n) + if b[:2] == b"\x1f\x8b": + try: + with gzip.open(path, "rb") as gz: + return gz.read(n) + except OSError: + return b + return b + + +def detect_kind(path: str | Path) -> str: + p = Path(path) + head = _head(p) + name = p.name.lower() + if head[:4] == b"IDAT": + return "idat" + if head[:3] == b"BPM": + return "bpm" + if head[:3] == b"gtc": + return "gtc" + if head.startswith(b"##fileformat=VCF"): + return "vcf" + if name.endswith(".egt") or (len(head) >= 5 and head[:4] == b"\x03\x00\x00\x00"): + return "egt" + text = head.decode("utf-8", errors="replace").lower() + if "[assay]" in text or ("[heading]" in text and "descriptor file name" in text): + return "csv_manifest" + if "sentrixbarcode_a" in text or "beadchipname" in text or "[data]" in text: + return "samplesheet" + return "unknown" diff --git a/src/simurg/formats/discovery.py b/src/simurg/formats/discovery.py new file mode 100644 index 0000000..ff6e3f1 --- /dev/null +++ b/src/simurg/formats/discovery.py @@ -0,0 +1,196 @@ +"""IDAT / GTC discovery and sample-set resolution — spec IN-016, IN-020 … IN-023, IN-030, CLI-006. + +Directory listing order never influences results: every path list is sorted bytewise (``os.fsencode``). +""" + +from __future__ import annotations + +import os +from collections import defaultdict +from dataclasses import dataclass, field +from pathlib import Path + +from simurg.core.errors import SimurgError +from simurg.formats.idat import parse_idat_filename, parse_run_time, read_idat +from simurg.formats.samplesheet import SampleSheet + + +@dataclass +class IdatPair: + barcode: str + position: str + grn: Path | None + red: Path | None + + @property + def sample_id(self) -> str: + return f"{self.barcode}_{self.position}" + + @property + def complete(self) -> bool: + return self.grn is not None and self.red is not None + + @property + def folder(self) -> Path: + p = self.grn or self.red + assert p is not None + return p.parent + + +@dataclass +class Discovery: + root: Path + pairs: dict[str, IdatPair] # sample_id -> pair (complete or not) + warnings: list[tuple[str, str]] = field(default_factory=list) + + +def _bytewise(paths: list[Path]) -> list[Path]: + return sorted(paths, key=lambda p: os.fsencode(str(p))) + + +def walk_files(root: Path, suffixes: tuple[str, ...]) -> list[Path]: + """Recursive listing following symlinks with cycle detection; sorted bytewise.""" + found: list[Path] = [] + seen_dirs: set[tuple[int, int]] = set() + stack = [root] + while stack: + d = stack.pop() + try: + st = d.stat() + except OSError: + continue + key = (st.st_dev, st.st_ino) + if key in seen_dirs: + continue + seen_dirs.add(key) + try: + entries = list(os.scandir(d)) + except OSError: + continue + for e in entries: + p = Path(e.path) + try: + if e.is_dir(follow_symlinks=True): + stack.append(p) + elif e.is_file(follow_symlinks=True) and e.name.lower().endswith(suffixes): + found.append(p) + except OSError: + continue + return _bytewise(found) + + +def _scan_time_key(p: Path) -> str: + try: + return parse_run_time(read_idat(p, load_arrays=False).scan_time_raw) + except SimurgError: + return "" + + +def discover_idats(root: str | Path, *, on_duplicate: str = "error") -> Discovery: + """Find ``__{Grn,Red}.idat[.gz]`` under ``root`` and pair them (IN-020 … IN-022).""" + rootp = Path(root) + if not rootp.is_dir(): + raise SimurgError("E_FILE_NOT_FOUND", f"IDAT folder not found: {rootp}", path=str(rootp)) + by_key: dict[tuple[str, str, str], list[Path]] = defaultdict(list) + for p in walk_files(rootp, (".idat", ".idat.gz")): + parsed = parse_idat_filename(p) + if parsed is None: + continue + by_key[parsed].append(p) + + warnings: list[tuple[str, str]] = [] + # group channel files by folder so that "the same sample in two folders" is a duplicate pair + per_sample: dict[tuple[str, str], dict[Path, dict[str, Path]]] = defaultdict(lambda: defaultdict(dict)) + for (barcode, pos, ch), paths in by_key.items(): + for p in paths: + per_sample[(barcode, pos)][p.parent][ch] = p + + pairs: dict[str, IdatPair] = {} + duplicates: dict[str, list[Path]] = {} + for (barcode, pos), folders in sorted(per_sample.items()): + candidates = [ + IdatPair(barcode, pos, chans.get("Grn"), chans.get("Red")) + for _, chans in sorted(folders.items(), key=lambda kv: os.fsencode(str(kv[0]))) + ] + sid = f"{barcode}_{pos}" + if len(candidates) > 1: + if on_duplicate == "error": + duplicates[sid] = [p for c in candidates for p in (c.grn, c.red) if p is not None] + continue + if on_duplicate == "newest": + # candidates are sorted bytewise; max() keeps the first (smallest path) among equal scan times + chosen = max(candidates, key=lambda c: _scan_time_key(c.grn or c.red)) # type: ignore[arg-type] + else: # "first": bytewise-smallest path + chosen = candidates[0] + warnings.append( + ( + "W_DUP_IDAT_RESOLVED", + f"{sid}: {len(candidates)} IDAT pairs found; using {chosen.folder} ({on_duplicate})", + ) + ) + else: + chosen = candidates[0] + if not chosen.complete: + have = "Grn" if chosen.grn else "Red" + warnings.append(("W_IDAT_UNPAIRED", f"{sid}: only the {have} IDAT was found in {chosen.folder}")) + pairs[sid] = chosen + if duplicates: + lines = [f"{sid}: " + ", ".join(str(p) for p in ps) for sid, ps in sorted(duplicates.items())] + raise SimurgError( + "E_DUP_IDAT", + f"{len(duplicates)} sample(s) have more than one IDAT pair:\n " + "\n ".join(lines), + hint="remove the extra copies or pass --on-duplicate-idat first|newest", + ) + return Discovery(rootp, pairs, warnings) + + +@dataclass +class ResolvedSample: + sample_id: str + sample_name: str # provisional: sheet Sample_Name or sample_id; IDAT field 406 is applied on read (IN-010) + pair: IdatPair | None + sheet_row: object | None + status: str # "READY" or "MISSING_INPUT" + + +def resolve_input_folder( + cli_folder: str | Path | None, sheet: SampleSheet | None +) -> tuple[Path, list[tuple[str, str]]]: + """Apply the input-set rule (CLI-006) and RootFolder precedence (IN-008).""" + notes: list[tuple[str, str]] = [] + if cli_folder is None and sheet is None: + raise SimurgError("E_NO_INPUT", "no input: give an IDAT folder, a sample sheet, or both") + if cli_folder is not None: + folder = Path(cli_folder) + if sheet is not None and sheet.root_folder is not None and sheet.root_folder.resolve() != folder.resolve(): + notes.append( + ("W_SHEET_ROOTFOLDER_OVERRIDDEN", f"CLI folder {folder} overrides sheet RootFolder {sheet.root_folder}") + ) + return folder, notes + assert sheet is not None + if sheet.root_folder is None: + raise SimurgError("E_NO_INPUT_FOLDER", f"{sheet.path}: sample sheet has no RootFolder and no folder was given") + return sheet.root_folder, notes + + +def resolve_samples(disc: Discovery, sheet: SampleSheet | None) -> tuple[list[ResolvedSample], list[tuple[str, str]]]: + """Match discovered IDAT pairs with sheet rows (IN-010, IN-016). Returns samples in a deterministic order.""" + notes: list[tuple[str, str]] = [] + out: list[ResolvedSample] = [] + if sheet is None: + for sid in sorted(disc.pairs): + found = disc.pairs[sid] + status = "READY" if found.complete else "MISSING_INPUT" + out.append(ResolvedSample(sid, sid, found, None, status)) + return out, notes + listed = set() + for row in sheet.rows: + listed.add(row.sample_id) + pair = disc.pairs.get(row.sample_id) + status = "READY" if pair is not None and pair.complete else "MISSING_INPUT" + name = row.sample_name or row.sample_id + out.append(ResolvedSample(row.sample_id, name, pair, row, status)) + extra = sorted(set(disc.pairs) - listed) + if extra: + notes.append(("I_UNLISTED_IDATS", f"{len(extra)} IDAT pair(s) not listed in the sample sheet are ignored")) + return out, notes diff --git a/src/simurg/formats/egt.py b/src/simurg/formats/egt.py new file mode 100644 index 0000000..66352c1 --- /dev/null +++ b/src/simurg/formats/egt.py @@ -0,0 +1,235 @@ +"""EGT cluster-file reader — plan §6.4, spec IN-044. + +Layout (little-endian, .NET strings):: + + int32 version (=3); string gencall_version, cluster_version, call_version, normalization_version, date_created + uint8 is_wgt (=1); string manifest_name; int32 data_block_version; string opa; int32 num_records + num_records × ClusterRecord (fixed size: 3 int32 + 12 float32 + extra floats) + num_records × ClusterScore (float sep, float total(GenTrain), float original, uint8 edited) + num_records × string genotype (unused) + num_records × string locus_name (join to the manifest BY NAME) + num_records × int32 address + num_records × int32[3] cluster counts (must equal the N in each record) + [data_block_version 9: num_records × float, meaning unknown] + +The number of extra floats per record differs between open readers (BeadArrayFiles: 2 for dbv<5, 15 for +dbv≥6; gtc2vcf: 0 for dbv 5, 15 for dbv≥7), and the intensity threshold is read for dbv 9 only +(BeadArrayFiles) or for dbv≥7 (gtc2vcf) — plan §20 U3. We try the candidate record sizes in order and accept +the first whose trailing cluster counts match the per-record N (a strong structural check). +""" + +from __future__ import annotations + +from dataclasses import dataclass +from pathlib import Path +from typing import Any + +import numpy as np +import numpy.typing as npt + +from simurg.core.errors import FormatError +from simurg.formats.binutil import BinaryReader, load_bytes + +_BASE_FIELDS = [ + ("n", " np.dtype[Any]: + fields: list[Any] = list(_BASE_FIELDS) + if extra: + fields.append(("extra", " list[int]: + if dbv < 5: + return [2, 0] + if dbv == 5: + return [0, 2, 15] + return [15, 0, 2] + + +@dataclass +class ClusterFile: + """Columnar EGT. Arrays with a trailing axis of 3 are ordered (AA, AB, BB).""" + + path: str + version: int + gencall_version: str + cluster_version: str + call_version: str + normalization_version: str + date_created: str + manifest_name: str + data_block_version: int + opa: str + names: list[str] + addresses: npt.NDArray[np.int32] + n: npt.NDArray[np.int32] # (records, 3) + theta_mean: npt.NDArray[np.float32] + theta_dev: npt.NDArray[np.float32] + r_mean: npt.NDArray[np.float32] + r_dev: npt.NDArray[np.float32] + intensity_threshold: npt.NDArray[np.float32] # NaN where not stored + cluster_sep: npt.NDArray[np.float32] + gentrain_score: npt.NDArray[np.float32] + original_score: npt.NDArray[np.float32] + edited: npt.NDArray[np.bool_] + extra_floats: npt.NDArray[np.float32] | None + trailer_floats: npt.NDArray[np.float32] | None + extra_per_record: int + + @property + def num_records(self) -> int: + return len(self.names) + + def name_index(self) -> dict[str, int]: + return {nm: i for i, nm in enumerate(self.names)} + + def mean_r(self) -> npt.NDArray[np.float64]: + """Count-weighted mean R — the expected R of intensity-only loci (plan §6.4).""" + w = self.n.astype(np.float64) + tot = w.sum(axis=1) + with np.errstate(invalid="ignore", divide="ignore"): + return np.asarray((w * self.r_mean).sum(axis=1) / tot, dtype=np.float64) + + def summary(self) -> dict[str, Any]: + thr = self.intensity_threshold + return { + "path": self.path, + "version": self.version, + "data_block_version": self.data_block_version, + "extra_floats_per_record": self.extra_per_record, + "num_records": self.num_records, + "manifest_name": self.manifest_name, + "gencall_version": self.gencall_version, + "cluster_version": self.cluster_version, + "call_version": self.call_version, + "normalization_version": self.normalization_version, + "date_created": self.date_created, + "intensity_threshold_nonzero": int(np.count_nonzero(np.nan_to_num(thr))), + "edited_records": int(self.edited.sum()), + "zero_gentrain": int((self.gentrain_score == 0).sum()), + } + + +def read_egt(path: str | Path) -> ClusterFile: + """Parse an EGT cluster file. Raises FormatError(E_EGT_LOAD).""" + spath = str(path) + buf = load_bytes(path) + r = BinaryReader(buf, "E_EGT_LOAD", spath) + version = r.i32("version") + if version != 3: + raise r.fail(f"unsupported cluster file version {version}") + gencall_version = r.string() + cluster_version = r.string() + call_version = r.string() + normalization_version = r.string() + date_created = r.string() + if (is_wgt := r.u8("is_wgt")) != 1: + raise r.fail(f"only WGT cluster files are supported (is_wgt={is_wgt})") + manifest_name = r.string() + dbv = r.i32("data block version") + if not 1 <= dbv <= 9: + raise r.fail(f"unsupported data block version {dbv}") + opa = r.string() + nrec = r.i32("number of records") + if nrec <= 0 or nrec > 50_000_000: + raise r.fail(f"implausible number of records {nrec}") + records_start = r.pos + + last_error: FormatError | None = None + for extra in _extra_candidates(dbv): + try: + return _parse_body( + r, + records_start, + nrec, + extra, + dbv, + spath, + version, + gencall_version, + cluster_version, + call_version, + normalization_version, + date_created, + manifest_name, + opa, + ) + except FormatError as e: + last_error = e + assert last_error is not None + raise last_error + + +def _parse_body( + r: BinaryReader, + start: int, + nrec: int, + extra: int, + dbv: int, + spath: str, + version: int, + gencall_version: str, + cluster_version: str, + call_version: str, + normalization_version: str, + date_created: str, + manifest_name: str, + opa: str, +) -> ClusterFile: + r.seek(start) + rec = r.array(_record_dtype(extra), nrec, "cluster records") + scores = r.array(_SCORE_DTYPE, nrec, "cluster scores") + for _ in range(nrec): + r.skip_string("genotype") + names = [r.string("locus name") for _ in range(nrec)] + addresses = np.array(r.array("= 4 * nrec: + trailer = np.array(r.array("= 7 and extra == 15: + thr = extra_arr[:, 0].copy() + else: + thr = np.full(nrec, np.nan, dtype=np.float32) + return ClusterFile( + path=spath, + version=version, + gencall_version=gencall_version, + cluster_version=cluster_version, + call_version=call_version, + normalization_version=normalization_version, + date_created=date_created, + manifest_name=manifest_name, + data_block_version=dbv, + opa=opa, + names=names, + addresses=addresses, + n=n, + theta_mean=np.array(rec["theta_mean"], dtype=np.float32), + theta_dev=np.array(rec["theta_dev"], dtype=np.float32), + r_mean=np.array(rec["r_mean"], dtype=np.float32), + r_dev=np.array(rec["r_dev"], dtype=np.float32), + intensity_threshold=thr, + cluster_sep=np.array(scores["cluster_sep"], dtype=np.float32), + gentrain_score=np.array(scores["total_score"], dtype=np.float32), + original_score=np.array(scores["original_score"], dtype=np.float32), + edited=np.array(scores["edited"] != 0), + extra_floats=extra_arr, + trailer_floats=trailer, + extra_per_record=extra, + ) diff --git a/src/simurg/formats/gtc.py b/src/simurg/formats/gtc.py new file mode 100644 index 0000000..0391b17 --- /dev/null +++ b/src/simurg/formats/gtc.py @@ -0,0 +1,548 @@ +"""GTC genotype file (versions 3–5) reader and writer — plan §6.5, spec IN-031/IN-032/IN-033, OUT-001. + +Layout:: + + char[3] "gtc"; uint8 version; int32 M; M × {uint16 id; uint32 offset} + (for ids 1, 2, 3 the "offset" IS the value) + +Authoritative description: ``docs/GTC_File_Format_v5.pdf`` in Illumina BeadArrayFiles (BSD-2); canonical write +order follows gtc2vcf ``gtc_write`` (MIT). The reader records the TOC order and the on-disk data order, so a +file can be re-written byte-for-byte (P1 acceptance: round-trip). +""" + +from __future__ import annotations + +import math +from dataclasses import dataclass, field +from pathlib import Path +from typing import Any + +import numpy as np +import numpy.typing as npt + +from simurg.formats.binutil import BinaryReader, BinaryWriter, load_bytes + +NUM_SNPS = 1 +PLOIDY = 2 +PLOIDY_TYPE = 3 +SAMPLE_NAME = 10 +SAMPLE_PLATE = 11 +SAMPLE_WELL = 12 +CLUSTER_FILE = 100 +SNP_MANIFEST = 101 +IMAGING_DATE = 200 +AUTOCALL_DATE = 201 +AUTOCALL_VERSION = 300 +NORMALIZATION_TRANSFORMS = 400 +CONTROLS_X = 500 +CONTROLS_Y = 501 +RAW_X = 1000 +RAW_Y = 1001 +GENOTYPES = 1002 +BASE_CALLS = 1003 +GENOTYPE_SCORES = 1004 +SCANNER_DATA = 1005 +CALL_RATE = 1006 +GENDER = 1007 +LOGR_DEV = 1008 +GC10 = 1009 +DX = 1010 +SAMPLE_DATA = 1011 +B_ALLELE_FREQS = 1012 +LOGR_RATIOS = 1013 +PERCENTILES_X = 1014 +PERCENTILES_Y = 1015 +SLIDE_IDENTIFIER = 1016 + +INLINE_IDS = (NUM_SNPS, PLOIDY, PLOIDY_TYPE) + +TOC_NAMES: dict[int, str] = { + NUM_SNPS: "num_snps", + PLOIDY: "ploidy", + PLOIDY_TYPE: "ploidy_type", + SAMPLE_NAME: "sample_name", + SAMPLE_PLATE: "sample_plate", + SAMPLE_WELL: "sample_well", + CLUSTER_FILE: "cluster_file", + SNP_MANIFEST: "snp_manifest", + IMAGING_DATE: "imaging_date", + AUTOCALL_DATE: "autocall_date", + AUTOCALL_VERSION: "autocall_version", + NORMALIZATION_TRANSFORMS: "normalization_transforms", + CONTROLS_X: "controls_x", + CONTROLS_Y: "controls_y", + RAW_X: "raw_x", + RAW_Y: "raw_y", + GENOTYPES: "genotypes", + BASE_CALLS: "base_calls", + GENOTYPE_SCORES: "genotype_scores", + SCANNER_DATA: "scanner_data", + CALL_RATE: "call_rate", + GENDER: "gender", + LOGR_DEV: "logr_dev", + GC10: "gc10", + DX: "dx", + SAMPLE_DATA: "sample_data", + B_ALLELE_FREQS: "b_allele_freqs", + LOGR_RATIOS: "logr_ratios", + PERCENTILES_X: "percentiles_x", + PERCENTILES_Y: "percentiles_y", + SLIDE_IDENTIFIER: "slide_identifier", +} + +# canonical order (gtc2vcf gtc_write); v3 omits PLOIDY, PLOIDY_TYPE and 1012..1016 +CANONICAL_ORDER_V5 = [ + NUM_SNPS, + PLOIDY, + PLOIDY_TYPE, + SAMPLE_NAME, + SAMPLE_PLATE, + SAMPLE_WELL, + CLUSTER_FILE, + SNP_MANIFEST, + IMAGING_DATE, + AUTOCALL_DATE, + AUTOCALL_VERSION, + NORMALIZATION_TRANSFORMS, + CONTROLS_X, + CONTROLS_Y, + RAW_X, + RAW_Y, + GENOTYPES, + BASE_CALLS, + GENOTYPE_SCORES, + SCANNER_DATA, + CALL_RATE, + GENDER, + LOGR_DEV, + GC10, + DX, + SAMPLE_DATA, + B_ALLELE_FREQS, + LOGR_RATIOS, + PERCENTILES_X, + PERCENTILES_Y, + SLIDE_IDENTIFIER, +] +_V4_PLUS = {PLOIDY, PLOIDY_TYPE, B_ALLELE_FREQS, LOGR_RATIOS, PERCENTILES_X, PERCENTILES_Y, SLIDE_IDENTIFIER} +CANONICAL_ORDER_V3 = [i for i in CANONICAL_ORDER_V5 if i not in _V4_PLUS] + +GENOTYPE_CODES = [ + "NC", + "AA", + "AB", + "BB", + "NULL", + "A", + "B", + "AAA", + "AAB", + "ABB", + "BBB", + "AAAA", + "AAAB", + "AABB", + "ABBB", + "BBBB", + "AAAAA", + "AAAAB", + "AAABB", + "AABBB", + "ABBBB", + "BBBBB", + "AAAAAA", + "AAAAAB", + "AAAABB", + "AAABBB", + "AABBBB", + "ABBBBB", + "BBBBBB", + "AAAAAAA", + "AAAAAAB", + "AAAAABB", + "AAAABBB", + "AAABBBB", + "AABBBBB", + "ABBBBBB", + "BBBBBBB", + "AAAAAAAA", + "AAAAAAAB", + "AAAAAABB", + "AAAAABBB", + "AAAABBBB", + "AAABBBBB", + "AABBBBBB", + "ABBBBBBB", + "BBBBBBBB", +] + +XFORM_DTYPE = np.dtype( + [ + ("version", " int: + return int.from_bytes(self.pmt_green, "little", signed=True) + + @property + def pmt_red_int(self) -> int: + return int.from_bytes(self.pmt_red, "little", signed=True) + + +@dataclass +class Gtc: + """In-memory GTC content. Arrays are per locus in BPM names order.""" + + version: int = 5 + num_snps: int = 0 + ploidy: int = 2 + ploidy_type: int = 1 + sample_name: str = "" + sample_plate: str = "" + sample_well: str = "" + cluster_file: str = "" + snp_manifest: str = "" + imaging_date: str = "" + autocall_date: str = "" + autocall_version: str = "" + normalization_transforms: npt.NDArray[Any] = field(default_factory=lambda: np.zeros(0, XFORM_DTYPE)) + controls_x: npt.NDArray[np.uint16] = field(default_factory=lambda: np.zeros(0, np.uint16)) + controls_y: npt.NDArray[np.uint16] = field(default_factory=lambda: np.zeros(0, np.uint16)) + raw_x: npt.NDArray[np.uint16] = field(default_factory=lambda: np.zeros(0, np.uint16)) + raw_y: npt.NDArray[np.uint16] = field(default_factory=lambda: np.zeros(0, np.uint16)) + genotypes: npt.NDArray[np.uint8] = field(default_factory=lambda: np.zeros(0, np.uint8)) + base_calls: npt.NDArray[np.uint8] = field(default_factory=lambda: np.zeros((0, 2), np.uint8)) + genotype_scores: npt.NDArray[np.float32] = field(default_factory=lambda: np.zeros(0, np.float32)) + scanner_data: ScannerData = field(default_factory=ScannerData) + call_rate: float = 0.0 + gender: bytes = b"U" + logr_dev: float = 0.0 + gc10: float = 0.0 + dx: int = 0 + sample_data: npt.NDArray[Any] = field(default_factory=lambda: np.zeros(1, SAMPLE_DATA_DTYPE)[0]) + b_allele_freqs: npt.NDArray[np.float32] = field(default_factory=lambda: np.zeros(0, np.float32)) + logr_ratios: npt.NDArray[np.float32] = field(default_factory=lambda: np.zeros(0, np.float32)) + percentiles_x: npt.NDArray[np.uint16] = field(default_factory=lambda: np.zeros(3, np.uint16)) + percentiles_y: npt.NDArray[np.uint16] = field(default_factory=lambda: np.zeros(3, np.uint16)) + slide_identifier: str = "" + # layout bookkeeping (set by the reader; None → canonical order) + toc_order: list[int] | None = None + data_order: list[int] | None = None + present: set[int] | None = None + path: str | None = None + + # -- convenience ------------------------------------------------------------------------------------ + @property + def gc50(self) -> float: + return float(self.sample_data["p50gc"]) + + @property + def num_calls(self) -> int: + return int(self.sample_data["num_calls"]) + + @property + def num_no_calls(self) -> int: + return int(self.sample_data["num_no_calls"]) + + @property + def num_intensity_only(self) -> int: + return int(self.sample_data["num_intensity_only"]) + + @property + def gender_char(self) -> str: + return self.gender.decode("latin-1") if self.gender else "" + + def genotype_strings(self) -> list[str]: + return [GENOTYPE_CODES[g] if g < len(GENOTYPE_CODES) else "?" for g in self.genotypes] + + def base_call_strings(self) -> list[str]: + """TOP-strand base calls; '-' for no-call (BeadArrayFiles ``get_base_calls`` for diploids).""" + return [bytes(bc).decode("latin-1") for bc in self.base_calls] + + def summary(self) -> dict[str, Any]: + def f(v: float) -> float | None: + return None if math.isnan(v) else round(float(v), 6) + + return { + "path": self.path, + "version": self.version, + "num_snps": self.num_snps, + "ploidy": self.ploidy, + "ploidy_type": self.ploidy_type, + "sample_name": self.sample_name, + "sample_plate": self.sample_plate, + "sample_well": self.sample_well, + "slide_identifier": self.slide_identifier, + "cluster_file": self.cluster_file, + "snp_manifest": self.snp_manifest, + "imaging_date": self.imaging_date, + "autocall_date": self.autocall_date, + "autocall_version": self.autocall_version, + "n_normalization_transforms": len(self.normalization_transforms), + "n_controls": len(self.controls_x), + "call_rate": f(self.call_rate), + "gender": self.gender_char, + "logr_dev": f(self.logr_dev), + "gc10": f(self.gc10), + "gc50": f(self.gc50), + "num_calls": self.num_calls, + "num_no_calls": self.num_no_calls, + "num_intensity_only": self.num_intensity_only, + "percentiles_x": [int(v) for v in self.percentiles_x], + "percentiles_y": [int(v) for v in self.percentiles_y], + "scanner": { + "name": self.scanner_data.name, + "version": self.scanner_data.version, + "user": self.scanner_data.user, + }, + "toc": [TOC_NAMES.get(i, str(i)) for i in (self.toc_order or [])], + } + + +def read_gtc(path: str | Path, *, check_write_complete: bool = True) -> Gtc: + """Parse a GTC file. Raises FormatError(E_GTC_CORRUPT) on any structural problem (IN-032).""" + spath = str(path) + buf = load_bytes(path) + r = BinaryReader(buf, "E_GTC_CORRUPT", spath) + if r.size < 8 or r.raw(3, "magic") != b"gtc": + raise r.fail("bad magic (not a GTC file)") + version = r.u8("version") + if version not in (3, 4, 5): + raise r.fail(f"unsupported GTC version {version}") + m = r.i32("TOC size") + if not 0 < m < 1024: + raise r.fail(f"implausible TOC size {m}") + toc: dict[int, int] = {} + order: list[int] = [] + for _ in range(m): + tid = r.u16("TOC id") + off = r.u32("TOC offset") + if tid in toc: + raise r.fail(f"duplicate TOC id {tid}") + toc[tid] = off + order.append(tid) + g = Gtc(version=version, toc_order=order, present=set(toc), path=spath) + for tid in INLINE_IDS: + if tid in toc: + setattr(g, TOC_NAMES[tid], _as_i32(toc[tid])) + if NUM_SNPS not in toc: + raise r.fail("missing NumSNPs entry") + n = g.num_snps + if n < 0: + raise r.fail(f"negative NumSNPs {n}") + + def at(tid: int) -> bool: + if tid not in toc: + return False + r.seek(toc[tid], TOC_NAMES.get(tid, str(tid))) + return True + + def arr(tid: int, dtype: Any, expect: int | None = None) -> npt.NDArray[Any]: + a = r.counted_array(dtype, TOC_NAMES[tid]) + if expect is not None and len(a) != expect: + raise r.fail(f"{TOC_NAMES[tid]} has {len(a)} entries, expected {expect}") + return np.array(a) + + for tid, attr in ( + (SAMPLE_NAME, "sample_name"), + (SAMPLE_PLATE, "sample_plate"), + (SAMPLE_WELL, "sample_well"), + (CLUSTER_FILE, "cluster_file"), + (SNP_MANIFEST, "snp_manifest"), + (IMAGING_DATE, "imaging_date"), + (AUTOCALL_DATE, "autocall_date"), + (AUTOCALL_VERSION, "autocall_version"), + (SLIDE_IDENTIFIER, "slide_identifier"), + ): + if at(tid): + setattr(g, attr, r.string(attr)) + if at(NORMALIZATION_TRANSFORMS): + g.normalization_transforms = arr(NORMALIZATION_TRANSFORMS, XFORM_DTYPE) + if at(CONTROLS_X): + g.controls_x = arr(CONTROLS_X, " int: + return v - (1 << 32) if v >= 1 << 31 else v + + +def _write_field(w: BinaryWriter, g: Gtc, tid: int) -> None: + n = g.num_snps + if tid in ( + SAMPLE_NAME, + SAMPLE_PLATE, + SAMPLE_WELL, + CLUSTER_FILE, + SNP_MANIFEST, + IMAGING_DATE, + AUTOCALL_DATE, + AUTOCALL_VERSION, + SLIDE_IDENTIFIER, + ): + w.string(getattr(g, TOC_NAMES[tid])) + elif tid == NORMALIZATION_TRANSFORMS: + xf = np.asarray(g.normalization_transforms, dtype=XFORM_DTYPE) + w.i32(len(xf)) + w.write(xf.tobytes()) + elif tid in (CONTROLS_X, CONTROLS_Y): + w.counted_array(getattr(g, TOC_NAMES[tid]), " None: + if len(getattr(g, TOC_NAMES[tid])) != n: + raise ValueError(f"{TOC_NAMES[tid]} length {len(getattr(g, TOC_NAMES[tid]))} != num_snps {n}") + + +def encode_gtc(g: Gtc, *, canonical: bool = False) -> bytes: + """Serialize. ``canonical=True`` ignores the recorded layout and uses the gtc2vcf order for ``g.version``.""" + if canonical or g.toc_order is None: + toc_order = list(CANONICAL_ORDER_V5 if g.version >= 4 else CANONICAL_ORDER_V3) + data_order = [t for t in toc_order if t not in INLINE_IDS] + else: + toc_order = list(g.toc_order) + data_order = list(g.data_order or [t for t in toc_order if t not in INLINE_IDS]) + header_len = 4 + 4 + 6 * len(toc_order) + body = BinaryWriter() + offsets: dict[int, int] = {} + for tid in data_order: + offsets[tid] = header_len + body.tell() + _write_field(body, g, tid) + w = BinaryWriter() + w.write(b"gtc") + w.u8(g.version) + w.i32(len(toc_order)) + for tid in toc_order: + w.u16(tid) + if tid in INLINE_IDS: + w.u32(getattr(g, TOC_NAMES[tid]) & 0xFFFFFFFF) + else: + w.u32(offsets[tid]) + w.write(body.getvalue()) + return w.getvalue() + + +def write_gtc(g: Gtc, path: str | Path, *, canonical: bool = False) -> None: + """Write atomically (spec CLI-011).""" + from simurg.core.atomic import atomic_write_bytes + + atomic_write_bytes(Path(path), encode_gtc(g, canonical=canonical)) + + +def validate_gtc(g: Gtc) -> list[str]: + """Light semantic checks used by ``simurg inspect`` (not raised).""" + problems: list[str] = [] + n = g.num_snps + for name in ("raw_x", "raw_y", "genotypes", "genotype_scores"): + if len(getattr(g, name)) not in (0, n): + problems.append(f"{name} length {len(getattr(g, name))} != num_snps {n}") + if len(g.genotypes) and int(g.genotypes.max(initial=0)) >= len(GENOTYPE_CODES): + problems.append("genotype code out of range") + return problems diff --git a/src/simurg/formats/idat.py b/src/simurg/formats/idat.py new file mode 100644 index 0000000..9ba8a72 --- /dev/null +++ b/src/simurg/formats/idat.py @@ -0,0 +1,343 @@ +"""IDAT (binary, version 3) reader — plan §6.1, spec IN-023/IN-024/IN-025. + +Layout (little-endian):: + + char[4] "IDAT"; int64 version (3); int32 n_fields; n_fields × {uint16 code; int64 offset} + +Field codes follow illuminaio (GPL, consulted as documentation only), gtc2vcf ``idat2gtc.c`` (MIT) and +BeadArrayFiles (BSD-2). The element type of the MidBlock (code 200) differs between those readers (int32 in +illuminaio, uint8 in gtc2vcf); we infer it from the space available before the next section (plan §20 U2). +""" + +from __future__ import annotations + +import datetime as _dt +import re +from dataclasses import dataclass, field +from pathlib import Path +from typing import Any + +import numpy as np +import numpy.typing as npt + +from simurg.core.errors import FormatError +from simurg.formats.binutil import BinaryReader, load_bytes + +MAGIC = b"IDAT" + +# field codes +N_SNPS = 1000 +ILLUMINA_ID = 102 +SD = 103 +MEAN = 104 +NBEADS = 107 +MID_BLOCK = 200 +RUN_INFO = 300 +RED_GREEN = 400 +MANIFEST = 401 +BARCODE = 402 +CHIP_TYPE = 403 +POSITION = 404 +BEADSET = 405 +SAMPLE_ID = 406 +DESCRIPTION = 407 +PLATE = 408 +WELL = 409 +SAMPLE_COUNT = 410 +VLN = 510 + +FIELD_NAMES: dict[int, str] = { + N_SNPS: "n_addresses", + ILLUMINA_ID: "address_ids", + SD: "sd", + MEAN: "mean", + NBEADS: "n_beads", + MID_BLOCK: "mid_block", + RUN_INFO: "run_info", + RED_GREEN: "red_green", + MANIFEST: "manifest", + BARCODE: "barcode", + CHIP_TYPE: "chip_type", + POSITION: "position", + BEADSET: "beadset", + SAMPLE_ID: "sample_id", + DESCRIPTION: "description", + PLATE: "plate", + WELL: "well", + SAMPLE_COUNT: "sample_count", + VLN: "vln", +} + +_STRING_FIELDS = { + MANIFEST: "manifest", + BARCODE: "barcode", + CHIP_TYPE: "chip_type", + POSITION: "position", + BEADSET: "beadset", + SAMPLE_ID: "sample_id", + DESCRIPTION: "description", + PLATE: "plate", + WELL: "well", + VLN: "vln", +} + +_NAME_RE = re.compile( + r"^(?P[^_/\\]+)_(?PR\d{2}C\d{2})_(?PGrn|Red)\.idat(\.gz)?$", re.IGNORECASE +) + +_DATE_FORMATS = ( + "%m/%d/%Y %I:%M:%S %p", + "%m/%d/%Y %H:%M:%S", + "%Y-%m-%dT%H:%M:%S", + "%Y-%m-%d %H:%M:%S", + "%d.%m.%Y %H:%M:%S", +) + + +@dataclass(frozen=True) +class RunInfoEntry: + run_time: str + block_type: str + block_pars: str + block_code: str + code_version: str + + +@dataclass +class Idat: + """One colour channel of one BeadChip section.""" + + path: str + version: int + toc: dict[int, int] + n_addresses: int + address_ids: npt.NDArray[np.int32] + mean: npt.NDArray[np.uint16] + sd: npt.NDArray[np.uint16] + n_beads: npt.NDArray[np.uint8] + mid_block: npt.NDArray[Any] | None = None + mid_block_itemsize: int | None = None + run_info: list[RunInfoEntry] = field(default_factory=list) + red_green: bytes = b"" + manifest: str = "" + barcode: str = "" + chip_type: str = "" + position: str = "" + beadset: str = "" + sample_id: str = "" + description: str = "" + plate: str = "" + well: str = "" + vln: str = "" + sample_count: int | None = None + unknown_fields: dict[int, int] = field(default_factory=dict) + + # -- derived metadata (IN-024) -------------------------------------------------------------------- + def scan_entry(self) -> RunInfoEntry | None: + scans = [r for r in self.run_info if r.block_type.strip().lower() == "scan"] + return scans[-1] if scans else None + + @property + def scan_time_raw(self) -> str: + e = self.scan_entry() + return e.run_time if e else "" + + @property + def scan_time_iso(self) -> str: + return parse_run_time(self.scan_time_raw) + + @property + def scanner_software(self) -> str: + e = self.scan_entry() + return e.code_version if e else "" + + @property + def scanner_pars(self) -> str: + e = self.scan_entry() + return e.block_pars if e else "" + + def metadata(self) -> dict[str, Any]: + return { + "path": self.path, + "version": self.version, + "n_addresses": self.n_addresses, + "barcode": self.barcode, + "position": self.position, + "chip_type": self.chip_type, + "manifest": self.manifest, + "beadset": self.beadset, + "sample_id": self.sample_id, + "description": self.description, + "plate": self.plate, + "well": self.well, + "sample_count": self.sample_count, + "vln": self.vln, + "red_green": self.red_green.hex(" "), + "mid_block_count": None if self.mid_block is None else len(self.mid_block), + "mid_block_itemsize": self.mid_block_itemsize, + "scan_time": self.scan_time_iso, + "scan_time_raw": self.scan_time_raw, + "scanner_software": self.scanner_software, + "toc_codes": sorted(self.toc), + "unknown_codes": sorted(self.unknown_fields), + "run_info": [r.__dict__ for r in self.run_info], + } + + +def parse_run_time(raw: str) -> str: + """Convert a RunInfo timestamp to ISO-8601 without consulting the process locale (EDGE-001). + + Returns ``""`` if the value is empty or cannot be parsed. + """ + s = raw.strip() + if not s: + return "" + for fmt in _DATE_FORMATS: + try: + return _dt.datetime.strptime(s, fmt).isoformat() + except ValueError: + continue + return "" + + +def parse_idat_filename(path: str | Path) -> tuple[str, str, str] | None: + """``__.idat[.gz]`` → (barcode, POSITION, 'Grn'|'Red'); None if not an IDAT name.""" + m = _NAME_RE.match(Path(path).name) + if not m: + return None + ch = "Grn" if m.group("channel").lower() == "grn" else "Red" + return m.group("barcode"), m.group("position").upper(), ch + + +def read_idat(path: str | Path, *, load_arrays: bool = True) -> Idat: + """Parse an IDAT file. Raises FormatError (E_IDAT_FORMAT / E_IDAT_ENCRYPTED).""" + spath = str(path) + buf = load_bytes(path) + r = BinaryReader(buf, "E_IDAT_FORMAT", spath) + if r.size < 16 or r.raw(4, "magic") != MAGIC: + raise r.fail("bad magic (not an IDAT file)") + version = r.i64("version") + # encrypted files carry version 1 in the low 32 bits followed by ciphertext (illuminaio reads int32) + if version & 0xFFFFFFFF == 1: + raise FormatError( + "E_IDAT_ENCRYPTED", + f"{spath}: encrypted (version 1) IDAT files are not supported", + path=spath, + hint="convert with an external tool first", + ) + if version != 3: + raise r.fail(f"unsupported IDAT version {version}") + n_fields = r.i32("TOC size") + if n_fields <= 0 or n_fields > 4096: + raise r.fail(f"implausible TOC size {n_fields}") + toc: dict[int, int] = {} + for _ in range(n_fields): + code = r.u16("TOC code") + off = r.i64("TOC offset") + if off < 0 or off >= r.size: + raise r.fail(f"TOC offset {off} for field {code} outside file (size {r.size})") + if code in toc: + raise r.fail(f"duplicate TOC field {code}") + toc[code] = off + + for required in (N_SNPS, ILLUMINA_ID, MEAN, SD, NBEADS): + if required not in toc: + raise r.fail(f"missing required field {required} ({FIELD_NAMES[required]})") + + r.seek(toc[N_SNPS]) + n = r.i32("number of addresses") + if n < 0: + raise r.fail(f"negative address count {n}") + + empty_i32 = np.empty(0, dtype=" int: + return next(o for o in offsets_sorted if o > off) + + for code, off in sorted(toc.items(), key=lambda kv: kv[1]): + r.seek(off) + if code in _STRING_FIELDS: + setattr(idat, _STRING_FIELDS[code], r.string(FIELD_NAMES[code])) + elif code == RUN_INFO: + m = r.i32("run info count") + if m < 0 or m > 100_000: + raise r.fail(f"implausible run info count {m}") + idat.run_info = [RunInfoEntry(*(r.string("run info") for _ in range(5))) for _ in range(m)] + elif code == RED_GREEN: + idat.red_green = r.raw(4, "RedGreen") + elif code == SAMPLE_COUNT: + idat.sample_count = r.i32("sample count") + elif code == MID_BLOCK: + if not load_arrays: + continue + count = r.i32("mid block count") + space = section_end(off) - off - 4 + if count < 0: + raise r.fail(f"negative mid block count {count}") + # element size inferred from the space before the next section (int32 vs uint8, plan §20 U2) + if count == space and count * 4 != space: + itemsize = 1 + elif count * 4 <= space: + itemsize = 4 + elif count <= space: + itemsize = 1 + else: + raise r.fail(f"mid block of {count} entries does not fit in {space} bytes") + idat.mid_block_itemsize = itemsize + idat.mid_block = r.array(" bool: + """Duplicate check; O(n) when IDs are already strictly increasing (the usual case).""" + if len(ids) < 2: + return False + d = np.diff(ids) + if (d > 0).all(): + return False + s = np.sort(ids) + return bool((s[1:] == s[:-1]).any()) + + +def channel_of(idat: Idat) -> str | None: + parsed = parse_idat_filename(idat.path) + return parsed[2] if parsed else None + + +def check_pair(grn: Idat, red: Idat) -> None: + """Red/Grn consistency (IN-023): same N and the same address set. Raises E_IDAT_PAIR_MISMATCH.""" + if grn.n_addresses != red.n_addresses: + raise FormatError( + "E_IDAT_PAIR_MISMATCH", + f"Grn has {grn.n_addresses} addresses, Red has {red.n_addresses}: {grn.path}, {red.path}", + ) + if len(grn.address_ids) and not np.array_equal(np.sort(grn.address_ids), np.sort(red.address_ids)): + raise FormatError("E_IDAT_PAIR_MISMATCH", f"address sets differ: {grn.path}, {red.path}") diff --git a/src/simurg/formats/samplesheet.py b/src/simurg/formats/samplesheet.py new file mode 100644 index 0000000..15cbe91 --- /dev/null +++ b/src/simurg/formats/samplesheet.py @@ -0,0 +1,359 @@ +"""Sample-sheet parser and validator — spec IN-001 … IN-017 (plan §6.6). + +A GenomeStudio-compatible CSV with optional ``[Header]``, ``[Manifests]`` (ignored) and ``[Data]`` sections. +All column matching uses Unicode case folding (locale independent, EDGE-001). +""" + +from __future__ import annotations + +import csv +import io +import re +from dataclasses import dataclass, field +from pathlib import Path + +from simurg.core.errors import SimurgError + +MAX_LEN = 500 +KEY_PAIRS = (("SentrixBarcode_A", "SentrixPosition_A"), ("beadChipName", "sampleSectionName")) +RECOGNISED = ( + "Sample_ID", + "Sample_Name", + "Sample_Plate", + "Sample_Well", + "Sample_Group", + "Gender", + "Sex", + "demographicSex", + "Parent1", + "Parent2", + "Replicate", +) +SEX_COLUMNS_PRIORITY = ("Sex", "Gender", "demographicSex") + +# Output columns of sample_summary (OUT-060) and their compatibility spellings (OUT-066): a sheet column may +# not collide with these case-insensitively (IN-005). Recognised input columns above are exempt. +RESERVED_OUTPUT_COLUMNS = ( + "sample_id", + "sample_name", + "sample_folder", + "status", + "fail_reason", + "autosomal_call_rate", + "call_rate", + "logr_dev", + "sex_est", + "sex_reported", + "sex_check", + "gc10", + "gc50", + "n_calls", + "n_nocalls", + "n_intensity_only", + "p05_x", + "p50_x", + "p95_x", + "p05_y", + "p50_y", + "p95_y", + "het_rate", + "contamination", + "functional_status", + "control_status", + "control_flag_codes", + "flags", + "Sample ID", + "Sample Name", + "Sample Folder", + "Autosomal Call Rate", + "Call Rate", + "Log R Ratio Std Dev", + "Sex Estimate", +) + +_BARCODE_RE = re.compile(r"^[0-9]{10,14}$") +_BARCODE_DECIMAL_RE = re.compile(r"^([0-9]{10,14})\.0+$") +_BARCODE_SCI_RE = re.compile(r"^[0-9][.,][0-9]+[eE][+]?[0-9]+$") +_POSITION_RE = re.compile(r"^[Rr][0-9]{2}[Cc][0-9]{2}$") +_NORM_RE = re.compile(r".*_norm_r$", re.IGNORECASE) +_NORM_COMPAT_RE = re.compile(r".* norm r$", re.IGNORECASE) + + +def fold(s: str) -> str: + return s.strip().casefold() + + +@dataclass +class SheetRow: + row_number: int # 1-based line number in the file + barcode: str + position: str + values: dict[str, str] # original header -> value (empty-header columns dropped) + sample_name: str | None = None + reported_sex: str = "U" + + @property + def sample_id(self) -> str: + return f"{self.barcode}_{self.position}" + + +@dataclass +class SampleSheet: + path: str + header: dict[str, str] + root_folder: Path | None + columns: list[str] + rows: list[SheetRow] + barcode_column: str + position_column: str + warnings: list[tuple[str, str]] = field(default_factory=list) + + @property + def has_sample_name(self) -> bool: + return any(fold(c) == fold("Sample_Name") for c in self.columns) + + def passthrough_columns(self) -> list[str]: + """Every non-key column in sheet order (IN-013).""" + keys = {fold(self.barcode_column), fold(self.position_column)} + return [c for c in self.columns if fold(c) not in keys] + + def by_sample_id(self) -> dict[str, SheetRow]: + return {r.sample_id: r for r in self.rows} + + +def _err(code: str, path: str, msg: str, hint: str | None = None) -> SimurgError: + return SimurgError(code, f"{path}: {msg}", hint=hint, path=path) + + +def normalise_sex(value: str) -> str: + v = fold(value) + if v in ("m", "male", "1"): + return "M" + if v in ("f", "female", "2"): + return "F" + return "U" + + +def validate_barcode(value: str, where: str, warnings: list[tuple[str, str]], path: str) -> str: + v = value.strip() + if _BARCODE_RE.match(v): + return v + m = _BARCODE_DECIMAL_RE.match(v) + if m: + warnings.append(("W_SHEET_BARCODE_REPAIRED", f"{where}: barcode {v!r} repaired to {m.group(1)}")) + return m.group(1) + hint = None + if _BARCODE_SCI_RE.match(v): + hint = "format the SentrixBarcode_A column as a number with 0 decimals (or as text) before saving" + raise _err("E_SHEET_BARCODE_FORMAT", path, f"{where}: invalid barcode {v!r}", hint) + + +def validate_position(value: str, where: str, path: str) -> str: + v = value.strip() + if not _POSITION_RE.match(v): + raise _err("E_SHEET_BARCODE_FORMAT", path, f"{where}: invalid position {v!r} (expected RxxCxx)") + return v.upper() + + +def _read_rows(path: Path) -> list[tuple[int, list[str]]]: + raw = path.read_bytes() + text = raw.decode("utf-8-sig", errors="replace") + stripped = text.lstrip() + if stripped.startswith("{") or (stripped.startswith("[") and stripped[1:].lstrip()[:1] in ("{", '"', "]")): + raise _err( + "E_SHEET_DELIMITER", + str(path), + "looks like JSON, not a CSV sample sheet", + "export the sheet as comma-separated CSV", + ) + rows: list[tuple[int, list[str]]] = [] + reader = csv.reader(io.StringIO(text, newline="")) + for row in reader: + rows.append((reader.line_num, row)) + return rows + + +def parse_samplesheet(path: str | Path, *, on_duplicate_sample: str = "error") -> SampleSheet: + """Parse and validate a sample sheet. Raises SimurgError with an E_SHEET_* code.""" + p = Path(path) + spath = str(p) + warnings: list[tuple[str, str]] = [] + rows = _read_rows(p) + + header_kv: dict[str, str] = {} + data_rows: list[tuple[int, list[str]]] = [] + section: str | None = None + saw_section = False + for line_no, row in rows: + if not row or all(not c.strip() for c in row): + continue + first = row[0].strip() + if first.startswith("[") and first.endswith("]") and first[1:-1].strip().replace(" ", "").isalnum(): + section = fold(first[1:-1]) + saw_section = True + continue + if not saw_section or section == "data": + data_rows.append((line_no, row)) + elif section == "header": + header_kv[first] = row[1].strip() if len(row) > 1 else "" + # [Manifests] and unknown sections are ignored + + if not data_rows: + raise _err("E_SHEET_EMPTY", spath, "no [Data] header row") + head_line, head = data_rows[0] + if len(head) == 1 and (";" in head[0] or "\t" in head[0]): + raise _err( + "E_SHEET_DELIMITER", + spath, + f"line {head_line}: header is not comma-delimited", + "save the sheet as comma-separated CSV (not semicolon or tab)", + ) + + # columns: drop empty headers (IN-004), check length (IN-006) and case-insensitive duplicates (IN-003) + col_idx: list[tuple[int, str]] = [] + empty_cols = 0 + seen: dict[str, str] = {} + for i, name in enumerate(head): + nm = name.strip() + if not nm: + empty_cols += 1 + continue + if len(nm) > MAX_LEN: + raise _err( + "E_SHEET_TOO_LONG", spath, f"column name at position {i + 1} has {len(nm)} characters (limit {MAX_LEN})" + ) + k = fold(nm) + if k in seen: + raise _err( + "E_SHEET_DUP_COLUMN", + spath, + f"Duplicate column '{nm}' (column names are case-insensitive). " + f"Rename or remove one of: {seen[k]}, {nm}", + ) + seen[k] = nm + col_idx.append((i, nm)) + if empty_cols: + warnings.append(("W_SHEET_EMPTY_COLUMN", f"{empty_cols} column(s) with an empty header ignored")) + columns = [nm for _, nm in col_idx] + folded = {fold(c): c for c in columns} + + key: tuple[str, str] | None = None + for bc, pc in KEY_PAIRS: + if fold(bc) in folded and fold(pc) in folded: + key = (folded[fold(bc)], folded[fold(pc)]) + break + if key is None: + raise _err( + "E_SHEET_NO_KEY_COLUMNS", + spath, + "missing SentrixBarcode_A/SentrixPosition_A (or beadChipName/sampleSectionName) columns", + ) + + recognised = {fold(c) for c in RECOGNISED} | {fold(key[0]), fold(key[1])} + reserved = {fold(c) for c in RESERVED_OUTPUT_COLUMNS} + for c in columns: + k = fold(c) + if k in recognised: + continue + if k in reserved or _NORM_RE.match(c) or _NORM_COMPAT_RE.match(c): + raise _err( + "E_SHEET_RESERVED", + spath, + f"column {c!r} is reserved for SimurgArray output", + "rename the column (or use --namespace-sheet-columns)", + ) + + out_rows: list[SheetRow] = [] + ids: dict[str, int] = {} + sex_cols = [folded[fold(c)] for c in SEX_COLUMNS_PRIORITY if fold(c) in folded] + name_col = folded.get(fold("Sample_Name")) + for line_no, row in data_rows[1:]: + values: dict[str, str] = {} + for i, nm in col_idx: + v = row[i] if i < len(row) else "" + if len(v) > MAX_LEN: + raise _err( + "E_SHEET_TOO_LONG", + spath, + f"line {line_no}, column {nm!r}: value has {len(v)} characters (limit {MAX_LEN})", + ) + values[nm] = v + where = f"line {line_no}" + barcode = validate_barcode(values[key[0]], where, warnings, spath) + position = validate_position(values[key[1]], where, spath) + sr = SheetRow(line_no, barcode, position, values) + if name_col is not None: + sr.sample_name = values[name_col].strip() + sexes = [normalise_sex(values[c]) for c in sex_cols] + known = [s for s in sexes if s != "U"] + if len(set(known)) > 1: + warnings.append(("W_SHEET_SEX_CONFLICT", f"{where}: reported-sex columns disagree; using {sex_cols[0]}")) + sr.reported_sex = sexes[0] if sexes else "U" + if sr.sample_id in ids: + if on_duplicate_sample == "drop": + warnings.append( + ( + "E_SHEET_DUP_SAMPLE", + f"{where}: duplicate sample {sr.sample_id} dropped (first at line {ids[sr.sample_id]})", + ) + ) + continue + raise _err( + "E_SHEET_DUP_SAMPLE", + spath, + f"{where}: sample {sr.sample_id} already listed at line {ids[sr.sample_id]}", + ) + ids[sr.sample_id] = line_no + out_rows.append(sr) + + if not out_rows: + raise _err("E_SHEET_EMPTY", spath, "no data rows") + + if name_col is not None: + names: dict[str, int] = {} + for sr in out_rows: + if not sr.sample_name: + raise _err( + "E_SHEET_NAME_EMPTY", + spath, + f"line {sr.row_number}: Sample_Name is empty (it must be set in every row when the column exists)", + ) + if sr.sample_name in ids and sr.sample_name != sr.sample_id: + raise _err( + "E_SHEET_NAME_COLLISION", + spath, + f"line {sr.row_number}: Sample_Name {sr.sample_name!r} equals the ID of another sample", + ) + if sr.sample_name in names: + raise _err( + "E_SHEET_NAME_COLLISION", + spath, + f"line {sr.row_number}: Sample_Name " + f"{sr.sample_name!r} already used at line {names[sr.sample_name]}", + ) + names[sr.sample_name] = sr.row_number + + root: Path | None = None + for k, v in header_kv.items(): + if fold(k) == "rootfolder" and v: + rp = Path(v) + root = rp if rp.is_absolute() else (p.parent / rp) + return SampleSheet( + path=spath, + header=header_kv, + root_folder=root, + columns=columns, + rows=out_rows, + barcode_column=key[0], + position_column=key[1], + warnings=warnings, + ) + + +_STEM_BAD = re.compile(r"[^A-Za-z0-9._-]") + + +def output_stem(name: str) -> str: + """File-name stem for a sample (IN-012).""" + s = _STEM_BAD.sub("_", name).lstrip(".") + return (s or "_")[:200] diff --git a/src/simurg/genotype/__init__.py b/src/simurg/genotype/__init__.py new file mode 100644 index 0000000..2eccfab --- /dev/null +++ b/src/simurg/genotype/__init__.py @@ -0,0 +1 @@ +"""SimurgArray genotype package.""" diff --git a/src/simurg/harmonize/__init__.py b/src/simurg/harmonize/__init__.py new file mode 100644 index 0000000..bb07a5b --- /dev/null +++ b/src/simurg/harmonize/__init__.py @@ -0,0 +1 @@ +"""SimurgArray harmonize package.""" diff --git a/src/simurg/pgx/__init__.py b/src/simurg/pgx/__init__.py new file mode 100644 index 0000000..066612c --- /dev/null +++ b/src/simurg/pgx/__init__.py @@ -0,0 +1 @@ +"""SimurgArray pgx package.""" diff --git a/src/simurg/pipeline/__init__.py b/src/simurg/pipeline/__init__.py new file mode 100644 index 0000000..2f4c9e2 --- /dev/null +++ b/src/simurg/pipeline/__init__.py @@ -0,0 +1 @@ +"""SimurgArray pipeline package.""" diff --git a/src/simurg/qc/__init__.py b/src/simurg/qc/__init__.py new file mode 100644 index 0000000..fa9a273 --- /dev/null +++ b/src/simurg/qc/__init__.py @@ -0,0 +1 @@ +"""SimurgArray qc package.""" diff --git a/src/simurg/signal/__init__.py b/src/simurg/signal/__init__.py new file mode 100644 index 0000000..611e4d1 --- /dev/null +++ b/src/simurg/signal/__init__.py @@ -0,0 +1 @@ +"""SimurgArray signal package.""" diff --git a/src/simurg/synth/__init__.py b/src/simurg/synth/__init__.py new file mode 100644 index 0000000..d41152a --- /dev/null +++ b/src/simurg/synth/__init__.py @@ -0,0 +1 @@ +"""SimurgArray synth package.""" diff --git a/src/simurg/synth/files.py b/src/simurg/synth/files.py new file mode 100644 index 0000000..fcf575b --- /dev/null +++ b/src/simurg/synth/files.py @@ -0,0 +1,247 @@ +"""Writers for small synthetic IDAT, BPM and EGT files (plan §16.4). + +They produce structurally valid files for tests and synthetic data sets. Layouts mirror the readers in +:mod:`simurg.formats`; they are not meant to reproduce vendor files byte-for-byte. +""" + +from __future__ import annotations + +from collections.abc import Sequence +from dataclasses import dataclass, field + +import numpy as np +import numpy.typing as npt + +from simurg.formats import idat as idatfmt +from simurg.formats.binutil import BinaryWriter + + +@dataclass +class SynthIdat: + address_ids: npt.NDArray[np.int32] + mean: npt.NDArray[np.uint16] + sd: npt.NDArray[np.uint16] | None = None + n_beads: npt.NDArray[np.uint8] | None = None + barcode: str = "200000000001" + position: str = "R01C01" + chip_type: str = "BeadChip 24x1x4" + sample_id: str = "" + plate: str = "" + well: str = "" + manifest: str = "" + run_info: list[tuple[str, str, str, str, str]] = field( + default_factory=lambda: [ + ("12/27/2023 5:52:21 PM", "Scan", "ScannerID=N0", "iScan Control Software", "4.2.1.729") + ] + ) + mid_block_itemsize: int = 4 + version: int = 3 + + +def encode_idat(s: SynthIdat) -> bytes: + n = len(s.address_ids) + sd = s.sd if s.sd is not None else np.full(n, 100, np.uint16) + nb = s.n_beads if s.n_beads is not None else np.full(n, 12, np.uint8) + sections: list[tuple[int, bytes]] = [] + + def sec(code: int, *parts: tuple[str, object]) -> None: + w = BinaryWriter() + for kind, val in parts: + if kind == "i32": + w.i32(int(val)) # type: ignore[call-overload] + elif kind == "str": + w.string(str(val)) + elif kind == "raw": + w.write(val) # type: ignore[arg-type] + else: + w.array(val, kind) # type: ignore[arg-type] + sections.append((code, w.getvalue())) + + mid_dtype = " bytes: + w = BinaryWriter() + w.write(b"BPM") + w.u8(1) + w.i32(bpm_version) + w.string(manifest_name) + w.string(control_config) + n = len(loci) + w.i32(n) + w.array(np.arange(1, n + 1), "= 6: + w.u8(0) + w.write(int(lc.exp_clusters).to_bytes(1, "little", signed=True)) + w.u8(lc.intensity_only) + w.u8(lc.assay_type) + if entry_version >= 7: + for f in (0.25, 0.25, 0.25, 0.25): + w.f32(f) + if entry_version >= 8: + w.string(lc.ref_strand) + w.i32(len(trailer)) + for line in trailer: + w.string(line) + return w.getvalue() + + +def encode_egt( + names: Sequence[str], + *, + theta_mean: npt.ArrayLike, + r_mean: npt.ArrayLike, + theta_dev: npt.ArrayLike | None = None, + r_dev: npt.ArrayLike | None = None, + n: npt.ArrayLike | None = None, + gentrain: npt.ArrayLike | None = None, + intensity_threshold: npt.ArrayLike | None = None, + data_block_version: int = 9, + manifest_name: str = "SYNTH-1v1_A1", + corrupt_counts: bool = False, +) -> bytes: + k = len(names) + tm = np.asarray(theta_mean, np.float32).reshape(k, 3) + rm = np.asarray(r_mean, np.float32).reshape(k, 3) + td = np.asarray(theta_dev if theta_dev is not None else np.full((k, 3), 0.03), np.float32).reshape(k, 3) + rd = np.asarray(r_dev if r_dev is not None else np.full((k, 3), 0.1), np.float32).reshape(k, 3) + nn = np.asarray(n if n is not None else np.full((k, 3), 100), np.int32).reshape(k, 3) + gt = np.asarray(gentrain if gentrain is not None else np.full(k, 0.8), np.float32) + thr = np.asarray(intensity_threshold if intensity_threshold is not None else np.zeros(k), np.float32) + w = BinaryWriter() + w.i32(3) + for s in ("< 6.1.3.0", "7.0.0", "7.0.0", "1.2.0", "1/1/2026 12:00:00 PM"): + w.string(s) + w.u8(1) + w.string(manifest_name) + w.i32(data_block_version) + w.string("") + w.i32(k) + for i in range(k): + w.array(nn[i], "= 7: + w.f32(float(thr[i])) + w.array(np.full(14, np.nan), "= 6: + w.array(np.zeros(15), " None: + missing = [p for p in paths if not p.exists()] + if missing: + pytest.skip(f"missing test data: {missing[0]} (run scripts/fetch_data.sh p0)") + + +@pytest.fixture(scope="session") +def beadarrayfiles(): # type: ignore[no-untyped-def] + """Illumina BeadArrayFiles (BSD-2) as an oracle.""" + path = REFS / "BeadArrayFiles" + if not (path / "module").is_dir(): + pytest.skip("refs/BeadArrayFiles not present (run scripts/bootstrap_refs.sh)") + sys.path.insert(0, str(path)) + import module # type: ignore[import-not-found] + + return module diff --git a/tests/oracle/test_p1_real_files.py b/tests/oracle/test_p1_real_files.py new file mode 100644 index 0000000..662b06d --- /dev/null +++ b/tests/oracle/test_p1_real_files.py @@ -0,0 +1,182 @@ +"""Phase P1 acceptance on real files (plan §18 P1): parsers agree field-by-field with Illumina's BeadArrayFiles +(BSD-2) and GTCs round-trip byte-identically. Needs ``scripts/fetch_data.sh p0``; GTC oracles are produced with +gtc2vcf ``idat2gtc`` (see docs/validation/P1.md). +""" + +from __future__ import annotations + +import numpy as np +import pytest +from tests.conftest import ( + BIOC_IDATS, + GDA_EGT, + GDA_SHEET, + GSA_BPM, + GSA_CSV, + GSA_EGT, + GSA_IDATS, + ORACLE_GTCS, + need, +) + +from simurg.core.errors import FormatError +from simurg.formats.bpm import read_bpm +from simurg.formats.csv_manifest import compare_manifests, read_csv_manifest +from simurg.formats.discovery import discover_idats +from simurg.formats.egt import read_egt +from simurg.formats.gtc import encode_gtc, read_gtc +from simurg.formats.idat import check_pair, read_idat +from simurg.formats.samplesheet import parse_samplesheet + +pytestmark = pytest.mark.data + + +@pytest.fixture(scope="module") +def gsa_bpm(): # type: ignore[no-untyped-def] + need(GSA_BPM) + return read_bpm(GSA_BPM) + + +def test_gsa_idats_all_parse_and_pair() -> None: + need(GSA_IDATS) + d = discover_idats(GSA_IDATS) + assert len(d.pairs) == 36 and all(p.complete for p in d.pairs.values()) + ref = None + for pair in d.pairs.values(): + g, r = read_idat(pair.grn), read_idat(pair.red) # type: ignore[arg-type] + check_pair(g, r) + assert g.barcode == pair.barcode and g.position == pair.position + assert g.scan_time_iso # every demo IDAT has a scan entry + ids = np.sort(g.address_ids) + if ref is None: + ref = ids + assert np.array_equal(ids, ref) + + +def test_bioc_idats() -> None: + need(BIOC_IDATS) + ok = { + "4019585376_B_Red.idat": 381079, + "4019585376_B_Red.idat.gz": 381079, + "5723646052_R02C02_Grn.idat": 622399, + "5723646052_R02C02_Grn.idat.gz": 622399, + } + for name, n in ok.items(): + assert read_idat(BIOC_IDATS / name).n_addresses == n + for name in ("4343238080_A_Grn.idat", "6016741005_A_Grn.idat"): + with pytest.raises(FormatError) as e: + read_idat(BIOC_IDATS / name) + assert e.value.code == "E_IDAT_ENCRYPTED" + + +def test_idat_addresses_cover_manifest(gsa_bpm) -> None: # type: ignore[no-untyped-def] # precursor of IN-027 + need(GSA_IDATS) + pair = next(iter(discover_idats(GSA_IDATS).pairs.values())) + ids = read_idat(pair.grn).address_ids # type: ignore[arg-type] + wanted = np.concatenate([gsa_bpm.address_a, gsa_bpm.address_b[gsa_bpm.address_b != 0]]) + assert np.isin(wanted, ids).mean() == 1.0 + + +def test_gtc_matches_beadarrayfiles_and_roundtrips(beadarrayfiles) -> None: # type: ignore[no-untyped-def] + need(ORACLE_GTCS) + files = sorted(ORACLE_GTCS.glob("*.gtc")) + assert files + GenotypeCalls = beadarrayfiles.GenotypeCalls + for f in files: + g = read_gtc(f) + assert encode_gtc(g) == f.read_bytes(), f"round-trip differs: {f}" + b = GenotypeCalls(str(f)) + assert np.array_equal(g.raw_x, b.get_raw_x_intensities()) + assert np.array_equal(g.raw_y, b.get_raw_y_intensities()) + assert g.genotypes.tolist() == list(b.get_genotypes()) + assert np.array_equal(g.genotype_scores, b.get_genotype_scores(), equal_nan=True) + assert np.array_equal(g.b_allele_freqs, b.get_ballele_freqs(), equal_nan=True) + assert np.array_equal(g.logr_ratios, b.get_logr_ratios(), equal_nan=True) + assert g.base_call_strings() == [x.decode() for x in b.get_base_calls()] + assert (g.call_rate, g.logr_dev, g.gc10, g.gc50) == ( + b.get_call_rate(), + b.get_logr_dev(), + b.get_gc10(), + b.get_gc50(), + ) + assert g.gender == b.get_gender() + assert (g.num_calls, g.num_no_calls, g.num_intensity_only) == ( + b.get_num_calls(), + b.get_num_no_calls(), + b.get_num_intensity_only(), + ) + assert g.percentiles_x.tolist() == list(b.get_percentiles_x()) + xf = [ + (t.offset_x, t.offset_y, t.scale_x, t.scale_y, t.shear, t.theta) for t in b.get_normalization_transforms() + ] + mine = [ + tuple(float(x[k]) for k in ("offset_x", "offset_y", "scale_x", "scale_y", "shear", "theta")) + for x in g.normalization_transforms + ] + assert mine == [tuple(float(v) for v in t) for t in xf] + + +def test_gtc_norm_transform_count_matches_bpm_bins(gsa_bpm) -> None: # type: ignore[no-untyped-def] + need(ORACLE_GTCS) + g = read_gtc(sorted(ORACLE_GTCS.glob("*.gtc"))[0]) + assert len(g.normalization_transforms) == len(np.unique(gsa_bpm.norm_bin_ids())) + assert g.num_snps == gsa_bpm.num_loci + + +@pytest.mark.slow +def test_bpm_matches_beadarrayfiles(gsa_bpm, beadarrayfiles) -> None: # type: ignore[no-untyped-def] + b = beadarrayfiles.BeadPoolManifest(str(GSA_BPM)) + m = gsa_bpm + assert m.names == b.names + assert m.snps == b.snps + assert m.chroms == b.chroms + assert m.map_infos.tolist() == b.map_infos + assert m.address_a.tolist() == b.addresses + assert m.assay_types.tolist() == b.assay_types + assert m.normalization_lookups().tolist() == b.normalization_lookups + rs = beadarrayfiles.RefStrand + assert [rs.from_string(s) for s in m.ref_strands] == b.ref_strands + assert m.manifest_name == b.manifest_name and m.control_config == b.control_config + + +@pytest.mark.slow +@pytest.mark.parametrize("egt_path", [GSA_EGT, GDA_EGT], ids=["gsa", "gda"]) +def test_egt_matches_beadarrayfiles(egt_path, beadarrayfiles) -> None: # type: ignore[no-untyped-def] + need(egt_path) + e = read_egt(egt_path) + with open(egt_path, "rb") as fh: + b = beadarrayfiles.ClusterFile.read_cluster_file(fh) + assert e.num_records == len(b.records) if hasattr(b, "records") else True + idx = e.name_index() + for name in e.names[:: max(1, e.num_records // 5000)]: # ~5000 evenly spaced records + rec = b.get_record(name) + i = idx[name] + for k, stats in enumerate((rec.aa_cluster_stats, rec.ab_cluster_stats, rec.bb_cluster_stats)): + assert e.n[i, k] == stats.N + assert e.theta_mean[i, k] == np.float32(stats.theta_mean) + assert e.theta_dev[i, k] == np.float32(stats.theta_dev) + assert e.r_mean[i, k] == np.float32(stats.r_mean) + assert e.r_dev[i, k] == np.float32(stats.r_dev) + assert e.intensity_threshold[i] == np.float32(rec.intensity_threshold) + assert e.gentrain_score[i] == np.float32(rec.cluster_score.total_score) + assert e.cluster_sep[i] == np.float32(rec.cluster_score.cluster_separation) + assert e.addresses[i] == rec.address + + +@pytest.mark.slow +def test_csv_manifest_consistent_with_bpm(gsa_bpm) -> None: # type: ignore[no-untyped-def] # IN-042 + need(GSA_CSV) + c = read_csv_manifest(GSA_CSV) + assert compare_manifests(gsa_bpm, c) == [] + # Assay type inferred from CSV sequences (fallback for BPM locus entries < v6, gtc2vcf rule) agrees with the + # BPM's explicit value for 99.71% of Infinium I loci on GSA v4; the BPM is authoritative (known_unknowns K-003). + inf1 = gsa_bpm.address_b != 0 + agree = (c.assay_types[inf1] == gsa_bpm.assay_types[inf1]).mean() + assert agree >= 0.995 + + +def test_gda_demo_sheet() -> None: + need(GDA_SHEET) + s = parse_samplesheet(GDA_SHEET) + assert len(s.rows) == 36 and s.has_sample_name + assert {r.values["Parent1"] for r in s.rows} >= {"NA12891"} diff --git a/tests/unit/test_binutil.py b/tests/unit/test_binutil.py new file mode 100644 index 0000000..47a6229 --- /dev/null +++ b/tests/unit/test_binutil.py @@ -0,0 +1,68 @@ +from __future__ import annotations + +import gzip + +import numpy as np +import pytest +from hypothesis import given +from hypothesis import strategies as st + +from simurg.core.errors import FormatError +from simurg.formats.binutil import BinaryReader, BinaryWriter, load_bytes + + +@given(st.text(max_size=400)) +def test_string_roundtrip(s: str) -> None: + w = BinaryWriter() + w.string(s) + r = BinaryReader(w.getvalue(), "E_GTC_CORRUPT") + assert r.string() == s + assert r.remaining() == 0 + + +@pytest.mark.parametrize("n", [0, 1, 127, 128, 16383, 16384, 300_000]) +def test_varint_boundaries(n: int) -> None: + w = BinaryWriter() + w.string(b"x" * n) + r = BinaryReader(w.getvalue(), "E_GTC_CORRUPT") + assert r.varint_length() == n + + +def test_invalid_utf8_is_replaced_not_raised() -> None: # IN-033, EDGE-002 + w = BinaryWriter() + w.string(b"ok\xff\xfe") + assert BinaryReader(w.getvalue(), "E_GTC_CORRUPT").string() == "ok��" + + +def test_string_past_eof_is_format_error() -> None: + buf = bytes([50]) + b"short" + with pytest.raises(FormatError) as e: + BinaryReader(buf, "E_GTC_CORRUPT", "x.gtc").string() + assert e.value.code == "E_GTC_CORRUPT" + + +def test_bad_varint_is_format_error() -> None: + with pytest.raises(FormatError): + BinaryReader(b"\xff\xff\xff\xff\xff\xff", "E_BPM_FORMAT").varint_length() + + +def test_array_bounds_and_seek() -> None: + r = BinaryReader(np.arange(4, dtype=" None: # type: ignore[no-untyped-def] + p = tmp_path / "a.bin.gz" + p.write_bytes(gzip.compress(b"hello")) + assert bytes(load_bytes(p)) == b"hello" + q = tmp_path / "b.bin" + q.write_bytes(b"world") + assert bytes(load_bytes(q)) == b"world" + e = tmp_path / "empty.bin" + e.write_bytes(b"") + assert bytes(load_bytes(e)) == b"" diff --git a/tests/unit/test_bpm_egt.py b/tests/unit/test_bpm_egt.py new file mode 100644 index 0000000..f0b39b7 --- /dev/null +++ b/tests/unit/test_bpm_egt.py @@ -0,0 +1,139 @@ +"""BPM and EGT readers on synthetic files (spec IN-040, IN-044; plan §6.2 normalization-bin rule).""" + +from __future__ import annotations + +import numpy as np +import pytest + +from simurg.core.errors import FormatError +from simurg.formats.bpm import infer_assay_type, parse_control_config, read_bpm +from simurg.formats.egt import read_egt +from simurg.synth.files import SynthLocus, encode_bpm, encode_egt + + +def _loci() -> list[SynthLocus]: + return [ + SynthLocus("rs1", "1", 100, "[A/G]", 11, 0, 0, norm_id=1), + SynthLocus("rs2", "2", 200, "[A/T]", 12, 13, 1, norm_id=1, ilmn_strand="BOT", ref_strand="-"), + SynthLocus("rs3", "X", 300, "[C/G]", 14, 15, 2, norm_id=99), + SynthLocus("rs4", "Y", 400, "[T/C]", 16, 0, 0, norm_id=2, intensity_only=1), + SynthLocus("rs5", "MT", 500, "[A/C]", 17, 0, 0, norm_id=99), + ] + + +def test_bpm_roundtrip(tmp_path) -> None: # type: ignore[no-untyped-def] + p = tmp_path / "m.bpm" + p.write_bytes(encode_bpm(_loci(), control_config="0000000001,Staining,Red,DNP (High)\n")) + m = read_bpm(p) + assert m.names == ["rs1", "rs2", "rs3", "rs4", "rs5"] + assert m.chroms == ["1", "2", "X", "Y", "MT"] + assert m.map_infos.tolist() == [100, 200, 300, 400, 500] + assert m.address_a.tolist() == [11, 12, 14, 16, 17] + assert m.address_b.tolist() == [0, 13, 15, 0, 0] + assert m.assay_types.tolist() == [0, 1, 2, 0, 0] + assert m.intensity_only.tolist() == [False, False, False, True, False] + assert m.ref_strands == ["+", "-", "+", "+", "+"] + assert m.ilmn_strands[1] == "BOT" + assert m.controls == [{"address": "0000000001", "category": "Staining", "color": "Red", "name": "DNP (High)"}] + assert m.header_lines == ["Illumina, Inc."] + + +def test_norm_bins_mod256_and_rank(tmp_path) -> None: # type: ignore[no-untyped-def] + # rs3: 99 + 100*2 = 299 -> 43 (byte wrap) ; rs2: 1 + 100 = 101 + p = tmp_path / "m.bpm" + p.write_bytes(encode_bpm(_loci())) + m = read_bpm(p) + assert m.norm_bin_ids().tolist() == [1, 101, 43, 2, 99] + # sorted distinct: [1, 2, 43, 99, 101] -> ranks + assert m.normalization_lookups().tolist() == [0, 4, 2, 1, 3] + + +def test_entries_mapped_by_name(tmp_path) -> None: # type: ignore[no-untyped-def] + p = tmp_path / "m.bpm" + p.write_bytes(encode_bpm(_loci(), entry_order=[4, 3, 2, 1, 0])) + m = read_bpm(p) + assert m.map_infos.tolist() == [100, 200, 300, 400, 500] + + +@pytest.mark.parametrize("entry_version", [4, 6, 7, 8]) +def test_entry_versions(tmp_path, entry_version) -> None: # type: ignore[no-untyped-def] + loci = [ + SynthLocus("rs1", source_seq="ACGT[A/G]TTGA", probe_a="ACGTA"), + SynthLocus("rs2", address_b=5, assay_type=1, source_seq="AAAT[A/T]AAAC", probe_a="ACGTA", probe_b="ACGTT"), + ] + p = tmp_path / "m.bpm" + p.write_bytes(encode_bpm(loci, entry_version=entry_version)) + m = read_bpm(p) + assert m.assay_types.tolist() == [0, 1] + assert (m.base_fractions is not None) == (entry_version >= 7) + assert m.has_ref_strand == (entry_version >= 8) + + +def test_bpm_errors(tmp_path) -> None: # type: ignore[no-untyped-def] + good = encode_bpm(_loci()) + cases = { + "magic": b"XPM" + good[3:], + "truncated": good[: len(good) // 2], + "bad_assay": encode_bpm([SynthLocus("rs1", address_b=0, assay_type=1)]), + "bad_norm": encode_bpm([SynthLocus("rs1", norm_id=150)]), + } + for label, data in cases.items(): + p = tmp_path / f"{label}.bpm" + p.write_bytes(data) + with pytest.raises(FormatError) as e: + read_bpm(p) + assert e.value.code == "E_BPM_FORMAT", label + + +def test_long_manifest_name(tmp_path) -> None: # type: ignore[no-untyped-def] # EDGE-006 + name = "X" * 200 + "_ünicode" + p = tmp_path / "long name (1).bpm" + p.write_bytes(encode_bpm(_loci(), manifest_name=name)) + assert read_bpm(p).manifest_name == name + + +def test_infer_assay_type() -> None: + assert infer_assay_type("ACGT", "", "AC[A/G]TT") == 0 + assert infer_assay_type("ACGA", "ACGT", "AAA[A/T]TTT") == 1 + assert infer_assay_type("ACGA", "ACGT", "CCC[C/G]GGG") == 2 + assert infer_assay_type("", "", "") == 255 + + +def test_control_config_parse() -> None: + txt = "0027630314:0027630314,Staining,Red,DNP (High)\n\nbad line\n1,Hyb,Green,Hyb (High), extra\n" + assert parse_control_config(txt)[1]["name"] == "Hyb (High),extra" + assert len(parse_control_config(txt)) == 2 + + +@pytest.mark.parametrize("dbv", [9, 8, 7, 6, 5, 4]) +def test_egt_versions(tmp_path, dbv) -> None: # type: ignore[no-untyped-def] + names = ["rs1", "rs2"] + tm = [[0.05, 0.5, 0.95], [0.1, 0.45, 0.9]] + rm = [[1.0, 1.1, 0.9], [0.8, 0.9, 1.0]] + p = tmp_path / "c.egt" + p.write_bytes(encode_egt(names, theta_mean=tm, r_mean=rm, intensity_threshold=[0.2, 0.0], data_block_version=dbv)) + e = read_egt(p) + assert e.names == names + np.testing.assert_allclose(e.theta_mean, tm, rtol=1e-6) + np.testing.assert_allclose(e.r_mean, rm, rtol=1e-6) + assert e.n.tolist() == [[100, 100, 100]] * 2 + if dbv >= 7: + assert e.intensity_threshold.tolist()[0] == pytest.approx(0.2) + else: + assert np.isnan(e.intensity_threshold).all() + np.testing.assert_allclose(e.mean_r(), np.mean(rm, axis=1), rtol=1e-6) + + +def test_egt_count_mismatch(tmp_path) -> None: # type: ignore[no-untyped-def] + p = tmp_path / "c.egt" + p.write_bytes(encode_egt(["rs1"], theta_mean=[[0.1, 0.5, 0.9]], r_mean=[[1, 1, 1]], corrupt_counts=True)) + with pytest.raises(FormatError) as e: + read_egt(p) + assert e.value.code == "E_EGT_LOAD" + + +def test_egt_bad_version(tmp_path) -> None: # type: ignore[no-untyped-def] + p = tmp_path / "c.egt" + p.write_bytes(b"\x04\x00\x00\x00" + b"\x00" * 50) + with pytest.raises(FormatError): + read_egt(p) diff --git a/tests/unit/test_cli_errors.py b/tests/unit/test_cli_errors.py new file mode 100644 index 0000000..ba53411 --- /dev/null +++ b/tests/unit/test_cli_errors.py @@ -0,0 +1,95 @@ +"""CLI behaviour (CLI-007, CLI-014, CLI-020) and the message catalogue (CLI-016).""" + +from __future__ import annotations + +import json +import subprocess +import sys +from pathlib import Path + +import numpy as np +from typer.testing import CliRunner + +from simurg.cli.main import app +from simurg.core.errors import CATALOGUE, render_markdown +from simurg.synth.files import SynthIdat, SynthLocus, encode_bpm, encode_egt, encode_idat + +ROOT = Path(__file__).resolve().parents[2] +runner = CliRunner() + + +def _files(tmp_path: Path) -> dict[str, Path]: + idat = tmp_path / "200000000001_R01C01_Grn.idat" + idat.write_bytes( + encode_idat(SynthIdat(address_ids=np.arange(1, 6, dtype=np.int32), mean=np.arange(5, dtype=np.uint16))) + ) + bpm = tmp_path / "m.bpm" + bpm.write_bytes(encode_bpm([SynthLocus("rs1"), SynthLocus("rs2", address_a=2)])) + egt = tmp_path / "c.egt" + egt.write_bytes(encode_egt(["rs1"], theta_mean=[[0.1, 0.5, 0.9]], r_mean=[[1, 1, 1]])) + return {"idat": idat, "bpm": bpm, "egt": egt} + + +def test_version() -> None: + r = runner.invoke(app, ["--version"]) + assert r.exit_code == 0 and r.stdout.startswith("SimurgArray ") + + +def test_inspect_json(tmp_path) -> None: # type: ignore[no-untyped-def] + f = _files(tmp_path) + r = runner.invoke(app, ["inspect", "--json", str(f["idat"]), str(f["bpm"]), str(f["egt"])]) + assert r.exit_code == 0, r.output + out = json.loads(r.stdout) + assert [o["kind"] for o in out] == ["idat", "bpm", "egt"] + assert out[1]["num_loci"] == 2 and out[2]["num_records"] == 1 + + +def test_inspect_missing_file_exit2() -> None: + r = runner.invoke(app, ["inspect", "/nonexistent/x.idat"]) + assert r.exit_code == 2 + assert "E_FILE_NOT_FOUND" in r.output + + +def test_inspect_partial_failure_exit3(tmp_path) -> None: # type: ignore[no-untyped-def] + f = _files(tmp_path) + bad = tmp_path / "bad.gtc" + bad.write_bytes(b"gtc\x05" + b"\x00" * 3) + r = runner.invoke(app, ["inspect", "--log-json", str(f["idat"]), str(bad)]) + assert r.exit_code == 3 + line = next(ln for ln in r.output.splitlines() if ln.startswith('{"schema"')) + rec = json.loads(line) + assert rec["schema"] == "simurg.log/1" and rec["code"] == "E_GTC_CORRUPT" and rec["level"] == "error" + + +def test_unknown_option_is_usage_error(tmp_path) -> None: # type: ignore[no-untyped-def] + r = runner.invoke(app, ["inspect", "--no-such-option", "x"]) + assert r.exit_code == 2 + + +def test_catalogue_integrity() -> None: + for code, d in CATALOGUE.items(): + assert code == d.code + assert code[:2] in ("E_", "W_", "I_") + assert d.level in ("error", "warning", "info") + + +def test_messages_doc_is_current() -> None: + doc = ROOT / "docs" / "messages.md" + assert doc.read_text(encoding="utf-8") == render_markdown(), "regenerate: python -m simurg.core.errors" + + +def test_locale_independent_output(tmp_path) -> None: # type: ignore[no-untyped-def] # EDGE-001 + f = _files(tmp_path) + outs = set() + for loc in ("C.UTF-8", "tr_TR.UTF-8", "de_DE.UTF-8"): + env = {"LC_ALL": loc, "PATH": "/usr/bin:/bin"} + r = subprocess.run( + [sys.executable, "-m", "simurg.cli.main", "inspect", "--json", str(f["bpm"])], + capture_output=True, + text=True, + env=env, + cwd=ROOT, + ) + assert r.returncode == 0, r.stderr + outs.add(r.stdout) + assert len(outs) == 1 diff --git a/tests/unit/test_discovery_atomic.py b/tests/unit/test_discovery_atomic.py new file mode 100644 index 0000000..d1d4748 --- /dev/null +++ b/tests/unit/test_discovery_atomic.py @@ -0,0 +1,120 @@ +"""IDAT discovery (IN-020 … IN-022, EDGE-005) and atomic writes / overwrite policy (CLI-010, CLI-011).""" + +from __future__ import annotations + +import os +from pathlib import Path + +import numpy as np +import pytest + +from simurg.core.atomic import atomic_open, atomic_write_bytes, check_outputs_absent +from simurg.core.errors import SimurgError +from simurg.formats import discovery as disc +from simurg.formats.samplesheet import parse_samplesheet +from simurg.synth.files import SynthIdat, encode_idat + + +def write_pair( + folder: Path, + barcode: str, + pos: str, + scan: str = "12/27/2023 5:52:21 PM", + channels: tuple[str, ...] = ("Grn", "Red"), +) -> None: + folder.mkdir(parents=True, exist_ok=True) + s = SynthIdat( + address_ids=np.arange(1, 11, dtype=np.int32), + mean=np.arange(10, dtype=np.uint16), + barcode=barcode, + position=pos, + run_info=[(scan, "Scan", "", "iScan", "1")], + ) + for ch in channels: + (folder / f"{barcode}_{pos}_{ch}.idat").write_bytes(encode_idat(s)) + + +def test_pairs_and_unpaired(tmp_path) -> None: # type: ignore[no-untyped-def] + write_pair(tmp_path / "a", "200000000001", "R01C01") + write_pair(tmp_path / "b" / "deep", "200000000001", "R02C01", channels=("Grn",)) + d = disc.discover_idats(tmp_path) + assert sorted(d.pairs) == ["200000000001_R01C01", "200000000001_R02C01"] + assert d.pairs["200000000001_R01C01"].complete + assert not d.pairs["200000000001_R02C01"].complete + assert [c for c, _ in d.warnings] == ["W_IDAT_UNPAIRED"] + + +def test_duplicate_pairs_error_lists_all_paths(tmp_path) -> None: # type: ignore[no-untyped-def] + write_pair(tmp_path / "run1", "200000000001", "R01C01") + write_pair(tmp_path / "run2", "200000000001", "R01C01") + with pytest.raises(SimurgError) as e: + disc.discover_idats(tmp_path) + assert e.value.code == "E_DUP_IDAT" + assert "run1" in e.value.message and "run2" in e.value.message + + +def test_duplicate_first_and_newest_deterministic(tmp_path, monkeypatch) -> None: # type: ignore[no-untyped-def] + write_pair(tmp_path / "z_old", "200000000001", "R01C01", scan="1/1/2020 1:00:00 PM") + write_pair(tmp_path / "a_new", "200000000001", "R01C01", scan="1/1/2024 1:00:00 PM") + first = disc.discover_idats(tmp_path, on_duplicate="first").pairs["200000000001_R01C01"] + newest = disc.discover_idats(tmp_path, on_duplicate="newest").pairs["200000000001_R01C01"] + assert first.folder.name == "a_new" # bytewise smallest path + assert newest.folder.name == "a_new" + # shuffled directory listing must not change anything + real_scandir = os.scandir + + def shuffled(path): # type: ignore[no-untyped-def] + return list(reversed(list(real_scandir(path)))) + + monkeypatch.setattr(os, "scandir", shuffled) + assert disc.discover_idats(tmp_path, on_duplicate="first").pairs["200000000001_R01C01"].folder.name == "a_new" + + +def test_symlink_loop(tmp_path) -> None: # type: ignore[no-untyped-def] + write_pair(tmp_path / "a", "200000000001", "R01C01") + (tmp_path / "a" / "loop").symlink_to(tmp_path, target_is_directory=True) + assert list(disc.discover_idats(tmp_path).pairs) == ["200000000001_R01C01"] + + +def test_resolve_with_sheet(tmp_path) -> None: # type: ignore[no-untyped-def] + write_pair(tmp_path / "idats", "200000000001", "R01C01") + write_pair(tmp_path / "idats", "200000000001", "R03C01") + sh = tmp_path / "s.csv" + sh.write_text( + "[Header]\nRootFolder,idats\n[Data]\nSentrixBarcode_A,SentrixPosition_A,Sample_Name\n" + "200000000001,R01C01,alpha\n200000000001,R02C01,beta\n" + ) + sheet = parse_samplesheet(sh) + folder, notes = disc.resolve_input_folder(None, sheet) + assert folder == tmp_path / "idats" and notes == [] + samples, notes = disc.resolve_samples(disc.discover_idats(folder), sheet) + assert [(s.sample_name, s.status) for s in samples] == [("alpha", "READY"), ("beta", "MISSING_INPUT")] + assert notes[0][0] == "I_UNLISTED_IDATS" + _, notes = disc.resolve_input_folder(tmp_path / "other", sheet) + assert notes[0][0] == "W_SHEET_ROOTFOLDER_OVERRIDDEN" + + +def test_input_rules() -> None: + with pytest.raises(SimurgError) as e: + disc.resolve_input_folder(None, None) + assert e.value.code == "E_NO_INPUT" + + +def test_atomic_write_and_failure(tmp_path) -> None: # type: ignore[no-untyped-def] + p = tmp_path / "out" / "x.bin" + atomic_write_bytes(p, b"abc") + assert p.read_bytes() == b"abc" + with pytest.raises(RuntimeError), atomic_open(tmp_path / "y.bin") as fh: + fh.write(b"partial") # type: ignore[arg-type] + raise RuntimeError("boom") + assert not (tmp_path / "y.bin").exists() + assert not (tmp_path / "y.bin.part").exists() + + +def test_overwrite_policy(tmp_path) -> None: # type: ignore[no-untyped-def] + p = tmp_path / "x" + p.write_text("1") + with pytest.raises(SimurgError) as e: + check_outputs_absent([p, tmp_path / "missing"], force=False) + assert e.value.code == "E_OUTPUT_EXISTS" + check_outputs_absent([p], force=True) diff --git a/tests/unit/test_gtc.py b/tests/unit/test_gtc.py new file mode 100644 index 0000000..7f244bf --- /dev/null +++ b/tests/unit/test_gtc.py @@ -0,0 +1,152 @@ +"""GTC codec: synthetic round-trips, property tests, corruption handling (spec IN-032, IN-033, EDGE-002).""" + +from __future__ import annotations + +import numpy as np +import pytest +from hypothesis import given, settings +from hypothesis import strategies as st + +from simurg.core.errors import FormatError +from simurg.formats.gtc import ( + CANONICAL_ORDER_V3, + CANONICAL_ORDER_V5, + XFORM_DTYPE, + Gtc, + ScannerData, + encode_gtc, + read_gtc, + write_gtc, +) + + +def make_gtc(n: int = 20, version: int = 5, seed: int = 0) -> Gtc: + rng = np.random.default_rng(seed) + xf = np.zeros(3, XFORM_DTYPE) + xf["version"] = 1 + xf["offset_x"] = rng.normal(size=3) + xf["theta"] = rng.normal(size=3) + g = Gtc( + version=version, + num_snps=n, + sample_name="NA12878 ünicode", + cluster_file="c.egt", + snp_manifest="m.bpm", + imaging_date="12/27/2023 5:53:34 PM", + autocall_date="01/01/26 12:00 PM", + autocall_version="simurg-0.1", + normalization_transforms=xf, + controls_x=rng.integers(0, 60000, 8).astype(np.uint16), + controls_y=rng.integers(0, 60000, 8).astype(np.uint16), + raw_x=rng.integers(0, 60000, n).astype(np.uint16), + raw_y=rng.integers(0, 60000, n).astype(np.uint16), + genotypes=rng.integers(0, 4, n).astype(np.uint8), + base_calls=np.frombuffer(b"AG" * n, np.uint8).reshape(n, 2).copy(), + genotype_scores=rng.random(n).astype(np.float32), + scanner_data=ScannerData("N103", b"\x01\x00\x00\x00", b"\x02\x00\x00\x00", "4.2.1", "user"), + call_rate=0.99, + gender=b"F", + logr_dev=0.12, + gc10=0.4, + b_allele_freqs=rng.random(n).astype(np.float32), + logr_ratios=rng.normal(size=n).astype(np.float32), + percentiles_x=np.array([1, 2, 3], np.uint16), + percentiles_y=np.array([4, 5, 6], np.uint16), + slide_identifier="200000000001", + ) + g.sample_data["p50gc"] = 0.7 + g.sample_data["num_calls"] = n - 2 + g.sample_data["num_no_calls"] = 2 + return g + + +@pytest.mark.parametrize("version", [3, 5]) +def test_roundtrip_canonical(tmp_path, version) -> None: # type: ignore[no-untyped-def] + g = make_gtc(version=version) + p = tmp_path / "s.gtc" + write_gtc(g, p) + h = read_gtc(p) + assert h.toc_order == (CANONICAL_ORDER_V5 if version == 5 else CANONICAL_ORDER_V3) + assert h.sample_name == g.sample_name + assert np.array_equal(h.raw_x, g.raw_x) and np.array_equal(h.genotypes, g.genotypes) + assert np.array_equal(h.genotype_scores, g.genotype_scores) + assert h.scanner_data == g.scanner_data + assert h.gc50 == pytest.approx(0.7) and h.num_calls == 18 + assert h.gender_char == "F" + assert h.base_call_strings()[0] == "AG" + assert encode_gtc(h) == p.read_bytes() # byte-identical round-trip + if version == 5: + assert np.array_equal(h.logr_ratios, g.logr_ratios) + + +def test_noncanonical_layout_preserved(tmp_path) -> None: # type: ignore[no-untyped-def] + g = make_gtc() + g.toc_order = list(reversed(CANONICAL_ORDER_V5)) + g.data_order = [t for t in CANONICAL_ORDER_V5 if t > 3][::-1] + data = encode_gtc(g) + p = tmp_path / "s.gtc" + p.write_bytes(data) + h = read_gtc(p) + assert h.toc_order == g.toc_order and h.data_order == g.data_order + assert encode_gtc(h) == data + + +@settings(max_examples=25, deadline=None) +@given(n=st.integers(0, 200), seed=st.integers(0, 10**6), name=st.text(max_size=300)) +def test_property_roundtrip(tmp_path_factory, n, seed, name) -> None: # type: ignore[no-untyped-def] + g = make_gtc(n=n, seed=seed) + g.sample_name = name + data = encode_gtc(g) + p = tmp_path_factory.mktemp("g") / "x.gtc" + p.write_bytes(data) + h = read_gtc(p) + assert h.sample_name == name and h.num_snps == n + assert encode_gtc(h) == data + + +def test_truncated_and_bad_magic(tmp_path) -> None: # type: ignore[no-untyped-def] + data = encode_gtc(make_gtc()) + for label, blob in { + "trunc": data[: len(data) // 3], + "magic": b"xtc" + data[3:], + "tiny": data[:5], + "version": data[:3] + b"\x09" + data[4:], + }.items(): + p = tmp_path / f"{label}.gtc" + p.write_bytes(blob) + with pytest.raises(FormatError) as e: + read_gtc(p) + assert e.value.code == "E_GTC_CORRUPT", label + + +def test_bad_string_length_is_clean_error(tmp_path) -> None: # type: ignore[no-untyped-def] # EDGE-002 + g = make_gtc() + data = bytearray(encode_gtc(g)) + header = 8 + 6 * len(CANONICAL_ORDER_V5) + data[header] = 0xFF # sample-name length prefix now runs past EOF + data[header + 1] = 0xFF + data[header + 2] = 0xFF + data[header + 3] = 0x7F + p = tmp_path / "s.gtc" + p.write_bytes(bytes(data)) + with pytest.raises(FormatError) as e: + read_gtc(p) + assert e.value.code == "E_GTC_CORRUPT" + + +def test_incomplete_write_detection(tmp_path) -> None: # type: ignore[no-untyped-def] + g = make_gtc() + g.toc_order = [t for t in CANONICAL_ORDER_V5 if t != 1016] + g.data_order = [t for t in g.toc_order if t > 3] + p = tmp_path / "s.gtc" + p.write_bytes(encode_gtc(g)) + with pytest.raises(FormatError): + read_gtc(p) + assert read_gtc(p, check_write_complete=False).num_snps == 20 + + +def test_length_mismatch_rejected_on_write() -> None: + g = make_gtc() + g.raw_x = g.raw_x[:-1] + with pytest.raises(ValueError): + encode_gtc(g) diff --git a/tests/unit/test_idat.py b/tests/unit/test_idat.py new file mode 100644 index 0000000..5ee180b --- /dev/null +++ b/tests/unit/test_idat.py @@ -0,0 +1,128 @@ +"""IDAT reader: synthetic files + malformed-input matrix (spec IN-023, EDGE-004, EDGE-033).""" + +from __future__ import annotations + +import gzip +import struct + +import numpy as np +import pytest + +from simurg.core.errors import FormatError +from simurg.formats.idat import check_pair, parse_idat_filename, parse_run_time, read_idat +from simurg.synth.files import SynthIdat, encode_idat + + +def _synth(n: int = 50, **kw) -> SynthIdat: # type: ignore[no-untyped-def] + rng = np.random.default_rng(1) + ids = np.sort(rng.choice(10_000_000, n, replace=False)).astype(np.int32) + return SynthIdat( + address_ids=ids, + mean=rng.integers(0, 60000, n).astype(np.uint16), + n_beads=rng.integers(0, 30, n).astype(np.uint8), + **kw, + ) + + +def test_roundtrip(tmp_path) -> None: # type: ignore[no-untyped-def] + s = _synth(sample_id="NA12878", plate="P1", well="A01") + p = tmp_path / "200000000001_R01C01_Grn.idat" + p.write_bytes(encode_idat(s)) + i = read_idat(p) + assert i.n_addresses == 50 + assert np.array_equal(i.address_ids, s.address_ids) + assert np.array_equal(i.mean, s.mean) + assert np.array_equal(i.n_beads, s.n_beads) + assert (i.barcode, i.position, i.sample_id, i.plate, i.well) == ("200000000001", "R01C01", "NA12878", "P1", "A01") + assert i.scan_time_iso == "2023-12-27T17:52:21" + assert i.scanner_software == "4.2.1.729" + assert i.mid_block_itemsize == 4 + + +def test_midblock_uint8_inferred(tmp_path) -> None: # type: ignore[no-untyped-def] + p = tmp_path / "x.idat" + p.write_bytes(encode_idat(_synth(n=40, mid_block_itemsize=1))) + assert read_idat(p).mid_block_itemsize == 1 + + +def test_gzip(tmp_path) -> None: # type: ignore[no-untyped-def] + p = tmp_path / "200000000001_R01C01_Red.idat.gz" + p.write_bytes(gzip.compress(encode_idat(_synth()))) + assert read_idat(p).n_addresses == 50 + + +def test_metadata_only(tmp_path) -> None: # type: ignore[no-untyped-def] + p = tmp_path / "x.idat" + p.write_bytes(encode_idat(_synth())) + i = read_idat(p, load_arrays=False) + assert i.n_addresses == 50 and len(i.mean) == 0 and i.barcode == "200000000001" + + +@pytest.mark.parametrize( + "mutate,code", + [ + (lambda b: b"XDAT" + b[4:], "E_IDAT_FORMAT"), # wrong magic + (lambda b: b[:4] + struct.pack(" None: # type: ignore[no-untyped-def] + p = tmp_path / "bad.idat" + p.write_bytes(mutate(encode_idat(_synth()))) + with pytest.raises(FormatError) as e: + read_idat(p) + assert e.value.code == code + + +def test_pair_mismatch(tmp_path) -> None: # type: ignore[no-untyped-def] + a, b = tmp_path / "a.idat", tmp_path / "b.idat" + a.write_bytes(encode_idat(_synth(n=50))) + b.write_bytes(encode_idat(_synth(n=49))) + with pytest.raises(FormatError) as e: + check_pair(read_idat(a), read_idat(b)) + assert e.value.code == "E_IDAT_PAIR_MISMATCH" + + +@pytest.mark.parametrize( + "name,expected", + [ + ("207895980019_R03C01_Grn.idat", ("207895980019", "R03C01", "Grn")), + ("207895980019_r03c01_red.idat.gz", ("207895980019", "R03C01", "Red")), + ("5723646052_R02C02_GRN.IDAT", ("5723646052", "R02C02", "Grn")), + ("4019585376_B_Red.idat", None), + ("notes.txt", None), + ], +) +def test_filename(name, expected) -> None: # type: ignore[no-untyped-def] + assert parse_idat_filename(name) == expected + + +@pytest.mark.parametrize( + "raw,iso", + [ + ("12/27/2023 5:53:34 PM", "2023-12-27T17:53:34"), + ("1/2/2020 12:00:01 AM", "2020-01-02T00:00:01"), + ("2020-01-02 13:14:15", "2020-01-02T13:14:15"), + ("", ""), + ("garbage", ""), + ], +) +def test_run_time(raw, iso) -> None: # type: ignore[no-untyped-def] + assert parse_run_time(raw) == iso + + +def test_duplicate_addresses_rejected(tmp_path) -> None: # type: ignore[no-untyped-def] + from simurg.formats.idat import has_duplicates + + assert not has_duplicates(np.array([1, 5, 9], np.int32)) + assert not has_duplicates(np.array([9, 1, 5], np.int32)) + assert has_duplicates(np.array([9, 1, 9], np.int32)) + s = _synth(n=5) + s.address_ids = np.array([5, 1, 5, 7, 8], np.int32) + p = tmp_path / "dup.idat" + p.write_bytes(encode_idat(s)) + with pytest.raises(FormatError): + read_idat(p) diff --git a/tests/unit/test_samplesheet.py b/tests/unit/test_samplesheet.py new file mode 100644 index 0000000..119a129 --- /dev/null +++ b/tests/unit/test_samplesheet.py @@ -0,0 +1,160 @@ +"""Sample-sheet matrix — spec EDGE-007 (IN-001 … IN-014).""" + +from __future__ import annotations + +from pathlib import Path + +import pytest + +from simurg.core.errors import SimurgError +from simurg.formats.samplesheet import output_stem, parse_samplesheet + +HEAD = "[Header]\nInvestigator Name,X\n[Manifests]\nA,foo.bpm\n[Data]\n" + + +def sheet(tmp_path: Path, body: str, name: str = "s.csv", raw: bytes | None = None) -> Path: + p = tmp_path / name + if raw is not None: + p.write_bytes(raw) + else: + p.write_text(body, encoding="utf-8") + return p + + +def code_of(tmp_path: Path, body: str) -> str: + with pytest.raises(SimurgError) as e: + parse_samplesheet(sheet(tmp_path, body)) + return e.value.code + + +def test_basic_sections_and_passthrough(tmp_path) -> None: # type: ignore[no-untyped-def] + s = parse_samplesheet( + sheet( + tmp_path, + HEAD + "Sample_ID,SentrixBarcode_A,SentrixPosition_A,Gender,Tissue\na,203790110013,r04c01,Female,blood\n", + ) + ) + assert s.rows[0].sample_id == "203790110013_R04C01" + assert s.rows[0].reported_sex == "F" + assert s.passthrough_columns() == ["Sample_ID", "Gender", "Tissue"] + assert s.header == {"Investigator Name": "X"} + + +def test_sectionless(tmp_path) -> None: # type: ignore[no-untyped-def] + s = parse_samplesheet(sheet(tmp_path, "SentrixBarcode_A,SentrixPosition_A\n203790110013,R01C01\n")) + assert len(s.rows) == 1 + + +def test_bom_crlf_quotes(tmp_path) -> None: # type: ignore[no-untyped-def] + raw = ('[Data]\r\nSentrixBarcode_A,SentrixPosition_A,Note\r\n203790110013,R01C01,"a, b"\r\n').encode() + s = parse_samplesheet(sheet(tmp_path, "", raw=raw)) + assert s.rows[0].values["Note"] == "a, b" + + +def test_trailing_empty_columns(tmp_path) -> None: # type: ignore[no-untyped-def] + s = parse_samplesheet(sheet(tmp_path, "[Data]\nSentrixBarcode_A,SentrixPosition_A,,\n203790110013,R01C01,,\n")) + assert s.columns == ["SentrixBarcode_A", "SentrixPosition_A"] + assert [c for c, _ in s.warnings] == ["W_SHEET_EMPTY_COLUMN"] + assert "" not in s.rows[0].values + + +def test_manifests_section_with_values_ignored(tmp_path) -> None: # type: ignore[no-untyped-def] + body = ( + "[Header],,\nDate,1,\n[Manifests],GSA.bpm,\nA,,\n[Data],,\nSentrixBarcode_A,SentrixPosition_A,x\n" + "203790110013,R01C01,1\n" + ) + s = parse_samplesheet(sheet(tmp_path, body)) + assert s.header == {"Date": "1"} and len(s.rows) == 1 + + +@pytest.mark.parametrize( + "body,code", + [ + ("[Data]\nSentrixBarcode_A,SentrixPosition_A,Sample_Name,sample_name\n1,2,3,4\n", "E_SHEET_DUP_COLUMN"), + ("[Data]\nSentrixBarcode_A,SentrixPosition_A,Sex Estimate\n203790110013,R01C01,M\n", "E_SHEET_RESERVED"), + ("[Data]\nSentrixBarcode_A,SentrixPosition_A,call_rate\n203790110013,R01C01,1\n", "E_SHEET_RESERVED"), + ("[Data]\nSentrixBarcode_A,SentrixPosition_A,x\n203790110013,R01C01," + "a" * 501 + "\n", "E_SHEET_TOO_LONG"), + ("[Data]\nSentrixBarcode_A,SentrixPosition_A\n2.04753E+11,R01C01\n", "E_SHEET_BARCODE_FORMAT"), + ("[Data]\nSentrixBarcode_A,SentrixPosition_A\n203790110013,R1C1\n", "E_SHEET_BARCODE_FORMAT"), + ("[Data]\nSentrixBarcode_A;SentrixPosition_A\n203790110013;R01C01\n", "E_SHEET_DELIMITER"), + ('{"samples": []}', "E_SHEET_DELIMITER"), + ("[Data]\nSample_ID,Foo\n1,2\n", "E_SHEET_NO_KEY_COLUMNS"), + ("[Data]\nSentrixBarcode_A,SentrixPosition_A\n", "E_SHEET_EMPTY"), + ("[Header]\nA,B\n", "E_SHEET_EMPTY"), + ( + "[Data]\nSentrixBarcode_A,SentrixPosition_A,Sample_Name\n203790110013,R01C01,a\n203790110013,R02C01,\n", + "E_SHEET_NAME_EMPTY", + ), + ( + "[Data]\nSentrixBarcode_A,SentrixPosition_A,Sample_Name\n203790110013,R01C01,203790110013_R02C01\n" + "203790110013,R02C01,b\n", + "E_SHEET_NAME_COLLISION", + ), + ( + "[Data]\nSentrixBarcode_A,SentrixPosition_A,Sample_Name\n203790110013,R01C01,a\n203790110013,R02C01,a\n", + "E_SHEET_NAME_COLLISION", + ), + ( + "[Data]\nSentrixBarcode_A,SentrixPosition_A\n203790110013,R01C01\n203790110013,r01c01\n", + "E_SHEET_DUP_SAMPLE", + ), + ], +) +def test_error_matrix(tmp_path, body, code) -> None: # type: ignore[no-untyped-def] + assert code_of(tmp_path, body) == code + + +def test_barcode_decimal_repaired(tmp_path) -> None: # type: ignore[no-untyped-def] + s = parse_samplesheet(sheet(tmp_path, "[Data]\nSentrixBarcode_A,SentrixPosition_A\n204753010023.0,R01C01\n")) + assert s.rows[0].barcode == "204753010023" + assert s.warnings[0][0] == "W_SHEET_BARCODE_REPAIRED" + + +def test_scientific_notation_hint(tmp_path) -> None: # type: ignore[no-untyped-def] + with pytest.raises(SimurgError) as e: + parse_samplesheet(sheet(tmp_path, '[Data]\nSentrixBarcode_A,SentrixPosition_A\n"2,04753E+11",R01C01\n')) + assert e.value.hint and "0 decimals" in e.value.hint + + +def test_alias_keys_and_rootfolder(tmp_path) -> None: # type: ignore[no-untyped-def] + s = parse_samplesheet( + sheet(tmp_path, "[Header]\nRootFolder,idats\n[Data]\nbeadChipName,sampleSectionName\n203790110013,R01C01\n") + ) + assert s.barcode_column == "beadChipName" + assert s.root_folder == tmp_path / "idats" + + +def test_sex_columns(tmp_path) -> None: # type: ignore[no-untyped-def] + s = parse_samplesheet( + sheet( + tmp_path, + "[Data]\nSentrixBarcode_A,SentrixPosition_A,Gender,Sex\n" + "203790110013,R01C01,Male,F\n203790110013,R02C01,,\n", + ) + ) + assert [r.reported_sex for r in s.rows] == ["F", "U"] + assert any(c == "W_SHEET_SEX_CONFLICT" for c, _ in s.warnings) + + +def test_turkish_case_folding_is_locale_independent(tmp_path) -> None: # type: ignore[no-untyped-def] # EDGE-001 + s = parse_samplesheet(sheet(tmp_path, "[Data]\nSENTRIXBARCODE_A,SentrixPosition_A,İl\n203790110013,R01C01,x\n")) + assert s.barcode_column == "SENTRIXBARCODE_A" + + +def test_ten_thousand_rows_fast(tmp_path) -> None: # type: ignore[no-untyped-def] # PERF-007 (parse part) + import time + + lines = ["[Data]", "SentrixBarcode_A,SentrixPosition_A,Sample_Name"] + for i in range(10_000): + lines.append(f"{200000000000 + i // 24},R{(i % 24) + 1:02d}C01,s{i}") + p = sheet(tmp_path, "\n".join(lines) + "\n") + t = time.perf_counter() + s = parse_samplesheet(p) + assert len(s.rows) == 10_000 + assert time.perf_counter() - t < 2.0 + + +def test_output_stem() -> None: + assert output_stem("NA 12878/ü") == "NA_12878__" + assert output_stem("..hidden") == "hidden" + assert len(output_stem("x" * 500)) == 200