Optional lastz scoring by divergence: a table, or per-pair last-train models (prototype)

src/cactus/cactus_progressive_config.xml

@@ -89,6 +89,15 @@
 	     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, 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 -->
 	<blast chunkSize="90000000"
 	       bigChunkSize="6000000000"
 		   overlapSize="10000"
@@ -124,6 +133,9 @@
 		   pickIngroupPrimaryAlignmentsSeparatelyToOutgroups="1"
 		   includeRootAsOutgroup="1"
 		   bigGenomeChunkScale="2"
+		   lastzTrain="1"
+		   lastzTrainGaps="1"
+		   lastzTrainCores="4"
 		   >
 
 		<!-- The following are parametrised to produce the same results as the default settings,
@@ -145,6 +157,13 @@
 				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; 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

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,69 @@ 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 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):'
+    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 values(key):
+        found = re.findall(r'# {}: ([-\d.eE+]+)'.format(key), txt)
+        if not found:
+            raise RuntimeError('no {} in {}'.format(key, train_file_path))
+        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(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(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
+    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.