P1: file parsers, sample sheet, discovery, CLI inspect
Readers for IDAT v3, BPM (locus entries v4-v8), CSV manifest, EGT (data blocks 5-9) and GTC v3-v5 (plus a layout-preserving GTC writer); sample-sheet validation (IN-001..IN-014), IDAT discovery and pairing (IN-020..IN-022), message catalogue and exit codes (CLI-014/016), JSON Lines logging (CLI-007/008), atomic writes (CLI-011), run manifest, and `simurg inspect`. Validation (docs/validation/P1.md): GTC, BPM and EGT readers agree exactly with Illumina BeadArrayFiles on GSA v4 / GDA v1 files; GTCs round-trip byte-identically; BPM and CSV manifests are consistent.
This commit is contained in:
kemal committed
1 parent
705138d70b
commit
e518ea03d5
40 files changed
+4742
No files matched your search
@@ -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 |
|
||||
@@ -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 |
|
||||
@@ -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.
|
||||
@@ -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"
|
||||
@@ -0,0 +1 @@
|
||||
"""SimurgArray cli package."""
|
||||
@@ -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()
|
||||
@@ -0,0 +1 @@
|
||||
"""SimurgArray cn package."""
|
||||
@@ -0,0 +1 @@
|
||||
"""SimurgArray core package."""
|
||||
@@ -0,0 +1,86 @@
|
||||
"""Atomic output writing and overwrite policy (spec CLI-010, CLI-011).
|
||||
|
||||
Every output is written to ``<name>.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
|
||||
@@ -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="")
|
||||
@@ -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()))
|
||||
@@ -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:<path>|profile:<name>|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")
|
||||
@@ -0,0 +1 @@
|
||||
"""SimurgArray formats package."""
|
||||
@@ -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("<b")
|
||||
_U8 = struct.Struct("<B")
|
||||
_I16 = struct.Struct("<h")
|
||||
_U16 = struct.Struct("<H")
|
||||
_I32 = struct.Struct("<i")
|
||||
_U32 = struct.Struct("<I")
|
||||
_I64 = struct.Struct("<q")
|
||||
_F32 = struct.Struct("<f")
|
||||
|
||||
GZIP_MAGIC = b"\x1f\x8b"
|
||||
|
||||
|
||||
def load_bytes(path: str | Path) -> 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)
|
||||
@@ -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("<i")
|
||||
_F32X4 = struct.Struct("<4f")
|
||||
|
||||
|
||||
@dataclass
|
||||
class Manifest:
|
||||
"""Columnar manifest; every list/array has one entry per locus in names order."""
|
||||
|
||||
path: str
|
||||
format: str # "bpm" or "csv"
|
||||
version: int
|
||||
manifest_name: str
|
||||
control_config: str
|
||||
names: list[str]
|
||||
ilmn_ids: list[str]
|
||||
snps: list[str]
|
||||
chroms: list[str]
|
||||
map_infos: npt.NDArray[np.int64]
|
||||
ilmn_strands: list[str]
|
||||
source_strands: list[str]
|
||||
ref_strands: list[str]
|
||||
ploidies: list[str]
|
||||
species: list[str]
|
||||
genome_builds: list[str]
|
||||
sources: list[str]
|
||||
source_versions: list[str]
|
||||
customer_strands: list[str]
|
||||
source_seqs: list[str]
|
||||
top_genomic_seqs: list[str]
|
||||
allele_a_probe_seqs: list[str]
|
||||
allele_b_probe_seqs: list[str]
|
||||
address_a: npt.NDArray[np.int32]
|
||||
address_b: npt.NDArray[np.int32]
|
||||
assay_types: npt.NDArray[np.uint8]
|
||||
exp_clusters: npt.NDArray[np.int8]
|
||||
intensity_only: npt.NDArray[np.bool_]
|
||||
norm_ids: npt.NDArray[np.uint8] # raw per-locus IDs (< 100)
|
||||
entry_versions: npt.NDArray[np.int8]
|
||||
base_fractions: npt.NDArray[np.float32] | None = None # (n, 4): A, C, T, G (entry v >= 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("<i4", n, "indexes")
|
||||
names = [r.string("locus name") for _ in range(n)]
|
||||
norm_ids = np.array(r.array("u1", n, "normalization ids"), dtype=np.uint8)
|
||||
if (norm_ids >= 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)
|
||||
@@ -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
|
||||
@@ -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"
|
||||
@@ -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 ``<barcode>_<pos>_{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
|
||||
@@ -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", "<i4", (3,)),
|
||||
("r_dev", "<f4", (3,)),
|
||||
("r_mean", "<f4", (3,)),
|
||||
("theta_dev", "<f4", (3,)),
|
||||
("theta_mean", "<f4", (3,)),
|
||||
]
|
||||
_SCORE_DTYPE = np.dtype([("cluster_sep", "<f4"), ("total_score", "<f4"), ("original_score", "<f4"), ("edited", "u1")])
|
||||
|
||||
|
||||
def _record_dtype(extra: int) -> np.dtype[Any]:
|
||||
fields: list[Any] = list(_BASE_FIELDS)
|
||||
if extra:
|
||||
fields.append(("extra", "<f4", (extra,)))
|
||||
return np.dtype(fields)
|
||||
|
||||
|
||||
def _extra_candidates(dbv: int) -> 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("<i4", nrec, "addresses"), dtype=np.int32)
|
||||
counts = r.array("<i4", nrec * 3, "cluster counts").reshape(nrec, 3)
|
||||
n = np.array(rec["n"], dtype=np.int32)
|
||||
if not np.array_equal(counts, n):
|
||||
raise r.fail(f"cluster counts do not match record N (assuming {extra} extra floats per record)")
|
||||
trailer = None
|
||||
if dbv == 9 and r.remaining() >= 4 * nrec:
|
||||
trailer = np.array(r.array("<f4", nrec, "trailer floats"), dtype=np.float32)
|
||||
if len(set(names)) != nrec:
|
||||
raise r.fail("duplicate locus names")
|
||||
extra_arr = np.array(rec["extra"], dtype=np.float32) if extra else None
|
||||
if extra_arr is not None and dbv >= 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,
|
||||
)
|
||||
@@ -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", "<i4"),
|
||||
("offset_x", "<f4"),
|
||||
("offset_y", "<f4"),
|
||||
("scale_x", "<f4"),
|
||||
("scale_y", "<f4"),
|
||||
("shear", "<f4"),
|
||||
("theta", "<f4"),
|
||||
("reserved", "<f4", (6,)),
|
||||
]
|
||||
)
|
||||
assert XFORM_DTYPE.itemsize == 52
|
||||
SAMPLE_DATA_DTYPE = np.dtype(
|
||||
[("p50gc", "<f4"), ("num_calls", "<i4"), ("num_no_calls", "<i4"), ("num_intensity_only", "<i4")]
|
||||
)
|
||||
|
||||
|
||||
@dataclass
|
||||
class ScannerData:
|
||||
name: str = ""
|
||||
pmt_green: bytes = b"\x00\x00\x00\x00" # 4 raw bytes: int32 in BeadArrayFiles, float in gtc2vcf
|
||||
pmt_red: bytes = b"\x00\x00\x00\x00"
|
||||
version: str = ""
|
||||
user: str = ""
|
||||
|
||||
@property
|
||||
def pmt_green_int(self) -> 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, "<u2")
|
||||
if at(CONTROLS_Y):
|
||||
g.controls_y = arr(CONTROLS_Y, "<u2")
|
||||
if at(RAW_X):
|
||||
g.raw_x = arr(RAW_X, "<u2", n)
|
||||
if at(RAW_Y):
|
||||
g.raw_y = arr(RAW_Y, "<u2", n)
|
||||
if at(GENOTYPES):
|
||||
g.genotypes = arr(GENOTYPES, "u1", n)
|
||||
if at(BASE_CALLS):
|
||||
cnt = r.i32("base calls count")
|
||||
if cnt != n:
|
||||
raise r.fail(f"base_calls has {cnt} entries, expected {n}")
|
||||
g.base_calls = np.array(r.array("u1", 2 * n, "base calls")).reshape(n, 2)
|
||||
if at(GENOTYPE_SCORES):
|
||||
g.genotype_scores = arr(GENOTYPE_SCORES, "<f4", n)
|
||||
if at(SCANNER_DATA):
|
||||
g.scanner_data = ScannerData(
|
||||
name=r.string("scanner name"),
|
||||
pmt_green=r.raw(4, "pmt green"),
|
||||
pmt_red=r.raw(4, "pmt red"),
|
||||
version=r.string("scanner version"),
|
||||
user=r.string("imaging user"),
|
||||
)
|
||||
if at(CALL_RATE):
|
||||
g.call_rate = r.f32("call rate")
|
||||
if at(GENDER):
|
||||
g.gender = r.raw(1, "gender")
|
||||
if at(LOGR_DEV):
|
||||
g.logr_dev = r.f32("logr dev")
|
||||
if at(GC10):
|
||||
g.gc10 = r.f32("gc10")
|
||||
if at(DX):
|
||||
g.dx = r.i32("dx")
|
||||
if at(SAMPLE_DATA):
|
||||
g.sample_data = np.array(r.array(SAMPLE_DATA_DTYPE, 1, "sample data"))[0]
|
||||
if at(B_ALLELE_FREQS):
|
||||
g.b_allele_freqs = arr(B_ALLELE_FREQS, "<f4", n)
|
||||
if at(LOGR_RATIOS):
|
||||
g.logr_ratios = arr(LOGR_RATIOS, "<f4", n)
|
||||
if at(PERCENTILES_X):
|
||||
g.percentiles_x = np.array(r.array("<u2", 3, "percentiles x"))
|
||||
if at(PERCENTILES_Y):
|
||||
g.percentiles_y = np.array(r.array("<u2", 3, "percentiles y"))
|
||||
unknown = [t for t in toc if t not in TOC_NAMES]
|
||||
if unknown:
|
||||
raise r.fail(f"unknown TOC ids {unknown}")
|
||||
g.data_order = sorted((t for t in toc if t not in INLINE_IDS), key=lambda t: toc[t])
|
||||
if check_write_complete:
|
||||
required = {3: SAMPLE_DATA, 4: LOGR_RATIOS, 5: SLIDE_IDENTIFIER}[version]
|
||||
if required not in toc:
|
||||
raise r.fail(f"GTC v{version} is incomplete (missing {TOC_NAMES[required]})")
|
||||
return g
|
||||
|
||||
|
||||
def _as_i32(v: int) -> 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]), "<u2")
|
||||
elif tid in (RAW_X, RAW_Y):
|
||||
_check_len(g, tid, n)
|
||||
w.counted_array(getattr(g, TOC_NAMES[tid]), "<u2")
|
||||
elif tid == GENOTYPES:
|
||||
_check_len(g, tid, n)
|
||||
w.counted_array(g.genotypes, "u1")
|
||||
elif tid == BASE_CALLS:
|
||||
bc = np.asarray(g.base_calls, dtype=np.uint8).reshape(-1, 2)
|
||||
if len(bc) != n:
|
||||
raise ValueError(f"base_calls has {len(bc)} rows, expected {n}")
|
||||
w.i32(n)
|
||||
w.write(bc.tobytes())
|
||||
elif tid in (GENOTYPE_SCORES, B_ALLELE_FREQS, LOGR_RATIOS):
|
||||
_check_len(g, tid, n)
|
||||
w.counted_array(getattr(g, TOC_NAMES[tid]), "<f4")
|
||||
elif tid == SCANNER_DATA:
|
||||
s = g.scanner_data
|
||||
w.string(s.name)
|
||||
w.write(s.pmt_green)
|
||||
w.write(s.pmt_red)
|
||||
w.string(s.version)
|
||||
w.string(s.user)
|
||||
elif tid in (CALL_RATE, LOGR_DEV, GC10):
|
||||
w.f32(getattr(g, TOC_NAMES[tid]))
|
||||
elif tid == GENDER:
|
||||
w.write(g.gender[:1] or b"U")
|
||||
elif tid == DX:
|
||||
w.i32(g.dx)
|
||||
elif tid == SAMPLE_DATA:
|
||||
w.write(np.asarray(g.sample_data, dtype=SAMPLE_DATA_DTYPE).tobytes())
|
||||
elif tid in (PERCENTILES_X, PERCENTILES_Y):
|
||||
w.array(getattr(g, TOC_NAMES[tid]), "<u2")
|
||||
else:
|
||||
raise ValueError(f"cannot write TOC id {tid}")
|
||||
|
||||
|
||||
def _check_len(g: Gtc, tid: int, n: int) -> 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
|
||||
@@ -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<barcode>[^_/\\]+)_(?P<position>R\d{2}C\d{2})_(?P<channel>Grn|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:
|
||||
"""``<barcode>_<RxxCxx>_<Grn|Red>.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="<i4")
|
||||
idat = Idat(
|
||||
path=spath,
|
||||
version=version,
|
||||
toc=dict(toc),
|
||||
n_addresses=n,
|
||||
address_ids=empty_i32,
|
||||
mean=np.empty(0, dtype="<u2"),
|
||||
sd=np.empty(0, dtype="<u2"),
|
||||
n_beads=np.empty(0, dtype="u1"),
|
||||
)
|
||||
|
||||
if load_arrays:
|
||||
r.seek(toc[ILLUMINA_ID])
|
||||
idat.address_ids = r.array("<i4", n, "address ids")
|
||||
r.seek(toc[SD])
|
||||
idat.sd = r.array("<u2", n, "SD")
|
||||
r.seek(toc[MEAN])
|
||||
idat.mean = r.array("<u2", n, "Mean")
|
||||
r.seek(toc[NBEADS])
|
||||
idat.n_beads = r.array("u1", n, "NBeads")
|
||||
if has_duplicates(idat.address_ids):
|
||||
raise r.fail("duplicate address IDs")
|
||||
|
||||
offsets_sorted = sorted(set(toc.values()) | {r.size})
|
||||
|
||||
def section_end(off: int) -> 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("<i4" if itemsize == 4 else "u1", count, "mid block")
|
||||
elif code in (N_SNPS, ILLUMINA_ID, SD, MEAN, NBEADS):
|
||||
continue
|
||||
else:
|
||||
idat.unknown_fields[code] = off
|
||||
return idat
|
||||
|
||||
|
||||
def has_duplicates(ids: npt.NDArray[np.int32]) -> 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}")
|
||||
@@ -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]
|
||||
@@ -0,0 +1 @@
|
||||
"""SimurgArray genotype package."""
|
||||
@@ -0,0 +1 @@
|
||||
"""SimurgArray harmonize package."""
|
||||
@@ -0,0 +1 @@
|
||||
"""SimurgArray pgx package."""
|
||||
@@ -0,0 +1 @@
|
||||
"""SimurgArray pipeline package."""
|
||||
@@ -0,0 +1 @@
|
||||
"""SimurgArray qc package."""
|
||||
@@ -0,0 +1 @@
|
||||
"""SimurgArray signal package."""
|
||||
@@ -0,0 +1 @@
|
||||
"""SimurgArray synth package."""
|
||||
@@ -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 = "<i4" if s.mid_block_itemsize == 4 else "u1"
|
||||
sec(idatfmt.N_SNPS, ("i32", n))
|
||||
sec(idatfmt.ILLUMINA_ID, ("<i4", s.address_ids))
|
||||
sec(idatfmt.SD, ("<u2", sd))
|
||||
sec(idatfmt.MEAN, ("<u2", s.mean))
|
||||
sec(idatfmt.NBEADS, ("u1", nb))
|
||||
sec(idatfmt.MID_BLOCK, ("i32", n), (mid_dtype, s.address_ids))
|
||||
run: list[tuple[str, object]] = [("i32", len(s.run_info))]
|
||||
run += [("str", part) for entry in s.run_info for part in entry]
|
||||
sec(idatfmt.RUN_INFO, *run)
|
||||
sec(idatfmt.RED_GREEN, ("raw", b"\x00\x00\x00\x00"))
|
||||
for code, val in (
|
||||
(idatfmt.MANIFEST, s.manifest),
|
||||
(idatfmt.BARCODE, s.barcode),
|
||||
(idatfmt.CHIP_TYPE, s.chip_type),
|
||||
(idatfmt.POSITION, s.position),
|
||||
(idatfmt.BEADSET, ""),
|
||||
(idatfmt.SAMPLE_ID, s.sample_id),
|
||||
(idatfmt.DESCRIPTION, ""),
|
||||
(idatfmt.PLATE, s.plate),
|
||||
(idatfmt.WELL, s.well),
|
||||
):
|
||||
sec(code, ("str", val))
|
||||
sec(idatfmt.SAMPLE_COUNT, ("i32", 1))
|
||||
|
||||
header = 4 + 8 + 4 + 10 * len(sections)
|
||||
out = BinaryWriter()
|
||||
out.write(b"IDAT")
|
||||
out.write(int(s.version).to_bytes(8, "little", signed=True))
|
||||
out.i32(len(sections))
|
||||
off = header
|
||||
for code, data in sections:
|
||||
out.u16(code)
|
||||
out.write(off.to_bytes(8, "little", signed=True))
|
||||
off += len(data)
|
||||
for _, data in sections:
|
||||
out.write(data)
|
||||
return out.getvalue()
|
||||
|
||||
|
||||
@dataclass
|
||||
class SynthLocus:
|
||||
name: str
|
||||
chrom: str = "1"
|
||||
pos: int = 1000
|
||||
snp: str = "[A/G]"
|
||||
address_a: int = 1
|
||||
address_b: int = 0
|
||||
assay_type: int = 0
|
||||
norm_id: int = 0
|
||||
ilmn_strand: str = "TOP"
|
||||
ref_strand: str = "+"
|
||||
source_strand: str = "TOP"
|
||||
intensity_only: int = 0
|
||||
exp_clusters: int = 3
|
||||
source_seq: str = ""
|
||||
probe_a: str = ""
|
||||
probe_b: str = ""
|
||||
ilmn_id: str | None = None
|
||||
|
||||
|
||||
def encode_bpm(
|
||||
loci: Sequence[SynthLocus],
|
||||
*,
|
||||
manifest_name: str = "SYNTH-1v1_A1",
|
||||
entry_version: int = 8,
|
||||
bpm_version: int = 4,
|
||||
control_config: str = "",
|
||||
trailer: Sequence[str] = ("Illumina, Inc.",),
|
||||
entry_order: Sequence[int] | None = None,
|
||||
) -> 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), "<i4")
|
||||
for lc in loci:
|
||||
w.string(lc.name)
|
||||
w.array([lc.norm_id for lc in loci], "u1")
|
||||
order = list(entry_order) if entry_order is not None else list(range(n))
|
||||
for i in order:
|
||||
lc = loci[i]
|
||||
w.i32(entry_version)
|
||||
w.string(lc.ilmn_id or f"{lc.name}-0_T_F_1")
|
||||
w.string(lc.name)
|
||||
for _ in range(3):
|
||||
w.string("")
|
||||
w.i32(i + 1)
|
||||
w.string("")
|
||||
w.string(lc.ilmn_strand)
|
||||
w.string(lc.snp)
|
||||
w.string(lc.chrom)
|
||||
w.string("diploid")
|
||||
w.string("Homo sapiens")
|
||||
w.string(str(lc.pos))
|
||||
w.string("") # TopGenomicSeq (v4 only)
|
||||
w.string("") # customer strand
|
||||
w.i32(lc.address_a)
|
||||
w.i32(lc.address_b)
|
||||
w.string(lc.probe_a if entry_version == 4 else "")
|
||||
w.string(lc.probe_b if entry_version == 4 else "")
|
||||
w.string("38")
|
||||
w.string("synthetic")
|
||||
w.string("0")
|
||||
w.string(lc.source_strand)
|
||||
w.string(lc.source_seq if entry_version == 4 else "")
|
||||
if entry_version >= 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], "<i4")
|
||||
w.array(rd[i], "<f4")
|
||||
w.array(rm[i], "<f4")
|
||||
w.array(td[i], "<f4")
|
||||
w.array(tm[i], "<f4")
|
||||
if data_block_version >= 7:
|
||||
w.f32(float(thr[i]))
|
||||
w.array(np.full(14, np.nan), "<f4")
|
||||
elif data_block_version >= 6:
|
||||
w.array(np.zeros(15), "<f4")
|
||||
elif data_block_version < 5:
|
||||
w.array(np.zeros(2), "<f4")
|
||||
for i in range(k):
|
||||
w.f32(0.5)
|
||||
w.f32(float(gt[i]))
|
||||
w.f32(float(gt[i]))
|
||||
w.u8(0)
|
||||
for _ in range(k):
|
||||
w.string("aa_ab_bb")
|
||||
for nm in names:
|
||||
w.string(nm)
|
||||
w.array(np.arange(1, k + 1), "<i4")
|
||||
counts = nn.copy()
|
||||
if corrupt_counts:
|
||||
counts[0, 0] += 1
|
||||
for i in range(k):
|
||||
w.array(counts[i], "<i4")
|
||||
if data_block_version == 9:
|
||||
w.array(np.zeros(k), "<f4")
|
||||
return w.getvalue()
|
||||
Whitespace-only changes.
@@ -0,0 +1,45 @@
|
||||
"""Shared test configuration.
|
||||
|
||||
Tests marked ``data`` need files fetched by ``scripts/fetch_data.sh p0``; tests marked ``oracle`` need
|
||||
reference repositories under ``refs/`` (``scripts/bootstrap_refs.sh``). Both are skipped when absent.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import sys
|
||||
from pathlib import Path
|
||||
|
||||
import pytest
|
||||
|
||||
ROOT = Path(__file__).resolve().parents[1]
|
||||
DATA = ROOT / "data" / "p0"
|
||||
REFS = ROOT / "refs"
|
||||
|
||||
GSA = DATA / "gsa_v4_pgx"
|
||||
GSA_BPM = GSA / "GSA-PGx-48v4-0_20079540_E2_bpm" / "GSA-PGx-48v4-0_20079540_E2.bpm"
|
||||
GSA_CSV = GSA / "GSA-PGx-48v4-0_20079540_E2_csv" / "GSA-PGx-48v4-0_20079540_E2.csv"
|
||||
GSA_EGT = GSA / "GSA-PGx-48v4-0_20079540_E2_ClusterFile" / "GSA-PGx-48v4-0_20079540_E2_ClusterFile.egt"
|
||||
GSA_IDATS = GSA / "idats_Demo_36" / "idats_Demo_36"
|
||||
GDA_EGT = DATA / "gda_v1" / "infinium-global-diversity-array-8-v1-0-D1-cluster-file" / "GDA-8v1-0_D1_ClusterFile.egt"
|
||||
GDA_BPM = DATA / "gda_v1" / "infinium-global-diversity-array-8-v1-0-D1-manifest-file-bpm" / "GDA-8v1-0_D1.bpm"
|
||||
GDA_SHEET = DATA / "gda_v1" / "infinium-global-diversity-array-8-v1-0_D1_SampleSheet_Demo_36.csv"
|
||||
BIOC_IDATS = DATA / "bioc" / "IlluminaDataTestFiles_1.50.0" / "IlluminaDataTestFiles" / "inst" / "extdata" / "idat"
|
||||
ORACLE_GTCS = DATA / "oracle" / "gsa_gtc_gtc2vcf"
|
||||
|
||||
|
||||
def need(*paths: Path) -> 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
|
||||
@@ -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"}
|
||||
@@ -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="<i4").tobytes(), "E_IDAT_FORMAT")
|
||||
assert r.array("<i4", 4).tolist() == [0, 1, 2, 3]
|
||||
r.seek(0)
|
||||
with pytest.raises(FormatError):
|
||||
r.array("<i4", 5)
|
||||
with pytest.raises(FormatError):
|
||||
r.seek(100)
|
||||
|
||||
|
||||
def test_load_bytes_gzip(tmp_path) -> 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""
|
||||
@@ -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)
|
||||
@@ -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
|
||||
@@ -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)
|
||||
@@ -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)
|
||||
@@ -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("<q", 4) + b[12:], "E_IDAT_FORMAT"), # unsupported version
|
||||
(lambda b: b[:4] + struct.pack("<q", 1 | (7 << 32)) + b[12:], "E_IDAT_ENCRYPTED"), # encrypted v1
|
||||
(lambda b: b[: len(b) // 2], "E_IDAT_FORMAT"), # truncated
|
||||
(lambda b: b[:18] + struct.pack("<q", 10**9) + b[26:], "E_IDAT_FORMAT"), # TOC offset outside file
|
||||
(lambda b: b[:8], "E_IDAT_FORMAT"), # tiny
|
||||
],
|
||||
)
|
||||
def test_malformed(tmp_path, mutate, code) -> 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)
|
||||
@@ -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
|
||||
Reference in new issue
Block a user