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
fb1ff51b76
commit
4c0ea5e04f
9 files changed
+970
-2
No files matched your search
@@ -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
|
||||
Reference in new issue
Block a user