From a4b96536f6b7cb8a328027a5eebc8dd7e7179d3a Mon Sep 17 00:00:00 2001 From: aaronkanzer <36093535+aaronkanzer@users.noreply.github.com> Date: Sat, 24 Oct 2026 19:56:34 -0400 Subject: [PATCH] feat: decode A3M multiple sequence alignments (0.8.0) The MSA format ColabFold and the MMseqs2 server write and AlphaFold, Boltz, Chai, Protenix and OpenFold read. Reports the query name and length, the depth within a 4 MiB read with truncated when the file goes on, ColabFold's # chain lengths and cardinalities, mean query coverage, rows with insertions, rows whose aligned length does not fit the query, HH-suite annotation rows, the hits' source databases and the first hit's header. The server's NUL-separated multi-query files list every query. Agrees with an independent parse on 188 real files from ColabFold, Protenix, Boltz, OpenFold3 and jackhmmer; 3000 mutants never raise. --- CHANGELOG.md | 13 ++ README.md | 14 +- pyproject.toml | 2 +- src/scigantic_headers/__init__.py | 3 +- src/scigantic_headers/a3m.py | 212 ++++++++++++++++++++++++++++++ tests/test_a3m.py | 165 +++++++++++++++++++++++ 6 files changed, 406 insertions(+), 3 deletions(-) create mode 100644 src/scigantic_headers/a3m.py create mode 100644 tests/test_a3m.py diff --git a/CHANGELOG.md b/CHANGELOG.md index e022a20..e70951e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,18 @@ # Changelog +## 0.8.0 + +- A3M: a decoder for multiple sequence alignments (`.a3m`, `.a3m.gz`), the + MSA format ColabFold and the MMseqs2 server write and AlphaFold, Boltz, + Chai, Protenix and OpenFold read. Reports the query name and length, the + depth within a 4 MiB read with `truncated` when the file goes on, ColabFold's + `#` chain lengths and cardinalities, mean query coverage, rows with + insertions, rows whose aligned length does not fit the query, HH-suite + annotation rows, the source databases of the hits and the first hit's + header. The server's NUL-separated multi-query files list every query. + Agrees with an independent parse on 188 real files from ColabFold, + Protenix, Boltz, OpenFold3 and jackhmmer. + ## 0.7.2 - FCS: the TEXT splitter now tells a legacy empty value from an escaped diff --git a/README.md b/README.md index 2f20401..49c9b2b 100644 --- a/README.md +++ b/README.md @@ -83,6 +83,17 @@ Decoders today: - **mmCIF / CIF**: structures (the modern PDB format). Entry id, title, keywords, experiment method, deposition date, and resolution (including the cryo-EM resolution item). +- **A3M**: multiple sequence alignments, the MSA format ColabFold and the + MMseqs2 server write and AlphaFold, Boltz, Chai, Protenix and OpenFold read. + Query name and length, depth (sequences in the read, `truncated` when the + file goes on past it), ColabFold's `#` chain lengths and copy counts, the + mean fraction of the query each row covers, rows with insertions, rows whose + aligned length does not fit the query (an un-realigned or mis-staged MSA), + HH-suite `ss_pred`/`ss_conf` annotation rows, the databases the hits came + from (UniRef, UniProtKB, UniParc, MGnify, SRA) and the first hit's header. + The server's NUL-separated multi-query files report every query. Checked + against an independent parse on 188 ColabFold, Protenix, Boltz, OpenFold3 + and jackhmmer files. - **GenBank**: annotated sequence records. LOCUS (name, length, molecule type, topology), definition, accession, version, and organism. - **GFF3 / GTF**: genome annotation. Version and a preview of the sources and @@ -275,6 +286,7 @@ detector-pixel-size / magnification pair are all handled. sam.py SAM alignment-header decode pdb.py PDB structure-header decode mmcif.py mmCIF / CIF structure-header decode + a3m.py A3M multiple-sequence-alignment decode genbank.py GenBank flat-file header decode gff.py GFF3 / GTF annotation-header decode bed.py BED interval-file decode @@ -297,7 +309,7 @@ detector-pixel-size / magnification pair are all handled. optics.py read_session_optics: find a data file's optics file cli.py scigantic-headers benchmark.py scigantic-headers-bench, measures the speed levers - tests/ pytest (454 tests, incl. fuzz, golden and sparse real-file fixtures) + tests/ pytest (507 tests, incl. fuzz, golden and sparse real-file fixtures) ## Robustness diff --git a/pyproject.toml b/pyproject.toml index 6a541b7..3d82b54 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "scigantic-headers" -version = "0.7.2" +version = "0.8.0" 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/__init__.py b/src/scigantic_headers/__init__.py index 4f1417a..f8f6cb4 100644 --- a/src/scigantic_headers/__init__.py +++ b/src/scigantic_headers/__init__.py @@ -11,6 +11,7 @@ from importlib.metadata import PackageNotFoundError, version as _pkg_version +from .a3m import decode_a3m from .acquisition import Finding, compare_acquisition from .bed import decode_bed from .czi import decode_czi @@ -85,7 +86,7 @@ "decode_nifti_header", "decode_npy_header", "decode_cryosparc_header", "parse_cryosparc_optics", "read_cryosparc_optics", "read_session_optics", "decode_parquet", "decode_fcs", "decode_mzml", - "decode_fastq", "decode_illumina_run", "decode_vcf", "decode_sam", "decode_pdb", + "decode_a3m", "decode_fastq", "decode_illumina_run", "decode_vcf", "decode_sam", "decode_pdb", "decode_mmcif", "decode_genbank", "decode_gff", "decode_bed", "decode_dicom", "decode_czi", "decode_hdf5", "decode_ics", "decode_klb", "decode_lif", "decode_nd2", "decode_obf", "decode_tiff", "parse_ome_xml", diff --git a/src/scigantic_headers/a3m.py b/src/scigantic_headers/a3m.py new file mode 100644 index 0000000..e811916 --- /dev/null +++ b/src/scigantic_headers/a3m.py @@ -0,0 +1,212 @@ +"""Decode an A3M multiple sequence alignment. + +A3M is the FASTA-like alignment format HH-suite defined and the MSA format the +structure predictors read: ColabFold and the MMseqs2 server write it; AlphaFold, +Boltz, Chai, Protenix and OpenFold take it as input. The first record is the +query. In every later record an upper-case letter or `-` is a match state +aligned to a query column and a lower-case letter is an insertion relative to +the query, which carries no column; dropping the lower-case letters leaves +every row exactly the query's length. + +Writers differ in the trimmings, and this reads each of them: + +- ColabFold's complex MSA starts with `#L1,L2C1,C2`: the chain lengths + and copy counts of the concatenated query (`#2311` for a monomer). +- The MMseqs2 / ColabFold server answers several queries in one file, blocks + separated by a NUL byte, each with its own `>101`, `>102` query record. +- HH-suite puts `>ss_pred` / `>ss_conf` annotation rows before the query. +- A sequence may be wrapped over several lines (reformat.pl) or be one line. + +Depth is the number that matters (prediction quality tracks it), so the record +count is reported for what the read covered, with `truncated` saying whether +the file went on past the read, together with the mean fraction of the query +each row covers and the rows whose aligned length does not match the query +(an MSA built for a different sequence, or a mis-staged file). +""" + +from __future__ import annotations + +import re +from typing import Dict, List, Optional, Tuple + +from .decoders import DecodedHeader, Read, register_decoder + +# A ColabFold uniref MSA for a 300-residue query is 2 to 3 MB; the environmental +# one can pass 10 MB. Counting past the read is not possible, so `truncated` +# says when the count is a lower bound. +A3M_READ_BYTES = 4 * 1024 * 1024 + +# HH-suite annotation rows: not sequences, excluded from the depth, listed. +_ANNOTATION = {b"ss_pred", b"ss_conf", b"ss_dssp", b"aa_dssp", b"aa_pred", b"sa_dssp"} +_LOWER = bytes(range(ord("a"), ord("z") + 1)) +_INSERT = _LOWER + b"." # a `.` is an insertion-column gap in a2m-style rows +_SEQ_CHARS = frozenset(b"ABCDEFGHIJKLMNOPQRSTUVWXYZ" + _LOWER + b"-.*") +_QUERIES_MAX = 32 +_INT_LIST = re.compile(rb"^\d+(,\d+)*$") +_HEADER_PREVIEW = 120 + +# Where a hit came from, read off its name. ColabFold's uniref.a3m is +# UniRef100_; its environmental MSA mixes MGnify proteins and SRA assemblies; +# jackhmmer / HHblits outputs carry UniProt `tr|` / `sp|` and UniRef90 names. +_DATABASES: List[Tuple[re.Pattern, str]] = [ + (re.compile(rb"^UniRef(100|90|50)_"), "UniRef%s"), + (re.compile(rb"^(?:tr|sp)\|"), "UniProtKB"), + (re.compile(rb"^UPI[0-9A-F]{10}"), "UniParc"), + (re.compile(rb"^MGY[PS]\d"), "MGnify"), + (re.compile(rb"^[SED]RR\d+"), "SRA"), +] + + +def _chain_spec(line: bytes) -> Tuple[Optional[List[int]], Optional[List[int]]]: + """ColabFold's `#101,1301,1` line: chain lengths and cardinalities, or + (None, None) when the comment is something else.""" + parts = line[1:].strip().split() + if len(parts) != 2 or not _INT_LIST.match(parts[0]) or not _INT_LIST.match(parts[1]): + return None, None + lengths = [int(x) for x in parts[0].split(b",")] + cards = [int(x) for x in parts[1].split(b",")] + if len(lengths) != len(cards) or min(lengths) <= 0 or min(cards) <= 0 or max(lengths + cards) >= 10 ** 9: + return None, None + return lengths, cards + + +def _records(block: bytes): + """(header, sequence) pairs of one block, plus its leading `#` comment + lines; None if the block is not alignment text.""" + comments: List[bytes] = [] + records: List[Tuple[bytes, bytearray]] = [] + for line in block.split(b"\n"): + line = line.rstrip(b" \t") + if not line: + continue + if line[:1] == b">": + records.append((line[1:], bytearray())) + elif not records: + if line[:1] == b"#": + comments.append(line) + continue + return None # sequence text before any header: not A3M + else: + records[-1][1].extend(line) + return comments, records + + +def _has_insertions(seq: bytes) -> bool: + return len(seq.translate(None, _LOWER)) != len(seq) + + +def _database(name: bytes) -> Optional[str]: + for pattern, tag in _DATABASES: + m = pattern.match(name) + if m: + return tag % m.group(1).decode("ascii") if "%s" in tag else tag + return None + + +def decode_a3m(data: bytes) -> Optional[DecodedHeader]: + if not data: + return None + truncated = len(data) >= A3M_READ_BYTES + blocks = [b for b in data.replace(b"\r", b"").split(b"\x00") if b.strip()] + if not blocks: + return None + + parsed = [] + for i, block in enumerate(blocks): + r = _records(block) + if r is None: + return None + comments, records = r + if truncated and i == len(blocks) - 1 and records: + records.pop() # the read ended inside this record + parsed.append((comments, records)) + + chain_lengths = chain_cards = None + for line in parsed[0][0]: + chain_lengths, chain_cards = _chain_spec(line) + if chain_lengths: + break + + queries = [] + first = None + for comments, records in parsed: + annotations: List[str] = [] + query = None + hits = [] + for header, seq in records: + name = header.split(None, 1)[0] if header.strip() else b"" + if name in _ANNOTATION: + annotations.append(name.decode("ascii")) + continue + if any(c not in _SEQ_CHARS for c in seq): + return None + if query is None: + query = (name, header, bytes(seq)) + else: + hits.append((name, header, bytes(seq))) + if query is None: + continue + query_len = len(query[2].translate(None, _INSERT)) + if query_len == 0: + continue + if first is None: + first = (query, hits, annotations, query_len) + queries.append({ + "name": query[0].decode("utf-8", "replace") or None, + "length": query_len, + "sequenceCount": 1 + len(hits), + }) + if first is None: + return None + + query, hits, annotations, query_len = first + with_insertions = 1 if _has_insertions(query[2]) else 0 + mismatches = 0 + coverage = 0.0 + covered = 0 + databases: Dict[str, int] = {} + for name, header, seq in hits: + match = seq.translate(None, _INSERT) + if _has_insertions(seq): + with_insertions += 1 + if len(match) != query_len: + mismatches += 1 + else: + coverage += 1.0 - match.count(b"-") / query_len + covered += 1 + tag = _database(name) + if tag: + databases[tag] = databases.get(tag, 0) + 1 + + first_hit = hits[0][1].decode("utf-8", "replace")[:_HEADER_PREVIEW] if hits else None + fields = { + "queryName": queries[0]["name"], + "queryLength": query_len, + "chainLengths": chain_lengths, + "chainCardinalities": chain_cards, + "sequenceCount": 1 + len(hits), + "meanQueryCoverage": round(coverage / covered, 4) if covered else None, + "sequencesWithInsertions": with_insertions, + "lengthMismatches": mismatches, + "annotationRows": annotations or None, + "databases": databases or None, + "firstHitHeader": first_hit, + "queryCount": len(queries), + "queries": queries[:_QUERIES_MAX], + "truncated": truncated, + } + + chains = "" + if chain_lengths and sum(chain_cards) > 1: + expanded = [str(n) for n, c in zip(chain_lengths, chain_cards) for _ in range(c)] + chains = " (%d chains: %s)" % (len(expanded), "+".join(expanded)) + depth = "%d%s sequences" % (1 + len(hits), "+" if truncated else "") + head = "A3M MSA, " if len(queries) == 1 else "A3M MSA, %d queries; first: " % len(queries) + return DecodedHeader( + format="a3m", + summary="%squery %d aa%s, %s" % (head, query_len, chains, depth), + fields=fields, + ) + + +register_decoder("a3m", decode_a3m, read=Read(leading=A3M_READ_BYTES)) # depth needs the rows diff --git a/tests/test_a3m.py b/tests/test_a3m.py new file mode 100644 index 0000000..cfa4efb --- /dev/null +++ b/tests/test_a3m.py @@ -0,0 +1,165 @@ +"""A3M multiple-sequence-alignment decoding.""" + +import json + +import pytest + +from scigantic_headers import decode_a3m, decode_bytes, has_decoder_for, is_decodable, read_for +from scigantic_headers.a3m import A3M_READ_BYTES + +# A ColabFold complex MSA: the `#` chain line, numeric query name, MMseqs2 +# tab-separated hit headers, a lower-case insertion and a gapped row. +COLABFOLD = ( + b"#5,4\t1,1\n" + b">101\n" + b"MKTAYWGHE\n" + b">UniRef100_A0A001\t75\t0.590\t6.662E-14\t2\t9\t9\t1419\t1427\t1504\n" + b"MKSAYWGHE\n" + b">UniRef100_A0A002\t60\t0.400\t1.0E-10\t1\t9\t9\t1\t9\t20\n" + b"MKTAYagWG--\n" + b">MGYP000000000001\t50\t0.3\t1.0E-5\t1\t5\t9\t1\t5\t5\n" + b"MKTAY----\n" +) + + +def test_colabfold_complex(): + h = decode_a3m(COLABFOLD) + assert h.format == "a3m" + f = h.fields + assert f["queryName"] == "101" + assert f["queryLength"] == 9 + assert f["chainLengths"] == [5, 4] + assert f["chainCardinalities"] == [1, 1] + assert f["sequenceCount"] == 4 # the query counts, as in every depth plot + assert f["sequencesWithInsertions"] == 1 + assert f["lengthMismatches"] == 0 + # coverage: 9/9, 7/9, 5/9 -> mean 0.7778 + assert f["meanQueryCoverage"] == pytest.approx(0.7778) + assert f["databases"] == {"UniRef100": 2, "MGnify": 1} + assert f["firstHitHeader"].startswith("UniRef100_A0A001\t75\t0.590") + assert f["queryCount"] == 1 + assert f["queries"] == [{"name": "101", "length": 9, "sequenceCount": 4}] + assert f["truncated"] is False + assert h.summary == "A3M MSA, query 9 aa (2 chains: 5+4), 4 sequences" + + +def test_homodimer_chain_line_expands_cardinality(): + h = decode_a3m(b"#9\t2\n>101\nMKTAYWGHE\n") + assert h.fields["chainLengths"] == [9] + assert h.fields["chainCardinalities"] == [2] + assert "(2 chains: 9+9)" in h.summary + + +def test_monomer_chain_line_is_not_announced(): + h = decode_a3m(b"#9\t1\n>101\nMKTAYWGHE\n") + assert h.fields["chainLengths"] == [9] + assert h.summary == "A3M MSA, query 9 aa, 1 sequences" + + +def test_server_multi_query_nul_blocks(): + # The MMseqs2 / ColabFold server answers two queries in one file, blocks + # separated by NUL, the file ending in one. + data = ( + b">101\nMKTAY\n>UniRef100_X\t1\t1\t1\t1\t5\t5\t1\t5\t5\nMKSAY\n\x00" + b">102\nWGHEQLN\n>UniRef100_Y\t1\t1\t1\t1\t7\t7\t1\t7\t7\nWGHE---\n" + b">UniRef100_Z\t1\t1\t1\t1\t7\t7\t1\t7\t7\nWGHEQLN\n\x00" + ) + h = decode_a3m(data) + f = h.fields + assert f["queryCount"] == 2 + assert f["queries"] == [ + {"name": "101", "length": 5, "sequenceCount": 2}, + {"name": "102", "length": 7, "sequenceCount": 3}, + ] + assert f["queryName"] == "101" and f["sequenceCount"] == 2 # top-level fields describe the first block + assert h.summary == "A3M MSA, 2 queries; first: query 5 aa, 2 sequences" + + +def test_hhsuite_annotation_rows_wrapped_lines_and_crlf(): + data = ( + b">ss_pred\r\nCCHHHHH\r\n" + b">ss_conf\r\n9876543\r\n" + b">d1abca_ b.1.1.1 (A:) some domain\r\n" + b"MKTA\r\nYWG\r\n" + b">tr|Q9XYZ1|Q9XYZ1_HUMAN Protein OS=Homo sapiens\r\n" + b"MKsaTAY\r\nWG\r\n" + ) + h = decode_a3m(data) + f = h.fields + assert f["annotationRows"] == ["ss_pred", "ss_conf"] + assert f["queryName"] == "d1abca_" + assert f["queryLength"] == 7 + assert f["sequenceCount"] == 2 # annotation rows are not sequences + assert f["sequencesWithInsertions"] == 1 + assert f["lengthMismatches"] == 0 + assert f["meanQueryCoverage"] == 1.0 + assert f["databases"] == {"UniProtKB": 1} + + +def test_length_mismatch_is_counted_not_hidden(): + # An MSA whose rows do not fit the query: the realignment step was skipped, + # or the file belongs to another sequence. Mismatched rows are counted and + # left out of the coverage mean. + data = b">query\nMKTAYWGHE\n>hit1\nMKTAYWGHE\n>hit2\nMKTAY\n>hit3\nMKTAYWGHEQQ\n" + f = decode_a3m(data).fields + assert f["sequenceCount"] == 4 + assert f["lengthMismatches"] == 2 + assert f["meanQueryCoverage"] == 1.0 + + +def test_query_only_has_no_coverage(): + f = decode_a3m(b">query \nGSHMTGKSVDMY\n").fields + assert f["queryName"] == "query" + assert f["sequenceCount"] == 1 + assert f["meanQueryCoverage"] is None + assert f["firstHitHeader"] is None + assert f["databases"] is None + + +def test_dot_is_an_insertion_gap_and_star_a_match_state(): + f = decode_a3m(b">q\nMKTAY*\n>h\nMK..TAY*\n").fields + assert f["queryLength"] == 6 + assert f["lengthMismatches"] == 0 + assert f["sequencesWithInsertions"] == 0 + + +def test_truncated_read_drops_the_cut_record_and_says_so(): + rows = b">q\n" + b"A" * 60 + b"\n" + hit = b">UniRef100_H\n" + b"A" * 60 + b"\n" + n = (A3M_READ_BYTES - len(rows)) // len(hit) + 2 + data = (rows + hit * n)[:A3M_READ_BYTES] + assert len(data) == A3M_READ_BYTES + f = decode_a3m(data).fields + assert f["truncated"] is True + # Every complete record counts; the one the read cut through does not. + whole = data.count(b"\n>") + assert f["sequenceCount"] in (whole, whole + 1) + assert f["lengthMismatches"] == 0 + assert decode_a3m(data).summary.endswith("+ sequences") + + +def test_rejects_non_alignment_text(): + assert decode_a3m(b"") is None + assert decode_a3m(b"MKTAYWGHE\n>q\nMKTAY\n") is None # sequence before any header + assert decode_a3m(b">q\nMKT AY WGHE\n") is None # not residue letters + assert decode_a3m(b">q\n123456\n") is None + assert decode_a3m(b"# just a comment\n") is None + assert decode_a3m(b">\n\n") is None # header with no residues + assert decode_a3m(b"\x00\x00") is None + + +def test_odd_chain_lines_are_ignored_not_fatal(): + assert decode_a3m(b"#hello world\n>q\nMKT\n").fields["chainLengths"] is None + assert decode_a3m(b"#1,2\t1\n>q\nMKT\n").fields["chainLengths"] is None # unequal lists + assert decode_a3m(b"#0\t1\n>q\nMKT\n").fields["chainLengths"] is None + + +def test_dispatch_and_read_size(): + assert has_decoder_for("msa/uniref.a3m") + assert is_decodable("msa/uniref.a3m.gz") # the reader inflates .gz before dispatch + assert decode_bytes("x.a3m", COLABFOLD).format == "a3m" + assert read_for("x.a3m").leading == A3M_READ_BYTES + + +def test_json_safe(): + json.dumps(decode_a3m(COLABFOLD).to_dict(), allow_nan=False)