[glennhickey]

Prototype for scoring lastz by divergence. lastz scores every pair with its default, HOXD70, whether the two genomes are 0.01 or 1.5 apart. Draft. Per-pair training is on by default (<blast lastzTrain="1">); the table is opt-in.

What it adds

A table of lastz scoring models by distance. Any number of these, in <blast>:

<lastzScoreModel maxDistance="0.08" matrix="A:A A:C ... T:T (16 scores)" gapOpen="744" gapExtend="22"/>

A pair is aligned with the model of smallest maxDistance that covers its distance (the same distance that picks its lastzArguments), and with HOXD70 past the last one. The model goes to lastz as a --scores file. lastz lets the command line override such a file, so each divergence class keeps its thresholds. Given HOXD70 itself as a model, lastz's output is byte-identical to the default.

Per-pair training (<blast lastzTrain>, on by default).

• When: each ancestor's blast step first trains a model with last-train on every pair of genomes it will align with lastz: the ingroup pairs, and each ingroup against each outgroup. Every lastz job for that pair then scores with it.
• Whole genomes: training uses the full sequences. Outgroups after the first are aligned against a trimmed ingroup (whatever earlier outgroups left), so training there would fit leftovers.
• Scale: each model is rescaled to HOXD70's λ, so lastz's thresholds and the chaining's, which are in HOXD70's units, keep their meaning.
• Fallback: if the smaller genome has under 500 kb of sequence, or training fails, the pair falls back to the table, else HOXD70.
• Options: lastzTrainGaps="0" keeps HOXD70's gap costs with the trained matrix; lastzTrainCores sets the threads, up to what Toil allows a job (--maxCores, or a single machine's cores).
• Repeat masking: training uses the inputs' own soft-masking (lastdb -R10, which lastal follows) rather than LAST's default of masking only simple repeats with tantan. Cactus's Red preprocessing provides that masking, the reconstructed ancestors keep it, and it's also what lastz itself doesn't seed in.

Training doesn't apply to KegAlign, minimap2 or FastGA; their alignments keep their default scoring, and KegAlign warns once per ancestor. A lastzScoreModel table with KegAlign still stops with an error.

Evolver mammals

Pooled leaf-pair precision / recall against the truth. bar is stock abPOA unless noted.

lastz scoring	P / R
HOXD70 (shipped)	0.816 / 0.740
table from an evolver last-train ladder (0.013–0.60)	0.829 / 0.743
same, matrices only (HOXD70 gaps)	0.826 / 0.744
same, gaps only (HOXD70 matrix)	0.821 / 0.740
per-pair training	0.830 / 0.744
per-pair training, matrix only	0.829 / 0.745
pecan per-pair models in bar (#2024), lastz HOXD70	0.841 / 0.755
pecan per-pair models + lastz table	0.854 / 0.759
pecan per-pair models + per-pair lastz training	0.856 / 0.760

• Gain: the trained matrices are what help, mostly by penalising mismatches, transitions especially, much harder than HOXD70 at short and moderate distances. The gain adds to bar's.
• Cost: ancestor recall falls 0.4–0.9 points at the higher ancestors, while ancestor precision and base accuracy rise everywhere. Matrices without the trained gaps lose about half as much recall.
• Primates: evolver primates don't move (0.9943/0.9863 to 0.9945/0.9863).
• Ancestor pairs: this is where per-pair models differ most from the table, because the table's distance adds the ancestors' heights. The two ancestors at the root trained to 93.6% identity; the table gave them its most distant model, fit for about 65%.
• Reproducibility: training is reproducible. The first ancestors' blast output was byte-identical across runs.

Cost on real genomes

Whole genomes from the 577-way tree, 4 threads. The mouse database was split into 4 GB volumes to fit a 31 GB machine.

pair (raw tree distance)	total CPU, LAST's default masking	with the inputs' masking (-R10)	model change	two equally good samples differ by
human vs mouse (0.33)	4,494 s	991 s	5.6%	8.3%
human vs fish (1.50)	1,144 s	753 s	21.7%	16.8%
two mouse assemblies (0.0003)	4,148 s	936 s	identity 99.35% → 100.00%	4.7%

"Model change" is the largest relative change among the eight distinct lastz parameters.

• Two mouse assemblies: with interspersed repeats left in, the default masking trained on alignments between repeat copies. The 100% with the inputs' masking is the right answer for two assemblies of one species.
• Evolver: accuracy was the same either way (0.830/0.744).
• For the 577-way: about 4,020 pairs, roughly 800–1,100 CPU-hours (about 0.3% of the 329,000 CPU-hours its lastz took). Training adds about 7 minutes before a mammal ancestor's lastz can start.

Ruled out, as the models moved well beyond that noise:

• Capping the database: 7–76% model change. The sampled query loses its true ortholog, so the model fits paralogs (close pairs) or junk (distant pairs).
• Sparse seeding (-uRY4): 11–66% change.

A smaller query sample stays within the noise but is no longer worth it, since the passes are now cheap.

LAST

With this PR, progressive cactus needs LAST by default. LAST is now built from a submodule (upstream 1651) rather than by downloadPangenomeTools, as of #2034, which is merged.

This uses only stock LAST: 1651 gives the same models. last-train is given --tolerance=1 when it has it (a patch we drafted, not in LAST). That saves a few passes, which now cost little.

Testing

The evolver runs were made on this code on top of the bar work in #2024. It was then rebased onto the current master, which has since gained #2029 (90 Mbp lastz chunks), #2030 (Red -sketch) and pangenome-only changes. The blast code (local_alignment.py, paf.py, common.py) is identical between the two branches.

CI caught one failure: testEvolverPrepareToil. cactus-prepare's Toil workflow runs each cactus-blast with --maxCores 1, and Toil refused the training job's 4 cores. Training now asks for no more cores than Toil allows a job, through a cactus_clamp_cores that works like cactus_clamp_memory. With that fix, the test's cactus-prepare-toil workflow ran to completion locally and passed the test's accuracy check. That run was on #2024's branch, which carries the same blast code. The test's last step, cactus-phast, wasn't run locally, as it needs phast tools that aren't installed here.

Not in this PR

• A table built from real genomes, as a fallback for pairs too small to train. Real-data validation of the defaults needs the 577-way pilot.
• Sharing one database across all of an ingroup's pairs, which would cut training CPU by roughly 2–3×.

🤖 Generated with Claude Code