P1: file parsers, sample sheet, discovery, CLI inspect
Readers for IDAT v3, BPM (locus entries v4-v8), CSV manifest, EGT (data blocks 5-9) and GTC v3-v5 (plus a layout-preserving GTC writer); sample-sheet validation (IN-001..IN-014), IDAT discovery and pairing (IN-020..IN-022), message catalogue and exit codes (CLI-014/016), JSON Lines logging (CLI-007/008), atomic writes (CLI-011), run manifest, and `simurg inspect`. Validation (docs/validation/P1.md): GTC, BPM and EGT readers agree exactly with Illumina BeadArrayFiles on GSA v4 / GDA v1 files; GTCs round-trip byte-identically; BPM and CSV manifests are consistent.
This commit is contained in:
kemal committed
1 parent
0a421dffc9
commit
5a20e7131c
40 files changed
+4742
No files matched your search
@@ -0,0 +1 @@
|
||||
"""SimurgArray formats package."""
|
||||
@@ -0,0 +1,208 @@
|
||||
"""Little-endian binary reading/writing helpers shared by the IDAT, BPM, EGT and GTC codecs.
|
||||
|
||||
All Illumina binary formats are little-endian. Strings are .NET ``BinaryWriter`` strings: a 7-bit varint
|
||||
(LEB128) length prefix followed by UTF-8 bytes (plan §6; BeadArrayFiles ``read_string``, BSD-2).
|
||||
Decoding is tolerant (spec IN-033): invalid UTF-8 becomes U+FFFD, while a length running past the end of the
|
||||
buffer is a :class:`~simurg.core.errors.FormatError`, never an uncaught exception.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import gzip
|
||||
import io
|
||||
import mmap
|
||||
import struct
|
||||
from pathlib import Path
|
||||
from typing import Any
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
from simurg.core.errors import FormatError
|
||||
|
||||
_I8 = struct.Struct("<b")
|
||||
_U8 = struct.Struct("<B")
|
||||
_I16 = struct.Struct("<h")
|
||||
_U16 = struct.Struct("<H")
|
||||
_I32 = struct.Struct("<i")
|
||||
_U32 = struct.Struct("<I")
|
||||
_I64 = struct.Struct("<q")
|
||||
_F32 = struct.Struct("<f")
|
||||
|
||||
GZIP_MAGIC = b"\x1f\x8b"
|
||||
|
||||
|
||||
def load_bytes(path: str | Path) -> bytes | mmap.mmap:
|
||||
"""Return the file content: memory-mapped when plain, decompressed when gzip."""
|
||||
p = Path(path)
|
||||
with p.open("rb") as fh:
|
||||
head = fh.read(2)
|
||||
if head == GZIP_MAGIC:
|
||||
fh.seek(0)
|
||||
with gzip.GzipFile(fileobj=fh) as gz:
|
||||
return gz.read()
|
||||
fh.seek(0, io.SEEK_END)
|
||||
if fh.tell() == 0:
|
||||
return b""
|
||||
return mmap.mmap(fh.fileno(), 0, access=mmap.ACCESS_READ)
|
||||
|
||||
|
||||
class BinaryReader:
|
||||
"""Bounds-checked sequential reader over an immutable buffer.
|
||||
|
||||
``err_code`` is the catalogue code raised on malformed input (e.g. ``E_IDAT_FORMAT``).
|
||||
"""
|
||||
|
||||
__slots__ = ("_buf", "_view", "err_code", "path", "pos", "size")
|
||||
|
||||
def __init__(self, buf: bytes | mmap.mmap | memoryview, err_code: str, path: str | None = None) -> None:
|
||||
self._buf = buf
|
||||
self._view = memoryview(buf)
|
||||
self.size = len(self._view)
|
||||
self.pos = 0
|
||||
self.err_code = err_code
|
||||
self.path = path
|
||||
|
||||
# -- errors ------------------------------------------------------------------------------------
|
||||
def fail(self, message: str) -> FormatError:
|
||||
where = f"{self.path}: " if self.path else ""
|
||||
return FormatError(self.err_code, f"{where}{message}", path=self.path)
|
||||
|
||||
def _need(self, n: int, what: str) -> int:
|
||||
if n < 0 or self.pos + n > self.size:
|
||||
raise self.fail(f"unexpected end of file reading {what} at offset {self.pos} (need {n} bytes)")
|
||||
start = self.pos
|
||||
self.pos += n
|
||||
return start
|
||||
|
||||
# -- positioning -------------------------------------------------------------------------------
|
||||
def seek(self, offset: int, what: str = "section") -> None:
|
||||
if offset < 0 or offset > self.size:
|
||||
raise self.fail(f"{what} offset {offset} outside file (size {self.size})")
|
||||
self.pos = offset
|
||||
|
||||
def remaining(self) -> int:
|
||||
return self.size - self.pos
|
||||
|
||||
# -- scalars -----------------------------------------------------------------------------------
|
||||
def _unpack(self, s: struct.Struct, what: str) -> Any:
|
||||
start = self._need(s.size, what)
|
||||
return s.unpack_from(self._view, start)[0]
|
||||
|
||||
def i8(self, what: str = "int8") -> int:
|
||||
return int(self._unpack(_I8, what))
|
||||
|
||||
def u8(self, what: str = "uint8") -> int:
|
||||
return int(self._unpack(_U8, what))
|
||||
|
||||
def i16(self, what: str = "int16") -> int:
|
||||
return int(self._unpack(_I16, what))
|
||||
|
||||
def u16(self, what: str = "uint16") -> int:
|
||||
return int(self._unpack(_U16, what))
|
||||
|
||||
def i32(self, what: str = "int32") -> int:
|
||||
return int(self._unpack(_I32, what))
|
||||
|
||||
def u32(self, what: str = "uint32") -> int:
|
||||
return int(self._unpack(_U32, what))
|
||||
|
||||
def i64(self, what: str = "int64") -> int:
|
||||
return int(self._unpack(_I64, what))
|
||||
|
||||
def f32(self, what: str = "float32") -> float:
|
||||
return float(self._unpack(_F32, what))
|
||||
|
||||
def raw(self, n: int, what: str = "bytes") -> bytes:
|
||||
start = self._need(n, what)
|
||||
return bytes(self._view[start : start + n])
|
||||
|
||||
# -- strings -----------------------------------------------------------------------------------
|
||||
def varint_length(self, what: str = "string length") -> int:
|
||||
""".NET 7-bit encoded int (at most 5 bytes)."""
|
||||
result = 0
|
||||
shift = 0
|
||||
for _ in range(5):
|
||||
b = self.u8(what)
|
||||
result |= (b & 0x7F) << shift
|
||||
if not b & 0x80:
|
||||
return result
|
||||
shift += 7
|
||||
raise self.fail(f"invalid 7-bit length prefix for {what} at offset {self.pos}")
|
||||
|
||||
def string(self, what: str = "string") -> str:
|
||||
n = self.varint_length(what)
|
||||
start = self._need(n, what)
|
||||
return bytes(self._view[start : start + n]).decode("utf-8", errors="replace")
|
||||
|
||||
def skip_string(self, what: str = "string") -> None:
|
||||
n = self.varint_length(what)
|
||||
self._need(n, what)
|
||||
|
||||
# -- arrays ------------------------------------------------------------------------------------
|
||||
def array(self, dtype: npt.DTypeLike, count: int, what: str = "array") -> npt.NDArray[Any]:
|
||||
"""Zero-copy little-endian array view (read-only)."""
|
||||
dt = np.dtype(dtype).newbyteorder("<")
|
||||
if count < 0:
|
||||
raise self.fail(f"negative element count {count} for {what}")
|
||||
start = self._need(count * dt.itemsize, what)
|
||||
return np.frombuffer(self._view, dtype=dt, count=count, offset=start)
|
||||
|
||||
def counted_array(self, dtype: npt.DTypeLike, what: str = "array") -> npt.NDArray[Any]:
|
||||
"""int32 element count followed by the elements."""
|
||||
return self.array(dtype, self.i32(f"{what} count"), what)
|
||||
|
||||
|
||||
class BinaryWriter:
|
||||
"""Little-endian writer into an in-memory buffer."""
|
||||
|
||||
def __init__(self) -> None:
|
||||
self._buf = io.BytesIO()
|
||||
|
||||
def tell(self) -> int:
|
||||
return self._buf.tell()
|
||||
|
||||
def getvalue(self) -> bytes:
|
||||
return self._buf.getvalue()
|
||||
|
||||
def write(self, b: bytes) -> None:
|
||||
self._buf.write(b)
|
||||
|
||||
def u8(self, v: int) -> None:
|
||||
self._buf.write(_U8.pack(v))
|
||||
|
||||
def i16(self, v: int) -> None:
|
||||
self._buf.write(_I16.pack(v))
|
||||
|
||||
def u16(self, v: int) -> None:
|
||||
self._buf.write(_U16.pack(v))
|
||||
|
||||
def i32(self, v: int) -> None:
|
||||
self._buf.write(_I32.pack(v))
|
||||
|
||||
def u32(self, v: int) -> None:
|
||||
self._buf.write(_U32.pack(v))
|
||||
|
||||
def f32(self, v: float) -> None:
|
||||
self._buf.write(_F32.pack(v))
|
||||
|
||||
def string(self, s: str | bytes) -> None:
|
||||
data = s.encode("utf-8") if isinstance(s, str) else s
|
||||
n = len(data)
|
||||
while True:
|
||||
b = n & 0x7F
|
||||
n >>= 7
|
||||
if n:
|
||||
self._buf.write(bytes([b | 0x80]))
|
||||
else:
|
||||
self._buf.write(bytes([b]))
|
||||
break
|
||||
self._buf.write(data)
|
||||
|
||||
def array(self, arr: npt.ArrayLike, dtype: npt.DTypeLike) -> None:
|
||||
self._buf.write(np.ascontiguousarray(arr, dtype=np.dtype(dtype).newbyteorder("<")).tobytes())
|
||||
|
||||
def counted_array(self, arr: npt.ArrayLike, dtype: npt.DTypeLike) -> None:
|
||||
a = np.asarray(arr)
|
||||
self.i32(len(a))
|
||||
self.array(a, dtype)
|
||||
@@ -0,0 +1,449 @@
|
||||
"""BPM (binary bead-pool manifest) reader — plan §6.2, spec IN-040.
|
||||
|
||||
Layout (little-endian, .NET strings)::
|
||||
|
||||
char[3] "BPM"; uint8 1
|
||||
int32 version (flag 0x1000 is XOR-ed off; 3..5 supported)
|
||||
string manifest_name
|
||||
string control_config (version > 1)
|
||||
int32 num_loci
|
||||
int32[num_loci] indexes
|
||||
string[num_loci] names
|
||||
uint8[num_loci] norm_ids (< 100)
|
||||
num_loci × LocusEntry
|
||||
[int32 m; m × string] trailer header lines (optional)
|
||||
|
||||
LocusEntry field order follows gtc2vcf ``locusentry_read`` (MIT) and BeadArrayFiles ``LocusEntry`` (BSD-2).
|
||||
Loci are returned in **names order**, which is the order of every per-locus array in a GTC file.
|
||||
|
||||
Normalization-bin rule (plan §6.2)::
|
||||
|
||||
norm_id' = (norm_id + 100 * assay_type) mod 256 # mimics the vendor's byte wrap (Omni5)
|
||||
lookup = rank of norm_id' among the sorted distinct values
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import struct
|
||||
from dataclasses import dataclass, field
|
||||
from pathlib import Path
|
||||
from typing import Any
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
from simurg.core.errors import FormatError
|
||||
from simurg.formats.binutil import BinaryReader, load_bytes
|
||||
|
||||
ASSAY_INF_II = 0
|
||||
ASSAY_INF_I_AT = 1 # red channel
|
||||
ASSAY_INF_I_GC = 2 # green channel
|
||||
|
||||
_I32 = struct.Struct("<i")
|
||||
_F32X4 = struct.Struct("<4f")
|
||||
|
||||
|
||||
@dataclass
|
||||
class Manifest:
|
||||
"""Columnar manifest; every list/array has one entry per locus in names order."""
|
||||
|
||||
path: str
|
||||
format: str # "bpm" or "csv"
|
||||
version: int
|
||||
manifest_name: str
|
||||
control_config: str
|
||||
names: list[str]
|
||||
ilmn_ids: list[str]
|
||||
snps: list[str]
|
||||
chroms: list[str]
|
||||
map_infos: npt.NDArray[np.int64]
|
||||
ilmn_strands: list[str]
|
||||
source_strands: list[str]
|
||||
ref_strands: list[str]
|
||||
ploidies: list[str]
|
||||
species: list[str]
|
||||
genome_builds: list[str]
|
||||
sources: list[str]
|
||||
source_versions: list[str]
|
||||
customer_strands: list[str]
|
||||
source_seqs: list[str]
|
||||
top_genomic_seqs: list[str]
|
||||
allele_a_probe_seqs: list[str]
|
||||
allele_b_probe_seqs: list[str]
|
||||
address_a: npt.NDArray[np.int32]
|
||||
address_b: npt.NDArray[np.int32]
|
||||
assay_types: npt.NDArray[np.uint8]
|
||||
exp_clusters: npt.NDArray[np.int8]
|
||||
intensity_only: npt.NDArray[np.bool_]
|
||||
norm_ids: npt.NDArray[np.uint8] # raw per-locus IDs (< 100)
|
||||
entry_versions: npt.NDArray[np.int8]
|
||||
base_fractions: npt.NDArray[np.float32] | None = None # (n, 4): A, C, T, G (entry v >= 7)
|
||||
header_lines: list[str] = field(default_factory=list)
|
||||
controls: list[dict[str, str]] = field(default_factory=list)
|
||||
heading: dict[str, str] = field(default_factory=dict)
|
||||
|
||||
@property
|
||||
def num_loci(self) -> int:
|
||||
return len(self.names)
|
||||
|
||||
@property
|
||||
def has_ref_strand(self) -> bool:
|
||||
return any(s for s in self.ref_strands)
|
||||
|
||||
def norm_bin_ids(self) -> npt.NDArray[np.uint8]:
|
||||
"""``(norm_id + 100 * assay_type) mod 256`` per locus."""
|
||||
return ((self.norm_ids.astype(np.int64) + 100 * self.assay_types.astype(np.int64)) % 256).astype(np.uint8)
|
||||
|
||||
def normalization_lookups(self) -> npt.NDArray[np.int32]:
|
||||
"""Index into the GTC's normalization-transform list for every locus."""
|
||||
ids = self.norm_bin_ids()
|
||||
uniq = np.unique(ids) # sorted
|
||||
return np.searchsorted(uniq, ids).astype(np.int32)
|
||||
|
||||
def name_index(self) -> dict[str, int]:
|
||||
return {n: i for i, n in enumerate(self.names)}
|
||||
|
||||
def summary(self) -> dict[str, Any]:
|
||||
at = self.assay_types
|
||||
return {
|
||||
"path": self.path,
|
||||
"format": self.format,
|
||||
"version": self.version,
|
||||
"manifest_name": self.manifest_name,
|
||||
"num_loci": self.num_loci,
|
||||
"entry_versions": sorted({int(v) for v in np.unique(self.entry_versions)}),
|
||||
"infinium_ii": int((at == ASSAY_INF_II).sum()),
|
||||
"infinium_i_at": int((at == ASSAY_INF_I_AT).sum()),
|
||||
"infinium_i_gc": int((at == ASSAY_INF_I_GC).sum()),
|
||||
"intensity_only": int(self.intensity_only.sum()),
|
||||
"normalization_bins": len(np.unique(self.norm_bin_ids())),
|
||||
"has_ref_strand": self.has_ref_strand,
|
||||
"n_controls": len(self.controls),
|
||||
"genome_builds": sorted(set(self.genome_builds) - {""}),
|
||||
"chromosomes": sorted(set(self.chroms), key=_chrom_key),
|
||||
}
|
||||
|
||||
|
||||
def _chrom_key(c: str) -> tuple[int, str]:
|
||||
return (int(c), "") if c.isdigit() else (1000, c)
|
||||
|
||||
|
||||
def parse_control_config(text: str) -> list[dict[str, str]]:
|
||||
"""Control-config lines ``addr[:addr..],Category,Color,Name`` → dicts."""
|
||||
out: list[dict[str, str]] = []
|
||||
for line in text.splitlines():
|
||||
parts = [p.strip() for p in line.split(",")]
|
||||
if len(parts) < 4 or not parts[0]:
|
||||
continue
|
||||
out.append({"address": parts[0], "category": parts[1], "color": parts[2], "name": ",".join(parts[3:])})
|
||||
return out
|
||||
|
||||
|
||||
def read_bpm(path: str | Path) -> Manifest:
|
||||
"""Parse a BPM file. Raises FormatError(E_BPM_FORMAT)."""
|
||||
spath = str(path)
|
||||
buf = load_bytes(path)
|
||||
r = BinaryReader(buf, "E_BPM_FORMAT", spath)
|
||||
if r.size < 8 or r.raw(3, "magic") != b"BPM":
|
||||
raise r.fail("bad magic (not a BPM file)")
|
||||
if (fmt := r.u8("format byte")) != 1:
|
||||
raise r.fail(f"unknown BPM format byte {fmt}")
|
||||
version = r.i32("version")
|
||||
if version & 0x1000:
|
||||
version ^= 0x1000
|
||||
if not 3 <= version <= 5:
|
||||
raise r.fail(f"unsupported BPM version {version}")
|
||||
manifest_name = r.string("manifest name")
|
||||
control_config = r.string("control config") if version > 1 else ""
|
||||
n = r.i32("number of loci")
|
||||
if n <= 0 or n > 50_000_000:
|
||||
raise r.fail(f"implausible number of loci {n}")
|
||||
r.array("<i4", n, "indexes")
|
||||
names = [r.string("locus name") for _ in range(n)]
|
||||
norm_ids = np.array(r.array("u1", n, "normalization ids"), dtype=np.uint8)
|
||||
if (norm_ids >= 100).any():
|
||||
raise r.fail("invalid normalization ID (>= 100)")
|
||||
name_to_idx: dict[str, int] = {}
|
||||
for i, nm in enumerate(names):
|
||||
if nm in name_to_idx:
|
||||
raise r.fail(f"locus name {nm!r} present more than once")
|
||||
name_to_idx[nm] = i
|
||||
|
||||
cols = _parse_locus_entries(r, n, name_to_idx)
|
||||
|
||||
header_lines: list[str] = []
|
||||
if r.remaining() >= 4:
|
||||
m = r.i32("trailer count")
|
||||
if 0 <= m <= 1_000_000:
|
||||
header_lines = [r.string("trailer line") for _ in range(m)]
|
||||
|
||||
return Manifest(
|
||||
path=spath,
|
||||
format="bpm",
|
||||
version=version,
|
||||
manifest_name=manifest_name,
|
||||
control_config=control_config,
|
||||
names=names,
|
||||
norm_ids=norm_ids,
|
||||
header_lines=header_lines,
|
||||
controls=parse_control_config(control_config),
|
||||
**cols,
|
||||
)
|
||||
|
||||
|
||||
def _parse_locus_entries(r: BinaryReader, n: int, name_to_idx: dict[str, int]) -> dict[str, Any]:
|
||||
"""Hot loop: parse ``n`` LocusEntry records (optimised pure Python; Rust in phase P9)."""
|
||||
view = memoryview(r._buf)
|
||||
size = r.size
|
||||
pos = r.pos
|
||||
intern: dict[bytes, str] = {}
|
||||
unpack_i32 = _I32.unpack_from
|
||||
unpack_4f = _F32X4.unpack_from
|
||||
|
||||
def fail(msg: str) -> FormatError:
|
||||
r.pos = pos
|
||||
return r.fail(msg)
|
||||
|
||||
def s_at(p: int) -> tuple[str, int]:
|
||||
# .NET 7-bit length-prefixed UTF-8 string, interned
|
||||
if p >= size:
|
||||
raise fail("unexpected end of file in locus entry")
|
||||
b = view[p]
|
||||
p += 1
|
||||
ln = b & 0x7F
|
||||
shift = 7
|
||||
while b & 0x80:
|
||||
if p >= size or shift > 28:
|
||||
raise fail("bad string length in locus entry")
|
||||
b = view[p]
|
||||
p += 1
|
||||
ln |= (b & 0x7F) << shift
|
||||
shift += 7
|
||||
end = p + ln
|
||||
if end > size:
|
||||
raise fail("string runs past end of file in locus entry")
|
||||
raw = bytes(view[p:end])
|
||||
s = intern.get(raw)
|
||||
if s is None:
|
||||
s = raw.decode("utf-8", errors="replace")
|
||||
if ln <= 64:
|
||||
intern[raw] = s
|
||||
return s, end
|
||||
|
||||
str_cols = [
|
||||
"ilmn_ids",
|
||||
"snps",
|
||||
"chroms",
|
||||
"ilmn_strands",
|
||||
"source_strands",
|
||||
"ref_strands",
|
||||
"ploidies",
|
||||
"species",
|
||||
"genome_builds",
|
||||
"sources",
|
||||
"source_versions",
|
||||
"customer_strands",
|
||||
"source_seqs",
|
||||
"top_genomic_seqs",
|
||||
"allele_a_probe_seqs",
|
||||
"allele_b_probe_seqs",
|
||||
]
|
||||
out: dict[str, list[str]] = {c: [""] * n for c in str_cols}
|
||||
map_infos = np.zeros(n, dtype=np.int64)
|
||||
address_a = np.zeros(n, dtype=np.int32)
|
||||
address_b = np.zeros(n, dtype=np.int32)
|
||||
assay = np.zeros(n, dtype=np.uint8)
|
||||
exp_clusters = np.zeros(n, dtype=np.int8)
|
||||
intensity_only = np.zeros(n, dtype=np.bool_)
|
||||
versions = np.zeros(n, dtype=np.int8)
|
||||
fracs = np.zeros((n, 4), dtype=np.float32)
|
||||
seen = np.zeros(n, dtype=np.bool_)
|
||||
any_v7 = False
|
||||
|
||||
o_ilmn, o_snp, o_chr = out["ilmn_ids"], out["snps"], out["chroms"]
|
||||
o_istr, o_sstr, o_rstr = out["ilmn_strands"], out["source_strands"], out["ref_strands"]
|
||||
o_ploidy, o_species, o_build = out["ploidies"], out["species"], out["genome_builds"]
|
||||
o_source, o_sver, o_cstr = out["sources"], out["source_versions"], out["customer_strands"]
|
||||
o_sseq, o_top, o_pa, o_pb = (
|
||||
out["source_seqs"],
|
||||
out["top_genomic_seqs"],
|
||||
out["allele_a_probe_seqs"],
|
||||
out["allele_b_probe_seqs"],
|
||||
)
|
||||
|
||||
for _ in range(n):
|
||||
if pos + 4 > size:
|
||||
raise fail("unexpected end of file at locus entry")
|
||||
ver = unpack_i32(view, pos)[0]
|
||||
pos += 4
|
||||
if ver not in (4, 6, 7, 8):
|
||||
raise fail(f"unsupported locus entry version {ver}")
|
||||
ilmn_id, pos = s_at(pos)
|
||||
name, pos = s_at(pos)
|
||||
_, pos = s_at(pos)
|
||||
_, pos = s_at(pos)
|
||||
_, pos = s_at(pos)
|
||||
if pos + 4 > size:
|
||||
raise fail("unexpected end of file at locus index")
|
||||
pos += 4 # index (1-based); BeadArrayFiles maps entries by name, as we do
|
||||
_, pos = s_at(pos)
|
||||
ilmn_strand, pos = s_at(pos)
|
||||
snp, pos = s_at(pos)
|
||||
chrom, pos = s_at(pos)
|
||||
ploidy, pos = s_at(pos)
|
||||
species, pos = s_at(pos)
|
||||
map_info, pos = s_at(pos)
|
||||
top_seq, pos = s_at(pos)
|
||||
cust_strand, pos = s_at(pos)
|
||||
if pos + 8 > size:
|
||||
raise fail("unexpected end of file at addresses")
|
||||
a_addr = unpack_i32(view, pos)[0]
|
||||
b_addr = unpack_i32(view, pos + 4)[0]
|
||||
pos += 8
|
||||
probe_a, pos = s_at(pos)
|
||||
probe_b, pos = s_at(pos)
|
||||
build, pos = s_at(pos)
|
||||
source, pos = s_at(pos)
|
||||
source_ver, pos = s_at(pos)
|
||||
source_strand, pos = s_at(pos)
|
||||
source_seq, pos = s_at(pos)
|
||||
at = -1
|
||||
exp_cl = 0
|
||||
int_only = 0
|
||||
if ver >= 6:
|
||||
if pos + 4 > size:
|
||||
raise fail("unexpected end of file at assay type")
|
||||
exp_cl = int.from_bytes(view[pos + 1 : pos + 2], "little", signed=True)
|
||||
int_only = view[pos + 2]
|
||||
at = view[pos + 3]
|
||||
pos += 4
|
||||
if at > 2:
|
||||
raise fail(f"invalid assay type {at}")
|
||||
if (b_addr == 0) != (at == 0):
|
||||
raise fail("assay type is inconsistent with address B")
|
||||
f4 = None
|
||||
if ver >= 7:
|
||||
if pos + 16 > size:
|
||||
raise fail("unexpected end of file at base fractions")
|
||||
f4 = unpack_4f(view, pos)
|
||||
pos += 16
|
||||
ref_strand = ""
|
||||
if ver >= 8:
|
||||
ref_strand, pos = s_at(pos)
|
||||
|
||||
idx = name_to_idx.get(name)
|
||||
if idx is None:
|
||||
raise fail(f"locus entry {name!r} not in the names table")
|
||||
if seen[idx]:
|
||||
raise fail(f"locus entry {name!r} appears twice")
|
||||
seen[idx] = True
|
||||
if at < 0:
|
||||
at = infer_assay_type(probe_a, probe_b, source_seq)
|
||||
o_ilmn[idx] = ilmn_id
|
||||
o_snp[idx] = snp
|
||||
o_chr[idx] = chrom
|
||||
o_istr[idx] = ilmn_strand
|
||||
o_sstr[idx] = source_strand
|
||||
o_rstr[idx] = ref_strand
|
||||
o_ploidy[idx] = ploidy
|
||||
o_species[idx] = species
|
||||
o_build[idx] = build
|
||||
o_source[idx] = source
|
||||
o_sver[idx] = source_ver
|
||||
o_cstr[idx] = cust_strand
|
||||
o_sseq[idx] = source_seq
|
||||
o_top[idx] = top_seq
|
||||
o_pa[idx] = probe_a
|
||||
o_pb[idx] = probe_b
|
||||
try:
|
||||
map_infos[idx] = int(map_info) if map_info else 0
|
||||
except ValueError:
|
||||
raise fail(f"non-integer MapInfo {map_info!r} for {name!r}") from None
|
||||
address_a[idx] = a_addr
|
||||
address_b[idx] = b_addr
|
||||
assay[idx] = at
|
||||
exp_clusters[idx] = exp_cl
|
||||
intensity_only[idx] = bool(int_only)
|
||||
versions[idx] = ver
|
||||
if f4 is not None:
|
||||
fracs[idx] = f4
|
||||
any_v7 = True
|
||||
|
||||
if not seen.all():
|
||||
raise fail(f"{int((~seen).sum())} loci in the names table have no locus entry")
|
||||
r.pos = pos
|
||||
return {
|
||||
**out,
|
||||
"map_infos": map_infos,
|
||||
"address_a": address_a,
|
||||
"address_b": address_b,
|
||||
"assay_types": assay,
|
||||
"exp_clusters": exp_clusters,
|
||||
"intensity_only": intensity_only,
|
||||
"entry_versions": versions,
|
||||
"base_fractions": fracs if any_v7 else None,
|
||||
}
|
||||
|
||||
|
||||
_IUPAC = {
|
||||
"A": 1,
|
||||
"C": 2,
|
||||
"G": 4,
|
||||
"T": 8,
|
||||
"R": 5,
|
||||
"Y": 10,
|
||||
"S": 6,
|
||||
"W": 9,
|
||||
"K": 12,
|
||||
"M": 3,
|
||||
"B": 14,
|
||||
"D": 13,
|
||||
"H": 11,
|
||||
"V": 7,
|
||||
"N": 15,
|
||||
}
|
||||
|
||||
|
||||
def infer_assay_type(probe_a: str, probe_b: str, source_seq: str) -> int:
|
||||
"""Assay type for locus entries older than v6 (gtc2vcf ``get_assay_type``, MIT).
|
||||
|
||||
Returns 0 for Infinium II, 1 for Infinium I red (A/T), 2 for Infinium I green (C/G), 255 if unknown.
|
||||
"""
|
||||
if not probe_a or not source_seq:
|
||||
return 255
|
||||
if not probe_b:
|
||||
return 0
|
||||
left = source_seq.find("[")
|
||||
right = source_seq.find("]")
|
||||
if left <= 0 or right < 0 or right + 1 >= len(source_seq):
|
||||
return 255
|
||||
tl = source_seq[left - 1].upper()
|
||||
tr = source_seq[right + 1].upper()
|
||||
if tl in "AT" and tr in "AT":
|
||||
return 1
|
||||
if tl in "CG" and tr in "CG":
|
||||
return 2
|
||||
i = 2
|
||||
while i <= min(len(probe_a), len(probe_b)) and not (
|
||||
_IUPAC.get(probe_a[-i].upper(), 0) & _IUPAC.get(probe_b[-i].upper(), 0)
|
||||
):
|
||||
i += 1
|
||||
if i > min(len(probe_a), len(probe_b)):
|
||||
return 255
|
||||
t = probe_a[-i].upper()
|
||||
if t in "CGS":
|
||||
return 1
|
||||
if t in "ATW":
|
||||
return 2
|
||||
rules = {
|
||||
("Y", "G"): 1,
|
||||
("Y", "T"): 1,
|
||||
("Y", "A"): 2,
|
||||
("K", "C"): 1,
|
||||
("K", "A"): 2,
|
||||
("M", "G"): 1,
|
||||
("M", "T"): 2,
|
||||
("R", "C"): 1,
|
||||
("R", "T"): 2,
|
||||
}
|
||||
return rules.get((t, tr), 255)
|
||||
@@ -0,0 +1,207 @@
|
||||
"""CSV manifest reader — plan §6.3, spec IN-041/IN-042.
|
||||
|
||||
Sections: ``[Heading]`` (key,value), ``[Assay]`` (header row + one row per locus), ``[Controls]``
|
||||
(``addr[:addr..],Category,Color,Name``). The CSV carries ``SourceSeq``/``TopGenomicSeq`` (needed for indels)
|
||||
but **no normalization IDs** — genotyping always needs the BPM.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import csv
|
||||
import io
|
||||
from pathlib import Path
|
||||
from typing import Any
|
||||
|
||||
import numpy as np
|
||||
|
||||
from simurg.core.errors import FormatError
|
||||
from simurg.formats.bpm import Manifest, infer_assay_type
|
||||
|
||||
REQUIRED_COLUMNS = ("IlmnID", "Name", "SNP", "AddressA_ID", "Chr", "MapInfo")
|
||||
_SEQUENCE_COLUMNS = ("SourceSeq", "TopGenomicSeq", "AlleleA_ProbeSeq", "AlleleB_ProbeSeq")
|
||||
|
||||
|
||||
def _open_text(path: Path) -> io.TextIOBase:
|
||||
raw = path.open("rb")
|
||||
head = raw.read(2)
|
||||
raw.seek(0)
|
||||
if head == b"\x1f\x8b":
|
||||
import gzip
|
||||
|
||||
return io.TextIOWrapper(gzip.GzipFile(fileobj=raw), encoding="utf-8-sig", errors="replace", newline="")
|
||||
return io.TextIOWrapper(raw, encoding="utf-8-sig", errors="replace", newline="")
|
||||
|
||||
|
||||
def _int_or_zero(s: str) -> int:
|
||||
s = s.strip()
|
||||
if not s:
|
||||
return 0
|
||||
try:
|
||||
return int(s)
|
||||
except ValueError:
|
||||
return int(float(s))
|
||||
|
||||
|
||||
def read_csv_manifest(path: str | Path, *, load_sequences: bool = True) -> Manifest:
|
||||
"""Parse a CSV manifest into the same columnar :class:`Manifest` as :func:`read_bpm`.
|
||||
|
||||
``load_sequences=False`` skips probe/source sequences to save memory. Raises E_CSV_MANIFEST_FORMAT.
|
||||
"""
|
||||
p = Path(path)
|
||||
spath = str(p)
|
||||
|
||||
def fail(msg: str) -> FormatError:
|
||||
return FormatError("E_CSV_MANIFEST_FORMAT", f"{spath}: {msg}", path=spath)
|
||||
|
||||
heading: dict[str, str] = {}
|
||||
controls: list[dict[str, str]] = []
|
||||
cols: dict[str, list[str]] = {}
|
||||
header: list[str] | None = None
|
||||
section = None
|
||||
intern: dict[str, str] = {}
|
||||
|
||||
def it(s: str) -> str:
|
||||
v = intern.get(s)
|
||||
if v is None:
|
||||
intern[s] = s
|
||||
v = s
|
||||
return v
|
||||
|
||||
with _open_text(p) as fh:
|
||||
for row in csv.reader(fh):
|
||||
if not row or all(not c.strip() for c in row):
|
||||
continue
|
||||
first = row[0].strip()
|
||||
if first.startswith("[") and first.endswith("]"):
|
||||
section = first[1:-1].strip().lower()
|
||||
continue
|
||||
if section == "heading" and len(row) >= 2:
|
||||
heading[first.rstrip()] = row[1].strip()
|
||||
elif section == "heading":
|
||||
continue
|
||||
elif section == "assay":
|
||||
if header is None:
|
||||
header = [c.strip() for c in row]
|
||||
missing = [c for c in REQUIRED_COLUMNS if c not in header]
|
||||
if missing:
|
||||
raise fail(f"[Assay] section lacks required columns: {', '.join(missing)}")
|
||||
keep = [c for c in header if load_sequences or c not in _SEQUENCE_COLUMNS]
|
||||
cols = {c: [] for c in keep}
|
||||
idx = [(c, header.index(c)) for c in keep]
|
||||
continue
|
||||
if len(row) < len(header):
|
||||
row = row + [""] * (len(header) - len(row))
|
||||
for c, i in idx:
|
||||
v = row[i]
|
||||
cols[c].append(
|
||||
v
|
||||
if c in ("IlmnID", "Name", "SourceSeq", "TopGenomicSeq", "AlleleA_ProbeSeq", "AlleleB_ProbeSeq")
|
||||
else it(v)
|
||||
)
|
||||
elif section == "controls" and len(row) >= 4:
|
||||
controls.append(
|
||||
{
|
||||
"address": first,
|
||||
"category": row[1].strip(),
|
||||
"color": row[2].strip(),
|
||||
"name": ",".join(row[3:]).strip(),
|
||||
}
|
||||
)
|
||||
if header is None:
|
||||
raise fail("no [Assay] section")
|
||||
n = len(cols["Name"])
|
||||
if n == 0:
|
||||
raise fail("[Assay] section has no rows")
|
||||
|
||||
def col(name: str) -> list[str]:
|
||||
return cols.get(name) or [""] * n
|
||||
|
||||
try:
|
||||
addr_a = np.array([_int_or_zero(v) for v in cols["AddressA_ID"]], dtype=np.int32)
|
||||
addr_b = np.array([_int_or_zero(v) for v in col("AddressB_ID")], dtype=np.int32)
|
||||
map_info = np.array([_int_or_zero(v) for v in cols["MapInfo"]], dtype=np.int64)
|
||||
except ValueError as e:
|
||||
raise fail(f"non-numeric address or MapInfo value ({e})") from None
|
||||
probe_a, probe_b, src = col("AlleleA_ProbeSeq"), col("AlleleB_ProbeSeq"), col("SourceSeq")
|
||||
assay = np.zeros(n, dtype=np.uint8)
|
||||
for k in np.flatnonzero(addr_b).tolist():
|
||||
assay[k] = infer_assay_type(probe_a[k], probe_b[k], src[k]) if load_sequences else 255
|
||||
exp = np.array([_int_or_zero(v) for v in col("Exp_Clusters")], dtype=np.int8)
|
||||
int_only = np.array([_int_or_zero(v) != 0 for v in col("Intensity_Only")], dtype=np.bool_)
|
||||
|
||||
return Manifest(
|
||||
path=spath,
|
||||
format="csv",
|
||||
version=0,
|
||||
manifest_name=heading.get("Descriptor File Name", p.name),
|
||||
control_config="",
|
||||
names=cols["Name"],
|
||||
ilmn_ids=cols["IlmnID"],
|
||||
snps=cols["SNP"],
|
||||
chroms=cols["Chr"],
|
||||
map_infos=map_info,
|
||||
ilmn_strands=col("IlmnStrand"),
|
||||
source_strands=col("SourceStrand"),
|
||||
ref_strands=col("RefStrand"),
|
||||
ploidies=col("Ploidy"),
|
||||
species=col("Species"),
|
||||
genome_builds=col("GenomeBuild"),
|
||||
sources=col("Source"),
|
||||
source_versions=col("SourceVersion"),
|
||||
customer_strands=[""] * n,
|
||||
source_seqs=src,
|
||||
top_genomic_seqs=col("TopGenomicSeq"),
|
||||
allele_a_probe_seqs=probe_a,
|
||||
allele_b_probe_seqs=probe_b,
|
||||
address_a=addr_a,
|
||||
address_b=addr_b,
|
||||
assay_types=assay,
|
||||
exp_clusters=exp,
|
||||
intensity_only=int_only,
|
||||
norm_ids=np.zeros(n, dtype=np.uint8),
|
||||
entry_versions=np.zeros(n, dtype=np.int8),
|
||||
controls=controls,
|
||||
heading=heading,
|
||||
)
|
||||
|
||||
|
||||
def assay_format(manifest: Manifest) -> str:
|
||||
"""Platform code from the ``Assay Format`` heading (e.g. ``Infinium EX`` → ``EX``; ``Infinium LCG`` → ``LCG``)."""
|
||||
v = manifest.heading.get("Assay Format", "")
|
||||
return v.split()[-1] if v else ""
|
||||
|
||||
|
||||
def compare_manifests(bpm: Manifest, csvm: Manifest, max_examples: int = 10) -> list[str]:
|
||||
"""BPM ↔ CSV consistency (spec IN-042). Returns human-readable differences (empty list = consistent)."""
|
||||
problems: list[str] = []
|
||||
if bpm.num_loci != csvm.num_loci:
|
||||
problems.append(f"locus count differs: BPM {bpm.num_loci}, CSV {csvm.num_loci}")
|
||||
bidx = bpm.name_index()
|
||||
cidx = csvm.name_index()
|
||||
only_b = [n for n in bpm.names if n not in cidx]
|
||||
only_c = [n for n in csvm.names if n not in bidx]
|
||||
if only_b:
|
||||
problems.append(f"{len(only_b)} loci only in BPM, e.g. {', '.join(only_b[:max_examples])}")
|
||||
if only_c:
|
||||
problems.append(f"{len(only_c)} loci only in CSV, e.g. {', '.join(only_c[:max_examples])}")
|
||||
fields: list[tuple[str, Any, Any]] = [
|
||||
("Chr", bpm.chroms, csvm.chroms),
|
||||
("MapInfo", bpm.map_infos, csvm.map_infos),
|
||||
("SNP", bpm.snps, csvm.snps),
|
||||
("AddressA", bpm.address_a, csvm.address_a),
|
||||
("AddressB", bpm.address_b, csvm.address_b),
|
||||
]
|
||||
for label, bcol, ccol in fields:
|
||||
bad: list[str] = []
|
||||
count = 0
|
||||
for name, i in bidx.items():
|
||||
j = cidx.get(name)
|
||||
if j is None:
|
||||
continue
|
||||
if bcol[i] != ccol[j]:
|
||||
count += 1
|
||||
if len(bad) < max_examples:
|
||||
bad.append(f"{name} ({bcol[i]} vs {ccol[j]})")
|
||||
if count:
|
||||
problems.append(f"{label} differs for {count} loci, e.g. {'; '.join(bad)}")
|
||||
return problems
|
||||
@@ -0,0 +1,42 @@
|
||||
"""File-type sniffing for ``simurg inspect`` (spec CLI-020)."""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import gzip
|
||||
from pathlib import Path
|
||||
|
||||
KINDS = ("idat", "bpm", "egt", "gtc", "csv_manifest", "samplesheet", "vcf", "unknown")
|
||||
|
||||
|
||||
def _head(path: Path, n: int = 4096) -> bytes:
|
||||
with path.open("rb") as fh:
|
||||
b = fh.read(n)
|
||||
if b[:2] == b"\x1f\x8b":
|
||||
try:
|
||||
with gzip.open(path, "rb") as gz:
|
||||
return gz.read(n)
|
||||
except OSError:
|
||||
return b
|
||||
return b
|
||||
|
||||
|
||||
def detect_kind(path: str | Path) -> str:
|
||||
p = Path(path)
|
||||
head = _head(p)
|
||||
name = p.name.lower()
|
||||
if head[:4] == b"IDAT":
|
||||
return "idat"
|
||||
if head[:3] == b"BPM":
|
||||
return "bpm"
|
||||
if head[:3] == b"gtc":
|
||||
return "gtc"
|
||||
if head.startswith(b"##fileformat=VCF"):
|
||||
return "vcf"
|
||||
if name.endswith(".egt") or (len(head) >= 5 and head[:4] == b"\x03\x00\x00\x00"):
|
||||
return "egt"
|
||||
text = head.decode("utf-8", errors="replace").lower()
|
||||
if "[assay]" in text or ("[heading]" in text and "descriptor file name" in text):
|
||||
return "csv_manifest"
|
||||
if "sentrixbarcode_a" in text or "beadchipname" in text or "[data]" in text:
|
||||
return "samplesheet"
|
||||
return "unknown"
|
||||
@@ -0,0 +1,196 @@
|
||||
"""IDAT / GTC discovery and sample-set resolution — spec IN-016, IN-020 … IN-023, IN-030, CLI-006.
|
||||
|
||||
Directory listing order never influences results: every path list is sorted bytewise (``os.fsencode``).
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import os
|
||||
from collections import defaultdict
|
||||
from dataclasses import dataclass, field
|
||||
from pathlib import Path
|
||||
|
||||
from simurg.core.errors import SimurgError
|
||||
from simurg.formats.idat import parse_idat_filename, parse_run_time, read_idat
|
||||
from simurg.formats.samplesheet import SampleSheet
|
||||
|
||||
|
||||
@dataclass
|
||||
class IdatPair:
|
||||
barcode: str
|
||||
position: str
|
||||
grn: Path | None
|
||||
red: Path | None
|
||||
|
||||
@property
|
||||
def sample_id(self) -> str:
|
||||
return f"{self.barcode}_{self.position}"
|
||||
|
||||
@property
|
||||
def complete(self) -> bool:
|
||||
return self.grn is not None and self.red is not None
|
||||
|
||||
@property
|
||||
def folder(self) -> Path:
|
||||
p = self.grn or self.red
|
||||
assert p is not None
|
||||
return p.parent
|
||||
|
||||
|
||||
@dataclass
|
||||
class Discovery:
|
||||
root: Path
|
||||
pairs: dict[str, IdatPair] # sample_id -> pair (complete or not)
|
||||
warnings: list[tuple[str, str]] = field(default_factory=list)
|
||||
|
||||
|
||||
def _bytewise(paths: list[Path]) -> list[Path]:
|
||||
return sorted(paths, key=lambda p: os.fsencode(str(p)))
|
||||
|
||||
|
||||
def walk_files(root: Path, suffixes: tuple[str, ...]) -> list[Path]:
|
||||
"""Recursive listing following symlinks with cycle detection; sorted bytewise."""
|
||||
found: list[Path] = []
|
||||
seen_dirs: set[tuple[int, int]] = set()
|
||||
stack = [root]
|
||||
while stack:
|
||||
d = stack.pop()
|
||||
try:
|
||||
st = d.stat()
|
||||
except OSError:
|
||||
continue
|
||||
key = (st.st_dev, st.st_ino)
|
||||
if key in seen_dirs:
|
||||
continue
|
||||
seen_dirs.add(key)
|
||||
try:
|
||||
entries = list(os.scandir(d))
|
||||
except OSError:
|
||||
continue
|
||||
for e in entries:
|
||||
p = Path(e.path)
|
||||
try:
|
||||
if e.is_dir(follow_symlinks=True):
|
||||
stack.append(p)
|
||||
elif e.is_file(follow_symlinks=True) and e.name.lower().endswith(suffixes):
|
||||
found.append(p)
|
||||
except OSError:
|
||||
continue
|
||||
return _bytewise(found)
|
||||
|
||||
|
||||
def _scan_time_key(p: Path) -> str:
|
||||
try:
|
||||
return parse_run_time(read_idat(p, load_arrays=False).scan_time_raw)
|
||||
except SimurgError:
|
||||
return ""
|
||||
|
||||
|
||||
def discover_idats(root: str | Path, *, on_duplicate: str = "error") -> Discovery:
|
||||
"""Find ``<barcode>_<pos>_{Grn,Red}.idat[.gz]`` under ``root`` and pair them (IN-020 … IN-022)."""
|
||||
rootp = Path(root)
|
||||
if not rootp.is_dir():
|
||||
raise SimurgError("E_FILE_NOT_FOUND", f"IDAT folder not found: {rootp}", path=str(rootp))
|
||||
by_key: dict[tuple[str, str, str], list[Path]] = defaultdict(list)
|
||||
for p in walk_files(rootp, (".idat", ".idat.gz")):
|
||||
parsed = parse_idat_filename(p)
|
||||
if parsed is None:
|
||||
continue
|
||||
by_key[parsed].append(p)
|
||||
|
||||
warnings: list[tuple[str, str]] = []
|
||||
# group channel files by folder so that "the same sample in two folders" is a duplicate pair
|
||||
per_sample: dict[tuple[str, str], dict[Path, dict[str, Path]]] = defaultdict(lambda: defaultdict(dict))
|
||||
for (barcode, pos, ch), paths in by_key.items():
|
||||
for p in paths:
|
||||
per_sample[(barcode, pos)][p.parent][ch] = p
|
||||
|
||||
pairs: dict[str, IdatPair] = {}
|
||||
duplicates: dict[str, list[Path]] = {}
|
||||
for (barcode, pos), folders in sorted(per_sample.items()):
|
||||
candidates = [
|
||||
IdatPair(barcode, pos, chans.get("Grn"), chans.get("Red"))
|
||||
for _, chans in sorted(folders.items(), key=lambda kv: os.fsencode(str(kv[0])))
|
||||
]
|
||||
sid = f"{barcode}_{pos}"
|
||||
if len(candidates) > 1:
|
||||
if on_duplicate == "error":
|
||||
duplicates[sid] = [p for c in candidates for p in (c.grn, c.red) if p is not None]
|
||||
continue
|
||||
if on_duplicate == "newest":
|
||||
# candidates are sorted bytewise; max() keeps the first (smallest path) among equal scan times
|
||||
chosen = max(candidates, key=lambda c: _scan_time_key(c.grn or c.red)) # type: ignore[arg-type]
|
||||
else: # "first": bytewise-smallest path
|
||||
chosen = candidates[0]
|
||||
warnings.append(
|
||||
(
|
||||
"W_DUP_IDAT_RESOLVED",
|
||||
f"{sid}: {len(candidates)} IDAT pairs found; using {chosen.folder} ({on_duplicate})",
|
||||
)
|
||||
)
|
||||
else:
|
||||
chosen = candidates[0]
|
||||
if not chosen.complete:
|
||||
have = "Grn" if chosen.grn else "Red"
|
||||
warnings.append(("W_IDAT_UNPAIRED", f"{sid}: only the {have} IDAT was found in {chosen.folder}"))
|
||||
pairs[sid] = chosen
|
||||
if duplicates:
|
||||
lines = [f"{sid}: " + ", ".join(str(p) for p in ps) for sid, ps in sorted(duplicates.items())]
|
||||
raise SimurgError(
|
||||
"E_DUP_IDAT",
|
||||
f"{len(duplicates)} sample(s) have more than one IDAT pair:\n " + "\n ".join(lines),
|
||||
hint="remove the extra copies or pass --on-duplicate-idat first|newest",
|
||||
)
|
||||
return Discovery(rootp, pairs, warnings)
|
||||
|
||||
|
||||
@dataclass
|
||||
class ResolvedSample:
|
||||
sample_id: str
|
||||
sample_name: str # provisional: sheet Sample_Name or sample_id; IDAT field 406 is applied on read (IN-010)
|
||||
pair: IdatPair | None
|
||||
sheet_row: object | None
|
||||
status: str # "READY" or "MISSING_INPUT"
|
||||
|
||||
|
||||
def resolve_input_folder(
|
||||
cli_folder: str | Path | None, sheet: SampleSheet | None
|
||||
) -> tuple[Path, list[tuple[str, str]]]:
|
||||
"""Apply the input-set rule (CLI-006) and RootFolder precedence (IN-008)."""
|
||||
notes: list[tuple[str, str]] = []
|
||||
if cli_folder is None and sheet is None:
|
||||
raise SimurgError("E_NO_INPUT", "no input: give an IDAT folder, a sample sheet, or both")
|
||||
if cli_folder is not None:
|
||||
folder = Path(cli_folder)
|
||||
if sheet is not None and sheet.root_folder is not None and sheet.root_folder.resolve() != folder.resolve():
|
||||
notes.append(
|
||||
("W_SHEET_ROOTFOLDER_OVERRIDDEN", f"CLI folder {folder} overrides sheet RootFolder {sheet.root_folder}")
|
||||
)
|
||||
return folder, notes
|
||||
assert sheet is not None
|
||||
if sheet.root_folder is None:
|
||||
raise SimurgError("E_NO_INPUT_FOLDER", f"{sheet.path}: sample sheet has no RootFolder and no folder was given")
|
||||
return sheet.root_folder, notes
|
||||
|
||||
|
||||
def resolve_samples(disc: Discovery, sheet: SampleSheet | None) -> tuple[list[ResolvedSample], list[tuple[str, str]]]:
|
||||
"""Match discovered IDAT pairs with sheet rows (IN-010, IN-016). Returns samples in a deterministic order."""
|
||||
notes: list[tuple[str, str]] = []
|
||||
out: list[ResolvedSample] = []
|
||||
if sheet is None:
|
||||
for sid in sorted(disc.pairs):
|
||||
found = disc.pairs[sid]
|
||||
status = "READY" if found.complete else "MISSING_INPUT"
|
||||
out.append(ResolvedSample(sid, sid, found, None, status))
|
||||
return out, notes
|
||||
listed = set()
|
||||
for row in sheet.rows:
|
||||
listed.add(row.sample_id)
|
||||
pair = disc.pairs.get(row.sample_id)
|
||||
status = "READY" if pair is not None and pair.complete else "MISSING_INPUT"
|
||||
name = row.sample_name or row.sample_id
|
||||
out.append(ResolvedSample(row.sample_id, name, pair, row, status))
|
||||
extra = sorted(set(disc.pairs) - listed)
|
||||
if extra:
|
||||
notes.append(("I_UNLISTED_IDATS", f"{len(extra)} IDAT pair(s) not listed in the sample sheet are ignored"))
|
||||
return out, notes
|
||||
@@ -0,0 +1,235 @@
|
||||
"""EGT cluster-file reader — plan §6.4, spec IN-044.
|
||||
|
||||
Layout (little-endian, .NET strings)::
|
||||
|
||||
int32 version (=3); string gencall_version, cluster_version, call_version, normalization_version, date_created
|
||||
uint8 is_wgt (=1); string manifest_name; int32 data_block_version; string opa; int32 num_records
|
||||
num_records × ClusterRecord (fixed size: 3 int32 + 12 float32 + extra floats)
|
||||
num_records × ClusterScore (float sep, float total(GenTrain), float original, uint8 edited)
|
||||
num_records × string genotype (unused)
|
||||
num_records × string locus_name (join to the manifest BY NAME)
|
||||
num_records × int32 address
|
||||
num_records × int32[3] cluster counts (must equal the N in each record)
|
||||
[data_block_version 9: num_records × float, meaning unknown]
|
||||
|
||||
The number of extra floats per record differs between open readers (BeadArrayFiles: 2 for dbv<5, 15 for
|
||||
dbv≥6; gtc2vcf: 0 for dbv 5, 15 for dbv≥7), and the intensity threshold is read for dbv 9 only
|
||||
(BeadArrayFiles) or for dbv≥7 (gtc2vcf) — plan §20 U3. We try the candidate record sizes in order and accept
|
||||
the first whose trailing cluster counts match the per-record N (a strong structural check).
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass
|
||||
from pathlib import Path
|
||||
from typing import Any
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
from simurg.core.errors import FormatError
|
||||
from simurg.formats.binutil import BinaryReader, load_bytes
|
||||
|
||||
_BASE_FIELDS = [
|
||||
("n", "<i4", (3,)),
|
||||
("r_dev", "<f4", (3,)),
|
||||
("r_mean", "<f4", (3,)),
|
||||
("theta_dev", "<f4", (3,)),
|
||||
("theta_mean", "<f4", (3,)),
|
||||
]
|
||||
_SCORE_DTYPE = np.dtype([("cluster_sep", "<f4"), ("total_score", "<f4"), ("original_score", "<f4"), ("edited", "u1")])
|
||||
|
||||
|
||||
def _record_dtype(extra: int) -> np.dtype[Any]:
|
||||
fields: list[Any] = list(_BASE_FIELDS)
|
||||
if extra:
|
||||
fields.append(("extra", "<f4", (extra,)))
|
||||
return np.dtype(fields)
|
||||
|
||||
|
||||
def _extra_candidates(dbv: int) -> list[int]:
|
||||
if dbv < 5:
|
||||
return [2, 0]
|
||||
if dbv == 5:
|
||||
return [0, 2, 15]
|
||||
return [15, 0, 2]
|
||||
|
||||
|
||||
@dataclass
|
||||
class ClusterFile:
|
||||
"""Columnar EGT. Arrays with a trailing axis of 3 are ordered (AA, AB, BB)."""
|
||||
|
||||
path: str
|
||||
version: int
|
||||
gencall_version: str
|
||||
cluster_version: str
|
||||
call_version: str
|
||||
normalization_version: str
|
||||
date_created: str
|
||||
manifest_name: str
|
||||
data_block_version: int
|
||||
opa: str
|
||||
names: list[str]
|
||||
addresses: npt.NDArray[np.int32]
|
||||
n: npt.NDArray[np.int32] # (records, 3)
|
||||
theta_mean: npt.NDArray[np.float32]
|
||||
theta_dev: npt.NDArray[np.float32]
|
||||
r_mean: npt.NDArray[np.float32]
|
||||
r_dev: npt.NDArray[np.float32]
|
||||
intensity_threshold: npt.NDArray[np.float32] # NaN where not stored
|
||||
cluster_sep: npt.NDArray[np.float32]
|
||||
gentrain_score: npt.NDArray[np.float32]
|
||||
original_score: npt.NDArray[np.float32]
|
||||
edited: npt.NDArray[np.bool_]
|
||||
extra_floats: npt.NDArray[np.float32] | None
|
||||
trailer_floats: npt.NDArray[np.float32] | None
|
||||
extra_per_record: int
|
||||
|
||||
@property
|
||||
def num_records(self) -> int:
|
||||
return len(self.names)
|
||||
|
||||
def name_index(self) -> dict[str, int]:
|
||||
return {nm: i for i, nm in enumerate(self.names)}
|
||||
|
||||
def mean_r(self) -> npt.NDArray[np.float64]:
|
||||
"""Count-weighted mean R — the expected R of intensity-only loci (plan §6.4)."""
|
||||
w = self.n.astype(np.float64)
|
||||
tot = w.sum(axis=1)
|
||||
with np.errstate(invalid="ignore", divide="ignore"):
|
||||
return np.asarray((w * self.r_mean).sum(axis=1) / tot, dtype=np.float64)
|
||||
|
||||
def summary(self) -> dict[str, Any]:
|
||||
thr = self.intensity_threshold
|
||||
return {
|
||||
"path": self.path,
|
||||
"version": self.version,
|
||||
"data_block_version": self.data_block_version,
|
||||
"extra_floats_per_record": self.extra_per_record,
|
||||
"num_records": self.num_records,
|
||||
"manifest_name": self.manifest_name,
|
||||
"gencall_version": self.gencall_version,
|
||||
"cluster_version": self.cluster_version,
|
||||
"call_version": self.call_version,
|
||||
"normalization_version": self.normalization_version,
|
||||
"date_created": self.date_created,
|
||||
"intensity_threshold_nonzero": int(np.count_nonzero(np.nan_to_num(thr))),
|
||||
"edited_records": int(self.edited.sum()),
|
||||
"zero_gentrain": int((self.gentrain_score == 0).sum()),
|
||||
}
|
||||
|
||||
|
||||
def read_egt(path: str | Path) -> ClusterFile:
|
||||
"""Parse an EGT cluster file. Raises FormatError(E_EGT_LOAD)."""
|
||||
spath = str(path)
|
||||
buf = load_bytes(path)
|
||||
r = BinaryReader(buf, "E_EGT_LOAD", spath)
|
||||
version = r.i32("version")
|
||||
if version != 3:
|
||||
raise r.fail(f"unsupported cluster file version {version}")
|
||||
gencall_version = r.string()
|
||||
cluster_version = r.string()
|
||||
call_version = r.string()
|
||||
normalization_version = r.string()
|
||||
date_created = r.string()
|
||||
if (is_wgt := r.u8("is_wgt")) != 1:
|
||||
raise r.fail(f"only WGT cluster files are supported (is_wgt={is_wgt})")
|
||||
manifest_name = r.string()
|
||||
dbv = r.i32("data block version")
|
||||
if not 1 <= dbv <= 9:
|
||||
raise r.fail(f"unsupported data block version {dbv}")
|
||||
opa = r.string()
|
||||
nrec = r.i32("number of records")
|
||||
if nrec <= 0 or nrec > 50_000_000:
|
||||
raise r.fail(f"implausible number of records {nrec}")
|
||||
records_start = r.pos
|
||||
|
||||
last_error: FormatError | None = None
|
||||
for extra in _extra_candidates(dbv):
|
||||
try:
|
||||
return _parse_body(
|
||||
r,
|
||||
records_start,
|
||||
nrec,
|
||||
extra,
|
||||
dbv,
|
||||
spath,
|
||||
version,
|
||||
gencall_version,
|
||||
cluster_version,
|
||||
call_version,
|
||||
normalization_version,
|
||||
date_created,
|
||||
manifest_name,
|
||||
opa,
|
||||
)
|
||||
except FormatError as e:
|
||||
last_error = e
|
||||
assert last_error is not None
|
||||
raise last_error
|
||||
|
||||
|
||||
def _parse_body(
|
||||
r: BinaryReader,
|
||||
start: int,
|
||||
nrec: int,
|
||||
extra: int,
|
||||
dbv: int,
|
||||
spath: str,
|
||||
version: int,
|
||||
gencall_version: str,
|
||||
cluster_version: str,
|
||||
call_version: str,
|
||||
normalization_version: str,
|
||||
date_created: str,
|
||||
manifest_name: str,
|
||||
opa: str,
|
||||
) -> ClusterFile:
|
||||
r.seek(start)
|
||||
rec = r.array(_record_dtype(extra), nrec, "cluster records")
|
||||
scores = r.array(_SCORE_DTYPE, nrec, "cluster scores")
|
||||
for _ in range(nrec):
|
||||
r.skip_string("genotype")
|
||||
names = [r.string("locus name") for _ in range(nrec)]
|
||||
addresses = np.array(r.array("<i4", nrec, "addresses"), dtype=np.int32)
|
||||
counts = r.array("<i4", nrec * 3, "cluster counts").reshape(nrec, 3)
|
||||
n = np.array(rec["n"], dtype=np.int32)
|
||||
if not np.array_equal(counts, n):
|
||||
raise r.fail(f"cluster counts do not match record N (assuming {extra} extra floats per record)")
|
||||
trailer = None
|
||||
if dbv == 9 and r.remaining() >= 4 * nrec:
|
||||
trailer = np.array(r.array("<f4", nrec, "trailer floats"), dtype=np.float32)
|
||||
if len(set(names)) != nrec:
|
||||
raise r.fail("duplicate locus names")
|
||||
extra_arr = np.array(rec["extra"], dtype=np.float32) if extra else None
|
||||
if extra_arr is not None and dbv >= 7 and extra == 15:
|
||||
thr = extra_arr[:, 0].copy()
|
||||
else:
|
||||
thr = np.full(nrec, np.nan, dtype=np.float32)
|
||||
return ClusterFile(
|
||||
path=spath,
|
||||
version=version,
|
||||
gencall_version=gencall_version,
|
||||
cluster_version=cluster_version,
|
||||
call_version=call_version,
|
||||
normalization_version=normalization_version,
|
||||
date_created=date_created,
|
||||
manifest_name=manifest_name,
|
||||
data_block_version=dbv,
|
||||
opa=opa,
|
||||
names=names,
|
||||
addresses=addresses,
|
||||
n=n,
|
||||
theta_mean=np.array(rec["theta_mean"], dtype=np.float32),
|
||||
theta_dev=np.array(rec["theta_dev"], dtype=np.float32),
|
||||
r_mean=np.array(rec["r_mean"], dtype=np.float32),
|
||||
r_dev=np.array(rec["r_dev"], dtype=np.float32),
|
||||
intensity_threshold=thr,
|
||||
cluster_sep=np.array(scores["cluster_sep"], dtype=np.float32),
|
||||
gentrain_score=np.array(scores["total_score"], dtype=np.float32),
|
||||
original_score=np.array(scores["original_score"], dtype=np.float32),
|
||||
edited=np.array(scores["edited"] != 0),
|
||||
extra_floats=extra_arr,
|
||||
trailer_floats=trailer,
|
||||
extra_per_record=extra,
|
||||
)
|
||||
@@ -0,0 +1,548 @@
|
||||
"""GTC genotype file (versions 3–5) reader and writer — plan §6.5, spec IN-031/IN-032/IN-033, OUT-001.
|
||||
|
||||
Layout::
|
||||
|
||||
char[3] "gtc"; uint8 version; int32 M; M × {uint16 id; uint32 offset}
|
||||
(for ids 1, 2, 3 the "offset" IS the value)
|
||||
|
||||
Authoritative description: ``docs/GTC_File_Format_v5.pdf`` in Illumina BeadArrayFiles (BSD-2); canonical write
|
||||
order follows gtc2vcf ``gtc_write`` (MIT). The reader records the TOC order and the on-disk data order, so a
|
||||
file can be re-written byte-for-byte (P1 acceptance: round-trip).
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import math
|
||||
from dataclasses import dataclass, field
|
||||
from pathlib import Path
|
||||
from typing import Any
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
from simurg.formats.binutil import BinaryReader, BinaryWriter, load_bytes
|
||||
|
||||
NUM_SNPS = 1
|
||||
PLOIDY = 2
|
||||
PLOIDY_TYPE = 3
|
||||
SAMPLE_NAME = 10
|
||||
SAMPLE_PLATE = 11
|
||||
SAMPLE_WELL = 12
|
||||
CLUSTER_FILE = 100
|
||||
SNP_MANIFEST = 101
|
||||
IMAGING_DATE = 200
|
||||
AUTOCALL_DATE = 201
|
||||
AUTOCALL_VERSION = 300
|
||||
NORMALIZATION_TRANSFORMS = 400
|
||||
CONTROLS_X = 500
|
||||
CONTROLS_Y = 501
|
||||
RAW_X = 1000
|
||||
RAW_Y = 1001
|
||||
GENOTYPES = 1002
|
||||
BASE_CALLS = 1003
|
||||
GENOTYPE_SCORES = 1004
|
||||
SCANNER_DATA = 1005
|
||||
CALL_RATE = 1006
|
||||
GENDER = 1007
|
||||
LOGR_DEV = 1008
|
||||
GC10 = 1009
|
||||
DX = 1010
|
||||
SAMPLE_DATA = 1011
|
||||
B_ALLELE_FREQS = 1012
|
||||
LOGR_RATIOS = 1013
|
||||
PERCENTILES_X = 1014
|
||||
PERCENTILES_Y = 1015
|
||||
SLIDE_IDENTIFIER = 1016
|
||||
|
||||
INLINE_IDS = (NUM_SNPS, PLOIDY, PLOIDY_TYPE)
|
||||
|
||||
TOC_NAMES: dict[int, str] = {
|
||||
NUM_SNPS: "num_snps",
|
||||
PLOIDY: "ploidy",
|
||||
PLOIDY_TYPE: "ploidy_type",
|
||||
SAMPLE_NAME: "sample_name",
|
||||
SAMPLE_PLATE: "sample_plate",
|
||||
SAMPLE_WELL: "sample_well",
|
||||
CLUSTER_FILE: "cluster_file",
|
||||
SNP_MANIFEST: "snp_manifest",
|
||||
IMAGING_DATE: "imaging_date",
|
||||
AUTOCALL_DATE: "autocall_date",
|
||||
AUTOCALL_VERSION: "autocall_version",
|
||||
NORMALIZATION_TRANSFORMS: "normalization_transforms",
|
||||
CONTROLS_X: "controls_x",
|
||||
CONTROLS_Y: "controls_y",
|
||||
RAW_X: "raw_x",
|
||||
RAW_Y: "raw_y",
|
||||
GENOTYPES: "genotypes",
|
||||
BASE_CALLS: "base_calls",
|
||||
GENOTYPE_SCORES: "genotype_scores",
|
||||
SCANNER_DATA: "scanner_data",
|
||||
CALL_RATE: "call_rate",
|
||||
GENDER: "gender",
|
||||
LOGR_DEV: "logr_dev",
|
||||
GC10: "gc10",
|
||||
DX: "dx",
|
||||
SAMPLE_DATA: "sample_data",
|
||||
B_ALLELE_FREQS: "b_allele_freqs",
|
||||
LOGR_RATIOS: "logr_ratios",
|
||||
PERCENTILES_X: "percentiles_x",
|
||||
PERCENTILES_Y: "percentiles_y",
|
||||
SLIDE_IDENTIFIER: "slide_identifier",
|
||||
}
|
||||
|
||||
# canonical order (gtc2vcf gtc_write); v3 omits PLOIDY, PLOIDY_TYPE and 1012..1016
|
||||
CANONICAL_ORDER_V5 = [
|
||||
NUM_SNPS,
|
||||
PLOIDY,
|
||||
PLOIDY_TYPE,
|
||||
SAMPLE_NAME,
|
||||
SAMPLE_PLATE,
|
||||
SAMPLE_WELL,
|
||||
CLUSTER_FILE,
|
||||
SNP_MANIFEST,
|
||||
IMAGING_DATE,
|
||||
AUTOCALL_DATE,
|
||||
AUTOCALL_VERSION,
|
||||
NORMALIZATION_TRANSFORMS,
|
||||
CONTROLS_X,
|
||||
CONTROLS_Y,
|
||||
RAW_X,
|
||||
RAW_Y,
|
||||
GENOTYPES,
|
||||
BASE_CALLS,
|
||||
GENOTYPE_SCORES,
|
||||
SCANNER_DATA,
|
||||
CALL_RATE,
|
||||
GENDER,
|
||||
LOGR_DEV,
|
||||
GC10,
|
||||
DX,
|
||||
SAMPLE_DATA,
|
||||
B_ALLELE_FREQS,
|
||||
LOGR_RATIOS,
|
||||
PERCENTILES_X,
|
||||
PERCENTILES_Y,
|
||||
SLIDE_IDENTIFIER,
|
||||
]
|
||||
_V4_PLUS = {PLOIDY, PLOIDY_TYPE, B_ALLELE_FREQS, LOGR_RATIOS, PERCENTILES_X, PERCENTILES_Y, SLIDE_IDENTIFIER}
|
||||
CANONICAL_ORDER_V3 = [i for i in CANONICAL_ORDER_V5 if i not in _V4_PLUS]
|
||||
|
||||
GENOTYPE_CODES = [
|
||||
"NC",
|
||||
"AA",
|
||||
"AB",
|
||||
"BB",
|
||||
"NULL",
|
||||
"A",
|
||||
"B",
|
||||
"AAA",
|
||||
"AAB",
|
||||
"ABB",
|
||||
"BBB",
|
||||
"AAAA",
|
||||
"AAAB",
|
||||
"AABB",
|
||||
"ABBB",
|
||||
"BBBB",
|
||||
"AAAAA",
|
||||
"AAAAB",
|
||||
"AAABB",
|
||||
"AABBB",
|
||||
"ABBBB",
|
||||
"BBBBB",
|
||||
"AAAAAA",
|
||||
"AAAAAB",
|
||||
"AAAABB",
|
||||
"AAABBB",
|
||||
"AABBBB",
|
||||
"ABBBBB",
|
||||
"BBBBBB",
|
||||
"AAAAAAA",
|
||||
"AAAAAAB",
|
||||
"AAAAABB",
|
||||
"AAAABBB",
|
||||
"AAABBBB",
|
||||
"AABBBBB",
|
||||
"ABBBBBB",
|
||||
"BBBBBBB",
|
||||
"AAAAAAAA",
|
||||
"AAAAAAAB",
|
||||
"AAAAAABB",
|
||||
"AAAAABBB",
|
||||
"AAAABBBB",
|
||||
"AAABBBBB",
|
||||
"AABBBBBB",
|
||||
"ABBBBBBB",
|
||||
"BBBBBBBB",
|
||||
]
|
||||
|
||||
XFORM_DTYPE = np.dtype(
|
||||
[
|
||||
("version", "<i4"),
|
||||
("offset_x", "<f4"),
|
||||
("offset_y", "<f4"),
|
||||
("scale_x", "<f4"),
|
||||
("scale_y", "<f4"),
|
||||
("shear", "<f4"),
|
||||
("theta", "<f4"),
|
||||
("reserved", "<f4", (6,)),
|
||||
]
|
||||
)
|
||||
assert XFORM_DTYPE.itemsize == 52
|
||||
SAMPLE_DATA_DTYPE = np.dtype(
|
||||
[("p50gc", "<f4"), ("num_calls", "<i4"), ("num_no_calls", "<i4"), ("num_intensity_only", "<i4")]
|
||||
)
|
||||
|
||||
|
||||
@dataclass
|
||||
class ScannerData:
|
||||
name: str = ""
|
||||
pmt_green: bytes = b"\x00\x00\x00\x00" # 4 raw bytes: int32 in BeadArrayFiles, float in gtc2vcf
|
||||
pmt_red: bytes = b"\x00\x00\x00\x00"
|
||||
version: str = ""
|
||||
user: str = ""
|
||||
|
||||
@property
|
||||
def pmt_green_int(self) -> int:
|
||||
return int.from_bytes(self.pmt_green, "little", signed=True)
|
||||
|
||||
@property
|
||||
def pmt_red_int(self) -> int:
|
||||
return int.from_bytes(self.pmt_red, "little", signed=True)
|
||||
|
||||
|
||||
@dataclass
|
||||
class Gtc:
|
||||
"""In-memory GTC content. Arrays are per locus in BPM names order."""
|
||||
|
||||
version: int = 5
|
||||
num_snps: int = 0
|
||||
ploidy: int = 2
|
||||
ploidy_type: int = 1
|
||||
sample_name: str = ""
|
||||
sample_plate: str = ""
|
||||
sample_well: str = ""
|
||||
cluster_file: str = ""
|
||||
snp_manifest: str = ""
|
||||
imaging_date: str = ""
|
||||
autocall_date: str = ""
|
||||
autocall_version: str = ""
|
||||
normalization_transforms: npt.NDArray[Any] = field(default_factory=lambda: np.zeros(0, XFORM_DTYPE))
|
||||
controls_x: npt.NDArray[np.uint16] = field(default_factory=lambda: np.zeros(0, np.uint16))
|
||||
controls_y: npt.NDArray[np.uint16] = field(default_factory=lambda: np.zeros(0, np.uint16))
|
||||
raw_x: npt.NDArray[np.uint16] = field(default_factory=lambda: np.zeros(0, np.uint16))
|
||||
raw_y: npt.NDArray[np.uint16] = field(default_factory=lambda: np.zeros(0, np.uint16))
|
||||
genotypes: npt.NDArray[np.uint8] = field(default_factory=lambda: np.zeros(0, np.uint8))
|
||||
base_calls: npt.NDArray[np.uint8] = field(default_factory=lambda: np.zeros((0, 2), np.uint8))
|
||||
genotype_scores: npt.NDArray[np.float32] = field(default_factory=lambda: np.zeros(0, np.float32))
|
||||
scanner_data: ScannerData = field(default_factory=ScannerData)
|
||||
call_rate: float = 0.0
|
||||
gender: bytes = b"U"
|
||||
logr_dev: float = 0.0
|
||||
gc10: float = 0.0
|
||||
dx: int = 0
|
||||
sample_data: npt.NDArray[Any] = field(default_factory=lambda: np.zeros(1, SAMPLE_DATA_DTYPE)[0])
|
||||
b_allele_freqs: npt.NDArray[np.float32] = field(default_factory=lambda: np.zeros(0, np.float32))
|
||||
logr_ratios: npt.NDArray[np.float32] = field(default_factory=lambda: np.zeros(0, np.float32))
|
||||
percentiles_x: npt.NDArray[np.uint16] = field(default_factory=lambda: np.zeros(3, np.uint16))
|
||||
percentiles_y: npt.NDArray[np.uint16] = field(default_factory=lambda: np.zeros(3, np.uint16))
|
||||
slide_identifier: str = ""
|
||||
# layout bookkeeping (set by the reader; None → canonical order)
|
||||
toc_order: list[int] | None = None
|
||||
data_order: list[int] | None = None
|
||||
present: set[int] | None = None
|
||||
path: str | None = None
|
||||
|
||||
# -- convenience ------------------------------------------------------------------------------------
|
||||
@property
|
||||
def gc50(self) -> float:
|
||||
return float(self.sample_data["p50gc"])
|
||||
|
||||
@property
|
||||
def num_calls(self) -> int:
|
||||
return int(self.sample_data["num_calls"])
|
||||
|
||||
@property
|
||||
def num_no_calls(self) -> int:
|
||||
return int(self.sample_data["num_no_calls"])
|
||||
|
||||
@property
|
||||
def num_intensity_only(self) -> int:
|
||||
return int(self.sample_data["num_intensity_only"])
|
||||
|
||||
@property
|
||||
def gender_char(self) -> str:
|
||||
return self.gender.decode("latin-1") if self.gender else ""
|
||||
|
||||
def genotype_strings(self) -> list[str]:
|
||||
return [GENOTYPE_CODES[g] if g < len(GENOTYPE_CODES) else "?" for g in self.genotypes]
|
||||
|
||||
def base_call_strings(self) -> list[str]:
|
||||
"""TOP-strand base calls; '-' for no-call (BeadArrayFiles ``get_base_calls`` for diploids)."""
|
||||
return [bytes(bc).decode("latin-1") for bc in self.base_calls]
|
||||
|
||||
def summary(self) -> dict[str, Any]:
|
||||
def f(v: float) -> float | None:
|
||||
return None if math.isnan(v) else round(float(v), 6)
|
||||
|
||||
return {
|
||||
"path": self.path,
|
||||
"version": self.version,
|
||||
"num_snps": self.num_snps,
|
||||
"ploidy": self.ploidy,
|
||||
"ploidy_type": self.ploidy_type,
|
||||
"sample_name": self.sample_name,
|
||||
"sample_plate": self.sample_plate,
|
||||
"sample_well": self.sample_well,
|
||||
"slide_identifier": self.slide_identifier,
|
||||
"cluster_file": self.cluster_file,
|
||||
"snp_manifest": self.snp_manifest,
|
||||
"imaging_date": self.imaging_date,
|
||||
"autocall_date": self.autocall_date,
|
||||
"autocall_version": self.autocall_version,
|
||||
"n_normalization_transforms": len(self.normalization_transforms),
|
||||
"n_controls": len(self.controls_x),
|
||||
"call_rate": f(self.call_rate),
|
||||
"gender": self.gender_char,
|
||||
"logr_dev": f(self.logr_dev),
|
||||
"gc10": f(self.gc10),
|
||||
"gc50": f(self.gc50),
|
||||
"num_calls": self.num_calls,
|
||||
"num_no_calls": self.num_no_calls,
|
||||
"num_intensity_only": self.num_intensity_only,
|
||||
"percentiles_x": [int(v) for v in self.percentiles_x],
|
||||
"percentiles_y": [int(v) for v in self.percentiles_y],
|
||||
"scanner": {
|
||||
"name": self.scanner_data.name,
|
||||
"version": self.scanner_data.version,
|
||||
"user": self.scanner_data.user,
|
||||
},
|
||||
"toc": [TOC_NAMES.get(i, str(i)) for i in (self.toc_order or [])],
|
||||
}
|
||||
|
||||
|
||||
def read_gtc(path: str | Path, *, check_write_complete: bool = True) -> Gtc:
|
||||
"""Parse a GTC file. Raises FormatError(E_GTC_CORRUPT) on any structural problem (IN-032)."""
|
||||
spath = str(path)
|
||||
buf = load_bytes(path)
|
||||
r = BinaryReader(buf, "E_GTC_CORRUPT", spath)
|
||||
if r.size < 8 or r.raw(3, "magic") != b"gtc":
|
||||
raise r.fail("bad magic (not a GTC file)")
|
||||
version = r.u8("version")
|
||||
if version not in (3, 4, 5):
|
||||
raise r.fail(f"unsupported GTC version {version}")
|
||||
m = r.i32("TOC size")
|
||||
if not 0 < m < 1024:
|
||||
raise r.fail(f"implausible TOC size {m}")
|
||||
toc: dict[int, int] = {}
|
||||
order: list[int] = []
|
||||
for _ in range(m):
|
||||
tid = r.u16("TOC id")
|
||||
off = r.u32("TOC offset")
|
||||
if tid in toc:
|
||||
raise r.fail(f"duplicate TOC id {tid}")
|
||||
toc[tid] = off
|
||||
order.append(tid)
|
||||
g = Gtc(version=version, toc_order=order, present=set(toc), path=spath)
|
||||
for tid in INLINE_IDS:
|
||||
if tid in toc:
|
||||
setattr(g, TOC_NAMES[tid], _as_i32(toc[tid]))
|
||||
if NUM_SNPS not in toc:
|
||||
raise r.fail("missing NumSNPs entry")
|
||||
n = g.num_snps
|
||||
if n < 0:
|
||||
raise r.fail(f"negative NumSNPs {n}")
|
||||
|
||||
def at(tid: int) -> bool:
|
||||
if tid not in toc:
|
||||
return False
|
||||
r.seek(toc[tid], TOC_NAMES.get(tid, str(tid)))
|
||||
return True
|
||||
|
||||
def arr(tid: int, dtype: Any, expect: int | None = None) -> npt.NDArray[Any]:
|
||||
a = r.counted_array(dtype, TOC_NAMES[tid])
|
||||
if expect is not None and len(a) != expect:
|
||||
raise r.fail(f"{TOC_NAMES[tid]} has {len(a)} entries, expected {expect}")
|
||||
return np.array(a)
|
||||
|
||||
for tid, attr in (
|
||||
(SAMPLE_NAME, "sample_name"),
|
||||
(SAMPLE_PLATE, "sample_plate"),
|
||||
(SAMPLE_WELL, "sample_well"),
|
||||
(CLUSTER_FILE, "cluster_file"),
|
||||
(SNP_MANIFEST, "snp_manifest"),
|
||||
(IMAGING_DATE, "imaging_date"),
|
||||
(AUTOCALL_DATE, "autocall_date"),
|
||||
(AUTOCALL_VERSION, "autocall_version"),
|
||||
(SLIDE_IDENTIFIER, "slide_identifier"),
|
||||
):
|
||||
if at(tid):
|
||||
setattr(g, attr, r.string(attr))
|
||||
if at(NORMALIZATION_TRANSFORMS):
|
||||
g.normalization_transforms = arr(NORMALIZATION_TRANSFORMS, XFORM_DTYPE)
|
||||
if at(CONTROLS_X):
|
||||
g.controls_x = arr(CONTROLS_X, "<u2")
|
||||
if at(CONTROLS_Y):
|
||||
g.controls_y = arr(CONTROLS_Y, "<u2")
|
||||
if at(RAW_X):
|
||||
g.raw_x = arr(RAW_X, "<u2", n)
|
||||
if at(RAW_Y):
|
||||
g.raw_y = arr(RAW_Y, "<u2", n)
|
||||
if at(GENOTYPES):
|
||||
g.genotypes = arr(GENOTYPES, "u1", n)
|
||||
if at(BASE_CALLS):
|
||||
cnt = r.i32("base calls count")
|
||||
if cnt != n:
|
||||
raise r.fail(f"base_calls has {cnt} entries, expected {n}")
|
||||
g.base_calls = np.array(r.array("u1", 2 * n, "base calls")).reshape(n, 2)
|
||||
if at(GENOTYPE_SCORES):
|
||||
g.genotype_scores = arr(GENOTYPE_SCORES, "<f4", n)
|
||||
if at(SCANNER_DATA):
|
||||
g.scanner_data = ScannerData(
|
||||
name=r.string("scanner name"),
|
||||
pmt_green=r.raw(4, "pmt green"),
|
||||
pmt_red=r.raw(4, "pmt red"),
|
||||
version=r.string("scanner version"),
|
||||
user=r.string("imaging user"),
|
||||
)
|
||||
if at(CALL_RATE):
|
||||
g.call_rate = r.f32("call rate")
|
||||
if at(GENDER):
|
||||
g.gender = r.raw(1, "gender")
|
||||
if at(LOGR_DEV):
|
||||
g.logr_dev = r.f32("logr dev")
|
||||
if at(GC10):
|
||||
g.gc10 = r.f32("gc10")
|
||||
if at(DX):
|
||||
g.dx = r.i32("dx")
|
||||
if at(SAMPLE_DATA):
|
||||
g.sample_data = np.array(r.array(SAMPLE_DATA_DTYPE, 1, "sample data"))[0]
|
||||
if at(B_ALLELE_FREQS):
|
||||
g.b_allele_freqs = arr(B_ALLELE_FREQS, "<f4", n)
|
||||
if at(LOGR_RATIOS):
|
||||
g.logr_ratios = arr(LOGR_RATIOS, "<f4", n)
|
||||
if at(PERCENTILES_X):
|
||||
g.percentiles_x = np.array(r.array("<u2", 3, "percentiles x"))
|
||||
if at(PERCENTILES_Y):
|
||||
g.percentiles_y = np.array(r.array("<u2", 3, "percentiles y"))
|
||||
unknown = [t for t in toc if t not in TOC_NAMES]
|
||||
if unknown:
|
||||
raise r.fail(f"unknown TOC ids {unknown}")
|
||||
g.data_order = sorted((t for t in toc if t not in INLINE_IDS), key=lambda t: toc[t])
|
||||
if check_write_complete:
|
||||
required = {3: SAMPLE_DATA, 4: LOGR_RATIOS, 5: SLIDE_IDENTIFIER}[version]
|
||||
if required not in toc:
|
||||
raise r.fail(f"GTC v{version} is incomplete (missing {TOC_NAMES[required]})")
|
||||
return g
|
||||
|
||||
|
||||
def _as_i32(v: int) -> int:
|
||||
return v - (1 << 32) if v >= 1 << 31 else v
|
||||
|
||||
|
||||
def _write_field(w: BinaryWriter, g: Gtc, tid: int) -> None:
|
||||
n = g.num_snps
|
||||
if tid in (
|
||||
SAMPLE_NAME,
|
||||
SAMPLE_PLATE,
|
||||
SAMPLE_WELL,
|
||||
CLUSTER_FILE,
|
||||
SNP_MANIFEST,
|
||||
IMAGING_DATE,
|
||||
AUTOCALL_DATE,
|
||||
AUTOCALL_VERSION,
|
||||
SLIDE_IDENTIFIER,
|
||||
):
|
||||
w.string(getattr(g, TOC_NAMES[tid]))
|
||||
elif tid == NORMALIZATION_TRANSFORMS:
|
||||
xf = np.asarray(g.normalization_transforms, dtype=XFORM_DTYPE)
|
||||
w.i32(len(xf))
|
||||
w.write(xf.tobytes())
|
||||
elif tid in (CONTROLS_X, CONTROLS_Y):
|
||||
w.counted_array(getattr(g, TOC_NAMES[tid]), "<u2")
|
||||
elif tid in (RAW_X, RAW_Y):
|
||||
_check_len(g, tid, n)
|
||||
w.counted_array(getattr(g, TOC_NAMES[tid]), "<u2")
|
||||
elif tid == GENOTYPES:
|
||||
_check_len(g, tid, n)
|
||||
w.counted_array(g.genotypes, "u1")
|
||||
elif tid == BASE_CALLS:
|
||||
bc = np.asarray(g.base_calls, dtype=np.uint8).reshape(-1, 2)
|
||||
if len(bc) != n:
|
||||
raise ValueError(f"base_calls has {len(bc)} rows, expected {n}")
|
||||
w.i32(n)
|
||||
w.write(bc.tobytes())
|
||||
elif tid in (GENOTYPE_SCORES, B_ALLELE_FREQS, LOGR_RATIOS):
|
||||
_check_len(g, tid, n)
|
||||
w.counted_array(getattr(g, TOC_NAMES[tid]), "<f4")
|
||||
elif tid == SCANNER_DATA:
|
||||
s = g.scanner_data
|
||||
w.string(s.name)
|
||||
w.write(s.pmt_green)
|
||||
w.write(s.pmt_red)
|
||||
w.string(s.version)
|
||||
w.string(s.user)
|
||||
elif tid in (CALL_RATE, LOGR_DEV, GC10):
|
||||
w.f32(getattr(g, TOC_NAMES[tid]))
|
||||
elif tid == GENDER:
|
||||
w.write(g.gender[:1] or b"U")
|
||||
elif tid == DX:
|
||||
w.i32(g.dx)
|
||||
elif tid == SAMPLE_DATA:
|
||||
w.write(np.asarray(g.sample_data, dtype=SAMPLE_DATA_DTYPE).tobytes())
|
||||
elif tid in (PERCENTILES_X, PERCENTILES_Y):
|
||||
w.array(getattr(g, TOC_NAMES[tid]), "<u2")
|
||||
else:
|
||||
raise ValueError(f"cannot write TOC id {tid}")
|
||||
|
||||
|
||||
def _check_len(g: Gtc, tid: int, n: int) -> None:
|
||||
if len(getattr(g, TOC_NAMES[tid])) != n:
|
||||
raise ValueError(f"{TOC_NAMES[tid]} length {len(getattr(g, TOC_NAMES[tid]))} != num_snps {n}")
|
||||
|
||||
|
||||
def encode_gtc(g: Gtc, *, canonical: bool = False) -> bytes:
|
||||
"""Serialize. ``canonical=True`` ignores the recorded layout and uses the gtc2vcf order for ``g.version``."""
|
||||
if canonical or g.toc_order is None:
|
||||
toc_order = list(CANONICAL_ORDER_V5 if g.version >= 4 else CANONICAL_ORDER_V3)
|
||||
data_order = [t for t in toc_order if t not in INLINE_IDS]
|
||||
else:
|
||||
toc_order = list(g.toc_order)
|
||||
data_order = list(g.data_order or [t for t in toc_order if t not in INLINE_IDS])
|
||||
header_len = 4 + 4 + 6 * len(toc_order)
|
||||
body = BinaryWriter()
|
||||
offsets: dict[int, int] = {}
|
||||
for tid in data_order:
|
||||
offsets[tid] = header_len + body.tell()
|
||||
_write_field(body, g, tid)
|
||||
w = BinaryWriter()
|
||||
w.write(b"gtc")
|
||||
w.u8(g.version)
|
||||
w.i32(len(toc_order))
|
||||
for tid in toc_order:
|
||||
w.u16(tid)
|
||||
if tid in INLINE_IDS:
|
||||
w.u32(getattr(g, TOC_NAMES[tid]) & 0xFFFFFFFF)
|
||||
else:
|
||||
w.u32(offsets[tid])
|
||||
w.write(body.getvalue())
|
||||
return w.getvalue()
|
||||
|
||||
|
||||
def write_gtc(g: Gtc, path: str | Path, *, canonical: bool = False) -> None:
|
||||
"""Write atomically (spec CLI-011)."""
|
||||
from simurg.core.atomic import atomic_write_bytes
|
||||
|
||||
atomic_write_bytes(Path(path), encode_gtc(g, canonical=canonical))
|
||||
|
||||
|
||||
def validate_gtc(g: Gtc) -> list[str]:
|
||||
"""Light semantic checks used by ``simurg inspect`` (not raised)."""
|
||||
problems: list[str] = []
|
||||
n = g.num_snps
|
||||
for name in ("raw_x", "raw_y", "genotypes", "genotype_scores"):
|
||||
if len(getattr(g, name)) not in (0, n):
|
||||
problems.append(f"{name} length {len(getattr(g, name))} != num_snps {n}")
|
||||
if len(g.genotypes) and int(g.genotypes.max(initial=0)) >= len(GENOTYPE_CODES):
|
||||
problems.append("genotype code out of range")
|
||||
return problems
|
||||
@@ -0,0 +1,343 @@
|
||||
"""IDAT (binary, version 3) reader — plan §6.1, spec IN-023/IN-024/IN-025.
|
||||
|
||||
Layout (little-endian)::
|
||||
|
||||
char[4] "IDAT"; int64 version (3); int32 n_fields; n_fields × {uint16 code; int64 offset}
|
||||
|
||||
Field codes follow illuminaio (GPL, consulted as documentation only), gtc2vcf ``idat2gtc.c`` (MIT) and
|
||||
BeadArrayFiles (BSD-2). The element type of the MidBlock (code 200) differs between those readers (int32 in
|
||||
illuminaio, uint8 in gtc2vcf); we infer it from the space available before the next section (plan §20 U2).
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import datetime as _dt
|
||||
import re
|
||||
from dataclasses import dataclass, field
|
||||
from pathlib import Path
|
||||
from typing import Any
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
from simurg.core.errors import FormatError
|
||||
from simurg.formats.binutil import BinaryReader, load_bytes
|
||||
|
||||
MAGIC = b"IDAT"
|
||||
|
||||
# field codes
|
||||
N_SNPS = 1000
|
||||
ILLUMINA_ID = 102
|
||||
SD = 103
|
||||
MEAN = 104
|
||||
NBEADS = 107
|
||||
MID_BLOCK = 200
|
||||
RUN_INFO = 300
|
||||
RED_GREEN = 400
|
||||
MANIFEST = 401
|
||||
BARCODE = 402
|
||||
CHIP_TYPE = 403
|
||||
POSITION = 404
|
||||
BEADSET = 405
|
||||
SAMPLE_ID = 406
|
||||
DESCRIPTION = 407
|
||||
PLATE = 408
|
||||
WELL = 409
|
||||
SAMPLE_COUNT = 410
|
||||
VLN = 510
|
||||
|
||||
FIELD_NAMES: dict[int, str] = {
|
||||
N_SNPS: "n_addresses",
|
||||
ILLUMINA_ID: "address_ids",
|
||||
SD: "sd",
|
||||
MEAN: "mean",
|
||||
NBEADS: "n_beads",
|
||||
MID_BLOCK: "mid_block",
|
||||
RUN_INFO: "run_info",
|
||||
RED_GREEN: "red_green",
|
||||
MANIFEST: "manifest",
|
||||
BARCODE: "barcode",
|
||||
CHIP_TYPE: "chip_type",
|
||||
POSITION: "position",
|
||||
BEADSET: "beadset",
|
||||
SAMPLE_ID: "sample_id",
|
||||
DESCRIPTION: "description",
|
||||
PLATE: "plate",
|
||||
WELL: "well",
|
||||
SAMPLE_COUNT: "sample_count",
|
||||
VLN: "vln",
|
||||
}
|
||||
|
||||
_STRING_FIELDS = {
|
||||
MANIFEST: "manifest",
|
||||
BARCODE: "barcode",
|
||||
CHIP_TYPE: "chip_type",
|
||||
POSITION: "position",
|
||||
BEADSET: "beadset",
|
||||
SAMPLE_ID: "sample_id",
|
||||
DESCRIPTION: "description",
|
||||
PLATE: "plate",
|
||||
WELL: "well",
|
||||
VLN: "vln",
|
||||
}
|
||||
|
||||
_NAME_RE = re.compile(
|
||||
r"^(?P<barcode>[^_/\\]+)_(?P<position>R\d{2}C\d{2})_(?P<channel>Grn|Red)\.idat(\.gz)?$", re.IGNORECASE
|
||||
)
|
||||
|
||||
_DATE_FORMATS = (
|
||||
"%m/%d/%Y %I:%M:%S %p",
|
||||
"%m/%d/%Y %H:%M:%S",
|
||||
"%Y-%m-%dT%H:%M:%S",
|
||||
"%Y-%m-%d %H:%M:%S",
|
||||
"%d.%m.%Y %H:%M:%S",
|
||||
)
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class RunInfoEntry:
|
||||
run_time: str
|
||||
block_type: str
|
||||
block_pars: str
|
||||
block_code: str
|
||||
code_version: str
|
||||
|
||||
|
||||
@dataclass
|
||||
class Idat:
|
||||
"""One colour channel of one BeadChip section."""
|
||||
|
||||
path: str
|
||||
version: int
|
||||
toc: dict[int, int]
|
||||
n_addresses: int
|
||||
address_ids: npt.NDArray[np.int32]
|
||||
mean: npt.NDArray[np.uint16]
|
||||
sd: npt.NDArray[np.uint16]
|
||||
n_beads: npt.NDArray[np.uint8]
|
||||
mid_block: npt.NDArray[Any] | None = None
|
||||
mid_block_itemsize: int | None = None
|
||||
run_info: list[RunInfoEntry] = field(default_factory=list)
|
||||
red_green: bytes = b""
|
||||
manifest: str = ""
|
||||
barcode: str = ""
|
||||
chip_type: str = ""
|
||||
position: str = ""
|
||||
beadset: str = ""
|
||||
sample_id: str = ""
|
||||
description: str = ""
|
||||
plate: str = ""
|
||||
well: str = ""
|
||||
vln: str = ""
|
||||
sample_count: int | None = None
|
||||
unknown_fields: dict[int, int] = field(default_factory=dict)
|
||||
|
||||
# -- derived metadata (IN-024) --------------------------------------------------------------------
|
||||
def scan_entry(self) -> RunInfoEntry | None:
|
||||
scans = [r for r in self.run_info if r.block_type.strip().lower() == "scan"]
|
||||
return scans[-1] if scans else None
|
||||
|
||||
@property
|
||||
def scan_time_raw(self) -> str:
|
||||
e = self.scan_entry()
|
||||
return e.run_time if e else ""
|
||||
|
||||
@property
|
||||
def scan_time_iso(self) -> str:
|
||||
return parse_run_time(self.scan_time_raw)
|
||||
|
||||
@property
|
||||
def scanner_software(self) -> str:
|
||||
e = self.scan_entry()
|
||||
return e.code_version if e else ""
|
||||
|
||||
@property
|
||||
def scanner_pars(self) -> str:
|
||||
e = self.scan_entry()
|
||||
return e.block_pars if e else ""
|
||||
|
||||
def metadata(self) -> dict[str, Any]:
|
||||
return {
|
||||
"path": self.path,
|
||||
"version": self.version,
|
||||
"n_addresses": self.n_addresses,
|
||||
"barcode": self.barcode,
|
||||
"position": self.position,
|
||||
"chip_type": self.chip_type,
|
||||
"manifest": self.manifest,
|
||||
"beadset": self.beadset,
|
||||
"sample_id": self.sample_id,
|
||||
"description": self.description,
|
||||
"plate": self.plate,
|
||||
"well": self.well,
|
||||
"sample_count": self.sample_count,
|
||||
"vln": self.vln,
|
||||
"red_green": self.red_green.hex(" "),
|
||||
"mid_block_count": None if self.mid_block is None else len(self.mid_block),
|
||||
"mid_block_itemsize": self.mid_block_itemsize,
|
||||
"scan_time": self.scan_time_iso,
|
||||
"scan_time_raw": self.scan_time_raw,
|
||||
"scanner_software": self.scanner_software,
|
||||
"toc_codes": sorted(self.toc),
|
||||
"unknown_codes": sorted(self.unknown_fields),
|
||||
"run_info": [r.__dict__ for r in self.run_info],
|
||||
}
|
||||
|
||||
|
||||
def parse_run_time(raw: str) -> str:
|
||||
"""Convert a RunInfo timestamp to ISO-8601 without consulting the process locale (EDGE-001).
|
||||
|
||||
Returns ``""`` if the value is empty or cannot be parsed.
|
||||
"""
|
||||
s = raw.strip()
|
||||
if not s:
|
||||
return ""
|
||||
for fmt in _DATE_FORMATS:
|
||||
try:
|
||||
return _dt.datetime.strptime(s, fmt).isoformat()
|
||||
except ValueError:
|
||||
continue
|
||||
return ""
|
||||
|
||||
|
||||
def parse_idat_filename(path: str | Path) -> tuple[str, str, str] | None:
|
||||
"""``<barcode>_<RxxCxx>_<Grn|Red>.idat[.gz]`` → (barcode, POSITION, 'Grn'|'Red'); None if not an IDAT name."""
|
||||
m = _NAME_RE.match(Path(path).name)
|
||||
if not m:
|
||||
return None
|
||||
ch = "Grn" if m.group("channel").lower() == "grn" else "Red"
|
||||
return m.group("barcode"), m.group("position").upper(), ch
|
||||
|
||||
|
||||
def read_idat(path: str | Path, *, load_arrays: bool = True) -> Idat:
|
||||
"""Parse an IDAT file. Raises FormatError (E_IDAT_FORMAT / E_IDAT_ENCRYPTED)."""
|
||||
spath = str(path)
|
||||
buf = load_bytes(path)
|
||||
r = BinaryReader(buf, "E_IDAT_FORMAT", spath)
|
||||
if r.size < 16 or r.raw(4, "magic") != MAGIC:
|
||||
raise r.fail("bad magic (not an IDAT file)")
|
||||
version = r.i64("version")
|
||||
# encrypted files carry version 1 in the low 32 bits followed by ciphertext (illuminaio reads int32)
|
||||
if version & 0xFFFFFFFF == 1:
|
||||
raise FormatError(
|
||||
"E_IDAT_ENCRYPTED",
|
||||
f"{spath}: encrypted (version 1) IDAT files are not supported",
|
||||
path=spath,
|
||||
hint="convert with an external tool first",
|
||||
)
|
||||
if version != 3:
|
||||
raise r.fail(f"unsupported IDAT version {version}")
|
||||
n_fields = r.i32("TOC size")
|
||||
if n_fields <= 0 or n_fields > 4096:
|
||||
raise r.fail(f"implausible TOC size {n_fields}")
|
||||
toc: dict[int, int] = {}
|
||||
for _ in range(n_fields):
|
||||
code = r.u16("TOC code")
|
||||
off = r.i64("TOC offset")
|
||||
if off < 0 or off >= r.size:
|
||||
raise r.fail(f"TOC offset {off} for field {code} outside file (size {r.size})")
|
||||
if code in toc:
|
||||
raise r.fail(f"duplicate TOC field {code}")
|
||||
toc[code] = off
|
||||
|
||||
for required in (N_SNPS, ILLUMINA_ID, MEAN, SD, NBEADS):
|
||||
if required not in toc:
|
||||
raise r.fail(f"missing required field {required} ({FIELD_NAMES[required]})")
|
||||
|
||||
r.seek(toc[N_SNPS])
|
||||
n = r.i32("number of addresses")
|
||||
if n < 0:
|
||||
raise r.fail(f"negative address count {n}")
|
||||
|
||||
empty_i32 = np.empty(0, dtype="<i4")
|
||||
idat = Idat(
|
||||
path=spath,
|
||||
version=version,
|
||||
toc=dict(toc),
|
||||
n_addresses=n,
|
||||
address_ids=empty_i32,
|
||||
mean=np.empty(0, dtype="<u2"),
|
||||
sd=np.empty(0, dtype="<u2"),
|
||||
n_beads=np.empty(0, dtype="u1"),
|
||||
)
|
||||
|
||||
if load_arrays:
|
||||
r.seek(toc[ILLUMINA_ID])
|
||||
idat.address_ids = r.array("<i4", n, "address ids")
|
||||
r.seek(toc[SD])
|
||||
idat.sd = r.array("<u2", n, "SD")
|
||||
r.seek(toc[MEAN])
|
||||
idat.mean = r.array("<u2", n, "Mean")
|
||||
r.seek(toc[NBEADS])
|
||||
idat.n_beads = r.array("u1", n, "NBeads")
|
||||
if has_duplicates(idat.address_ids):
|
||||
raise r.fail("duplicate address IDs")
|
||||
|
||||
offsets_sorted = sorted(set(toc.values()) | {r.size})
|
||||
|
||||
def section_end(off: int) -> int:
|
||||
return next(o for o in offsets_sorted if o > off)
|
||||
|
||||
for code, off in sorted(toc.items(), key=lambda kv: kv[1]):
|
||||
r.seek(off)
|
||||
if code in _STRING_FIELDS:
|
||||
setattr(idat, _STRING_FIELDS[code], r.string(FIELD_NAMES[code]))
|
||||
elif code == RUN_INFO:
|
||||
m = r.i32("run info count")
|
||||
if m < 0 or m > 100_000:
|
||||
raise r.fail(f"implausible run info count {m}")
|
||||
idat.run_info = [RunInfoEntry(*(r.string("run info") for _ in range(5))) for _ in range(m)]
|
||||
elif code == RED_GREEN:
|
||||
idat.red_green = r.raw(4, "RedGreen")
|
||||
elif code == SAMPLE_COUNT:
|
||||
idat.sample_count = r.i32("sample count")
|
||||
elif code == MID_BLOCK:
|
||||
if not load_arrays:
|
||||
continue
|
||||
count = r.i32("mid block count")
|
||||
space = section_end(off) - off - 4
|
||||
if count < 0:
|
||||
raise r.fail(f"negative mid block count {count}")
|
||||
# element size inferred from the space before the next section (int32 vs uint8, plan §20 U2)
|
||||
if count == space and count * 4 != space:
|
||||
itemsize = 1
|
||||
elif count * 4 <= space:
|
||||
itemsize = 4
|
||||
elif count <= space:
|
||||
itemsize = 1
|
||||
else:
|
||||
raise r.fail(f"mid block of {count} entries does not fit in {space} bytes")
|
||||
idat.mid_block_itemsize = itemsize
|
||||
idat.mid_block = r.array("<i4" if itemsize == 4 else "u1", count, "mid block")
|
||||
elif code in (N_SNPS, ILLUMINA_ID, SD, MEAN, NBEADS):
|
||||
continue
|
||||
else:
|
||||
idat.unknown_fields[code] = off
|
||||
return idat
|
||||
|
||||
|
||||
def has_duplicates(ids: npt.NDArray[np.int32]) -> bool:
|
||||
"""Duplicate check; O(n) when IDs are already strictly increasing (the usual case)."""
|
||||
if len(ids) < 2:
|
||||
return False
|
||||
d = np.diff(ids)
|
||||
if (d > 0).all():
|
||||
return False
|
||||
s = np.sort(ids)
|
||||
return bool((s[1:] == s[:-1]).any())
|
||||
|
||||
|
||||
def channel_of(idat: Idat) -> str | None:
|
||||
parsed = parse_idat_filename(idat.path)
|
||||
return parsed[2] if parsed else None
|
||||
|
||||
|
||||
def check_pair(grn: Idat, red: Idat) -> None:
|
||||
"""Red/Grn consistency (IN-023): same N and the same address set. Raises E_IDAT_PAIR_MISMATCH."""
|
||||
if grn.n_addresses != red.n_addresses:
|
||||
raise FormatError(
|
||||
"E_IDAT_PAIR_MISMATCH",
|
||||
f"Grn has {grn.n_addresses} addresses, Red has {red.n_addresses}: {grn.path}, {red.path}",
|
||||
)
|
||||
if len(grn.address_ids) and not np.array_equal(np.sort(grn.address_ids), np.sort(red.address_ids)):
|
||||
raise FormatError("E_IDAT_PAIR_MISMATCH", f"address sets differ: {grn.path}, {red.path}")
|
||||
@@ -0,0 +1,359 @@
|
||||
"""Sample-sheet parser and validator — spec IN-001 … IN-017 (plan §6.6).
|
||||
|
||||
A GenomeStudio-compatible CSV with optional ``[Header]``, ``[Manifests]`` (ignored) and ``[Data]`` sections.
|
||||
All column matching uses Unicode case folding (locale independent, EDGE-001).
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import csv
|
||||
import io
|
||||
import re
|
||||
from dataclasses import dataclass, field
|
||||
from pathlib import Path
|
||||
|
||||
from simurg.core.errors import SimurgError
|
||||
|
||||
MAX_LEN = 500
|
||||
KEY_PAIRS = (("SentrixBarcode_A", "SentrixPosition_A"), ("beadChipName", "sampleSectionName"))
|
||||
RECOGNISED = (
|
||||
"Sample_ID",
|
||||
"Sample_Name",
|
||||
"Sample_Plate",
|
||||
"Sample_Well",
|
||||
"Sample_Group",
|
||||
"Gender",
|
||||
"Sex",
|
||||
"demographicSex",
|
||||
"Parent1",
|
||||
"Parent2",
|
||||
"Replicate",
|
||||
)
|
||||
SEX_COLUMNS_PRIORITY = ("Sex", "Gender", "demographicSex")
|
||||
|
||||
# Output columns of sample_summary (OUT-060) and their compatibility spellings (OUT-066): a sheet column may
|
||||
# not collide with these case-insensitively (IN-005). Recognised input columns above are exempt.
|
||||
RESERVED_OUTPUT_COLUMNS = (
|
||||
"sample_id",
|
||||
"sample_name",
|
||||
"sample_folder",
|
||||
"status",
|
||||
"fail_reason",
|
||||
"autosomal_call_rate",
|
||||
"call_rate",
|
||||
"logr_dev",
|
||||
"sex_est",
|
||||
"sex_reported",
|
||||
"sex_check",
|
||||
"gc10",
|
||||
"gc50",
|
||||
"n_calls",
|
||||
"n_nocalls",
|
||||
"n_intensity_only",
|
||||
"p05_x",
|
||||
"p50_x",
|
||||
"p95_x",
|
||||
"p05_y",
|
||||
"p50_y",
|
||||
"p95_y",
|
||||
"het_rate",
|
||||
"contamination",
|
||||
"functional_status",
|
||||
"control_status",
|
||||
"control_flag_codes",
|
||||
"flags",
|
||||
"Sample ID",
|
||||
"Sample Name",
|
||||
"Sample Folder",
|
||||
"Autosomal Call Rate",
|
||||
"Call Rate",
|
||||
"Log R Ratio Std Dev",
|
||||
"Sex Estimate",
|
||||
)
|
||||
|
||||
_BARCODE_RE = re.compile(r"^[0-9]{10,14}$")
|
||||
_BARCODE_DECIMAL_RE = re.compile(r"^([0-9]{10,14})\.0+$")
|
||||
_BARCODE_SCI_RE = re.compile(r"^[0-9][.,][0-9]+[eE][+]?[0-9]+$")
|
||||
_POSITION_RE = re.compile(r"^[Rr][0-9]{2}[Cc][0-9]{2}$")
|
||||
_NORM_RE = re.compile(r".*_norm_r$", re.IGNORECASE)
|
||||
_NORM_COMPAT_RE = re.compile(r".* norm r$", re.IGNORECASE)
|
||||
|
||||
|
||||
def fold(s: str) -> str:
|
||||
return s.strip().casefold()
|
||||
|
||||
|
||||
@dataclass
|
||||
class SheetRow:
|
||||
row_number: int # 1-based line number in the file
|
||||
barcode: str
|
||||
position: str
|
||||
values: dict[str, str] # original header -> value (empty-header columns dropped)
|
||||
sample_name: str | None = None
|
||||
reported_sex: str = "U"
|
||||
|
||||
@property
|
||||
def sample_id(self) -> str:
|
||||
return f"{self.barcode}_{self.position}"
|
||||
|
||||
|
||||
@dataclass
|
||||
class SampleSheet:
|
||||
path: str
|
||||
header: dict[str, str]
|
||||
root_folder: Path | None
|
||||
columns: list[str]
|
||||
rows: list[SheetRow]
|
||||
barcode_column: str
|
||||
position_column: str
|
||||
warnings: list[tuple[str, str]] = field(default_factory=list)
|
||||
|
||||
@property
|
||||
def has_sample_name(self) -> bool:
|
||||
return any(fold(c) == fold("Sample_Name") for c in self.columns)
|
||||
|
||||
def passthrough_columns(self) -> list[str]:
|
||||
"""Every non-key column in sheet order (IN-013)."""
|
||||
keys = {fold(self.barcode_column), fold(self.position_column)}
|
||||
return [c for c in self.columns if fold(c) not in keys]
|
||||
|
||||
def by_sample_id(self) -> dict[str, SheetRow]:
|
||||
return {r.sample_id: r for r in self.rows}
|
||||
|
||||
|
||||
def _err(code: str, path: str, msg: str, hint: str | None = None) -> SimurgError:
|
||||
return SimurgError(code, f"{path}: {msg}", hint=hint, path=path)
|
||||
|
||||
|
||||
def normalise_sex(value: str) -> str:
|
||||
v = fold(value)
|
||||
if v in ("m", "male", "1"):
|
||||
return "M"
|
||||
if v in ("f", "female", "2"):
|
||||
return "F"
|
||||
return "U"
|
||||
|
||||
|
||||
def validate_barcode(value: str, where: str, warnings: list[tuple[str, str]], path: str) -> str:
|
||||
v = value.strip()
|
||||
if _BARCODE_RE.match(v):
|
||||
return v
|
||||
m = _BARCODE_DECIMAL_RE.match(v)
|
||||
if m:
|
||||
warnings.append(("W_SHEET_BARCODE_REPAIRED", f"{where}: barcode {v!r} repaired to {m.group(1)}"))
|
||||
return m.group(1)
|
||||
hint = None
|
||||
if _BARCODE_SCI_RE.match(v):
|
||||
hint = "format the SentrixBarcode_A column as a number with 0 decimals (or as text) before saving"
|
||||
raise _err("E_SHEET_BARCODE_FORMAT", path, f"{where}: invalid barcode {v!r}", hint)
|
||||
|
||||
|
||||
def validate_position(value: str, where: str, path: str) -> str:
|
||||
v = value.strip()
|
||||
if not _POSITION_RE.match(v):
|
||||
raise _err("E_SHEET_BARCODE_FORMAT", path, f"{where}: invalid position {v!r} (expected RxxCxx)")
|
||||
return v.upper()
|
||||
|
||||
|
||||
def _read_rows(path: Path) -> list[tuple[int, list[str]]]:
|
||||
raw = path.read_bytes()
|
||||
text = raw.decode("utf-8-sig", errors="replace")
|
||||
stripped = text.lstrip()
|
||||
if stripped.startswith("{") or (stripped.startswith("[") and stripped[1:].lstrip()[:1] in ("{", '"', "]")):
|
||||
raise _err(
|
||||
"E_SHEET_DELIMITER",
|
||||
str(path),
|
||||
"looks like JSON, not a CSV sample sheet",
|
||||
"export the sheet as comma-separated CSV",
|
||||
)
|
||||
rows: list[tuple[int, list[str]]] = []
|
||||
reader = csv.reader(io.StringIO(text, newline=""))
|
||||
for row in reader:
|
||||
rows.append((reader.line_num, row))
|
||||
return rows
|
||||
|
||||
|
||||
def parse_samplesheet(path: str | Path, *, on_duplicate_sample: str = "error") -> SampleSheet:
|
||||
"""Parse and validate a sample sheet. Raises SimurgError with an E_SHEET_* code."""
|
||||
p = Path(path)
|
||||
spath = str(p)
|
||||
warnings: list[tuple[str, str]] = []
|
||||
rows = _read_rows(p)
|
||||
|
||||
header_kv: dict[str, str] = {}
|
||||
data_rows: list[tuple[int, list[str]]] = []
|
||||
section: str | None = None
|
||||
saw_section = False
|
||||
for line_no, row in rows:
|
||||
if not row or all(not c.strip() for c in row):
|
||||
continue
|
||||
first = row[0].strip()
|
||||
if first.startswith("[") and first.endswith("]") and first[1:-1].strip().replace(" ", "").isalnum():
|
||||
section = fold(first[1:-1])
|
||||
saw_section = True
|
||||
continue
|
||||
if not saw_section or section == "data":
|
||||
data_rows.append((line_no, row))
|
||||
elif section == "header":
|
||||
header_kv[first] = row[1].strip() if len(row) > 1 else ""
|
||||
# [Manifests] and unknown sections are ignored
|
||||
|
||||
if not data_rows:
|
||||
raise _err("E_SHEET_EMPTY", spath, "no [Data] header row")
|
||||
head_line, head = data_rows[0]
|
||||
if len(head) == 1 and (";" in head[0] or "\t" in head[0]):
|
||||
raise _err(
|
||||
"E_SHEET_DELIMITER",
|
||||
spath,
|
||||
f"line {head_line}: header is not comma-delimited",
|
||||
"save the sheet as comma-separated CSV (not semicolon or tab)",
|
||||
)
|
||||
|
||||
# columns: drop empty headers (IN-004), check length (IN-006) and case-insensitive duplicates (IN-003)
|
||||
col_idx: list[tuple[int, str]] = []
|
||||
empty_cols = 0
|
||||
seen: dict[str, str] = {}
|
||||
for i, name in enumerate(head):
|
||||
nm = name.strip()
|
||||
if not nm:
|
||||
empty_cols += 1
|
||||
continue
|
||||
if len(nm) > MAX_LEN:
|
||||
raise _err(
|
||||
"E_SHEET_TOO_LONG", spath, f"column name at position {i + 1} has {len(nm)} characters (limit {MAX_LEN})"
|
||||
)
|
||||
k = fold(nm)
|
||||
if k in seen:
|
||||
raise _err(
|
||||
"E_SHEET_DUP_COLUMN",
|
||||
spath,
|
||||
f"Duplicate column '{nm}' (column names are case-insensitive). "
|
||||
f"Rename or remove one of: {seen[k]}, {nm}",
|
||||
)
|
||||
seen[k] = nm
|
||||
col_idx.append((i, nm))
|
||||
if empty_cols:
|
||||
warnings.append(("W_SHEET_EMPTY_COLUMN", f"{empty_cols} column(s) with an empty header ignored"))
|
||||
columns = [nm for _, nm in col_idx]
|
||||
folded = {fold(c): c for c in columns}
|
||||
|
||||
key: tuple[str, str] | None = None
|
||||
for bc, pc in KEY_PAIRS:
|
||||
if fold(bc) in folded and fold(pc) in folded:
|
||||
key = (folded[fold(bc)], folded[fold(pc)])
|
||||
break
|
||||
if key is None:
|
||||
raise _err(
|
||||
"E_SHEET_NO_KEY_COLUMNS",
|
||||
spath,
|
||||
"missing SentrixBarcode_A/SentrixPosition_A (or beadChipName/sampleSectionName) columns",
|
||||
)
|
||||
|
||||
recognised = {fold(c) for c in RECOGNISED} | {fold(key[0]), fold(key[1])}
|
||||
reserved = {fold(c) for c in RESERVED_OUTPUT_COLUMNS}
|
||||
for c in columns:
|
||||
k = fold(c)
|
||||
if k in recognised:
|
||||
continue
|
||||
if k in reserved or _NORM_RE.match(c) or _NORM_COMPAT_RE.match(c):
|
||||
raise _err(
|
||||
"E_SHEET_RESERVED",
|
||||
spath,
|
||||
f"column {c!r} is reserved for SimurgArray output",
|
||||
"rename the column (or use --namespace-sheet-columns)",
|
||||
)
|
||||
|
||||
out_rows: list[SheetRow] = []
|
||||
ids: dict[str, int] = {}
|
||||
sex_cols = [folded[fold(c)] for c in SEX_COLUMNS_PRIORITY if fold(c) in folded]
|
||||
name_col = folded.get(fold("Sample_Name"))
|
||||
for line_no, row in data_rows[1:]:
|
||||
values: dict[str, str] = {}
|
||||
for i, nm in col_idx:
|
||||
v = row[i] if i < len(row) else ""
|
||||
if len(v) > MAX_LEN:
|
||||
raise _err(
|
||||
"E_SHEET_TOO_LONG",
|
||||
spath,
|
||||
f"line {line_no}, column {nm!r}: value has {len(v)} characters (limit {MAX_LEN})",
|
||||
)
|
||||
values[nm] = v
|
||||
where = f"line {line_no}"
|
||||
barcode = validate_barcode(values[key[0]], where, warnings, spath)
|
||||
position = validate_position(values[key[1]], where, spath)
|
||||
sr = SheetRow(line_no, barcode, position, values)
|
||||
if name_col is not None:
|
||||
sr.sample_name = values[name_col].strip()
|
||||
sexes = [normalise_sex(values[c]) for c in sex_cols]
|
||||
known = [s for s in sexes if s != "U"]
|
||||
if len(set(known)) > 1:
|
||||
warnings.append(("W_SHEET_SEX_CONFLICT", f"{where}: reported-sex columns disagree; using {sex_cols[0]}"))
|
||||
sr.reported_sex = sexes[0] if sexes else "U"
|
||||
if sr.sample_id in ids:
|
||||
if on_duplicate_sample == "drop":
|
||||
warnings.append(
|
||||
(
|
||||
"E_SHEET_DUP_SAMPLE",
|
||||
f"{where}: duplicate sample {sr.sample_id} dropped (first at line {ids[sr.sample_id]})",
|
||||
)
|
||||
)
|
||||
continue
|
||||
raise _err(
|
||||
"E_SHEET_DUP_SAMPLE",
|
||||
spath,
|
||||
f"{where}: sample {sr.sample_id} already listed at line {ids[sr.sample_id]}",
|
||||
)
|
||||
ids[sr.sample_id] = line_no
|
||||
out_rows.append(sr)
|
||||
|
||||
if not out_rows:
|
||||
raise _err("E_SHEET_EMPTY", spath, "no data rows")
|
||||
|
||||
if name_col is not None:
|
||||
names: dict[str, int] = {}
|
||||
for sr in out_rows:
|
||||
if not sr.sample_name:
|
||||
raise _err(
|
||||
"E_SHEET_NAME_EMPTY",
|
||||
spath,
|
||||
f"line {sr.row_number}: Sample_Name is empty (it must be set in every row when the column exists)",
|
||||
)
|
||||
if sr.sample_name in ids and sr.sample_name != sr.sample_id:
|
||||
raise _err(
|
||||
"E_SHEET_NAME_COLLISION",
|
||||
spath,
|
||||
f"line {sr.row_number}: Sample_Name {sr.sample_name!r} equals the ID of another sample",
|
||||
)
|
||||
if sr.sample_name in names:
|
||||
raise _err(
|
||||
"E_SHEET_NAME_COLLISION",
|
||||
spath,
|
||||
f"line {sr.row_number}: Sample_Name "
|
||||
f"{sr.sample_name!r} already used at line {names[sr.sample_name]}",
|
||||
)
|
||||
names[sr.sample_name] = sr.row_number
|
||||
|
||||
root: Path | None = None
|
||||
for k, v in header_kv.items():
|
||||
if fold(k) == "rootfolder" and v:
|
||||
rp = Path(v)
|
||||
root = rp if rp.is_absolute() else (p.parent / rp)
|
||||
return SampleSheet(
|
||||
path=spath,
|
||||
header=header_kv,
|
||||
root_folder=root,
|
||||
columns=columns,
|
||||
rows=out_rows,
|
||||
barcode_column=key[0],
|
||||
position_column=key[1],
|
||||
warnings=warnings,
|
||||
)
|
||||
|
||||
|
||||
_STEM_BAD = re.compile(r"[^A-Za-z0-9._-]")
|
||||
|
||||
|
||||
def output_stem(name: str) -> str:
|
||||
"""File-name stem for a sample (IN-012)."""
|
||||
s = _STEM_BAD.sub("_", name).lstrip(".")
|
||||
return (s or "_")[:200]
|
||||
Reference in new issue
Block a user