Stop a 'within' mutation from steering a neighboring deletion in NORMALIZE

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<string> seq_ids = verify_ref_seq_info.seq_ids();
-		vector<string> 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<string>::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<string> seq_ids = verify_ref_seq_info.seq_ids();
+
+			for (vector<string>::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).

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=<MOB id>' 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);

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

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

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=<MOB id>'
+# 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}

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

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}