Decode A3M multiple sequence alignments (0.8.0)

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

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 <file|url|--dir>
       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
 

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"

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",

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,L2<TAB>C1,C2`: the chain lengths
+  and copy counts of the concatenated query (`#231<TAB>1` 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,130<TAB>1,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