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:
kemal committed
1 parent
a7a235beb5
commit
87f9b3499c
9 files changed
+970
-2
No files matched your search
+50
-1
@@ -1,6 +1,6 @@
|
||||
"""``simurg`` command-line interface (spec §2, plan §14).
|
||||
|
||||
Phase P1 provides ``simurg version`` and ``simurg inspect``. Later phases add the other commands of CLI-001.
|
||||
Available: ``version``, ``inspect`` (P1), ``ingest`` (P2). Later phases add the other commands of CLI-001.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
@@ -203,5 +203,54 @@ def inspect_cmd(
|
||||
run_guarded(rep, body)
|
||||
|
||||
|
||||
# --------------------------------------------------------------------------------------------------------
|
||||
@app.command("ingest")
|
||||
def ingest_cmd(
|
||||
bpm: Annotated[Path, typer.Option("--bpm", help="BPM manifest (required).")],
|
||||
out: Annotated[Path, typer.Option("--out", help="Output folder (required).")],
|
||||
idats: Annotated[Path | None, typer.Option("--idats", help="Folder searched recursively for IDAT pairs.")] = None,
|
||||
sheet: Annotated[Path | None, typer.Option("--sheet", help="Sample sheet (CSV).")] = None,
|
||||
csv: Annotated[Path | None, typer.Option("--csv", help="CSV manifest, checked against the BPM.")] = None,
|
||||
threads: Annotated[int, typer.Option("--threads", min=0, max=1024, help="Worker processes (0 = auto).")] = 0,
|
||||
force: Annotated[bool, typer.Option("--force", help="Replace existing outputs.")] = False,
|
||||
resume: Annotated[bool, typer.Option("--resume", help="Keep outputs whose inputs are unchanged.")] = False,
|
||||
layout: Annotated[str, typer.Option("--layout", help="nested|flat")] = "nested",
|
||||
on_duplicate_idat: Annotated[str, typer.Option("--on-duplicate-idat", help="error|first|newest")] = "error",
|
||||
log_json: Annotated[bool, typer.Option("--log-json", "--json-log", help="JSON Lines log on stderr.")] = False,
|
||||
debug: Annotated[bool, typer.Option("--debug", help="Debug logging with tracebacks.")] = False,
|
||||
log_level: Annotated[str, typer.Option("--log-level", help="debug|info|warning|error")] = "info",
|
||||
quiet: Annotated[bool, typer.Option("--quiet", "--no-progress", help="Suppress info messages.")] = False,
|
||||
) -> None:
|
||||
"""IDAT pairs + BPM → canonical per-sample raw intensities (Parquet). Phase P2."""
|
||||
from simurg.pipeline.ingest import IngestOptions, run_ingest
|
||||
|
||||
if layout not in ("nested", "flat"):
|
||||
raise typer.BadParameter("--layout must be nested or flat")
|
||||
if on_duplicate_idat not in ("error", "first", "newest"):
|
||||
raise typer.BadParameter("--on-duplicate-idat must be error, first or newest")
|
||||
rep = Reporter(
|
||||
command="ingest",
|
||||
level="debug" if debug else log_level,
|
||||
json_stderr=log_json,
|
||||
debug=debug,
|
||||
quiet=quiet,
|
||||
log_file=out / "logs" / "ingest.log",
|
||||
)
|
||||
opts = IngestOptions(
|
||||
out=out,
|
||||
bpm=bpm,
|
||||
idats=idats,
|
||||
sheet=sheet,
|
||||
csv=csv,
|
||||
threads=threads,
|
||||
force=force,
|
||||
resume=resume,
|
||||
layout=layout,
|
||||
on_duplicate_idat=on_duplicate_idat,
|
||||
argv=sys.argv[1:],
|
||||
)
|
||||
run_guarded(rep, lambda: run_ingest(opts, rep))
|
||||
|
||||
|
||||
if __name__ == "__main__": # pragma: no cover
|
||||
app()
|
||||
@@ -72,6 +72,9 @@ CATALOGUE: dict[str, MessageDef] = dict(
|
||||
_m("E_IDAT_ENCRYPTED", "error", "Encrypted (version 1) IDATs are not supported", "IN-023"),
|
||||
_m("E_IDAT_PAIR_MISMATCH", "error", "Red and Grn IDATs of a sample do not match", "IN-023"),
|
||||
_m("W_IDAT_NAME_MISMATCH", "warning", "IDAT barcode/position disagrees with its file name", "IN-023"),
|
||||
_m("E_IDAT_MANIFEST_MISMATCH", "error", "IDATs and manifest belong to different products", "IN-027"),
|
||||
_m("W_IDAT_LOW_ADDRESS_COVERAGE", "warning", "Some manifest addresses are absent from the IDAT", "IN-027"),
|
||||
_m("W_IDAT_MANIFEST_NAME", "warning", "IDAT names a different manifest than the BPM", "IN-027"),
|
||||
_m("W_IDAT_NO_SCAN_TIME", "warning", "IDAT has no scan entry in its run info", "IN-024"),
|
||||
_m("I_UNLISTED_IDATS", "info", "IDAT pairs not listed in the sample sheet are ignored", "IN-016"),
|
||||
_m("E_DUP_GTC", "error", "Several GTC files resolve to the same sample", "IN-030"),
|
||||
|
||||
@@ -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
|
||||
|
||||
import struct
|
||||
from dataclasses import dataclass, field
|
||||
from dataclasses import dataclass, field, replace
|
||||
from pathlib import Path
|
||||
from typing import Any
|
||||
|
||||
@@ -103,6 +103,37 @@ class Manifest:
|
||||
def name_index(self) -> dict[str, int]:
|
||||
return {n: i for i, n in enumerate(self.names)}
|
||||
|
||||
def lite(self) -> Manifest:
|
||||
"""Copy keeping only what intensity assembly needs (cheap to pickle to worker processes)."""
|
||||
n = self.num_loci
|
||||
empty = [""] * 0
|
||||
return replace(
|
||||
self,
|
||||
names=list(self.names),
|
||||
ilmn_ids=empty,
|
||||
snps=empty,
|
||||
chroms=empty,
|
||||
ilmn_strands=empty,
|
||||
source_strands=empty,
|
||||
ref_strands=empty,
|
||||
ploidies=empty,
|
||||
species=empty,
|
||||
genome_builds=empty,
|
||||
sources=empty,
|
||||
source_versions=empty,
|
||||
customer_strands=empty,
|
||||
source_seqs=empty,
|
||||
top_genomic_seqs=empty,
|
||||
allele_a_probe_seqs=empty,
|
||||
allele_b_probe_seqs=empty,
|
||||
map_infos=np.zeros(0, np.int64),
|
||||
base_fractions=None,
|
||||
header_lines=[],
|
||||
heading={},
|
||||
exp_clusters=self.exp_clusters[:n],
|
||||
intensity_only=self.intensity_only[:n],
|
||||
)
|
||||
|
||||
def summary(self) -> dict[str, Any]:
|
||||
at = self.assay_types
|
||||
return {
|
||||
|
||||
@@ -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)
|
||||
Reference in new issue
Block a user