Fix abPOA's cg_backtrack crash with large trained scores; divide trained models by their common factor

Makefile

@@ -93,6 +93,7 @@ all_libs.blastLib: all_libs.api
 testModules = \
     progressive/outgroupTest.py \
     hal/cactus_validateTest.py \
+    paf/pafTest.py \
     preprocessor/cactus_preprocessorTest.py \
     preprocessor/checkPreprocessedSequenceTest.py \
     preprocessor/lastzRepeatMasking/cactus_lastzRepeatMaskTest.py \

ReleaseNotes.md

@@ -1,5 +1,6 @@
 # Unreleased
 
+- Fix abPOA exiting with `Error in cg_backtrack` when its scores are large, which crashed a CHM13 HPRC build on chr2, whose trained scoring model learned a gap extension of 2. abPOA now picks 16- or 32-bit scores from bounds on the scores, checks the alignments it cannot prove correct, and redoes those at 32 bits. Trained scoring models are also divided by the common factor of their scores, which gives the same alignments with smaller numbers and keeps abPOA in range, and an unusable trained model is rejected and falls back as if training had failed.
 - The GAF overlap filter now resolves overlapping minigraph mappings exactly, per query segment, instead of deleting whole mappings (`gaffilter -x` from cactus-gfa-tools, `GAFOverlapFilterExact="1"` in `<graphmap>`, on by default): each segment goes to the mapping that beats every other one claiming it, newly anchored sequence must be 98% identical and may not place reference sequence its haplotype already covers, and new junctions that cannot be bounded rearrangements are broken by a 31 kb query gap. On a whole-genome CHM13 30-way pangenome it had 83k (3.5%) fewer dipcall errors (without normalization) over 14 samples, 13 of them better (the other is worse on chrY's ampliconic repeats), left HG002 against GIAB v5.0q unchanged, and represented more real inversions as two-sided snarls; on chr16 alone it anchored 15.3 Mb more (15 Mb less clipping). A review log of junctions between pieces placed apart on the reference, with how many haplotypes share each, is written to `<paf>.junctions.tsv`. `GAFOverlapFilterExact="0"` (or a config without it) gives the previous whole-mapping deletion. `GAFOverlapFilterTrim` (`gaffilter -t`) has been removed, and a config that sets it to 1 is an error. See "Overlapping Mappings" in [doc/pangenome.md](doc/pangenome.md).
 - Fix minigraph mappings with wrong base-level alignments: minigraph (since r508, so every release since 2.1.0) printed the `cg` CIGAR of every GAF record that followed a query's first reverse-strand stable-coordinate record with its operations in reverse order, while the strand, path, coordinates and counts stayed correct. `gaf2paf` turned these into PAF anchors that pinch query bases to the wrong graph bases. Most such records are short pericentromeric alignments that the 250 kb `minGAFBlockLength` drops (about 0.01% of what `cactus-align` received in a CHM13 30-way test), but the GAFs themselves were affected much more (27% of records there). Graph construction was not affected. minigraph is now pinned at [c59d03d](https://github.com/glennhickey/minigraph/commit/c59d03d), which resets the flag for each record. GAFs from earlier runs that are reused with `--inGAF` still carry the bug; map those genomes again.
 - `rgfa-zip` (new in cactus-gfa-tools, and installed by `downloadPangenomeTools`) now rewrites each minigraph graph after construction, on by default (`zipAlleles="1"` in `<graphmap>`; `zipAlleles="0"` turns it off). Alleles minigraph stored as novel sequence are aligned, on both strands, to the reference window they bypass, and confident copies (5 kb at 95% gap-compressed identity, uniquely placed, by default) are replaced by the sequence they duplicate, inversions included. Graphmap, the split, `cactus-align` and the join all see the zipped graph, and `<x>.sv.unzipped.gfa.gz` and its report `<x>.sv.zip.tsv` are kept beside `<x>.sv.gfa.gz`. `zipOptions` is passed through, and defaults to `rgfa-zip`'s own thresholds plus `--no-alt`, which zips onto the reference only: without `--no-alt` a second pass zips alleles onto parallel alleles with the same bounding nodes, which on a CHM13 30-way pangenome was a small but consistent loss. `zipWalks="gaf"` first maps every genome to the unzipped graph as graphmap does, so that the zip works from observed haplotype walks, kept as `<x>.sv.unzipped.gaf.gz`. Under `--mgSplit` it runs per chromosome, on each chromosome's whole-reference graph before the prune (see below). It always aligns with the `minimap2` cactus ships, and cannot be combined with `--inGAF` (without `--remap`) when adding genomes to a graph, so that now needs `zipAlleles="0"`; resuming from a graph with nothing to add is unaffected. A config without `zipAlleles`, such as one made from an older default, leaves the zip off.

bar/tests/poaBarTest.c

@@ -482,12 +482,71 @@ void test_make_oriented_partial_order_alignment(CuTest *testCase) {
     }
 }
 
+/*
+ * abPOA with the scores of the HPRC oct9 CHM13 chr2 model (a learned gap extension of 2, scaled 30x): the abPOA
+ * cactus pinned before could not hold them in its 16-bit DP, whose "minus infinity" floor came out above real
+ * scores, and aborted with "Error in cg_backtrack" or misaligned silently.  The same model with every score
+ * divided by div is the same model, and must give the same alignment.
+ */
+static PoaParameters *chr2_model_params(int div) {
+    static const int chr2_mat[25] = { 150, -960, -780, -1020, -780,  -960, 210, -900, -780, -780,
+                                     -780, -900, 210, -960, -780,  -1020, -780, -960, 150, -780,
+                                     -780, -780, -780, -780, 150 };
+    PoaParameters *poaParams = st_calloc(1, sizeof(PoaParameters));
+    poaParams->engine = BASE_ALIGNER_ABPOA;
+    abpoa_para_t *abpt = abpoa_init_para();
+    abpt->wb = 1000; abpt->wf = 0.1; // the shipped band, under which every short window is a full-band DP
+    abpt->gap_open1 = 1290 / div; abpt->gap_ext1 = 60 / div; abpt->gap_open2 = 3870 / div; abpt->gap_ext2 = 2 / div;
+    abpt->progressive_poa = 1;
+    abpoa_post_set_para(abpt);
+    // as abpoaParamaters_constructFromCactusParams overrides the matrix
+    abpt->use_score_matrix = 1; abpt->max_mat = abpt->min_mis = 0;
+    for (int i = 0; i < 25; i++) {
+        abpt->mat[i] = chr2_mat[i] / div;
+        if (abpt->mat[i] > abpt->max_mat) abpt->max_mat = abpt->mat[i];
+        if (-abpt->mat[i] > abpt->min_mis) abpt->min_mis = -abpt->mat[i];
+    }
+    poaParams->abpt = abpt;
+    return poaParams;
+}
+
+void test_trained_scores_beyond_16_bits(CuTest *testCase) {
+    PoaParameters *full = chr2_model_params(1), *half = chr2_model_params(2);
+    for (int64_t test = 0; test < 200; test++) {
+        char *parent_string = getRandomACGTSequence(st_randomInt(1, 150));
+        int64_t seq_no = st_randomInt(2, 20);
+        char **seqs = st_malloc(sizeof(char *) * seq_no), **seqs_half = st_malloc(sizeof(char *) * seq_no);
+        int *seq_lens = st_malloc(sizeof(int) * seq_no), *seq_lens_half = st_malloc(sizeof(int) * seq_no);
+        for (int64_t i = 0; i < seq_no; i++) {
+            seqs[i] = evolveSequence(parent_string);
+            seqs_half[i] = stString_copy(seqs[i]);
+            seq_lens[i] = seq_lens_half[i] = strlen(seqs[i]);
+        }
+        // the msas take ownership of their sequences
+        Msa *msa = msa_make_partial_order_alignment(seqs, seq_lens, seq_no, 10000, 5000, 1.0, full);
+        Msa *msa_half = msa_make_partial_order_alignment(seqs_half, seq_lens_half, seq_no, 10000, 5000, 1.0, half);
+        int64_t lengths[seq_no];
+        validate_msa(testCase, msa, lengths);
+        CuAssertIntEquals(testCase, msa_half->column_no, msa->column_no);
+        for (int64_t i = 0; i < seq_no; i++) {
+            CuAssertTrue(testCase, lengths[i] == seq_lens[i]);
+            char *r = msa_row_string(msa, i), *r_half = msa_row_string(msa_half, i);
+            CuAssertStrEquals(testCase, r_half, r);
+            free(r); free(r_half);
+        }
+        msa_destruct(msa); msa_destruct(msa_half);
+        free(parent_string);
+    }
+    poaParameters_destruct(full); poaParameters_destruct(half);
+}
+
 CuSuite* poaBarAlignerTestSuite(void) {
     CuSuite* suite = CuSuiteNew();
     SUITE_ADD_TEST(suite, test_shipped_config_is_loadable);
     SUITE_ADD_TEST(suite, test_baseAligner_selection);
     SUITE_ADD_TEST(suite, test_make_partial_order_alignment);
     SUITE_ADD_TEST(suite, test_make_oriented_partial_order_alignment);
+    SUITE_ADD_TEST(suite, test_trained_scores_beyond_16_bits);
     SUITE_ADD_TEST(suite, test_make_consistent_partial_order_alignments_two_ends);
     SUITE_ADD_TEST(suite, test_make_flower_alignment_poa);
     SUITE_ADD_TEST(suite, test_get_adjacency_string_and_overlap_bounded);

src/cactus/paf/last_scoring.py

@@ -58,6 +58,12 @@ def parse_train_file(train_file_path):
         for key in ['GAP-OPEN', 'GAP-EXTEND', 'A', 'C', 'G', 'T']:
             if key not in score_dict:
                 raise RuntimeError('Information for {} not parsed from {}'.format(key, train_file_path))
+        # apply_long_gap() needs an extension of at least 1 to scale by, and a model with a negative gap
+        # open or a match score of 0 or less is not one the aligners can use.  Raising here makes
+        # last_train() fall back to the defaults (or a borrowed model) rather than hand such a model on.
+        if score_dict['GAP-EXTEND'] < 1 or score_dict['GAP-OPEN'] < 0 or min(score_dict[b][b] for b in 'ACGT') <= 0:
+            raise RuntimeError('Unusable scoring model in {}: gap open {}, gap extend {}, matches {}'.format(
+                train_file_path, score_dict['GAP-OPEN'], score_dict['GAP-EXTEND'], [score_dict[b][b] for b in 'ACGT']))
         for i in ['A', 'C', 'G', 'T']:
             rci = { 'A' : 'T', 'C' : 'G', 'G' : 'C', 'T' : 'A' }[i]
             for j in ['A', 'C', 'G', 'T']:
@@ -90,10 +96,52 @@ def apply_long_gap(score_dict, open_factor, extend_factor):
 
     score_dict['GAP-OPEN-2'] = score_dict['GAP-OPEN'] * open_factor
     score_dict['GAP-EXTEND-2'] = max(1, score_dict['GAP-EXTEND'] // extend_factor)
-
-def apply_scores_to_config(score_dict, config_xml):
+    return reduce_scores(score_dict)
+
+def reduce_scores(score_dict):
+    """ divide every score of a model (substitutions and both gap pieces) by their greatest common
+    divisor, and return it.  The aligners' choices, ties included, are unchanged by a common positive
+    factor, so the alignments are identical with smaller numbers.  apply_long_gap() scales by
+    extend_factor (30) to make room for a long-gap extension 30 times cheaper, and a learned extension
+    that shares a factor with it leaves every score that many times larger than needed: a learned
+    extension of 2 reached abPOA as 60 (GapOpen 1290; GapExtend 60; GapOpen2 3870; GapExtend2 2), more
+    than the abPOA cactus pinned before could hold in its 16-bit DP (the HPRC oct9 CHM13 chr2 model
+    aborted with "Error in cg_backtrack", and misaligned windows that did not abort).  Divided by 2 it
+    is 645, 30, 1935, 1.  Extensions of 1 (all but one model so far) are left alone. """
+    bases = ['A', 'C', 'G', 'T']
+    gaps = ['GAP-OPEN', 'GAP-EXTEND', 'GAP-OPEN-2', 'GAP-EXTEND-2']
+    d = 0
+    for v in [score_dict[i][j] for i in bases for j in bases] + [score_dict[k] for k in gaps]:
+        d = math.gcd(d, abs(v))
+    if d > 1:
+        for i in bases:
+            for j in bases:
+                score_dict[i][j] //= d
+        for k in gaps:
+            score_dict[k] //= d
+    return max(d, 1)
+
+def abpoa_exact_depth(matrix, gap_open1, gap_ext1, gap_open2, gap_ext2, window):
+    """ the deepest graph (rows on a path) for which abPOA can prove, up front, that a window-long
+    sequence aligns exactly in 32 bits, in its worst case (AVX-512 lanes, a band narrower than the
+    window).  This mirrors simd_abpoa_score_range() in cactus's abPOA: deeper graphs still align
+    exactly, but abPOA has to check every score as it goes, and aborts with a message rather than
+    give a wrong alignment if the scores ever outgrow 32 bits.  Trained models allow over 100
+    window-lengths of depth; much less means the scores are far larger than they need to be.
+    Path scores are left out, as cactus does not use them. """
+    M, X = max(0, max(matrix)), max(0, -min(matrix))
+    convex = gap_open1 > 0 and gap_open2 > 0  # abPOA's affine and linear modes use one piece
+    oe = [gap_open1 + gap_ext1] + ([gap_open2 + gap_ext2] if convex else [])
+    OE, EM, lanes = max(oe), max(gap_ext1, gap_ext2) if convex else gap_ext1, 16
+    floor = -2**31 + (window + 2) * X + OE + lanes // 2 * EM
+    ztop = floor + (window + 1) * M + M + 2 * X + 2 * OE
+    # every real score is at least -max(X, min(oe)) * (window + depth)
+    return max(0, (-ztop - 1) // max(1, X, min(oe)) - window)
+
+def apply_scores_to_config(score_dict, config_xml, label=None):
     """ load the score dict into the config.  since last won't train long gaps,
-    should be specified by the parameters.  note that they are disabled with 0, untouched with None"""
+    should be specified by the parameters.  note that they are disabled with 0, untouched with None.
+    label (the chromosome, say) only goes into the log, so a model can be traced back to its job """
 
     bar_node = findRequiredNode(config_xml, "bar")
     poa_node = findRequiredNode(bar_node, "poa")
@@ -103,10 +151,16 @@ def apply_scores_to_config(score_dict, config_xml):
     # that's applied to the trained gap parameters
     long_gap_open_factor = int(poa_node.attrib['partialOrderAlignmentTrainedGapOpen2Factor'])
     long_gap_extend_factor = int(poa_node.attrib['partialOrderAlignmentTrainedGapExtension2Factor'])
-    # the open factor was being passed the extend factor, so TrainedGapOpen2Factor was read and
-    # then discarded.  Both default to 3, so the shipped config is unaffected; setting them to
-    # different values silently did the wrong thing.
-    apply_long_gap(score_dict, long_gap_open_factor, long_gap_extend_factor)
+    # the open factor used to be passed the extend factor, so TrainedGapOpen2Factor was read and
+    # then discarded.  That went unnoticed while both defaulted to 3; the defaults are now 3 for the
+    # open and 30 for the extension.
+    trained_ext = score_dict['GAP-EXTEND']
+    divisor = apply_long_gap(score_dict, long_gap_open_factor, long_gap_extend_factor)
+    if divisor > 1:
+        RealtimeLogger.info('Trained scoring model{} learned a gap extension of {}: its scores, multiplied by {} for the '
+                            'long-gap piece, are all divided by their common factor {} (the alignments are the same, '
+                            'the numbers smaller)'.format(' for ' + label if label else '', trained_ext,
+                                                          long_gap_extend_factor, divisor))
 
     poa_node.attrib['partialOrderAlignmentGapOpenPenalty1'] = str(score_dict['GAP-OPEN'])
     poa_node.attrib['partialOrderAlignmentGapExtensionPenalty1'] = str(score_dict['GAP-EXTEND'])
@@ -139,6 +193,16 @@ def apply_scores_to_config(score_dict, config_xml):
 
     poa_node.attrib['partialOrderAlignmentSubMatrix'] = score_string
 
+    window = int(poa_node.attrib.get('partialOrderAlignmentWindow', 10000))
+    depth = abpoa_exact_depth([int(v) for v in score_string.split()], score_dict['GAP-OPEN'], score_dict['GAP-EXTEND'],
+                              score_dict['GAP-OPEN-2'], score_dict['GAP-EXTEND-2'], window)
+    if depth < 50 * window:
+        RealtimeLogger.warning('Trained scoring model{} is very large (gap open {}, gap extend {}, matrix {}): abPOA can only '
+                               'guarantee exact {} bp alignments to graphs up to {} deep, and will check deeper ones as it goes, '
+                               'aborting if they overflow.  Consider --scoresFile or lastTrain="0"'.format(
+                                   ' for ' + label if label else '', score_dict['GAP-OPEN'], score_dict['GAP-EXTEND'],
+                                   score_string, window, depth))
+
     # minipoa gets the learned gaps verbatim, as a single affine piece: last-train fits a single
     # affine gap model, and the GAP-OPEN-2/GAP-EXTEND-2 pair above is a synthesised second piece
     # that exists only because abPOA is unstable without one, so it would be wrong to hand it on.
@@ -156,12 +220,14 @@ def apply_scores_to_config(score_dict, config_xml):
             minipoa_node.attrib['minipoaGapOpenPenalty'],
             minipoa_node.attrib['minipoaGapExtensionPenalty']))
 
-    RealtimeLogger.info("Overriding abPOA scores with trained values: GapOpen {}; GapExtend {}; GapOpen2 {}; GapExtend2 {}; SubMatrix {}".format(
+    # the line keeps its old prefix, which logs are grepped for; the label and scaling go on the end
+    RealtimeLogger.info("Overriding abPOA scores with trained values: GapOpen {}; GapExtend {}; GapOpen2 {}; GapExtend2 {}; SubMatrix {} ({}learned gap extension {}; all scores multiplied by {}{})".format(
         poa_node.attrib['partialOrderAlignmentGapOpenPenalty1'],
         poa_node.attrib['partialOrderAlignmentGapExtensionPenalty1'],
         poa_node.attrib['partialOrderAlignmentGapOpenPenalty2'],
         poa_node.attrib['partialOrderAlignmentGapExtensionPenalty2'],
-        poa_node.attrib['partialOrderAlignmentSubMatrix']))
+        poa_node.attrib['partialOrderAlignmentSubMatrix'], label + ': ' if label else '', trained_ext, long_gap_extend_factor,
+        ', then divided by their common factor {}'.format(divisor) if divisor > 1 else ''))
 
 def minigraph_wfa_penalties(score_dict):
     """ convert a last-train model into minigraph's base-alignment penalties, as (x, o, e) for