Repository: cogent3/cogent3
Pull request #2108 (2024-11-13 01:55:24)

More performance optimisations
## Summary by Sourcery

Enhancements:
- Optimize performance by applying Numba JIT compilation to several functions, including _gap_spans, seq_to_align_index, align_to_seq_index, decompose_gapped_seq, and compose_gapped_seq.


Patch:
diff --git a/src/cogent3/core/location.py b/src/cogent3/core/location.py
index b27cb761a..d8c1b4f23 100644
--- a/src/cogent3/core/location.py
+++ b/src/cogent3/core/location.py
@@ -8,6 +8,7 @@ from bisect import bisect_left, bisect_right
 from functools import total_ordering
 from typing import Any, Iterator, Optional, Sequence, Union
 
+import numba
 import numpy
 from numpy.typing import NDArray
 
@@ -962,13 +963,14 @@ def spans_to_gap_coords(
     return numpy.array(gap_pos, dtype=dtype), numpy.array(cum_lengths, dtype=dtype)
 
 
+@numba.jit
 def _gap_spans(
     gap_pos: IntArrayTypes, cum_gap_lengths: IntArrayTypes
-) -> tuple[IntArrayTypes, IntArrayTypes]:
+) -> tuple[IntArrayTypes, IntArrayTypes]:  # pragma: no cover
     """returns 1D arrays in alignment coordinates of
     gap start, gap stop"""
     if not len(gap_pos):
-        r = numpy.array([], dtype=_DEFAULT_GAP_DTYPE)
+        r = numpy.empty((0,), dtype=gap_pos.dtype)
         return r, r
 
     ends = gap_pos + cum_gap_lengths
@@ -1077,6 +1079,125 @@ def _input_vals_neg_step(
     return (0, 0, 1) if start < stop else (start, stop, step)
 
 
+@numba.jit
+def seq_to_align_index(
+    gap_pos: IntArrayTypes,
+    cum_lengths: IntArrayTypes,
+    parent_length: int,
+    num_gaps: int,
+    seq_index: int,
+    slice_stop: bool = False,
+) -> int:  # pragma: no cover
+    """convert a sequence index into an alignment index
+
+    Parameters
+    ----------
+    gap_pos
+        array of gap positions in sequence coordinates
+    cum_lengths
+        array of cumulative gap lengths
+    parent_length
+        length of parent sequence (i.e. aligned sequence without gaps)
+    num_gaps
+        the number of gaps
+    seq_index
+        coordinate on the sequence, must be < parent_length
+    slice_stop
+        set to True if the index is to be the end of an alignment slice.
+        In that case, and if seq_index is in gap_pos then it returns
+        the first alignment index of the gap run.
+    """
+    # NOTE I explicitly cast all returned values to python int's due to
+    # need for json serialisation, which does not support numpy int classes
+    assert len(gap_pos) == len(cum_lengths) == num_gaps
+    if seq_index < 0:
+        seq_index += parent_length
+
+    if seq_index < 0:
+        raise IndexError(f"{seq_index} negative seq_index beyond limit ")
+
+    if not num_gaps or seq_index < gap_pos[0]:
+        return int(seq_index)
+
+    # if slice_stop, check if the seq_index corresponds to a gap position
+    match = seq_index == gap_pos
+    if slice_stop and match.any():
+        # if so, we return the alignment coord for the first gap position
+        (idx,) = numpy.where(match)[0]
+        if idx:
+            gap_len = cum_lengths[idx] - cum_lengths[idx - 1]
+        else:
+            gap_len = cum_lengths[idx]
+        gap_end = gap_pos[idx] + cum_lengths[idx]
+        return int(gap_end - gap_len)
+
+    if seq_index >= gap_pos[-1]:
+        return int(seq_index + cum_lengths[-1])
+
+    # find gap position before seq_index
+    index = numpy.searchsorted(gap_pos, seq_index, side="left")
+    if seq_index < gap_pos[index]:
+        gap_lengths = cum_lengths[index - 1] if index else 0
+    else:
+        gap_lengths = cum_lengths[index]
+
+    return int(seq_index + gap_lengths)
+
+
+@numba.jit
+def align_to_seq_index(
+    gap_pos: IntArrayTypes,
+    cum_lengths: IntArrayTypes,
+    len_aligned: int,
+    num_gaps: int,
+    align_index: int,
+) -> int:  # pragma: no cover
+    """converts alignment index to sequence index
+
+    Parameters
+    ----------
+    gap_pos
+        array of gap positions in sequence coordinates
+    cum_lengths
+        array of cumulative gap lengths
+    len_aligned
+        length of aligned sequence
+    num_gaps
+        the number of gaps
+    align_index
+        coordinate on the alignment, must be < len_aligned
+    """
+    assert len(gap_pos) == len(cum_lengths) == num_gaps
+    # NOTE I explicitly cast all returned values to python int's due to
+    # need for json serialisation, which does not support numpy int classes
+    if align_index < 0:
+        align_index = len_aligned + align_index
+    if align_index < 0:
+        raise IndexError(f"{align_index} align_index beyond limit")
+
+    if not num_gaps or align_index < gap_pos[0]:
+        return align_index
+
+    # these are alignment indices for gaps
+    gap_starts, gap_ends = _gap_spans(gap_pos, cum_lengths)
+    if align_index >= gap_ends[-1]:
+        return int(align_index - cum_lengths[-1])
+
+    index = numpy.searchsorted(gap_ends, align_index, side="left")
+    if align_index < gap_starts[index]:
+        # before the gap at index
+        return int(align_index - cum_lengths[index - 1])
+
+    if align_index == gap_ends[index]:
+        # after the gap at index
+        return int(align_index - cum_lengths[index])
+
+    if gap_starts[index] <= align_index < gap_ends[index]:
+        # within the gap at index
+        # so the gap insertion position is the sequence position
+        return int(gap_pos[index])
+
+
 @dataclasses.dataclass
 class IndelMap(MapABC):
     """store locations of deletions in a Aligned sequence
@@ -1093,7 +1214,7 @@ class IndelMap(MapABC):
         if ``True``, returns new instance with terminal gaps indicated as
         unknown character '?'
     parent_length
-        length of parent sequence (i.e. aligned sequence with gaps)
+        length of parent sequence (i.e. aligned sequence without gaps)
     """
 
     # gap data is gap positions, gap lengths on input, stored
@@ -1122,6 +1243,9 @@ class IndelMap(MapABC):
                 f"gap position {self.gap_pos[-1]} outside parent_length {self.parent_length}"
             )
 
+        # force all to int32
+        self.gap_pos = self.gap_pos.astype(numpy.int32)
+        self.cum_gap_lengths = self.cum_gap_lengths.astype(numpy.int32)
         # make gap array immutable
         self.gap_pos.flags.writeable = False
         self.cum_gap_lengths.flags.writeable = False
@@ -1387,73 +1511,30 @@ class IndelMap(MapABC):
             In that case, and if seq_index is in gap_pos then it returns
             the first alignment index of the gap run.
         """
-        cum_lengths = self.cum_gap_lengths
-        gap_pos = self.gap_pos
+        index = seq_to_align_index(
+            self.gap_pos,
+            self.cum_gap_lengths,
+            self.parent_length,
+            self.num_gaps,
+            seq_index,
+            slice_stop,
+        )
         # NOTE I explicitly cast all returned values to python int's due to
         # need for json serialisation, which does not support numpy int classes
-        if seq_index < 0:
-            seq_index += self.parent_length
-
-        if seq_index < 0:
-            raise IndexError(f"{seq_index} negative seq_index beyond limit ")
-
-        if not self.num_gaps or seq_index < gap_pos[0]:
-            return int(seq_index)
-
-        # if stop_index, check if the seq_index corresponds to a gap position
-        if slice_stop and (match := seq_index == gap_pos).any():
-            # if so, we return the alignment coord for the first gap position
-            (idx,) = numpy.where(match)[0]
-            if idx:
-                gap_len = cum_lengths[idx] - cum_lengths[idx - 1]
-            else:
-                gap_len = cum_lengths[idx]
-            gap_end = gap_pos[idx] + cum_lengths[idx]
-            return int(gap_end - gap_len)
-
-        if seq_index >= gap_pos[-1]:
-            return int(seq_index + cum_lengths[-1])
-
-        # find gap position before seq_index
-        index = numpy.searchsorted(gap_pos, seq_index, side="left")
-        if seq_index < gap_pos[index]:
-            gap_lengths = cum_lengths[index - 1] if index else 0
-        else:
-            gap_lengths = cum_lengths[index]
-
-        return int(seq_index + gap_lengths)
+        return int(index)
 
     def get_seq_index(self, align_index: int) -> int:
         """converts alignment index to sequence index"""
+        val = align_to_seq_index(
+            self.gap_pos,
+            self.cum_gap_lengths,
+            len(self),
+            self.num_gaps,
+            align_index,
+        )
         # NOTE I explicitly cast all returned values to python int's due to
         # need for json serialisation, which does not support numpy int classes
-        if align_index < 0:
-            align_index = len(self) + align_index
-        if align_index < 0:
-            raise IndexError(f"{align_index} align_index beyond limit")
-
-        if not self.num_gaps or align_index < self.gap_pos[0]:
-            return align_index
-
-        # these are alignment indices for gaps
-        cum_lengths = self.cum_gap_lengths
-        gap_starts, gap_ends = _gap_spans(self.gap_pos, cum_lengths)
-        if align_index >= gap_ends[-1]:
-            return int(align_index - cum_lengths[-1])
-
-        index = numpy.searchsorted(gap_ends, align_index, side="left")
-        if align_index < gap_starts[index]:
-            # before the gap at index
-            return int(align_index - cum_lengths[index - 1])
-
-        if align_index == gap_ends[index]:
-            # after the gap at index
-            return int(align_index - cum_lengths[index])
-
-        if gap_starts[index] <= align_index < gap_ends[index]:
-            # within the gap at index
-            # so the gap insertion position is the sequence position
-            return int(self.gap_pos[index])
+        return int(val)
 
     def __len__(self) -> int:
         length_gaps = self.cum_gap_lengths[-1] if self.num_gaps else 0
diff --git a/src/cogent3/core/new_alignment.py b/src/cogent3/core/new_alignment.py
index b32f10a9a..2b68e8909 100644
--- a/src/cogent3/core/new_alignment.py
+++ b/src/cogent3/core/new_alignment.py
@@ -18,7 +18,13 @@ import numpy
 
 from cogent3 import get_app
 from cogent3._version import __version__
-from cogent3.core import new_alphabet, new_genetic_code, new_moltype, new_sequence
+from cogent3.core import (
+    location,
+    new_alphabet,
+    new_genetic_code,
+    new_moltype,
+    new_sequence,
+)
 from cogent3.core.annotation import Feature
 from cogent3.core.annotation_db import (
     BasicAnnotationDb,
@@ -848,10 +854,7 @@ class SequenceCollection:
         as_array
             if True, sequences are returned as numpy arrays, otherwise as strings
         """
-        return {
-            name: (numpy.array(self.seqs[name]) if as_array else str(self.seqs[name]))
-            for name in self.names
-        }
+        return {s.name: (numpy.array(s) if as_array else str(s)) for s in self.seqs}
 
     def to_rich_dict(self) -> dict[str, str | dict[str, str]]:
         """returns a json serialisable dict
@@ -2540,7 +2543,7 @@ def _(
 
 
 @singledispatch
-def seq_to_gap_coords(
+def decompose_gapped_seq(
     seq: typing.union[StrORBytesORArray, new_sequence.Sequence],
     *,
     alphabet: new_alphabet.AlphabetABC,
@@ -2549,11 +2552,13 @@ def seq_to_gap_coords(
     Takes a sequence with (or without) gaps and returns an ungapped sequence
     and a map of the position and length of gaps in the original parent sequence
     """
-    raise NotImplementedError(f"seq_to_gap_coords not implemented for type {type(seq)}")
+    raise NotImplementedError(
+        f"decompose_gapped_seq not implemented for type {type(seq)}"
+    )
 
 
 @numba.jit
-def gapped_seq_to_ungapped_and_gaps(
+def decompose_gapped_seq_array(
     seq: numpy.ndarray,
     gap_index: int,
 ) -> tuple[numpy.ndarray, numpy.ndarray]:  # pragma: no cover
@@ -2573,7 +2578,7 @@ def gapped_seq_to_ungapped_and_gaps(
 
     Notes
     -----
-    being called by seq_to_gap_coords
+    being called by decompose_gapped_seq
     """
     seqlen = len(seq)
     working = numpy.empty((seqlen, numpy.int64(2)), dtype=numpy.int64)
@@ -2615,39 +2620,70 @@ def gapped_seq_to_ungapped_and_gaps(
     return ungapped, gap_coords
 
 
-@seq_to_gap_coords.register
+@numba.jit
+def compose_gapped_seq(
+    ungapped_seq: numpy.ndarray, gaps: numpy.ndarray, gap_index: int
+) -> numpy.ndarray:  # pragma: no cover
+    """reconstruct a gapped sequence from an ungapped sequence and gap data"""
+    if not len(gaps):
+        return ungapped_seq
+
+    gapped_len = len(ungapped_seq) + gaps[-1, 1]
+
+    gapped_seq = numpy.empty(gapped_len, dtype=ungapped_seq.dtype)
+
+    pos = 0
+    ungapped_pos = 0
+    prev_gap_len = 0
+    for gap_pos, cum_gap_len in gaps:
+        gap_len = cum_gap_len - prev_gap_len
+        prev_gap_len = cum_gap_len
+
+        gapped_seq[pos : pos + gap_pos - ungapped_pos] = ungapped_seq[
+            ungapped_pos:gap_pos
+        ]
+        pos += gap_pos - ungapped_pos
+        ungapped_pos = gap_pos
+
+        gapped_seq[pos : pos + gap_len] = gap_index
+        pos += gap_len
+
+    gapped_seq[pos:] = ungapped_seq[ungapped_pos:]
+
+    return gapped_seq
+
+
+@decompose_gapped_seq.register
 def _(
     seq: numpy.ndarray,
     *,
     alphabet: new_alphabet.AlphabetABC,
 ) -> tuple[numpy.ndarray, numpy.ndarray]:
-    return gapped_seq_to_ungapped_and_gaps(
-        seq.astype(alphabet.dtype), alphabet.gap_index
-    )
+    return decompose_gapped_seq_array(seq.astype(alphabet.dtype), alphabet.gap_index)
 
 
-@seq_to_gap_coords.register
+@decompose_gapped_seq.register
 def _(
     seq: str, *, alphabet: new_alphabet.AlphabetABC
 ) -> tuple[numpy.ndarray, numpy.ndarray]:
     if not alphabet.is_valid(seq):
         raise new_alphabet.AlphabetError(f"Sequence is invalid for alphabet {alphabet}")
 
-    return seq_to_gap_coords(alphabet.to_indices(seq), alphabet=alphabet)
+    return decompose_gapped_seq(alphabet.to_indices(seq), alphabet=alphabet)
 
 
-@seq_to_gap_coords.register
+@decompose_gapped_seq.register
 def _(
     seq: bytes, *, alphabet: new_alphabet.AlphabetABC
 ) -> tuple[numpy.ndarray, numpy.ndarray]:
-    return seq_to_gap_coords(seq.decode("utf-8"), alphabet=alphabet)
+    return decompose_gapped_seq(seq.decode("utf-8"), alphabet=alphabet)
 
 
-@seq_to_gap_coords.register
+@decompose_gapped_seq.register
 def _(
     seq: new_sequence.Sequence, *, alphabet: new_alphabet.AlphabetABC
 ) -> tuple[numpy.ndarray, numpy.ndarray]:
-    return seq_to_gap_coords(numpy.array(seq), alphabet=alphabet)
+    return decompose_gapped_seq(numpy.array(seq), alphabet=alphabet)
 
 
 class Aligned:
@@ -3030,7 +3066,7 @@ class AlignedSeqsData(AlignedSeqsDataABC):
         seqs = {}
         gaps = {}
         for name, seq in data.items():
-            seq, gap_map = seq_to_gap_coords(seq, alphabet=alphabet)
+            seq, gap_map = decompose_gapped_seq(seq, alphabet=alphabet)
             seq = alphabet.to_indices(seq)
             seq.flags.writeable = False
             gap_map.flags.writeable = False
@@ -3099,7 +3135,7 @@ class AlignedSeqsData(AlignedSeqsDataABC):
         seqs = {}
         gaps = {}
         for name, seq in zip(names, data):
-            seq, gap_map = seq_to_gap_coords(seq, alphabet=alphabet)
+            seq, gap_map = decompose_gapped_seq(seq, alphabet=alphabet)
             seq.flags.writeable = False
             gap_map.flags.writeable = False
             seqs[name], gaps[name] = seq, gap_map
@@ -3191,64 +3227,17 @@ class AlignedSeqsData(AlignedSeqsDataABC):
         """Return sequence data corresponding to seqid as an array of indices.
         start/stop are in alignment coordinates. Includes gaps.
         """
-
-        # refactor: design
-        # should this pass the data sliced data to the ADV which then applies the step???
-
-        step = 1 if step is None else step
-        func = _input_vals_pos_step if step > 0 else _input_vals_neg_step
-        start, stop, step = func(self.align_len, start, stop, step)
-        seq = self._seqs[seqid]
-
-        # if there's no gaps, just slice the sequence
-        if len(seq) == self.align_len:
-            return seq[start:stop:step]
-
-        unknown = self.alphabet.gap_index
-
-        # there is only gaps
-        if len(seq) == 0:
-            return numpy.full(
-                int(-((stop - start) // -step)), unknown, dtype=numpy.uint8
-            )
-
+        # reminder: design
+        # this is recreating the orignal seq, so will be memory and compute inefficient
+        # better to build the smallest gapped seq represented by the start/stop
+        # determine the seq coordinates for slicing the ungapped sequence then apply the step
+        # note with the latter approach, if the step is negative, need to apply the step as
+        # a positive integer than return the reversed array
+        gap_index = self.alphabet.gap_index
         gaps = self._gaps[seqid]
-        indel_map = IndelMap(
-            gap_pos=gaps[:, 0], cum_gap_lengths=gaps[:, 1], parent_length=len(seq)
-        )
-
-        if step < 0:
-            seq = seq[::-1]
-            indel_map = indel_map[::-1]
-            start = -start - 1
-            stop = -stop - 1
-            step = -step
-
-        # convert from alignment coordinates to sequences coordinates
-        # selecting only the indices that are in the frame of the stride
-        # refactor: design
-        # only select indices if we have values for start/stop/step
-        indices = []
-        for span in indel_map.nongap():
-            indices += [
-                indel_map.get_seq_index(x)
-                for x in range(span.start, span.end)
-                if start <= x < stop and (x - start) % step == 0
-            ]
-
-        seq_sliced = seq[indices]
-        map_sliced = indel_map[start:stop:step]
-        concat = numpy.array([], dtype=numpy.uint8)
-
-        # iterate through spans
-        for span in map_sliced.spans:
-            if span.lost:
-                seg = numpy.full(span.length, unknown, dtype=numpy.uint8)
-            else:
-                seg = seq_sliced[span.start : span.end]
-            concat = numpy.concatenate((concat, seg))
-
-        return concat
+        seq = self._seqs[seqid]
+        gapped = compose_gapped_seq(seq, gaps, gap_index)
+        return gapped[start:stop:step]
 
     def get_seq_str(
         self,
@@ -3318,8 +3307,12 @@ class AlignedSeqsData(AlignedSeqsDataABC):
         step: OptInt = None,
     ) -> tuple[dict, dict]:
         seqs = {}
+        # redesign
+        # if gaps exist, don't go via gapped seq
+        # convert alignment coords into sequence coords using the location.align_to_seq_index function
+        # this means we will need to convert coordinates to a plus strand slice
         for name in name_map.values():
-            # use get_gapped_seq_array because we want to slice in alignment coordinates
+            # using get_gapped_seq_array so we can slice in alignment coordinates
             seq = self.get_gapped_seq_array(
                 seqid=name, start=start, stop=stop, step=step
             )
@@ -3358,7 +3351,7 @@ class AlignedSeqsData(AlignedSeqsDataABC):
         new_seqs = {}
         new_gaps = {}
         for name, seq in seqs.items():
-            seq, gap_map = seq_to_gap_coords(seq, alphabet=self.alphabet)
+            seq, gap_map = decompose_gapped_seq(seq, alphabet=self.alphabet)
             seq = self.alphabet.to_indices(seq)
             seq.flags.writeable = False
             gap_map.flags.writeable = False
diff --git a/src/cogent3/core/sequence.py b/src/cogent3/core/sequence.py
index 058502587..dcb0026f0 100644
--- a/src/cogent3/core/sequence.py
+++ b/src/cogent3/core/sequence.py
@@ -31,6 +31,7 @@ from numpy import (
     array,
     compress,
     floating,
+    int32,
     integer,
     issubdtype,
     logical_not,


Test patch:
diff --git a/tests/test_core/test_new_alignment.py b/tests/test_core/test_new_alignment.py
index 987a2a7f7..ac7515ca1 100644
--- a/tests/test_core/test_new_alignment.py
+++ b/tests/test_core/test_new_alignment.py
@@ -2732,10 +2732,10 @@ def gap_seqs():
 
 
 @pytest.mark.parametrize("i", range(3))
-def test_seq_to_gap_coords_sequences(gap_seqs, i, dna_alphabet):
+def test_decompose_gapped_seq_sequences(gap_seqs, i, dna_alphabet):
     seq, gap_coords = gap_seqs[i]
     dna = new_moltype.get_moltype("dna")
-    got_ungapped, got_map = new_alignment.seq_to_gap_coords(
+    got_ungapped, got_map = new_alignment.decompose_gapped_seq(
         dna.make_seq(seq=seq), alphabet=dna_alphabet
     )
     expect = dna_alphabet.to_indices(seq.replace("-", ""))
@@ -2744,18 +2744,20 @@ def test_seq_to_gap_coords_sequences(gap_seqs, i, dna_alphabet):
 
 
 @pytest.mark.parametrize("i", range(3))
-def test_seq_to_gap_coords_str(gap_seqs, i, dna_alphabet):
+def test_decompose_gapped_seq_str(gap_seqs, i, dna_alphabet):
     seq, gap_coords = gap_seqs[i]
-    got_ungapped, got_map = new_alignment.seq_to_gap_coords(seq, alphabet=dna_alphabet)
+    got_ungapped, got_map = new_alignment.decompose_gapped_seq(
+        seq, alphabet=dna_alphabet
+    )
     expect = dna_alphabet.to_indices(seq.replace("-", ""))
     assert numpy.array_equal(got_ungapped, expect)
     assert numpy.array_equal(got_map, gap_coords)
 
 
-def test_seq_to_gap_coords_str_all_gaps(dna_alphabet):
+def test_decompose_gapped_seq_str_all_gaps(dna_alphabet):
     parent_seq = "-----"
     expect_gaplen = numpy.array([len(parent_seq)])
-    got_ungap, got_map = new_alignment.seq_to_gap_coords(
+    got_ungap, got_map = new_alignment.decompose_gapped_seq(
         parent_seq, alphabet=dna_alphabet
     )
     expect = numpy.array([])
@@ -2763,9 +2765,9 @@ def test_seq_to_gap_coords_str_all_gaps(dna_alphabet):
     assert got_map[:, 1] == expect_gaplen
 
 
-def test_seq_to_gap_coords_str_no_gaps(dna_alphabet):
+def test_decompose_gapped_seq_str_no_gaps(dna_alphabet):
     parent_seq = "ACTGC"
-    got_ungap, got_map = new_alignment.seq_to_gap_coords(
+    got_ungap, got_map = new_alignment.decompose_gapped_seq(
         parent_seq, alphabet=dna_alphabet
     )
     expect = dna_alphabet.to_indices(parent_seq)
@@ -2773,18 +2775,18 @@ def test_seq_to_gap_coords_str_no_gaps(dna_alphabet):
     assert got_map.size == 0
 
 
-def test_seq_to_gap_coords_arr_all_gaps(dna_alphabet):
+def test_decompose_gapped_seq_arr_all_gaps(dna_alphabet):
     parent_seq = dna_alphabet.to_indices("-----")
-    got_ungap, got_map = new_alignment.seq_to_gap_coords(
+    got_ungap, got_map = new_alignment.decompose_gapped_seq(
         parent_seq, alphabet=dna_alphabet
     )
     assert got_ungap.size == 0
     assert numpy.array_equal(got_map, numpy.array([[0, 5]]))
 
 
-def test_seq_to_gap_coords_arr_no_gaps(dna_alphabet):
+def test_decompose_gapped_seq_arr_no_gaps(dna_alphabet):
     parent_seq = dna_alphabet.to_indices("ACTGC")
-    got_ungap, got_empty_arr = new_alignment.seq_to_gap_coords(
+    got_ungap, got_empty_arr = new_alignment.decompose_gapped_seq(
         parent_seq, alphabet=dna_alphabet
     )
     assert numpy.array_equal(got_ungap, parent_seq)
@@ -2792,28 +2794,30 @@ def test_seq_to_gap_coords_arr_no_gaps(dna_alphabet):
 
 
 @pytest.mark.parametrize("i", range(3))
-def test_seq_to_gap_coords_arr(gap_seqs, i, dna_alphabet):
+def test_decompose_gapped_seq_arr(gap_seqs, i, dna_alphabet):
     seq, gap_coords = gap_seqs[i]
     seq = dna_alphabet.to_indices(seq)
-    got_ungapped, got_map = new_alignment.seq_to_gap_coords(seq, alphabet=dna_alphabet)
+    got_ungapped, got_map = new_alignment.decompose_gapped_seq(
+        seq, alphabet=dna_alphabet
+    )
     expect = seq[seq != 4]  # gap_char = 4
     assert numpy.array_equal(got_ungapped, expect)
     assert numpy.array_equal(got_map, gap_coords)
 
 
 @pytest.mark.parametrize("i", range(3))
-def test_seq_to_gap_coords_arr_dispatch_equal(gap_seqs, i, dna_alphabet):
-    """seq_to_gap_coords should return the same gap coords when input is a string/array/bytes"""
+def test_decompose_gapped_seq_arr_dispatch_equal(gap_seqs, i, dna_alphabet):
+    """decompose_gapped_seq should return the same gap coords when input is a string/array/bytes"""
     seq_str, _ = gap_seqs[i]
     seq_array = dna_alphabet.to_indices(seq_str)
     seq_bytes = dna_alphabet.array_to_bytes(seq_array)
-    seq_from_str, gaps_from_str = new_alignment.seq_to_gap_coords(
+    seq_from_str, gaps_from_str = new_alignment.decompose_gapped_seq(
         seq_str, alphabet=dna_alphabet
     )
-    seq_from_array, gaps_from_arr = new_alignment.seq_to_gap_coords(
+    seq_from_array, gaps_from_arr = new_alignment.decompose_gapped_seq(
         seq_array, alphabet=dna_alphabet
     )
-    seq_from_bytes, gaps_from_bytes = new_alignment.seq_to_gap_coords(
+    seq_from_bytes, gaps_from_bytes = new_alignment.decompose_gapped_seq(
         seq_bytes, alphabet=dna_alphabet
     )
     assert numpy.array_equal(gaps_from_str, gaps_from_arr)
@@ -2983,11 +2987,11 @@ def test_aligned_seqs_data_init_gapped(
     }
 
     seq_data = {
-        name: new_alignment.seq_to_gap_coords(seq, alphabet=dna_alphabet)[0]
+        name: new_alignment.decompose_gapped_seq(seq, alphabet=dna_alphabet)[0]
         for name, seq in typed_data.items()
     }
     gap_data = {
-        name: new_alignment.seq_to_gap_coords(seq, alphabet=dna_alphabet)[1]
+        name: new_alignment.decompose_gapped_seq(seq, alphabet=dna_alphabet)[1]
         for name, seq in typed_data.items()
     }
     asd = new_alignment.AlignedSeqsData(
@@ -3008,11 +3012,11 @@ def test_aligned_seqs_data_unequal_seqlens_raises(data_type, dna_alphabet, dna_m
         _ = new_alignment.AlignedSeqsData.from_seqs(data=data, alphabet=dna_alphabet)
     # directly creating an AlignedSeqsData object should also raise an error
     seq_data = {
-        name: new_alignment.seq_to_gap_coords(seq, alphabet=dna_alphabet)[0]
+        name: new_alignment.decompose_gapped_seq(seq, alphabet=dna_alphabet)[0]
         for name, seq in data.items()
     }
     gap_data = {
-        name: new_alignment.seq_to_gap_coords(seq, alphabet=dna_alphabet)[1]
+        name: new_alignment.decompose_gapped_seq(seq, alphabet=dna_alphabet)[1]
         for name, seq in data.items()
     }
     with pytest.raises(ValueError):
@@ -4676,3 +4680,45 @@ def test_indexing_seqs_iter(names_seqs, func):
     got = list(map(str, obj.seqs))
     expect = list(seqs)
     assert got == expect
+
+
+@pytest.mark.parametrize(
+    "gapped_seq",
+    (
+        "ACGGCTA",  # No gaps
+        "AC--GT-AGC",  # With gaps
+        "------",  # All gaps
+        # Edge cases
+        "-A-C-G-T",
+        "A-C-G-T-",
+    ),
+)
+def test_gapped_seq_round_trip(gapped_seq):
+    gapped_seq = new_moltype.DNA.gapped_alphabet.to_indices(gapped_seq)
+    # split into components
+    ungapped_seq, gaps = new_alignment.decompose_gapped_seq_array(gapped_seq, 4)
+
+    # Recreate the gapped sequence from the ungapped sequence and gaps
+    recreated_gapped_seq = new_alignment.compose_gapped_seq(ungapped_seq, gaps, 4)
+
+    # Test the output of gapped_seq_from_components against the original sequence
+    numpy.testing.assert_array_equal(recreated_gapped_seq, gapped_seq)
+
+
+def test_asd_get_gapped_seq(aligned_dict, dna_alphabet):
+    a_slice = slice(None, 4, 2)
+    seqid = "seq4"
+    orig_seq = aligned_dict[seqid][a_slice]
+    orig_array = dna_alphabet.to_indices(orig_seq)
+    orig_revd = orig_array[::-1]
+
+    asd = new_alignment.AlignedSeqsData.from_seqs(
+        data=aligned_dict, alphabet=new_moltype.DNA.most_degen_alphabet()
+    )
+    view = asd.get_view(seqid)
+    fwd = view[a_slice]
+    rev = fwd[::-1]
+    fwd_gapped = fwd.gapped_array_value
+    numpy.testing.assert_array_equal(fwd_gapped, orig_array)
+    rev_gapped = rev.gapped_array_value
+    numpy.testing.assert_array_equal(rev_gapped, orig_revd)