diff --git a/Makefile b/Makefile index 3683fa7fc..e026a3300 100644 --- a/Makefile +++ b/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 \ diff --git a/ReleaseNotes.md b/ReleaseNotes.md index 5ff4dda60..32f851c74 100644 --- a/ReleaseNotes.md +++ b/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 ``, 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 `.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 ``; `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 `.sv.unzipped.gfa.gz` and its report `.sv.zip.tsv` are kept beside `.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 `.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. diff --git a/bar/tests/poaBarTest.c b/bar/tests/poaBarTest.c index 163a9f242..0737a7797 100644 --- a/bar/tests/poaBarTest.c +++ b/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); diff --git a/src/cactus/paf/last_scoring.py b/src/cactus/paf/last_scoring.py index d006cf190..96d6ac8f3 100644 --- a/src/cactus/paf/last_scoring.py +++ b/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 diff --git a/src/cactus/paf/pafTest.py b/src/cactus/paf/pafTest.py index a1f3a4f29..ac11ad598 100644 --- a/src/cactus/paf/pafTest.py +++ b/src/cactus/paf/pafTest.py @@ -8,18 +8,24 @@ import unittest from sonLib.bioio import newickTreeParser -from paf import * +from cactus.paf.paf import * +from cactus.paf.last_scoring import apply_long_gap, apply_scores_to_config, abpoa_exact_depth, parse_train_file +from toil.realtimeLogger import RealtimeLogger +import logging +import os +import tempfile +import xml.etree.ElementTree as ET class TestCase(unittest.TestCase): def setUp(self): self.tree_string = "((simCow:0.1,simDog:0.2)anc1:0.3,(simChimp:0.4,simHuman:0.5)anc2:0.2)anc0:0.1;" self.tree = newickTreeParser(self.tree_string) - def test_get_leaf_event_pairs(self): - self.assertEqual(6, len(list(get_leaf_event_pairs(self.tree)))) - distances = get_distances(self.tree) + def test_get_event_pairs(self): leaves = get_leaves(self.tree) - for species_a, species_b, distance in get_leaf_event_pairs(self.tree): + self.assertEqual(6, len(list(get_event_pairs(self.tree, leaves)))) + distances = get_distances(self.tree) + for species_a, species_b, distance in get_event_pairs(self.tree, leaves): self.assertTrue(species_a in leaves) self.assertTrue(species_b in leaves) self.assertEqual(distances[(species_a, species_b)], distance) @@ -46,5 +52,102 @@ def test_get_subtree_nodes(self): def test_get_leaves(self): self.assertEqual({ "simCow", "simDog", "simChimp", "simHuman" }, { i.iD for i in get_leaves(self.tree) }) +class LastScoringTest(unittest.TestCase): + """ how a last-train model is scaled for abPOA (cactus.paf.last_scoring) """ + ROWS = {'A': [5, -32, -26, -34], 'C': [-32, 7, -30, -26], 'G': [-26, -30, 7, -32], 'T': [-34, -26, -32, 5]} + + def model(self, gap_open, gap_extend): + # the HPRC oct9 CHM13 chr2 model's matrix (that chromosome learned an extension of 2, the rest 1) + d = {i: dict(zip('ACGT', self.ROWS[i])) for i in 'ACGT'} + d['GAP-OPEN'], d['GAP-EXTEND'] = gap_open, gap_extend + return d + + @staticmethod + def gaps(d): + return d['GAP-OPEN'], d['GAP-EXTEND'], d['GAP-OPEN-2'], d['GAP-EXTEND-2'] + + @staticmethod + def config(): + return ET.fromstring('') + + @staticmethod + def write_train(path, gap_open, gap_extend, rows): + with open(path, 'w') as f: + f.write('#last -a {0}\n#last -A {0}\n#last -b {1}\n#last -B {1}\n'.format(gap_open, gap_extend)) + f.write(' A C G T\n') + for i in 'ACGT': + f.write('{} {}\n'.format(i, ' '.join(str(v) for v in rows[i]))) + + def test_extension_one_unchanged(self): + # the usual model: scaled 30x and nothing else, as before + d = self.model(43, 1) + self.assertEqual(apply_long_gap(d, 3, 30), 1) + self.assertEqual(self.gaps(d), (1290, 30, 3870, 1)) + self.assertEqual(d['A']['C'], -960) + + def test_extension_two_reduced(self): + # scaled 30x, every score is even: halved, the same model with abPOA's usual extension of 30 + d = self.model(43, 2) + self.assertEqual(apply_long_gap(d, 3, 30), 2) + self.assertEqual(self.gaps(d), (645, 30, 1935, 1)) + self.assertEqual([d['A'][j] for j in 'ACGT'], [75, -480, -390, -510]) + + def test_reduction_is_exact(self): + # whatever the extension, the result is the 30x model divided by a common factor + for ext in range(1, 9): + d, ref = self.model(43, ext), self.model(43, ext) + div = apply_long_gap(d, 3, 30) + self.assertEqual(div, {1: 1, 2: 2, 3: 3, 4: 2, 5: 5, 6: 6, 7: 1, 8: 2}[ext]) + self.assertEqual(self.gaps(d), (43 * 30 // div, ext * 30 // div, 43 * 90 // div, ext // div)) + for i in 'ACGT': + for j in 'ACGT': + self.assertEqual(d[i][j] * div, ref[i][j] * 30) + + def test_config_gets_reduced_scores(self): + root = self.config() + apply_scores_to_config(self.model(43, 2), root, label='chr2') + poa, minipoa = root.find('bar/poa'), root.find('bar/minipoa') + self.assertEqual([poa.attrib['partialOrderAlignmentGapOpenPenalty1'], poa.attrib['partialOrderAlignmentGapExtensionPenalty1'], + poa.attrib['partialOrderAlignmentGapOpenPenalty2'], poa.attrib['partialOrderAlignmentGapExtensionPenalty2']], + ['645', '30', '1935', '1']) + self.assertEqual(poa.attrib['partialOrderAlignmentSubMatrix'].split()[:5], ['75', '-480', '-390', '-510', '-390']) + self.assertEqual((minipoa.attrib['minipoaGapOpenPenalty'], minipoa.attrib['minipoaGapExtensionPenalty']), ('645', '30')) + + def test_config_log(self): + # dividing a model is exact, so nothing is warned about, and the override line keeps the prefix logs are grepped for + RealtimeLogger.getLogger() # outside a workflow its first use sets it to CRITICAL, which would undo assertLogs + with self.assertLogs('toil-rt', level='INFO') as logs: + apply_scores_to_config(self.model(43, 2), self.config(), label='chr2') + self.assertEqual([r.getMessage() for r in logs.records if r.levelno >= logging.WARNING], []) + override = [r.getMessage() for r in logs.records if r.getMessage().startswith( + 'Overriding abPOA scores with trained values: GapOpen 645; GapExtend 30; GapOpen2 1935; GapExtend2 1; SubMatrix 75 ')] + self.assertEqual(len(override), 1) + self.assertIn('chr2: learned gap extension 2; all scores multiplied by 30, then divided by their common factor 2', override[0]) + + def test_exact_depth(self): + # abPOA's own numbers for a trained model at the 10 kb window: far beyond any window's graph + d = self.model(47, 1) + apply_long_gap(d, 3, 30) + mat = [d[i][j] for i in 'ACGT' for j in 'ACGT'] + [-990, 210, -990, 210, 150] + self.assertGreater(abpoa_exact_depth(mat, *self.gaps(d), 10000), 100 * 10000) + # four times the scores leaves less than 50 windows of depth, which is warned about + self.assertLess(abpoa_exact_depth([4 * v for v in mat], *[4 * g for g in self.gaps(d)], 10000), 50 * 10000) + # convex gaps whose two pieces have the same open + extension (110): the larger extension (60) still counts, + # as in abPOA's simd_abpoa_score_range(), which gives this depth + self.assertEqual(abpoa_exact_depth([10 if i % 6 == 0 else -20 for i in range(25)], 100, 10, 50, 60, 10000), 19509843) + + def test_unusable_model_rejected(self): + # an extension below 1, a negative open, or any match score of 0 or less (here A and T, the others positive) + zero_at = {i: [0 if i == j and i in 'AT' else v for j, v in zip('ACGT', self.ROWS[i])] for i in 'ACGT'} + with tempfile.TemporaryDirectory() as tmp: + path = os.path.join(tmp, 'model.train') + self.write_train(path, 43, 1, self.ROWS) + self.assertEqual(parse_train_file(path)['GAP-EXTEND'], 1) + for gap_open, gap_extend, rows in [(43, 0, self.ROWS), (-1, 1, self.ROWS), (43, 1, zero_at)]: + self.write_train(path, gap_open, gap_extend, rows) + with self.assertRaises(RuntimeError): + parse_train_file(path) + if __name__ == '__main__': unittest.main() \ No newline at end of file diff --git a/src/cactus/setup/cactus_align.py b/src/cactus/setup/cactus_align.py index bce3e5f3b..cce9676bc 100644 --- a/src/cactus/setup/cactus_align.py +++ b/src/cactus/setup/cactus_align.py @@ -482,7 +482,7 @@ def cactus_align(job, config_wrapper, mc_tree, input_seq_map, input_seq_id_map, scores_path = os.path.join(job.fileStore.getLocalTempDir(), 'last.train') job.fileStore.readGlobalFile(scores_id, scores_path) score_dict = parse_train_file(scores_path) - apply_scores_to_config(score_dict, config_wrapper.xmlRoot) + apply_scores_to_config(score_dict, config_wrapper.xmlRoot, label=chrom_name) # unzip the input sequences and enforce unique header prefixes sanitize_job = head_job.addChildJobFn(sanitize_fasta_headers, input_seq_id_map, pangenome=doVG or doGFA or do_filter_paf, walltime=cactus_walltime()) diff --git a/submodules/abPOA b/submodules/abPOA index 36110ac80..12f99c4bc 160000 --- a/submodules/abPOA +++ b/submodules/abPOA @@ -1 +1 @@ -Subproject commit 36110ac805f6f2d250621f8a08e5ef092e09d48b +Subproject commit 12f99c4bc2a9f92546d583e111bd15245a21da23