Validation
CRAIC makes two claims that can be tested rather than asserted:
- that its per-column reliability score identifies the columns that are actually wrong, and
- that masking the flagged columns leaves an alignment more accurate than removing the same number of columns by a trivial rule.
Both are measured by the benchmark harness that ships with the software
(benchmarks/), against two independent sources of ground truth. This page
reports what it found. The numbers are regenerated by the commands given at the
bottom, so they can be checked.
How accuracy is measured
SP (sum-of-pairs) is the fraction of the reference's homologous residue pairs the alignment recovers — recall. TC (total column) is the fraction of the reference's columns reproduced exactly, which is much harsher: one misplaced residue loses the column. Reliability AUC is the ROC area under the curve of the reliability score predicting which columns are actually correct; 0.5 is chance, and this is the number that decides whether the reliability claim is worth anything.
A scope caveat that matters for comparing against published tables: BAliBASE SP
and TC figures in the literature are computed over the reference's core
blocks only — the regions the structural alignment is confident about — using
the official bali_score. An all-column figure is a different and generally
lower quantity. Both are reported below and each is labelled.
Real data: BAliBASE 3
Every family was re-aligned with CRAIC's built-in engine and with MAFFT 7.505,
MUSCLE 5.1, Clustal Omega 1.2.4, ProbCons 1.12 and PRANK v.170427, then scored
two ways: with the official bali_score over core blocks, and with CRAIC's own
all-column metrics (which is what the reliability AUC and the masking controls
need, since they are per-column quantities).
Only families that every aligner completed are included, so the means compare like with like. 163 families qualify, covering RV11 (76 of 76) and RV12 (87 of 88 — one family's dense pairwise posteriors exceed the memory budget, and it is skipped for every engine alike).
RV20, RV30, RV40 and RV50 are not reported. The consistency transformation does not scale to families of 30–450 sequences (see Limitations), so those sets can only ever be partially covered, and a partial set is worse than an absent one: which families drop out is decided by CRAIC's memory ceiling rather than by anything about the aligners, so the survivors are a biased sample. RV11 and RV12 are reported because they are complete.
Figures below are means ± standard error.
Alignment accuracy, official core-block bali_score
These are the numbers comparable with published BAliBASE tables.
| aligner | SP | TC |
|---|---|---|
| MUSCLE 5.1 | 0.856 ± 0.012 | 0.712 ± 0.020 |
| ProbCons 1.12 | 0.850 ± 0.013 | 0.705 ± 0.021 |
| CRAIC built-in | 0.827 ± 0.015 | 0.674 ± 0.022 |
| MAFFT 7.505 | 0.818 ± 0.016 | 0.667 ± 0.022 |
| Clustal Omega 1.2.4 | 0.766 ± 0.018 | 0.593 ± 0.024 |
| PRANK v.170427 | 0.714 ± 0.018 | 0.506 ± 0.023 |
Paired Wilcoxon tests against the built-in engine, on SP:
| aligner | difference | better / worse | p |
|---|---|---|---|
| MUSCLE | +0.029 | 107 / 40 | 3.7 × 10⁻⁹ |
| ProbCons | +0.023 | 99 / 52 | 4.6 × 10⁻⁶ |
| MAFFT | −0.009 | 61 / 92 | 0.12 |
| Clustal Omega | −0.061 | 33 / 128 | 2.3 × 10⁻¹² |
| PRANK | −0.113 | 5 / 157 | 1.1 × 10⁻²⁷ |
Three things in that table are worth drawing out.
ProbCons is there on purpose. CRAIC's built-in engine follows ProbCons's method, so a table reporting it against MAFFT and MUSCLE alone would invite the reader to conclude that posterior-consistency alignment is the weaker approach. ProbCons itself scores 0.023 above the built-in engine, so most of that difference is this implementation rather than the method — the built-in engine estimates its pair-HMM parameters per dataset where ProbCons EM-trains them once, and that is the substantive remaining difference between them.
The built-in engine and MAFFT are indistinguishable here (p = 0.12). It is behind MUSCLE and ProbCons, ahead of Clustal Omega, and the ordering is not uniform across the sets: it is fourth on BB11 and BB12 and third on BBS1, where it beats MAFFT by 0.042. Per-set figures are in the supplementary tables.
PRANK was re-run in 0.5.10
CRAIC 0.5.9 passed PRANK -gaprate=0.025 -gapext=0.5 for every dataset.
PRANK's own defaults for protein are 0.005 / 0.5, so the first version of
this benchmark aligned every family with five times PRANK's gap-opening
rate. The PRANK rows here are re-run with PRANK's own defaults (SP 0.686 →
0.714, TC 0.470 → 0.506); every other aligner's rows are unchanged, and so
are the conclusions drawn from them.
PRANK's position is a property of the benchmark, not a verdict on PRANK. It is phylogeny-aware and declines to over-align — it treats insertions as insertions rather than forcing them into shared columns — while reference benchmarks reward maximal alignment. Scoring it this way largely measures that design choice; it is included because CRAIC wraps it and omitting an engine the workbench offers would be a selection the reader could not see.
Per-set accuracy
| aligner | BB11 SP | BB12 SP | BBS1 SP |
|---|---|---|---|
| n=38 | n=43 | n=82 | |
| MUSCLE | 0.685 | 0.942 | 0.889 |
| ProbCons | 0.670 | 0.941 | 0.886 |
| CRAIC built-in | 0.613 | 0.926 | 0.874 |
| MAFFT | 0.651 | 0.938 | 0.832 |
| Clustal Omega | 0.590 | 0.905 | 0.775 |
| PRANK | 0.490 | 0.851 | 0.745 |
No method reconstructs these alignments
The premise the workbench is built on is measurable, and it is not a matter of degree. Over the same 163 families:
| aligner | families reconstructed exactly (TC = 1.000) |
|---|---|
| MUSCLE | 6 / 163 (3.7%) |
| ProbCons | 6 / 163 (3.7%) |
| CRAIC built-in | 6 / 163 (3.7%) |
| MAFFT | 5 / 163 (3.1%) |
| Clustal Omega | 0 / 163 (0.0%) |
| PRANK | 0 / 163 (0.0%) |
| any method | 10 / 163 (6.1%) |
| all six methods | 0 / 163 (0.0%) |
94% of families are reconstructed incorrectly by every aligner tested, and there is not one family that all six get right. Nor is the best method knowable in advance: MUSCLE produces the best alignment for 88 of the 163 families, leaving 75 won by one of the other five. An oracle allowed to pick the best method per family would reach TC 0.739 against MUSCLE's 0.712 — an improvement unavailable in practice, and one that still leaves a quarter of columns misplaced.
These are capable tools; each recovers the large majority of homologous pairs. The point is that no automatic procedure available today returns an alignment that can be treated as settled, and that the errors differ from method to method — so deciding which parts to rely on remains a judgement a person has to make. That is the case for putting the uncertainty in front of them.
Reliability: does the score find the wrong columns?
This is the claim that matters. The score used here is the consistency component alone (the benchmark was run without the perturbation ensemble, which is an order of magnitude more expensive); the combined score adds perturbation on top of it.
These are CRAIC's own all-column metrics, not the core-block figures above, since the AUC is a per-column quantity and needs every column.
| alignment built by | SP | TC | reliability AUC |
|---|---|---|---|
| MUSCLE | 0.740 | 0.466 | 0.833 ± 0.010 |
| ProbCons | 0.732 | 0.460 | 0.845 ± 0.009 |
| CRAIC built-in | 0.710 | 0.444 | 0.843 ± 0.008 |
| MAFFT | 0.700 | 0.431 | 0.874 ± 0.007 |
| Clustal Omega | 0.665 | 0.399 | 0.876 ± 0.009 |
| PRANK | 0.596 | 0.332 | 0.905 ± 0.007 |
An AUC of 0.83–0.90 means that, picking a correct column and an incorrect one at random, the score ranks them the right way round between five and nine times in ten. It works on alignments CRAIC did not build as well as on its own, which is the property that matters for a curation tool.
This column is not a ranking of the aligners — it is very close to the reverse of one. Across the six engines, reliability AUC and alignment accuracy are almost perfectly anticorrelated (Spearman ρ = −0.94, p = 0.005): the AUC is highest on PRANK's alignments, which are the least accurate, and lowest on MUSCLE's, which are the most. That is expected rather than perverse. A less accurate alignment contains more genuinely wrong columns, so there is more signal for any classifier to separate; a nearly-correct alignment offers few errors to find and every mistake costs proportionally more. The column says how well the score identifies wrong columns given an alignment. It says nothing whatever about which alignment to prefer, and reading it as if it did would invert the answer.
Masking, against controls
Masking raises mean accuracy for any column score correlated with difficulty, so the all-column/kept-column contrast is not evidence of anything. Both controls therefore retain the same number of columns: a matched random mask, and the trivial gap-fraction rule. Only the margin over those is evidence.
| alignment built by | kept columns | matched random | gap-fraction rule | z |
|---|---|---|---|---|
| PRANK | 0.846 | 0.592 | 0.671 | +10.8 |
| Clustal Omega | 0.842 | 0.615 | 0.710 | +9.8 |
| MUSCLE | 0.837 | 0.671 | 0.733 | +8.1 |
| ProbCons | 0.832 | 0.666 | 0.733 | +8.1 |
| CRAIC built-in | 0.829 | 0.645 | 0.728 | +8.0 |
| MAFFT | 0.828 | 0.634 | 0.728 | +9.5 |
A random mask is an unbiased subsample, so it recovers the all-column mean rather than improving on it; its job here is to fix the baseline and supply the spread the z is measured in. The gap-fraction rule is the real competitor: it costs nothing to compute and already reaches 0.71–0.73, which is over half the distance from random to the reliability score. The reliability score reaches 0.83–0.84.
The informative margin is therefore the last ≈ 0.10 (0.10–0.18 across the six engines), not the ≈ 0.19 over random and not what an uncontrolled before-and-after comparison would suggest. Gap fraction is one line of code, so the pair-HMM, the consistency transformation and the perturbation ensemble have to beat it to earn the compute they cost. They do, consistently and across every engine — but by about a third of what the naive comparison would let you claim.
One caveat on the z column: it is computed against the random-mask null only, so it tests the ≈ 0.19 margin rather than the ≈ 0.10 one that matters.
Simulated data
Simulation gives per-residue ground truth, so it can measure things a structural reference cannot — in particular how accuracy and reliability behave as a function of divergence, and whether masking helps the downstream tree.
Sequences are evolved down random binary trees with substitutions, insertions and deletions and among-site rate variation (gamma, α = 0.8): 360 datasets over 8 and 16 taxa, root length 400 nt, six divergence levels (0.1 to 1.5) and three indel rates, 10 replicates per condition.
| aligner | SP | TC | reliability AUC | kept | random | gap rule | z |
|---|---|---|---|---|---|---|---|
| CRAIC built-in | 0.539 | 0.355 | 0.875 | 0.893 | 0.461 | 0.592 | +10.6 |
| MAFFT | 0.608 | 0.379 | 0.893 | 0.843 | 0.533 | 0.548 | +11.3 |
| MUSCLE | 0.663 | 0.387 | 0.858 | 0.884 | 0.581 | 0.669 | +8.5 |
| Clustal Omega | 0.549 | 0.261 | 0.943 | 0.836 | 0.529 | 0.532 | +12.2 |
The reliability AUC is higher here than on BAliBASE (0.86–0.94 against 0.83–0.90), which is what you would expect and is worth saying out loud: the simulated sequences are generated under a model closer to the pair-HMM that computes the score than real protein families are. The BAliBASE figures are the honest ones; these are the upper bound.
Where alignment breaks
Averages hide the interesting part. Accuracy does not decay gracefully with divergence — it falls off a cliff between 0.25 and 0.5, and every method falls off it together:
| divergence | CRAIC | MAFFT | MUSCLE | Clustal Omega |
|---|---|---|---|---|
| 0.1 | 0.982 | 0.965 | 0.964 | 0.860 |
| 0.25 | 0.905 | 0.887 | 0.892 | 0.746 |
| 0.5 | 0.598 | 0.680 | 0.730 | 0.594 |
| 0.75 | 0.366 | 0.509 | 0.590 | 0.473 |
| 1.0 | 0.249 | 0.396 | 0.490 | 0.382 |
| 1.5 | 0.131 | 0.214 | 0.314 | 0.241 |
Two things follow. The built-in engine is competitive at low divergence and falls behind faster than the others as divergence rises, which is the expected signature of an engine without empirically-tuned parameters — another reason to use MAFFT or MUSCLE on hard data. And past a divergence of about 0.5 every aligner is recovering less than three-quarters of the true homologies, which is precisely the regime where an interactive view of which columns are wrong stops being a nicety.
The reliability score, by contrast, keeps working across the whole range (AUC 0.83–0.90 on CRAIC's own alignments) — it does not degrade where the alignment does.
The honest negative result
Masking improves per-column accuracy substantially and consistently. It does not improve neighbour-joining tree recovery. It makes it worse, at every divergence above the trivial, and increasingly so — Robinson–Foulds distance to the true topology, CRAIC alignments, default 0.5 threshold:
| divergence | RF, full alignment | RF, masked |
|---|---|---|
| 0.1 | 0.00 | 0.00 |
| 0.25 | 0.04 | 0.15 |
| 0.5 | 0.16 | 0.47 |
| 0.75 | 0.35 | 0.65 |
| 1.0 | 0.52 | 0.73 |
| 1.5 | 0.65 | 0.84 |
The reason is not mysterious: at high divergence the threshold discards most of the columns, and with them the phylogenetic signal. A more accurate alignment of a quarter of the data is a worse basis for a tree than a sloppier alignment of all of it. This is consistent with Tan et al. (2015), who found that automated alignment filtering frequently worsens single-gene phylogenetic inference.
We report it prominently because it bears directly on how the tool should be used. The reliability signal is for locating and interrogating ambiguous regions — and, if you must act on it, for curating them by hand. It is not an automatic filter to put in front of a tree search, and CRAIC's masking slider should not be treated as one.
Masking residues instead of columns
Masking a column throws away its correctly aligned residues along with the wrong ones, and residue-level filters such as Divvier (Ali et al. 2019) and CLOAK (Wheeler et al. 2026) are reported to lose less signal. CRAIC 0.5.10 can mask individual residues, so the same simulations were used to ask whether masking the residues the reliability score flags does better than masking the columns. For each CRAIC alignment the tree was built from the full alignment; with columns below 0.5 removed; with the lowest-scoring residues masked as missing data — exactly as many as the column mask removed, so both masks remove the same amount; with that many residues masked at random; and with every residue below 0.5 masked. Neighbour-joining skips masked residues pairwise, like gaps; IQ-TREE 2 (JC+G4, the simulating model), run at divergences 0.25 and 0.5, handles them site by site. Datasets in which any mask left too little data for a tree are excluded from every comparison (n below, of 60 per divergence).
| divergence | tree | removed | full | columns | residues, lowest | residues, random | residues < 0.5 | n |
|---|---|---|---|---|---|---|---|---|
| 0.1 | NJ | 2% | 0.00 | 0.01 | 0.01 | 0.01 | 0.01 | 60 |
| 0.25 | NJ | 25% | 0.04 | 0.16 | 0.22 | 0.21 | 0.22 | 60 |
| 0.5 | NJ | 70% | 0.14 | 0.47 | 0.59 | 0.63 | 0.57 | 55 |
| 0.75 | NJ | 89% | 0.28 | 0.62 | 0.76 | 0.84 | 0.68 | 51 |
| 1.0 | NJ | 96% | 0.45 | 0.75 | 0.87 | 0.92 | 0.76 | 37 |
| 1.5 | NJ | 99% | 0.55 | 0.81 | 0.93 | 0.97 | 0.79 | 20 |
| 0.25 | ML | 25% | 0.02 | 0.12 | 0.19 | 0.18 | 0.18 | 60 |
| 0.5 | ML | 70% | 0.12 | 0.47 | 0.57 | 0.61 | 0.53 | 52 |
Masking residues did not help. At the same amount removed, masking the lowest-scoring residues gave worse trees than masking columns (mean Robinson–Foulds difference +0.09 with neighbour-joining, better in 18 datasets and worse in 111, Wilcoxon p = 3 × 10⁻¹⁶; +0.09 with maximum likelihood, 6 better and 40 worse, p = 5 × 10⁻⁷), and only a little better than masking the same number of residues at random (−0.03, p = 0.005, and −0.01, p = 0.19). Masking every residue below 0.5 was also worse than masking columns. The flagged residues are not spread evenly across sequences, so the sequences that are hardest to align are left with the least data. At this threshold the residues the score flags carry phylogenetic signal as well as error. Gentler filtering, removing far fewer residues, has been reported to improve gene trees (Wheeler et al. 2026) and was not tested here.
These trees come from a re-run of the simulation with 0.5.10, whose built-in engine gives slightly different alignments from those behind the tables above (SP within about 0.01), so the full and column-masked values differ a little from the previous table's.
So residue masks are a curation tool — for a residue you have looked at and judged wrong — not a better automatic filter.
What this does not show
- Nothing here validates the perturbation component separately against GUIDANCE, because CRAIC's perturbation score is deliberately not equivalent to it (16 replicates of CRAIC's own engine, not ~100 of the aligner under test). See Limitations.
- Large families are excluded. The consistency-based analysis is O(N³) in sequences and O(N²L²) in memory; the benchmark skips any family whose pairwise posteriors would exceed a stated budget, and skipped families are excluded from every aligner's mean so the comparison stays fair. The number skipped is reported by the harness.
- The built-in aligner is not being advanced as competitive. It is included because a reader is entitled to know how the always-available engine performs, and the answer is: adequately, and behind MAFFT and MUSCLE.
Reproducing this
You need the aligners on your PATH and BAliBASE 3 unpacked. The compiled Rust
core makes the whole thing roughly ten times faster but changes no results.
# 1. simulation arm — 360 datasets
python benchmarks/run_benchmark.py --sim --taxa 8 16 --length 400 --reps 10 \
--divergence 0.1 0.25 0.5 0.75 1.0 1.5 --indel 0.5 2.0 5.0 \
--rate-alpha 0.8 --perturbation --out sim.csv
# 2. real-data arm — BAliBASE 3, the two reference sets reported above.
# --effort max runs the built-in engine with iterative refinement; leaving it
# at the default would measure the engine with a stage switched off.
# --engine-timeout 0 disables the per-aligner cap: PRANK needs several minutes
# on the larger families and a cap records our impatience as its failure.
python benchmarks/run_benchmark.py --protein \
--reference "bb3_release/RV1*/*.xml" --effort max \
--write-alignments aligned_out --max-mem-gb 3 --engine-timeout 0 \
--out balibase.csv
# 3. official core-block scores for the alignments written in step 2
python benchmarks/score_balibase.py aligned_out \
--reference-dir bb3_release --bali-score ./bali_score \
-o balibase_official.csv
Building the official scorer is described in the header of
benchmarks/score_balibase.py: BAliBASE ships its source, but the bundled
libexpat.a is a 2008 Alpha binary that will not link on any current machine, so
it has to be built against the system expat. Verify it by scoring a reference
against itself, which must give exactly 1.000 / 1.000.