diff --git a/CHANGELOG.md b/CHANGELOG.md index 12e2e75..e022a20 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,31 @@ # Changelog +## 0.7.2 + +- FCS: the TEXT splitter now tells a legacy empty value from an escaped + delimiter by position. BD CellQuest and other FCS 2.0 / 3.0 writers put an + empty value down as a doubled delimiter (`KEY\\NEXT\value\`); read as an + escape, KEY swallowed NEXT and every later pair swapped keyword and value. + A keyword or value may not begin with the delimiter (FCS 3.1 3.2.7), so a + doubled delimiter straight after a keyword is an empty value, and inside a + value it is still a literal (`$P3F/488//10/` is `488/10`). On flowio's own + `data1.fcs` this reads the 4 empty CellQuest pairs that flowio 1.4.0 itself + misreads; the other 31 corpus files give identical keywords. A keyword that + contains the delimiter is the one layout this cannot read, and it costs + only that pair. +- FCS: a supplemental TEXT segment that lies within the leading read is parsed + too, so `$CYT` and `$PnS` placed there (optional keywords may be) are found. + Primary TEXT keywords win. +- FCS: a channel without `$PnN` takes its `$PnS` as the name; when `$PAR` is + missing or unusable, channels are listed from the `$PnN` keywords present. +- FCS: values that are not UTF-8 (CellQuest's Latin-1 `\xaa`) decode as + Latin-1 instead of becoming replacement characters; `$DATATYPE` and `$CYT` + are stripped and `$DATATYPE` upper-cased; a HEADER offset with a stray byte + after the digits is read; `multipleDatasets` is set when `$NEXTDATA` points + at a second dataset (Coulter LMD files hold an FCS 2.0 and an FCS 3.0 copy). +- FCS: the splitter works on `bytes.split` runs instead of a per-byte loop + (1 MiB of TEXT in 13 ms instead of ~300 ms). + ## 0.7.1 Fixes found by cross-checking these decoders against other readers. diff --git a/pyproject.toml b/pyproject.toml index c943efb..6a541b7 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "scigantic-headers" -version = "0.7.1" +version = "0.7.2" description = "Decode scientific file headers (cryo-EM, light microscopy and light-sheet imaging incl. TIFF/OME-TIFF, CZI, HDF5/Imaris/BigDataViewer, LIF, ND2, Zarr/N5; tabular, flow cytometry, mass spec, sequencing, structural biology) into typed fields. Pure, zero-dependency, with bounded readers that follow internal pointers without downloading the file." readme = "README.md" requires-python = ">=3.9" diff --git a/src/scigantic_headers/fcs.py b/src/scigantic_headers/fcs.py index 1e17422..87b4f24 100644 --- a/src/scigantic_headers/fcs.py +++ b/src/scigantic_headers/fcs.py @@ -15,9 +15,18 @@ The TEXT segment is a run of delimiter-separated keyword/value pairs. Its first byte is the delimiter; the rest is `$KEY value $KEY ...`. -A doubled delimiter is an escaped literal delimiter inside a value. The keywords -this reads: $PAR (parameters), $TOT (events), $DATATYPE, $MODE, $CYT (cytometer), -and $PnN / $PnS (per-channel short and long names). +A doubled delimiter inside a value is an escaped literal delimiter (FCS 3.1 +section 3.2.7). The keywords this reads: $PAR (parameters), $TOT (events), +$DATATYPE, $MODE, $CYT (cytometer), $NEXTDATA, $BEGINSTEXT/$ENDSTEXT, and +$PnN / $PnS (per-channel short and long names). + +Legacy writers (BD CellQuest and other FCS 2.0 / 3.0 vendors) write an empty +value as a doubled delimiter: `KEY\\\\NEXTKEY\\value\\`. Read strictly, that +doubled delimiter is a literal, KEY's value swallows NEXTKEY, and every later +pair swaps keyword and value. The two readings can be told apart by position: +a keyword or value may not begin with the delimiter (section 3.2.7), so a +doubled delimiter that directly follows a keyword cannot start its value and +must be the empty value. Inside a value it is a literal. See `_split`. The TEXT segment sits near the start but can run well past 1 KiB, so this decoder registers a larger leading read (see the register_decoder call below). @@ -25,55 +34,129 @@ from __future__ import annotations -from typing import List, Optional +import re +from typing import Dict, List, Optional from .decoders import DecodedHeader, Read, register_decoder _FCS_VERSIONS = (b"FCS2.0", b"FCS3.0", b"FCS3.1", b"FCS3.2") _HEADER_MIN = 58 # version + reserved + six 8-byte offset fields _DATATYPE = {"I": "int", "F": "float32", "D": "float64", "A": "ascii"} +_MAX_PAR = 100000 +_PNN = re.compile(r"^\$P(\d+)N$") def _offset(data: bytes, start: int) -> int: - """One 8-byte ASCII offset field. -1 if it is not a non-negative integer.""" + """One 8-byte ASCII offset field. -1 if it holds no non-negative integer. + Blank means 0 (the spec's value for a segment that is absent or too far + into the file for the HEADER). A stray trailing byte after the digits, seen + in some LMD exports, is ignored.""" s = data[start:start + 8].strip() if not s: return 0 try: return int(s) except ValueError: - return -1 + m = re.match(rb"\d+", s) + return int(m.group()) if m else -1 def _split(body: bytes, delim: int) -> List[bytes]: - """Split the TEXT body on the delimiter, treating a doubled delimiter as an - escaped literal delimiter within a value (the FCS convention).""" + """Split a TEXT body (the bytes after the leading delimiter) into alternating + keyword and value fields. + + Delimiters come in runs. A run that follows a value is `run // 2` literal + delimiters, plus one separator if the run is odd (spec escaping: a doubled + delimiter is a literal). A run of two or more that follows a keyword is the + keyword's separator and an empty value: a value may not begin with the + delimiter, so it cannot be a literal there. That reads spec-escaped values + such as `$P3F/488//10/` (value `488/10`) and legacy empty values such as + `KEY\\\\NEXT\\v\\` (KEY empty, NEXT = v) from the same text. A keyword that + itself contains an escaped delimiter (a `CD GAG/100` label written with `/` + as the delimiter) reads as two keywords around an empty value; that misreads + only that pair, where the strict reading of a legacy empty value shifts + every pair after it.""" + d = bytes((delim,)) + pieces = body.lstrip(d).split(d) # a leading delimiter run is junk: no field may start with one out: List[bytes] = [] - cur = bytearray() - i, n = 0, len(body) - while i < n: - if body[i] == delim: - if i + 1 < n and body[i + 1] == delim: # escaped literal delimiter - cur.append(delim) - i += 2 - continue - out.append(bytes(cur)) - cur = bytearray() - i += 1 + cur = bytearray(pieces[0]) + npieces = len(pieces) + k = 1 + while k < npieces: + j = k + while j < npieces and not pieces[j]: + j += 1 + if j < npieces: + run = j - k + 1 # delimiters between cur and pieces[j] + nxt = pieces[j] else: - cur.append(body[i]) - i += 1 + run = j - k # delimiters after cur at the end of the body + nxt = b"" + while run: + if len(out) % 2 == 0: # cur is a keyword + out.append(bytes(cur)) + cur = bytearray() + if run >= 2: # keyword, then an empty value (legacy writers) + out.append(b"") + run -= 2 + else: + run -= 1 + else: # cur is a value: pairs are literals, an odd one ends it + cur.extend(d * (run // 2)) + if run % 2: + out.append(bytes(cur)) + cur = bytearray() + run = 0 + cur.extend(nxt) + k = j + 1 if cur: out.append(bytes(cur)) return out +def _text(v: bytes) -> str: + """A keyword value as text. FCS 3.1 values are UTF-8; older files from + Windows instruments carry Latin-1 bytes (CellQuest's '\\xaa'), which are + kept rather than replaced.""" + try: + return v.decode("utf-8") + except UnicodeDecodeError: + return v.decode("latin-1") + + +def _pairs(body: bytes, delim: int, *, drop_last: bool = False) -> Dict[str, str]: + """Parse a TEXT body into a dict. Keys are upper-cased and stripped, empty + keys are skipped, and a repeated key keeps its last value (as flowio and + fcsparser do). drop_last discards the final field, which may have been cut + by a short read.""" + parts = _split(body, delim) + if drop_last and parts: + parts = parts[:-1] + kv: Dict[str, str] = {} + for i in range(0, len(parts) - 1, 2): + key = parts[i].decode("ascii", "replace").strip().upper() + if key: + kv[key] = _text(parts[i + 1]) + return kv + + +def _as_int(kv: Dict[str, str], key: str) -> Optional[int]: + v = kv.get(key) + if v is None: + return None + try: + return int(v.strip()) + except ValueError: + return None + + def decode_fcs(data: bytes) -> Optional[DecodedHeader]: """Decode an FCS file from its leading bytes (version, HEADER offsets, and the TEXT segment). Returns None if the magic is wrong. A TEXT segment longer than `data` (many parameters, a large spillover matrix) is decoded as far as it goes: the last part may be cut mid-value, so it is dropped, and - textTruncated says keywords may be missing.""" + textTruncated says keywords may be missing. A supplemental TEXT segment + that lies within `data` is read too; primary TEXT keywords win.""" if len(data) < _HEADER_MIN or data[:6] not in _FCS_VERSIONS: return None version = data[:6].decode("ascii") @@ -89,56 +172,63 @@ def decode_fcs(data: bytes) -> Optional[DecodedHeader]: truncated = True text_end = len(data) - 1 - body = data[text_begin + 1:text_end + 1] # skip the leading delimiter byte delim = data[text_begin] - parts = _split(body, delim) - if truncated and parts: - parts = parts[:-1] - - kv = {} - it = iter(parts) - for key in it: - try: - val = next(it) - except StopIteration: - break - kv[key.decode("ascii", "replace").strip().upper()] = val.decode("utf-8", "replace") - - def as_int(key: str) -> Optional[int]: - v = kv.get(key) - if v is None: - return None - try: - return int(v.strip()) - except ValueError: - return None - - npar = as_int("$PAR") - ntot = as_int("$TOT") - + kv = _pairs(data[text_begin + 1:text_end + 1], delim, drop_last=truncated) + + # Supplemental TEXT: optional keywords only ($CYT, $PnS, spillover...), same + # delimiter (3.2.5). Some writers start it with the delimiter, some do not. + stext_begin = _as_int(kv, "$BEGINSTEXT") or 0 + stext_end = _as_int(kv, "$ENDSTEXT") or 0 + if ( + 0 < stext_begin < stext_end < len(data) + and (stext_begin, stext_end) != (text_begin, text_end) + and not truncated + ): + seg = data[stext_begin:stext_end + 1] + if seg[:1] == bytes((delim,)): + seg = seg[1:] + for key, val in _pairs(seg, delim).items(): + kv.setdefault(key, val) # primary TEXT wins + + npar = _as_int(kv, "$PAR") + ntot = _as_int(kv, "$TOT") + + # Channel count: $PAR, or when it is missing or unusable, the highest $PnN + # index present. A channel without $PnN (seen in FCS 2.0) takes its $PnS. + count = npar if npar is not None and 0 < npar <= _MAX_PAR else 0 + if count == 0: + for key in kv: + m = _PNN.match(key) + if m: + count = max(count, min(int(m.group(1)), _MAX_PAR)) channels = [] - if npar is not None and 0 < npar <= 100000: - for i in range(1, npar + 1): - name = kv.get("$P%dN" % i) - if name is None: + for i in range(1, count + 1): + name = kv.get("$P%dN" % i) + long = kv.get("$P%dS" % i) + if name is None: + if not long: continue - ch = {"name": name} - long = kv.get("$P%dS" % i) - if long: - ch["label"] = long - channels.append(ch) + name, long = long, None + ch = {"name": name} + if long: + ch["label"] = long + channels.append(ch) + datatype = (kv.get("$DATATYPE") or "").strip().upper() fields = { "version": version, "numParameters": npar, "numEvents": ntot, - "dataType": _DATATYPE.get(kv.get("$DATATYPE", ""), kv.get("$DATATYPE")), + "dataType": _DATATYPE.get(datatype, kv.get("$DATATYPE")), "mode": kv.get("$MODE"), "channels": channels, } - cyt = kv.get("$CYT") + cyt = (kv.get("$CYT") or "").strip() if cyt: fields["cytometer"] = cyt + next_data = _as_int(kv, "$NEXTDATA") + if next_data: + fields["multipleDatasets"] = True # a Coulter LMD holds an FCS 2.0 and an FCS 3.0 copy if truncated: fields["textTruncated"] = True diff --git a/tests/test_fcs.py b/tests/test_fcs.py index 25a5ca1..375b577 100644 --- a/tests/test_fcs.py +++ b/tests/test_fcs.py @@ -235,3 +235,173 @@ def test_matches_flowio_channel_and_event_counts(): # channel short names ($PnN) agree ref_names = [fd.channels[i]["pnn"] for i in range(1, fd.channel_count + 1)] assert [c["name"] for c in h.fields["channels"]] == ref_names + + +# ── TEXT segment edge cases: legacy empty values, escapes, encodings, STEXT ─── +# +# These build the TEXT bytes by hand so the exact delimiter runs are under test. + +def raw_fcs(text, version=b"FCS3.1", tail=b""): + """HEADER + the given TEXT bytes (+ tail), with no DATA segment.""" + text_begin = 58 + text_end = text_begin + len(text) - 1 + header = version + b" " + b"%8d%8d" % (text_begin, text_end) + b"%8d" % 0 * 4 + assert len(header) == 58 + return header + text + tail + + +def test_legacy_empty_value_does_not_shift_later_pairs(): + # BD CellQuest (FCS 2.0) writes an empty value as a doubled delimiter: + # `KEY\\NEXT\value\`. Read as an escaped literal, KEY's value would swallow + # NEXT and every later keyword would become a value. Position decides: a + # doubled delimiter straight after a keyword is an empty value. + text = ( + b"\\$PAR\\2\\$TOT\\5\\$DATATYPE\\I\\$MODE\\L" + b"\\&5Data File Prefix Part #1\\\\&6Data File Prefix Part #2\\\\" + b"&8Acquisition Doc.\\LYMPH SUBSET ACQ\\$CYT\\FACScan\\$P1N\\FSC-H\\$P2N\\SSC-H" + b"\\&13Analysis Doc.\\\\" + ) + h = decode_fcs(raw_fcs(text, version=b"FCS2.0")) + assert h.fields["cytometer"] == "FACScan" + assert [c["name"] for c in h.fields["channels"]] == ["FSC-H", "SSC-H"] + assert h.fields["numEvents"] == 5 + + +def test_doubled_delimiter_inside_value_is_a_literal(): + # Spec escaping (3.2.7): `$P3F/488//10/` is the filter "488/10". A value may + # also end in a literal (`C:/data/`, a run of three), and four in a row are + # two literals (a Miltenyi writer's `561////10 nm`). + text = ( + b"/$PAR/2/$TOT/1/$P1N/FSC/$P1S/488//10/$P2N/SSC/$P2S/561////10 nm/" + b"$CYT/C://data///$MODE/L/" + ) + h = decode_fcs(raw_fcs(text)) + assert h.fields["channels"][0]["label"] == "488/10" + assert h.fields["channels"][1]["label"] == "561//10 nm" + assert h.fields["cytometer"] == "C:/data/" + assert h.fields["mode"] == "L" + + +def test_legacy_empty_value_and_escaped_literal_in_one_text(): + text = b"/$PAR/1/$TOT/3/NOTE//$P1N/FSC/$P1S/a//b/$CYT/Cyto/" + h = decode_fcs(raw_fcs(text)) + assert h.fields["channels"] == [{"name": "FSC", "label": "a/b"}] + assert h.fields["cytometer"] == "Cyto" + assert h.fields["numEvents"] == 3 + + +def test_split_round_trips_random_texts(): + # Property: any keywords without the delimiter and any values that do not + # start with it (the spec's rule), empty values written the legacy way and + # literals doubled, come back exactly. + import random + + from scigantic_headers.fcs import _split + + rng = random.Random(2026) + for delim in (b"/", b"\\", b"|", b"\x0c"): + alphabet = [b"a", b"B", b"1", b" ", b"$", delim] + for _ in range(300): + pairs = [] + for _ in range(rng.randint(1, 12)): + key = b"".join(rng.choice(alphabet[:-1]) for _ in range(rng.randint(1, 6))) + if rng.random() < 0.3: + val = b"" + else: + val = rng.choice(alphabet[:-1]) + b"".join( + rng.choice(alphabet) for _ in range(rng.randint(0, 8)) + ) + pairs.append((key, val)) + body = b"".join(k + delim + v.replace(delim, delim * 2) + delim for k, v in pairs) + got = _split(body, delim[0]) + assert got == [x for kv in pairs for x in kv], (delim, pairs, got) + + +def test_latin1_value_is_kept_not_replaced(): + # FCS 2.0/3.0 TEXT from Windows instruments is not always UTF-8. + text = b"/$PAR/1/$TOT/1/$P1N/FSC/$CYT/CELLQuest\xaa 3.3/" + h = decode_fcs(raw_fcs(text, version=b"FCS2.0")) + assert h.fields["cytometer"] == "CELLQuest\xaa 3.3" + assert h.fields["cytometer"].encode("latin-1") == b"CELLQuest\xaa 3.3" + + +def test_missing_par_still_lists_channels_from_pnn(): + text = b"/$TOT/4/$P1N/FSC/$P2N/SSC/$P3N/FL1/$DATATYPE/F/" + h = decode_fcs(raw_fcs(text)) + assert h.fields["numParameters"] is None + assert [c["name"] for c in h.fields["channels"]] == ["FSC", "SSC", "FL1"] + assert "? parameters" in h.summary + + +def test_channel_without_pnn_falls_back_to_pns(): + text = b"/$PAR/2/$TOT/1/$P1S/Forward Scatter/$P2N/SSC/$P2S/Side/" + h = decode_fcs(raw_fcs(text)) + assert h.fields["channels"] == [{"name": "Forward Scatter"}, {"name": "SSC", "label": "Side"}] + + +def _with_stext(stext, leading_delim=True, offset_from_text_end=0, begin_override=None, end_override=None): + """A file whose primary TEXT points at a supplemental TEXT placed after it.""" + primary = b"/$PAR/2/$TOT/1/$P1N/FSC/$P2N/SSC/$BEGINSTEXT/%08d/$ENDSTEXT/%08d/" + length = len(primary % (0, 0)) + begin = 58 + length + offset_from_text_end + seg = (b"/" if leading_delim else b"") + stext + end = begin + len(seg) - 1 + text = primary % (begin_override if begin_override is not None else begin, + end_override if end_override is not None else end) + return raw_fcs(text, tail=b"\x00" * offset_from_text_end + seg) + + +@pytest.mark.parametrize("leading_delim", [True, False]) +def test_supplemental_text_is_read_and_primary_wins(leading_delim): + # Optional keywords may live in the supplemental TEXT. $PAR is given twice; + # the primary value (2) wins over the supplemental one (99). + h = decode_fcs(_with_stext(b"$CYT/Aurora/$P1S/Forward/$PAR/99/", leading_delim=leading_delim, offset_from_text_end=16)) + assert h.fields["cytometer"] == "Aurora" + assert h.fields["channels"][0] == {"name": "FSC", "label": "Forward"} + assert h.fields["numParameters"] == 2 + + +def test_supplemental_text_outside_the_read_is_ignored(): + blob = _with_stext(b"$CYT/Aurora/", offset_from_text_end=16) + cut = blob[:blob.index(b"\x00" * 16)] # everything up to the STEXT + h = decode_fcs(cut) + assert h.fields["numParameters"] == 2 + assert "cytometer" not in h.fields + + +def test_supplemental_text_with_bad_offsets_is_ignored(): + h = decode_fcs(_with_stext(b"$CYT/Aurora/", begin_override=500, end_override=400)) + assert h is not None and "cytometer" not in h.fields + h = decode_fcs(_with_stext(b"$CYT/Aurora/", begin_override=99999990, end_override=99999999)) + assert h is not None and "cytometer" not in h.fields + + +def test_next_dataset_sets_multiple_datasets(): + # A Coulter .lmd carries an FCS 2.0 and an FCS 3.0 dataset; $NEXTDATA points + # at the second. A single-dataset file has no flag at all. + assert decode_fcs(make_fcs(extra=[("$NEXTDATA", "363602")])).fields["multipleDatasets"] is True + assert "multipleDatasets" not in decode_fcs(make_fcs()).fields + + +def test_keyword_values_are_normalized(): + h = decode_fcs(make_fcs(extra=[("$DATATYPE", " f "), ("$CYT", " Attune NxT "), ("$TOT", " 42 ")])) + assert h.fields["dataType"] == "float32" + assert h.fields["cytometer"] == "Attune NxT" + assert h.fields["numEvents"] == 42 + + +def test_header_offset_with_trailing_byte(): + # Some LMD exports leave a stray byte after the digits of an offset field. + blob = bytearray(make_fcs()) + assert blob[10:18] == b" 58" + blob[10:18] = b" 58\x00" + h = decode_fcs(bytes(blob)) + assert h is not None and h.fields["numParameters"] == 4 + + +def test_leading_delimiter_run_in_body_is_skipped(): + # A TEXT that opens with two delimiters (an empty first keyword) still + # yields the pairs after it. + h = decode_fcs(raw_fcs(b"//$PAR/1/$TOT/9/$P1N/FSC/")) + assert h.fields["numEvents"] == 9 + assert h.fields["channels"] == [{"name": "FSC"}]