Repository: cogent3/cogent3
Pull request #946 (2021-09-20 00:48:13)

significant refactor to align_to_ref
Strange bugs were difficult to debug because the code was not structured in a way that facilitated troubleshooting. This patch makes changes to simplify debugging and also fixes the recently identified bug.


Patch:
diff --git a/src/cogent3/app/align.py b/src/cogent3/app/align.py
index 3c5718db9..d64e43cb8 100644
--- a/src/cogent3/app/align.py
+++ b/src/cogent3/app/align.py
@@ -1,5 +1,7 @@
 import warnings
 
+from bisect import bisect_left
+
 from cogent3 import make_tree
 from cogent3.align import (
     global_pairwise,
@@ -9,12 +11,7 @@ from cogent3.align import (
 from cogent3.align.progressive import TreeAlign
 from cogent3.app import dist
 from cogent3.core.alignment import Aligned, Alignment
-from cogent3.core.location import (
-    _gap_insertion_data,
-    _gap_pos_to_map,
-    _interconvert_seq_aln_coords,
-    _merged_gaps,
-)
+from cogent3.core.location import gap_coords_to_map
 from cogent3.core.moltype import get_moltype
 from cogent3.evolve.models import get_model
 
@@ -38,40 +35,294 @@ __email__ = "Gavin.Huttley@anu.edu.au"
 __status__ = "Alpha"
 
 
-def _map_ref_gaps_to_seq(all_ref_seq, curr_ref, other):
-    """inject ref gaps missing from curr_ref into other map
+class _GapOffset:
+    """computes sum of gap lengths preceding a position. Acts like a dict
+    for getting the offset for an integer key with the __getitem__ returning
+    the offset.
+    If your coordinate is an alignment position, set invert=True.
+    Examples
+    --------
+    From sequence coordinate to an alignment coordinate
+
+    >>> seq2aln = _GapOffset({1:3, 7:1})
+    >>> seq_pos = 2
+    >>> aln_pos = seq_pos + seq2aln[seq_pos]
+    >>> aln_pos
+    5
+
+    From alignment coordinate to a sequence coordinate
+
+    >>> aln2seq = _GapOffset({1:3, 7:1}, invert=True)
+    >>> seq_pos = aln_pos - aln2seq[seq_pos]
+    >>> seq_pos
+    2
+    """
+
+    def __init__(self, gaps_lengths, invert=False):
+        """
+        Parameters
+        ----------
+        gaps_lengths : dict
+            {pos: length, ...} where pos is a gap insert position and length
+            how long the gap is.
+        invert : bool
+            if True, query keys are taken as being in alignment coordinates
+        """
+        offset = 0
+        min_val = None
+        result = {}
+        k = -1
+        for k, l in sorted(gaps_lengths.items()):
+            if invert:
+                result[k + offset + l] = offset + l
+                result[k + offset] = offset
+            else:
+                result[k] = offset
+
+            offset += l
+            if min_val is None:
+                min_val = k
+
+        self._store = result
+        self.min_pos = min_val
+        self.max_pos = k + offset if invert else k
+        self.total = offset
+        self._ordered = None
+        self._invert = invert
+
+    def __repr__(self):
+        return repr(self._store)
+
+    def __str__(self):
+        return str(self._store)
+
+    def __getitem__(self, k):
+        if not self._store:
+            return 0
+
+        if k in self._store:
+            return self._store[k]
+
+        if k < self.min_pos:
+            return 0
+
+        if k > self.max_pos:
+            return self.total
+
+        if self._ordered is None:
+            self._ordered = sorted(self._store)
+
+        # k is definitely bounded by min and max positions here
+        i = bisect_left(self._ordered, k)
+        pos = self._ordered[i]
+        if self._invert:
+            pos = pos if pos in [k, 0] else self._ordered[i - 1]
+        return self._store[pos]
+
+
+def _gap_union(seqs) -> dict:
+    """returns the union of all gaps in seqs"""
+    seq_name = None
+    all_gaps = {}
+    for seq in seqs:
+        if not isinstance(seq, Aligned):
+            raise TypeError(f"must be Aligned instances, not {type(seq)}")
+        if seq_name is None:
+            seq_name = seq.name
+        if seq.name != seq_name:
+            raise ValueError("all sequences must have the same name")
+
+        gaps_lengths = dict(seq.map.get_gap_coordinates())
+        all_gaps = _merged_gaps(all_gaps, gaps_lengths)
+    return all_gaps
+
+
+def _gap_difference(seq_gaps: dict, union_gaps: dict) -> tuple:
+    """
+
+    Parameters
+    ----------
+    seq_gaps
+        {gap pos: length, } of a sequence used to generate the gap union
+    union_gaps
+        {gap pos: maximum length, ...} derived from the same seq aligned
+        to different sequences
 
     Returns
     -------
-    Map
+    gaps missing from seq_gaps, seq_gaps that overlap with union gaps
     """
-    all_gap_pos, _ = _gap_insertion_data(all_ref_seq)
-    c_gpos, c_offsets = _gap_insertion_data(curr_ref)
-    o_gpos, o_offsets = _gap_insertion_data(other)
-    unique_ref_gaps = set(all_gap_pos) - set(c_gpos)
-    if not unique_ref_gaps:  # other seq map unchanged
-        return other.map
-
-    # we add the offset from the full ref up to the new gap
-    new_gaps = []
-    for unique in unique_ref_gaps:
-        # convert ref seq coord into current alignment coord
-        aln_pos = _interconvert_seq_aln_coords(
-            c_gpos, c_offsets, unique[0], seq_pos=True
-        )
-        # convert current alignment coord into query seq coord
-        seq_pos = _interconvert_seq_aln_coords(
-            o_gpos, o_offsets, aln_pos, seq_pos=False
+    missing = {}
+    overlapping = {}
+    for position, length in union_gaps.items():
+        if position not in seq_gaps:
+            missing[position] = length
+        elif seq_gaps[position] != length:
+            overlapping[position] = length - seq_gaps[position]
+    return missing, overlapping
+
+
+def _merged_gaps(a_gaps: dict, b_gaps: dict) -> dict:
+    """merges gaps that occupy same position
+
+    Parameters
+    ----------
+    a_gaps, b_gaps
+        [(gap position, length),...]
+
+    Returns
+    -------
+    Merged places as {gap position: length, ...}
+
+    Notes
+    -----
+    If a_gaps and b_gaps are from the same underlying sequence, set
+    function to 'max'. Use 'sum' when the gaps derive from different
+    sequences.
+    """
+
+    if not a_gaps:
+        return b_gaps
+
+    if not b_gaps:
+        return a_gaps
+
+    places = set(a_gaps) | set(b_gaps)
+    return {
+        place: max(
+            a_gaps.get(place, 0),
+            b_gaps.get(place, 0),
         )
-        new_gaps.append((seq_pos, unique[1]))
+        for place in places
+    }
+
+
+def _subset_gaps_to_align_coords(
+    subset_gaps: dict, orig_gaps: dict, seq_2_aln: _GapOffset
+) -> dict:
+    """compute alignment coords of subset gaps
 
-    # because we combine gap positions from two different sequences,
-    # the reference and the other, we need to merge overlapping gaps by
-    # summation, rather than max
-    o_gpos = _merged_gaps(o_gpos, new_gaps, function="sum")
-    gap_pos, gap_lengths = list(zip(*o_gpos))
-    seq = other.data
-    return _gap_pos_to_map(gap_pos, gap_lengths, len(seq))
+    Parameters
+    ----------
+    subset_gaps : dict
+        {position: length delta} lengths are the adjusted gap lengths
+    orig_gaps : dict
+        {position: orig length} the original gap lengths are from a pairwise
+        alignment
+    seq_2_aln : dict
+        {seq position: alignment position, ...}
+
+    Returns
+    -------
+    dict
+        {alignment position + orig length: length delta, ...}
+
+    Notes
+    -----
+    """
+    result = {}
+    for p in subset_gaps:
+        offset = seq_2_aln[p]
+        result[offset + p + orig_gaps[p]] = subset_gaps[p]
+
+    return result
+
+
+def _combined_refseq_gaps(seq_gaps: dict, union_gaps: dict) -> dict:
+    # takes union gaps and refseq gaps, converts into diffs and
+    # subset diffs
+    seq2aln = _GapOffset(seq_gaps)
+    diff_gaps, subset_gaps = _gap_difference(seq_gaps, union_gaps)
+    align_coord_gaps = _subset_gaps_to_align_coords(subset_gaps, seq_gaps, seq2aln)
+    align_coord_gaps.update({p + seq2aln[p]: diff_gaps[p] for p in diff_gaps})
+    return align_coord_gaps
+
+
+def _gaps_for_injection(other_seq_gaps: dict, refseq_gaps: dict, seqlen: int) -> dict:
+    """projects refseq aligned gaps into otherseq
+
+    Parameters
+    ----------
+    other_seq_gaps : dict
+        {gap in other seq position: gap length}
+    refseq_gaps : dict
+        {gap as alignment position: gap length}
+    seqlen : int
+        length of sequence being injected into
+
+    Returns
+    -------
+    dict
+        {gap in other seq position: gap length}
+    """
+    aln2seq = _GapOffset(other_seq_gaps, invert=True)
+    # to inject a gap means to convert it from alignment coordinates into
+    # sequence coordinates
+    # we probably need to include the refseq gap union because we need to
+    # establish whether a refseq gap overlaps with a gap in other seq
+    # and
+    all_gaps = {}
+    all_gaps.update(other_seq_gaps)
+    for gap_pos, gap_length in sorted(refseq_gaps.items()):
+        offset = aln2seq[gap_pos]
+        gap_pos -= offset
+        gap_pos = min(seqlen, gap_pos)
+        if gap_pos < 0:
+            raise ValueError(
+                f"computed gap_pos {gap_pos} < 0, correct reference sequence?"
+            )
+        if gap_pos in all_gaps:
+            gap_length += all_gaps[gap_pos]
+
+        all_gaps[gap_pos] = gap_length
+
+    return all_gaps
+
+
+def pairwise_to_multiple(pwise, ref_seq, moltype, info=None):
+    """
+    turns pairwise alignments to a reference into a multiple alignment
+
+    Parameters
+    ----------
+    pwise
+        Series of pairwise alignments to ref_seq as
+        [(non-refseq name, aligned pair), ...]
+    ref_seq
+        The sequence common in all pairwise alignments
+    moltype
+        molecular type for the returned alignment
+    info
+        info object
+
+    Returns
+    -------
+    ArrayAlign
+    """
+    if not hasattr(ref_seq, "name"):
+        raise TypeError(f"ref_seq must be a cogent3 sequence, not {type(ref_seq)}")
+
+    refseqs = [s for _, aln in pwise for s in aln.seqs if s.name == ref_seq.name]
+    ref_gaps = _gap_union(refseqs)
+
+    m = gap_coords_to_map(ref_gaps, len(ref_seq))
+    aligned = [Aligned(m, ref_seq)]
+    for other_name, aln in pwise:
+        curr_ref = aln.named_seqs[ref_seq.name]
+        curr_ref_gaps = dict(curr_ref.map.get_gap_coordinates())
+        other_seq = aln.named_seqs[other_name]
+        other_gaps = dict(other_seq.map.get_gap_coordinates())
+        diff_gaps = _combined_refseq_gaps(curr_ref_gaps, ref_gaps)
+        inject = _gaps_for_injection(other_gaps, diff_gaps, len(other_seq.data))
+        if inject:
+            m = gap_coords_to_map(inject, len(other_seq.data))
+            other_seq = Aligned(m, other_seq.data)
+
+        aligned.append(other_seq)
+    # default to ArrayAlign
+    return Alignment(aligned, moltype=moltype, info=info).to_type(
+        array_align=True, moltype=moltype
+    )
 
 
 class align_to_ref(ComposableSeq):
@@ -154,7 +405,6 @@ class align_to_ref(ComposableSeq):
             seqs = seqs.to_moltype(self._moltype)
 
         ref_seq = seqs.get_seq(self._ref_name)
-        ref_gaps = []
         kwargs = self._kwargs.copy()
         pwise = []
         for seq in seqs.seqs:
@@ -163,24 +413,8 @@ class align_to_ref(ComposableSeq):
 
             aln = global_pairwise(ref_seq, seq, **kwargs).to_type(array_align=False)
             pwise.append(((seq.name, aln)))
-            gap_pos, _ = _gap_insertion_data(aln.named_seqs[ref_seq.name])
-            ref_gaps = _merged_gaps(ref_gaps, gap_pos, function="max")
-
-        gap_pos, gap_lengths = list(zip(*ref_gaps)) if ref_gaps else ([], [])
-        m = _gap_pos_to_map(gap_pos, gap_lengths, len(ref_seq))
-        aligned = [Aligned(m, ref_seq)]
-        for other_name, aln in pwise:
-            curr_ref = aln.named_seqs[ref_seq.name]
-            other_seq = aln.named_seqs[other_name]
-            m = _map_ref_gaps_to_seq(aligned[0], curr_ref, other_seq)
-            aligned.append(
-                other_seq if m is other_seq.map else Aligned(m, other_seq.data)
-            )
 
-        # default to ArrayAlign
-        return Alignment(aligned, moltype=self._moltype, info=seqs.info).to_type(
-            array_align=True, moltype=self._moltype
-        )
+        return pairwise_to_multiple(pwise, ref_seq, self._moltype, info=seqs.info)
 
 
 class progressive_align(ComposableSeq):
diff --git a/src/cogent3/core/location.py b/src/cogent3/core/location.py
index 7d2c0025f..00c7540fb 100644
--- a/src/cogent3/core/location.py
+++ b/src/cogent3/core/location.py
@@ -713,6 +713,18 @@ class Map(object):
         spans = [s.reversed_relative_to(self.parent_length) for s in self.spans]
         return Map(spans=spans, parent_length=self.parent_length)
 
+    def get_gap_coordinates(self):
+        """returns [(gap pos, gap length), ...]"""
+        gap_pos = []
+        for i, span in enumerate(self.spans):
+            if not span.lost:
+                continue
+
+            pos = self.spans[i - 1].end if i else 0
+            gap_pos.append((pos, len(span)))
+
+        return gap_pos
+
     def gaps(self):
         """The gaps (lost spans) in this map"""
         locations = []
@@ -1048,98 +1060,32 @@ def RangeFromString(string, delimiter=","):
     return result
 
 
-def _gap_insertion_data(seq):
-    """compute gap position, length and cumulative offsets
-
-    Parameters
-    ----------
-    seq
-        a cogent3 annotatable sequence
-
-    Returns
-    -------
-    [(seq position, gap length), ...], [cum sum gap length, ...]
-
-    Notes
-    -----
-    The sequence position is in unaligned sequence coordinates. offsets are
-    calculated as the cumulative sum of gap lengths. The offset
-    plus the sequence position gives the alignment coordinate for a gap.
+def gap_coords_to_map(gaps_lengths: dict, seq_length: int) -> Map:
     """
-    gap_pos = []
-    offsets = []
-    offset = 0
-    for i, span in enumerate(seq.map.spans):
-        if not span.lost:
-            continue
-        pos = seq.map.spans[i - 1].end if i else 0
-        gap_pos.append((pos, len(span)))
-        offsets.append(offset)
-        offset += span.length
-
-    return gap_pos, offsets
-
-
-def _merged_gaps(a_gaps, b_gaps, function="max"):
-    """merges gaps that occupy same position
-
     Parameters
     ----------
-    a_gaps, b_gaps
-        [(gap position, length),...]
-    function
-        When a and b contain a gap at the same position, function is applied
-        to the gap lengths. Valid values are either 'max' or 'sum'.
+    gaps_lengths
+        {gap insertion pos: gap length, ...}
+    seq_length : int
+        length of unaligned sequence
 
     Returns
     -------
-    Merged places as [(gap position, length),...]
-
-    Notes
-    -----
-    If a_gaps and b_gaps are from the same underlying sequence, set
-    function to 'max'. Use 'sum' when the gaps derive from different
-    sequences.
+    Map
     """
-    if function.lower() not in "maxsum":
-        raise ValueError(f"{function} not allowed, choose either 'sum' or 'max'")
 
-    function = max if function.lower() == "max" else sum
-
-    if not a_gaps:
-        return b_gaps
-
-    if not b_gaps:
-        return a_gaps
-
-    places = sorted(a_gaps + b_gaps)
-
-    positions, lengths = [places[0][0]], [places[0][1]]
-    for i, (pos, l) in enumerate(places[1:], 1):
-        if positions[-1] == pos:
-            lengths[-1] = function([lengths[-1], l])
-            continue
-
-        positions.append(pos)
-        lengths.append(l)
-    return list(zip(positions, lengths))
-
-
-def _gap_pos_to_map(gap_pos, gap_lengths, seq_length):
-    """[(pos, gap length), ...]"""
-
-    if not gap_pos:
+    if not gaps_lengths:
         return Map([(0, seq_length)], parent_length=seq_length)
 
     spans = []
     last = pos = 0
-    for i, pos in enumerate(gap_pos):
+    for pos in sorted(gaps_lengths):
         if pos > seq_length:
             raise ValueError(
                 f"cannot have gap at position {pos} beyond seq_length= {seq_length}"
             )
 
-        gap = LostSpan(length=gap_lengths[i])
+        gap = LostSpan(length=gaps_lengths[pos])
         spans.extend([gap] if pos == 0 else [Span(last, pos), gap])
         last = pos
 
@@ -1147,40 +1093,3 @@ def _gap_pos_to_map(gap_pos, gap_lengths, seq_length):
         spans.append(Span(last, seq_length))
 
     return Map(spans=spans, parent_length=seq_length)
-
-
-def _interconvert_seq_aln_coords(gaps, offsets, pos, seq_pos=True):
-    """converts sequence position to an alignment position
-
-    Parameters
-    ----------
-    gaps
-        series of [(seq pos, length), ..]
-    offsets
-        the offset of seq pos, which is basically the sum of all lengths up
-        to the previous gap
-    p
-        the sequence coordinate to convert
-    seq_pos : bool
-        whether pos is in sequence coordinates. If False, means it is in
-        alignment coordinates.
-    Returns
-    -------
-    alignment coordinate
-    """
-    if pos < 0:
-        raise ValueError(f"negative value {pos}")
-
-    if not gaps or pos == 0:
-        return pos
-
-    offset = 0
-    for i, (p, len) in enumerate(gaps):
-        assert p >= 0 and len > 0
-        if p + offset >= pos:
-            break
-        offset += len
-    assert offset < pos, f"calculated offset {offset} greater than align pos {p}"
-    if not seq_pos:  # we subtract the gap length total to get to seq coords
-        offset = -offset
-    return pos + offset


Test patch:
diff --git a/tests/test_app/test_align.py b/tests/test_app/test_align.py
index 50bd897ad..a8f8da9cb 100644
--- a/tests/test_app/test_align.py
+++ b/tests/test_app/test_align.py
@@ -9,18 +9,18 @@ from cogent3 import (
 )
 from cogent3.align.align import make_generic_scoring_dict
 from cogent3.app import align as align_app
-from cogent3.app.align import _map_ref_gaps_to_seq
-from cogent3.app.composable import NotCompleted
-from cogent3.core.alignment import Aligned, Alignment
-from cogent3.core.location import (
-    LostSpan,
-    Map,
-    Span,
-    _gap_insertion_data,
-    _gap_pos_to_map,
-    _interconvert_seq_aln_coords,
+from cogent3.app.align import (
+    _combined_refseq_gaps,
+    _gap_difference,
+    _gap_union,
+    _GapOffset,
+    _gaps_for_injection,
     _merged_gaps,
+    pairwise_to_multiple,
 )
+from cogent3.app.composable import NotCompleted
+from cogent3.core.alignment import Aligned
+from cogent3.core.location import gap_coords_to_map
 
 
 __author__ = "Gavin Huttley"
@@ -65,6 +65,27 @@ _codon_models = [
 ]
 
 
+def make_pairwise(data, refseq_name, moltype="dna", array_align=False):
+    """returns series of refseq, [(n, pwise aln),..]. All alignments are to ref_seq"""
+    aln = make_aligned_seqs(
+        data,
+        array_align=array_align,
+        moltype=moltype,
+    )
+    refseq = aln.get_seq(refseq_name)
+    pwise = [
+        (n, aln.take_seqs([refseq_name, n]).omit_gap_pos())
+        for n in aln.names
+        if n != refseq_name
+    ]
+    return refseq, pwise
+
+
+def make_aligned(gaps_lengths, seq, name="seq1"):
+    seq = seq.moltype.make_seq(seq, name=name)
+    return Aligned(gap_coords_to_map(gaps_lengths, len(seq)), seq)
+
+
 class RefalignmentTests(TestCase):
     seqs = make_unaligned_seqs(_seqs, moltype=DNA)
 
@@ -97,124 +118,16 @@ class RefalignmentTests(TestCase):
         got = aligner(self.seqs)
         self.assertEqual(got.moltype.label, "dna")
 
-    def test__gap_insertion_data(self):
-        """identifies gap locations and lengths"""
-        seq = DNA.make_seq("AACCCCCCGGG")
-        # no gaps
-        m = _gap_pos_to_map([], [], len(seq))
-        got = _gap_insertion_data(Aligned(m, seq))
-        expect = [], []
-        self.assertEqual(got, expect)
-        # one gap
-        gap_positions = [2]
-        gap_lengths = [3]
-        m = _gap_pos_to_map(gap_positions, gap_lengths, len(seq))
-        got = _gap_insertion_data(Aligned(m, seq))
-        expect = list(zip(gap_positions, gap_lengths)), [0]
-        self.assertEqual(got, expect)
-        # two gaps
-        gap_positions = [2, 8]
-        gap_lengths = [3, 1]
-        m = _gap_pos_to_map(gap_positions, gap_lengths, len(seq))
-        got = _gap_insertion_data(Aligned(m, seq))
-        expect = list(zip(gap_positions, gap_lengths)), [0, 3]
-        self.assertEqual(got, expect)
-
-    def test__merged_gaps(self):
+    def test_merged_gaps(self):
         """correctly merges gaps"""
-        with self.assertRaises(ValueError):
-            _merged_gaps([], [], function="blah")
-
-        a = [(2, 3), (4, 9)]
-        b = [(2, 6), (8, 5)]
+        a = dict([(2, 3), (4, 9)])
+        b = dict([(2, 6), (8, 5)])
         # omitting one just returns the other
-        self.assertIs(_merged_gaps(a, []), a)
-        self.assertIs(_merged_gaps([], b), b)
-        # specifying 'sum' means (2, 9)
-        got = _merged_gaps(a, b, function="sum")
-        self.assertEqual(got, [(2, 9), (4, 9), (8, 5)])
-        # specifying 'max' means (2, 6)
-        got = _merged_gaps(a, b, function="max")
+        self.assertIs(_merged_gaps(a, {}), a)
+        self.assertIs(_merged_gaps({}, b), b)
+        got = _merged_gaps(a, b)
         self.assertEqual(got, [(2, 6), (4, 9), (8, 5)])
 
-    def test__gap_pos_to_map(self):
-        """correctly converts a gap positions to a Map"""
-        gap_positions = [2, 9]
-        gap_lengths = [3, 1]
-        map = _gap_pos_to_map(gap_positions, gap_lengths, 20)
-        self.assertTrue(all(map.spans[i].lost for i in (1, 3)))
-        self.assertEqual(len(map), 24)
-        # no gaps
-        map = _gap_pos_to_map([], [], 20)
-        self.assertEqual(len(map), 20)
-        self.assertEqual(len(map.spans), 1)
-        # gap at start
-        gap_positions = [0]
-        gap_lengths = [3]
-        map = _gap_pos_to_map(gap_positions, gap_lengths, 20)
-        self.assertTrue(map.spans[0].lost and not map.spans[1].lost)
-        self.assertEqual(len(map), 23)
-        self.assertEqual(len(map.spans), 2)
-        # gap at end
-        gap_positions = [20]
-        gap_lengths = [3]
-        map = _gap_pos_to_map(gap_positions, gap_lengths, 20)
-        self.assertTrue(map.spans[-1].lost and not map.spans[0].lost)
-        self.assertEqual(len(map), 23)
-        self.assertEqual(len(map.spans), 2)
-        # fail if pos beyond sequence
-        with self.assertRaises(ValueError):
-            _gap_pos_to_map([40], [2], 20)
-
-    def test__map_ref_gaps_to_seq(self):
-        """correctly handle case where ref and curr ref are equal"""
-        aln = make_aligned_seqs(
-            {
-                "Ref": "CAG---GAGAACAGAAACCCAT--TACTCACT",
-                "Qu2": "CAGCATGAGAACAGAAACCCGT--TA---ACT",
-            },
-            array_align=False,
-            moltype="dna",
-        )
-        ref = aln.named_seqs["Ref"]
-        other = aln.named_seqs["Qu2"]
-        got = _map_ref_gaps_to_seq(ref, ref, other)
-        self.assertIs(got, other.map)
-
-    def test_seq_2_aln_coords(self):
-        """correctly converts sequence coordinates to alignment coordinates"""
-        got = _interconvert_seq_aln_coords([], [], 20, seq_pos=True)
-        assert got == 20
-
-        # aligned seq has one gap after
-        got = _interconvert_seq_aln_coords([(22, 3)], [0], 20, seq_pos=True)
-        assert got == 20
-
-        # aligned seq has one 2bp gap before
-        got = _interconvert_seq_aln_coords([(12, 2)], [0], 20, seq_pos=True)
-        assert got == 22
-
-        # aligned seq has two 2bp gap before
-        got = _interconvert_seq_aln_coords([(10, 2), (12, 2)], [0], 20, seq_pos=True)
-        assert got == 24
-
-    def test_aln_2_seq_coord(self):
-        """correctly converts alignment coordinates to sequence coordinates"""
-        got = _interconvert_seq_aln_coords([], [], 20, seq_pos=False)
-        assert got == 20
-
-        # query has one gap after
-        got = _interconvert_seq_aln_coords([(22, 3)], [0], 20, seq_pos=False)
-        assert got == 20
-
-        # query has one 2bp gap before
-        got = _interconvert_seq_aln_coords([(12, 2)], [0], 20, seq_pos=False)
-        assert got == 18
-
-        # query has two 2bp gap before
-        got = _interconvert_seq_aln_coords([(10, 2), (12, 2)], [0], 20, seq_pos=False)
-        assert got == 16
-
     def test_aln_to_ref_known(self):
         """correctly recapitulates known case"""
         orig = make_aligned_seqs(
@@ -235,6 +148,175 @@ class RefalignmentTests(TestCase):
         aln = aligner(orig.degap())
         self.assertEqual(aln.to_dict(), expect)
 
+    def test_gap_union(self):
+        """correctly identifies the union of all gaps"""
+        # fails if not all sequences same
+        seq = DNA.make_seq("AACCCGTT")
+        all_gaps = dict([(0, 3), (2, 1), (5, 3), (6, 3)])
+        final_seq = make_aligned(all_gaps, seq)
+        gap_sets = [
+            dict([(5, 1), (6, 3)]),
+            dict([(2, 1), (5, 3)]),
+            dict([(2, 1), (5, 1), (6, 2)]),
+            dict([(0, 3)]),
+        ]
+        seqs = [make_aligned(gaps, seq) for gaps in gap_sets]
+        got = _gap_union(seqs)
+        self.assertEqual(got, dict(all_gaps))
+
+        # must all be Aligned instances
+        with self.assertRaises(TypeError):
+            _gap_union(seqs + ["GGGGGGGG"])
+
+        # must all have the same name
+        with self.assertRaises(ValueError):
+            _gap_union(seqs + [make_aligned({}, seq, name="blah")])
+
+    def test_gap_difference(self):
+        """correctly identifies the difference in gaps"""
+        seq = DNA.make_seq("AACCCGTT")
+        all_gaps = dict([(0, 3), (2, 1), (5, 3), (6, 3)])
+        gap_sets = [
+            dict([(5, 1), (6, 3)]),
+            dict([(2, 1), (5, 3)]),
+            dict([(2, 1), (5, 1), (6, 2)]),
+            dict([(0, 3)]),
+        ]
+        seqs = [make_aligned(gaps, seq) for gaps in gap_sets]
+        union = _gap_union(seqs)
+        expects = [
+            [dict([(0, 3), (2, 1)]), dict([(5, 2)])],
+            [dict([(0, 3), (6, 3)]), {}],
+            [dict([(0, 3)]), dict([(5, 2), (6, 1)])],
+            [dict([(2, 1), (5, 3), (6, 3)]), {}],
+        ]
+        for seq, (plain, overlap) in zip(seqs, expects):
+            seq_gaps = dict(seq.map.get_gap_coordinates())
+            got_plain, got_overlap = _gap_difference(seq_gaps, union)
+            self.assertEqual(got_plain, dict(plain))
+            self.assertEqual(got_overlap, dict(overlap))
+
+    def test_merged_gaps(self):
+        """correctly handles gap values"""
+        a_gaps = {0: 2}
+        b_gaps = {2: 2}
+        self.assertEqual(_merged_gaps(a_gaps, {}), a_gaps)
+        self.assertEqual(_merged_gaps({}, b_gaps), b_gaps)
+
+    def test_combined_refseq_gaps(self):
+        union = dict([(0, 3), (2, 1), (5, 3), (6, 3)])
+        gap_sets = [
+            [(5, 1), (6, 3)],
+            [(2, 1), (5, 3)],
+            [(2, 1), (5, 1), (6, 2)],
+            [(0, 3)],
+        ]
+        # for subset gaps, their alignment position is the
+        # offset + their position + their gap length
+        expects = [
+            dict([(6, 2), (0, 3), (2, 1)]),
+            dict([(0, 3), (10, 3)]),
+            dict([(0, 3), (5 + 1 + 1, 2), (6 + 2 + 2, 1)]),
+            dict([(2 + 3, 1), (5 + 3, 3), (6 + 3, 3)]),
+        ]
+        for i, gap_set in enumerate(gap_sets):
+            got = _combined_refseq_gaps(dict(gap_set), union)
+            self.assertEqual(got, expects[i])
+
+        # if union gaps equals ref gaps
+        got = _combined_refseq_gaps({2: 2}, {2: 2})
+        self.assertEqual(got, {})
+
+    def test_gaps_for_injection(self):
+        # for gaps before any otherseq gaps, alignment coord is otherseq coord
+        oseq_gaps = {2: 1, 6: 2}
+        rseq_gaps = {0: 3}
+        expect = {0: 3, 2: 1, 6: 2}
+        seqlen = 50
+        got = _gaps_for_injection(oseq_gaps, rseq_gaps, seqlen)
+        self.assertEqual(got, expect)
+        # for gaps after otherseq gaps seq coord is align coord minus gap
+        # length totals
+        got = _gaps_for_injection(oseq_gaps, {4: 3}, seqlen)
+        expect = {2: 1, 3: 3, 6: 2}
+        self.assertEqual(got, expect)
+        got = _gaps_for_injection(oseq_gaps, {11: 3}, seqlen)
+        expect = {2: 1, 6: 2, 8: 3}
+        self.assertEqual(got, expect)
+        # gaps beyond sequence length added to end of sequence
+        got = _gaps_for_injection({2: 1, 6: 2}, {11: 3, 8: 3}, 7)
+        expect = {2: 1, 6: 2, 7: 6}
+        self.assertEqual(got, expect)
+
+    def test_pairwise_to_multiple(self):
+        """the standalone function constructs a multiple alignment"""
+        expect = {
+            "Ref": "CAG---GAGAACAGAAACCCAT--TACTCACT",
+            "Qu1": "CAG---GAGAACAG---CCCGTGTTACTCACT",
+            "Qu2": "CAGCATGAGAACAGAAACCCGT--TA---ACT",
+            "Qu3": "CAGCATGAGAACAGAAACCCGT----CTCACT",
+            "Qu7": "CAG---GA--ACAGA--CCCGT--TA---ACT",
+            "Qu4": "CAGCATGAGAACAGAAACCCGTGTTACTCACT",
+            "Qu5": "CAG---GAGAACAG---CCCAT--TACTCACT",
+            "Qu6": "CAG---GA-AACAG---CCCAT--TACTCACT",
+        }
+        aln = make_aligned_seqs(expect, moltype="dna").omit_gap_pos()
+        expect = aln.to_dict()
+        for refseq_name in ["Qu3"]:
+            refseq, pwise = make_pairwise(expect, refseq_name)
+            got = pairwise_to_multiple(pwise, ref_seq=refseq, moltype=refseq.moltype)
+            self.assertEqual(len(got), len(aln))
+            orig = dict(pwise)
+            _, pwise = make_pairwise(got.to_dict(), refseq_name)
+            got = dict(pwise)
+            # should be able to recover the original pairwise alignments
+            for key, value in got.items():
+                self.assertEqual(value.to_dict(), orig[key].to_dict(), msg=refseq_name)
+
+            with self.assertRaises(TypeError):
+                pairwise_to_multiple(pwise, "ACGG", DNA)
+
+    def test_pairwise_to_multiple_2(self):
+        """correctly handle alignments with gaps beyond end of query"""
+        # cogent3.core.alignment.DataError: Not all sequences are the same length:
+        # max is 425, min is 419
+        def make_pwise(data, ref_name):
+            result = []
+            for n, seqs in data.items():
+                result.append(
+                    [n, make_aligned_seqs(data=seqs, moltype="dna", array_align=False)]
+                )
+            ref_seq = result[0][1].get_seq(ref_name)
+            return result, ref_seq
+
+        pwise = {
+            "Platypus": {
+                "Opossum": "-----------------GTGC------GAT-------------------------------CCAAAAACCTGTGTC--ACCGT--------GCC----CAGAGCCTCC----CTCAGGCCGCTCGGGGAG---TG-------GCCCCCCG--GC-GGAGGGCAGGGATGGGGAGT-AGGGGTGGCAGTC----GGAACTGGAAGAGCTT-TACAAACC---------GA--------------------GGCT-AGAGGGTC-TGCTTAC-------TTTTTACCTTGG------------GTTTG-CCAGGAGGTAG----------AGGATGA-----------------CTAC--ATCAAG----AGC------------TGGG-------------",
+                "Platypus": "CAGGATGACTACATCAAGAGCTGGGAAGATAACCAGCAAGGAGATGAAGCTCTGGACACTACCAAAGACCCCTGCCAGAACGTGAAGTGCAGCCGACACAAGGTCTGCATCGCTCAGGGCTACCAGAGAGCCATGTGTATCAGCCGCAAGAAGCTGGAGCACAGGATCAAGCAGCCAGCCCTGAAACTCCATGGAAACAGAGAGAGCTTCTGCAAGCCTTGTCACATGACCCAGCTGGCCTCTGTCTGCGGCTCGGACGGACACACTTACAGCTCCGTGTGCAAACTGGAGCAGCAGGCCTGTCTGACCAGCAAGCAGCTGACAGTCAAGTGTGAAGGCCAGTGCCCGTGCCCCACCGATCATGTTCCAGCCTCCACCGCTGATGGAAAACAAGAGACCT",
+            },
+            "Wombat": {
+                "Opossum": "GTGCGATCCAAAAACCTGTGTCACCGTGCCCAGAGCCTCCCTCAGGCCGCTCGG-GGAGTGGCCCCCCGGCGGAGGGCAGGGATGGGGAGTAGGGGTGGCAGTCGGAACTGGAAGAGCTTTACAAACCGAGGCTAGAGGGTCTGCTTACTTTTTACCTTGG------GTTT--GC-CAGGA---GGT----AGAGGATGACTACATCAAGAGCTGGG---------------------------",
+                "Wombat": "--------CA----------TCACCGC-CCCTGCACC---------CGGCTCGGCGGAGGGGGATTCTAA-GGGGGTCAAGGATGGCGAG-ACCCCTGGCAATTTCA--TGGAGGA------CGAGCAATGGCT-----GTC-GTCCATCTCCCAGTATAGCGGCAAGATCAAGCACTGGAACCGCTTCCGAGACGATGACTACATCAAGAGCTGGGAGGACAGTCAGCAAGGAGATGAAGCGC",
+            },
+        }
+        pwise, ref_seq = make_pwise(pwise, "Opossum")
+        aln = pairwise_to_multiple(pwise, ref_seq, ref_seq.moltype)
+        self.assertNotIsInstance(aln, NotCompleted)
+
+        pwise = {
+            "Platypus": {
+                "Opossum": "-----------------GTGC------GAT-------------------------------CCAAAAACCTGTGTC",
+                "Platypus": "CAGGATGACTACATCAAGAGCTGGGAAGATAACCAGCAAGGAGATGAAGCTCTGGACACTACCAAAGACCCCTGCC",
+            },
+            "Wombat": {
+                "Opossum": "GTGCGATCCAAAAACCTGTGTC",
+                "Wombat": "--------CA----------TC",
+            },
+        }
+        pwise, ref_seq = make_pwise(pwise, "Opossum")
+        aln = pairwise_to_multiple(pwise, ref_seq, ref_seq.moltype)
+        self.assertNotIsInstance(aln, NotCompleted)
+
 
 class ProgressiveAlignment(TestCase):
     seqs = make_unaligned_seqs(_seqs, moltype=DNA)
@@ -363,5 +445,45 @@ class ProgressiveAlignment(TestCase):
         self.assertEqual(len(aln), 14)
 
 
+class GapOffsetTests(TestCase):
+    def test_empty(self):
+        """create an empty offset"""
+        goff = _GapOffset({})
+        for i in range(4):
+            self.assertEqual(goff[i], 0)
+
+        goff = _GapOffset({}, invert=True)
+        for i in range(4):
+            self.assertEqual(goff[i], 0)
+
+    def test_repr_str(self):
+        """repr and str work"""
+        goff = _GapOffset({}, invert=True)
+        for func in (str, repr):
+            self.assertEqual(func(goff), "{}")
+
+    def test_gap_offset(self):
+        goff = _GapOffset({1: 2, 3: 4})
+        self.assertEqual(goff.min_pos, 1)
+        self.assertEqual(goff.max_pos, 3)
+        self.assertEqual(goff.total, 6)
+        self.assertEqual(goff[0], 0)
+        self.assertEqual(goff[1], 0)
+        self.assertEqual(goff[2], 2)
+        self.assertEqual(goff[3], 2)
+        self.assertEqual(goff[4], 6)
+
+    def test_gap_offset_invert(self):
+        aln2seq = _GapOffset({2: 1, 5: 2, 7: 2}, invert=True)
+        self.assertEqual(aln2seq._store, {3: 1, 2: 0, 8: 3, 6: 1, 12: 5, 10: 3})
+        self.assertEqual(aln2seq.max_pos, 12)
+        self.assertEqual(aln2seq.min_pos, 2)
+        self.assertEqual(aln2seq[11], 3)
+        seq2aln = _GapOffset({2: 1, 5: 2, 7: 2})
+        for seq_pos in range(20):
+            aln_pos = seq_pos + seq2aln[seq_pos]
+            self.assertEqual(aln_pos - aln2seq[aln_pos], seq_pos)
+
+
 if __name__ == "__main__":
     main()
diff --git a/tests/test_core/test_location.py b/tests/test_core/test_location.py
index 590afd7fb..a21219949 100644
--- a/tests/test_core/test_location.py
+++ b/tests/test_core/test_location.py
@@ -4,7 +4,14 @@
 """
 from unittest import TestCase, main
 
-from cogent3.core.location import Map, Range, RangeFromString, Span
+from cogent3 import DNA
+from cogent3.core.location import (
+    Map,
+    Range,
+    RangeFromString,
+    Span,
+    gap_coords_to_map,
+)
 
 
 __author__ = "Rob Knight"
@@ -528,6 +535,38 @@ class MapTests(TestCase):
         coords = map.get_coordinates()
         self.assertEqual(coords, spans)
 
+    def test_get_gap_coords(self):
+        """returns gap start and lengths"""
+        m, seq = DNA.make_seq("-AC--GT-TTA--").parse_out_gaps()
+        got = m.get_gap_coordinates()
+        self.assertEqual(dict(got), {0: 1, 2: 2, 4: 1, 7: 2})
+
+    def test_gap_coords_to_map(self):
+        """construct a Map from coordinates of gap alone"""
+        m, seq = DNA.make_seq("-AC--GT-TTA--").parse_out_gaps()
+        gap_coords = {0: 1, 2: 2, 4: 1, 7: 2}
+        seqlen = 70
+        got = gap_coords_to_map(gap_coords, seqlen)
+        self.assertEqual(len(got), seqlen + sum(gap_coords.values()))
+
+        gap_coords = {5: 2, 17: 3, 10: 2}
+        seqlen = 20
+        got = gap_coords_to_map(gap_coords, seqlen)
+        self.assertEqual(len(got), sum(gap_coords.values()) + seqlen)
+
+        # roundtrip from Map.get_gap_coordinates()
+        self.assertEqual(dict(got.get_gap_coordinates()), gap_coords)
+
+        # and no gaps
+        m, seq = DNA.make_seq("ACGTTTA").parse_out_gaps()
+        got = gap_coords_to_map({}, len(seq))
+        self.assertEqual(len(got), len(m))
+        self.assertEqual(got.get_coordinates(), m.get_coordinates())
+
+        # and gaps outside sequence
+        with self.assertRaises(ValueError):
+            got = gap_coords_to_map({20: 1}, len(seq))
+
 
 # run the following if invoked from command-line
 if __name__ == "__main__":