Repository: cogent3/cogent3
Pull request #1817 (2024-04-21 22:45:20)

Refactor of map
IndelMap now uses numpy arrays and the API is reduced compared to the previous Map implementation.


Patch:
diff --git a/changelog.d/20240329_152534_Gavin.Huttley.md b/changelog.d/20240329_152534_Gavin.Huttley.md
index 779132349..b8b15800d 100644
--- a/changelog.d/20240329_152534_Gavin.Huttley.md
+++ b/changelog.d/20240329_152534_Gavin.Huttley.md
@@ -22,12 +22,11 @@ Uncomment the section that is right (remove the HTML comment wrapper).
   length is set to None. Think of this tree as the true "evolutionary tree".
   Thanks to Von Bing Yap for suggesting this!
 
-<!--
 ### BUG
 
-- A bullet item for the BUG category.
+- Alignment.quick_tree() default settings referred to a distance calculator
+  that has been removed. Now updated to the pdist calculator.
 
--->
 <!--
 ### DOC
 
diff --git a/requirements.txt b/requirements.txt
index 03b4568a2..8b70ef5b3 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,13 +1,13 @@
 chardet==5.2.0
-pillow==10.2.0
+pillow==10.3.0
 psutil==5.9.8
 numpy==1.26.4
 numba>0.48.0;python_version<'3.9'
 numba>0.53; python_version>='3.9'
 numba>0.54; python_version>='3.10'
 numba>=0.57.0; python_version=='3.11'
-scipy==1.12.0
-nox==2024.3.2
+scipy==1.13.0
+nox==2024.4.15
 plotly==5.20.0
 typing_extensions>=3.7
 .[dev]
diff --git a/src/cogent3/align/pairwise.py b/src/cogent3/align/pairwise.py
index c517092e1..2214f0edb 100644
--- a/src/cogent3/align/pairwise.py
+++ b/src/cogent3/align/pairwise.py
@@ -13,6 +13,7 @@ import numpy
 
 from cogent3.align.traceback import alignment_traceback, map_traceback
 from cogent3.core.alignment import Aligned
+from cogent3.core.location import IndelMap
 from cogent3.evolve.likelihood_tree import LikelihoodTreeEdge
 from cogent3.util.misc import ascontiguousarray
 
@@ -477,8 +478,20 @@ class AlignablePOG(_Alignable):
         aligneds = []
         for dim, child in enumerate(children):
             for seq_name, aligned in child.aligneds:
-                aligned = aligned.remapped_to((maps[dim] * word_length).inverse())
-                aligneds.append((seq_name, aligned))
+                # if word_length != 1 then maps are forced to have
+                # sequences lengths that are modulo word_length
+                new_map = aligned.map.merge_maps(
+                    maps[dim] * word_length,
+                    parent_length=maps[dim].parent_length * word_length,
+                )
+                # Likewise, if the data is not modulo word_length,
+                # it is trimmed to match
+                data = (
+                    aligned.data
+                    if new_map.parent_length == len(aligned.data)
+                    else aligned.data[: new_map.parent_length]
+                )
+                aligneds.append((seq_name, Aligned(new_map, data)))
         return aligneds
 
     def backward(self):
@@ -526,7 +539,11 @@ class AlignableSeq(_Alignable):
         _Alignable.__init__(self, leaf)
         if hasattr(leaf, "sequence"):
             self.seq = leaf.sequence
-            aligned = Aligned([(0, len(self.seq))], self.seq, len(self.seq))
+            seqlen = len(self.seq)
+            imap = IndelMap.from_aligned_segments(
+                locations=[(0, seqlen)], aligned_length=seqlen
+            )
+            aligned = Aligned(imap, self.seq)
             self.aligneds = [(self.leaf.edge_name, aligned)]
         self.max_preds = 1
         self._pog = None
diff --git a/src/cogent3/align/traceback.py b/src/cogent3/align/traceback.py
index 16a8d0ceb..4a0ca3dba 100644
--- a/src/cogent3/align/traceback.py
+++ b/src/cogent3/align/traceback.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 """Conversion of dynamic program results ("arrays of arrows") into gap vectors,
 gapped sequences or Cogent Alignment objects"""
 
@@ -52,18 +51,17 @@ def gap_traceback(aligned_positions):
         if consuming[dimension]:
             gv.append(a)
         gap_vectors[dimension] = [(gv[i], gv[i + 1]) for i in range(0, len(gv), 2)]
-    return (starts, ends, gap_vectors, a)
+    return starts, ends, gap_vectors, a
 
 
 def map_traceback(aligned_positions):
     # using IndelMap's to keep track of gaps for indel alignment
-    (starts, ends, gap_vectors, alignment_len) = gap_traceback(aligned_positions)
-    # print 'gv', gap_vectors
+    starts, ends, gap_vectors, alignment_len = gap_traceback(aligned_positions)
     maps = [
-        IndelMap(locations=gv, parent_length=alignment_len).inverse()
+        IndelMap.from_aligned_segments(locations=gv, aligned_length=alignment_len)
         for gv in gap_vectors
     ]
-    return (starts, ends, maps)
+    return starts, ends, maps
 
 
 def alignment_traceback(seqs, aligned_positions, word_length):
diff --git a/src/cogent3/app/align.py b/src/cogent3/app/align.py
index 94a2200be..04a2abdd4 100644
--- a/src/cogent3/app/align.py
+++ b/src/cogent3/app/align.py
@@ -171,7 +171,7 @@ def _merged_gaps(a_gaps: dict, b_gaps: dict) -> dict:
     function to 'max'. Use 'sum' when the gaps derive from different
     sequences.
     """
-
+    # todo convert to using IndelMap functions
     if not a_gaps:
         return b_gaps
 
@@ -251,7 +251,8 @@ def _gaps_for_injection(other_seq_gaps: dict, refseq_gaps: dict, seqlen: int) ->
     # 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
+    # todo convert these functions to using IndelMap and the numpy set
+    #  operation functions
     all_gaps = {}
     all_gaps.update(other_seq_gaps)
     for gap_pos, gap_length in sorted(refseq_gaps.items()):
@@ -304,8 +305,7 @@ def pairwise_to_multiple(pwise, ref_seq, moltype, info=None):
         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:
+        if inject := _gaps_for_injection(other_gaps, diff_gaps, len(other_seq.data)):
             m = gap_coords_to_map(inject, len(other_seq.data))
             other_seq = Aligned(m, other_seq.data)
 
diff --git a/src/cogent3/core/alignment.py b/src/cogent3/core/alignment.py
index 3ba2a0f6a..e34cf8b79 100644
--- a/src/cogent3/core/alignment.py
+++ b/src/cogent3/core/alignment.py
@@ -30,7 +30,7 @@ import warnings
 
 from collections import Counter, defaultdict
 from copy import deepcopy
-from functools import total_ordering
+from functools import singledispatchmethod, total_ordering
 from itertools import combinations
 from typing import Any, Iterable, Iterator, List, Optional, Tuple, Union
 
@@ -2136,11 +2136,10 @@ class Aligned:
     """One sequence in an alignment, a map between alignment coordinates and
     sequence coordinates"""
 
-    def __init__(self, map, data, length=None):
-        # Unlike the normal map constructor, here we take a list of pairs of
-        # alignment coordinates, NOT a list of pairs of sequence coordinates
-        if isinstance(map, list):
-            map = IndelMap(locations=map, parent_length=length).inverse()
+    @c3warn.deprecated_args(
+        "2024.9", reason="now requires an IndelMap", discontinued=["length"]
+    )
+    def __init__(self, map, data):
         self.map = map
         self.data = data
         if hasattr(data, "info"):
@@ -2188,12 +2187,8 @@ class Aligned:
             if strand == "-" or exclude_annotations:
                 new_seq.annotation_db = None
             else:
-                new_seq.annotation_offset = self.map.start
                 new_seq.annotation_db = db
-            new_map = self.map.zeroed()
-        else:
-            new_map = self.map
-
+        new_map = self.map
         return self.__class__(new_map, new_seq)
 
     def __repr__(self):
@@ -2243,19 +2238,52 @@ class Aligned:
             map, seq = seq.parse_out_gaps()
         return Aligned(map, seq)
 
-    def __getitem__(self, slice):
+    @singledispatchmethod
+    def __getitem__(self, span: int):
+        raise NotImplementedError(f"unexpected {type(span)}")
+
+    @__getitem__.register
+    def _(self, span: int):
+        return self[span : span + 1]
+
+    @__getitem__.register
+    def _(self, span: FeatureMap):
+        # we assume the feature map is in align coordinates
+        start, end = span.start, span.end
+        if span.useful and start > end:
+            empty = numpy.array([], dtype=self.map.gap_pos.dtype)
+            im = IndelMap(gap_pos=empty, cum_gap_lengths=empty, parent_length=0)
+            data = self.data[:0]
+        elif span.useful and len(list(span.spans)) == 1:
+            im = self.map[start:end]
+            seq_start = self.map.get_seq_index(start)
+            seq_end = self.map.get_seq_index(end)
+            data = self.data[seq_start:seq_end]
+        elif not span.useful:
+            im = self.map[start:end]
+            data = self.data[:0]
+        else:
+            # multiple spans
+            align_coords = span.get_coordinates()
+            im = self.map.joined_segments(align_coords)
+            seq_map = self.map.make_seq_feature_map(span, include_gaps=False)
+            data = self.data.gapped_by_map(seq_map)
+
+        return Aligned(im, data)
+
+    @__getitem__.register
+    def _(self, span: slice):
         # todo we need to get the sequence coordinates that slice corresponds to
-        #  so we can update the self.data, plus we will need to zero new_map
-        new_map = self.map[slice]
-        data = (
-            self.data[new_map.start : new_map.end] if new_map.useful else self.data[:0]
-        )
-        if new_map.useful and new_map.start > new_map.end:
+        #  so we can update the self.data
+        new_map = self.map[span]
+        seq_start = self.map.get_seq_index(span.start or 0)
+        seq_end = self.map.get_seq_index(span.stop or len(self))
+        data = self.data[seq_start:seq_end] if new_map.useful else self.data[:0]
+        if new_map.useful and seq_start > seq_end:
             # For now, a reverse slice means we should have an empty sequence
             # todo modify this clause if a negative step is ever allowed
             new_map = new_map.__class__(locations=(), parent_length=len(self.data))
-        elif new_map.useful:
-            new_map = new_map.zeroed()
+
         return Aligned(new_map, data)
 
     def rc(self):
@@ -2274,12 +2302,9 @@ class Aligned:
         return Aligned(self.map, self.data.to_dna())
 
     def to_rich_dict(self):
-        coords = self.map.get_covering_span().get_coordinates()
-        if len(coords) != 1:
-            raise NotImplementedError
         data = dict(version=__version__, type=get_object_provenance(self))
         # we are resetting the stored data to begin at 0
-        data["map_init"] = self.map.zeroed().to_rich_dict()
+        data["map_init"] = self.map.to_rich_dict()
         data["seq_init"] = self.data.to_rich_dict(exclude_annotations=True)
         return data
 
@@ -2318,7 +2343,7 @@ class Aligned:
     def make_feature(self, feature: FeatureDataType, alignment: "Alignment") -> Feature:
         """returns a feature, not written into annotation_db"""
         annot = self.data.make_feature(feature)
-        inverted = self.map.inverse()
+        inverted = self.map.to_feature_map().inverse()
         # todo should indicate whether tidy or not
         return annot.remapped_to(alignment, inverted)
 
@@ -3285,7 +3310,7 @@ class AlignmentI(object):
     @extend_docstring_from(distance_matrix, pre=False)
     def quick_tree(
         self,
-        calc="percent",
+        calc="pdist",
         bootstrap=None,
         drop_invalid=False,
         show_progress=False,
@@ -4721,18 +4746,18 @@ class Alignment(AlignmentI, SequenceCollection):
 
     def _mapped(self, slicemap):
         align = []
-        for name in self.names:
-            align.append((name, self.named_seqs[name][slicemap]))
+        for seq in self.seqs:
+            sliced = seq[slicemap]
+            align.append((sliced.name, sliced))
         return self.__class__(moltype=self.moltype, data=align, info=self.info)
 
     def gapped_by_map(self, keep, **kwargs):
         # keep is a Map
+        # seqs = [seq[keep] for seq in self.seqs]
         seqs = []
-        for seq_name in self.names:
-            aligned = self.named_seqs[seq_name]
-            seqmap = aligned.map[keep]
-            seq = aligned.data.gapped_by_map(seqmap)
-            seqs.append((seq_name, seq))
+        for seq in self.seqs:
+            selected = seq[keep]
+            seqs.append(selected)
         return self.__class__(moltype=self.moltype, data=seqs, **kwargs)
 
     def get_projected_feature(self, *, seqid: str, feature: Feature) -> Feature:
@@ -4834,8 +4859,9 @@ class Alignment(AlignmentI, SequenceCollection):
             return None
 
         locations = [(gv[i], gv[i + 1]) for i in range(0, len(gv), 2)]
+        # these are alignment coordinate locations
+        keep = FeatureMap.from_locations(locations=locations, parent_length=len(self))
 
-        keep = IndelMap(locations=locations, parent_length=len(self))
         return self.gapped_by_map(keep, info=self.info)
 
     def get_seq(self, seqname):
@@ -5231,7 +5257,9 @@ class Alignment(AlignmentI, SequenceCollection):
         # there's no sequence to bind to, the feature is directly on self
         revd = feature.pop("strand", None) == "-"
         feature["strand"] = "-" if revd else "+"
-        fmap = FeatureMap(parent_length=len(self), locations=feature.pop("spans"))
+        fmap = FeatureMap.from_locations(
+            locations=feature.pop("spans"), parent_length=len(self)
+        )
         if revd:
             fmap = fmap.nucleic_reversed()
         return Feature(parent=self, map=fmap, **feature)
@@ -5367,7 +5395,7 @@ class Alignment(AlignmentI, SequenceCollection):
             if feature["seqid"]:
                 raise RuntimeError(f"{on_alignment=} {feature=}")
             if seq_map is None:
-                seq_map = self.seqs[0].map
+                seq_map = self.seqs[0].map.to_feature_map()
                 *_, strand = self.seqs[0].data.parent_coordinates()
 
             spans = numpy.array(feature["spans"])
diff --git a/src/cogent3/core/annotation_db.py b/src/cogent3/core/annotation_db.py
index 60bd064ed..83fb68885 100644
--- a/src/cogent3/core/annotation_db.py
+++ b/src/cogent3/core/annotation_db.py
@@ -333,7 +333,7 @@ def _matching_conditions(
         sql.append(" AND ".join(conds))
         vals = tuple(vals)
 
-    if isinstance(start, int) and isinstance(stop, int):
+    if start is not None and stop is not None:
         if allow_partial:
             # allow matches that overlap the segment
             cond = [
@@ -347,11 +347,11 @@ def _matching_conditions(
             # only matches within bounds
             cond = f"start >= {start} AND stop <= {stop}"
         sql.append(f"({cond})")
-    elif isinstance(start, int):
+    elif start is not None:
         # if query has no stop, then any feature containing start
         cond = f"(start <= {start} AND {start} < stop)"
         sql.append(f"({cond})")
-    elif isinstance(stop, int):
+    elif stop is not None:
         # if query has no start, then any feature containing stop
         cond = f"(start <= {stop} AND {stop} < stop)"
         sql.append(f"({cond})")
diff --git a/src/cogent3/core/location.py b/src/cogent3/core/location.py
index 476b93e7c..4f8f274c6 100644
--- a/src/cogent3/core/location.py
+++ b/src/cogent3/core/location.py
@@ -21,7 +21,9 @@ former contains a and b while the latter contains c.
 """
 import copy
 import dataclasses
+import functools
 import inspect
+import json
 
 from abc import ABC, abstractmethod
 from bisect import bisect_left, bisect_right
@@ -29,10 +31,13 @@ from functools import total_ordering
 from itertools import chain
 from typing import Iterator, List, Optional, Sequence, Tuple, Union
 
-from numpy import array, ndarray
+import numpy
+
+from numpy.typing import NDArray
 
 from cogent3._version import __version__
 from cogent3.util import warning as c3warn
+from cogent3.util.deserialise import register_deserialiser
 from cogent3.util.misc import (
     ClassChecker,
     ConstrainedList,
@@ -54,21 +59,6 @@ def _norm_index(i, length, default):
     return min(max(i, 0), length)
 
 
-def _norm_slice(index, length):
-    """_norm_slice(slice(1, -2, 3), 10) -> (1,8,3)"""
-    if isinstance(index, slice):
-        start = _norm_index(index.start, length, 0)
-        end = _norm_index(index.stop, length, length)
-        return (start, end, index.step)
-    else:
-        start = index
-        if start < 0:
-            start += length
-        if start >= length:
-            raise IndexError(index)
-        return (start, start + 1, 1)
-
-
 def as_map(slice, length, cls):
     """Take anything that might be used as a subscript: Integer, Slice,
     or MapABC, and return cls."""
@@ -78,13 +68,15 @@ def as_map(slice, length, cls):
         for i in slice:
             spans.extend(as_map(i, length, cls).spans)
         return cls(spans=spans, parent_length=length)
-    elif isinstance(slice, (FeatureMap, IndelMap)):
+    elif isinstance(slice, (FeatureMap, IndelMap, Map)):
         return slice
     else:
         lo, hi, step = _norm_slice(slice, length)
         assert (step or 1) == 1
         # since we disallow step, a reverse slice means an empty series
         locations = [] if lo > hi else [(lo, hi)]
+        # IndelMap has this class method
+        cls = getattr(cls, "from_locations", cls)
         return cls(locations=locations, parent_length=length)
 
 
@@ -802,7 +794,7 @@ class Map:  # pragma: no cover
         """returns dicts for contained spans [dict(), ..]"""
         spans = [s.to_rich_dict() for s in self.spans]
         data = copy.deepcopy(self._serialisable)
-        data.pop("locations")
+        data.pop("locations", None)
         data["spans"] = spans
         data["type"] = get_object_provenance(self)
         data["version"] = __version__
@@ -839,7 +831,7 @@ class Map:  # pragma: no cover
 
         return zeroed
 
-    T = Union[ndarray, int]
+    T = Union[numpy.ndarray, int]
 
     def absolute_position(self, rel_pos: T) -> T:
         """converts rel_pos into an absolute position
@@ -848,7 +840,9 @@ class Map:  # pragma: no cover
         ------
         raises ValueError if rel_pos < 0
         """
-        check = array([rel_pos], dtype=int) if isinstance(rel_pos, int) else rel_pos
+        check = (
+            numpy.array([rel_pos], dtype=int) if isinstance(rel_pos, int) else rel_pos
+        )
         if check.min() < 0:
             raise ValueError(f"must positive, not {rel_pos=}")
 
@@ -865,7 +859,9 @@ class Map:  # pragma: no cover
         ------
         raises ValueError if abs_pos < 0
         """
-        check = array([abs_pos], dtype=int) if isinstance(abs_pos, int) else abs_pos
+        check = (
+            numpy.array([abs_pos], dtype=int) if isinstance(abs_pos, int) else abs_pos
+        )
         if check.min() < 0:
             raise ValueError(f"must positive, not {abs_pos=}")
         return abs_pos - self.start
@@ -1123,10 +1119,6 @@ class MapABC(ABC):
     def __add__(self, other):
         ...
 
-    @abstractmethod
-    def gaps(self):
-        ...
-
     @abstractmethod
     def nongap(self):
         ...
@@ -1143,8 +1135,18 @@ class MapABC(ABC):
     def to_rich_dict(self):
         ...
 
+    @classmethod
+    def from_locations(cls, locations, parent_length, **kwargs):
+        if len(locations):
+            spans = _spans_from_locations(locations, parent_length=parent_length)
+        else:
+            spans = ()
+
+        return cls.from_spans(spans=spans, parent_length=parent_length, **kwargs)
+
+    @classmethod
     @abstractmethod
-    def inverse(self):
+    def from_spans(cls, spans, parent_length, **kwargs):
         ...
 
 
@@ -1195,6 +1197,70 @@ def _spans_from_locations(locations, parent_length) -> Tuple[Union[Span, _LostSp
     return tuple(spans)
 
 
+O = tuple[numpy.ndarray, Sequence]
+
+
+def spans_to_gap_coords(
+    indel_spans: list[Union[Span, LostSpan]], dtype: Union[type, str] = numpy.int32
+) -> tuple[numpy.ndarray, numpy.ndarray]:
+    """returns coordinates of sequence gaps
+
+    Parameters
+    ----------
+    indel_map
+        old style Map object
+    dtype
+        string or numpy type, default to 32-bit integer
+
+    Returns
+    -------
+    numpy.array([gap pos,...]), numpy.array([cum gap length,...]),
+    """
+    gap_pos = []
+    cum_lengths = []
+    cum_length = 0
+    for i, span in enumerate(indel_spans):
+        if not span.lost:
+            continue
+
+        pos = indel_spans[i - 1].end if i else 0
+        cum_length += len(span)
+        gap_pos.append(pos)
+        cum_lengths.append(cum_length)
+
+    return numpy.array(gap_pos, dtype=dtype), numpy.array(cum_lengths, dtype=dtype)
+
+
+def _gap_spans(
+    gap_pos: NDArray[int], cum_gap_lengths: NDArray[int]
+) -> tuple[NDArray[int], NDArray[int]]:
+    """returns 1D arrays in alignment coordinates of
+    gap start, gap stop"""
+    if not len(gap_pos):
+        r = numpy.array([], dtype=gap_pos.dtype)
+        return r, r
+
+    ends = gap_pos + cum_gap_lengths
+    starts = gap_pos.copy()
+    starts[1:] += cum_gap_lengths[:-1]
+
+    return starts, ends
+
+
+def _update_lengths(result_pos, result_lengths, gap_pos, gap_lengths):
+    """modifies result_lengths in place with gap_lengths
+    where elements in gap_pos occur in result_pos
+
+    Notes
+    -----
+    result_pos is a superset of gap_pos
+    """
+    _, result_indices, other_indices = numpy.intersect1d(
+        result_pos, gap_pos, assume_unique=True, return_indices=True
+    )
+    result_lengths[result_indices] += gap_lengths[other_indices]
+
+
 T = Union[List[int], Tuple[int]]
 
 
@@ -1202,146 +1268,496 @@ T = Union[List[int], Tuple[int]]
 class IndelMap(MapABC):
     """store locations of deletions in a Aligned sequence"""
 
-    spans: dataclasses.InitVar[Optional[tuple]] = ()
-    locations: dataclasses.InitVar[Optional[Sequence[T]]] = None
-    termini_unknown: dataclasses.InitVar[bool] = dataclasses.field(default=False)
+    # gap data is gap positions, gap lengths on input, stored
+    gap_pos: numpy.ndarray
+    cum_gap_lengths: Optional[numpy.ndarray] = None
+    gap_lengths: dataclasses.InitVar[Optional[numpy.ndarray]] = None
+    termini_unknown: bool = False
     parent_length: int = 0
-    start: Optional[int] = dataclasses.field(init=False, default=0)
-    end: Optional[int] = dataclasses.field(init=False, default=0)
-    length: int = dataclasses.field(init=False, default=0)
     _serialisable: dict = dataclasses.field(init=False, repr=False)
-    _spans: Tuple[Union[Span, _LostSpan, TerminalPadding]] = dataclasses.field(
-        default=(), init=False
-    )
+    num_gaps: int = dataclasses.field(init=False, repr=False, default=0)
 
-    def __post_init__(self, spans, locations, termini_unknown):
-        if locations is not None and len(locations):
-            spans = _spans_from_locations(locations, parent_length=self.parent_length)
-        elif isinstance(spans, property):
-            # This clause is due a known issue with dataclasses.
-            # As we have a spans property, the default spans value is
-            # ignored, so we have to check for its value being a property
-            # and then set the default value here
-            spans = ()
+    def __post_init__(self, gap_lengths):
+        assert gap_lengths is None or self.cum_gap_lengths is None
+        if gap_lengths is not None:
+            self.cum_gap_lengths = gap_lengths.cumsum()
 
-        spans = tuple(spans)
-        last_not_lost = None
-        start, end = None, None
-        for i, span in enumerate(spans):
-            self.length += len(span)
-            if span.lost:
-                continue
-            elif start is None:
-                # this ugly logic because we're using spans!
-                start = span.end if span.reverse else span.start
+        if len(self.gap_pos) != len(self.cum_gap_lengths):
+            raise ValueError(
+                f"length of gap_ pos {len(self.gap_pos)} != "
+                f"length of gap lengths {len(self.cum_gap_lengths)}"
+            )
 
-            last_not_lost = i
+        self.num_gaps = self.gap_pos.shape[0]
+        if self.num_gaps and self.gap_pos[-1] > self.parent_length:
+            raise ValueError(
+                f"gap position {self.gap_pos[-1]} outside parent_length {self.parent_length}"
+            )
 
-        if last_not_lost is not None:
-            span = spans[last_not_lost]
-            end = span.start if span.reverse else span.end
+        # make gap array immutable
+        self.gap_pos.flags.writeable = False
+        self.cum_gap_lengths.flags.writeable = False
+        self._serialisable.pop("gap_lengths", None)
 
-        if start is None:
-            start = 0
+    @classmethod
+    def from_spans(cls, spans, parent_length, termini_unknown=False):
+        gap_pos, cum_lengths = spans_to_gap_coords(spans)
+        return cls(
+            gap_pos=gap_pos,
+            cum_gap_lengths=cum_lengths,
+            parent_length=parent_length,
+            termini_unknown=termini_unknown,
+        )
 
-        if end is None:
-            end = self.parent_length
+    @classmethod
+    def from_aligned_segments(
+        cls, locations: list[tuple[int, int]], aligned_length: int
+    ):
+        """
+        converts coordinates from aligned segments into IndelMap for ungapped sequence
+
+        Parameters
+        ----------
+        locations
+            list of ungapped segment in alignment coordinates
+        aligned_length
+            length of the alignment
+        """
+        if not locations or (
+            len(locations) == 1
+            and locations[0][0] == 0
+            and locations[0][-1] == aligned_length
+        ):
+            empty = numpy.array([], dtype=int)
+            return cls(
+                gap_pos=empty,
+                cum_gap_lengths=empty.copy(),
+                parent_length=aligned_length,
+            )
 
-        self.start, self.end = start, end
+        if locations[0][0] != 0:
+            # starts with a gap
+            locations = [(0, 0)] + locations
+        if locations[-1][1] < aligned_length:
+            # ends with a gap
+            locations += [(aligned_length, aligned_length)]
+
+        locations = numpy.array(locations, dtype=numpy.int32).flatten()[1:-1]
+        gap_coords = locations.reshape((locations.shape[0] // 2, 2))
+        gap_ends, gap_starts = gap_coords[:, ::-1].T
+        gap_lengths = gap_ends - gap_starts
+        cum_lens = gap_lengths.cumsum()
+        # convert to sequence coords
+        gap_pos = gap_starts.copy()
+        gap_pos[1:] -= cum_lens[:-1]
+        seq_length = aligned_length - cum_lens[-1]
+
+        return cls(
+            gap_pos=gap_pos,
+            cum_gap_lengths=cum_lens,
+            parent_length=seq_length,
+        )
 
-        if termini_unknown:
-            spans = list(spans)
-            if spans[0].lost:
-                spans[0] = TerminalPadding(spans[0].length)
-            if spans[-1].lost:
-                spans[-1] = TerminalPadding(spans[-1].length)
+    @functools.singledispatchmethod
+    def __getitem__(self, item):
+        raise NotImplementedError(f"cannot slice using {type(item)}")
+
+    @__getitem__.register
+    def _(self, item: int):
+        return self[item : item + 1]
+
+    @__getitem__.register
+    def _(self, item: slice):
+        # we're assuming that this gap object is associated with a sequence
+        # that will also be sliced. Hence, we need to shift the gap insertion
+        # positions relative to this newly sliced sequence.
+        if item.step is not None:
+            raise NotImplementedError(
+                f"{type(self).__name__!r} does not yet support strides"
+            )
+        zero_array = numpy.array([], dtype=self.gap_pos.dtype)
+        start = item.start or 0
+        stop = item.stop or len(self)
+
+        # convert negative indices
+        start = start if start >= 0 else len(self) + start
+        stop = stop if stop >= 0 else len(self) + stop
+        if min((start, stop)) < 0:
+            raise IndexError(f"one of adjusted {start, stop} is < 0")
+
+        if start >= stop:
+            # standard slice behaviour without negative step
+            return self.__class__(
+                gap_pos=zero_array.copy(),
+                cum_gap_lengths=zero_array.copy(),
+                parent_length=0,
+            )
 
-        self._spans = tuple(spans)
+        # we address three easy cases:
+        # 1 - no gaps; 2 - slice before first gap; 3 - after last gap
+        no_gaps = self.__class__(
+            gap_pos=zero_array.copy(),
+            cum_gap_lengths=zero_array.copy(),
+            parent_length=stop - start,
+        )
+        if not self.num_gaps:
+            return no_gaps
+
+        first_gap = self.gap_pos[0]
+        last_gap = self.gap_pos[-1] + self.cum_gap_lengths[-1]
+        if stop < first_gap or start >= last_gap:
+            return no_gaps
+
+        gap_starts, gap_ends = _gap_spans(self.gap_pos, self.cum_gap_lengths)
+        gap_pos = self.gap_pos.copy()
+        cum_lengths = self.cum_gap_lengths.copy()
+        # we find where the slice starts
+        l = numpy.searchsorted(gap_ends, start, side="left")
+        if gap_starts[l] <= start < gap_ends[l] and stop <= gap_ends[l]:
+            # entire span is within a single gap
+            # pos now 0
+            gap_pos = numpy.array([0], dtype=self.gap_pos.dtype)
+            cum_lengths = cum_lengths[l : l + 1]
+            cum_lengths[0] = stop - start
+            return self.__class__(
+                gap_pos=gap_pos, cum_gap_lengths=cum_lengths, parent_length=0
+            )
 
-    def __getitem__(self, slice):
-        # A possible shorter map at the same level
-        new_map = as_map(slice, len(self), self.__class__)
-        new_parts = []
-        for span in new_map.spans:
-            # we reset tidy start / end to false to avoid changes to
-            # default behaviour of tidy on IndelMap (which ignores these
-            # anyway)
-            old_tidy = span.tidy_start, span.tidy_end
-            span.tidy_start, span.tidy_end = False, False
-            new_parts.extend(span.remap_with(self))
-            span.tidy_start, span.tidy_end = old_tidy
-        return self.__class__(spans=new_parts, parent_length=self.parent_length)
+        lengths = self.get_gap_lengths()
+        if start < first_gap:
+            # start is before the first gap, we don't slice or shift
+            shift = start
+            begin = 0
+        elif gap_starts[l] <= start < gap_ends[l]:
+            # start is within a gap
+            # so the absolute gap_pos value remains unchanged, but we shorten
+            # the gap length
+            begin = l
+            begin_diff = start - gap_starts[l]  # if l else self.gap_pos[l]
+            lengths[l] -= begin_diff
+            shift = (start - cum_lengths[l - 1] - begin_diff) if l else gap_pos[0]
+        elif start == gap_ends[l]:
+            # at gap boundary, so beginning of non-gapped segment
+            # no adjustment to gap lengths
+            begin = l + 1
+            shift = start - cum_lengths[l]
+        else:
+            # not within a gap
+            begin = l
+            shift = start - cum_lengths[l - 1] if l else start
+
+        # start search for stop from l index
+        r = numpy.searchsorted(gap_ends[l:], stop, side="right") + l
+        if r == self.num_gaps:
+            # stop is after last gap
+            end = r
+        elif gap_starts[r] < stop <= gap_ends[r]:
+            # within gap
+            end = r + 1
+            end_diff = gap_ends[r] - stop
+            lengths[r] -= end_diff
+        else:
+            end = r
+
+        pos_result = gap_pos[begin:end]
+        pos_result -= shift
+        lengths = lengths[begin:end]
+        parent_length = self.get_seq_index(stop) - self.get_seq_index(start)
+
+        return self.__class__(
+            gap_pos=pos_result,
+            gap_lengths=lengths,
+            parent_length=parent_length,
+        )
+
+    def get_align_index(self, seq_index: int, slice_stop: bool = False) -> int:
+        """convert a sequence index into an alignment index
+
+        Parameters
+        ----------
+        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
+        if seq_index < 0:
+            seq_index += self.parent_length
+
+        if seq_index < 0:
+            raise NotImplementedError(f"{seq_index} negative seq_index beyond limit ")
+
+        if not self.num_gaps or seq_index < self.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 == self.gap_pos).any():
+            # if so, we return the alignment coord for the first gap position
+            (idx,) = numpy.where(match)[0]
+            if idx:
+                gap_len = self.cum_gap_lengths[idx] - self.cum_gap_lengths[idx - 1]
+            else:
+                gap_len = self.cum_gap_lengths[idx]
+            gap_end = self.gap_pos[idx] + self.cum_gap_lengths[idx]
+            return int(gap_end - gap_len)
+
+        cum_gap_lengths = self.cum_gap_lengths
+        gap_pos = self.gap_pos
+
+        if seq_index >= gap_pos[-1]:
+            return int(seq_index + cum_gap_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_gap_lengths[index - 1] if index else 0
+        else:
+            gap_lengths = cum_gap_lengths[index]
+
+        return int(seq_index + gap_lengths)
+
+    def get_seq_index(self, align_index: int) -> int:
+        """converts alignment index to sequence 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 NotImplementedError(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])
 
     def __len__(self):
-        return self.length
+        length_gaps = self.cum_gap_lengths[-1] if self.num_gaps else 0
+        return int(self.parent_length + length_gaps)
 
     def __add__(self, other):
-        if other.parent_length != self.parent_length:
-            raise ValueError("Those maps belong to different sequences")
+        # what was the purpose of this method? The code seems designed for
+        # combining Maps from the same parent sequence, which is a union rather
+        # than addition
+        # I'm revising this to be suitable for the concatenation of two aligned
+        # sequences
+        gap_pos = self.gap_pos.tolist() + (self.parent_length + other.gap_pos).tolist()
+        gap_pos = numpy.array(gap_pos, dtype=self.gap_pos.dtype)
+
+        cum_length = self.cum_gap_lengths[-1] if self.num_gaps else 0
+        cum_gap_lengths = (
+            self.cum_gap_lengths.tolist()
+            + (cum_length + other.cum_gap_lengths).tolist()
+        )
+        cum_gap_lengths = numpy.array(cum_gap_lengths, dtype=self.cum_gap_lengths.dtype)
+
         return self.__class__(
-            spans=self._spans + tuple(other.spans), parent_length=self.parent_length
+            gap_pos=gap_pos,
+            cum_gap_lengths=cum_gap_lengths,
+            parent_length=self.parent_length + other.parent_length,
         )
 
     def __mul__(self, scale):
-        new_parts = [span * scale for span in self._spans]
-        return self.__class__(spans=new_parts, parent_length=self.parent_length * scale)
+        # could be used for going from amino-acid alignment to codon alignment
+        gap_pos = self.gap_pos * scale
+        cum_gap_lengths = self.cum_gap_lengths * scale
+        return self.__class__(
+            gap_pos=gap_pos,
+            cum_gap_lengths=cum_gap_lengths,
+            parent_length=self.parent_length * scale,
+        )
 
     def __repr__(self):
-        return f"{self._spans!r}/{self.parent_length}"
+        gap_data = numpy.array([self.gap_pos, self.cum_gap_lengths]).T
+        return f"{gap_data.tolist()!r}/{self.parent_length}"
+
+    def get_gap_lengths(self) -> NDArray[int]:
+        lengths = self.cum_gap_lengths.copy()
+        lengths[1:] = numpy.diff(lengths)
+        return lengths
 
     @property
     def offsets(self):
-        return [0] + array([s.length for s in self._spans[:-1]]).cumsum().tolist()
+        # offsets are the aligned indices for every starting point of a segment
+        # when we encounter a gap, we include that position and the end of that gap
+
+        starts, ends = _gap_spans(self.gap_pos, self.cum_gap_lengths)
+        return numpy.array([starts, ends]).T.flatten()[:-1].tolist()
+
+    def nongap(self) -> Iterator[Union[Span, LostSpan]]:
+        """ungappeed segments in this map in aligned coordinates"""
+        # we want to know the coordinates of the ungapped segments on
+        # the aligned sequence. The gap_pos attribute is in sequence
+        # coordinates
+        prev_pos = 0
+        for i, pos in enumerate(self.gap_pos):
+            if pos == 0:
+                # we start with a gap
+                prev_pos = pos
+                continue
 
-    def gaps(self):
-        """The gaps (lost spans) in this map"""
-        locations = []
-        offset = 0
-        for s in self._spans:
-            if s.lost:
-                locations.append((offset, offset + s.length))
-            offset += s.length
-        return self.__class__(locations=locations, parent_length=len(self))
+            cum_length = 0 if i == 0 else self.cum_gap_lengths[i - 1]
+            start = 0 if i == 0 else prev_pos
+            start += cum_length
+            end = self.gap_pos[i] + cum_length
+            yield Span(start, end)
+            prev_pos = pos
 
-    def nongap(self):
-        """ungappeed segments in this map"""
-        locations = []
-        offset = 0
-        for s in self._spans:
-            if not s.lost:
-                locations.append((offset, offset + s.length))
-            offset += s.length
-        return self.__class__(locations=locations, parent_length=len(self))
+        if self.num_gaps and self.gap_pos[-1] + self.cum_gap_lengths[-1] < len(self):
+            yield Span(prev_pos, len(self))
 
     @property
-    def spans(self) -> Iterator[Span]:
+    def spans(self) -> Iterator[Union[Span, LostSpan]]:
         """generator of spans"""
-        yield from self._spans
+        if not self.num_gaps:
+            yield Span(0, self.parent_length)
+            return
+
+        for i, pos in enumerate(self.gap_pos):
+            cum_length = self.cum_gap_lengths[i]
+            if pos == 0:
+                cls = TerminalPadding if self.termini_unknown else LostSpan
+                yield cls(cum_length)
+                continue
+
+            if i == 0:
+                start = 0
+                prev_length = 0
+            else:
+                start = self.gap_pos[i - 1]
+                prev_length = self.cum_gap_lengths[i - 1]
+            yield Span(start, pos)
+            cls = (
+                TerminalPadding
+                if self.termini_unknown and i == self.num_gaps - 1
+                else LostSpan
+            )
+            yield cls(cum_length - prev_length)
+
+        if self.num_gaps and self.gap_pos[-1] < self.parent_length:
+            yield Span(self.gap_pos[-1], self.parent_length)
 
     @property
     def complete(self):
         """whether any span represents a gap"""
-        return not any(span.lost for span in self._spans)
+        return self.num_gaps != 0 and self.useful
 
     @property
     def useful(self):
-        return not all(span.lost for span in self._spans)
+        return self.parent_length != 0
 
     def get_coordinates(self):
-        """returns span coordinates as [(start, end), ...]"""
-        return [(s.start, s.end) for s in self._spans if not s.lost]
+        """returns sequence coordinates of ungapped segments
+
+        Returns
+        -------
+        [(start, end), ...]
+        """
+        coords = []
+        last = 0
+        for pos, cum_length in zip(self.gap_pos, self.cum_gap_lengths):
+            pos, cum_length = int(pos), int(cum_length)
+
+            if not pos:
+                continue
+            coords.append((last, pos))
+            last = pos
+
+        if (
+            self.num_gaps
+            and self.gap_pos[-1] + self.cum_gap_lengths[-1] < self.parent_length
+        ):
+            coords.append((last, self.parent_length))
+        return coords
 
     def get_gap_coordinates(self):
         """returns [(gap pos, gap length), ...]"""
-        gap_pos = []
-        for i, span in enumerate(self._spans):
-            if not span.lost:
-                continue
+        cum_lengths = self.cum_gap_lengths.copy()
+        diffs = numpy.diff(cum_lengths)
+        cum_lengths[1:] = diffs
+        return numpy.array([self.gap_pos, cum_lengths]).T.tolist()
+
+    def merge_maps(self, other, parent_length: Optional[int] = None):
+        """merge gaps of other with self
+
+        Parameters
+        ----------
+        indel_map
+            instance for same sequence
+        parent_length
+            overrides property
+        """
+        unique_pos = numpy.union1d(self.gap_pos, other.gap_pos)
+        gap_lengths = numpy.zeros(unique_pos.shape, dtype=self.cum_gap_lengths.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)
+        parent_length = parent_length or self.parent_length
+        return self.__class__(
+            gap_pos=unique_pos,
+            gap_lengths=gap_lengths,
+            parent_length=parent_length,
+        )
 
-            pos = self._spans[i - 1].end if i else 0
-            gap_pos.append((pos, len(span)))
+    def joined_segments(self, coords: list[tuple[int, int]]):
+        """returns new map with disjoint gapped segments joined
 
-        return gap_pos
+        Parameters
+        ----------
+        coords
+            sequence insert gap coordinates [(gap start, gap end), ...]
+        Returns
+        -------
+
+        """
+        coords = list(sorted(coords))
+        # using a dict here because joining can produce a gap merge
+        gaps = {}
+        cum_length = 0
+        cum_parent_length = 0
+        dtype = self.gap_pos.dtype
+        for start, end in coords:
+            im = self[start:end]
+            for i in range(im.num_gaps):
+                pos = im.gap_pos[i] + cum_parent_length
+                if pos in gaps:
+                    gaps[pos] += im.cum_gap_lengths[i]
+                else:
+                    gaps[pos] = im.cum_gap_lengths[i] + cum_length
+            cum_parent_length += im.parent_length
+            if im.num_gaps:
+                cum_length += im.cum_gap_lengths[-1]
+        gap_pos = numpy.empty(len(gaps), dtype=dtype)
+        cum_lengths = numpy.empty(len(gaps), dtype=dtype)
+        for i, (pos, length) in enumerate(sorted(gaps.items())):
+            gap_pos[i] = pos
+            cum_lengths[i] = length
+        return self.__class__(
+            gap_pos=gap_pos,
+            cum_gap_lengths=cum_lengths,
+            parent_length=cum_parent_length,
+        )
 
     def nucleic_reversed(self):
         """map for a sequence that has itself been reversed and complemented
@@ -1350,32 +1766,34 @@ class IndelMap(MapABC):
         -----
         discards reverse attribute on both spans and self
         """
-        spans = []
-        parent_length = self.parent_length
-        for s in self._spans:
-            if not s.lost:
-                start = parent_length - s.end
-                assert start >= 0
-                end = start + s.length
-                s = Span(start=start, end=end)
-            spans.append(s)
+        new_pos = self.gap_pos.copy()
+        lengths = self.get_gap_lengths()
 
-        spans.reverse()
-        return self.__class__(spans=spans, parent_length=self.parent_length)
+        if len(new_pos):
+            new_pos = self.parent_length - new_pos
+            new_pos = new_pos[::-1]
+            lengths = lengths[::-1]
+
+        return self.__class__(
+            gap_pos=new_pos,
+            gap_lengths=lengths,
+            parent_length=self.parent_length,
+        )
 
     def to_rich_dict(self):
         """returns dicts for contained spans [dict(), ..]"""
-        spans = [s.to_rich_dict() for s in self._spans]
         # exclude spans from deep copy since being overwritten
-        data = copy.deepcopy(
-            {k: v for k, v in self._serialisable.items() if k != "spans"}
-        )
-        data.pop("locations")
-        data["spans"] = spans
+        data = copy.deepcopy(dict(self._serialisable.items()))
         data["type"] = get_object_provenance(self)
         data["version"] = __version__
+        data["gap_pos"] = self.gap_pos.tolist()
+        data["cum_gap_lengths"] = self.cum_gap_lengths.tolist()
+        data["parent_length"] = int(self.parent_length)
         return data
 
+    def to_json(self):
+        return json.dumps(self.to_rich_dict())
+
     @classmethod
     def from_rich_dict(cls, map_element):
         from cogent3.util.deserialise import _get_class
@@ -1383,158 +1801,53 @@ class IndelMap(MapABC):
         map_element.pop("version", None)
         type_ = map_element.pop("type")
         assert _get_class(type_) == cls
-        spans = []
-        for element in map_element["spans"]:
-            element.pop("version", None)
-            klass = _get_class(element.pop("type"))
-            instance = klass(**element)
-            spans.append(instance)
+        map_element["gap_pos"] = numpy.array(map_element["gap_pos"])
+        map_element["cum_gap_lengths"] = numpy.array(map_element["cum_gap_lengths"])
 
-        map_element["spans"] = spans
         return cls(**map_element)
 
     def with_termini_unknown(self):
         """returns new instance with terminal gaps indicated as unknown"""
         return self.__class__(
-            spans=self._spans[:],
+            gap_pos=self.gap_pos.copy(),
+            cum_gap_lengths=self.cum_gap_lengths.copy(),
             parent_length=self.parent_length,
             termini_unknown=True,
         )
 
-    def inverse(self):
-        """returns instance with coordinates updated for aligned, unaligned"""
-        # is this only required for parse_out_gaps?
-        # NO also used in cogent3.align code
-
-        # can't work if there are overlaps in the map
-        # tidy ends don't survive inversion
-        if self.parent_length is None:
-            raise ValueError("Uninvertable. parent length not known")
-
-        cum_posn = 0
-        temp = []
-        for span in self._spans:
-            if not span.lost:
-                if span.reverse:
-                    temp.append(
-                        (span.start, span.end, cum_posn + span.length, cum_posn)
-                    )
-                else:
-                    temp.append(
-                        (span.start, span.end, cum_posn, cum_posn + span.length)
-                    )
-            cum_posn += span.length
-
-        temp.sort()
-        new_spans = []
-        last_start = 0
-        for start, end, cum_start, cum_end in temp:
-            if start > last_start:
-                new_spans.append(LostSpan(start - last_start))
-            elif start < last_start:
-                raise ValueError(f"Uninvertable. Overlap: {start} < {last_start}")
-
-            # we force tidy_<start/end> to be same as self, attribute has no meaning
-            # for IndelMap, but retained for compatability for now
-            new_spans.append(
-                Span(
-                    cum_start,
-                    cum_end,
-                    reverse=cum_start > cum_end,
-                )
-            )
-            last_start = end
-
-        if self.parent_length > last_start:
-            new_spans.append(LostSpan(self.parent_length - last_start))
-
-        return self.__class__(spans=new_spans, parent_length=len(self))
-
-    T = Union[ndarray, int]
-
-    def absolute_position(self, rel_pos: T) -> T:
-        """converts rel_pos into an absolute position
-
-        Raises
-        ------
-        raises ValueError if rel_pos < 0
-        """
-        check = array([rel_pos], dtype=int) if isinstance(rel_pos, int) else rel_pos
-        if check.min() < 0:
-            raise ValueError(f"must positive, not {rel_pos=}")
-
-        if len(self) == self.parent_length:
-            # handle case of reversed here?
-            return rel_pos
-
-        return self.start + rel_pos
+    def to_feature_map(self):
+        """returns a Map type, suited to Features"""
+        return FeatureMap(spans=list(self.spans), parent_length=self.parent_length)
 
-    def relative_position(self, abs_pos: T) -> T:
-        """converts abs_pos into an relative position
+    def make_seq_feature_map(self, align_feature_map, include_gaps=True):
+        """converts align_feature_map to a FeatureMap with sequence coordinates
 
-        Raises
-        ------
-        raises ValueError if abs_pos < 0
+        Parameters
+        ----------
+        align_feature_map
+            with alignment coordinates
+        include_gaps
+            whether to include gaps from self as LostSpan's
         """
-        check = array([abs_pos], dtype=int) if isinstance(abs_pos, int) else abs_pos
-        if check.min() < 0:
-            raise ValueError(f"must positive, not {abs_pos=}")
-        return abs_pos - self.start
-
-    def get_covering_span(self):
-        span = self.start, self.end
-        return self.__class__(locations=[span], parent_length=self.parent_length)
-
-    def zeroed(self):
-        """returns a new instance with the first span starting at 0
-
-        Note
-        ----
-
-        Useful when an annotated sequence is sliced, but the connection to
-        the original parent is being deliberately broken as in the
-        Sequence.deepcopy(sliced=True) case.
-        """
-        # todo is this really required, best ifn we can rely on SeqView to
-        #  store all relationship to underlying sequence
-        min_val = min(self.start, self.end)
+        gap_spans = {}
+        if include_gaps and self.num_gaps:
+            last = 0
+            for pos, cum_length in zip(self.gap_pos, self.cum_gap_lengths):
+                gap_spans[pos] = LostSpan(cum_length - last)
+                last = cum_length
         spans = []
-        for span in self._spans:
+        for span in align_feature_map.spans:
             if span.lost:
                 spans.append(span)
                 continue
-            kwargs = span.to_rich_dict()
-            del kwargs["version"]
-            del kwargs["type"]
-            kwargs["start"] = kwargs["start"] - min_val
-            kwargs["end"] = kwargs["end"] - min_val
-            spans.append(Span(**kwargs))
-
-        kwargs = self.to_rich_dict()
-        del kwargs["version"]
-        del kwargs["type"]
-        kwargs["spans"] = spans
-        kwargs["parent_length"] = abs(self.start - self.end)
-        return self.__class__(**kwargs)
 
-    def to_feature_map(self):
-        """returns a Map type, suited to Features"""
-        spans = []
-        for span in self._spans:
-            if span.lost:
-                spans.append(span)
-                continue
-            kwargs = span.to_rich_dict()
-            del kwargs["version"]
-            del kwargs["type"]
-            spans.append(Span(**kwargs))
+            start = self.get_seq_index(span.start)
+            end = self.get_seq_index(span.end)
+            if lost := gap_spans.pop(span.start, None):
+                spans.append(lost)
+            spans.append(Span(start, end))
 
-        kwargs = self.to_rich_dict()
-        del kwargs["version"]
-        del kwargs["type"]
-        del kwargs["termini_unknown"]  # not used for FeatureMap
-        kwargs["spans"] = spans
-        return FeatureMap(**kwargs)
+        return FeatureMap(spans=spans, parent_length=align_feature_map.parent_length)
 
 
 @dataclasses.dataclass
@@ -1542,7 +1855,6 @@ class FeatureMap(MapABC):
     """A map holds a list of spans."""
 
     spans: dataclasses.InitVar[Optional[tuple]] = ()
-    locations: dataclasses.InitVar[Optional[Sequence[T]]] = None
     parent_length: int = 0
     offsets: list[int] = dataclasses.field(init=False, repr=False)
     useful: bool = dataclasses.field(init=False, repr=False, default=False)
@@ -1552,12 +1864,10 @@ class FeatureMap(MapABC):
         default=(), init=False
     )
 
-    def __post_init__(self, spans, locations):
+    def __post_init__(self, spans):
         assert self.parent_length is not None
-
-        if locations is not None and len(locations):
-            spans = _spans_from_locations(locations, parent_length=self.parent_length)
-        elif isinstance(spans, property):
+        self.parent_length = int(self.parent_length)
+        if isinstance(spans, property):
             # This clause is due a known issue with dataclasses.
             # As we have a spans property, the default spans value is
             # ignored, so we have to check for its value being a property
@@ -1585,6 +1895,10 @@ class FeatureMap(MapABC):
         self._spans = tuple(spans)
         self.length = posn
 
+    @classmethod
+    def from_spans(cls, spans, parent_length):
+        return cls(spans=spans, parent_length=parent_length)
+
     def __len__(self):
         return self.length
 
@@ -1620,17 +1934,11 @@ class FeatureMap(MapABC):
     def spans(self):
         yield from self._spans
 
-    def _with_termini_unknown(self):
-        return self.__class__(
-            self,
-            spans=self.spans[:],
-            parent_length=self.parent_length,
-            termini_unknown=True,
-        )
-
     def get_covering_span(self):
         span = (self.start, self.end)
-        return self.__class__(locations=[span], parent_length=self.parent_length)
+        return self.__class__.from_locations(
+            locations=[span], parent_length=self.parent_length
+        )
 
     def covered(self):
         """>>> Map([(10,20), (15, 25), (80, 90)]).covered().spans
@@ -1659,7 +1967,9 @@ class FeatureMap(MapABC):
             last_x = x
             last_y = y
         assert y == 0
-        return self.__class__(locations=result, parent_length=self.parent_length)
+        return self.__class__.from_locations(
+            locations=result, parent_length=self.parent_length
+        )
 
     def nucleic_reversed(self):
         """map for a sequence that has itself been reversed and complemented
@@ -1702,7 +2012,9 @@ class FeatureMap(MapABC):
             if s.lost:
                 locations.append((offset, offset + s.length))
             offset += s.length
-        return self.__class__(locations=locations, parent_length=len(self))
+        return self.__class__.from_locations(
+            locations=locations, parent_length=len(self)
+        )
 
     def shadow(self):
         """The 'negative' map of the spans not included in this map"""
@@ -1715,7 +2027,7 @@ class FeatureMap(MapABC):
             if not s.lost:
                 locations.append((offset, offset + s.length))
             offset += s.length
-        return self.__class__(locations=locations, parent_length=len(self))
+        return _spans_from_locations(locations=locations, parent_length=len(self))
 
     def without_gaps(self):
         return self.__class__(
@@ -1772,35 +2084,6 @@ class FeatureMap(MapABC):
 
         return self.__class__(spans=new_spans, parent_length=len(self))
 
-    def _inverse(self):
-        # can't work if there are overlaps in the map
-        # tidy ends don't survive inversion
-        if self.parent_length is None:
-            raise ValueError("Uninvertable. parent length not known")
-        posn = 0
-        temp = []
-        for span in self.spans:
-            if not span.lost:
-                if span.reverse:
-                    temp.append((span.start, span.end, posn + span.length, posn))
-                else:
-                    temp.append((span.start, span.end, posn, posn + span.length))
-            posn += span.length
-
-        temp.sort()
-        new_spans = []
-        last_hi = 0
-        for lo, hi, start, end in temp:
-            if lo > last_hi:
-                new_spans.append(LostSpan(lo - last_hi))
-            elif lo < last_hi:
-                raise ValueError(f"Uninvertable. Overlap: {lo} < {last_hi}")
-            new_spans.append(Span(start, end, reverse=start > end))
-            last_hi = hi
-        if self.parent_length > last_hi:
-            new_spans.append(LostSpan(self.parent_length - last_hi))
-        return self.__class__(spans=new_spans, parent_length=len(self))
-
     def get_coordinates(self):
         """returns span coordinates as [(v1, v2), ...]
 
@@ -1813,12 +2096,33 @@ class FeatureMap(MapABC):
         """returns dicts for contained spans [dict(), ..]"""
         spans = [s.to_rich_dict() for s in self.spans]
         data = copy.deepcopy(self._serialisable)
-        data.pop("locations")
+        data.pop("locations", None)
         data["spans"] = spans
         data["type"] = get_object_provenance(self)
         data["version"] = __version__
+        data["parent_length"] = int(self.parent_length)
         return data
 
+    def to_json(self):
+        return json.dumps(self.to_rich_dict())
+
+    @classmethod
+    def from_rich_dict(cls, map_element):
+        from cogent3.util.deserialise import _get_class
+
+        map_element.pop("version", None)
+        type_ = map_element.pop("type")
+        assert _get_class(type_) == cls
+        spans = []
+        for element in map_element.pop("spans"):
+            element.pop("version", None)
+            klass = _get_class(element.pop("type"))
+            instance = klass(**element)
+            spans.append(instance)
+
+        map_element["spans"] = spans
+        return cls(**map_element)
+
     def zeroed(self):
         """returns a new instance with the first span starting at 0
 
@@ -1850,7 +2154,7 @@ class FeatureMap(MapABC):
 
         return zeroed
 
-    T = Union[ndarray, int]
+    T = Union[numpy.ndarray, int]
 
     def absolute_position(self, rel_pos: T) -> T:
         """converts rel_pos into an absolute position
@@ -1859,7 +2163,9 @@ class FeatureMap(MapABC):
         ------
         raises ValueError if rel_pos < 0
         """
-        check = array([rel_pos], dtype=int) if isinstance(rel_pos, int) else rel_pos
+        check = (
+            numpy.array([rel_pos], dtype=int) if isinstance(rel_pos, int) else rel_pos
+        )
         if check.min() < 0:
             raise ValueError(f"must positive, not {rel_pos=}")
 
@@ -1876,13 +2182,15 @@ class FeatureMap(MapABC):
         ------
         raises ValueError if abs_pos < 0
         """
-        check = array([abs_pos], dtype=int) if isinstance(abs_pos, int) else abs_pos
+        check = (
+            numpy.array([abs_pos], dtype=int) if isinstance(abs_pos, int) else abs_pos
+        )
         if check.min() < 0:
             raise ValueError(f"must positive, not {abs_pos=}")
         return abs_pos - self.start
 
 
-def gap_coords_to_map(gaps_lengths: dict, seq_length: int) -> FeatureMap:
+def gap_coords_to_map(gaps_lengths: dict, seq_length: int) -> IndelMap:
     """
     Parameters
     ----------
@@ -1897,21 +2205,53 @@ def gap_coords_to_map(gaps_lengths: dict, seq_length: int) -> FeatureMap:
     """
 
     if not gaps_lengths:
-        return FeatureMap(locations=[(0, seq_length)], parent_length=seq_length)
+        gap_pos = numpy.array([], dtype=int)
+        lengths = gap_pos.copy()
+    else:
+        gap_pos, lengths = list(zip(*sorted(gaps_lengths.items())))
+        gap_pos = numpy.array(gap_pos, dtype=int)
+        lengths = numpy.array(lengths, dtype=int)
+
+    return IndelMap(gap_pos=gap_pos, gap_lengths=lengths, parent_length=seq_length)
+
+
+@register_deserialiser(get_object_provenance(IndelMap))
+def deserialise_indelmap(data: dict) -> IndelMap:
+    return IndelMap.from_rich_dict(data)
+
+
+@register_deserialiser(get_object_provenance(FeatureMap))
+def deserialise_featuremap(data: dict) -> FeatureMap:
+    return FeatureMap.from_rich_dict(data)
+
+
+@functools.singledispatch
+def _norm_slice(index, length):
+    """_norm_slice(slice(1, -2, 3), 10) -> (1,8,3)"""
+    start = index
+    if start < 0:
+        start += length
+    if start >= length:
+        raise IndexError(index)
+    return start, start + 1, 1
+
+
+@_norm_slice.register
+def _(index: slice, length):
+    start = _norm_index(index.start, length, 0)
+    end = _norm_index(index.stop, length, length)
+    return start, end, index.step
 
-    spans = []
-    last = pos = 0
-    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=gaps_lengths[pos])
-        spans.extend([gap] if pos == 0 else [Span(last, pos), gap])
-        last = pos
+@_norm_slice.register
+def _(index: Span, length):
+    start = _norm_index(index.start, length, 0)
+    end = _norm_index(index.end, length, length)
+    return start, end, None
 
-    if pos < seq_length:
-        spans.append(Span(last, seq_length))
 
-    return FeatureMap(spans=spans, parent_length=seq_length)
+@_norm_slice.register
+def _(index: FeatureMap, length):
+    start = _norm_index(index.start, length, 0)
+    end = _norm_index(index.end, length, length)
+    return start, end, None
diff --git a/src/cogent3/core/sequence.py b/src/cogent3/core/sequence.py
index 022fd5a57..2b28e59e7 100644
--- a/src/cogent3/core/sequence.py
+++ b/src/cogent3/core/sequence.py
@@ -143,7 +143,7 @@ class SequenceI(object):
             version=__version__,
         )
         if hasattr(self, "annotation_offset"):
-            offset = self._seq.parent_start
+            offset = int(self._seq.parent_start)
             data.update(dict(annotation_offset=offset))
 
         if (
@@ -1060,7 +1060,7 @@ class Sequence(SequenceI):
                 continue
             new_spans.append(new.tolist())
 
-        fmap = FeatureMap(locations=new_spans, parent_length=len(self))
+        fmap = FeatureMap.from_locations(locations=new_spans, parent_length=len(self))
         if pre or post:
             # create a lost span to represent the segment missing from
             # the instance
@@ -1256,7 +1256,7 @@ class Sequence(SequenceI):
                 seqid=self.name,
                 name=None,
                 biotype=None,
-                map=FeatureMap(locations=[], parent_length=len(self)),
+                map=FeatureMap.from_locations(locations=[], parent_length=len(self)),
                 strand="+",
             )
         else:
@@ -1459,18 +1459,27 @@ class Sequence(SequenceI):
 
     def parse_out_gaps(self):
         """returns Map corresponding to gap locations and ungapped Sequence"""
-        gapless = []
-        segments = []
-        nongap = re.compile(f"([^{re.escape('-')}]+)")
-        for match in nongap.finditer(str(self)):
-            segments.append(match.span())
-            gapless.append(match.group())
-        map = IndelMap(locations=segments, parent_length=len(self)).inverse()
+        gap = re.compile(f"[{re.escape(self.moltype.gap)}]+")
+        seq = str(self)
+        gap_pos = []
+        cum_lengths = []
+        for match in gap.finditer(seq):
+            pos = match.start()
+            gap_pos.append(pos)
+            cum_lengths.append(match.end() - pos)
+
+        gap_pos = array(gap_pos)
+        cum_lengths = array(cum_lengths).cumsum()
+        gap_pos[1:] = gap_pos[1:] - cum_lengths[:-1]
+
         seq = self.__class__(
-            "".join(gapless), name=self.get_name(), info=self.info, preserve_case=True
+            gap.sub("", seq), name=self.get_name(), info=self.info, preserve_case=True
+        )
+        indel_map = IndelMap(
+            gap_pos=gap_pos, cum_gap_lengths=cum_lengths, parent_length=len(seq)
         )
         seq.annotation_db = self.annotation_db
-        return map, seq
+        return indel_map, seq
 
     def replace(self, oldchar, newchar):
         """return new instance with oldchar replaced by newchar"""
@@ -1728,7 +1737,7 @@ class NucleicAcidSequence(Sequence):
         if not gc.is_stop(end):
             return self
 
-        if not len(m.gaps()):
+        if not m.num_gaps:
             # has zero length if no gaps
             return self[:-3]
 
@@ -1933,7 +1942,8 @@ class SeqView:
 
     @offset.setter
     def offset(self, value: int):
-        self._offset = value or 0
+        value = value or 0
+        self._offset = int(value)
 
     @property
     def seqid(self) -> str:
@@ -2329,7 +2339,7 @@ class SeqView:
             start, stop = self.start, self.stop
 
         data["init_args"]["seq"] = self.seq[start:stop]
-        data["init_args"]["offset"] = self.parent_start
+        data["init_args"]["offset"] = int(self.parent_start)
         data["init_args"]["seqid"] = self.seqid
         return data
 
diff --git a/src/cogent3/draw/dotplot.py b/src/cogent3/draw/dotplot.py
index 00c67c48e..c08795777 100644
--- a/src/cogent3/draw/dotplot.py
+++ b/src/cogent3/draw/dotplot.py
@@ -27,14 +27,9 @@ def suitable_threshold(window, desired_probability):
     return matches
 
 
-def len_seq(span):
-    """length of a Annotatable map object"""
-    return len(span.nongap())
-
-
 def not_gap(span):
     """whether a span corresponds to a non-gap"""
-    return len(span.gaps()) == 0
+    return span.num_gaps == 0
 
 
 def _convert_input(seq, moltype):
@@ -77,20 +72,23 @@ def get_align_coords(map1, map2, aligned=False) -> MatchedSeqPaths:
         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 = len_seq(map1[:i])
-            start_y = len_seq(map2[:i])
+            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, len_seq(map1[:i]) - 1),
-                    segment(start_y, len_seq(map2[:i]) - 1),
+                    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, len_seq(map1) - 1), segment(start_y, len_seq(map2) - 1))
+            (
+                segment(start_x, map1.parent_length - 1),
+                segment(start_y, map2.parent_length - 1),
+            )
         )
 
     return paths
diff --git a/src/cogent3/parse/cigar.py b/src/cogent3/parse/cigar.py
index 6c6dc28f6..0abfabc2d 100644
--- a/src/cogent3/parse/cigar.py
+++ b/src/cogent3/parse/cigar.py
@@ -49,7 +49,7 @@ def cigar_to_map(cigar_text):
             posn += n
         else:
             spans.append(LostSpan(n))
-    return IndelMap(spans=spans, parent_length=posn)
+    return IndelMap.from_spans(spans=spans, parent_length=posn)
 
 
 def aligned_from_cigar(cigar_text, seq, moltype=DNA):
@@ -62,16 +62,18 @@ def aligned_from_cigar(cigar_text, seq, moltype=DNA):
 
 def _slice_by_aln(map, left, right):
     slicemap = map[left:right]
-    location = [slicemap.start, slicemap.end] if hasattr(slicemap, "start") else []
+    location = [
+        map.get_seq_index(left),
+        map.get_seq_index(right),
+    ]
     return slicemap, location
 
 
 def _slice_by_seq(map, start, end):
-    re_map = map.inverse()
-    slicemap = re_map[start:end]
-    aln_start, aln_end = slicemap.start, slicemap.end
-    new_map = map[aln_start:aln_end]
-    return new_map, [aln_start, aln_end]
+    # start, end are in sequence coords
+    aln_start = map.get_align_index(start)
+    aln_end = map.get_align_index(end, slice_stop=True)
+    return map[aln_start:aln_end], [aln_start, aln_end]
 
 
 def _remap(map):
@@ -136,10 +138,11 @@ def CigarParser(
         ]:
             m, seq_loc = slice_cigar(cigars[seqname], aln_loc[0], aln_loc[1])
             if seq_loc:
-                seq = seqs[seqname]
+                start, stop = seq_loc
+                seq = seqs[seqname][start:stop]
                 if isinstance(seq, str):
                     seq = moltype.make_seq(seq)
-                data[seqname] = seq[seq_loc[0] : seq_loc[1]].gapped_by_map(m)
+                data[seqname] = seq.gapped_by_map(m)
             else:
                 data[seqname] = DNA.make_seq("-" * (aln_loc[1] - aln_loc[0]))
     return make_aligned_seqs(data)
diff --git a/src/cogent3/parse/gbseq.py b/src/cogent3/parse/gbseq.py
index 3e3476599..a917dbd2c 100644
--- a/src/cogent3/parse/gbseq.py
+++ b/src/cogent3/parse/gbseq.py
@@ -64,7 +64,9 @@ def GbSeqXmlParser(doc):
 
         seq = alphabet.make_seq(raw_string, name=name)
 
-        feat = location.FeatureMap(locations=[(0, len(seq))], parent_length=len(seq))
+        feat = location.FeatureMap.from_locations(
+            locations=[(0, len(seq))], parent_length=len(seq)
+        )
         seq.add_feature(biotype="source", name=name, spans=feat.get_coordinates())
 
         organism = str(
diff --git a/src/cogent3/util/parallel.py b/src/cogent3/util/parallel.py
index a1251a094..cf81811ac 100644
--- a/src/cogent3/util/parallel.py
+++ b/src/cogent3/util/parallel.py
@@ -216,9 +216,8 @@ def _as_completed_mpi(f, s, max_workers, if_serial, chunksize=None):
 
 def _as_completed_mproc(f, s, max_workers):
     """multiprocess version of as_completed"""
-    if not max_workers:
+    if not max_workers or max_workers > multiprocessing.cpu_count():
         max_workers = multiprocessing.cpu_count() - 1
-    assert max_workers < multiprocessing.cpu_count()
 
     f = PicklableAndCallable(f)
 


Test patch:
diff --git a/tests/test_app/test_align.py b/tests/test_app/test_align.py
index 9e53d9b4e..288a84230 100644
--- a/tests/test_app/test_align.py
+++ b/tests/test_app/test_align.py
@@ -10,6 +10,7 @@ from cogent3 import (
     DNA,
     get_app,
     get_moltype,
+    load_aligned_seqs,
     make_aligned_seqs,
     make_tree,
     make_unaligned_seqs,
@@ -148,7 +149,7 @@ class RefalignmentTests(TestCase):
         )
         expect = orig.to_dict()
         aligner = align_app.align_to_ref(ref_seq="Ref")
-        aln = aligner(orig.degap())
+        aln = aligner.main(orig.degap())
         self.assertEqual(aln.to_dict(), expect)
 
     def test_gap_union(self):
@@ -726,3 +727,31 @@ def test_smith_waterman_raises(seqs):
     coll = make_unaligned_seqs(data=[seqs.get_seq("Human")], moltype="dna")
     aln = aligner(coll)
     assert isinstance(aln, NotCompleted)
+
+
+def test_aln_two():
+    """correctly recapitulates known case"""
+    orig = make_aligned_seqs(
+        {
+            "Ref": "CAGGAGAACAGAAACCCATTACTCACT",
+            "Qu7": "CAGGA--ACAGA--CCCGTTA---ACT",
+        },
+        moltype="dna",
+    )
+    expect = orig.to_dict()
+    aligner = align_app.align_to_ref(ref_seq="Ref")
+    seqs = orig.degap()
+    aln = aligner.main(seqs)
+    assert aln.to_dict() == expect
+
+
+def test_codon_incomplete(DATA_DIR):
+    names = ["FlyingFox", "DogFaced", "FreeTaile"]
+    aln = load_aligned_seqs(DATA_DIR / "brca1.fasta", moltype="dna")
+    seqs = aln.take_seqs(names)[2700:3000].degap()
+    aligner = align_app.progressive_align("codon")
+    aln = aligner(seqs)
+    assert aln  # will fail if aln is a NotCompleted instance
+    # 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}
diff --git a/tests/test_core/test_alignment.py b/tests/test_core/test_alignment.py
index 84c79fa00..06fe5fd72 100644
--- a/tests/test_core/test_alignment.py
+++ b/tests/test_core/test_alignment.py
@@ -455,54 +455,6 @@ class SequenceCollectionBaseTests(object):
         align = align1 + align2
         self.assertEqual(align.info["key"], "foo")
 
-    def test_add_seqs(self):
-        """add_seqs should return an alignment with the new sequences appended or inserted"""
-        data = [("name1", "AAA"), ("name2", "AAA"), ("name3", "AAA"), ("name4", "AAA")]
-        data1 = [("name1", "AAA"), ("name2", "AAA")]
-        data2 = [("name3", "AAA"), ("name4", "AAA")]
-        data3 = [("name5", "BBB"), ("name6", "CCC")]
-        aln = self.Class(data)
-        aln3 = self.Class(data3)
-
-        out_aln = aln.add_seqs(aln3)
-        # test append at the end
-        self.assertEqual(str(out_aln), str(self.Class(data + data3)))
-
-        out_aln = aln.add_seqs(aln3, before_name="name3")
-        self.assertEqual(
-            str(out_aln), str(self.Class(data1 + data3 + data2))
-        )  # test insert before
-
-        out_aln = aln.add_seqs(aln3, after_name="name2")
-        self.assertEqual(
-            str(out_aln), str(self.Class(data1 + data3 + data2))
-        )  # test insert after
-
-        out_aln = aln.add_seqs(aln3, before_name="name1")
-        # test if insert before first seq works
-        self.assertEqual(str(out_aln), str(self.Class(data3 + data)))
-
-        out_aln = aln.add_seqs(aln3, after_name="name4")
-        # test if insert after last seq works
-        self.assertEqual(str(out_aln), str(self.Class(data + data3)))
-
-        self.assertRaises(
-            ValueError, aln.add_seqs, aln3, before_name="name5"
-        )  # wrong after/before name
-        self.assertRaises(
-            ValueError, aln.add_seqs, aln3, after_name="name5"
-        )  # wrong after/before name
-
-        if isinstance(aln, Alignment) or isinstance(aln, ArrayAlignment):
-            self.assertRaises((DataError, ValueError), aln.add_seqs, aln3 + aln3)
-        else:
-            exp = set([seq for name, seq in data])
-            exp.update([seq + seq for name, seq in data3])
-            got = set()
-            for seq in aln.add_seqs(aln3 + aln3).seqs:
-                got.update([str(seq).strip()])
-            self.assertEqual(got, exp)
-
     def test_add_seqs_info(self):
         """add_seqs should preserve info attribute"""
         data = [("name1", "AAA"), ("name2", "AAA"), ("name3", "AAA"), ("name4", "AAA")]
@@ -823,23 +775,6 @@ class SequenceCollectionBaseTests(object):
                 dict(num_seqs=10, num_pos=60, ref_name="longest", wrap=60),
             )
 
-    def test_set_repr_policy_valid_input(self):
-        """repr_policy should be set to new values"""
-        seqs = self.Class({"a": "AAAAA", "b": "AAA--"})
-        seqs.set_repr_policy(num_seqs=5, num_pos=40, ref_name="a", wrap=10)
-        self.assertEqual(
-            seqs._repr_policy, dict(num_seqs=5, num_pos=40, ref_name="a", wrap=10)
-        )
-
-        if self.Class == SequenceCollection:
-            # this class cannot slice
-            return
-
-        # should persist in slicing
-        self.assertEqual(
-            seqs[:2]._repr_policy, dict(num_seqs=5, num_pos=40, ref_name="a", wrap=10)
-        )
-
     def test_set_wrap_affects_repr_html(self):
         """the wrap argument affects the number of columns"""
         if self.Class == SequenceCollection:
@@ -1070,93 +1005,6 @@ class AlignmentBaseTests(SequenceCollectionBaseTests):
             {"a": "AA", "b": "A-", "c": "AA"},
         )
 
-    def test_no_degenerates(self):
-        """no_degenerates correctly excludes columns containing IUPAC ambiguity codes"""
-        data = {
-            "s1": "AAA CCC GGG TTT".replace(" ", ""),
-            "s2": "CCC GGG T-T AAA".replace(" ", ""),
-            "s3": "GGR YTT AAA CCC".replace(" ", ""),
-        }
-        aln = self.Class(data=data, moltype=DNA)
-
-        # motif length of 1, defaults - no gaps allowed
-        result = aln.no_degenerates().to_dict()
-        expect = {
-            "s1": "AA CC GG TTT".replace(" ", ""),
-            "s2": "CC GG TT AAA".replace(" ", ""),
-            "s3": "GG TT AA CCC".replace(" ", ""),
-        }
-        self.assertEqual(result, expect)
-
-        # allow gaps
-        result = aln.no_degenerates(allow_gap=True).to_dict()
-        expect = {
-            "s1": "AA CC GGG TTT".replace(" ", ""),
-            "s2": "CC GG T-T AAA".replace(" ", ""),
-            "s3": "GG TT AAA CCC".replace(" ", ""),
-        }
-        self.assertEqual(result, expect)
-
-        # motif length of 3, defaults - no gaps allowed
-        result = aln.no_degenerates(motif_length=3).to_dict()
-        expect = {
-            "s1": "TTT".replace(" ", ""),
-            "s2": "AAA".replace(" ", ""),
-            "s3": "CCC".replace(" ", ""),
-        }
-        self.assertEqual(result, expect)
-
-        # allow gaps
-        result = aln.no_degenerates(motif_length=3, allow_gap=True).to_dict()
-        expect = {
-            "s1": "GGG TTT".replace(" ", ""),
-            "s2": "T-T AAA".replace(" ", ""),
-            "s3": "AAA CCC".replace(" ", ""),
-        }
-        self.assertEqual(result, expect)
-
-        # raises ValueError if a default moltype -- with no
-        # degen characters -- is used
-        aln = self.Class(data=data)
-        self.assertRaises(ValueError, aln.no_degenerates)
-
-    def test_omit_gap_pos(self):
-        """Alignment omit_gap_pos should return alignment w/o positions of gaps"""
-        aln = self.end_gaps
-        # first, check behavior when we're just acting on the cols (and not
-        # trying to delete the naughty seqs).
-
-        # default should strip out cols that are 100% gaps
-        result = aln.omit_gap_pos()
-        self.assertEqual(result.to_dict(), {"a": "-ABC", "b": "CBA-", "c": "-DEF"})
-        # if allowed_gap_frac is 1, shouldn't delete anything
-        self.assertEqual(
-            aln.omit_gap_pos(1).to_dict(),
-            {"a": "--A-BC-", "b": "-CB-A--", "c": "--D-EF-"},
-        )
-        # if allowed_gap_frac is 0, should strip out any cols containing gaps
-        self.assertEqual(
-            aln.omit_gap_pos(0).to_dict(), {"a": "AB", "b": "BA", "c": "DE"}
-        )
-        # intermediate numbers should work as expected
-        self.assertEqual(
-            aln.omit_gap_pos(0.4).to_dict(), {"a": "ABC", "b": "BA-", "c": "DEF"}
-        )
-        self.assertEqual(
-            aln.omit_gap_pos(0.7).to_dict(), {"a": "-ABC", "b": "CBA-", "c": "-DEF"}
-        )
-
-        # when we increase the number of sequences to 6, more differences
-        # start to appear.
-        new_aln_data = aln.named_seqs.copy()
-        new_aln_data["d"] = "-------"
-        new_aln_data["e"] = "XYZXYZX"
-        new_aln_data["f"] = "AB-CDEF"
-        aln = self.Class(new_aln_data)
-
-        # if no gaps are allowed, we get None
-        self.assertEqual(aln.omit_gap_pos(0), None)
-
     def test_omit_gap_pos2(self):
         """consistency with different motif_length values"""
         data = {
@@ -1416,41 +1264,6 @@ class AlignmentBaseTests(SequenceCollectionBaseTests):
         self.assertEqual(aln.variable_positions(include_gap_motif=True), [])
         self.assertEqual(aln.variable_positions(include_gap_motif=False), [])
 
-    def test_to_type(self):
-        """correctly interconvert between alignment types"""
-        new_seqs = {"seq1": "ACGTACGTA", "seq2": "ACCGAA---", "seq3": "ACGTACGTT"}
-        array_align = self.Class == ArrayAlignment
-        # when array_align arg matches instance class, no conversion
-        # and get back self
-        aln = self.Class(data=new_seqs)
-        new = aln.to_type(array_align=array_align)
-        self.assertEqual(id(aln), id(new))
-
-        # when array_align arg does not match, should get back the opposite type
-        new = aln.to_type(array_align=not array_align)
-        self.assertFalse(isinstance(new, self.Class))
-
-        # we should be able to specify moltype and alignment
-        new = aln.to_type(array_align=not array_align, moltype=DNA)
-        self.assertEqual(new.to_dict(), new_seqs)
-        # and translate
-        self.assertEqual(
-            new.get_translation().to_dict(),
-            {"seq1": "TYV", "seq3": "TYV", "seq2": "TE-"},
-        )
-
-        # should work on ArrayAlign when just moltype changes
-        new_seqs = {"seq1": "ACGTACGTA", "seq2": "ACCGAA---"}
-        aln = self.Class(data=new_seqs)
-        new = aln.to_type(array_align=array_align, moltype=DNA)
-        new = new.no_degenerates()  # this should not fail!
-        self.assertEqual(len(new), len(aln) - 3)
-
-        # should correctly apply to existing moltype
-        aln = self.Class(data=new_seqs, moltype=DNA)
-        new = aln.to_type(array_align=not array_align)
-        self.assertEqual(aln.moltype, new.moltype)
-
     def test_to_type_info(self):
         """interconverting between alignment types preserves info attribute"""
         new_seqs = {"seq1": "ACGTACGTA", "seq2": "ACCGAA---", "seq3": "ACGTACGTT"}
@@ -1665,19 +1478,6 @@ class AlignmentBaseTests(SequenceCollectionBaseTests):
         # but setting drop_invalid=False allows calc
         dists = aln.distance_matrix(calc="paralinear", drop_invalid=True)
 
-    def test_quick_tree(self):
-        """quick tree method returns tree"""
-        aln = self.Class(self.brca1_data, moltype=DNA)
-        aln = aln.take_seqs(aln.names[:5])[:100]
-        tree = aln.quick_tree(calc="hamming", show_progress=False)
-        # bootstrap
-        tree = aln.quick_tree(calc="hamming", bootstrap=2, show_progress=False)
-        self.assertEqual(set(tree.get_tip_names()), set(aln.names))
-        for edge in tree.preorder():
-            if edge.is_root():
-                continue
-            self.assertIsInstance(edge.params["support"], float)
-
     def test_get_gapped_seq(self):
         """Alignment.get_gapped_seq should return seq, with gaps"""
         aln = self.Class({"seq1": "--TTT?", "seq2": "GATC??"})
@@ -2179,6 +1979,61 @@ class IntegrationTests(TestCase):
         self.assertEqual(self.r1.name, "x")
 
 
+@pytest.mark.parametrize(
+    "raw_seq,coords",
+    (("ACGGTAAAG", ((2, 4), (5, 8))), ("CCC---CCC", ((0, 3), (6, 9)))),
+)
+def test_featuremap_slice_aligned(raw_seq, coords):
+    from cogent3.core.alignment import Aligned
+    from cogent3.core.location import FeatureMap, Span
+
+    im, seq = DNA.make_seq(raw_seq).parse_out_gaps()
+    ia = Aligned(im, seq)
+    length = len(raw_seq)
+    fmap = FeatureMap(spans=[Span(s, e) for s, e in coords], parent_length=length)
+    expect = "".join(raw_seq[s:e] for s, e in fmap.get_coordinates())
+    got = ia[fmap]
+    assert str(got) == expect
+
+
+@pytest.mark.parametrize("cls", (ArrayAlignment, Alignment))
+def test_to_type(cls):
+    """correctly interconvert between alignment types"""
+    new_seqs = {"seq1": "ACGTACGTA", "seq2": "ACCGAA---", "seq3": "ACGTACGTT"}
+    array_align = cls == ArrayAlignment
+    # when array_align arg matches instance class, no conversion
+    # and get back self
+    aln = cls(data=new_seqs)
+    new = aln.to_type(array_align=array_align)
+    assert id(aln) == id(new)
+
+    # when array_align arg does not match, should get back the opposite type
+    new = aln.to_type(array_align=not array_align)
+    assert not isinstance(new, cls)
+
+    # we should be able to specify moltype and alignment
+    new = aln.to_type(array_align=not array_align, moltype=DNA)
+    assert new.to_dict() == new_seqs
+    # and translate
+    assert new.get_translation().to_dict() == {
+        "seq1": "TYV",
+        "seq3": "TYV",
+        "seq2": "TE-",
+    }
+
+    # should work on ArrayAlign when just moltype changes
+    new_seqs = {"seq1": "ACGTACGTA", "seq2": "ACCGAA---"}
+    aln = cls(data=new_seqs)
+    new = aln.to_type(array_align=array_align, moltype=DNA)
+    new = new.no_degenerates()  # this should not fail!
+    assert len(new) == len(aln) - 3
+
+    # should correctly apply to existing moltype
+    aln = cls(data=new_seqs, moltype=DNA)
+    new = aln.to_type(array_align=not array_align)
+    assert aln.moltype == new.moltype
+
+
 @pytest.mark.parametrize(
     "moltype,array_align",
     tuple(itertools.product(["rna", "dna", "protein"], [True, False])),
@@ -2961,7 +2816,8 @@ def test_get_seq_with_sliced_rced_aln_multiple_spans(name):
     }
     aln = make_aligned_seqs(data=seqs, moltype="dna", array_align=False)
     start, stop = 1, 10
-    a1 = aln[start:stop].rc()
+    a1 = aln[start:stop]
+    a1 = a1.rc()
     got = str(a1.get_seq(name))
     dna = get_moltype("dna")
     expect = dna.complement(seqs[name][start:stop].replace("-", ""))[::-1]
@@ -3613,7 +3469,202 @@ def test_positions(cls):
     assert list(r.positions) == expect
 
 
+@pytest.mark.parametrize("cls", (Alignment, ArrayAlignment, SequenceCollection))
+def test_set_repr_policy_valid_input(cls):
+    """repr_policy should be set to new values"""
+    seqs = cls({"a": "AAAAA", "b": "AAA--"})
+    seqs.set_repr_policy(num_seqs=5, num_pos=40, ref_name="a", wrap=10)
+    assert seqs._repr_policy == dict(num_seqs=5, num_pos=40, ref_name="a", wrap=10)
+
+
+@pytest.mark.parametrize("cls", (Alignment, ArrayAlignment))
+def test_set_repr_policy_valid_input_slices(cls):
+    """repr_policy should be set to new values"""
+    seqs = cls({"a": "AAAAA", "b": "AAA--"})
+    seqs.set_repr_policy(num_seqs=5, num_pos=40, ref_name="a", wrap=10)
+    # should persist in slicing
+    sliced = seqs[:2]
+    assert sliced._repr_policy == dict(num_seqs=5, num_pos=40, ref_name="a", wrap=10)
+
+
 def test_array_align_error_with_mixed_length():
     data = dict(s1="ACGG", s2="A-G")
     with pytest.raises(ValueError, match=".* not all the same length.*"):
         make_aligned_seqs(data=data)
+
+
+@pytest.mark.parametrize("cls", (Alignment, ArrayAlignment, SequenceCollection))
+def test_add_seqs(cls):
+    """add_seqs should return an alignment with the new sequences appended or inserted"""
+    data = [("name1", "AAA"), ("name2", "AAA"), ("name3", "AAA"), ("name4", "AAA")]
+    data1 = [("name1", "AAA"), ("name2", "AAA")]
+    data2 = [("name3", "AAA"), ("name4", "AAA")]
+    data3 = [("name5", "BBB"), ("name6", "CCC")]
+    aln = cls(data)
+    aln3 = cls(data3)
+
+    out_aln = aln.add_seqs(aln3)
+    # test append at the end
+    assert str(out_aln) == str(cls(data + data3))
+
+    out_aln = aln.add_seqs(aln3, before_name="name3")
+    assert str(out_aln) == str(cls(data1 + data3 + data2))
+
+    # test insert before
+
+    out_aln = aln.add_seqs(aln3, after_name="name2")
+    assert str(out_aln) == str(cls(data1 + data3 + data2))  # test insert after
+
+    out_aln = aln.add_seqs(aln3, before_name="name1")
+    # test if insert before first seq works
+    assert str(out_aln) == str(cls(data3 + data))
+
+    out_aln = aln.add_seqs(aln3, after_name="name4")
+    # test if insert after last seq works
+    assert str(out_aln) == str(cls(data + data3))
+
+    with pytest.raises(ValueError):
+        # wrong after/before name
+        aln.add_seqs(aln3, before_name="name5")
+
+    with pytest.raises(ValueError):
+        # wrong after/before name
+        aln.add_seqs(aln3, after_name="name5")
+
+    if isinstance(aln, Alignment) or isinstance(aln, ArrayAlignment):
+        with pytest.raises((DataError, ValueError)):
+            aln.add_seqs(aln3 + aln3)
+    else:
+        exp = set([seq for name, seq in data])
+        exp.update([seq + seq for name, seq in data3])
+        got = set()
+        for seq in aln.add_seqs(aln3 + aln3).seqs:
+            got.update([str(seq).strip()])
+        assert got == exp
+
+
+@pytest.mark.parametrize("cls", (Alignment, ArrayAlignment))
+def test_omit_gap_pos(cls):
+    """Alignment omit_gap_pos should return alignment w/o positions of gaps"""
+    aln = cls({"a": "--A-BC-", "b": "-CB-A--", "c": "--D-EF-"}, names=["a", "b", "c"])
+    # first, check behavior when we're just acting on the cols (and not
+    # trying to delete the naughty seqs).
+
+    # default should strip out cols that are 100% gaps
+    result = aln.omit_gap_pos()
+    assert result.to_dict() == {"a": "-ABC", "b": "CBA-", "c": "-DEF"}
+    # if allowed_gap_frac is 1, shouldn't delete anything
+    assert aln.omit_gap_pos(1).to_dict() == {
+        "a": "--A-BC-",
+        "b": "-CB-A--",
+        "c": "--D-EF-",
+    }
+
+    # if allowed_gap_frac is 0, should strip out any cols containing gaps
+    assert aln.omit_gap_pos(0).to_dict() == {"a": "AB", "b": "BA", "c": "DE"}
+    # intermediate numbers should work as expected
+    assert aln.omit_gap_pos(0.4).to_dict() == {"a": "ABC", "b": "BA-", "c": "DEF"}
+    assert aln.omit_gap_pos(0.7).to_dict() == {"a": "-ABC", "b": "CBA-", "c": "-DEF"}
+
+    # when we increase the number of sequences to 6, more differences
+    # start to appear.
+    new_aln_data = aln.named_seqs.copy()
+    new_aln_data["d"] = "-------"
+    new_aln_data["e"] = "XYZXYZX"
+    new_aln_data["f"] = "AB-CDEF"
+    aln = cls(new_aln_data)
+
+    # if no gaps are allowed, we get None
+    assert aln.omit_gap_pos(0) is None
+
+
+@pytest.mark.parametrize("cls", (Alignment, ArrayAlignment))
+def test_no_degenerates(cls):
+    """no_degenerates correctly excludes columns containing IUPAC ambiguity codes"""
+    data = {
+        "s1": "AAA CCC GGG TTT".replace(" ", ""),
+        "s2": "CCC GGG T-T AAA".replace(" ", ""),
+        "s3": "GGR YTT AAA CCC".replace(" ", ""),
+    }
+    aln = cls(data=data, moltype=DNA)
+
+    # motif length of 1, defaults - no gaps allowed
+    result = aln.no_degenerates().to_dict()
+    expect = {
+        "s1": "AA CC GG TTT".replace(" ", ""),
+        "s2": "CC GG TT AAA".replace(" ", ""),
+        "s3": "GG TT AA CCC".replace(" ", ""),
+    }
+    assert result == expect
+
+    # allow gaps
+    result = aln.no_degenerates(allow_gap=True).to_dict()
+    expect = {
+        "s1": "AA CC GGG TTT".replace(" ", ""),
+        "s2": "CC GG T-T AAA".replace(" ", ""),
+        "s3": "GG TT AAA CCC".replace(" ", ""),
+    }
+    assert result == expect
+
+    # motif length of 3, defaults - no gaps allowed
+    result = aln.no_degenerates(motif_length=3, allow_gap=False).to_dict()
+    expect = {
+        "s1": "TTT".replace(" ", ""),
+        "s2": "AAA".replace(" ", ""),
+        "s3": "CCC".replace(" ", ""),
+    }
+    assert result == expect
+
+    # allow gaps
+    result = aln.no_degenerates(motif_length=3, allow_gap=True).to_dict()
+    expect = {
+        "s1": "GGG TTT".replace(" ", ""),
+        "s2": "T-T AAA".replace(" ", ""),
+        "s3": "AAA CCC".replace(" ", ""),
+    }
+    assert result == expect
+
+
+@pytest.mark.parametrize("cls", (Alignment, ArrayAlignment))
+@pytest.mark.parametrize("moltype", ("bytes", "text"))
+def test_no_degenerates_invalid_moltype(cls, moltype):
+    # raises ValueError if a default moltype -- with no
+    # degen characters -- is used
+    data = {
+        "s1": "AAA CCC GGG TTT".replace(" ", ""),
+        "s2": "CCC GGG T-T AAA".replace(" ", ""),
+        "s3": "GGR YTT AAA CCC".replace(" ", ""),
+    }
+
+    aln = cls(data=data, moltype=moltype)
+    with pytest.raises(ValueError):
+        aln.no_degenerates()
+
+
+@pytest.mark.parametrize("cls", (Alignment, ArrayAlignment))
+@pytest.mark.parametrize("calc", ("hamming", None))
+def test_quick_tree(cls, calc, brca1_data):
+    """quick tree method returns tree"""
+    aln = cls(brca1_data, moltype=DNA)[:100]
+    aln = aln.take_seqs(["Human", "Rhesus", "HowlerMon", "Galago", "Mouse"])
+    kwargs = dict(show_progress=False)
+    if calc:
+        kwargs["calc"] = calc
+
+    # bootstrap
+    tree = aln.quick_tree(bootstrap=2, **kwargs)
+    assert set(tree.get_tip_names()) == set(aln.names)
+    types = {
+        type(float(edge.params["support"]))
+        for edge in tree.preorder()
+        if not edge.is_root()
+    }
+    assert types == {float}
+
+
+@pytest.mark.parametrize("raw", ("-AAAGGGGGAACCCT", "AAAGGGGGAACCCT"))
+def test_slice_aligned(raw):
+    imap, seq = DNA.make_seq(raw, name="x").parse_out_gaps()
+    al = Aligned(imap, seq)
+    sliced = al[:-3]
+    assert str(sliced) == raw[:-3]
diff --git a/tests/test_core/test_annotation_db.py b/tests/test_core/test_annotation_db.py
index 05eb53188..a0910bb13 100644
--- a/tests/test_core/test_annotation_db.py
+++ b/tests/test_core/test_annotation_db.py
@@ -413,7 +413,9 @@ def test_feature_nucleic():
     #                              111111
     #                    0123456789012345
     seq = make_seq("AACCTTTGGGGAATTT", moltype="dna")
-    mmap = loc.FeatureMap(locations=[(4, 7), (10, 12)], parent_length=len(seq))
+    mmap = loc.FeatureMap.from_locations(
+        locations=[(4, 7), (10, 12)], parent_length=len(seq)
+    )
     expect = seq[mmap]
 
     rcseq = seq.rc()
diff --git a/tests/test_core/test_features.py b/tests/test_core/test_features.py
index 5c1952762..40290d0e6 100755
--- a/tests/test_core/test_features.py
+++ b/tests/test_core/test_features.py
@@ -1,5 +1,6 @@
 from unittest import TestCase
 
+import numpy
 import pytest
 
 from cogent3 import ASCII, DNA, get_moltype, make_aligned_seqs
@@ -135,7 +136,9 @@ def test_feature_residue():
     db.add_feature(seqid="x", biotype="exon", name="ex1", spans=[(0, 4)])
     aln_exons = list(aln.get_features(seqid="x", biotype="exon"))
     assert "biotype='exon', name='ex1', map=[0:1, 2:5]/10" in str(aln_exons)
-    assert aln_exons[0].get_slice().to_dict() == dict(x="CCCC", y="----")
+    exon = aln_exons[0]
+    exon_seq = exon.get_slice()
+    assert exon_seq.to_dict() == dict(x="CCCC", y="----")
     # Feature.as_one_span(), is applied to the exon that
     # straddles the gap in x. The result is we preserve that feature.
     exon_full_aln = aln_exons[0].as_one_span()
@@ -621,6 +624,21 @@ def test_feature_out_range():
     assert not f
 
 
+@pytest.mark.parametrize("cast", (list, numpy.array))
+def test_search_with_ints(cast):
+    """searching for features with numpy ints should work"""
+    start, stop = cast([2, 5])
+    seq = DNA.make_seq("AAAGGGGGAACCCT", name="x")
+    db = GffAnnotationDb()
+    db.add_feature(seqid="x", biotype="exon", name="E1", spans=[(3, 8)])
+    db.add_feature(seqid="x", biotype="exon", name="E2", spans=[(10, 13)])
+    seq.annotation_db = db
+    feats = list(
+        seq.get_features(biotype="exon", allow_partial=True, start=start, stop=stop)
+    )
+    assert len(feats) == 1
+
+
 def test_roundtripped_alignment_with_slices():
     """Sliced Alignment with annotations roundtrips correctly"""
     # annotations just on member sequences
@@ -632,9 +650,15 @@ def test_roundtripped_alignment_with_slices():
     db.add_feature(seqid="x", biotype="exon", name="E2", spans=[(10, 13)])
     aln.annotation_db = db
     # at the alignment level
+    sl = aln.seqs[0][:-3]
+    assert str(sl) == "-AAAGGGGGAACCCT"[:-3]
     sub_aln = aln[:-3]
+    feats = list(sub_aln.get_features(biotype="exon", allow_partial=True))
+    assert len(feats) == 2
     new = deserialise_object(sub_aln.to_json())
-    gf1, gf2 = list(new.get_features(biotype="exon", allow_partial=True))
+    feats = list(new.get_features(biotype="exon", allow_partial=True))
+    assert len(feats) == 2
+    gf1, gf2 = feats
     assert gf1.get_slice().to_dict() == {"x": "GGGGG", "y": "--TTT"}
     assert gf2.get_slice().to_dict() == {"x": "C", "y": "G"}
 
diff --git a/tests/test_core/test_location.py b/tests/test_core/test_location.py
index 61273e9d3..ece8c9a84 100644
--- a/tests/test_core/test_location.py
+++ b/tests/test_core/test_location.py
@@ -1,10 +1,12 @@
 """Unit tests for Span classes.
 """
+from itertools import combinations
 from unittest import TestCase
 
+import numpy
 import pytest
 
-from cogent3 import DNA
+from cogent3 import DNA, make_seq
 from cogent3.core.location import (
     FeatureMap,
     IndelMap,
@@ -269,42 +271,15 @@ class MapTests(TestCase):
         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))
-
-
-@pytest.mark.parametrize("cls", (IndelMap, FeatureMap))
-def test_map_plus_position(cls):
+def test_map_plus_position():
     # seq is 9 long
     # plus coords  012345678
     # +slice         **
     # plus seq     AAACCCTGG
 
     # orig = Aligned(*DNA.make_seq("AAACCCTGG", name="a").parse_out_gaps())
-    orig = cls(locations=[(0, 9)], parent_length=9)
+    orig = FeatureMap.from_locations(locations=[(0, 9)], parent_length=9)
     assert orig.absolute_position(2) == 2
     assert orig.absolute_position(6) == 6
 
@@ -329,22 +304,19 @@ def test_map_plus_position(cls):
 
 
 def test_indel_map_useful_complete():
-    im = IndelMap(spans=[LostSpan(3)], parent_length=0)
+    im = IndelMap.from_spans(spans=[LostSpan(3)], parent_length=0)
     assert not im.useful
     assert not im.complete
-    assert len(im) == im.length == 3
+    assert len(im) == 3
 
 
-@pytest.mark.parametrize("cls", (FeatureMap, IndelMap))
-def test_map_nucleic_reversed(cls):
+def test_map_nucleic_reversed():
     expect = [(0, 9)]
     # seq is 9 long
     # plus coords  012345678
-    # +slice         **
     # plus seq     AAACCCTGG
 
-    # orig = Aligned(*DNA.make_seq("AAACCCTGG", name="a").parse_out_gaps())
-    orig = cls(locations=[(0, 9)], parent_length=9)
+    orig = FeatureMap.from_locations(locations=[(0, 9)], parent_length=9)
     # minus coords 012345678
     # rel coord      01234
     # -slice         *****
@@ -361,7 +333,7 @@ def test_coordinate(cls):
     #                   01   2 345
     seq = DNA.make_seq("AC---G-TAA--")
     m, _ = seq.parse_out_gaps()
-    m = cls(spans=tuple(m.spans), parent_length=m.parent_length)
+    m = cls.from_spans(spans=tuple(m.spans), parent_length=m.parent_length)
     got = m.get_coordinates()
     assert got == [(0, 2), (2, 3), (3, 6)]
 
@@ -370,19 +342,18 @@ def test_coordinate(cls):
 def test_gap_coordinate(cls):
     seq = DNA.make_seq("AC---G-TAA--")
     m, _ = seq.parse_out_gaps()
-    m = cls(spans=tuple(m.spans), parent_length=m.parent_length)
+    m = cls.from_spans(spans=tuple(m.spans), parent_length=m.parent_length)
     got = m.get_gap_coordinates()
-    assert got == [(2, 3), (3, 1), (6, 2)]
+    assert [tuple(c) for c in got] == [(2, 3), (3, 1), (6, 2)]
 
 
-@pytest.mark.parametrize("cls", (FeatureMap, IndelMap))
-def test_gaps(cls):
+def test_gaps():
     # returns spans corresponding to position on "aligned" seq of gaps
     #                   000000000011
     #                   012345678901
     seq = DNA.make_seq("AC---G-TAA--")
     m, s = seq.parse_out_gaps()
-    m = cls(spans=tuple(m.spans), parent_length=m.parent_length)
+    m = FeatureMap.from_spans(spans=tuple(m.spans), parent_length=m.parent_length)
     got = [(g.start, g.end) for g in tuple(m.gaps().spans)]
     assert got == [(2, 5), (6, 7), (10, 12)]
 
@@ -394,41 +365,76 @@ def test_nongap(cls):
     #                   012345678901
     seq = DNA.make_seq("AC---G-TAA--")
     m, _ = seq.parse_out_gaps()
-    m = cls(spans=tuple(m.spans), parent_length=m.parent_length)
+    m = cls.from_spans(spans=tuple(m.spans), parent_length=m.parent_length)
 
-    got = [(g.start, g.end) for g in m.nongap().spans]
+    got = [(s.start, s.end) for s in m.nongap()]
     assert got == [(0, 2), (5, 6), (7, 10)]
 
 
+def test_spans_gen():
+    # returns spans corresponding to position on "aligned" seq of nongaps
+    #                   000000000011
+    #                   012345678901
+    seq = DNA.make_seq("AC---G-TAA--")
+    expect = [Span(0, 2), LostSpan(3), Span(2, 3), LostSpan(1), Span(3, 6), LostSpan(2)]
+    m, _ = seq.parse_out_gaps()
+    gap_data = numpy.array([(2, 3), (3, 1), (6, 2)])
+    pos, lengths = gap_data.T
+    im = IndelMap(gap_pos=pos, gap_lengths=lengths, parent_length=m.parent_length)
+    got = list(im.spans)
+    assert got == expect
+
+
+def test_spans_gap_start():
+    # returns spans corresponding to position on "aligned" seq of nongaps
+    #                   000000000011
+    #                   012345678901
+    seq = DNA.make_seq("---TAA")
+    expect = [LostSpan(3), Span(0, 3)]
+    im, _ = seq.parse_out_gaps()
+    got = list(im.spans)
+    assert got == expect
+
+
+def test_spans_gen_nogaps():
+    # returns spans corresponding to position on "aligned" seq of nongaps
+    #                   000000000011
+    #                   012345678901
+    seq = DNA.make_seq("ACGTAA")
+    m, _ = seq.parse_out_gaps()
+    spans = list(m.spans)
+    assert len(spans) == 1
+    assert len(spans[0]) == len(seq)
+
+
 def test_round_trip_rich_dict():
     seq = DNA.make_seq("AC---G-TAA--")
-    m, s = seq.parse_out_gaps()
-    # reversed() reverses the order of spans, but keeps their coordinates
-    # differs from nucleic reversed, which computes a new relative position
-    im = IndelMap(spans=m.spans, parent_length=m.parent_length)
+    im, s = seq.parse_out_gaps()
     rd = im.to_rich_dict()
     got = IndelMap.from_rich_dict(rd)
     assert im is not got
     assert got.to_rich_dict() == im.to_rich_dict()
 
 
-def test_serialisable_attr():
-    im = IndelMap(locations=[(0, 2)], parent_length=20)
-    set_vals = {"locations": [(0, 2)], "parent_length": 20}
-    got = {k: im._serialisable[k] for k in set_vals}
-    assert got == set_vals
-
-
 def test_terminal_unknown():
-    # span idx          01 2  345  6
-    seq = DNA.make_seq("-AC---G-TAA--")
-    m, _ = seq.parse_out_gaps()
+    # seq  idx   01   2 345  6
+    #           -AC---G-TAA--
+    # aligned seq length is 13
+    gap_data = numpy.array([[0, 1], [2, 3], [3, 1], [6, 2]])
+    gap_pos, gap_lengths = gap_data.T
+
+    m = IndelMap(
+        gap_pos=gap_pos.copy(), gap_lengths=gap_lengths.copy(), parent_length=6
+    )
     # not unknown, by default
     m_spans = tuple(m.spans)
     assert m_spans[0].lost and not isinstance(m_spans[0], TerminalPadding)
     # use the constructor arg
     m = IndelMap(
-        spans=tuple(m.spans), parent_length=m.parent_length, termini_unknown=True
+        gap_pos=gap_pos.copy(),
+        gap_lengths=gap_lengths.copy(),
+        parent_length=6,
+        termini_unknown=True,
     )
     m_spans = tuple(m.spans)
     assert isinstance(m_spans[0], TerminalPadding)
@@ -437,9 +443,8 @@ def test_terminal_unknown():
     assert m_spans[4].lost and not isinstance(m_spans[2], TerminalPadding)
 
     # use the method
-    m, _ = seq.parse_out_gaps()
     m = IndelMap(
-        spans=tuple(m.spans), parent_length=m.parent_length
+        gap_pos=gap_pos.copy(), gap_lengths=gap_lengths.copy(), parent_length=6
     ).with_termini_unknown()
     m_spans = tuple(m.spans)
     assert isinstance(m_spans[0], TerminalPadding)
@@ -448,25 +453,26 @@ def test_terminal_unknown():
     assert not isinstance(m_spans[2], TerminalPadding)
 
     # no gaps, no effect
-    seq = DNA.make_seq("ACGTAA")
-    m, _ = seq.parse_out_gaps()
     # use the constructor arg
+    empty = numpy.array([], dtype=int)
     m = IndelMap(
-        spans=tuple(m.spans), parent_length=m.parent_length, termini_unknown=True
+        gap_pos=empty.copy(),
+        gap_lengths=empty.copy(),
+        parent_length=6,
+        termini_unknown=True,
     )
     m_spans = tuple(m.spans)
-    assert not isinstance(m_spans[0], TerminalPadding)
+    assert len(m_spans) == 1 and not m_spans[0].lost
     # use the method
     m = IndelMap(
-        spans=tuple(m.spans), parent_length=m.parent_length
+        gap_pos=empty.copy(), gap_lengths=empty.copy(), parent_length=6
     ).with_termini_unknown()
     m_spans = tuple(m.spans)
-    assert not isinstance(m_spans[0], TerminalPadding)
+    assert len(m_spans) == 1 and not m_spans[0].lost
 
 
-@pytest.mark.parametrize("cls", (FeatureMap, IndelMap))
-def test_map_inverse(cls):
-    m = cls(locations=[(0, 2), (4, 6)], parent_length=6)
+def test_featuremap_inverse():
+    m = FeatureMap.from_locations(locations=[(0, 2), (4, 6)], parent_length=6)
     assert len(m) == 4
     mi = m.inverse()
     assert len(mi) == 6
@@ -479,6 +485,67 @@ def test_map_inverse(cls):
     assert got == expect
 
 
+def test_indelmap_from_aligned_segments():
+    locations = [(0, 2), (4, 6)]
+    im = IndelMap.from_aligned_segments(locations=locations, aligned_length=6)
+    assert len(im) == 6
+    expected_length = sum(e - s for s, e in locations)
+    assert im.parent_length == expected_length
+    im_spans = tuple(im.spans)
+    assert im_spans[1].lost and len(im_spans[1]) == 2
+
+
+def test_indelmap_from_aligned_segments2():
+    locations = [(0, 5), (7, 12), (14, 21), (24, 27)]
+    im = IndelMap.from_aligned_segments(locations=locations, aligned_length=27)
+    expected_length = sum(e - s for s, e in locations)
+    assert im.parent_length == expected_length
+
+
+def test_indelmap_merge():
+    im1 = IndelMap(
+        gap_pos=numpy.array([], dtype=int),
+        gap_lengths=numpy.array([], dtype=int),
+        parent_length=4,
+    )
+    im2 = IndelMap(
+        gap_pos=numpy.array([0], dtype=int),
+        gap_lengths=numpy.array([2], dtype=int),
+        parent_length=4,
+    )
+
+    assert len(im1) == 4
+    assert im1.parent_length == 4
+    inv = im1.merge_maps(im2)
+    assert inv.parent_length == 4
+    assert len(inv) == 6
+
+    fmap1 = im1.to_feature_map()
+    fmap2 = im2.to_feature_map()
+    got = fmap2.inverse()[fmap1.inverse()].inverse()
+    assert inv.get_gap_coordinates() == [
+        list(coords) for coords in got.get_gap_coordinates()
+    ]
+
+
+def test_indelmap_merge_parent_length():
+    im1 = IndelMap(
+        gap_pos=numpy.array([], dtype=int),
+        gap_lengths=numpy.array([], dtype=int),
+        parent_length=4,
+    )
+    im2 = IndelMap(
+        gap_pos=numpy.array([0], dtype=int),
+        gap_lengths=numpy.array([2], dtype=int),
+        parent_length=4,
+    )
+    # providing a value for parent_length overrides standard
+    im3 = im1.merge_maps(im2)
+    ov = im1.merge_maps(im2, parent_length=20)
+    assert ov.parent_length != im2.parent_length
+    assert ov.parent_length == 20
+
+
 @pytest.mark.parametrize("cls", (FeatureMap, IndelMap))
 def test_map_offsets(cls):
     # offsets are absolute starts of spans
@@ -490,9 +557,8 @@ def test_map_offsets(cls):
     assert got == [0, 1, 3, 6, 7, 8, 11]
 
 
-@pytest.mark.parametrize("cls", (FeatureMap, IndelMap))
-def test_map_indexed(cls):
-    m = cls(locations=[(0, 2), (4, 6)], parent_length=6).inverse()
+def test_map_indexed():
+    m = FeatureMap.from_locations(locations=[(0, 2), (4, 6)], parent_length=6).inverse()
     indexed = m[2]
     assert len(indexed) == 1
 
@@ -500,90 +566,474 @@ def test_map_indexed(cls):
 def test_compare_map_indexed():
     from cogent3.core.alignment import Aligned
 
-    seq = DNA.make_seq("--AC-GTAA--".replace("-", ""))
-    spans = [LostSpan(2), Span(0, 2), LostSpan(2), Span(2, 6), LostSpan(2)]
-    kwargs = dict(spans=spans, parent_length=len(seq))
-    mm = FeatureMap(**kwargs)
-    im = IndelMap(**kwargs)
-    ma = Aligned(mm, seq)
+    raw_seq = "--AC-GTAA--"
+    im, seq = DNA.make_seq(raw_seq).parse_out_gaps()
     ia = Aligned(im, seq)
-    first = ia[0]
-    assert first == "-"
-    length = len(ma)
+    length = len(raw_seq)
     got = [str(ia[i]) for i in range(length)]
-    expect = [str(ma[i]) for i in range(length)]
+    expect = list("--AC-GTAA--")
     assert got == expect
 
 
 @pytest.mark.parametrize("slice_it", (True, False))
-def test_indel_map_zeroed(slice_it):
+def test_feature_map_zeroed(slice_it):
     spans = [LostSpan(2), Span(2, 4), LostSpan(2), Span(4, 8), LostSpan(2)]
     kwargs = dict(spans=spans, parent_length=6)
 
-    mm = FeatureMap(**kwargs)
-    if slice_it:
-        mm = mm[:6]
-    mm_zeroed = mm.zeroed()
-
-    im = IndelMap(**kwargs)
+    fm = FeatureMap(**kwargs)
     if slice_it:
-        im = im[:6]
+        fm = fm[3:6]
 
-    im_zeroed = im.zeroed()
-    assert im_zeroed.get_coordinates() == mm_zeroed.get_coordinates()
-    assert im_zeroed.parent_length == mm_zeroed.parent_length
+    fm_zeroed = fm.zeroed()
+    assert fm.start > 0
+    assert fm_zeroed.start == 0
 
 
 def test_indelmap_to_feature_map():
     spans = [LostSpan(2), Span(2, 4), LostSpan(2), Span(4, 8), LostSpan(2)]
-    kwargs = dict(spans=spans, parent_length=6)
-    im = IndelMap(**kwargs)
+    kwargs = dict(spans=spans, parent_length=8)
+    im = IndelMap.from_spans(**kwargs)
     mm = im.to_feature_map()
     assert mm.get_coordinates() == im.get_coordinates()
 
 
 def test_indelmap_nucleic_reversed():
-    spans = [LostSpan(2), Span(2, 4), LostSpan(2), Span(4, 8), LostSpan(2)]
-    kwargs = dict(spans=spans, parent_length=12)
-    orig = IndelMap(**kwargs)
+    #   01  2345
+    # --AC--GGGG--
+    orig, s = DNA.make_seq("--AC--GGGG--").parse_out_gaps()
+    plus_coords = [(0, 2), (2, 6)]
+    assert orig.get_coordinates() == plus_coords
     rev = orig.nucleic_reversed()
-    rev_spans = tuple(rev.spans)
-    assert rev_spans[1].reverse == rev_spans[3].reverse == False
-    old = orig.nucleic_reversed()
-    old_spans = tuple(old.spans)
-    assert rev.get_coordinates() == old.get_coordinates()
-
-
-def test_indelmap_with_reverse_span():
-    spans = [
-        LostSpan(2),
-        Span(8, 4, reverse=True),
-        LostSpan(2),
-        Span(4, 2, reverse=True),
-        LostSpan(2),
-    ]
-    imap = IndelMap(spans=spans, parent_length=12)
-    assert (imap.start, imap.end) == (8, 2)
-
-
-def test_indelmap_no_gaps():
-    imap = IndelMap(locations=(), parent_length=6)
-    gaps = imap.gaps()
-    assert not gaps
+    assert rev.num_gaps == orig.num_gaps
+    minus_coords = [(0, 4), (4, 6)]
+    assert rev.get_coordinates() == minus_coords
+    # reversing again returns original
+    back = rev.nucleic_reversed()
+    assert orig.get_coordinates() == back.get_coordinates()
 
 
 def test_get_coords():
     """get_coordinates should return raw coordinates matching input"""
     spans = [(0, 9), (20, 32)]
-    fmap = FeatureMap(locations=spans, parent_length=100)
+    fmap = FeatureMap.from_locations(locations=spans, parent_length=100)
     coords = fmap.get_coordinates()
     assert coords == spans
 
 
+def test_indel_map_get_coords():
+    """get_coordinates should return raw coordinates matching input"""
+    imap = IndelMap(
+        gap_pos=numpy.array([9]), gap_lengths=numpy.array([10]), parent_length=20
+    )
+    locations = [(0, 9), (9, 20)]
+    coords = imap.get_coordinates()
+    assert coords == locations
+
+
+def test_indel_map_get_gap_coords():
+    gap_data = numpy.array([(2, 3), (3, 1), (6, 2)])
+    gap_pos, lengths = gap_data[:, 0], gap_data[:, 1]
+    m = IndelMap(gap_pos=gap_pos, gap_lengths=lengths, parent_length=10)
+    got = m.get_gap_coordinates()
+    assert got == gap_data.tolist()
+
+
+@pytest.mark.parametrize(
+    "coords", ([(0, 3), (7, 11)], [(0, 3)], [(2, 4), (6, 10)], [(4, 6)])
+)
+def test_indelmap_joined_segments(coords):
+    raw = "--AC--GGGG--"
+    expect, _ = DNA.make_seq("".join(raw[s:e] for s, e in coords)).parse_out_gaps()
+    imap, s = DNA.make_seq(raw).parse_out_gaps()
+    got = imap.joined_segments(coords)
+    assert got.gap_pos.tolist() == expect.gap_pos.tolist()
+    assert got.cum_gap_lengths.tolist() == expect.cum_gap_lengths.tolist()
+
+
+def test_indelmap_slice_terminating():
+    raw = "-CB-A--"
+    start, end = 4, 6
+    expect, _ = DNA.make_seq(raw[start:end]).parse_out_gaps()
+    imap, s = DNA.make_seq(raw).parse_out_gaps()
+    got = imap[start:end]
+    assert got.gap_pos.tolist() == expect.gap_pos.tolist()
+    assert got.cum_gap_lengths.tolist() == expect.cum_gap_lengths.tolist()
+
+
+def test_indelmap_slice_innner_gap():
+    start, end = 4, 6
+    raw = "--AC--GGGG--"
+    expect, _ = DNA.make_seq(raw[start:end]).parse_out_gaps()
+    imap, _ = DNA.make_seq(raw).parse_out_gaps()
+    imap = imap[start:end]
+    assert imap.gap_pos.tolist() == expect.gap_pos.tolist()
+
+
+def test_indelmap_slice_cum_length():
+    start, end = 7, 11
+    raw = "--AC--GGGG--"
+    expect, _ = DNA.make_seq(raw[start:end]).parse_out_gaps()
+    imap, _ = DNA.make_seq(raw).parse_out_gaps()
+    imap = imap[start:end]
+    assert imap.gap_pos.tolist() == expect.gap_pos.tolist()
+    assert imap.cum_gap_lengths.tolist() == expect.cum_gap_lengths.tolist()
+
+
 def test_get_coords_invalid_order():
     """get_coordinates should return raw coordinates matching input"""
 
     # should work for reversed Maps too
     spans = [(32, 20), (9, 0)]
     with pytest.raises(ValueError):
-        FeatureMap(locations=spans, parent_length=100)
+        FeatureMap.from_locations(locations=spans, parent_length=100)
+
+
+def test_indelmap_slice_gap_into_seq():
+    pos, lengths = numpy.array([[3, 1], [7, 1]]).T
+    gaps = IndelMap(gap_pos=pos, gap_lengths=lengths, parent_length=7)
+    # example seq
+    # 012345678
+    # 012 3456
+    # ACA-CGAC-
+    # slicing from aligned coords 3-5, gives '-C'
+    got = gaps[3:5]
+    assert got.gap_pos.tolist() == [0]
+    assert got.cum_gap_lengths.tolist() == [1]
+
+
+def test_indelmap_slice_one_gap():
+    pos, lengths = numpy.array([[3, 6]]).T
+    gaps = IndelMap(gap_pos=pos, gap_lengths=lengths, parent_length=24)
+    # example seq
+    #           1
+    # 01234567890
+    # 012      3456789
+    # GTG------GTAGAAGTTCCAAATAATGAA
+    # slicing from aligned coords 3-5, gives 'TG------G'
+    # in sliced seq, gap at 2
+    got = gaps[1:10]
+    assert got.gap_pos.tolist() == [2]
+    assert got.cum_gap_lengths.tolist() == [6]
+
+
+@pytest.mark.parametrize(
+    "data",
+    (
+        "AB---CD--EF",
+        "---ABCD--EF",
+        "ABCD---EF--",
+        "-----ABCDEF",
+        "ABCDEF-----",
+        "-ABCDEF----",
+        "-A-B-C-D-EF",
+        "A-B-C-D-EF-",
+    ),
+)
+@pytest.mark.parametrize("index", range(4, 6))  # the ungapped sequence is 6 long
+def test_gapped_convert_seq2aln(data, index):
+    # converting a sequence index to alignment index
+    ungapped = data.replace("-", "")
+    seq = make_seq(data, moltype="text")
+    gaps, s = seq.parse_out_gaps()
+    idx = gaps.get_align_index(index)
+    assert data[idx] == ungapped[index]
+
+
+@pytest.mark.parametrize(
+    "data",
+    (
+        "AB---CD--EF",
+        "---ABCD--EF",
+        "ABCD---EF--",
+        "-----ABCDEF",
+        "ABCDEF-----",
+        "-ABCDEF----",
+        "-A-B-C-D-EF",
+        "A-B-C-D-EF-",
+    ),
+)
+@pytest.mark.parametrize(
+    "start,end", list(combinations(range(6), 2))
+)  # the ungapped sequence is 6 long
+def test_indelmap_align_index_slice_stop(data, start, end):
+    # converting a sequence index to alignment index
+    ungapped = data.replace("-", "")
+    seq = make_seq(data, moltype="text")
+    gaps, s = seq.parse_out_gaps()
+    stop = gaps.get_align_index(end, slice_stop=True)
+    assert data[stop - 1] == ungapped[end - 1]
+
+
+@pytest.mark.parametrize(
+    "data",
+    (
+        "AB---CD--EF",
+        "---ABCD--EF",
+        "ABCD---EF--",
+        "-----ABCDEF",
+        "ABCDEF-----",
+        "-ABCDEF----",
+        "-A-B-C-D-EF",
+        "A-B-C-D-EF-",
+    ),
+)
+@pytest.mark.parametrize("index", range(6))  # the ungapped sequence is 6 long
+def test_gapped_convert_seq2aln2seq(data, index):
+    # round tripping seq to alignment to seq indices
+    gaps, s = make_seq(data, moltype="text").parse_out_gaps()
+    align_index = gaps.get_align_index(index)
+    got = gaps.get_seq_index(align_index)
+    assert got == index
+
+
+@pytest.mark.parametrize("expect", range(10))
+def test_indelmap_get_seq_index_negative(expect):
+    parent_length = 10
+    gap_pos = [0]
+    gap_lengths = [3]
+    imap = IndelMap(
+        gap_pos=numpy.array(gap_pos),
+        gap_lengths=numpy.array(gap_lengths),
+        parent_length=parent_length,
+    )
+    neg_index = expect - parent_length
+    got = imap.get_seq_index(neg_index)
+    assert got == expect
+
+
+@pytest.mark.parametrize("expect", range(10))
+def test_indelmap_get_align_index_negative(expect):
+    parent_length = 10
+    gap_pos = [0]
+    gap_lengths = [3]
+    gap_length = gap_lengths[0]
+    imap = IndelMap(
+        gap_pos=numpy.array(gap_pos),
+        gap_lengths=numpy.array(gap_lengths),
+        parent_length=parent_length,
+    )
+    neg_index = expect + gap_length - len(imap)
+    got = imap.get_seq_index(neg_index)
+    assert got == expect
+
+
+@pytest.mark.parametrize(
+    "data",
+    (
+        "AB--CDE-FG",
+        "--ABC-DEFG",
+        "AB--CDE-FG--",
+        "ABCDE--FG---",
+        "-----ABCDEFG",
+        "-A-B-C-D-E-F-G-",
+    ),
+)
+@pytest.mark.parametrize("seq_index", range(7))
+def test_gapped_convert_aln2seq_nongap_char(data, seq_index):
+    # test alignment indexes when aligned position is NOT a gap
+    ungapped = "ABCDEFG"
+    align_index = data.find(ungapped[seq_index])
+    gaps, seq = make_seq(data, moltype="text").parse_out_gaps()
+    idx = gaps.get_seq_index(align_index)
+    assert idx == seq_index
+
+
+def _find_nth_gap_index(data: str, n: int) -> int:
+    num = -1
+    for i, c in enumerate(data):
+        if c == "-":
+            num += 1
+        if num == n:
+            return i
+    raise ValueError(f"{data=}, {n=}")
+
+
+def _get_expected_seqindex(data: str, align_index: int) -> int:
+    # compute the expected seqindex
+    refseq = data.replace("-", "")
+    got = data[align_index:].lstrip("-")
+    return refseq.find(got[0]) if got else len(refseq)
+
+
+@pytest.mark.parametrize(
+    "data",
+    (
+        "AB-----CDE-F--G",
+        "----ABC-DEFG---",
+        "AB--CDE-FG-----",
+        "ABCDE--FG------",
+        "--------ABCDEFG",
+        "-A-B-C-D-E-F-G-",
+    ),
+)
+@pytest.mark.parametrize("gap_number", range(8))
+def test_gapped_convert_aln2seq_gapchar(data, gap_number):
+    # test alignment indexes when aligned position IS a gap
+    # in this case we expect the position of the next non-gap
+    # to be the result
+    # find the alignment index corresponding to the
+    align_index = _find_nth_gap_index(data, gap_number)
+    assert data[align_index] == "-", (data, gap_number)
+    # find nearest non-gap
+    seq_index = _get_expected_seqindex(data, align_index)
+    gaps, seq = make_seq(data, moltype="text").parse_out_gaps()
+    idx = gaps.get_seq_index(align_index)
+    assert idx == seq_index
+
+
+def test_gapped_convert_aln2seq_invalid():
+    gaps, seq = make_seq("AC--GTA-TG", moltype="dna").parse_out_gaps()
+    with pytest.raises(NotImplementedError):
+        # absolute value of negative indices must be < seq length
+        gaps.get_seq_index(-100)
+
+
+@pytest.mark.parametrize(
+    "invalid_slice",
+    (slice(None, None, -1), slice(None, None, 2)),
+)
+def test_gap_pos_invalid_slice(invalid_slice):
+    pos, lengths = numpy.array([[1, 3]], dtype=numpy.int32).T
+    gp = IndelMap(gap_pos=pos, gap_lengths=lengths, parent_length=20)
+    with pytest.raises(NotImplementedError):
+        gp[invalid_slice]
+
+
+@pytest.mark.parametrize(
+    "slice",
+    (
+        slice(3, 7),
+        slice(20, None),
+    ),
+)
+def test_no_gaps_in_slice(slice):
+    # aligned length is 25
+    seq_length = 20
+    gap_length = 5
+    pos, lengths = numpy.array([[10, gap_length]], dtype=numpy.int32).T
+    gp = IndelMap(gap_pos=pos, gap_lengths=lengths, parent_length=seq_length)
+    got = gp[slice]
+    assert not got.num_gaps
+    start = slice.start or 0
+    stop = slice.stop or (seq_length + gap_length)
+    assert len(got) == stop - start
+
+
+def test_len_gapped():
+    seq_length = 20
+    gap_length = 5
+    pos, lengths = numpy.array([[10, gap_length]], dtype=numpy.int32).T
+    gp = IndelMap(gap_pos=pos, gap_lengths=lengths, parent_length=seq_length)
+    assert len(gp) == (seq_length + gap_length)
+
+
+def test_all_gaps_in_slice():
+    # slicing GapPositions
+    # sample seq 1
+    data = "AC--GTA-TG"
+    gp, seq = make_seq(data, moltype="dna").parse_out_gaps()
+    sl = slice(1, 9)
+
+    got = gp[sl]
+    expect_gaps, expect_seq = make_seq(data[sl], moltype="dna").parse_out_gaps()
+    assert got.get_gap_coordinates() == expect_gaps.get_gap_coordinates()
+    assert got.parent_length == 5
+
+
+@pytest.mark.parametrize(
+    "data",
+    (
+        "----GTA-TG",
+        "AC--GTA---",
+        "AC--GTA-TG",
+        "A-C-G-T-A-",
+        "-A-C-G-T-A",
+        "ACGTAACGTA",
+        "----------",
+    ),
+)
+@pytest.mark.parametrize(
+    "slice",
+    [slice(i, j) for i, j in combinations(range(10), 2)],
+)
+def test_variant_slices(data, slice):
+    gaps, _ = make_seq(data, moltype="dna").parse_out_gaps()
+    got = gaps[slice]
+
+    expect_gaps, expect_seq = make_seq(data[slice], moltype="dna").parse_out_gaps()
+    assert got.parent_length == len(expect_seq)
+    assert (got.gap_pos == expect_gaps.gap_pos).all()
+    assert (got.cum_gap_lengths == expect_gaps.cum_gap_lengths).all()
+
+
+def test_indelmap_joined():
+    pos = numpy.array([0, 1])
+    cum_len = numpy.array([1, 5])
+    imap = IndelMap(gap_pos=pos, cum_gap_lengths=cum_len, parent_length=5)
+    fmap = FeatureMap.from_locations(locations=[(0, 1), (2, 5)], parent_length=10)
+    coords = fmap.get_coordinates()
+    got = imap.joined_segments(coords)
+    assert got.num_gaps == 1
+    assert got.gap_pos[0] == 0
+    assert got.cum_gap_lengths[0] == (1 + 5 - 2)
+
+
+def test_indel_map_sliced_with_negative():
+    imap = IndelMap(
+        gap_pos=numpy.array([0]), cum_gap_lengths=numpy.array([1]), parent_length=14
+    )
+    got = imap[:-3]
+    assert got.parent_length == 14 - 3
+
+
+def test_indelmap_roundtrip_json():
+    from cogent3.util.deserialise import deserialise_object
+
+    imap = IndelMap(
+        gap_pos=numpy.array([0, 9]),
+        cum_gap_lengths=numpy.array([1, 3]),
+        parent_length=14,
+    )
+    got = deserialise_object(imap.to_json())
+    assert (got.gap_pos == imap.gap_pos).all()
+    assert (got.cum_gap_lengths == imap.cum_gap_lengths).all()
+    assert got.parent_length == imap.parent_length
+
+
+def test_featuremap_roundtrip_json():
+    from cogent3.util.deserialise import deserialise_object
+
+    fmap = FeatureMap.from_locations(
+        locations=[[0, 9], [20, 23]],
+        parent_length=140,
+    )
+    got = deserialise_object(fmap.to_json())
+    coords = fmap.get_coordinates()
+    assert coords == [(0, 9), (20, 23)]
+    assert got.parent_length == fmap.parent_length == 140
+
+
+def test_gap_coords_to_map():
+    """construct a Map from coordinates of gap alone"""
+    gap_coords = {0: 1, 2: 2, 4: 1, 7: 2}
+    seqlen = 70
+    got = gap_coords_to_map(gap_coords, seqlen)
+    assert 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)
+    assert len(got) == sum(gap_coords.values()) + seqlen
+
+    # roundtrip from Map.get_gap_coordinates()
+    assert 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))
+    assert len(got) == len(m)
+    assert got.get_coordinates() == m.get_coordinates()
+
+    # and gaps outside sequence
+    with pytest.raises(ValueError):
+        gap_coords_to_map({20: 1}, len(seq))
diff --git a/tests/test_core/test_maps.py b/tests/test_core/test_maps.py
index 54c3e4829..5eca67498 100644
--- a/tests/test_core/test_maps.py
+++ b/tests/test_core/test_maps.py
@@ -9,7 +9,7 @@ class MapTest(unittest.TestCase):
 
     def test_spans(self):
         # a simple two part map of length 10
-        map = FeatureMap(locations=[(0, 5), (5, 10)], parent_length=10)
+        map = FeatureMap.from_locations(locations=[(0, 5), (5, 10)], parent_length=10)
         # try different spans on the above map
         for (start, end), expected in [
             ((0, 4), "[0:4]"),
diff --git a/tests/test_draw/test_dotplot.py b/tests/test_draw/test_dotplot.py
index 38c8d6a63..19180a7c0 100644
--- a/tests/test_draw/test_dotplot.py
+++ b/tests/test_draw/test_dotplot.py
@@ -6,7 +6,6 @@ from cogent3.draw.dotplot import (
     Dotplot,
     _convert_input,
     get_align_coords,
-    len_seq,
     not_gap,
 )
 
@@ -15,7 +14,7 @@ class TestUtilFunctions(TestCase):
     def test_len_seq(self):
         """returns length of sequence minus gaps"""
         m, seq = DNA.make_seq("ACGGT--A").parse_out_gaps()
-        self.assertEqual(len_seq(m), 6)
+        self.assertEqual(m.parent_length, 6)
 
     def test_not_gap(self):
         """distinguishes Map instances that include gap or not"""
diff --git a/tests/test_parse/test_cigar.py b/tests/test_parse/test_cigar.py
index 5f09808ff..e8547afb4 100644
--- a/tests/test_parse/test_cigar.py
+++ b/tests/test_parse/test_cigar.py
@@ -49,7 +49,7 @@ class TestCigar(unittest.TestCase):
                 slicealn1 = self.seq[loc1[0] : loc1[1]].gapped_by_map(map1)
                 assert ori1 == slicealn1
             else:
-                assert map1.length == len(ori1)
+                assert len(map1) == len(ori1)
 
             # test by_align = False
             map2, loc2 = slice_cigar(self.cigar_text, start, end, by_align=False)
@@ -84,7 +84,7 @@ class TestCigar(unittest.TestCase):
                 start=start,
                 end=end,
             )
-            assert cmp_aln.to_dict() == slice_aln.to_dict()
+            assert cmp_aln.to_dict() == slice_aln.to_dict(), (start, end)
 
 
 if __name__ == "__main__":