Repository: cogent3/cogent3
Pull request #2259 (2025-03-22 00:17:38)

BUG: fixed a significant error in usage of IndelMap.merge_maps
[FIXED] The alignment code applies a merge operation to IndelMap objects
    progressively as it aligns sequences. This means the incoming IndelMap
    object is in alignment coordinates, not the original sequence
    coordinates. Introduced the aligned_indices parameter to IndelMap.merge_maps
    to handle this case.

    This will have caused significant errors in the output of the cogent3
    progressive aligners.

## Summary by Sourcery

Fixes a bug in the progressive aligner related to merging gaps in IndelMap objects during progressive alignment. Introduces the `aligned_indices` parameter to `IndelMap.merge_maps` to correctly handle alignment coordinates. Adds new tests to verify the fix.

Bug Fixes:
- Fixes a major bug in the cogent3 progressive aligner where gaps were being incorrectly merged during progressive steps. Users are advised to re-run any analyses that used the progressive aligner. This does not affect creating alignment objects from existing data, e.g. alignments produced by other algorithms.
- Fixes incorrect merging of gaps during progressive alignment by introducing the aligned_indices parameter to IndelMap.merge_maps to handle alignment coordinates correctly.

Enhancements:
- IndelMap.merge_maps(other: "IndelMap", alignd_indices: bool) handles case where `other.gap_pos` are in alignment coordinates.

Tests:
- Adds tests for indelmap merge with ending gap and aligned indices.


Patch:
diff --git a/changelog.d/20250314_131212_Gavin.Huttley.md b/changelog.d/20250314_131212_Gavin.Huttley.md
index d06efdc24..53ec4981f 100644
--- a/changelog.d/20250314_131212_Gavin.Huttley.md
+++ b/changelog.d/20250314_131212_Gavin.Huttley.md
@@ -8,6 +8,7 @@ Uncomment the section that is right (remove the HTML comment wrapper).
 
 - First time contributions from @c-jonesy, @Yixuan567, @Qinzi27
   and @je634 👏🚀🎉
+- @GavinHuttley bug fixes and maintenance
 
 ### ENH
 
@@ -30,11 +31,20 @@ Uncomment the section that is right (remove the HTML comment wrapper).
   For example, `use_hook="my_package"` would invoke the one defined above. Set `use_hook="cogent3"`
   to force using the cogent3 default.
 
+- IndelMap.merge_maps(other: "IndelMap", alignd_indices: bool) handles case where
+  `other.gap_pos` are in alignment coordinates. This is essential for the case
+  of progressive alignments where the alignment coordinates are sequentially
+  adjusted as more profiles are aligned. This is the key change that fixes the
+  bug in the progressive aligner.
 
 ### BUG
 
 - Dotplot on new_type Alignments now show the alignment path 🎉
-
+- Fixed a MAJOR bug affecting the cogent3 progressive aligner. Gaps
+  were being incorrectly merged during the progressive steps. Users
+  are advised to re-run any analyses that used the progressive aligner.
+  This DOES NOT affect creating alignment objects from existing data,
+  e.g. alignments produced by other algorithms.
 
 <!--
 ### DOC
@@ -54,5 +64,3 @@ Uncomment the section that is right (remove the HTML comment wrapper).
 - new_alignment.Alignment.quick_tree() no longer supports doing bootstraps.
   This is because bootstrapping is not quick! If you want this, post a feature
   request and we will implement a separate app for this.
-
-
diff --git a/src/cogent3/align/pairwise.py b/src/cogent3/align/pairwise.py
index 3e4420792..33ad70352 100644
--- a/src/cogent3/align/pairwise.py
+++ b/src/cogent3/align/pairwise.py
@@ -494,6 +494,7 @@ class AlignablePOG(_Alignable):
                 new_map = aligned.map.merge_maps(
                     maps[dim] * word_length,
                     parent_length=maps[dim].parent_length * word_length,
+                    aligned_indices=True,
                 )
                 # Likewise, if the data is not modulo word_length,
                 # it is trimmed to match
diff --git a/src/cogent3/core/location.py b/src/cogent3/core/location.py
index f014f5449..4f759c2c2 100644
--- a/src/cogent3/core/location.py
+++ b/src/cogent3/core/location.py
@@ -11,6 +11,7 @@ from typing import Any, NoReturn, Optional, Union
 
 import numba
 import numpy
+import typing_extensions
 from numpy.typing import NDArray
 
 from cogent3._version import __version__
@@ -1383,6 +1384,7 @@ class IndelMap(MapABC):
         self,
         other: "IndelMap",
         parent_length: int | None = None,
+        aligned_indices: bool = False,
     ) -> "IndelMap":
         """merge gaps of other with self
 
@@ -1392,13 +1394,21 @@ class IndelMap(MapABC):
             instance for same sequence
         parent_length
             overrides property
+        aligned_indices
+            other.gap_pos are alignment indices of self
         """
-        unique_pos = numpy.union1d(self.gap_pos, other.gap_pos)
+        if aligned_indices:
+            other_pos, other_lengths = self._aligned_to_seq_map(other)
+        else:
+            other_pos = other.gap_pos
+            other_lengths = other.get_gap_lengths()
+
+        unique_pos = numpy.union1d(self.gap_pos, other_pos)
         gap_lengths = numpy.zeros(unique_pos.shape, dtype=_DEFAULT_GAP_DTYPE)
         self_lengths = self.get_gap_lengths()
-        other_lengths = other.get_gap_lengths()
+
         _update_lengths(unique_pos, gap_lengths, self.gap_pos, self_lengths)
-        _update_lengths(unique_pos, gap_lengths, other.gap_pos, other_lengths)
+        _update_lengths(unique_pos, gap_lengths, other_pos, other_lengths)
         parent_length = parent_length or self.parent_length
         return self.__class__(
             gap_pos=unique_pos,
@@ -1406,6 +1416,22 @@ class IndelMap(MapABC):
             parent_length=parent_length,
         )
 
+    def _aligned_to_seq_map(
+        self,
+        other: typing_extensions.Self,
+    ) -> tuple[numpy.ndarray[int], numpy.ndarray[int]]:
+        # converts the gap pos of other into seq indices,
+        # returning unique positions and their summed lengths
+        other_pos = numpy.array(
+            [self.get_seq_index(p) for p in other.gap_pos],
+            dtype=_DEFAULT_GAP_DTYPE,
+        )
+        other_lengths = other.get_gap_lengths()
+        # reduce to unique positions
+        other_pos, unique_indices = numpy.unique(other_pos, return_index=True)
+        other_lengths = numpy.add.reduceat(other_lengths, unique_indices)
+        return other_pos, other_lengths
+
     def joined_segments(self, coords: SeqCoordTypes) -> "IndelMap":
         """returns new map with disjoint gapped segments joined
 


Test patch:
diff --git a/tests/test_app/test_align.py b/tests/test_app/test_align.py
index 9da8d46ea..ecf672da8 100644
--- a/tests/test_app/test_align.py
+++ b/tests/test_app/test_align.py
@@ -803,3 +803,29 @@ def test_codon_incomplete(DATA_DIR):
     # now make sure the resulting ungapped sequences are modulo 3
     seqs = aln.degap().to_dict().values()
     assert {len(s) % 3 for s in seqs} == {0}
+
+
+def test_align_with_gap_merging_case():
+    tree = "((RoundEare:0.165,FalseVamp:0.061):0.057,LeafNose:0.098,Sloth:0.044);"
+    aln = make_aligned_seqs(
+        {
+            "LeafNose": "AGGGTCACCAAAGACGGGCCAGGAAATTAGAGTCCTCAGAAGAGGA-------------",
+            "FalseVamp": "AGGGTCACCGAAGA---------------------------------------------",
+            "RoundEare": "GGGGTCACCCCAGAAGGGCCAGGAAACTCGAGTCCTCAGAAGAGGA-------------",
+            "Sloth": "--------AAAAGACGGGCCAGGAAATTAGACTCCTCAGAAGAGGATGTGTCTAGTGAG",
+        },
+        moltype="dna",
+    )
+    seqs = aln.degap()
+    # trigger an alignment with many gaps using a lower
+    # gap insertion penalty than the gap extension penalty
+    aligner = get_app(
+        "progressive_align",
+        "nucleotide",
+        indel_rate=1e-2,
+        indel_length=1e-9,
+        guide_tree=tree,
+    )
+    aln = aligner.main(seqs)
+    # does not fail
+    assert len(aln)
diff --git a/tests/test_core/test_location.py b/tests/test_core/test_location.py
index 9ba149bfb..1328e8c62 100644
--- a/tests/test_core/test_location.py
+++ b/tests/test_core/test_location.py
@@ -1519,3 +1519,44 @@ def test_indelmap_scaled(codon_and_aa_maps):
     amap_to_cmap = amap.scaled(3)
     numpy.testing.assert_allclose(cmap_to_amap.array, amap.array)
     numpy.testing.assert_allclose(amap_to_cmap.array, cmap.array)
+
+
+def test_indelmap_merge_ending_gap():
+    im1 = IndelMap(
+        gap_pos=numpy.array([5], dtype=int),
+        gap_lengths=numpy.array([6], dtype=int),
+        parent_length=5,
+    )
+    im2 = IndelMap(
+        gap_pos=numpy.array([11], dtype=int),
+        gap_lengths=numpy.array([2], dtype=int),
+        parent_length=11,
+    )
+    # providing a value for parent_length overrides standard
+    ov = im1.merge_maps(im2, parent_length=im2.parent_length, aligned_indices=True)
+    assert ov.gap_pos.tolist() == [5]
+    assert ov.cum_gap_lengths.tolist() == [8]
+
+
+def test_indel_map_merge_aligned_indices():
+    im1 = IndelMap(
+        gap_pos=numpy.array([0, 4, 7], dtype=int),
+        gap_lengths=numpy.array([14, 9, 2], dtype=int),
+        parent_length=14,
+    )
+    # in im2, the first and second gaps fall within the first
+    # gap of im1 producing redundant seq indices that have to be
+    # merged and their individual lengths summed
+    im2 = IndelMap(
+        gap_pos=numpy.array([3, 6, 31], dtype=int),
+        gap_lengths=numpy.array([1, 2, 1], dtype=int),
+        parent_length=46,
+    )
+    merged = im1.merge_maps(im2, parent_length=46, aligned_indices=True)
+    exp = IndelMap(
+        gap_pos=numpy.array([0, 4, 7], dtype=int),
+        gap_lengths=numpy.array([17, 9, 3], dtype=int),
+        parent_length=46,
+    )
+
+    numpy.testing.assert_allclose(merged.array, exp.array)