From c4b2d91441525da3138bda38f323e995c9def27a Mon Sep 17 00:00:00 2001 From: Jeffrey Barrick Date: Sat, 22 Aug 2026 18:55:16 -0400 Subject: [PATCH] Stop a 'within' mutation from steering a neighboring deletion in NORMALIZE A mutation tagged within= has a POSITION in a pseudo-coordinate space -- the MOB's position plus an offset into the newly inserted element -- so it occupies zero reference bases while freely colliding with real reference coordinates downstream of the MOB. normalize_and_annotate_tandem_repeat_mutations skipped such a mutation from being shifted, but its 'goto next_mutation' still assigned it to last_mut, so that pseudo-coordinate bounded the next deletion's shift. On a real curated LTEE clone this moved DEL 1270157 535 to 1270063, 94 bases outside its equivalence window [1270157, 1270234], and the output then failed APPLY because the deletion had swallowed the MOB it was backing off from. The sister clone, whose only difference is not carrying that within SNP, normalized the same authored line to 1270234 -- one shared ancestral event with two coordinates, which reads downstream as a homoplasy that never happened. Four changes: - A 'within' mutation is now skipped entirely and never becomes last_mut. The mutation it is within is itself in the list and supplies the real barrier at real coordinates. 'no_normalize' keeps the old path: those positions are genuine reference coordinates. - The DEL and INS back-off fires only when the shift is what created the overlap. The comment always said "did we get shifted into the position of the next mutation?", but the test was plain interval overlap, evaluated even when nothing moved. An interval that already covered last_mut in the input got pushed LEFT of where the author put it, by a distance that tracks last_mut's coordinate. Requiring the pre-shift interval to end before last_mut also bounds the result from below. - The APPLY self-check applies the normalized diff instead of re-parsing the input file. It was comparing the input against itself, so "Failed APPLY test" was unreachable and the non-zero exit added alongside it could never fire. It now also re-validates the output, which is what names this failure precisely: a shift can make the diff ambiguous even when the two applied sequences still agree. - NORMALIZE validates its input the way APPLY does. Normalizing an ambiguous diff is not meaningful -- the orderings give different genomes and normalization has to pick one -- and accepting it silently is how a bad coordinate got minted from a bad input. -x skips both checks. Adds the first two tests for NORMALIZE: gdtools_normalize_1 pins the within case (it emits DEL 19800, 230 bases left, without these changes) and gdtools_normalize_2 asserts the ambiguous input is rejected. Co-Authored-By: Claude Opus 5 (1M context) --- src/breseq/gdtools_cmdline.cpp | 53 ++++++++++++++++++++------- src/breseq/mutation_predictor.cpp | 40 +++++++++++++++++--- tests/gdtools_normalize_1/expected.gd | 6 +++ tests/gdtools_normalize_1/input.gd | 5 +++ tests/gdtools_normalize_1/testcmd.sh | 30 +++++++++++++++ tests/gdtools_normalize_2/input.gd | 4 ++ tests/gdtools_normalize_2/testcmd.sh | 22 +++++++++++ 7 files changed, 140 insertions(+), 20 deletions(-) create mode 100644 tests/gdtools_normalize_1/expected.gd create mode 100644 tests/gdtools_normalize_1/input.gd create mode 100755 tests/gdtools_normalize_1/testcmd.sh create mode 100644 tests/gdtools_normalize_2/input.gd create mode 100755 tests/gdtools_normalize_2/testcmd.sh diff --git a/src/breseq/gdtools_cmdline.cpp b/src/breseq/gdtools_cmdline.cpp index 299ac465..b4f3e4ed 100644 --- a/src/breseq/gdtools_cmdline.cpp +++ b/src/breseq/gdtools_cmdline.cpp @@ -1517,7 +1517,7 @@ int do_normalize_gd(int argc, char* argv[]) options("reference,r" , "File containing reference sequences in GenBank, GFF3, or FASTA format. Option may be provided multiple times for multiple files (REQUIRED)"); options("reassign-ids,s" , "reassign ids to lowest numbers possible.", TAKES_NO_ARGUMENT); options("repeat-adjacent,a" , "mark repeat-region adjacent, mediated, and between mutations.", TAKES_NO_ARGUMENT); - options("dont-check-apply,x" , "skip step that checks consistency of normalize using APPLY.", TAKES_NO_ARGUMENT); + options("dont-check-apply,x" , "skip both the check that the input is valid against the reference sequences and the step that checks consistency of normalize using APPLY.", TAKES_NO_ARGUMENT); const int32_t kDistanceToRepeat = 20; @@ -1581,7 +1581,16 @@ int do_normalize_gd(int argc, char* argv[]) ref_seq_info.LoadFiles(reference_file_names); cGenomeDiff gd(input); - + + // Normalizing an ambiguous Genome Diff is not meaningful: when one mutation overlaps bases + // another one deletes or duplicates, the two can be applied in either order to different + // results, and normalization has to pick one. APPLY rejects these outright, so NORMALIZE + // does too rather than silently shifting coordinates based on a reading APPLY will not share. + // '--dont-check-apply' skips this along with the sequence check below. + if (!options.count("dont-check-apply")) { + gd.valid_with_reference_sequences(ref_seq_info); + } + Settings settings; cReferenceSequences new_ref_seq_info; @@ -1629,28 +1638,44 @@ int do_normalize_gd(int argc, char* argv[]) gd.reassign_unique_ids(); } + uout("Writing output Genome Diff file", options["output"]); + gd.write(options["output"]); + bool apply_test_failed = false; if (!options.count("dont-check-apply")) { uout("Using APPLY to check that normalization didn't change the mutated sequence."); cReferenceSequences verify_ref_seq_info = cReferenceSequences::deep_copy(ref_seq_info); - cGenomeDiff verify_gd(input); // must load new copy or positions will be shifted by apply_to_sequences - verify_gd.apply_to_sequences(ref_seq_info, verify_ref_seq_info, false, kDistanceToRepeat, settings.size_cutoff_AMP_becomes_INS_DEL_mutation); - vector seq_ids = verify_ref_seq_info.seq_ids(); - vector new_seq_ids = new_ref_seq_info.seq_ids(); + // Read back what we just wrote rather than reusing 'gd': apply_to_sequences shifts + // positions in place, and cGenomeDiff holds its entries by shared_ptr so copying it would + // alias them. Loading the output file also means the check covers exactly the bytes the + // caller gets. Note this must be the NORMALIZED diff -- applying 'input' here compares + // the input against itself, which is a tautology that can never fail. + cGenomeDiff verify_gd(options["output"]); + + // The shift can move a mutation onto bases another mutation deletes or duplicates, which + // makes the output ambiguous even when the two applied sequences still agree. The input + // was already checked for this above, so any such error is one normalization introduced. + cFileParseErrors verify_parse_errors = verify_gd.valid_with_reference_sequences(verify_ref_seq_info, true); + if (verify_parse_errors._errors.size()) { + WARN("Failed APPLY test. NORMALIZE produced a Genome Diff that is no longer valid against the reference sequences."); + verify_parse_errors.print_errors(false); + apply_test_failed = true; + } else { + verify_gd.apply_to_sequences(ref_seq_info, verify_ref_seq_info, false, kDistanceToRepeat, settings.size_cutoff_AMP_becomes_INS_DEL_mutation); - for (vector::const_iterator it = seq_ids.begin(); it != seq_ids.end(); it++) - { - if (new_ref_seq_info[*it].m_fasta_sequence.get_sequence() != verify_ref_seq_info[*it].m_fasta_sequence.get_sequence()) { - WARN("Failed APPLY test. Discrepancies between sequences produced before and after NORMALIZE. Check ordering of mutations."); - apply_test_failed = true; + vector seq_ids = verify_ref_seq_info.seq_ids(); + + for (vector::const_iterator it = seq_ids.begin(); it != seq_ids.end(); it++) + { + if (new_ref_seq_info[*it].m_fasta_sequence.get_sequence() != verify_ref_seq_info[*it].m_fasta_sequence.get_sequence()) { + WARN("Failed APPLY test. Discrepancies between sequences produced before and after NORMALIZE. Check ordering of mutations."); + apply_test_failed = true; + } } } } - uout("Writing output Genome Diff file", options["output"]); - gd.write(options["output"]); - // A failed self-consistency check means NORMALIZE altered the mutated sequence, so the // output cannot be trusted. Return non-zero even though output was written (the caller // discards the declared output on a non-zero exit). diff --git a/src/breseq/mutation_predictor.cpp b/src/breseq/mutation_predictor.cpp index 778f72fe..a48b5360 100644 --- a/src/breseq/mutation_predictor.cpp +++ b/src/breseq/mutation_predictor.cpp @@ -2593,9 +2593,16 @@ namespace breseq { // We are still potentially in danger of doing the wrong thing here, // because a mutation could be applied only after one with the 'before' tag, making the shift // incorrect. So, setting 'no_normalize' is an out that can be used. - - if (mut.entry_exists("within") || mut.entry_exists("no_normalize") ) goto next_mutation; - + + // A 'within' mutation is skipped entirely -- it is not shifted, and it does not become + // 'last_mut'. Its POSITION is not a reference coordinate: for 'within=' it is the + // MOB's position plus an offset into the newly inserted element, so it names bases that do + // not exist in the reference and cannot legitimately bound a neighbor's shift. The mutation + // it is within is itself in this list and provides the real barrier at its real coordinates. + if (mut.entry_exists("within")) continue; + + if (mut.entry_exists("no_normalize")) goto next_mutation; + if (mut._type == INS) { int32_t size = mut["new_seq"].size(); @@ -2621,7 +2628,13 @@ namespace breseq { // Did we get shifted into the position of the next mutation? Then back off // Note: We don't do this with converted AMPs as this creates problems (an insertion within them can shift their position) - if (last_mut && (mut[SEQ_ID] == last_mut->get(SEQ_ID)) && (mut.get_reference_coordinate_end() >= last_mut->get_reference_coordinate_start()) && !mut.entry_exists("_was_AMP")) { + // As in the DEL case below, only back off when the shift is what created the overlap -- + // an insertion that already sat at or past last_mut in the input is not ours to move. + bool ins_overlap_existed_before_shift = + last_mut && (mut[SEQ_ID] == last_mut->get(SEQ_ID)) + && (cReferenceCoordinate(position, insert_position) >= last_mut->get_reference_coordinate_start()); + + if (last_mut && !ins_overlap_existed_before_shift && (mut[SEQ_ID] == last_mut->get(SEQ_ID)) && (mut.get_reference_coordinate_end() >= last_mut->get_reference_coordinate_start()) && !mut.entry_exists("_was_AMP")) { // The position of this insert mutation should be one before the mutation, unless it is another INS, // In the INS case, we need to properly update all of the insert positions @@ -2703,10 +2716,25 @@ namespace breseq { // Begin consensus mode shifting of coordinates ----> if (!settings.polymorphism_prediction) { + // Remember where the author put it, so we can tell a shift-induced overlap from one + // that was already in the input. + int32_t position_before_shift = position; + normalizeDELposition(ref_seq_info[mut["seq_id"]], mut, repeat_unit_sequence); - + // Did we get shifted into the position of the next mutation? Then back off. - if (last_mut && !gd.applied_before_id(last_mut->_id, mut._id) && (mut[SEQ_ID] == last_mut->get(SEQ_ID)) && (mut.get_reference_coordinate_end() >= last_mut->get_reference_coordinate_start())) { + // Only if the shift is what created the overlap: a deletion that already covered + // last_mut in the input is not ours to move. Backing off from that pushes the deletion + // LEFT of where the author put it, onto bases it was never equivalent to, and the + // distance moved tracks last_mut's coordinate -- so the same authored deletion lands + // somewhere different in every sample that happens to carry a mutation inside it. + // Requiring the original interval to end before last_mut also bounds the result from + // below: position_before_shift <= last_mut position - size whenever this fires. + bool overlap_existed_before_shift = + last_mut && (mut[SEQ_ID] == last_mut->get(SEQ_ID)) + && (cReferenceCoordinate(position_before_shift + size - 1) >= last_mut->get_reference_coordinate_start()); + + if (last_mut && !overlap_existed_before_shift && !gd.applied_before_id(last_mut->_id, mut._id) && (mut[SEQ_ID] == last_mut->get(SEQ_ID)) && (mut.get_reference_coordinate_end() >= last_mut->get_reference_coordinate_start())) { // The position of this insert mutation should be as many bases as the // deletion is long before the next mutation mut["position"] = s(n(last_mut->get("position")) - size); diff --git a/tests/gdtools_normalize_1/expected.gd b/tests/gdtools_normalize_1/expected.gd new file mode 100644 index 00000000..aab1e9ef --- /dev/null +++ b/tests/gdtools_normalize_1/expected.gd @@ -0,0 +1,6 @@ +#=GENOME_DIFF 1.0 +#=TITLE input +#=REFSEQ REL606.fragment.gbk +MOB 1 . REL606-5 20000 IS1 1 9 repeat_size=768 +DEL 2 . REL606-5 20032 500 +SNP 3 . REL606-5 20300 A within=1 diff --git a/tests/gdtools_normalize_1/input.gd b/tests/gdtools_normalize_1/input.gd new file mode 100644 index 00000000..c5429530 --- /dev/null +++ b/tests/gdtools_normalize_1/input.gd @@ -0,0 +1,5 @@ +#=GENOME_DIFF 1.0 +#=REFSEQ REL606.fragment.gbk +MOB 1 . REL606-5 20000 IS1 1 9 +DEL 2 . REL606-5 20030 500 +SNP 3 . REL606-5 20300 A within=1 diff --git a/tests/gdtools_normalize_1/testcmd.sh b/tests/gdtools_normalize_1/testcmd.sh new file mode 100755 index 00000000..072e98a3 --- /dev/null +++ b/tests/gdtools_normalize_1/testcmd.sh @@ -0,0 +1,30 @@ +#!/bin/bash + +SELF=`dirname ${BASH_SOURCE}` +. ${SELF}/../common.sh + +# NORMALIZE must not let a 'within' mutation steer a neighboring deletion's coordinate. +# +# input.gd puts a SNP inside the new IS1 copy created by the MOB (within=1). For 'within=' +# the POSITION field is not a reference coordinate -- it is the MOB's position plus an offset into +# the inserted element -- so that SNP occupies no reference bases and cannot legitimately bound the +# DEL. It nonetheless falls inside the DEL's reference interval (20030..20529), which used to send +# the DEL to 20300-500 = 19800: 230 bases left of where the author put it, deleting sequence it was +# never equivalent to, and tracking the SNP's coordinate so the same authored DEL landed elsewhere +# in every sample that happened to carry a mutation inside it. +# +# Expected: the DEL only right-shifts within its equivalence window (20030 -> 20032), and NORMALIZE +# exits 0 -- the post-normalization validity check would otherwise reject the output. + +CURRENT_OUTPUTS[0]="${SELF}/output.gd" +EXPECTED_OUTPUTS[0]="${SELF}/expected.gd" + +TESTCMD="\ + ${GDTOOLS} \ + NORMALIZE \ + -o ${SELF}/output.gd \ + -r ${DATADIR}/REL606/REL606.fragment.gbk \ + ${SELF}/input.gd \ + " + +do_test $1 ${SELF} diff --git a/tests/gdtools_normalize_2/input.gd b/tests/gdtools_normalize_2/input.gd new file mode 100644 index 00000000..1494ff5d --- /dev/null +++ b/tests/gdtools_normalize_2/input.gd @@ -0,0 +1,4 @@ +#=GENOME_DIFF 1.0 +#=REFSEQ REL606.fragment.gbk +DEL 1 . REL606-5 20030 500 +SNP 2 . REL606-5 20300 A diff --git a/tests/gdtools_normalize_2/testcmd.sh b/tests/gdtools_normalize_2/testcmd.sh new file mode 100755 index 00000000..f59ac2d0 --- /dev/null +++ b/tests/gdtools_normalize_2/testcmd.sh @@ -0,0 +1,22 @@ +#!/bin/bash + +SELF=`dirname ${BASH_SOURCE}` +. ${SELF}/../common.sh + +# Failure test: NORMALIZE must reject an ambiguous Genome Diff, as APPLY already does. +# +# input.gd has a plain SNP sitting on bases the DEL removes, with no 'within' or 'before' to say +# which applies first. The two orderings give different genomes, so there is no single coordinate +# to normalize to. NORMALIZE used to accept this and shift the DEL anyway, which produced a +# silently different deletion; it now fails with the same error APPLY gives. +EXPECTED_EXIT_CODE=1 + +TESTCMD="\ + ${GDTOOLS} \ + NORMALIZE \ + -o ${SELF}/output.gd \ + -r ${DATADIR}/REL606/REL606.fragment.gbk \ + ${SELF}/input.gd \ + " + +do_test $1 ${SELF}