From b969be14e57a36cc0fc529844c9b440ac0305d56 Mon Sep 17 00:00:00 2001
From: Glenn Hickey <glenn.hickey@gmail.com>
Date: Fri, 2 Oct 2026 09:52:40 -0400
Subject: [PATCH 1/6] blast: optional lastz scoring models, chosen by
 divergence

lastz scored every pair with its default HOXD70, however diverged.
<blast> can now hold any number of <lastzScoreModel maxDistance= matrix=
gapOpen= gapExtend=>.  A pair is aligned with the model of smallest
maxDistance its distance is within -- the distance that picks its
lastzArguments -- written to a scores file and passed as --scores, and
with HOXD70 past them all.  lastz lets the command line override a
scores file, so each divergence class keeps its thresholds.  None are
configured by default, so nothing changes unless a config adds them.
KegAlign does not take them, and a run that asks for both stops rather
than ignoring the models.

HOXD70 given as a model gives byte-identical lastz output to the
default.  On evolver mammals, with bar stock, a table of last-train
models trained on evolver pairs from 0.013 to 0.6 apart, rescaled to
HOXD70's lambda, raised pooled leaf-pair precision from 0.816 to 0.829
and recall from 0.740 to 0.743; ancestor precision rose and recall fell,
by up to 0.9 points.  Evolver primates were unchanged.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/cactus/cactus_progressive_config.xml |  6 ++++
 src/cactus/paf/local_alignment.py        | 37 ++++++++++++++++++++++++
 2 files changed, 43 insertions(+)

diff --git a/src/cactus/cactus_progressive_config.xml b/src/cactus/cactus_progressive_config.xml
index 9cbe53932..92cdb1639 100644
--- a/src/cactus/cactus_progressive_config.xml
+++ b/src/cactus/cactus_progressive_config.xml
@@ -145,6 +145,12 @@
 				five="--step=2 --ambiguous=iupac,100,100 --ydrop=4000 --hspthresh=2400 --gappedthresh=2600 --queryhspbest=300000"
 				default="--step=1 --ambiguous=iupac,100,100 --ydrop=4000 --hspthresh=2200 --gappedthresh=2400 --queryhspbest=300000"
 		/>
+		<!-- lastz scores by divergence, in place of its default HOXD70 (none by default).  Any number of
+		     <lastzScoreModel maxDistance="D" matrix="M" gapOpen="O" gapExtend="E"/>, where M is 16 scores, rows and
+		     columns A, C, G, T, and a gap of length k costs O + k*E, as with lastz's gap penalties.  A pair is
+		     aligned with the model of smallest maxDistance its distance (the one that picks its lastzArguments above)
+		     is within, and with HOXD70 when past them all.  Keep a model on HOXD70's scale (its lambda): the
+		     thresholds above, and the chaining's, are in HOXD70's units.  Not supported with KegAlign. -->
 		<!-- Estimated runtime, in seconds, of one lastz job in each divergence class, for a
 		     pair of full-sized (chunkSize, 90 Mbp) chunks.  These are the p99 of 827k lastz runs from
 		     the VGP 577-way alignment (90 Mbp chunks, queryhspbest 100000, so they predate 300000), halved for the 2x lastz speedup since; Cactus scales
diff --git a/src/cactus/paf/local_alignment.py b/src/cactus/paf/local_alignment.py
index 64e2dd3e7..cf62fcab5 100755
--- a/src/cactus/paf/local_alignment.py
+++ b/src/cactus/paf/local_alignment.py
@@ -39,6 +39,32 @@ def get_divergence_class(distance, params):
     return "default"
 
 
+def get_lastz_score_model(distance, params):
+    """ The <blast><lastzScoreModel> lastz aligns a pair at the given distance with: the one with the
+    smallest maxDistance the distance is within, or None, leaving lastz its default (HOXD70), when there
+    are none or the distance is past them all.  The distance is the one the divergence class is chosen by. """
+    models = params.find("blast").findall("lastzScoreModel")
+    for model in sorted(models, key=lambda m: float(m.attrib["maxDistance"])):
+        if distance <= float(model.attrib["maxDistance"]):
+            return model
+    return None
+
+
+def write_lastz_scores(model, path):
+    """ A <blast><lastzScoreModel> as a lastz scores file: its substitution matrix (rows and columns A, C,
+    G, T) and gap penalties, a gap of length k costing gapOpen + k * gapExtend, as with lastz's --gap.
+    Anything else the file could set is left to the command line. """
+    matrix = [int(s) for s in model.attrib["matrix"].split()]
+    if len(matrix) != 16:
+        raise RuntimeError("lastzScoreModel matrix has {} scores, not 16".format(len(matrix)))
+    with open(path, 'w') as scores_file:
+        scores_file.write('gap_open_penalty = {}\n'.format(int(model.attrib["gapOpen"])))
+        scores_file.write('gap_extend_penalty = {}\n'.format(int(model.attrib["gapExtend"])))
+        scores_file.write('  A C G T\n')
+        for i, base in enumerate('ACGT'):
+            scores_file.write('{} {}\n'.format(base, ' '.join(str(s) for s in matrix[4 * i:4 * i + 4])))
+
+
 def get_lastz_walltime(distance, params, chunk_a, chunk_b):
     """ Estimated seconds for one lastz/kegalign job on this pair of chunks.
 
@@ -85,9 +111,20 @@ def run_lastz(job, name_A, genome_A, name_B, genome_B, distance, params):
     lastz_divergence_node = lastz_params_node.find("kegalignArguments" if gpu else "lastzArguments")
     divergences = params.find("constants").find("divergences")
     lastz_params = lastz_divergence_node.attrib[get_divergence_class(distance, params)]
+    score_model = get_lastz_score_model(distance, params)
+    if score_model is not None:
+        if gpu:
+            raise RuntimeError("<blast><lastzScoreModel> is not supported with KegAlign")
+        # cactus_call runs lastz in work_dir
+        write_lastz_scores(score_model, os.path.join(work_dir, 'lastz.scores'))
+        lastz_params += ' --scores=lastz.scores'
     if not getOptionalAttrib(divergences, 'useDefault', typeFn=bool, default=False):
         logger.info("For distance {} for genomes {}, {} using {} lastz parameters".format(distance, genome_A,
                                                                                           genome_B, lastz_params))
+    if score_model is not None:
+        logger.info("Scoring with the lastzScoreModel for distances up to {}: {}, gap {}+{}k".format(
+            score_model.attrib["maxDistance"], score_model.attrib["matrix"], score_model.attrib["gapOpen"],
+            score_model.attrib["gapExtend"]))
     if gpu:
         lastz_bin = 'run_kegalign'
         suffix_a, suffix_b = '', ''

From 6cfb490aa0e61b10931bf4a9e9a196909754980b Mon Sep 17 00:00:00 2001
From: Glenn Hickey <glenn.hickey@gmail.com>
Date: Fri, 2 Oct 2026 12:14:39 -0400
Subject: [PATCH 2/6] blast: optionally train lastz's scores on each pair it
 aligns

With <blast lastzTrain="1">, 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 -- and every
lastz job for the pair scores with it, rescaled to HOXD70's lambda so
that lastz's thresholds and the chaining's mean what they did.  Training
uses the whole genomes: the ingroup aligned to an outgroup after the
first is only what earlier outgroups left of it.  A pair whose smaller
genome has under 500kb, or whose training fails, falls back to the
lastzScoreModel table, else HOXD70.  lastzTrainGaps="0" keeps HOXD70's
gap penalties with the trained matrix; lastzTrainCores sets the threads.
last-train is given --tolerance=1 when it has it (a proposed LAST patch).

last_scoring.py gets the conversion: last-train's fine-scale model,
which parse_train_file's integer one is too coarse to rescale, put on
HOXD70's scale.  It gives the evolver ladder's models exactly as the
table experiments had them.

On evolver mammals all 18 pairs trained, about 9s each.  Per-pair
models matched or beat the table built from the same evolver pairs:
pooled leaf-pair precision/recall 0.830/0.744 against the table's
0.829/0.743 with stock bar, and 0.856/0.760 against 0.854/0.759 with
pecan's per-pair models (HOXD70 with stock bar: 0.816/0.740).  Where
they differ most is pairs with a reconstructed ancestor: the two
ancestors at the root trained to 92.9% identity, where the table, going
by a distance with both ancestors' heights added, gave them its most
distant model, for ~57%.  Training is reproducible: the first ancestors'
blast output was byte-identical between two runs.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/cactus/cactus_progressive_config.xml |  15 +-
 src/cactus/paf/last_scoring.py           |  62 +++++++++
 src/cactus/paf/local_alignment.py        | 170 +++++++++++++++++++----
 3 files changed, 219 insertions(+), 28 deletions(-)

diff --git a/src/cactus/cactus_progressive_config.xml b/src/cactus/cactus_progressive_config.xml
index 92cdb1639..7f0992654 100644
--- a/src/cactus/cactus_progressive_config.xml
+++ b/src/cactus/cactus_progressive_config.xml
@@ -89,6 +89,13 @@
 	     an outgroup rather than an ingroup, so that its coverage to its descendants is not starved during chaining. Set
 	     to 0 to restore the old behaviour of treating the pinned root as an ingroup -->
 	<!-- leftAlign: move each indel in the lastz alignments along its repeat without changing the alignment's score (paffy left_align). "none", or "canonical" (default): at whichever end makes the repeat read smaller than its reverse complement, which gives every pairwise alignment through a repeat the same placement whatever the pair or strand, so they agree when pinched together. "left" (left on the target) only makes alignments agree when no sequence is aligned on its reverse strand, and makes things worse otherwise; it is here for comparison. lastz only -->
+	<!-- lastzTrain Train lastz's scores on each pair of genomes it aligns with last-train (needs LAST), on the whole
+	     genomes, and align the pair with them, rescaled to HOXD70's units (see lastzScoreModel below).  A pair that
+	     cannot be trained (the smaller genome under 500kb, or last-train failing) falls back to the lastzScoreModel
+	     table, else HOXD70.  Not supported with KegAlign -->
+	<!-- lastzTrainGaps With lastzTrain, use the trained gap penalties; 0 keeps HOXD70's (400 and 30) with the trained
+	     substitution scores -->
+	<!-- lastzTrainCores Cores for each pair's training -->
 	<blast chunkSize="90000000"
 	       bigChunkSize="6000000000"
 		   overlapSize="10000"
@@ -124,6 +131,9 @@
 		   pickIngroupPrimaryAlignmentsSeparatelyToOutgroups="1"
 		   includeRootAsOutgroup="1"
 		   bigGenomeChunkScale="2"
+		   lastzTrain="0"
+		   lastzTrainGaps="1"
+		   lastzTrainCores="4"
 		   >
 
 		<!-- The following are parametrised to produce the same results as the default settings,
@@ -149,8 +159,9 @@
 		     <lastzScoreModel maxDistance="D" matrix="M" gapOpen="O" gapExtend="E"/>, where M is 16 scores, rows and
 		     columns A, C, G, T, and a gap of length k costs O + k*E, as with lastz's gap penalties.  A pair is
 		     aligned with the model of smallest maxDistance its distance (the one that picks its lastzArguments above)
-		     is within, and with HOXD70 when past them all.  Keep a model on HOXD70's scale (its lambda): the
-		     thresholds above, and the chaining's, are in HOXD70's units.  Not supported with KegAlign. -->
+		     is within, and with HOXD70 when past them all; with lastzTrain, a pair's trained model comes first and
+		     this table is its fallback.  Keep a model on HOXD70's scale (its lambda): the thresholds above, and the
+		     chaining's, are in HOXD70's units.  Not supported with KegAlign. -->
 		<!-- Estimated runtime, in seconds, of one lastz job in each divergence class, for a
 		     pair of full-sized (chunkSize, 90 Mbp) chunks.  These are the p99 of 827k lastz runs from
 		     the VGP 577-way alignment (90 Mbp chunks, queryhspbest 100000, so they predate 300000), halved for the 2x lastz speedup since; Cactus scales
diff --git a/src/cactus/paf/last_scoring.py b/src/cactus/paf/last_scoring.py
index d006cf190..e63bf8758 100644
--- a/src/cactus/paf/last_scoring.py
+++ b/src/cactus/paf/last_scoring.py
@@ -14,6 +14,7 @@
 from sonLib.bioio import newickTreeParser
 from toil.realtimeLogger import RealtimeLogger
 import os
+import re
 import math
 from cactus.shared.common import cactus_call, getOptionalAttrib
 from cactus.shared.common import cactus_clamp_memory
@@ -192,6 +193,67 @@ def minigraph_wfa_penalties(score_dict):
             match, -mismatch, score_dict['GAP-OPEN'], score_dict['GAP-EXTEND']))
     return max(1, round(x / e)), round(o / e), 1
 
+# lastz's default scores, HOXD70, and the base frequencies its scale is measured with.  A model trained for
+# lastz is put on the same scale (lambda): lastz's thresholds, and the chaining's, are in HOXD70's units.
+HOXD70 = [[91, -114, -31, -123], [-114, 100, -125, -31], [-31, -125, 100, -114], [-123, -31, -114, 91]]
+HOXD70_FREQS = [0.26585, 0.23415, 0.23415, 0.26585]
+
+def implied_lambda(matrix, freqs):
+    """ the lambda with sum_ij f_i f_j exp(lambda s_ij) = 1: the scale of a score matrix, in nats per unit """
+    def excess(lam):
+        return sum(freqs[i] * freqs[j] * math.exp(lam * matrix[i][j]) for i in range(4) for j in range(4)) - 1
+    # last-train's matrices are at ~80 units per nat, its integer ones at ~4: lambda ~0.01 to ~0.25
+    lo, hi = 1e-6, 1.0
+    if not (excess(lo) < 0 < excess(hi)):
+        raise RuntimeError('score matrix {} has no positive lambda'.format(matrix))
+    for _ in range(200):
+        mid = (lo + hi) / 2
+        if excess(mid) < 0:
+            lo = mid
+        else:
+            hi = mid
+    return (lo + hi) / 2
+
+def parse_train_model(train_file_path):
+    """ last-train's final model at the scale it trains at (parse_train_file reads the integer one it ends
+    with, whose scores are too coarse to rescale): matrix rows and columns A, C, G, T, the gap existence and
+    extension costs, the query's base frequencies and the substitution identity.  The training scale differs
+    from one run to the next (last-train raises it 10% at a time when rounding gets too coarse), so the
+    numbers are only comparable between runs once rescaled. """
+    with open(train_file_path) as train_file:
+        txt = train_file.read()
+    header = '# score matrix (query letters = columns, reference letters = rows):'
+    blocks = txt.split(header)
+    if len(blocks) < 3 or '#last -t' not in txt:
+        raise RuntimeError('{} is not a finished last-train run'.format(train_file_path))
+    rows = {}
+    for line in blocks[-2].splitlines()[2:6]:
+        toks = line.split()
+        rows[toks[1]] = [int(x) for x in toks[2:6]]
+    def last_value(key):
+        values = re.findall(r'# {}: ([-\d.eE+]+)'.format(key), txt)
+        if not values:
+            raise RuntimeError('no {} in {}'.format(key, train_file_path))
+        return values[-1]
+    return {'matrix': [rows[b] for b in 'ACGT'],
+            'open': int(last_value('delExistCost')), 'extend': int(last_value('delExtendCost')),
+            'freqs': [float(x) / 100 for x in re.findall(r'# qry letter %: (.*)', txt)[-1].split()],
+            'identity': float(last_value('substitution percent identity'))}
+
+def lastz_scores_from_train(train_file_path, trained_gaps=True):
+    """ the model last-train fitted, as lastz scores on HOXD70's scale: a dict of a <blast><lastzScoreModel>'s
+    matrix, gapOpen and gapExtend (see local_alignment.write_lastz_scores), plus the trained identity.  Without
+    trained_gaps the gap penalties are HOXD70's, 400 and 30, as lastz has them by default. """
+    model = parse_train_model(train_file_path)
+    scale = implied_lambda(model['matrix'], model['freqs']) / implied_lambda(HOXD70, HOXD70_FREQS)
+    matrix = [int(round(s * scale)) for row in model['matrix'] for s in row]
+    if any(matrix[5 * i] <= 0 for i in range(4)):
+        raise RuntimeError('trained matrix {} has a non-positive match score'.format(matrix))
+    return {'matrix': ' '.join(str(s) for s in matrix),
+            'gapOpen': str(int(round(model['open'] * scale))) if trained_gaps else '400',
+            'gapExtend': str(max(1, int(round(model['extend'] * scale)))) if trained_gaps else '30',
+            'identity': model['identity']}
+
 def last_train(job, config, seq_order, seq_id_map, ref_name=None):
     """ run last_train on a pair of fasta files, using the first as the database.
 
diff --git a/src/cactus/paf/local_alignment.py b/src/cactus/paf/local_alignment.py
index cf62fcab5..2adefc705 100755
--- a/src/cactus/paf/local_alignment.py
+++ b/src/cactus/paf/local_alignment.py
@@ -20,6 +20,7 @@
 import copy
 from Bio import SeqIO
 from cactus.paf.paf import get_event_pairs, get_leaves, get_node, get_distances
+from cactus.paf.last_scoring import lastz_scores_from_train
 from cactus.shared.common import cactus_call, getOptionalAttrib, zip_gz, cactus_walltime
 from cactus.preprocessor.checkUniqueHeaders import sanitize_fasta_headers
 from cactus.preprocessor.unmasking import unmask_contigs_all
@@ -51,20 +52,96 @@ def get_lastz_score_model(distance, params):
 
 
 def write_lastz_scores(model, path):
-    """ A <blast><lastzScoreModel> as a lastz scores file: its substitution matrix (rows and columns A, C,
-    G, T) and gap penalties, a gap of length k costing gapOpen + k * gapExtend, as with lastz's --gap.
-    Anything else the file could set is left to the command line. """
-    matrix = [int(s) for s in model.attrib["matrix"].split()]
+    """ A lastz scoring model -- a <blast><lastzScoreModel>'s attributes, or a trained model -- as a lastz
+    scores file: its substitution matrix (rows and columns A, C, G, T) and gap penalties, a gap of length k
+    costing gapOpen + k * gapExtend, as with lastz's --gap.  Anything else the file could set is left to the
+    command line. """
+    matrix = [int(s) for s in model["matrix"].split()]
     if len(matrix) != 16:
-        raise RuntimeError("lastzScoreModel matrix has {} scores, not 16".format(len(matrix)))
+        raise RuntimeError("lastz scoring matrix has {} scores, not 16".format(len(matrix)))
     with open(path, 'w') as scores_file:
-        scores_file.write('gap_open_penalty = {}\n'.format(int(model.attrib["gapOpen"])))
-        scores_file.write('gap_extend_penalty = {}\n'.format(int(model.attrib["gapExtend"])))
+        scores_file.write('gap_open_penalty = {}\n'.format(int(model["gapOpen"])))
+        scores_file.write('gap_extend_penalty = {}\n'.format(int(model["gapExtend"])))
         scores_file.write('  A C G T\n')
         for i, base in enumerate('ACGT'):
             scores_file.write('{} {}\n'.format(base, ' '.join(str(s) for s in matrix[4 * i:4 * i + 4])))
 
 
+def lastz_train_enabled(params):
+    """ Whether to train lastz's scores on each pair of genomes it aligns: <blast lastzTrain> """
+    lastz_params_node = params.find("blast")
+    if not getOptionalAttrib(lastz_params_node, 'lastzTrain', typeFn=bool, default=False):
+        return False
+    if lastz_params_node.attrib["mapper"] != "lastz":
+        logger.warning("<blast lastzTrain> is ignored: it is for lastz, and the mapper is {}".format(
+            lastz_params_node.attrib["mapper"]))
+        return False
+    if getOptionalAttrib(lastz_params_node, 'gpu', typeFn=int, default=0):
+        raise RuntimeError("<blast lastzTrain> is not supported with KegAlign")
+    return True
+
+
+def lastz_pair_key(event_a, event_b):
+    """ The key of a pair of genomes' trained lastz model, whichever way round they are aligned """
+    return ' '.join(sorted([event_a, event_b]))
+
+
+def last_train_has_tolerance():
+    """ Whether last-train has --tolerance, which stops training once no score moves by more than a given
+    percentage, rather than waiting for the scores to repeat exactly (a proposed patch to LAST) """
+    try:
+        return '--tolerance' in cactus_call(parameters=['last-train', '--help'], check_output=True)
+    except Exception:
+        return False
+
+
+def train_lastz_score_model(job, event_a, genome_a, event_b, genome_b, params):
+    """ A scoring model for lastz to align a pair of genomes with, trained on the pair with last-train, as
+    lastz_scores_from_train makes it -- or None, leaving the pair to the lastzScoreModel table or HOXD70,
+    when the smaller genome has under 500kb to train on or training fails. """
+    # lastdb is most of the work, so the smaller genome is the database, judged by file size, a gzipped fasta
+    # (as an ancestor's is) holding about four times its size
+    genomes = []
+    for event, genome in (event_a, genome_a), (event_b, genome_b):
+        with job.fileStore.readGlobalFileStream(genome) as stream:
+            gzipped = stream.read(2) == b'\x1f\x8b'
+        genomes.append((genome.size * (4 if gzipped else 1), event, genome, gzipped))
+    genomes.sort(key=lambda g: (g[0], g[1]))
+    db_event, query_event = genomes[0][1], genomes[1][1]
+    work_dir = job.fileStore.getLocalTempDir()
+    local_fastas = []
+    for name, (_, _, genome, gzipped) in zip(('db', 'query'), genomes):
+        # last-train goes by the name to tell a gzipped fasta
+        local_fastas.append(name + ('.fa.gz' if gzipped else '.fa'))
+        job.fileStore.readGlobalFile(genome, os.path.join(work_dir, local_fastas[-1]))
+    trained_gaps = getOptionalAttrib(params.find("blast"), 'lastzTrainGaps', typeFn=bool, default=True)
+
+    # lastdb -c and -C2 are LAST's recipe for training on genomes: -c keeps masked letters out of the seeds, and
+    # -C2 drops gapless alignments nested in two others
+    train_cmd = ['last-train', '--revsym', '--matsym', '--gapsym', '-C2', '-P', str(job.cores)]
+    if last_train_has_tolerance():
+        train_cmd.append('--tolerance=1')
+    train_cmd += ['db', local_fastas[1]]
+    train_file = os.path.join(work_dir, 'pair.train')
+    try:
+        cactus_call(parameters=['lastdb', '-c', '-P', str(job.cores), 'db', local_fastas[0]], work_dir=work_dir)
+        with open(os.path.join(work_dir, 'db.prj')) as prj:
+            db_letters = int(next(line for line in prj if line.startswith('numofletters=')).split('=')[1])
+        if db_letters < 500000:
+            RealtimeLogger.warning("Not training lastz scores for {} and {}: {} has {} bases, too few to train on".format(
+                event_a, event_b, db_event, db_letters))
+            return None
+        cactus_call(parameters=train_cmd, outfile=train_file, work_dir=work_dir)
+        model = lastz_scores_from_train(train_file, trained_gaps)
+    except Exception as e:
+        RealtimeLogger.warning("Training lastz scores for {} and {} failed, so the pair falls back to the "
+                               "lastzScoreModel table or HOXD70: {}".format(event_a, event_b, e))
+        return None
+    RealtimeLogger.info("Trained lastz scores for {} and {} ({} against {}): identity {:.1f}%, matrix {}, gap {}+{}k".format(
+        event_a, event_b, query_event, db_event, model['identity'], model['matrix'], model['gapOpen'], model['gapExtend']))
+    return model
+
+
 def get_lastz_walltime(distance, params, chunk_a, chunk_b):
     """ Estimated seconds for one lastz/kegalign job on this pair of chunks.
 
@@ -93,7 +170,8 @@ def get_lastz_walltime(distance, params, chunk_a, chunk_b):
     return base
 
 
-def run_lastz(job, name_A, genome_A, name_B, genome_B, distance, params):
+def run_lastz(job, name_A, genome_A, name_B, genome_B, distance, params, trained_score_model=None):
+    """ trained_score_model is the pair's model from train_lastz_score_model, if <blast lastzTrain> made one """
     # Create a local temporary file to put the alignments in.
     work_dir = job.fileStore.getLocalTempDir()
     alignment_file = os.path.join(work_dir, '{}_{}.paf'.format(name_A, name_B))
@@ -111,10 +189,16 @@ def run_lastz(job, name_A, genome_A, name_B, genome_B, distance, params):
     lastz_divergence_node = lastz_params_node.find("kegalignArguments" if gpu else "lastzArguments")
     divergences = params.find("constants").find("divergences")
     lastz_params = lastz_divergence_node.attrib[get_divergence_class(distance, params)]
-    score_model = get_lastz_score_model(distance, params)
+    # the pair's trained model if it has one, else the table's for its distance, else lastz's own HOXD70
+    score_model, score_model_source = trained_score_model, 'trained on this pair'
+    if score_model is None:
+        table_model = get_lastz_score_model(distance, params)
+        if table_model is not None:
+            score_model = table_model.attrib
+            score_model_source = 'the lastzScoreModel for distances up to {}'.format(table_model.attrib["maxDistance"])
     if score_model is not None:
         if gpu:
-            raise RuntimeError("<blast><lastzScoreModel> is not supported with KegAlign")
+            raise RuntimeError("lastz scoring models are not supported with KegAlign")
         # cactus_call runs lastz in work_dir
         write_lastz_scores(score_model, os.path.join(work_dir, 'lastz.scores'))
         lastz_params += ' --scores=lastz.scores'
@@ -122,9 +206,8 @@ def run_lastz(job, name_A, genome_A, name_B, genome_B, distance, params):
         logger.info("For distance {} for genomes {}, {} using {} lastz parameters".format(distance, genome_A,
                                                                                           genome_B, lastz_params))
     if score_model is not None:
-        logger.info("Scoring with the lastzScoreModel for distances up to {}: {}, gap {}+{}k".format(
-            score_model.attrib["maxDistance"], score_model.attrib["matrix"], score_model.attrib["gapOpen"],
-            score_model.attrib["gapExtend"]))
+        logger.info("Scoring with {}: {}, gap {}+{}k".format(score_model_source, score_model["matrix"],
+                                                            score_model["gapOpen"], score_model["gapExtend"]))
     if gpu:
         lastz_bin = 'run_kegalign'
         suffix_a, suffix_b = '', ''
@@ -529,7 +612,8 @@ def chunked_alignment_walltime(genome_a, genome_b):
     return cactus_walltime(CHUNK_SECS_PER_GB * seq_size / 1e9, io_bytes=2 * seq_size)
 
 
-def make_chunked_alignments(job, event_a, genome_a, event_b, genome_b, distance, params):
+def make_chunked_alignments(job, event_a, genome_a, event_b, genome_b, distance, params, score_model=None):
+    """ score_model is the pair's trained lastz model, if <blast lastzTrain> made one """
     lastz_params_node = params.find("blast")
     gpu = getOptionalAttrib(lastz_params_node, 'gpu', typeFn=int, default=0)
     fastga = getOptionalAttrib(lastz_params_node, 'mapper', typeFn=str) == 'fastga'        
@@ -562,6 +646,7 @@ def make_chunks(genome):
         for j, chunk_b in enumerate(chunks_b):
             mappers = { "lastz":run_lastz, "minimap2":run_minimap2, "fastga":run_fastga}
             mappingFn = mappers[params.find("blast").attrib["mapper"]]
+            mapper_args = {'trained_score_model': score_model} if score_model is not None and mappingFn is run_lastz else {}
             memory = lastz_memory if lastz_memory else max(200000000, 15*(chunk_a.size+chunk_b.size))
             chunked_alignment_files.append(job.addChildJobFn(mappingFn, '{}_{}'.format(event_a, i), chunk_a,
                                                              '{}_{}'.format(event_b, j), chunk_b, distance, params,
@@ -569,7 +654,8 @@ def make_chunks(genome):
                                                              disk=max(4*(chunk_a.size+chunk_b.size), memory),
                                                              memory=cactus_clamp_memory(memory),
                                                              accelerators=accelerators,
-                                                             walltime=cactus_walltime(get_lastz_walltime(distance, params, chunk_a, chunk_b))).rv())
+                                                             walltime=cactus_walltime(get_lastz_walltime(distance, params, chunk_a, chunk_b)),
+                                                             **mapper_args).rv())
 
     dechunk_batch_size = getOptionalAttrib(lastz_params_node, 'dechunkBatchSize', typeFn=int, default=1e9)
     # combine_chunks either dechunks the chunk pafs itself or fans them out in batches, so it never
@@ -591,10 +677,13 @@ def invert_alignments(job, alignment_file):
     return job.fileStore.writeGlobalFile(inverted_alignment_file_local)
 
 
-def make_ingroup_to_outgroup_alignments_0(job, ingroup_event, outgroup_events, event_names_to_sequences, distances, params):
+def make_ingroup_to_outgroup_alignments_0(job, ingroup_event, outgroup_events, event_names_to_sequences, distances, params,
+                                          score_models=None):
+    """ score_models holds the trained lastz models, by lastz_pair_key, if <blast lastzTrain> made them """
     # Generate the alignments fle
     alignment_file = job.addChildJobFn(make_ingroup_to_outgroup_alignments_1, ingroup_event, outgroup_events,
-                                            event_names_to_sequences, distances, params, walltime=cactus_walltime()).rv()
+                                            event_names_to_sequences, distances, params, walltime=cactus_walltime(),
+                                            score_models=score_models).rv()
 
     # Invert the final alignment so that the query is the outgroup and the target is the ingroup.
     # paffy invert holds the alignment and its inverse on local disk at once. The alignment's size
@@ -612,7 +701,8 @@ def make_ingroup_to_outgroup_alignments_0(job, ingroup_event, outgroup_events, e
                                                          io_bytes=4 * alignment_seq_size)).rv()
 
 
-def make_ingroup_to_outgroup_alignments_1(job, ingroup_event, outgroup_events, event_names_to_sequences, distances, params):
+def make_ingroup_to_outgroup_alignments_1(job, ingroup_event, outgroup_events, event_names_to_sequences, distances, params,
+                                          score_models=None):
     #  a job should never set its own follow-on, so we hang everything off root_job here to encapsulate
     root_job = Job(walltime=cactus_walltime())
     job.addChild(root_job)
@@ -620,6 +710,8 @@ def make_ingroup_to_outgroup_alignments_1(job, ingroup_event, outgroup_events, e
     #  align ingroup to first outgroup to produce paf alignments
     outgroup = outgroup_events[0] # The first outgroup
     logger.info("Building alignment between ingroup event: {} and outgroup event: {}".format(ingroup_event.iD, outgroup.iD))
+    # the model trained on the whole ingroup, though only what earlier outgroups left of it is aligned here
+    score_model = score_models.get(lastz_pair_key(ingroup_event.iD, outgroup.iD)) if score_models else None
     alignment = root_job.addChildJobFn(make_chunked_alignments,
                                        outgroup.iD, event_names_to_sequences[outgroup.iD],
                                        ingroup_event.iD, event_names_to_sequences[ingroup_event.iD], distances[ingroup_event, outgroup], params,
@@ -627,14 +719,15 @@ def make_ingroup_to_outgroup_alignments_1(job, ingroup_event, outgroup_events, e
                                                                        event_names_to_sequences[ingroup_event.iD], params),
                                        disk=4*(event_names_to_sequences[ingroup_event.iD].size+event_names_to_sequences[outgroup.iD].size),
                                        walltime=chunked_alignment_walltime(event_names_to_sequences[outgroup.iD],
-                                                                           event_names_to_sequences[ingroup_event.iD])).rv()
+                                                                           event_names_to_sequences[ingroup_event.iD]),
+                                       score_model=score_model).rv()
 
     #  post process the alignments and recursively generate alignments to remaining outgroups
     # Memory: paffy to_bed creates SequenceCountArray (2 bytes per base) for query (ingroup) sequences
     ingroup_size = event_names_to_sequences[ingroup_event.iD].size
     outgroup_size = event_names_to_sequences[outgroup.iD].size
     return root_job.addFollowOnJobFn(make_ingroup_to_outgroup_alignments_2, alignment, ingroup_event, outgroup_events[1:],
-                                     event_names_to_sequences, distances, params,
+                                     event_names_to_sequences, distances, params, score_models=score_models,
                                      disk=4*(ingroup_size + outgroup_size),
                                      memory=cactus_clamp_memory(4*ingroup_size + 2*outgroup_size),
                                      # _2 runs paffy to_bed over the ingroup-to-outgroup paf and then
@@ -646,7 +739,7 @@ def make_ingroup_to_outgroup_alignments_1(job, ingroup_event, outgroup_events, e
 
 
 def make_ingroup_to_outgroup_alignments_2(job, alignments, ingroup_event, outgroup_events,
-                                          event_names_to_sequences, distances, params):
+                                          event_names_to_sequences, distances, params, score_models=None):
     
     # a job should never set its own follow-on, so we hang everything off root_job here to encapsulate
     root_job = Job(walltime=cactus_walltime())
@@ -681,7 +774,8 @@ def make_ingroup_to_outgroup_alignments_2(job, alignments, ingroup_event, outgro
 
     # recursively make alignments with the remaining outgroups
     alignments2 = root_job.addChildJobFn(make_ingroup_to_outgroup_alignments_1, ingroup_event, outgroup_events,
-                                         event_names_to_sequences, distances, params, walltime=cactus_walltime()).rv()
+                                         event_names_to_sequences, distances, params, walltime=cactus_walltime(),
+                                         score_models=score_models).rv()
 
     # _3 with has_resources unset only re-schedules itself once the two pafs have resolved, so it
     # is coordination; the instance that does the merging is sized below
@@ -1018,10 +1112,31 @@ def make_paf_alignments(job, event_tree_string, event_names_to_sequences, ancest
         unmask_job = root_job.addChildJobFn(unmask_contigs_all, input_sequence_map, ingroups, params, walltime=cactus_walltime())
         for i,ingroup in enumerate(ingroups):
             event_names_to_sequences[ingroup] = unmask_job.rv(i)
+
+    # train lastz's scores on each pair of genomes it will align (see <blast lastzTrain>), on the whole genomes
+    # (the trimmed ingroups aligned to later outgroups are what earlier ones left) and alongside the unmasking,
+    # which only changes case: last-train masks repeats its own way
+    score_models = {}
+    if lastz_train_enabled(params):
+        train_cores = getOptionalAttrib(lastz_params_node, 'lastzTrainCores', typeFn=int, default=4)
+        pairs = [(ingroup, ingroup2) for ingroup, ingroup2, _ in get_event_pairs(ancestor_event, ingroup_events)] + \
+                [(ingroup, outgroup) for ingroup in ingroup_events for outgroup in outgroup_events]
+        for event_a, event_b in pairs:
+            genome_a, genome_b = input_sequence_map[event_a.iD], input_sequence_map[event_b.iD]
+            # lastdb takes the smaller genome, and its index is several times the sequence (more for a gzipped one)
+            db_size = min(genome_a.size, genome_b.size)
+            score_models[lastz_pair_key(event_a.iD, event_b.iD)] = root_job.addChildJobFn(
+                train_lastz_score_model, event_a.iD, genome_a, event_b.iD, genome_b, params,
+                cores=train_cores,
+                memory=cactus_clamp_memory(int(2e9 + 20 * db_size)),
+                disk=int(2e9 + 25 * db_size + 2 * max(genome_a.size, genome_b.size)),
+                walltime=cactus_walltime(1800 + 7200 * db_size / 1e9, io_bytes=genome_a.size + genome_b.size)).rv()
+
+    if unmask_job is not None or score_models:
         new_root_job = Job(walltime=cactus_walltime())
         root_job.addFollowOn(new_root_job)
         root_job = new_root_job
-    
+
     # for each pair of ingroups make alignments
     ingroup_alignments = []
     # for better logs
@@ -1039,7 +1154,8 @@ def make_paf_alignments(job, event_tree_string, event_names_to_sequences, ancest
                                                                                          input_sequence_map[ingroup2.iD], params),
                                                          disk=2*total_sequence_size,
                                                          walltime=chunked_alignment_walltime(input_sequence_map[ingroup.iD],
-                                                                                             input_sequence_map[ingroup2.iD])).rv())
+                                                                                             input_sequence_map[ingroup2.iD]),
+                                                         score_model=score_models.get(lastz_pair_key(ingroup.iD, ingroup2.iD))).rv())
         ingroup_alignment_names.append('{}-{}_vs_{}'.format(ancestor_event_string, ingroup.iD, ingroup2.iD))
 
     distances = get_distances(event_tree)  # Distances between all pairs of nodes
@@ -1052,7 +1168,8 @@ def make_paf_alignments(job, event_tree_string, event_names_to_sequences, ancest
     # for each ingroup make alignments to the outgroups
     if int(params.find("blast").attrib["trimIngroups"]):  # Trim the ingroup sequences
         outgroup_alignments = [root_job.addChildJobFn(make_ingroup_to_outgroup_alignments_0, ingroup, outgroup_events,
-                                                      dict(event_names_to_sequences), distances, params, walltime=cactus_walltime()).rv()
+                                                      dict(event_names_to_sequences), distances, params, walltime=cactus_walltime(),
+                                                      score_models=score_models).rv()
                                 for ingroup in ingroup_events] if len(outgroup_events) > 0 else []
     else:
         outgroup_alignments = [root_job.addChildJobFn(make_chunked_alignments,
@@ -1064,7 +1181,8 @@ def make_paf_alignments(job, event_tree_string, event_names_to_sequences, ancest
                                                                                       input_sequence_map[outgroup.iD], params),
                                                       disk=2*total_sequence_size,
                                                       walltime=chunked_alignment_walltime(input_sequence_map[ingroup.iD],
-                                                                                          input_sequence_map[outgroup.iD])).rv()
+                                                                                          input_sequence_map[outgroup.iD]),
+                                                      score_model=score_models.get(lastz_pair_key(ingroup.iD, outgroup.iD))).rv()
                                for ingroup in ingroup_events for outgroup in outgroup_events]
     # for better logs
     outgroup_alignment_names = ['{}-og_{}'.format(ancestor_event_string, i) for i in range(len(outgroup_alignments))]

From a933213fb949efaf168a7aea437ae7b0a58f17c8 Mon Sep 17 00:00:00 2001
From: Glenn Hickey <glenn.hickey@gmail.com>
Date: Fri, 2 Oct 2026 16:46:05 -0400
Subject: [PATCH 3/6] blast: size lastz training by bases, and log the identity
 it fitted

The training job's memory, disk and walltime were asked for by file size,
which for a gzipped fasta (as an ancestor's is) is about a quarter of its
bases, and for a plain one counted the bases four times over in memory.
They are now by estimated bases, a gzipped fasta counting four per byte:
lastdb's database is about 4.6 bytes per base (6.3 GB for a 1.4 Gb fish
genome) and lastdb, and lastal in training, hold it all.  The same
estimate picks the smaller genome as the database.  On evolver mammals
the blast output was byte-identical to before.

The identity logged with a trained model was the last one last-train
prints, which is the one its final integer scores imply: their rounding
moved it from the fitted 69.9% to 76.2% for human against mouse.  It is
now the last pass's, from the alignments it counted.  The scores lastz
gets are unchanged.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/cactus/paf/last_scoring.py    | 20 +++++++++++---------
 src/cactus/paf/local_alignment.py | 31 +++++++++++++++++++++----------
 2 files changed, 32 insertions(+), 19 deletions(-)

diff --git a/src/cactus/paf/last_scoring.py b/src/cactus/paf/last_scoring.py
index e63bf8758..1e810944d 100644
--- a/src/cactus/paf/last_scoring.py
+++ b/src/cactus/paf/last_scoring.py
@@ -217,9 +217,9 @@ def excess(lam):
 def parse_train_model(train_file_path):
     """ last-train's final model at the scale it trains at (parse_train_file reads the integer one it ends
     with, whose scores are too coarse to rescale): matrix rows and columns A, C, G, T, the gap existence and
-    extension costs, the query's base frequencies and the substitution identity.  The training scale differs
-    from one run to the next (last-train raises it 10% at a time when rounding gets too coarse), so the
-    numbers are only comparable between runs once rescaled. """
+    extension costs, the query's base frequencies and the identity of the last pass's alignments.  The
+    training scale differs from one run to the next (last-train raises it 10% at a time when rounding gets
+    too coarse), so the numbers are only comparable between runs once rescaled. """
     with open(train_file_path) as train_file:
         txt = train_file.read()
     header = '# score matrix (query letters = columns, reference letters = rows):'
@@ -230,15 +230,17 @@ def parse_train_model(train_file_path):
     for line in blocks[-2].splitlines()[2:6]:
         toks = line.split()
         rows[toks[1]] = [int(x) for x in toks[2:6]]
-    def last_value(key):
-        values = re.findall(r'# {}: ([-\d.eE+]+)'.format(key), txt)
-        if not values:
+    def values(key):
+        found = re.findall(r'# {}: ([-\d.eE+]+)'.format(key), txt)
+        if not found:
             raise RuntimeError('no {} in {}'.format(key, train_file_path))
-        return values[-1]
+        return found
+    # the last identity printed is the one the final integer scores imply, whose rounding can move it by
+    # several points; the one before is the last pass's, from the alignments it counted
     return {'matrix': [rows[b] for b in 'ACGT'],
-            'open': int(last_value('delExistCost')), 'extend': int(last_value('delExtendCost')),
+            'open': int(values('delExistCost')[-1]), 'extend': int(values('delExtendCost')[-1]),
             'freqs': [float(x) / 100 for x in re.findall(r'# qry letter %: (.*)', txt)[-1].split()],
-            'identity': float(last_value('substitution percent identity'))}
+            'identity': float(values('substitution percent identity')[-2])}
 
 def lastz_scores_from_train(train_file_path, trained_gaps=True):
     """ the model last-train fitted, as lastz scores on HOXD70's scale: a dict of a <blast><lastzScoreModel>'s
diff --git a/src/cactus/paf/local_alignment.py b/src/cactus/paf/local_alignment.py
index 2adefc705..2f7c4e869 100755
--- a/src/cactus/paf/local_alignment.py
+++ b/src/cactus/paf/local_alignment.py
@@ -81,6 +81,16 @@ def lastz_train_enabled(params):
     return True
 
 
+def fasta_is_gzipped(job, genome):
+    with job.fileStore.readGlobalFileStream(genome) as stream:
+        return stream.read(2) == b'\x1f\x8b'
+
+
+def fasta_letters_estimate(job, genome):
+    """ About how many bases a fasta file holds: its size, or four times that when it is gzipped """
+    return genome.size * (4 if fasta_is_gzipped(job, genome) else 1)
+
+
 def lastz_pair_key(event_a, event_b):
     """ The key of a pair of genomes' trained lastz model, whichever way round they are aligned """
     return ' '.join(sorted([event_a, event_b]))
@@ -99,13 +109,10 @@ def train_lastz_score_model(job, event_a, genome_a, event_b, genome_b, params):
     """ A scoring model for lastz to align a pair of genomes with, trained on the pair with last-train, as
     lastz_scores_from_train makes it -- or None, leaving the pair to the lastzScoreModel table or HOXD70,
     when the smaller genome has under 500kb to train on or training fails. """
-    # lastdb is most of the work, so the smaller genome is the database, judged by file size, a gzipped fasta
-    # (as an ancestor's is) holding about four times its size
+    # lastdb is most of the work, so the smaller genome is the database (an ancestor's fasta is gzipped)
     genomes = []
     for event, genome in (event_a, genome_a), (event_b, genome_b):
-        with job.fileStore.readGlobalFileStream(genome) as stream:
-            gzipped = stream.read(2) == b'\x1f\x8b'
-        genomes.append((genome.size * (4 if gzipped else 1), event, genome, gzipped))
+        genomes.append((fasta_letters_estimate(job, genome), event, genome, fasta_is_gzipped(job, genome)))
     genomes.sort(key=lambda g: (g[0], g[1]))
     db_event, query_event = genomes[0][1], genomes[1][1]
     work_dir = job.fileStore.getLocalTempDir()
@@ -1121,16 +1128,20 @@ def make_paf_alignments(job, event_tree_string, event_names_to_sequences, ancest
         train_cores = getOptionalAttrib(lastz_params_node, 'lastzTrainCores', typeFn=int, default=4)
         pairs = [(ingroup, ingroup2) for ingroup, ingroup2, _ in get_event_pairs(ancestor_event, ingroup_events)] + \
                 [(ingroup, outgroup) for ingroup in ingroup_events for outgroup in outgroup_events]
+        letters = {event.iD: fasta_letters_estimate(job, input_sequence_map[event.iD])
+                   for pair in pairs for event in pair}
         for event_a, event_b in pairs:
             genome_a, genome_b = input_sequence_map[event_a.iD], input_sequence_map[event_b.iD]
-            # lastdb takes the smaller genome, and its index is several times the sequence (more for a gzipped one)
-            db_size = min(genome_a.size, genome_b.size)
+            # lastdb takes the smaller genome.  Its database is about 4.6 bytes per base (6.3 GB for a 1.4 Gb fish
+            # genome), and lastdb, and lastal in training, hold it all.  With 4 threads lastdb took 6 minutes on the
+            # fish genome and 34 on a mouse one (split into volumes), and training took 2 to 17 minutes.
+            db_letters = min(letters[event_a.iD], letters[event_b.iD])
             score_models[lastz_pair_key(event_a.iD, event_b.iD)] = root_job.addChildJobFn(
                 train_lastz_score_model, event_a.iD, genome_a, event_b.iD, genome_b, params,
                 cores=train_cores,
-                memory=cactus_clamp_memory(int(2e9 + 20 * db_size)),
-                disk=int(2e9 + 25 * db_size + 2 * max(genome_a.size, genome_b.size)),
-                walltime=cactus_walltime(1800 + 7200 * db_size / 1e9, io_bytes=genome_a.size + genome_b.size)).rv()
+                memory=cactus_clamp_memory(int(1e9 + 5.5 * db_letters)),
+                disk=int(2e9 + 5 * db_letters + genome_a.size + genome_b.size),
+                walltime=cactus_walltime(1800 + 1500 * db_letters / 1e9, io_bytes=genome_a.size + genome_b.size)).rv()
 
     if unmask_job is not None or score_models:
         new_root_job = Job(walltime=cactus_walltime())

From 64fdba8cbf9d4ac2b18769a5c4660cb83d28cf7e Mon Sep 17 00:00:00 2001
From: Glenn Hickey <glenn.hickey@gmail.com>
Date: Fri, 2 Oct 2026 19:29:59 -0400
Subject: [PATCH 4/6] blast: train lastz's scores with the inputs' own repeat
 masking

lastdb -R10 in place of LAST's default -R01: the inputs' soft-masking,
which cactus's Red preprocessing puts there and the reconstructed
ancestors keep, masks repeats from seeding, as it does for lastz, and
lastal follows lastdb's -R.  The default uppercased the inputs and
masked only simple repeats, with tantan, which was 29% of lastdb's time,
and left interspersed repeats seeding.

On whole genomes (4 threads) human against mouse took 991 CPU-s in all
rather than 4494, the training passes 135 rather than 1916, and the model
moved 5.6%, inside the 8.3% two equally good samples differ by; human
against a fish took 753 rather than 1144, moving 21.7% against 16.8%.
Two mouse assemblies trained to 100.00% identity rather than 99.35%:
the default had been fitting alignments between repeat copies.  On
evolver mammals the alignment was as accurate as before (pooled
leaf-pair precision/recall 0.830/0.744 both ways).

LAST 1651, as cactus installs it, gives the same models as the patched
last-train with --tolerance, in a few more passes that now cost little.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/cactus/paf/local_alignment.py | 9 +++++++--
 1 file changed, 7 insertions(+), 2 deletions(-)

diff --git a/src/cactus/paf/local_alignment.py b/src/cactus/paf/local_alignment.py
index 2f7c4e869..b8dca6f60 100755
--- a/src/cactus/paf/local_alignment.py
+++ b/src/cactus/paf/local_alignment.py
@@ -124,14 +124,19 @@ def train_lastz_score_model(job, event_a, genome_a, event_b, genome_b, params):
     trained_gaps = getOptionalAttrib(params.find("blast"), 'lastzTrainGaps', typeFn=bool, default=True)
 
     # lastdb -c and -C2 are LAST's recipe for training on genomes: -c keeps masked letters out of the seeds, and
-    # -C2 drops gapless alignments nested in two others
+    # -C2 drops gapless alignments nested in two others.  -R10 makes the masked letters the inputs' soft-masking
+    # (cactus's Red), which masks repeats from seeding as it does for lastz, in place of LAST's default of masking
+    # only simple repeats with tantan (lastal follows lastdb's -R).  That makes lastdb ~3x and training ~10x cheaper
+    # on a mammal genome (human against mouse: 991 rather than 4494 CPU-s in all), and keeps interspersed repeats
+    # out of the model: two mouse assemblies trained to 99.35% identity without, from alignments between repeat
+    # copies, and 100.00% with.
     train_cmd = ['last-train', '--revsym', '--matsym', '--gapsym', '-C2', '-P', str(job.cores)]
     if last_train_has_tolerance():
         train_cmd.append('--tolerance=1')
     train_cmd += ['db', local_fastas[1]]
     train_file = os.path.join(work_dir, 'pair.train')
     try:
-        cactus_call(parameters=['lastdb', '-c', '-P', str(job.cores), 'db', local_fastas[0]], work_dir=work_dir)
+        cactus_call(parameters=['lastdb', '-c', '-R10', '-P', str(job.cores), 'db', local_fastas[0]], work_dir=work_dir)
         with open(os.path.join(work_dir, 'db.prj')) as prj:
             db_letters = int(next(line for line in prj if line.startswith('numofletters=')).split('=')[1])
         if db_letters < 500000:

From a7a797b76718b65e2138e6e854f1aa69ed3f47b9 Mon Sep 17 00:00:00 2001
From: Glenn Hickey <glenn.hickey@gmail.com>
Date: Mon, 5 Oct 2026 16:42:32 -0400
Subject: [PATCH 5/6] blast: train lastz's scores by default

<blast lastzTrain> is now on.  It passes over the mappers it does not
apply to rather than stopping them: minimap2 and FastGA quietly, KegAlign
with a warning (it used to be an error).  A pair it cannot train, as when
LAST is missing, falls back to the lastzScoreModel table or HOXD70 as
before.  On evolver mammals it raised pooled leaf-pair precision from
0.816 to 0.830 and recall from 0.740 to 0.744; on the 577-way it would
cost about 0.3% of the lastz CPU.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/cactus/cactus_progressive_config.xml | 8 +++++---
 src/cactus/paf/local_alignment.py        | 8 ++++----
 2 files changed, 9 insertions(+), 7 deletions(-)

diff --git a/src/cactus/cactus_progressive_config.xml b/src/cactus/cactus_progressive_config.xml
index 7f0992654..8f4a4389b 100644
--- a/src/cactus/cactus_progressive_config.xml
+++ b/src/cactus/cactus_progressive_config.xml
@@ -91,8 +91,10 @@
 	<!-- leftAlign: move each indel in the lastz alignments along its repeat without changing the alignment's score (paffy left_align). "none", or "canonical" (default): at whichever end makes the repeat read smaller than its reverse complement, which gives every pairwise alignment through a repeat the same placement whatever the pair or strand, so they agree when pinched together. "left" (left on the target) only makes alignments agree when no sequence is aligned on its reverse strand, and makes things worse otherwise; it is here for comparison. lastz only -->
 	<!-- lastzTrain Train lastz's scores on each pair of genomes it aligns with last-train (needs LAST), on the whole
 	     genomes, and align the pair with them, rescaled to HOXD70's units (see lastzScoreModel below).  A pair that
-	     cannot be trained (the smaller genome under 500kb, or last-train failing) falls back to the lastzScoreModel
-	     table, else HOXD70.  Not supported with KegAlign -->
+	     cannot be trained (the smaller genome under 500kb, or last-train failing, as without LAST) falls back to the
+	     lastzScoreModel table, else HOXD70.  It does not apply to KegAlign, minimap2 or FastGA.  On evolver mammals
+	     it raised pooled leaf-pair precision from 0.816 to 0.830 and recall from 0.740 to 0.744, for about 0.3% of
+	     the 577-way's lastz CPU; set it to 0 for lastz's HOXD70 everywhere -->
 	<!-- lastzTrainGaps With lastzTrain, use the trained gap penalties; 0 keeps HOXD70's (400 and 30) with the trained
 	     substitution scores -->
 	<!-- lastzTrainCores Cores for each pair's training -->
@@ -131,7 +133,7 @@
 		   pickIngroupPrimaryAlignmentsSeparatelyToOutgroups="1"
 		   includeRootAsOutgroup="1"
 		   bigGenomeChunkScale="2"
-		   lastzTrain="0"
+		   lastzTrain="1"
 		   lastzTrainGaps="1"
 		   lastzTrainCores="4"
 		   >
diff --git a/src/cactus/paf/local_alignment.py b/src/cactus/paf/local_alignment.py
index b8dca6f60..5934f813f 100755
--- a/src/cactus/paf/local_alignment.py
+++ b/src/cactus/paf/local_alignment.py
@@ -68,16 +68,16 @@ def write_lastz_scores(model, path):
 
 
 def lastz_train_enabled(params):
-    """ Whether to train lastz's scores on each pair of genomes it aligns: <blast lastzTrain> """
+    """ Whether to train lastz's scores on each pair of genomes it aligns: <blast lastzTrain>, on by default, so
+    it passes over the mappers it does not apply to, minimap2, FastGA and KegAlign, rather than stopping them """
     lastz_params_node = params.find("blast")
     if not getOptionalAttrib(lastz_params_node, 'lastzTrain', typeFn=bool, default=False):
         return False
     if lastz_params_node.attrib["mapper"] != "lastz":
-        logger.warning("<blast lastzTrain> is ignored: it is for lastz, and the mapper is {}".format(
-            lastz_params_node.attrib["mapper"]))
         return False
     if getOptionalAttrib(lastz_params_node, 'gpu', typeFn=int, default=0):
-        raise RuntimeError("<blast lastzTrain> is not supported with KegAlign")
+        logger.warning("<blast lastzTrain> does not apply to KegAlign, whose alignments keep its default scores")
+        return False
     return True
 
 

From 2988906c5a5858e3e952d2c39dd1e9083dc19d91 Mon Sep 17 00:00:00 2001
From: Glenn Hickey <glenn.hickey@gmail.com>
Date: Tue, 6 Oct 2026 07:37:11 -0400
Subject: [PATCH 6/6] blast: ask for no more cores for training than Toil
 allows

Toil refuses a job that asks for more cores than --maxCores, or on a
single machine than it has, so lastz training's 4 (lastzTrainCores)
failed every cactus-blast run with fewer -- as cactus-prepare-toil's,
whose blast jobs run cactus-blast with --maxCores 1.  cactus now keeps
that limit in CACTUS_MAX_CORES, as it does --maxMemory's, and
cactus_clamp_cores holds a job to it.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/cactus/paf/local_alignment.py |  4 ++--
 src/cactus/shared/common.py       | 12 ++++++++++++
 2 files changed, 14 insertions(+), 2 deletions(-)

diff --git a/src/cactus/paf/local_alignment.py b/src/cactus/paf/local_alignment.py
index 5934f813f..8eaabd844 100755
--- a/src/cactus/paf/local_alignment.py
+++ b/src/cactus/paf/local_alignment.py
@@ -25,7 +25,7 @@
 from cactus.preprocessor.checkUniqueHeaders import sanitize_fasta_headers
 from cactus.preprocessor.unmasking import unmask_contigs_all
 from cactus.preprocessor.cactus_preprocessor import clean_if_different
-from cactus.shared.common import cactus_clamp_memory
+from cactus.shared.common import cactus_clamp_memory, cactus_clamp_cores
 
 def get_divergence_class(distance, params):
     """ The <constants><divergences> bucket a pair at the given distance falls in ("one" through
@@ -1143,7 +1143,7 @@ def make_paf_alignments(job, event_tree_string, event_names_to_sequences, ancest
             db_letters = min(letters[event_a.iD], letters[event_b.iD])
             score_models[lastz_pair_key(event_a.iD, event_b.iD)] = root_job.addChildJobFn(
                 train_lastz_score_model, event_a.iD, genome_a, event_b.iD, genome_b, params,
-                cores=train_cores,
+                cores=cactus_clamp_cores(train_cores),
                 memory=cactus_clamp_memory(int(1e9 + 5.5 * db_letters)),
                 disk=int(2e9 + 5 * db_letters + genome_a.size + genome_b.size),
                 walltime=cactus_walltime(1800 + 1500 * db_letters / 1e9, io_bytes=genome_a.size + genome_b.size)).rv()
diff --git a/src/cactus/shared/common.py b/src/cactus/shared/common.py
index 05244c851..bdcaacf0c 100644
--- a/src/cactus/shared/common.py
+++ b/src/cactus/shared/common.py
@@ -229,6 +229,13 @@ def cactus_override_toil_options(options):
     os.environ['CACTUS_MAX_MEMORY'] = str(max_mem)
     os.environ['CACTUS_DEFAULT_MEMORY'] = str(human2bytes(str(options.defaultMemory)) if options.defaultMemory else 2**31)
 
+    # and the most cores a job can ask for: Toil refuses a job that asks for more than --maxCores, or on a single
+    # machine than the machine has (see cactus_clamp_cores)
+    max_cores = sys.maxsize if options.maxCores in (None, float('inf')) else int(options.maxCores)
+    if options.batchSystem.lower() in ['single_machine', 'singleMachine']:
+        max_cores = min(max_cores, cactus_cpu_count())
+    os.environ['CACTUS_MAX_CORES'] = str(max(1, max_cores))
+
     # store the walltime knobs here so cactus_walltime() can get at them without carrying
     # options around, and so they propagate to workers just like the memory limits above
     os.environ['CACTUS_WALLTIME_FACTOR'] = str(getattr(options, 'walltimeFactor', WALLTIME_FACTOR))
@@ -261,6 +268,11 @@ def cactus_override_toil_options(options):
     except:
         pass
 
+def cactus_clamp_cores(cores):
+    """ use the environment variable from --maxCores (and a single machine's cores) to clamp a job's cores, as
+    cactus_clamp_memory does its memory: Toil refuses a job that asks for more """
+    return max(1, min(int(os.environ.get('CACTUS_MAX_CORES', cores)), int(cores)))
+
 def cactus_clamp_memory(memory_bytes):
     """ use the environment variables from --maxMemory and --defaultMemory to clamp a given memory value """
     return max(min(int(os.environ['CACTUS_MAX_MEMORY']), int(memory_bytes)), int(os.environ['CACTUS_DEFAULT_MEMORY']))