Cigar tools: fix edge cases of strip

Author: Donaim

micall/core/contig_stitcher.py

@@ -103,7 +103,7 @@ def lstrip_query(self) -> 'AlignedContig':
         alignment. The CIGAR alignment is also updated to reflect the trimming.
         """
 
-        alignment = self.alignment.lstrip_query()
+        alignment = self.alignment.lstrip_query().lstrip_reference()
         q_remainder, query = self.cut_query(alignment.q_st - 0.5)
         alignment = alignment.translate(0, -1 * alignment.q_st)
         result = AlignedContig.make(query, alignment, self.strand)
@@ -120,7 +120,7 @@ def rstrip_query(self) -> 'AlignedContig':
         alignment. The CIGAR alignment is also updated to reflect the trimming.
         """
 
-        alignment = self.alignment.rstrip_query()
+        alignment = self.alignment.rstrip_query().rstrip_reference()
         query, q_remainder = self.cut_query(alignment.q_ei + 0.5)
         result = AlignedContig.make(query, alignment, self.strand)
         logger.debug("Doing rstrip of %r resulted in %r, so %s (len %s) became %s (len %s)",

micall/tests/test_cigar_tools.py

@@ -343,9 +343,11 @@ def test_cigar_hit_ref_cut_add_prop_exhaustive(hit, cut_point):
     ('4D6I5M@1->1', '4D5M@7->1'),
     ('4I6D5M@1->1', '6D5M@5->1'),
     ('6I4D@1->1', '4D@7->1'),
-    ('6D4I@1->1', '6D4I@1->1'),
-    ('4D6I@1->1', '4D6I@1->1'),
+    ('6D4I@1->1', '6D@5->1'),
+    ('4D6I@1->1', '4D@7->1'),
     ('4I6D@1->1', '6D@5->1'),
+    ('4I@1->1', '@5->1'),
+    ('4D@1->1', '4D@1->1'),
     ('@1->1', '@1->1'),
 ]
 
@@ -364,9 +366,11 @@ def test_cigar_hit_ref_cut_add_prop_exhaustive(hit, cut_point):
     ('5M6D4I@1->1', '5M6D@1->1'),
     ('5M6I4D@1->1', '5M4D@1->1'),
     ('6D4I@1->1', '6D@1->1'),
-    ('6I4D@1->1', '6I4D@1->1'),
-    ('4I6D@1->1', '4I6D@1->1'),
+    ('6I4D@1->1', '4D@1->1'),
+    ('4I6D@1->1', '6D@1->1'),
     ('4D6I@1->1', '4D@1->1'),
+    ('4I@1->1', '@1->1'),
+    ('4D@1->1', '4D@1->1'),
     ('@1->1', '@1->1'),
 ]
 
@@ -384,10 +388,12 @@ def test_cigar_hit_ref_cut_add_prop_exhaustive(hit, cut_point):
     ('3I2D3I2D5M@1->1', '6I5M@1->5'),
     ('4D6I5M@1->1', '6I5M@1->5'),
     ('4I6D5M@1->1', '4I5M@1->7'),
-    ('6I4D@1->1', '6I4D@1->1'),
+    ('6I4D@1->1', '6I@1->5'),
     ('6D4I@1->1', '4I@1->7'),
     ('4D6I@1->1', '6I@1->5'),
-    ('4I6D@1->1', '4I6D@1->1'),
+    ('4I6D@1->1', '4I@1->7'),
+    ('4I@1->1', '4I@1->1'),
+    ('4D@1->1', '@1->5'),
     ('@1->1', '@1->1'),
 ]
 
@@ -405,10 +411,12 @@ def test_cigar_hit_ref_cut_add_prop_exhaustive(hit, cut_point):
     ('5M2D3I2D3I@1->1', '5M6I@1->1'),
     ('5M6D4I@1->1', '5M4I@1->1'),
     ('5M6I4D@1->1', '5M6I@1->1'),
-    ('6D4I@1->1', '6D4I@1->1'),
+    ('6D4I@1->1', '4I@1->1'),
     ('6I4D@1->1', '6I@1->1'),
     ('4I6D@1->1', '4I@1->1'),
-    ('4D6I@1->1', '4D6I@1->1'),
+    ('4D6I@1->1', '6I@1->1'),
+    ('4I@1->1', '4I@1->1'),
+    ('4D@1->1', '@1->1'),
     ('@1->1', '@1->1'),
 ]
 
@@ -493,8 +501,12 @@ def test_cigar_hit_reference_strip_is_idempotent(hit):
 def test_cigar_hit_reference_strips_are_commutative(hit):
     hit = parsed_hit(hit)
 
-    assert hit.rstrip_reference().lstrip_reference() \
-        == hit.lstrip_reference().rstrip_reference()
+    if len(hit.cigar.coordinate_mapping.ref_to_query) > 0:
+        assert hit.rstrip_reference().lstrip_reference() \
+            == hit.lstrip_reference().rstrip_reference()
+    else:
+        assert hit.rstrip_reference().lstrip_reference().cigar \
+            == hit.lstrip_reference().rstrip_reference().cigar
 
 
 @pytest.mark.parametrize('hit, expected', lstrip_query_cases)
@@ -540,13 +552,13 @@ def test_cigar_hit_query_strip_combines_with_add(hit):
 
 
 @pytest.mark.parametrize('hit', strip_prop_cases_all)
-def test_cigar_hit_query_strip_never_crashes(hit):
+def test_cigar_hit_strips_work_together(hit):
     hit = parsed_hit(hit)
 
-    hit.rstrip_query().lstrip_query()
-    hit.lstrip_query().rstrip_query()
-    hit.lstrip_query().lstrip_query()
-    hit.rstrip_query().rstrip_query()
+    rstrip = str(hit.rstrip_query().rstrip_reference().cigar)
+    assert not rstrip.endswith("I") and not rstrip.endswith("D")
+    lstrip = str(hit.lstrip_query().lstrip_reference().cigar)
+    assert not lstrip.startswith("I") and not lstrip.startswith("D")
 
 
 @pytest.mark.parametrize('hit', strip_prop_cases_all)
@@ -570,8 +582,12 @@ def test_cigar_hit_query_strip_is_idempotent(hit):
 def test_cigar_hit_query_strips_are_commutative(hit):
     hit = parsed_hit(hit)
 
-    assert hit.rstrip_query().lstrip_query() \
-        == hit.lstrip_query().rstrip_query()
+    if len(hit.cigar.coordinate_mapping.ref_to_query) > 0:
+        assert hit.rstrip_query().lstrip_query() \
+            == hit.lstrip_query().rstrip_query()
+    else:
+        assert hit.rstrip_query().lstrip_query().cigar \
+            == hit.lstrip_query().rstrip_query().cigar
 
 
 @pytest.mark.parametrize('hit, cut_point', [(x[0], x[1]) for x in cigar_hit_ref_cut_cases

micall/utils/cigar_tools.py

@@ -229,7 +229,7 @@ def slice_operations(self, start_inclusive, end_noninclusive) -> 'Cigar':
     def lstrip_reference(self) -> 'Cigar':
         """ Return a copy of the Cigar with leading (unmatched) reference elements removed. """
 
-        min_r = min(self.coordinate_mapping.ref_to_query.keys(), default=0)
+        min_r = min(self.coordinate_mapping.ref_to_query.keys(), default=None)
         min_op = self.coordinate_mapping.ref_to_op.get(min_r, float("inf"))
 
         ops = [(1, op) for i, (op, ref_pointer, query_pointer)
@@ -241,9 +241,8 @@ def lstrip_reference(self) -> 'Cigar':
     def rstrip_reference(self) -> 'Cigar':
         """ Return a copy of the Cigar with trailing (unmatched) reference elements removed. """
 
-        max_r = max(self.coordinate_mapping.ref_to_query.keys(),
-                    default=len(self.coordinate_mapping.ref_to_op) - 1)
-        max_op = self.coordinate_mapping.ref_to_op.get(max_r, float("inf"))
+        max_r = max(self.coordinate_mapping.ref_to_query.keys(), default=None)
+        max_op = self.coordinate_mapping.ref_to_op.get(max_r, float("-inf"))
 
         ops = [(1, op) for i, (op, ref_pointer, query_pointer)
                in enumerate(self.iterate_operations_with_pointers())
@@ -254,7 +253,7 @@ def rstrip_reference(self) -> 'Cigar':
     def lstrip_query(self) -> 'Cigar':
         """ Return a copy of the Cigar with leading (unmatched) query elements removed. """
 
-        min_q = min(self.coordinate_mapping.query_to_ref.keys(), default=0)
+        min_q = min(self.coordinate_mapping.query_to_ref.keys(), default=None)
         min_op = self.coordinate_mapping.query_to_op.get(min_q, float("inf"))
 
         ops = [(1, op) for i, (op, ref_pointer, query_pointer)
@@ -266,9 +265,8 @@ def lstrip_query(self) -> 'Cigar':
     def rstrip_query(self) -> 'Cigar':
         """ Return a copy of the Cigar with trailing (unmatched) query elements removed. """
 
-        max_q = max(self.coordinate_mapping.query_to_ref.keys(),
-                    default=len(self.coordinate_mapping.query_to_op) - 1)
-        max_op = self.coordinate_mapping.query_to_op.get(max_q, float("inf"))
+        max_q = max(self.coordinate_mapping.query_to_ref.keys(), default=None)
+        max_op = self.coordinate_mapping.query_to_op.get(max_q, float("-inf"))
 
         ops = [(1, op) for i, (op, ref_pointer, query_pointer)
                in enumerate(self.iterate_operations_with_pointers())