Repository: cogent3/cogent3
Pull request #1824 (2024-04-24 11:24:08)

MAINT: use IndelMap for dotplotting
[CHANGED] reimplement the get_align_coords() function, now simpler
    and more robust

[NEW] added IndelMap.get_gap_align_coordinates() returns a 2D array
    of gap start, end in alignment coordinates. Used by get_align_coords().


Patch:
diff --git a/requirements.txt b/requirements.txt
index 8b70ef5b3..091cbe1fd 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -8,6 +8,6 @@ numba>0.54; python_version>='3.10'
 numba>=0.57.0; python_version=='3.11'
 scipy==1.13.0
 nox==2024.4.15
-plotly==5.20.0
+plotly==5.21.0
 typing_extensions>=3.7
 .[dev]
diff --git a/src/cogent3/core/location.py b/src/cogent3/core/location.py
index 4f8f274c6..a19bf5170 100644
--- a/src/cogent3/core/location.py
+++ b/src/cogent3/core/location.py
@@ -1697,6 +1697,11 @@ class IndelMap(MapABC):
         cum_lengths[1:] = diffs
         return numpy.array([self.gap_pos, cum_lengths]).T.tolist()
 
+    def get_gap_align_coordinates(self) -> NDArray[int]:
+        """returns [(gap start, gap end), ...] in alignment indices"""
+        starts, ends = _gap_spans(self.gap_pos, self.cum_gap_lengths)
+        return numpy.array([starts, ends]).T
+
     def merge_maps(self, other, parent_length: Optional[int] = None):
         """merge gaps of other with self
 
diff --git a/src/cogent3/draw/dotplot.py b/src/cogent3/draw/dotplot.py
index c08795777..6c0ecf73f 100644
--- a/src/cogent3/draw/dotplot.py
+++ b/src/cogent3/draw/dotplot.py
@@ -1,3 +1,5 @@
+import numpy
+
 from cogent3.align.pycompare import (
     MatchedSeqPaths,
     SeqKmers,
@@ -27,9 +29,23 @@ def suitable_threshold(window, desired_probability):
     return matches
 
 
-def not_gap(span):
-    """whether a span corresponds to a non-gap"""
-    return span.num_gaps == 0
+def _ungapped_spans(
+    starts_ends: list[list[int, int]], align_length: int, aligned=False
+) -> numpy.ndarray:
+    """returns numpy array of [(non gap start, non gap stop), ...}"""
+    if len(starts_ends) == 0:
+        data = [(0, align_length)] if aligned else []
+        return numpy.array(data, dtype=int)
+
+    if starts_ends[0][0] != 0:
+        # does not start with a gap, so we adjust to get an ungapped span at
+        # the start
+        starts_ends = [(0, 0)] + starts_ends
+    if starts_ends[-1][1] != align_length:
+        # does not end with a gap, so we adjust to get an ungapped span at
+        # the end
+        starts_ends.append((align_length, align_length))
+    return numpy.array(starts_ends).flatten()[1:-1].reshape((len(starts_ends) - 1, 2))
 
 
 def _convert_input(seq, moltype):
@@ -54,43 +70,24 @@ def get_align_coords(map1, map2, aligned=False) -> MatchedSeqPaths:
     """sequence coordinates of aligned segments"""
     from cogent3.align.pycompare import segment
 
-    paths = MatchedSeqPaths("Alignment")
-    if not_gap(map1) and not_gap(map2):
-        # no gaps
-        if aligned:
-            assert len(map1) == len(map2), "Aligned sequences inconsistent length"
-            paths[0].append((segment(0, len(map1)), segment(0, len(map2))))
-        return paths
-
     assert len(map1) == len(map2), "aligned sequences not equal length"
-    # diagonals are places where both sequences are NOT gaps
-    # so we record start of a diagonal and when we hit a 'gap'
-    # in either sequence, we have the end of the diagonal
-
-    start_x = start_y = None
-    for i in range(len(map1)):
-        x_not_gap = not_gap(map1[i])
-        y_not_gap = not_gap(map2[i])
-        if x_not_gap and y_not_gap and start_x is None:
-            start_x = map1[:i].parent_length
-            start_y = map2[:i].parent_length
-        elif (not x_not_gap or not y_not_gap) and start_x is not None:
-            paths[start_y - start_x].append(
-                (
-                    segment(start_x, map1[:i].parent_length - 1),
-                    segment(start_y, map2[:i].parent_length - 1),
-                )
-            )
-            start_x = start_y = None
-
-    if start_x is not None:
-        paths[start_y - start_x].append(
-            (
-                segment(start_x, map1.parent_length - 1),
-                segment(start_y, map2.parent_length - 1),
-            )
-        )
 
+    # we get the gap coordinates in alignment indices for both sequences
+    # sorting this allows us to trivially identify the alignment indices of
+    # ungapped segments,
+    # e.g. [(gap start 1, gap end 1), (gap start 2, gap end 2),...
+    # then [(gap end 1, gap start 2), ...
+    # is what we need to identify the segments for plotting
+    starts_ends = map1.get_gap_align_coordinates().tolist()
+    starts_ends.extend(map2.get_gap_align_coordinates().tolist())
+    starts_ends.sort()
+    starts_ends = _ungapped_spans(starts_ends, len(map1), aligned=aligned)
+    paths = MatchedSeqPaths("Alignment")
+    for align_start, align_end in starts_ends:
+        # ends are inclusive for plotting
+        s1 = segment(map1.get_seq_index(align_start), map1.get_seq_index(align_end) - 1)
+        s2 = segment(map2.get_seq_index(align_start), map2.get_seq_index(align_end) - 1)
+        paths.append(s1, s2)
     return paths
 
 


Test patch:
diff --git a/.github/workflows/testing_develop.yml b/.github/workflows/testing_develop.yml
index 4445bfdeb..1135ecc86 100644
--- a/.github/workflows/testing_develop.yml
+++ b/.github/workflows/testing_develop.yml
@@ -15,7 +15,7 @@ jobs:
 
     strategy:
       matrix:
-        os: [ubuntu-latest, macos-latest, windows-latest]
+        os: [ubuntu-latest, macos-13, windows-latest]
         python-version: ["3.9", "3.10", "3.11", "3.12"]
 
     steps:
diff --git a/tests/test_draw/test_dotplot.py b/tests/test_draw/test_dotplot.py
index 19180a7c0..334961bf9 100644
--- a/tests/test_draw/test_dotplot.py
+++ b/tests/test_draw/test_dotplot.py
@@ -1,13 +1,11 @@
 from unittest import TestCase
 
+import numpy
+
 from cogent3 import DNA, make_unaligned_seqs
 from cogent3.core.alignment import Aligned, ArrayAlignment
-from cogent3.draw.dotplot import (
-    Dotplot,
-    _convert_input,
-    get_align_coords,
-    not_gap,
-)
+from cogent3.core.location import IndelMap
+from cogent3.draw.dotplot import Dotplot, _convert_input, get_align_coords
 
 
 class TestUtilFunctions(TestCase):
@@ -16,12 +14,6 @@ class TestUtilFunctions(TestCase):
         m, seq = DNA.make_seq("ACGGT--A").parse_out_gaps()
         self.assertEqual(m.parent_length, 6)
 
-    def test_not_gap(self):
-        """distinguishes Map instances that include gap or not"""
-        m, seq = DNA.make_seq("ACGGT--A").parse_out_gaps()
-        self.assertTrue(not_gap(m[0]))
-        self.assertFalse(not_gap(m[5]))
-
     def test_convert_input(self):
         """converts data for dotplotting"""
         m, seq = DNA.make_seq("ACGGT--A").parse_out_gaps()
@@ -56,7 +48,8 @@ class TestUtilFunctions(TestCase):
         m1, seq1 = DNA.make_seq("ACGGTTTA").parse_out_gaps()
         m2, seq2 = DNA.make_seq("GGGGTTTA").parse_out_gaps()
         paths = get_align_coords(m1, m2, aligned=True)
-        self.assertEqual(paths.get_coords(), ([0, len(seq1)], [0, len(seq1)]))
+        # display ranges are inclusive, thus length - 1
+        self.assertEqual(paths.get_coords(), ([0, len(seq1) - 1], [0, len(seq1) - 1]))
 
         # raises an exception if the Aligned seqs are different lengths
         m1, seq1 = DNA.make_seq("ACGGTTTA").parse_out_gaps()
@@ -66,6 +59,13 @@ class TestUtilFunctions(TestCase):
 
     def test_display2d(self):
         """correctly constructs a Display2d"""
+        #           111111
+        # 0123456789012345
+        #            11111
+        #  012345678901234
+        # -TGATGTAAGGTAGTT
+        # CTGG---AAG---GGT
+        # 0123   456   789
         # check alignment coords are correct
         dp = Dotplot("-TGATGTAAGGTAGTT", "CTGG---AAG---GGT", is_aligned=True, window=5)
         expect = [0, 2, None, 6, 8, None, 12, 14], [1, 3, None, 4, 6, None, 7, 9]
@@ -79,7 +79,7 @@ class TestUtilFunctions(TestCase):
         self.assertEqual(traces[-1].name, "Alignment")
         self.assertEqual(traces[0].name, "+ strand")
         # check + strand has integers/float/None
-        expect = {int, float, type(None)}
+        expect = {int, float, type(None), numpy.int64, numpy.int32}
         for trace in traces:
             for axis in "xy":
                 got = {type(v) for v in trace[axis]}
@@ -146,3 +146,23 @@ class TestUtilFunctions(TestCase):
         )
         dp = seqs.dotplot("seq1", "seq3", title="")
         self.assertEqual(dp.figure.layout.title, "")
+
+
+def test_aligned_path():
+    imap1 = IndelMap(
+        gap_pos=numpy.array([4, 5, 6, 8, 10]),
+        cum_gap_lengths=numpy.array([6, 10, 12, 14, 15]),
+        parent_length=10,
+    )
+    imap2 = IndelMap(
+        gap_pos=numpy.array([], dtype=int),
+        cum_gap_lengths=numpy.array([], dtype=int),
+        parent_length=25,
+    )
+
+    coords = get_align_coords(imap1, imap2)
+
+    assert coords.get_coords() == (
+        [0, 3, None, 4, 4, None, 5, 5, None, 6, 7, None, 8, 9],
+        [0, 3, None, 10, 10, None, 15, 15, None, 18, 19, None, 22, 23],
+    )