Python API

CRAIC is a library as well as a workbench. Everything the GUI does is available from Python, which is the supported way to script it or to build on it.

from craic import io, progressive, evaluate
from craic.ambiguity import reliability
from craic.domain import Alphabet

records = io.read_records("seqs.fasta")
aln = progressive.align(records, Alphabet.PROTEIN)

report = reliability.analyse(aln, do_perturbation=True)
print(report.col_combined)                  # per-column reliability, 0..1

truth = io.load_alignment("truth.fasta", alphabet=Alphabet.PROTEIN)
acc = evaluate.compare_to_reference(aln, truth)
print(acc.sp, acc.tc, acc.reliability_auc(report.col_combined))

Nothing below imports Qt, so all of it works headless.

---

The canonical model

craic.domain

Core domain model.

A single canonical representation is stored at the sequences' native level (nucleotide rows for DNA/RNA, amino-acid rows for protein). Codon and amino-acid views are derived projections of a coding nucleotide alignment, so the GUI can switch levels live without ever holding three copies of the data.

Alphabet

Bases: Enum

Source code in craic/domain.py

class Alphabet(Enum):
    DNA = "dna"
    RNA = "rna"
    PROTEIN = "protein"

    @property
    def is_nucleotide(self) -> bool:
        return self in (Alphabet.DNA, Alphabet.RNA)

Level

Bases: Enum

Source code in craic/domain.py

class Level(Enum):
    NT = "nt"
    CODON = "codon"
    AA = "aa"

Alignment dataclass

Source code in craic/domain.py

@dataclass
class Alignment:
    ids: List[str]
    rows: List[str]
    alphabet: Alphabet
    coding: Optional[CodingSpec] = None
    meta: Dict[str, object] = field(default_factory=dict)

    def __post_init__(self) -> None:
        if not self.rows:
            return
        # Canonicalise gap characters once, here, so that every residue-index
        # walk downstream can test ``ch != "-"``. See GAP_CHARS.
        self.rows = [normalise_gaps(r) for r in self.rows]
        width = len(self.rows[0])
        if any(len(r) != width for r in self.rows):
            raise ValueError("alignment rows are not equal length")
        if len(self.ids) != len(self.rows):
            raise ValueError("ids and rows count mismatch")

    # -- basic shape ------------------------------------------------------- #
    @property
    def n_seqs(self) -> int:
        return len(self.rows)

    @property
    def length(self) -> int:
        return len(self.rows[0]) if self.rows else 0

    def column(self, j: int) -> str:
        return "".join(r[j] for r in self.rows)

    def gap_fraction(self, j: int) -> float:
        col = self.column(j)
        return sum(c in _GAP for c in col) / len(col)

    # -- levels ------------------------------------------------------------ #
    def available_levels(self) -> List[Level]:
        if self.alphabet == Alphabet.PROTEIN:
            return [Level.AA]
        if self.coding is not None:
            return [Level.NT, Level.CODON, Level.AA]
        return [Level.NT]

    def n_codons(self) -> int:
        if not self.coding:
            return 0
        return max(0, (self.length - self.coding.frame) // 3)

    def _codon_cols(self, c: int) -> Tuple[int, int, int]:
        start = self.coding.frame + 3 * c
        return start, start + 1, start + 2

    def is_codon_aware(self) -> bool:
        """True if no codon contains a partial gap (frame is intact)."""
        if not self.coding:
            return False
        for r in self.rows:
            for c in range(self.n_codons()):
                a, b, d = self._codon_cols(c)
                cod = r[a] + r[b] + r[d]
                gaps = cod.count("-")
                if gaps not in (0, 3):
                    return False
        return True

    def frame_break_columns(self) -> List[int]:
        """Codon columns where at least one row has a partial-gap codon."""
        if not self.coding:
            return []
        out = []
        for c in range(self.n_codons()):
            a, b, d = self._codon_cols(c)
            broken = any(
                0 < (r[a] + r[b] + r[d]).count("-") < 3 for r in self.rows
            )
            if broken:
                out.append(c)
        return out

    def internal_stops(self) -> List[Tuple[int, int]]:
        """``(row, codon)`` for every stop codon that is not the last codon of
        its sequence.

        A stop in the middle of a coding sequence usually means one of three
        things, and the user needs to be told which is being assumed: the wrong
        genetic code (mitochondrial and Mollicute genes read TGA as tryptophan,
        the standard code reads it as a stop), the wrong reading frame, or a
        pseudogene. One check catches all three, which is why it is worth making
        even though CRAIC cannot tell them apart.

        Trailing gaps are ignored, so a short sequence padded out to the
        alignment width is not reported as stopping early.
        """
        if not self.coding:
            return []
        out: List[Tuple[int, int]] = []
        for ri, row in enumerate(self.rows):
            last = -1
            for c in range(self.n_codons()):
                a, b, d = self._codon_cols(c)
                if (row[a] + row[b] + row[d]) != "---":
                    last = c
            for c in range(self.n_codons()):
                if c >= last:
                    break                      # the final codon may legitimately stop
                a, b, d = self._codon_cols(c)
                if translate_codon(row[a] + row[b] + row[d], self.coding.table) == "*":
                    out.append((ri, c))
        return out

    def display(self, level: Level) -> "DisplayMatrix":
        """Project to a renderable matrix of string cells.

        Returns ``cells[row][col]`` plus, for each display column, the starting
        nucleotide index so the GUI can map selections back to nt coordinates.
        """
        if level == Level.NT or self.alphabet == Alphabet.PROTEIN:
            cells = [list(r) for r in self.rows]
            nt_start = list(range(self.length))
            return DisplayMatrix(cells, nt_start, span=1)

        if not self.coding:
            raise ValueError("codon/aa view requires a CodingSpec")

        nc = self.n_codons()
        nt_start = [self.coding.frame + 3 * c for c in range(nc)]
        cells: List[List[str]] = []
        for r in self.rows:
            row_cells = []
            for c in range(nc):
                a, b, d = self._codon_cols(c)
                cod = r[a] + r[b] + r[d]
                if level == Level.CODON:
                    row_cells.append(cod)
                else:  # AA
                    row_cells.append(translate_codon(cod, self.coding.table))
            cells.append(row_cells)
        return DisplayMatrix(cells, nt_start, span=3)

    def translated_rows(self) -> List[str]:
        dm = self.display(Level.AA)
        return ["".join(c) for c in dm.cells]

    def amino_acid_alignment(self) -> "Alignment":
        """A protein Alignment of the translated rows, for amino-acid-level tools.

        For protein data this is the alignment itself; for a coding nucleotide
        alignment each amino-acid column corresponds to one codon column.
        """
        if self.alphabet == Alphabet.PROTEIN:
            return self
        if not self.coding:
            raise ValueError("amino-acid view requires a CodingSpec")
        return Alignment(list(self.ids), self.translated_rows(), Alphabet.PROTEIN)

    def codon_col(self, nt_col: int) -> int:
        """Map a nucleotide column to its codon / amino-acid column index."""
        if not self.coding:
            return nt_col
        return (nt_col - self.coding.frame) // 3

    # -- io helpers -------------------------------------------------------- #
    def slice_columns(self, start: int, stop: int) -> "Alignment":
        return Alignment(
            list(self.ids),
            [r[start:stop] for r in self.rows],
            self.alphabet,
            coding=None,
            meta=dict(self.meta),
        )

    @classmethod
    def from_records(
        cls,
        records: Sequence[Tuple[str, str]],
        alphabet: Optional[Alphabet] = None,
        coding: Optional[CodingSpec] = None,
    ) -> "Alignment":
        ids = [r[0] for r in records]
        rows = [r[1] for r in records]
        alph = alphabet or detect_alphabet(rows)
        return cls(ids, rows, alph, coding=coding)

is_codon_aware

is_codon_aware() -> bool

True if no codon contains a partial gap (frame is intact).

Source code in craic/domain.py

def is_codon_aware(self) -> bool:
    """True if no codon contains a partial gap (frame is intact)."""
    if not self.coding:
        return False
    for r in self.rows:
        for c in range(self.n_codons()):
            a, b, d = self._codon_cols(c)
            cod = r[a] + r[b] + r[d]
            gaps = cod.count("-")
            if gaps not in (0, 3):
                return False
    return True

frame_break_columns

frame_break_columns() -> List[int]

Codon columns where at least one row has a partial-gap codon.

Source code in craic/domain.py

def frame_break_columns(self) -> List[int]:
    """Codon columns where at least one row has a partial-gap codon."""
    if not self.coding:
        return []
    out = []
    for c in range(self.n_codons()):
        a, b, d = self._codon_cols(c)
        broken = any(
            0 < (r[a] + r[b] + r[d]).count("-") < 3 for r in self.rows
        )
        if broken:
            out.append(c)
    return out

internal_stops

internal_stops() -> List[Tuple[int, int]]

(row, codon) for every stop codon that is not the last codon of its sequence.

A stop in the middle of a coding sequence usually means one of three things, and the user needs to be told which is being assumed: the wrong genetic code (mitochondrial and Mollicute genes read TGA as tryptophan, the standard code reads it as a stop), the wrong reading frame, or a pseudogene. One check catches all three, which is why it is worth making even though CRAIC cannot tell them apart.

Trailing gaps are ignored, so a short sequence padded out to the alignment width is not reported as stopping early.

Source code in craic/domain.py

def internal_stops(self) -> List[Tuple[int, int]]:
    """``(row, codon)`` for every stop codon that is not the last codon of
    its sequence.

    A stop in the middle of a coding sequence usually means one of three
    things, and the user needs to be told which is being assumed: the wrong
    genetic code (mitochondrial and Mollicute genes read TGA as tryptophan,
    the standard code reads it as a stop), the wrong reading frame, or a
    pseudogene. One check catches all three, which is why it is worth making
    even though CRAIC cannot tell them apart.

    Trailing gaps are ignored, so a short sequence padded out to the
    alignment width is not reported as stopping early.
    """
    if not self.coding:
        return []
    out: List[Tuple[int, int]] = []
    for ri, row in enumerate(self.rows):
        last = -1
        for c in range(self.n_codons()):
            a, b, d = self._codon_cols(c)
            if (row[a] + row[b] + row[d]) != "---":
                last = c
        for c in range(self.n_codons()):
            if c >= last:
                break                      # the final codon may legitimately stop
            a, b, d = self._codon_cols(c)
            if translate_codon(row[a] + row[b] + row[d], self.coding.table) == "*":
                out.append((ri, c))
    return out

display

display(level: Level) -> 'DisplayMatrix'

Project to a renderable matrix of string cells.

Returns cells[row][col] plus, for each display column, the starting nucleotide index so the GUI can map selections back to nt coordinates.

Source code in craic/domain.py

def display(self, level: Level) -> "DisplayMatrix":
    """Project to a renderable matrix of string cells.

    Returns ``cells[row][col]`` plus, for each display column, the starting
    nucleotide index so the GUI can map selections back to nt coordinates.
    """
    if level == Level.NT or self.alphabet == Alphabet.PROTEIN:
        cells = [list(r) for r in self.rows]
        nt_start = list(range(self.length))
        return DisplayMatrix(cells, nt_start, span=1)

    if not self.coding:
        raise ValueError("codon/aa view requires a CodingSpec")

    nc = self.n_codons()
    nt_start = [self.coding.frame + 3 * c for c in range(nc)]
    cells: List[List[str]] = []
    for r in self.rows:
        row_cells = []
        for c in range(nc):
            a, b, d = self._codon_cols(c)
            cod = r[a] + r[b] + r[d]
            if level == Level.CODON:
                row_cells.append(cod)
            else:  # AA
                row_cells.append(translate_codon(cod, self.coding.table))
        cells.append(row_cells)
    return DisplayMatrix(cells, nt_start, span=3)

amino_acid_alignment

amino_acid_alignment() -> 'Alignment'

A protein Alignment of the translated rows, for amino-acid-level tools.

For protein data this is the alignment itself; for a coding nucleotide alignment each amino-acid column corresponds to one codon column.

Source code in craic/domain.py

def amino_acid_alignment(self) -> "Alignment":
    """A protein Alignment of the translated rows, for amino-acid-level tools.

    For protein data this is the alignment itself; for a coding nucleotide
    alignment each amino-acid column corresponds to one codon column.
    """
    if self.alphabet == Alphabet.PROTEIN:
        return self
    if not self.coding:
        raise ValueError("amino-acid view requires a CodingSpec")
    return Alignment(list(self.ids), self.translated_rows(), Alphabet.PROTEIN)

codon_col

codon_col(nt_col: int) -> int

Map a nucleotide column to its codon / amino-acid column index.

Source code in craic/domain.py

def codon_col(self, nt_col: int) -> int:
    """Map a nucleotide column to its codon / amino-acid column index."""
    if not self.coding:
        return nt_col
    return (nt_col - self.coding.frame) // 3

CodingSpec dataclass

How a nucleotide alignment maps to codons.

Source code in craic/domain.py

@dataclass(frozen=True)
class CodingSpec:
    """How a nucleotide alignment maps to codons."""

    frame: int = 0   # leading alignment columns to skip before the first codon
    table: int = 1   # NCBI genetic-code table id

genetic_codes

genetic_codes() -> List[Tuple[int, str]]

Every NCBI translation table Biopython knows, as (id, name).

The list is not hard-coded here because it is not ours to curate: NCBI adds tables, and anything shipped as a fixed list goes stale silently. Table 1 is forced to the front so the standard code is the first thing offered.

Source code in craic/domain.py

def genetic_codes() -> List[Tuple[int, str]]:
    """Every NCBI translation table Biopython knows, as ``(id, name)``.

    The list is not hard-coded here because it is not ours to curate: NCBI adds
    tables, and anything shipped as a fixed list goes stale silently. Table 1 is
    forced to the front so the standard code is the first thing offered.
    """
    from Bio.Data import CodonTable

    codes = []
    for tid, tbl in sorted(CodonTable.unambiguous_dna_by_id.items()):
        name = next((n for n in tbl.names if n and not n.startswith("SGC")), f"table {tid}")
        codes.append((tid, name))
    return sorted(codes, key=lambda c: (c[0] != 1, c[0]))

code_name

code_name(table_id: int) -> str

The NCBI name for a table id, or a bare label if it is unknown.

Source code in craic/domain.py

def code_name(table_id: int) -> str:
    """The NCBI name for a table id, or a bare label if it is unknown."""
    for tid, name in genetic_codes():
        if tid == table_id:
            return name
    return f"table {table_id}"

translate_codon

translate_codon(codon: str, table_id: int = 1) -> str

Translate one (possibly gapped) codon to a single residue.

'---' -> '-' (clean gap), any partial gap -> 'X' (frameshift / indel ambiguity, surfaced rather than hidden), unknown -> 'X', stop -> '*'.

Source code in craic/domain.py

def translate_codon(codon: str, table_id: int = 1) -> str:
    """Translate one (possibly gapped) codon to a single residue.

    '---' -> '-' (clean gap), any partial gap -> 'X' (frameshift / indel
    ambiguity, surfaced rather than hidden), unknown -> 'X', stop -> '*'.
    """
    codon = codon.upper().replace("U", "T")
    if codon == "---":
        return "-"
    if "-" in codon:
        return "X"
    if table_id not in _TABLE_CACHE:
        _TABLE_CACHE[table_id] = _build_table(table_id)
    return _TABLE_CACHE[table_id].get(codon, "X")

residue_index

residue_index(row: str) -> List[int]

Per column, the index of its residue in the ungapped sequence, or -1.

Source code in craic/domain.py

def residue_index(row: str) -> List[int]:
    """Per column, the index of its residue in the ungapped sequence, or -1."""
    out: List[int] = []
    r = 0
    for ch in row:
        if ch == "-":
            out.append(-1)
        else:
            out.append(r)
            r += 1
    return out

column_index

column_index(row: str) -> Dict[int, int]

The inverse of :func:residue_index: residue index -> column.

Source code in craic/domain.py

def column_index(row: str) -> Dict[int, int]:
    """The inverse of :func:`residue_index`: residue index -> column."""
    return {r: c for c, r in enumerate(residue_index(row)) if r >= 0}

membership

membership(aln: 'Alignment', by: str = 'index') -> Dict[Tuple[Any, int], int]

(sequence, residue) -> column for every non-gap cell.

by="index" keys on the sequence's position, by="id" on its identifier. Both are needed: scores computed per row want the position, while anything comparing two alignments has to match on the name.

Source code in craic/domain.py

def membership(aln: "Alignment", by: str = "index") -> Dict[Tuple[Any, int], int]:
    """``(sequence, residue) -> column`` for every non-gap cell.

    ``by="index"`` keys on the sequence's position, ``by="id"`` on its
    identifier. Both are needed: scores computed per row want the position,
    while anything comparing two alignments has to match on the name.
    """
    if by not in ("index", "id"):
        raise ValueError(f"by must be 'index' or 'id', not {by!r}")
    keys = range(aln.n_seqs) if by == "index" else aln.ids
    mem: Dict[Tuple[Any, int], int] = {}
    for key, row in zip(keys, aln.rows):
        for c, r in enumerate(residue_index(row)):
            if r >= 0:
                mem[(key, r)] = c
    return mem

normalise_gaps

normalise_gaps(row: str) -> str

Map every accepted gap character onto the canonical -.

Source code in craic/domain.py

def normalise_gaps(row: str) -> str:
    """Map every accepted gap character onto the canonical ``-``."""
    return row.translate(_GAP_NORMALISE)

ungap

ungap(row: str) -> str

The residues of a row, with every gap character removed.

Source code in craic/domain.py

def ungap(row: str) -> str:
    """The residues of a row, with every gap character removed."""
    return normalise_gaps(row).replace("-", "")

match_ids

match_ids(reference: Alignment, other: Alignment) -> Dict[str, str]

Map each reference id onto the corresponding id in other.

Several aligners (MAFFT among them) truncate FASTA identifiers at the first whitespace. Matching on the raw string then misses every sequence, and because a miss increments the pair denominator without incrementing the numerator, every column silently scores zero agreement — a wrong result that looks exactly like a real one. So: try the ids as given, then try them truncated, and if neither matches, say so.

Source code in craic/domain.py

def match_ids(reference: Alignment, other: Alignment) -> Dict[str, str]:
    """Map each reference id onto the corresponding id in ``other``.

    Several aligners (MAFFT among them) truncate FASTA identifiers at the first
    whitespace. Matching on the raw string then misses every sequence, and
    because a miss increments the pair denominator without incrementing the
    numerator, *every column silently scores zero agreement* — a wrong result
    that looks exactly like a real one. So: try the ids as given, then try them
    truncated, and if neither matches, say so.
    """
    other_ids = set(other.ids)
    if set(reference.ids) <= other_ids:
        return {rid: rid for rid in reference.ids}

    by_short: Dict[str, List[str]] = {}
    for oid in other.ids:
        by_short.setdefault(_short(oid), []).append(oid)
    mapped = {}
    for rid in reference.ids:
        cands = by_short.get(_short(rid), [])
        if len(cands) == 1:
            mapped[rid] = cands[0]
    if len(mapped) == len(reference.ids):
        return mapped

    missing = [rid for rid in reference.ids if rid not in mapped]
    raise IdMismatch(
        f"{len(missing)} of {len(reference.ids)} sequence ids in the reference "
        f"have no unique counterpart in the other alignment "
        f"(first unmatched: {missing[0]!r}); the two alignments cannot be compared"
    )

---

Reading and writing

craic.io

Reading and writing alignments / unaligned record sets.

read_records

read_records(path: str, fmt: Optional[str] = None) -> List[Record]

Source code in craic/io.py

def read_records(path: str, fmt: Optional[str] = None) -> List[Record]:
    from Bio import SeqIO

    fmt = fmt or guess_format(path)
    return [(r.id, str(r.seq)) for r in SeqIO.parse(path, fmt)]

load_alignment

load_alignment(path: str, fmt: Optional[str] = None, alphabet: Optional[Alphabet] = None, coding: Optional[CodingSpec] = None) -> Alignment

Source code in craic/io.py

def load_alignment(
    path: str,
    fmt: Optional[str] = None,
    alphabet: Optional[Alphabet] = None,
    coding: Optional[CodingSpec] = None,
) -> Alignment:
    recs = read_records(path, fmt)
    if not recs:
        return Alignment([], [], alphabet or Alphabet.PROTEIN, coding=coding)
    alph = alphabet or detect_alphabet([s for _, s in recs])
    groups = load_groups(path)
    annotations = load_annotations(path)
    if looks_aligned(recs):
        aln = Alignment.from_records(recs, alphabet=alph, coding=coding)
        if groups:
            aln.meta["groups"] = groups
        if annotations:
            aln.meta["annotations"] = annotations
        return aln
    # Not aligned: keep as a degenerate single-column-less alignment is wrong;
    # callers that load unaligned data should run an engine. We still return a
    # padded alignment so the viewer can show the raw sequences.
    width = max(len(s) for _, s in recs)
    padded = [(i, s.ljust(width, "-")) for i, s in recs]
    aln = Alignment.from_records(padded, alphabet=alph, coding=coding)
    aln.meta["unaligned"] = True
    if groups:
        aln.meta["groups"] = groups
    if annotations:
        aln.meta["annotations"] = annotations
    return aln

write_fasta

write_fasta(path: str, ids: Sequence[str], rows: Sequence[str], wrap: int = 60) -> None

Source code in craic/io.py

def write_fasta(path: str, ids: Sequence[str], rows: Sequence[str], wrap: int = 60) -> None:
    with open(path, "w") as fh:
        for name, seq in zip(ids, rows):
            fh.write(f">{name}\n")
            if wrap and wrap > 0:
                for k in range(0, len(seq), wrap):
                    fh.write(seq[k : k + wrap] + "\n")
            else:
                fh.write(seq + "\n")

---

The acceleration layer

The single primitive everything else is built on. The Rust core and the NumPy mirror satisfy the same contract and are cross-validated to ~1e-9; callers cannot tell which answered.

craic.accel

Acceleration layer: one pair-HMM posterior primitive, Rust or NumPy.

Design notes

• The compiled craic_accel extension (src/lib.rs) is used when present.
• A NumPy fallback mirrors it line-for-line so the workbench runs anywhere and so the two implementations can be cross-validated in the test-suite (Liskov: callers cannot tell which one answered).
• Emission models are built here, not in Rust, so the core stays biology-agnostic. DNA/RNA use an identity-driven model; protein uses BLOSUM45 converted to a proper joint probability via Robinson background frequencies.

backend

backend() -> str

Source code in craic/accel.py

def backend() -> str:
    return "rust" if HAVE_RUST else "numpy"

pair_posteriors

pair_posteriors(a_idx: Sequence[int], b_idx: Sequence[int], k: int, joint_flat: Sequence[float], bg: Sequence[float], delta: float = 0.02, epsilon: float = 0.5) -> np.ndarray

Dispatch to Rust if available, else the NumPy mirror. Returns (m, n).

delta and epsilon are validated here, in the one place both backends go through. Out-of-range values are not a curiosity: delta >= 0.5 makes log(1 - 2*delta) undefined, which the Rust core turns into NaN that propagates silently through every posterior and out into every reliability score, while the NumPy path raises. A silently wrong number is worse than a loud failure, so both now fail loudly.

Source code in craic/accel.py

def pair_posteriors(
    a_idx: Sequence[int],
    b_idx: Sequence[int],
    k: int,
    joint_flat: Sequence[float],
    bg: Sequence[float],
    delta: float = 0.02,
    epsilon: float = 0.5,
) -> np.ndarray:
    """Dispatch to Rust if available, else the NumPy mirror. Returns (m, n).

    ``delta`` and ``epsilon`` are validated here, in the one place both backends
    go through. Out-of-range values are not a curiosity: ``delta >= 0.5`` makes
    ``log(1 - 2*delta)`` undefined, which the Rust core turns into NaN that
    propagates silently through every posterior and out into every reliability
    score, while the NumPy path raises. A silently wrong number is worse than a
    loud failure, so both now fail loudly.
    """
    _check_gap_params(delta, epsilon)
    if HAVE_RUST:
        post, m, n = _rust.pair_posteriors(
            list(a_idx), list(b_idx), k, list(joint_flat), list(bg), delta, epsilon
        )
        return np.asarray(post, dtype=float).reshape(m, n)
    return _pp_numpy(a_idx, b_idx, k, joint_flat, bg, delta, epsilon)

---

The built-in aligner

craic.progressive

A small, self-contained progressive aligner (the always-available engine).

ProbCons-style, built entirely on CRAIC's pair-HMM posteriors:

1. estimate the emission divergence and gap parameters from the data (no hard-coded assumption that sequences are ~90% identical),
2. compute all pairwise posterior "match" matrices,
3. apply the ProbCons consistency transformation (re-estimate each pairwise posterior using every third sequence), and
4. align progressively by pure maximum-expected-accuracy decoding of the transformed posteriors -- gap placement falls out of the model, so there are no ad-hoc gap constants.

It is still deliberately modest (no sequence weighting, no iterative refinement) and is meant to always be available, not to beat MAFFT; but it is genuine posterior decoding rather than a heuristic. The pair-HMM inner loop runs in the Rust core when built, with the NumPy fallback otherwise.

align

align(records: Sequence[Record], alphabet: Alphabet, model: Optional[EmissionModel] = None, delta: Optional[float] = None, epsilon: Optional[float] = None, effort: str = 'med', consistency_iters: Optional[int] = None, refine_iters: Optional[int] = None, guide: Optional[str] = None, estimate: bool = True, consistency_mem_gb: float = 1.0, progress=None, cancelled=None, guide_seed: Optional[int] = None, sparse: Optional[bool] = None, matrix: Optional[str] = None) -> Alignment

Progressive multiple alignment by posterior decoding, tiered by memory.

Small/medium families get the full ProbCons consistency transformation (accurate). Families whose dense pairwise posteriors would exceed consistency_mem_gb fall back automatically to streaming plain-MEA: each pairwise posterior is computed on demand and discarded, so peak memory is a single L x L matrix and the aligner scales to large inputs (less accurate; for large jobs an external engine is still preferable). progress(done, total) is called before each merge and cancelled() is polled throughout. matrix names the protein substitution matrix (default accel.PROTEIN_MATRIX); it is ignored for nucleotides and when model is given.

Source code in craic/progressive.py

def align(
    records: Sequence[Record],
    alphabet: Alphabet,
    model: Optional[EmissionModel] = None,
    delta: Optional[float] = None,
    epsilon: Optional[float] = None,
    effort: str = "med",
    consistency_iters: Optional[int] = None,
    refine_iters: Optional[int] = None,
    guide: Optional[str] = None,
    estimate: bool = True,
    consistency_mem_gb: float = 1.0,
    progress=None,
    cancelled=None,
    guide_seed: Optional[int] = None,
    sparse: Optional[bool] = None,
    matrix: Optional[str] = None,
) -> Alignment:
    """Progressive multiple alignment by posterior decoding, tiered by memory.

    Small/medium families get the full ProbCons **consistency transformation**
    (accurate). Families whose dense pairwise posteriors would exceed
    ``consistency_mem_gb`` fall back automatically to **streaming plain-MEA**: each
    pairwise posterior is computed on demand and discarded, so peak memory is a
    single L x L matrix and the aligner scales to large inputs (less accurate; for
    large jobs an external engine is still preferable). ``progress(done, total)``
    is called before each merge and ``cancelled()`` is polled throughout.
    ``matrix`` names the protein substitution matrix (default
    ``accel.PROTEIN_MATRIX``); it is ignored for nucleotides and when ``model``
    is given.
    """
    preset = _EFFORT.get(effort, _EFFORT["med"])
    consistency_iters = preset["consistency_iters"] if consistency_iters is None else consistency_iters
    refine_iters = preset["refine_iters"] if refine_iters is None else refine_iters
    guide = preset["guide"] if guide is None else guide

    ids = [r[0] for r in records]
    seqs = [r[1].replace("-", "") for r in records]
    if not seqs:
        raise ValueError("no sequences to align")
    if len(set(ids)) != len(ids):
        # Rows are reordered back onto the input order through a dict keyed on
        # id, and the sandbox splices by id too, so duplicates silently drop
        # sequences rather than failing.
        dup = next(i for i in ids if ids.count(i) > 1)
        raise ValueError(f"duplicate sequence id {dup!r}; ids must be unique")
    if len(seqs) == 1:
        return Alignment(list(ids), list(seqs), alphabet)

    kind = "protein" if alphabet == Alphabet.PROTEIN else "dna"
    matrix = matrix or accel.PROTEIN_MATRIX
    n = len(seqs)
    maxlen = max(len(s) for s in seqs)
    will_sparsify = sparse_available() and maxlen >= SPARSE_MIN_LEN
    use_consistency = (consistency_iters > 0
                       and _posterior_gb(n, maxlen, will_sparsify) <= consistency_mem_gb)

    if model is None and estimate:
        # streaming plain-MEA pilot to read off divergence + gap rates (memory-safe)
        pm = emission_model(kind, matrix=matrix)

        def pilot_post(x, y):
            return accel.posterior_matrix(seqs[x], seqs[y], pm, _DEFAULT_DELTA, _DEFAULT_EPSILON)

        pilot = [r for _, r in _progressive(seqs, pilot_post, cancelled=cancelled)]
        theta = round(min(0.95, max(0.55, _estimate_identity(pilot))), 2)
        gdelta, gepsilon = _estimate_gaps(pilot)
        model = emission_model(kind, theta, matrix)
        if delta is None:
            delta = gdelta
        if epsilon is None:
            epsilon = gepsilon

    if model is None:
        model = emission_model(kind, matrix=matrix)
    if delta is None:
        delta = _DEFAULT_DELTA
    if epsilon is None:
        epsilon = _DEFAULT_EPSILON

    if use_consistency:
        raw = _all_pairs_posteriors(seqs, model, delta, epsilon)
        D = _posterior_distances(seqs, raw) if guide == "posterior" else None
        P = consistency_transform(raw, n, consistency_iters, cancelled, sparse=sparse)

        def post(x, y):
            return _post(P, x, y)
    else:
        D = None                                     # streaming -> cheap k-mer guide tree
        def post(x, y):
            return accel.posterior_matrix(seqs[x], seqs[y], model, delta, epsilon)

    final = _progressive(seqs, post, dist=D, guide_seed=guide_seed,
                         progress=progress, cancelled=cancelled)
    if refine_iters > 0 and use_consistency:
        rng = np.random.default_rng(0 if guide_seed is None else guide_seed)
        final = _refine(final, post, refine_iters, rng, cancelled)

    order = {i: k for k, i in enumerate(ids)}
    pairs = sorted(((ids[gi], row) for gi, row in final), key=lambda p: order[p[0]])
    return Alignment([p[0] for p in pairs], [p[1] for p in pairs], alphabet)

estimate_params

estimate_params(rows: Sequence[str], alphabet: Alphabet)

Emission model and gap parameters read off a committed alignment.

Returns (model, delta, epsilon). This is the single place that turns an alignment into pair-HMM parameters, so that a column is always scored under the model the data implies rather than under a fixed default. The reliability scores and the aligner's own pilot pass both go through it.

Note that for protein data emission_model is currently independent of the estimated identity — one matrix, BLOSUM45, is used at every divergence — so theta affects nucleotide models only. The gap parameters are estimated for both alphabets. The clamp on theta below means it carries no information under 55% identity in any case, which is where most protein reference data sits; keying the protein matrix to divergence would have to read the raw identity, not this.

Source code in craic/progressive.py

def estimate_params(rows: Sequence[str], alphabet: Alphabet):
    """Emission model and gap parameters read off a committed alignment.

    Returns ``(model, delta, epsilon)``. This is the single place that turns an
    alignment into pair-HMM parameters, so that a column is always *scored*
    under the model the data implies rather than under a fixed default. The
    reliability scores and the aligner's own pilot pass both go through it.

    Note that for protein data ``emission_model`` is currently independent of
    the estimated identity — one matrix, BLOSUM45, is used at every divergence —
    so ``theta`` affects nucleotide models only. The gap parameters are
    estimated for both alphabets. The clamp on ``theta`` below means it carries
    no information under 55% identity in any case, which is where most protein
    reference data sits; keying the protein matrix to divergence would have to
    read the raw identity, not this.
    """
    kind = "protein" if alphabet == Alphabet.PROTEIN else "dna"
    theta = round(min(0.95, max(0.55, _estimate_identity(rows))), 2)
    delta, epsilon = _estimate_gaps(rows)
    return emission_model(kind, theta), delta, epsilon

consistency_transform

consistency_transform(P, n, iters=2, cancelled=None, sparse: Optional[bool] = None, threshold: float = SPARSE_THRESHOLD)

ProbCons consistency transformation: re-estimate each pairwise posterior P(x,y) using every third sequence z, P'(x,y) = (1/n) sum_z P(x,z) P(z,y).

This is the most expensive thing CRAIC does — O(n^3) matrix products of L x L matrices per iteration — and it is what sets the largest family the workbench can handle with consistency turned on. Posteriors are overwhelmingly near-zero, so when SciPy is available the matrices are thresholded at threshold and the products run sparse, which is how ProbCons itself does it. sparse=False forces the dense implementation; the two are compared directly in the test-suite rather than assumed equivalent.

Dropping entries below threshold does change the numbers slightly. It is a deliberate approximation, not a refactor, which is why the dense path is kept and why the tests check that the resulting alignments agree.

Source code in craic/progressive.py

def consistency_transform(P, n, iters=2, cancelled=None,
                          sparse: Optional[bool] = None,
                          threshold: float = SPARSE_THRESHOLD):
    """ProbCons consistency transformation: re-estimate each pairwise posterior
    P(x,y) using every third sequence z, P'(x,y) = (1/n) sum_z P(x,z) P(z,y).

    This is the most expensive thing CRAIC does — ``O(n^3)`` matrix products of
    ``L x L`` matrices per iteration — and it is what sets the largest family the
    workbench can handle with consistency turned on. Posteriors are overwhelmingly
    near-zero, so when SciPy is available the matrices are thresholded at
    ``threshold`` and the products run sparse, which is how ProbCons itself does
    it. ``sparse=False`` forces the dense implementation; the two are compared
    directly in the test-suite rather than assumed equivalent.

    Dropping entries below ``threshold`` does change the numbers slightly. It is
    a deliberate approximation, not a refactor, which is why the dense path is
    kept and why the tests check that the resulting *alignments* agree.
    """
    if sparse is None:
        # Auto: sparse only where it actually pays. See SPARSE_MIN_LEN.
        biggest = max((max(v.shape) for v in P.values()), default=0)
        sparse = sparse_available() and biggest >= SPARSE_MIN_LEN
    if sparse and not sparse_available():
        raise RuntimeError("the sparse consistency transform needs SciPy "
                           "(pip install scipy), or pass sparse=False")
    if sparse:
        P = {k: _sparsify(v, threshold) for k, v in P.items()}

    for _ in range(max(0, iters)):
        new = {}
        for x in range(n):
            for y in range(x + 1, n):
                if cancelled is not None and cancelled():
                    raise Cancelled()
                acc = 2.0 * P[(x, y)]                    # z = x and z = y (identity) terms
                for z in range(n):
                    if z == x or z == y:
                        continue
                    a = P[(x, z)] if x < z else P[(z, x)].T
                    b = P[(z, y)] if z < y else P[(y, z)].T
                    acc = acc + a @ b
                acc = acc / float(n)
                # Re-sparsify each iteration: a sum of sparse products fills in,
                # so without this the second pass is effectively dense again and
                # the memory saving evaporates. ProbCons does the same.
                new[(x, y)] = _sparsify(acc.toarray(), threshold) if sparse else acc
        P = new
    return P

sparse_available

sparse_available() -> bool

Whether SciPy is present, which is what the sparse transform needs.

Source code in craic/progressive.py

def sparse_available() -> bool:
    """Whether SciPy is present, which is what the sparse transform needs."""
    try:
        import scipy.sparse  # noqa: F401
    except Exception:
        return False
    return True

---

Reliability

craic.ambiguity.reliability

Per-column / per-residue reliability + live masking.

Two complementary, well-grounded signals:

• consistency (TCS / heads-style): for each alignment column, average the pair-HMM posterior P(res_i ~ res_j) over the residue pairs that column asserts are homologous. High = the unaligned-sequence evidence supports the column. Uses the CRAIC core directly.
• perturbation (guide-tree / gap-regime sensitivity, in the spirit of GUIDANCE but not equivalent to it): re-align the same sequences under an ensemble of perturbed guide trees across a few gap regimes, and measure what fraction of each residue's asserted homologies survive. See :func:perturbation for what the ensemble does and does not sample.

Both live in [0, 1]; combined averages them. A threshold turns either into a column mask that the viewer previews live before you commit.

Two properties this module is careful about, because getting either wrong produces a plausible-looking number that is wrong:

• Both scores are computed under parameters estimated from the alignment being scored (:func:craic.progressive.estimate_params), not under a fixed default. Scoring a 40%-identity family under a 90%-identity model makes every column look unreliable and biases masking toward destroying exactly the divergent data that masking decisions are about.
• A failed ensemble is an error, not a score of zero. If no replicate completes, :func:perturbation raises :class:EnsembleFailure rather than returning all-nan, which downstream would be indistinguishable from "every column is maximally unreliable" and would mask the whole alignment away.

EnsembleFailure

Bases: RuntimeError

Every replicate in the perturbation ensemble failed.

Raised rather than returning an all-nan score, which a threshold would read as "mask everything".

Source code in craic/ambiguity/reliability.py

class EnsembleFailure(RuntimeError):
    """Every replicate in the perturbation ensemble failed.

    Raised rather than returning an all-nan score, which a threshold would read
    as "mask everything".
    """

analyse

analyse(aln: Alignment, do_perturbation: bool = True, model: Optional[EmissionModel] = None, n_replicates: int = 16) -> Reliability

Source code in craic/ambiguity/reliability.py

def analyse(
    aln: Alignment,
    do_perturbation: bool = True,
    model: Optional[accel.EmissionModel] = None,
    n_replicates: int = 16,
) -> Reliability:
    col_c, cell_c = consistency(aln, model=model)
    stats: Dict[str, int] = {}
    if do_perturbation:
        col_p, cell_p = perturbation(aln, n_replicates=n_replicates, stats=stats)
    else:
        col_p = np.full(aln.length, np.nan)
        cell_p = np.full((aln.n_seqs, aln.length), np.nan)
    stack = np.vstack([col_c, col_p])
    combined = _col_nanmean(stack)
    return Reliability(
        aln.length, col_c, col_p, combined, cell_c, cell_p,
        n_replicates_ok=stats.get("n_ok", 0),
        n_replicates_failed=stats.get("n_failed", 0),
    )

consistency

consistency(aln: Alignment, model: Optional[EmissionModel] = None, max_pairs: int = 300, seed: int = 0, delta: Optional[float] = None, epsilon: Optional[float] = None) -> Tuple[np.ndarray, np.ndarray]

Return (col_score[length], cell_score[n, length] with nan at gaps).

model, delta and epsilon default to the values estimated from aln itself rather than to library defaults; pass them explicitly only to score under a model of your own choosing.

Source code in craic/ambiguity/reliability.py

def consistency(
    aln: Alignment,
    model: Optional[accel.EmissionModel] = None,
    max_pairs: int = 300,
    seed: int = 0,
    delta: Optional[float] = None,
    epsilon: Optional[float] = None,
) -> Tuple[np.ndarray, np.ndarray]:
    """Return (col_score[length], cell_score[n, length] with nan at gaps).

    ``model``, ``delta`` and ``epsilon`` default to the values estimated from
    ``aln`` itself rather than to library defaults; pass them explicitly only
    to score under a model of your own choosing.
    """
    est_model, est_delta, est_epsilon = progressive.estimate_params(aln.rows, aln.alphabet)
    if model is None:
        model = est_model
    if delta is None:
        delta = est_delta
    if epsilon is None:
        epsilon = est_epsilon

    n = aln.n_seqs
    L = aln.length
    maps = [domain.residue_index(r) for r in aln.rows]

    pairs = [(i, j) for i in range(n) for j in range(i + 1, n)]
    if len(pairs) > max_pairs:
        rng = np.random.default_rng(seed)
        pairs = [pairs[k] for k in rng.choice(len(pairs), max_pairs, replace=False)]

    # accumulate per-cell sums and counts
    cell_sum = np.zeros((n, L))
    cell_cnt = np.zeros((n, L))
    col_sum = np.zeros(L)
    col_cnt = np.zeros(L)

    for i, j in pairs:
        P = accel.posterior_matrix(aln.rows[i], aln.rows[j], model, delta, epsilon)
        if P.size == 0:
            continue
        mi, mj = maps[i], maps[j]
        for c in range(L):
            ri, rj = mi[c], mj[c]
            if ri < 0 or rj < 0:
                continue
            p = float(P[ri, rj])
            col_sum[c] += p
            col_cnt[c] += 1
            cell_sum[i, c] += p
            cell_cnt[i, c] += 1
            cell_sum[j, c] += p
            cell_cnt[j, c] += 1

    col = np.divide(col_sum, col_cnt, out=np.full(L, np.nan), where=col_cnt > 0)
    cell = np.divide(cell_sum, cell_cnt, out=np.full((n, L), np.nan), where=cell_cnt > 0)
    return col, cell

perturbation

perturbation(aln: Alignment, alphabet: Optional[Alphabet] = None, n_replicates: int = 16, seed: int = 0, deltas: Tuple[float, ...] = _PERTURB_DELTAS, stats: Optional[Dict[str, int]] = None) -> Tuple[np.ndarray, np.ndarray]

Return (col_score[length], cell_score[n, length] with nan at gaps).

Re-align the same sequences under n_replicates perturbed guide trees (cycling through the gap regimes in deltas) and, for each residue, score the fraction of its reference-asserted homologies that survive across the ensemble. A residue with no asserted partners (a singleton column) has no residue-pair evidence and is left nan rather than scored a free 1.0.

What this is. A sensitivity analysis of one engine: how stable are the homologies this alignment asserts, when the guide tree is perturbed and the gap model is varied? Low scores mark residues whose placement is an artefact of a particular tree or gap cost.

What this is not. It is not GUIDANCE. GUIDANCE bootstraps alignment columns to build perturbed guide trees and re-aligns with the same aligner that produced the reference, over ~100 replicates. Here the replicates come from CRAIC's built-in engine under plain posterior decoding (no consistency pass), which is a weaker aligner than the one that may have produced the reference and a different aligner if the reference came from MAFFT or PRANK. In that case the score conflates genuine alignment uncertainty with systematic between-method difference. Read it as a stability probe, and read the disagreement map for between-method evidence.

Raises :class:EnsembleFailure if no replicate completes. If some fail, a warning is issued and the surviving replicates are used; analyse records the count on the report.

Source code in craic/ambiguity/reliability.py

def perturbation(
    aln: Alignment,
    alphabet: Optional[Alphabet] = None,
    n_replicates: int = 16,
    seed: int = 0,
    deltas: Tuple[float, ...] = _PERTURB_DELTAS,
    stats: Optional[Dict[str, int]] = None,
) -> Tuple[np.ndarray, np.ndarray]:
    """Return (col_score[length], cell_score[n, length] with nan at gaps).

    Re-align the same sequences under ``n_replicates`` perturbed guide trees
    (cycling through the gap regimes in ``deltas``) and, for each residue, score
    the fraction of its reference-asserted homologies that survive across the
    ensemble. A residue with no asserted partners (a singleton column) has no
    residue-pair evidence and is left nan rather than scored a free 1.0.

    **What this is.** A sensitivity analysis of one engine: how stable are the
    homologies this alignment asserts, when the guide tree is perturbed and the
    gap model is varied? Low scores mark residues whose placement is an artefact
    of a particular tree or gap cost.

    **What this is not.** It is not GUIDANCE. GUIDANCE bootstraps alignment
    columns to build perturbed guide trees and re-aligns with the *same* aligner
    that produced the reference, over ~100 replicates. Here the replicates come
    from CRAIC's built-in engine under plain posterior decoding (no consistency
    pass), which is a weaker aligner than the one that may have produced the
    reference and a *different* aligner if the reference came from MAFFT or
    PRANK. In that case the score conflates genuine alignment uncertainty with
    systematic between-method difference. Read it as a stability probe, and read
    the disagreement map for between-method evidence.

    Raises :class:`EnsembleFailure` if no replicate completes. If some fail, a
    warning is issued and the surviving replicates are used; ``analyse`` records
    the count on the report.
    """
    alphabet = alphabet or aln.alphabet
    seqs = [domain.ungap(r) for r in aln.rows]
    records = [(str(i), s) for i, s in enumerate(seqs)]

    model, _est_delta, _est_epsilon = progressive.estimate_params(aln.rows, alphabet)

    wanted = max(1, n_replicates)
    alts: List[Alignment] = []
    first_error: Optional[BaseException] = None
    for i in range(wanted):
        try:
            alts.append(progressive.align(
                records, alphabet, model=model,
                delta=deltas[i % len(deltas)], epsilon=0.5,
                estimate=False, consistency_iters=0, guide_seed=seed + i + 1))
        except Exception as exc:          # noqa: BLE001 - recorded, then re-raised if total
            if first_error is None:
                first_error = exc
    n_failed = wanted - len(alts)
    if not alts:
        raise EnsembleFailure(
            f"all {wanted} perturbation replicates failed; "
            f"first error: {first_error!r}"
        ) from first_error
    if n_failed:
        warnings.warn(
            f"{n_failed} of {wanted} perturbation replicates failed; "
            "scores are based on the remainder",
            RuntimeWarning,
            stacklevel=2,
        )
    if stats is not None:
        stats["n_ok"], stats["n_failed"] = len(alts), n_failed

    n = aln.n_seqs
    L = aln.length

    # reference membership uses *positional* sequence index; alts use str(index)
    def mem_by_index(a: Alignment) -> Dict[ResKey, int]:
        mem: Dict[ResKey, int] = {}
        for row_id, row in zip(a.ids, a.rows):
            si = int(row_id)
            r = 0
            for c, ch in enumerate(row):
                if ch != "-":
                    mem[(si, r)] = c
                    r += 1
        return mem

    ref_mem = {(si, r): c for (si, r), c in domain.membership(aln).items()}
    alt_mems = [mem_by_index(a) for a in alts]

    # partners asserted by the reference, per residue
    partners: Dict[ResKey, List[ResKey]] = defaultdict(list)
    maps = [domain.residue_index(r) for r in aln.rows]
    for c in range(L):
        present = [(si, maps[si][c]) for si in range(n) if maps[si][c] >= 0]
        for a_idx in range(len(present)):
            for b_idx in range(len(present)):
                if a_idx != b_idx:
                    partners[present[a_idx]].append(present[b_idx])

    res_score: Dict[ResKey, float] = {}
    for x, plist in partners.items():
        if not alt_mems or not plist:
            continue  # no ensemble, or a singleton column: no residue-pair evidence -> nan
        tot = ok = 0
        for y in plist:
            for mem in alt_mems:
                tot += 1
                cx, cy = mem.get(x), mem.get(y)
                if cx is not None and cx == cy:
                    ok += 1
        if tot:
            res_score[x] = ok / tot

    cell = np.full((n, L), np.nan)
    for (si, r), score in res_score.items():
        c = ref_mem.get((si, r))
        if c is not None:
            cell[si, c] = score
    col = _col_nanmean(cell)
    return col, cell

keep_mask

keep_mask(scores: ndarray, threshold: float, unscored: str = 'drop') -> np.ndarray

Columns to keep: scores >= threshold.

This is the one definition of a reliability mask in CRAIC. It previously existed twice with opposite handling of unscoreable (nan) columns — the GUI kept them, the benchmark dropped them — so the mask a user exported was not the mask the benchmark evaluated.

A nan column is one with no residue-pair evidence at all (a singleton column, or one whose ensemble evidence is missing). unscored="drop" is the default and the conservative choice for downstream phylogenetics: a column that could not be assessed is not evidence. unscored="keep" retains them.

Raises ValueError if every column is unscoreable, which means the reliability analysis failed rather than that the alignment is worthless.

Source code in craic/ambiguity/reliability.py

def keep_mask(
    scores: np.ndarray,
    threshold: float,
    unscored: str = "drop",
) -> np.ndarray:
    """Columns to keep: ``scores >= threshold``.

    This is the *one* definition of a reliability mask in CRAIC. It previously
    existed twice with opposite handling of unscoreable (nan) columns — the GUI
    kept them, the benchmark dropped them — so the mask a user exported was not
    the mask the benchmark evaluated.

    A nan column is one with no residue-pair evidence at all (a singleton
    column, or one whose ensemble evidence is missing). ``unscored="drop"`` is
    the default and the conservative choice for downstream phylogenetics: a
    column that could not be assessed is not evidence. ``unscored="keep"``
    retains them.

    Raises ``ValueError`` if *every* column is unscoreable, which means the
    reliability analysis failed rather than that the alignment is worthless.
    """
    if unscored not in ("drop", "keep"):
        raise ValueError(f"unscored must be 'drop' or 'keep', not {unscored!r}")
    scores = np.asarray(scores, dtype=float)
    if scores.size and bool(np.all(np.isnan(scores))):
        raise ValueError(
            "every column is unscored, so no threshold is meaningful; "
            "the reliability analysis did not produce usable scores"
        )
    fill = 0.0 if unscored == "drop" else 1.0
    return np.nan_to_num(scores, nan=fill) >= threshold

---

Accuracy against a known answer

craic.evaluate

Accuracy of an alignment against a known-true reference.

Sum-of-pairs and total-column scores, per-column correctness, and the ROC AUC of a score predicting that correctness. These used to live in benchmarks/ where only the benchmark harness could reach them; they are here so that the GUI's truth mode, the command line and the benchmark all compute accuracy the same way. A number shown to a student and a number printed in a paper being produced by two different implementations is exactly the drift worth designing out.

A "reference" here is a trusted alignment of the same sequences — a structural reference such as BAliBASE, or the known-true alignment of a simulated dataset.

Accuracy dataclass

How an inferred alignment compares with a trusted reference.

Source code in craic/evaluate.py

@dataclass
class Accuracy:
    """How an inferred alignment compares with a trusted reference."""

    sp: float                      #: recall of the reference's homologous pairs
    precision: float               #: fraction of asserted pairs that are true
    tc: float                      #: fraction of reference columns reproduced exactly
    col_correct: np.ndarray        #: per inferred column, fraction of its pairs that are true
    n_scorable: int                #: columns that assert at least one pair

    @property
    def correct_mask(self) -> np.ndarray:
        """Columns that are entirely correct (nan columns count as not-correct,
        since they assert nothing to be correct about)."""
        return np.nan_to_num(self.col_correct, nan=0.0) >= 0.999

    def reliability_auc(self, scores: np.ndarray) -> float:
        """How well a per-column score predicts which columns are actually right."""
        scores = np.asarray(scores, dtype=float)
        ok = ~np.isnan(scores) & ~np.isnan(self.col_correct)
        if ok.sum() < 4:
            return float("nan")
        return auc(scores[ok], self.col_correct[ok] >= 0.999)

correct_mask property

correct_mask: ndarray

Columns that are entirely correct (nan columns count as not-correct, since they assert nothing to be correct about).

reliability_auc

reliability_auc(scores: ndarray) -> float

How well a per-column score predicts which columns are actually right.

Source code in craic/evaluate.py

def reliability_auc(self, scores: np.ndarray) -> float:
    """How well a per-column score predicts which columns are actually right."""
    scores = np.asarray(scores, dtype=float)
    ok = ~np.isnan(scores) & ~np.isnan(self.col_correct)
    if ok.sum() < 4:
        return float("nan")
    return auc(scores[ok], self.col_correct[ok] >= 0.999)

ColumnTruth dataclass

What the reference says about one inferred column.

agree are the residues the reference groups together and the alignment got right; intruders are residues placed here that belong elsewhere; and missing are residues that belong here but were placed elsewhere. A correct column has neither of the latter two.

Source code in craic/evaluate.py

@dataclass
class ColumnTruth:
    """What the reference says about one inferred column.

    ``agree`` are the residues the reference groups together and the alignment
    got right; ``intruders`` are residues placed here that belong elsewhere; and
    ``missing`` are residues that belong here but were placed elsewhere. A
    correct column has neither of the latter two.
    """

    col: int
    true_col: Optional[int]        #: the reference column this one mostly represents
    agree: List[Placement] = field(default_factory=list)
    intruders: List[Placement] = field(default_factory=list)
    missing: List[Placement] = field(default_factory=list)

    @property
    def correct(self) -> bool:
        return not self.intruders and not self.missing

    def describe(self) -> str:
        """The same thing in words, for the command line and for tooltips."""
        if self.true_col is None:
            return f"Column {self.col + 1} is empty."
        def names(ps):
            return ", ".join(f"{p.seq_id}/{p.char}{p.residue + 1}" for p in ps)
        lines = [f"Column {self.col + 1} — reference column {self.true_col + 1}"]
        if self.agree:
            lines.append(f"Correctly grouped: {names(self.agree)}")
        if self.intruders:
            for p in self.intruders:
                where = ("no true partners in this alignment" if p.offset is None
                         else f"belongs {abs(p.offset)} column"
                              f"{'s' if abs(p.offset) != 1 else ''} to the "
                              f"{'right' if p.offset > 0 else 'left'}, "
                              f"in column {p.target_col + 1}")
                lines.append(f"Does not belong here: {p.seq_id}/{p.char}"
                             f"{p.residue + 1} — {where}")
        if self.missing:
            for p in self.missing:
                lines.append(f"Should be here but is not: {p.seq_id}/{p.char}"
                             f"{p.residue + 1} — currently in column {p.col + 1}")
        if self.correct:
            lines.append("This column is exactly right.")
        return "\n".join(lines)

describe

describe() -> str

The same thing in words, for the command line and for tooltips.

Source code in craic/evaluate.py

def describe(self) -> str:
    """The same thing in words, for the command line and for tooltips."""
    if self.true_col is None:
        return f"Column {self.col + 1} is empty."
    def names(ps):
        return ", ".join(f"{p.seq_id}/{p.char}{p.residue + 1}" for p in ps)
    lines = [f"Column {self.col + 1} — reference column {self.true_col + 1}"]
    if self.agree:
        lines.append(f"Correctly grouped: {names(self.agree)}")
    if self.intruders:
        for p in self.intruders:
            where = ("no true partners in this alignment" if p.offset is None
                     else f"belongs {abs(p.offset)} column"
                          f"{'s' if abs(p.offset) != 1 else ''} to the "
                          f"{'right' if p.offset > 0 else 'left'}, "
                          f"in column {p.target_col + 1}")
            lines.append(f"Does not belong here: {p.seq_id}/{p.char}"
                         f"{p.residue + 1} — {where}")
    if self.missing:
        for p in self.missing:
            lines.append(f"Should be here but is not: {p.seq_id}/{p.char}"
                         f"{p.residue + 1} — currently in column {p.col + 1}")
    if self.correct:
        lines.append("This column is exactly right.")
    return "\n".join(lines)

Placement dataclass

One residue, where it sits, and where the reference says it belongs.

Source code in craic/evaluate.py

@dataclass
class Placement:
    """One residue, where it sits, and where the reference says it belongs."""

    seq: int                       #: row index in the inferred alignment
    seq_id: str                    #: that row's identifier
    residue: int                   #: index in the ungapped sequence, 0-based
    char: str                      #: the residue itself
    true_col: int                  #: the column the reference puts it in
    col: int                       #: the column it currently occupies
    target_col: Optional[int] = None
    """The inferred column holding most of its true partners — where it would
    have to move to join them. ``None`` if it has no true partners to join."""

    @property
    def offset(self) -> Optional[int]:
        """Columns to the right (negative: to the left) it would have to move."""
        return None if self.target_col is None else self.target_col - self.col

target_col class-attribute instance-attribute

target_col: Optional[int] = None

The inferred column holding most of its true partners — where it would have to move to join them. None if it has no true partners to join.

offset property

offset: Optional[int]

Columns to the right (negative: to the left) it would have to move.

compare_to_reference

compare_to_reference(aln: Alignment, reference: Alignment, core: Optional[ndarray] = None) -> Accuracy

Score aln against a trusted reference alignment of the same sequences.

See :func:matched_rows on how the two are paired up.

Source code in craic/evaluate.py

def compare_to_reference(aln: Alignment, reference: Alignment,
                         core: Optional[np.ndarray] = None) -> Accuracy:
    """Score ``aln`` against a trusted ``reference`` alignment of the same sequences.

    See :func:`matched_rows` on how the two are paired up.
    """
    true_rows = matched_rows(aln, reference)
    inf_rows = list(aln.rows)
    sp, prec = sp_score(true_rows, inf_rows, core)
    tc = tc_score(true_rows, inf_rows, core)
    frac = col_correct_fraction(inf_rows, true_rows)
    return Accuracy(sp=sp, precision=prec, tc=tc, col_correct=frac,
                    n_scorable=int((~np.isnan(frac)).sum()))

matched_rows

matched_rows(aln: Alignment, reference: Alignment) -> List[str]

The reference's rows, reordered onto aln's own sequence order.

The two are matched on sequence id (tolerating the whitespace truncation some aligners apply), so that a reference whose sequences are listed differently — which is usual — is still usable. If the two do not describe the same sequences, or describe them with different residues, that is an error rather than a very low score: a reference for the wrong data would otherwise read as an alignment that is entirely wrong.

Source code in craic/evaluate.py

def matched_rows(aln: Alignment, reference: Alignment) -> List[str]:
    """The reference's rows, reordered onto ``aln``'s own sequence order.

    The two are matched on sequence id (tolerating the whitespace truncation some
    aligners apply), so that a reference whose sequences are listed differently —
    which is usual — is still usable. If the two do not describe the same
    sequences, or describe them with different residues, that is an error rather
    than a very low score: a reference for the wrong data would otherwise read as
    an alignment that is entirely wrong.
    """
    mapping = match_ids(aln, reference)                     # raises IdMismatch
    by_id = dict(zip(reference.ids, reference.rows))
    true_rows = [by_id[mapping[rid]] for rid in aln.ids]

    for rid, inferred, truth in zip(aln.ids, aln.rows, true_rows):
        if ungap(inferred).upper() != ungap(truth).upper():
            raise IdMismatch(
                f"sequence {rid!r} differs between the alignment and the reference "
                f"({len(ungap(inferred))} vs {len(ungap(truth))} residues); "
                "the reference must be an alignment of the same sequences"
            )
    return true_rows

sp_score

sp_score(true_rows: Sequence[str], inf_rows: Sequence[str], core: Optional[ndarray] = None) -> Tuple[float, float]

(recall, precision) of homologous residue pairs. Recall is the SP score.

core is a boolean mask over the reference columns. Published BAliBASE SP figures are restricted to the reference's core blocks, so a score computed over all columns is not comparable with them — pass the core mask when one is available, and label the result accordingly when it is not.

Source code in craic/evaluate.py

def sp_score(true_rows: Sequence[str], inf_rows: Sequence[str],
             core: Optional[np.ndarray] = None) -> Tuple[float, float]:
    """(recall, precision) of homologous residue pairs. Recall is the SP score.

    ``core`` is a boolean mask over the *reference* columns. Published BAliBASE
    SP figures are restricted to the reference's core blocks, so a score
    computed over all columns is not comparable with them — pass the core mask
    when one is available, and label the result accordingly when it is not.
    """
    T, I = _pairs(true_rows, core), _pairs(inf_rows)
    inter = len(T & I)
    recall = inter / len(T) if T else 1.0
    prec = inter / len(I) if I else 1.0
    return recall, prec

tc_score

tc_score(true_rows: Sequence[str], inf_rows: Sequence[str], core: Optional[ndarray] = None) -> float

Total-column score: fraction of reference columns reproduced exactly.

See :func:sp_score on core.

Source code in craic/evaluate.py

def tc_score(true_rows: Sequence[str], inf_rows: Sequence[str],
             core: Optional[np.ndarray] = None) -> float:
    """Total-column score: fraction of reference columns reproduced exactly.

    See :func:`sp_score` on ``core``.
    """
    T, I = _columns(true_rows, core), _columns(inf_rows)
    return len(T & I) / len(T) if T else 1.0

col_correct_fraction

col_correct_fraction(inf_rows: Sequence[str], true_rows: Sequence[str]) -> np.ndarray

Per inferred column, the fraction of its asserted pairs that are true.

nan for a column that asserts no pair at all (a singleton column): it makes no claim, so it can be neither right nor wrong. Those columns are not counted as errors anywhere downstream.

Source code in craic/evaluate.py

def col_correct_fraction(inf_rows: Sequence[str], true_rows: Sequence[str]) -> np.ndarray:
    """Per inferred column, the fraction of its asserted pairs that are true.

    ``nan`` for a column that asserts no pair at all (a singleton column): it
    makes no claim, so it can be neither right nor wrong. Those columns are not
    counted as errors anywhere downstream.
    """
    true_pairs = _pairs(true_rows)
    ridx = [residue_index(r) for r in inf_rows]
    out = []
    for c in range(len(inf_rows[0])):
        present = [(si, ridx[si][c]) for si in range(len(inf_rows)) if inf_rows[si][c] != "-"]
        pairs = [(present[a], present[b])
                 for a in range(len(present)) for b in range(a + 1, len(present))]
        out.append(np.mean([p in true_pairs for p in pairs]) if pairs else np.nan)
    return np.array(out)

cell_correct_fraction

cell_correct_fraction(inf_rows: Sequence[str], true_rows: Sequence[str], grid: Optional[ndarray] = None) -> np.ndarray

Per residue, the fraction of its partners in that column that are homologous.

The per-cell counterpart of :func:col_correct_fraction, and exactly consistent with it: a column's score is the mean of its cells' scores. One misplaced residue in an otherwise sound column scores 0 while its neighbours stay near 1, which is what makes the culprit visible rather than merely the damage.

nan for a gap, and for a residue with no partners at all — a residue alone in its column asserts nothing, so it can be neither right nor wrong.

Source code in craic/evaluate.py

def cell_correct_fraction(inf_rows: Sequence[str], true_rows: Sequence[str],
                          grid: Optional[np.ndarray] = None) -> np.ndarray:
    """Per residue, the fraction of its partners in that column that are homologous.

    The per-cell counterpart of :func:`col_correct_fraction`, and exactly
    consistent with it: a column's score is the mean of its cells' scores. One
    misplaced residue in an otherwise sound
    column scores 0 while its neighbours stay near 1, which is what makes the
    culprit visible rather than merely the damage.

    ``nan`` for a gap, and for a residue with no partners at all — a residue
    alone in its column asserts nothing, so it can be neither right nor wrong.
    """
    grid = true_column_grid(inf_rows, true_rows) if grid is None else grid
    out = np.full(grid.shape, np.nan, dtype=float)
    for c in range(grid.shape[1]):
        col = grid[:, c]
        present = np.flatnonzero(col >= 0)
        if present.size < 2:
            continue
        counts = Counter(int(v) for v in col[present])
        for s in present:
            out[s, c] = (counts[int(col[s])] - 1) / (present.size - 1)
    return out

true_column_grid

true_column_grid(inf_rows: Sequence[str], true_rows: Sequence[str]) -> np.ndarray

(n_seqs, n_cols) of true-alignment columns, -1 where a cell is a gap.

Entry [s, c] is the column the reference puts the residue that the inferred alignment placed at row s, column c. Two residues are homologous exactly when their entries are equal, which is the whole of the comparison — every function below is a different way of reading this grid.

Source code in craic/evaluate.py

def true_column_grid(inf_rows: Sequence[str], true_rows: Sequence[str]) -> np.ndarray:
    """``(n_seqs, n_cols)`` of true-alignment columns, ``-1`` where a cell is a gap.

    Entry ``[s, c]`` is the column the *reference* puts the residue that the
    inferred alignment placed at row ``s``, column ``c``. Two residues are
    homologous exactly when their entries are equal, which is the whole of the
    comparison — every function below is a different way of reading this grid.
    """
    n, L = len(inf_rows), len(inf_rows[0])
    grid = np.full((n, L), -1, dtype=int)
    for s, (inferred, truth) in enumerate(zip(inf_rows, true_rows)):
        col_of = column_index(truth)                   # residue index -> true column
        for c, r in enumerate(residue_index(inferred)):
            if r >= 0:
                grid[s, c] = col_of[r]
    return grid

explain_column

explain_column(inf_rows: Sequence[str], true_rows: Sequence[str], col: int, ids: Optional[Sequence[str]] = None, grid: Optional[ndarray] = None) -> ColumnTruth

Account for one inferred column against the reference.

Pass grid from :func:true_column_grid to explain many columns without rebuilding the lookup each time.

Source code in craic/evaluate.py

def explain_column(inf_rows: Sequence[str], true_rows: Sequence[str], col: int,
                   ids: Optional[Sequence[str]] = None,
                   grid: Optional[np.ndarray] = None) -> ColumnTruth:
    """Account for one inferred column against the reference.

    Pass ``grid`` from :func:`true_column_grid` to explain many columns without
    rebuilding the lookup each time.
    """
    grid = true_column_grid(inf_rows, true_rows) if grid is None else grid
    n = len(inf_rows)
    ids = list(ids) if ids is not None else [f"seq_{i + 1}" for i in range(n)]
    ridx = [residue_index(r) for r in inf_rows]

    def placement(s: int, c: int, true_col: int, target: Optional[int] = None) -> Placement:
        return Placement(seq=s, seq_id=ids[s], residue=ridx[s][c], char=inf_rows[s][c],
                         true_col=true_col, col=c, target_col=target)

    here = grid[:, col]
    present = [int(s) for s in np.flatnonzero(here >= 0)]
    if not present:
        return ColumnTruth(col=col, true_col=None)

    # The reference column this one mostly represents. Ties go to the leftmost,
    # so the answer never depends on dictionary order.
    counts = Counter(int(here[s]) for s in present)
    best = max(counts.values())
    true_col = min(t for t, k in counts.items() if k == best)

    def column_of(s: int, t: int) -> Optional[int]:
        """Which inferred column holds sequence ``s``'s residue from true column ``t``."""
        found = np.flatnonzero(grid[s] == t)
        return int(found[0]) if found.size else None

    def target_for(s: int, t: int) -> Optional[int]:
        """Where this residue's true partners currently sit, by majority."""
        elsewhere = Counter()
        for other in range(n):
            if other == s:
                continue
            c = column_of(other, t)
            if c is not None:
                elsewhere[c] += 1
        if not elsewhere:
            return None
        most = max(elsewhere.values())
        return min(c for c, k in elsewhere.items() if k == most)

    out = ColumnTruth(col=col, true_col=true_col)
    for s in present:
        t = int(here[s])
        if t == true_col:
            out.agree.append(placement(s, col, t))
        else:
            out.intruders.append(placement(s, col, t, target_for(s, t)))

    for s in range(n):
        if here[s] == true_col:
            continue
        c = column_of(s, true_col)
        if c is not None:
            out.missing.append(placement(s, c, true_col))
    return out

auc

auc(score: ndarray, positive: ndarray) -> float

ROC AUC: probability a random positive outranks a random negative.

Ties count as half, via midranks, so a score that cannot discriminate at all returns 0.5 rather than an order-dependent number.

Source code in craic/evaluate.py

def auc(score: np.ndarray, positive: np.ndarray) -> float:
    """ROC AUC: probability a random positive outranks a random negative.

    Ties count as half, via midranks, so a score that cannot discriminate at all
    returns 0.5 rather than an order-dependent number.
    """
    score = np.asarray(score, dtype=float)
    positive = np.asarray(positive, dtype=bool)
    n_pos, n_neg = int(positive.sum()), int((~positive).sum())
    if n_pos == 0 or n_neg == 0:
        return float("nan")
    ranks = _midranks(score)
    return float((ranks[positive].sum() - n_pos * (n_pos + 1) / 2.0) / (n_pos * n_neg))

---

Simulation

craic.simulate

Ground-truth sequence simulator.

Used by the benchmark harness and by the workbench itself: a dataset whose true alignment is known exactly is the only way to show, rather than assert, what a reliability score is doing — which makes it as much a teaching device as a benchmarking one.

Evolves DNA down a random binary tree with JC69 substitutions and insertions / deletions, tracking the TRUE alignment (which residues are homologous) and the TRUE topology. Divergence is controlled by the branch-length range and indel load by indel_rate. Because the homology of every residue is known exactly, this is the clean way to evaluate a per-column reliability score.

A caveat that belongs in any write-up using these numbers. The default generative model here — JC69, uniform base frequencies, one rate for every site, geometric indel lengths — is close to the model CRAIC's pair-HMM assumes, so CRAIC's inference model is nearly well-specified on this data while MAFFT and MUSCLE carry empirical matrices tuned for real sequences. That is a structural advantage to CRAIC, not a bug, but it must be disclosed. rate_alpha turns on gamma-distributed among-site rate variation and indel_zipf gives indel lengths a heavy tail, both of which break the correspondence; run the sweep with and without them and report both.

simulate

simulate(taxa: int = 8, root_len: int = 200, seed: int = 0, indel_rate: float = 2.0, bmin: float = 0.04, bmax: float = 0.28, rate_alpha: float = 0.0, indel_zipf: float = 0.0) -> dict

Return dict(names, seqs=[(name, ungapped)], true_rows=[gapped], tree).

rate_alpha > 0 draws each site's relative substitution rate from a gamma(alpha, 1/alpha), i.e. among-site rate variation with mean 1; smaller alpha means more heterogeneity. indel_zipf > 1 replaces the geometric indel-length distribution with a Zipf one. Both default to off, which reproduces earlier results exactly; both should be exercised before claiming a result generalises — see the module docstring.

Source code in craic/simulate.py

def simulate(taxa: int = 8, root_len: int = 200, seed: int = 0,
             indel_rate: float = 2.0, bmin: float = 0.04, bmax: float = 0.28,
             rate_alpha: float = 0.0, indel_zipf: float = 0.0) -> dict:
    """Return dict(names, seqs=[(name, ungapped)], true_rows=[gapped], tree).

    ``rate_alpha`` > 0 draws each site's relative substitution rate from a
    gamma(alpha, 1/alpha), i.e. among-site rate variation with mean 1; smaller
    alpha means more heterogeneity. ``indel_zipf`` > 1 replaces the geometric
    indel-length distribution with a Zipf one. Both default to off, which
    reproduces earlier results exactly; both should be exercised before claiming
    a result generalises — see the module docstring.
    """
    rng = np.random.default_rng(seed)
    tree = random_tree(taxa, rng, bmin, bmax)
    ctr = _Ctr()
    root = [[ctr.new(float(i)), rng.choice(list(BASES))] for i in range(root_len)]
    rates = None
    if rate_alpha and rate_alpha > 0:
        rates = {i: float(rng.gamma(rate_alpha, 1.0 / rate_alpha)) for i, _ in root}
    leaves: Dict[str, list] = {}

    def rec(node, sites):
        if not node["children"]:
            leaves[node["name"]] = sites
            return
        for child, bl in node["children"]:
            kids = _evolve(sites, bl, rng, ctr, indel_rate, rates, indel_zipf)
            if rates is not None:
                for i, _ in kids:                  # newly inserted sites need a rate
                    rates.setdefault(i, float(rng.gamma(rate_alpha, 1.0 / rate_alpha)))
            rec(child, kids)

    rec(tree, root)
    names = sorted(leaves, key=lambda n: int(n[1:]))
    allids = sorted({i for s in leaves.values() for i, _ in s}, key=lambda i: (ctr.pos[i], i))
    colidx = {i: c for c, i in enumerate(allids)}
    rows = []
    for nm in names:
        row = ["-"] * len(allids)
        for i, b in leaves[nm]:
            row[colidx[i]] = b
        rows.append("".join(row))
    seqs = [(nm, "".join(b for _, b in leaves[nm])) for nm in names]
    return dict(names=names, seqs=seqs, true_rows=rows, tree=tree)

---

Model-free trimming

Faithful reimplementations of the published column-selection rules, offered as complementary signals to the model-based reliability score.

craic.ambiguity.trimming

Model-free column trimming and sequence-outlier detection.

Faithful reimplementations of the algorithms behind widely-used tools, offered as complements to CRAIC's model-based posterior reliability (which is derived from the pair-HMM). These use only gaps and residue conservation -- no probabilistic model -- so they are cheap and always available:

• gap_score / gap_threshold_mask -- MSA_trimmer-style gappiness.
• similarity_score -- per-column conservation.
• gappyout_mask / trimal_strict_mask -- trimAl-style automatic trimming.
• gblocks_mask -- Gblocks-style conserved blocks.
• sequence_scores / outlier_sequences -- EvalMSA / OD-seq-style outliers.

These are reimplementations of the published algorithms, not the original code, and are not guaranteed identical to the reference tools; for authoritative filtering, run the originals (trimAl, Gblocks, EvalMSA).

gap_score

gap_score(aln: Alignment) -> np.ndarray

Fraction of non-gap residues per column (trimAl gap score). 1 = no gaps.

Source code in craic/ambiguity/trimming.py

def gap_score(aln: Alignment) -> np.ndarray:
    """Fraction of non-gap residues per column (trimAl gap score). 1 = no gaps."""
    M = _char_matrix(aln)
    return (M != "-").mean(axis=0)

similarity_score

similarity_score(aln: Alignment) -> np.ndarray

Mean pairwise residue similarity per column in [0, 1] (trimAl-style). Gaps excluded; a column with 0/1 residues scores 0/1.

Source code in craic/ambiguity/trimming.py

def similarity_score(aln: Alignment) -> np.ndarray:
    """Mean pairwise residue similarity per column in [0, 1] (trimAl-style).
    Gaps excluded; a column with 0/1 residues scores 0/1."""
    E, sim = _encode(aln)
    n, L = E.shape
    out = np.zeros(L)
    for c in range(L):
        col = E[:, c]
        present = col[col >= 0]
        if present.size < 2:
            out[c] = 1.0 if present.size == 1 else 0.0
            continue
        sub = sim[np.ix_(present, present)]
        iu = np.triu_indices(present.size, 1)
        out[c] = float(sub[iu].mean())
    return out

---

Sessions

craic.session

The working session as a document.

A curation session is more than its alignment. It is also the reference alignment truth mode is scoring against, the column annotations, the named sequence groups, the provenance log of how the alignment was arrived at, and the view state that makes any of it legible. Exporting FASTA or NEXUS keeps the residues and silently discards all of that — fine as an export, useless as a way to put work down and pick it up again.

So a session is its own small document: readable JSON, no dependencies, and versioned so that an old file can be recognised rather than misread. The alignment can still be saved to any standard format at any time; the two are separate acts, because they answer different questions ("give this to another program" versus "let me carry on tomorrow").

The module is deliberately free of Qt, so that the command line, the tests and the crash-recovery autosave can all use it without a display.

Session dataclass

Everything needed to resume a curation session.

Source code in craic/session.py

@dataclass
class Session:
    """Everything needed to resume a curation session."""

    alignment: Optional[Alignment] = None
    #: the trusted alignment truth mode scores against, and what to call it
    reference: Optional[Alignment] = None
    reference_label: str = ""
    annotations: List[dict] = field(default_factory=list)
    groups: Dict[str, List[str]] = field(default_factory=dict)
    #: where the alignment itself was last read from or written to
    source_path: str = ""
    #: view state — which overlay, the mask threshold, the viewing level
    view: Dict[str, Any] = field(default_factory=dict)
    #: set by the autosave; an explicit save clears it
    autosaved: bool = False
    saved_at: float = 0.0
    craic_version: str = ""

    # -- serialisation ----------------------------------------------------- #
    def to_dict(self) -> Dict[str, Any]:
        return {
            "format": "craic-session",
            "format_version": FORMAT_VERSION,
            "craic_version": self.craic_version or __version__,
            "saved_at": self.saved_at or time.time(),
            "autosaved": bool(self.autosaved),
            "source_path": self.source_path,
            "alignment": _alignment_to_dict(self.alignment) if self.alignment else None,
            "reference": _alignment_to_dict(self.reference) if self.reference else None,
            "reference_label": self.reference_label,
            "annotations": list(self.annotations),
            "groups": {k: list(v) for k, v in self.groups.items()},
            "view": dict(self.view),
        }

    @classmethod
    def from_dict(cls, d: Dict[str, Any]) -> "Session":
        if not isinstance(d, dict) or d.get("format") != "craic-session":
            raise SessionError("this file is not a CRAIC session")
        version = d.get("format_version")
        if not isinstance(version, int) or version > FORMAT_VERSION:
            raise SessionError(
                f"this session was written by a newer CRAIC (format {version}); "
                f"this one understands up to {FORMAT_VERSION}")
        return cls(
            alignment=_alignment_from_dict(d.get("alignment")),
            reference=_alignment_from_dict(d.get("reference")),
            reference_label=d.get("reference_label", "") or "",
            annotations=list(d.get("annotations") or []),
            groups={k: list(v) for k, v in (d.get("groups") or {}).items()},
            source_path=d.get("source_path", "") or "",
            view=dict(d.get("view") or {}),
            autosaved=bool(d.get("autosaved")),
            saved_at=float(d.get("saved_at") or 0.0),
            craic_version=d.get("craic_version", "") or "",
        )

    # -- files -------------------------------------------------------------- #
    def save(self, path: str) -> None:
        """Write the session, atomically.

        Via a temporary file and a rename, because the autosave writes this
        repeatedly and a crash during the write is exactly the moment the file
        matters most — a half-written recovery file is worse than none.
        """
        self.saved_at = time.time()
        self.craic_version = __version__
        tmp = path + ".part"
        directory = os.path.dirname(os.path.abspath(path))
        os.makedirs(directory, exist_ok=True)
        with open(tmp, "w", encoding="utf-8") as fh:
            json.dump(self.to_dict(), fh, indent=1)
        os.replace(tmp, path)

    @classmethod
    def load(cls, path: str) -> "Session":
        try:
            with open(path, encoding="utf-8") as fh:
                data = json.load(fh)
        except OSError as exc:
            raise SessionError(f"cannot read {path}: {exc}") from exc
        except ValueError as exc:
            raise SessionError(f"{os.path.basename(path)} is not valid JSON: {exc}") from exc
        return cls.from_dict(data)

    # -- description -------------------------------------------------------- #
    def describe(self) -> str:
        """One-line summary, for a recovery prompt or the command line."""
        if self.alignment is None:
            return "an empty session"
        name = os.path.basename(self.source_path) if self.source_path else "unsaved sequences"
        bits = [f"{name}: {self.alignment.n_seqs} sequences x "
                f"{self.alignment.length} columns"]
        if self.reference is not None:
            bits.append(f"reference: {self.reference_label or 'loaded'}")
        history = (self.alignment.meta or {}).get("history") or []
        if history:
            bits.append(f"{len(history)} steps of history")
        return "; ".join(bits)

save

save(path: str) -> None

Write the session, atomically.

Via a temporary file and a rename, because the autosave writes this repeatedly and a crash during the write is exactly the moment the file matters most — a half-written recovery file is worse than none.

Source code in craic/session.py

def save(self, path: str) -> None:
    """Write the session, atomically.

    Via a temporary file and a rename, because the autosave writes this
    repeatedly and a crash during the write is exactly the moment the file
    matters most — a half-written recovery file is worse than none.
    """
    self.saved_at = time.time()
    self.craic_version = __version__
    tmp = path + ".part"
    directory = os.path.dirname(os.path.abspath(path))
    os.makedirs(directory, exist_ok=True)
    with open(tmp, "w", encoding="utf-8") as fh:
        json.dump(self.to_dict(), fh, indent=1)
    os.replace(tmp, path)

describe

describe() -> str

One-line summary, for a recovery prompt or the command line.

Source code in craic/session.py

def describe(self) -> str:
    """One-line summary, for a recovery prompt or the command line."""
    if self.alignment is None:
        return "an empty session"
    name = os.path.basename(self.source_path) if self.source_path else "unsaved sequences"
    bits = [f"{name}: {self.alignment.n_seqs} sequences x "
            f"{self.alignment.length} columns"]
    if self.reference is not None:
        bits.append(f"reference: {self.reference_label or 'loaded'}")
    history = (self.alignment.meta or {}).get("history") or []
    if history:
        bits.append(f"{len(history)} steps of history")
    return "; ".join(bits)

SessionError

Bases: ValueError

A session file could not be read as one.

Source code in craic/session.py

class SessionError(ValueError):
    """A session file could not be read as one."""

is_session_file

is_session_file(path: str) -> bool

Whether path looks like a session document, cheaply.

Source code in craic/session.py

def is_session_file(path: str) -> bool:
    """Whether ``path`` looks like a session document, cheaply."""
    if not path.endswith(SUFFIX) and not path.endswith(".json"):
        return False
    try:
        with open(path, encoding="utf-8") as fh:
            return '"craic-session"' in fh.read(400)
    except OSError:
        return False