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