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}