Repository: cogent3/cogent3
Pull request #1367 (2023-04-17 04:38:27)

Annotation refactor


Patch:
diff --git a/src/cogent3/__init__.py b/src/cogent3/__init__.py
index e7d00edbe..3de4e33b5 100644
--- a/src/cogent3/__init__.py
+++ b/src/cogent3/__init__.py
@@ -17,7 +17,9 @@ from cogent3.core.alignment import (
     Sequence,
     SequenceCollection,
 )
+from cogent3.core.annotation_db import load_annotations
 from cogent3.core.genetic_code import available_codes, get_code
+
 # note that moltype has to be imported last, because it sets the moltype in
 # the objects created by the other modules.
 from cogent3.core.moltype import (
@@ -25,8 +27,6 @@ from cogent3.core.moltype import (
     DNA,
     PROTEIN,
     RNA,
-    STANDARD_CODON,
-    CodonAlphabet,
     available_moltypes,
     get_moltype,
 )
@@ -234,6 +234,66 @@ def _load_seqs(file_format, filename, fmt, kw, parser_kw):
 
 
 def load_seq(
+    filename: Union[str, pathlib.Path],
+    annotation_path=None,
+    format: Optional[str] = None,
+    moltype: Optional[str] = None,
+    label_to_name: Optional[Callable] = None,
+    parser_kw: dict = None,
+    info: dict = None,
+    **kw,
+) -> Sequence:
+    """
+    loads unaligned sequences from file
+
+    Parameters
+    ----------
+    filename : str
+        path to sequence file
+    format : str
+        sequence file format, if not specified tries to guess from the path suffix
+    moltype : str
+        the moltype, eg DNA, PROTEIN, 'dna', 'protein'
+    label_to_name : callable
+        function for converting original name into another name.
+    parser_kw : dict
+        optional arguments for the parser
+    info : dict
+        a dict from which to make an info object
+    **kw
+        other keyword arguments passed to SequenceCollection
+
+    Notes
+    -----
+    Returns **one** sequence from a file. Use load_aligned_seqs or
+    load_unaligned_seqs to get a collection.
+
+    Returns
+    -------
+    ``Sequence``
+    """
+    info = info or {}
+    info["source"] = str(filename)
+    file_format, _ = get_format_suffixes(filename)
+    if file_format == "json":
+        seq = load_from_json(filename, (Sequence,))
+        seq.name = label_to_name(seq.name) if label_to_name else seq.name
+        return seq
+
+    data = _load_seqs(file_format, filename, format, kw, parser_kw)
+    name, seq = data[0]
+    name = label_to_name(name) if label_to_name else name
+    result = make_seq(seq, name, moltype=moltype)
+    result.info.update(info)
+
+    if "gb" in file_format:
+        result.annotation_db = load_annotations(filename, seqids=[name])
+    elif annotation_path is not None:
+        result.annotation_db = load_annotations(annotation_path, seqids=[name])
+    return result
+
+
+def load_seq_old(
     filename: Union[str, pathlib.Path],
     format: Optional[str] = None,
     moltype: Optional[str] = None,
diff --git a/src/cogent3/core/alignment.py b/src/cogent3/core/alignment.py
index f17f59bfd..5aa546167 100644
--- a/src/cogent3/core/alignment.py
+++ b/src/cogent3/core/alignment.py
@@ -52,6 +52,7 @@ from numpy.random import choice, permutation, randint
 import cogent3  # will use to get at cogent3.parse.fasta.MinimalFastaParser,
 
 from cogent3.core.annotation import Map, _Annotatable
+from cogent3.core.annotation_db import GenbankAnnotationDb, load_annotations
 from cogent3.core.genetic_code import get_code
 from cogent3.core.info import Info as InfoClass
 from cogent3.core.profile import PSSM, MotifCountsArray
@@ -509,6 +510,19 @@ class _SequenceCollectionBase:
         self._set_additional_attributes(curr_seqs)
 
         self._repr_policy = dict(num_seqs=10, num_pos=60, ref_name="longest", wrap=60)
+        self._annotation_db = None
+
+    @property
+    def annotation_db(self):
+        return self._annotation_db
+
+    @annotation_db.setter
+    def annotation_db(self, value):
+        from cogent3.core.annotation_db import SupportsFeatures
+
+        if not isinstance(value, SupportsFeatures):
+            raise TypeError
+        self._annotation_db = value
 
     def __str__(self):
         """Returns self in FASTA-format, respecting name order."""
@@ -1588,8 +1602,7 @@ class _SequenceCollectionBase:
 
     def to_dna(self):
         """returns copy of self as an alignment of DNA moltype seqs"""
-        make_seq = cogent3.DNA.make_seq
-        data = [make_seq(s, s.name) for s in self.seqs]
+        data = [s.to_moltype("dna") for s in self.seqs]
         new = self.__class__(
             data=data, moltype=cogent3.DNA, name=self.name, info=self.info
         )
@@ -1599,8 +1612,7 @@ class _SequenceCollectionBase:
 
     def to_rna(self):
         """returns copy of self as an alignment of RNA moltype seqs"""
-        make_seq = cogent3.RNA.make_seq
-        data = [make_seq(s, s.name) for s in self.seqs]
+        data = [s.to_moltype("rna") for s in self.seqs]
         new = self.__class__(
             data=data, moltype=cogent3.RNA, name=self.name, info=self.info
         )
@@ -1610,8 +1622,7 @@ class _SequenceCollectionBase:
 
     def to_protein(self):
         """returns copy of self as an alignment of PROTEIN moltype seqs"""
-        make_seq = cogent3.PROTEIN.make_seq
-        data = [make_seq(s, s.name) for s in self.seqs]
+        data = [s.to_moltype("protein") for s in self.seqs]
         new = self.__class__(
             data=data, moltype=cogent3.PROTEIN, name=self.name, info=self.info
         )
@@ -1630,7 +1641,7 @@ class _SequenceCollectionBase:
         return rc
 
     def reverse_complement(self):
-        """Returns the reverse complement alignment. A synonymn for rc."""
+        """Returns the reverse complement alignment. A synonym for rc."""
         return self.rc()
 
     def pad_seqs(self, pad_length=None, **kwargs):
@@ -1999,30 +2010,29 @@ class SequenceCollection(_SequenceCollectionBase):
             if name in self.named_seqs:
                 self.named_seqs[name].copy_annotations(seq)
 
-    def annotate_from_gff(self, f):
-        """Copies annotations from gff-format file to self.
+    def annotate_from_gff(self, f: os.PathLike, seq_ids: [str]):
+        """copies annotations from a gff file to a sequence in self
 
-        Matches by name of sequence. This method accepts string path
-        or pathlib.Path or file-like object (e.g. StringIO)
+        Parameters
+        ----------
+        f : path to gff annotation file.
+        seq_name : names of seqs to be annotated.
+        does not support setting offset, set offset directly on sequences with seq.annotation_offset = offset
 
-        Skips sequences in the file that are not in self.
         """
+        if isinstance(seq_ids, str):
+            seq_ids = [seq_ids]
 
-        seq_dict = {}
+        if self.annotation_db is None:
+            self.annotation_db = load_annotations(f, seq_ids)
+        elif isinstance(self.annotation_db, GenbankAnnotationDb):
+            raise ValueError("GenbankAnnotationDb already attached")
+        else:
+            self.annotation_db = load_annotations(f, seq_ids, db=self.annotation_db._db)
 
-        for gff_dict in gff_parser(f):
-            if gff_dict["SeqID"] not in self.named_seqs:
-                continue
-            seq_id = gff_dict["SeqID"]
-            if seq_id not in seq_dict.keys():
-                seq_dict[seq_id] = [gff_dict]
-                continue
-            seq_dict[seq_id].append(gff_dict)
-        for seq_id in seq_dict.keys():
-            seq = self.named_seqs[seq_id]
-            if not hasattr(seq, "annotations"):
-                seq = seq.data
-            seq.annotate_from_gff(seq_dict[seq_id], pre_parsed=True)
+        # add reference to this annotation_db on each sequence that was provided in seq_ids
+        for seq in seq_ids:
+            self.get_seq(seq).annotation_db = self.annotation_db
 
 
 @total_ordering
@@ -2179,8 +2189,24 @@ class Aligned(object):
         return Aligned(map[self.map.inverse()].inverse(), self.data)
 
     def get_annotations_matching(self, alignment, annotation_type="*", **kwargs):
-        for annot in self.data.get_annotations_matching(
-            annotation_type=annotation_type, **kwargs
+
+        from cogent3.util.warning import deprecated
+
+        deprecated(
+            "method",
+            "get_annotations_matching",
+            "get_features_matching",
+            "2023.3",
+        )
+
+        for annot in self.data.get_features_matching(
+            feature_type=annotation_type, **kwargs
+        ):
+            yield annot.remapped_to(alignment, self.map.inverse())
+
+    def get_features_matching(self, alignment, annotation_type="*", **kwargs):
+        for annot in self.data.get_features_matching(
+            feature_type=annotation_type, **kwargs
         ):
             yield annot.remapped_to(alignment, self.map.inverse())
 
@@ -2313,12 +2339,12 @@ class AlignmentI(object):
             col_lookup = dict.fromkeys(cols)
             for name, seq in list(self.named_seqs.items()):
                 result[name] = make_seq(
-                    [seq[i] for i in range(len(seq)) if i not in col_lookup]
+                    [str(seq[i]) for i in range(len(seq)) if i not in col_lookup]
                 )
         # otherwise, just get the requested indices
         else:
             for name, seq in list(self.named_seqs.items()):
-                result[name] = make_seq([seq[i] for i in cols])
+                result[name] = make_seq("".join(str(seq[i]) for i in cols))
         return self.__class__(result, names=self.names, info=self.info)
 
     def get_position_indices(self, f, native=False, negate=False):
@@ -4389,11 +4415,24 @@ class Alignment(_Annotatable, AlignmentI, SequenceCollection):
                 aligned_seqs.append(self._seq_to_aligned(s, n))
         self.named_seqs = dict(list(zip(names, aligned_seqs)))
         self.seq_data = self._seqs = aligned_seqs
+        self._annotation_db = None
+
+    @property
+    def annotation_db(self):
+        return self._annotation_db
+
+    @annotation_db.setter
+    def annotation_db(self, value):
+        from cogent3.core.annotation_db import SupportsFeatures
+
+        if not isinstance(value, SupportsFeatures):
+            raise TypeError
+        self._annotation_db = value
 
     def _coerce_seqs(self, seqs, is_array):
-        if not min(
-            [isinstance(seq, _Annotatable) or isinstance(seq, Aligned) for seq in seqs]
-        ):
+        if any(isinstance(seq, ArraySequence) for seq in seqs):
+            seqs = [self.moltype.make_seq(str(seq), name=seq.name) for seq in seqs]
+        elif not any(isinstance(seq, (_Annotatable, Aligned)) for seq in seqs):
             seqs = list(map(self.moltype.make_seq, seqs))
         return seqs
 
@@ -4404,6 +4443,27 @@ class Alignment(_Annotatable, AlignmentI, SequenceCollection):
         ).parse_out_gaps()
         return Aligned(map, seq)
 
+    def __getitem__(self, index):
+
+        if hasattr(index, "get_slice"):
+            return index.get_slice()
+
+        if isinstance(index, Map):
+            new = self._mapped(index)
+
+        elif isinstance(index, (int, slice)):
+            seqs = [s[index] for s in self.seqs]
+
+            new = self.__class__(seqs, name=self.name, info=self.info)
+
+        if self.annotation_db is not None:
+            new.annotation_db = self.annotation_db
+
+        if hasattr(self, "_repr_policy"):
+            new._repr_policy.update(self._repr_policy)
+
+        return new
+
     def __repr__(self):
         seqs = []
         limit = 10
@@ -4443,12 +4503,12 @@ class Alignment(_Annotatable, AlignmentI, SequenceCollection):
         return annot.remapped_to(target_aligned.data, target_aligned.map)
 
     def get_projected_annotations(self, seq_name, *args):
-        aln_annots = self.get_annotations_matching(*args)
+        aln_annots = self.get_features_matching(*args)
         return [self.project_annotation(seq_name, a) for a in aln_annots]
 
     def get_annotations_from_seq(self, seq_name, annotation_type="*", **kwargs):
         aligned = self.named_seqs[seq_name]
-        return aligned.get_annotations_matching(
+        return aligned.get_features_matching(
             self, annotation_type=annotation_type, **kwargs
         )
 
@@ -4588,7 +4648,7 @@ class Alignment(_Annotatable, AlignmentI, SequenceCollection):
                 if g == tgp:
                     combo.append(gap)
                 else:
-                    combo.append(s)
+                    combo.append(str(s))
             result[name] = combo
         return Alignment(result, alphabet=self.alphabet.with_gap_motif())
 
diff --git a/src/cogent3/core/annotation.py b/src/cogent3/core/annotation.py
index b42ab2470..b5abac1cf 100644
--- a/src/cogent3/core/annotation.py
+++ b/src/cogent3/core/annotation.py
@@ -13,7 +13,7 @@ from .location import Map, as_map
 
 __author__ = "Peter Maxwell and Gavin Huttley"
 __copyright__ = "Copyright 2007-2022, The Cogent Project"
-__credits__ = ["Peter Maxwell", "Gavin Huttley"]
+__credits__ = ["Peter Maxwell", "Gavin Huttley", "Katherine Caley"]
 __license__ = "BSD-3"
 __version__ = "2023.2.12a1"
 __maintainer__ = "Gavin Huttley"
@@ -21,7 +21,104 @@ __email__ = "gavin.huttley@anu.edu.au"
 __status__ = "Production"
 
 
-class _Annotatable:
+class _AnnotationMixin:
+    def add_feature(self, biotype, name, spans):
+        feat = Feature(self, biotype, name, spans)
+
+        if feat.parent is not self:
+            raise ValueError("doesn't belong here")
+
+        if feat.attached:
+            raise ValueError("already attached")
+
+        if self.annotations is self.__class__.annotations:
+            self.annotations = []
+
+        self.annotations.extend([feat])
+
+        feat.attached = True
+
+        return feat
+
+    def get_features_matching(self, feature_type, name=None, extend_query=False):
+        """
+
+        Parameters
+        ----------
+        feature_type : string
+            name of the feature type. Wild-cards allowed.
+        name : string
+            name of the instance. Wild-cards allowed.
+        extend_query : boolean
+            queries sub-annotations if True
+        Returns
+        -------
+        list of AnnotatableFeatures
+        """
+        result = []
+        if len(self.annotations) == 0:
+            return result
+        for annotation in self.annotations:
+            if fnmatch(annotation.type, feature_type) and (
+                name is None or fnmatch(annotation.name, name)
+            ):
+                result.append(annotation)
+            if extend_query:
+                result.extend(
+                    annotation.get_features_matching(
+                        feature_type, name, extend_query=extend_query
+                    )
+                )
+        return result
+
+    def get_drawables(self):
+        """returns a dict of drawables, keyed by type"""
+        result = defaultdict(list)
+        for a in self.annotations:
+            result[a.type].append(a.get_drawable())
+        return result
+
+    def get_drawable(self, width=600, vertical=False):
+        """returns Drawable instance"""
+        from cogent3.draw.drawable import Drawable
+
+        drawables = self.get_drawables()
+        if not drawables:
+            return None
+        # we order by tracks
+        top = 0
+        space = 0.25
+        annotes = []
+        for feature_type in drawables:
+            new_bottom = top + space
+            for i, annott in enumerate(drawables[feature_type]):
+                annott.shift(y=new_bottom - annott.bottom)
+                if i > 0:
+                    annott._showlegend = False
+                annotes.append(annott)
+
+            top = annott.top
+
+        top += space
+        height = max((top / len(self)) * width, 300)
+        xaxis = dict(range=[0, len(self)], zeroline=False, showline=True)
+        yaxis = dict(range=[0, top], visible=False, zeroline=True, showline=True)
+
+        if vertical:
+            all_traces = [t.T.as_trace() for t in annotes]
+            width, height = height, width
+            xaxis, yaxis = yaxis, xaxis
+        else:
+            all_traces = [t.as_trace() for t in annotes]
+
+        drawer = Drawable(
+            title=self.name, traces=all_traces, width=width, height=height
+        )
+        drawer.layout.update(xaxis=xaxis, yaxis=yaxis)
+        return drawer
+
+
+class _Annotatable(_AnnotationMixin):
     # default
     annotations = ()
 
@@ -96,61 +193,54 @@ class _Annotatable:
     def _mapped(self, map):
         raise NotImplementedError
 
-    def get_drawables(self):
-        """returns a dict of drawables, keyed by type"""
-        result = defaultdict(list)
-        for a in self.annotations:
-            result[a.type].append(a.get_drawable())
-        return result
-
     def add_annotation(self, klass, *args, **kw):
+        from cogent3.core.sequence import Sequence
+
+        if isinstance(self, Sequence):
+            from cogent3.util.warning import deprecated
+
+            deprecated(
+                "method",
+                "Sequence.add_annotation",
+                "_Annotatable.add_annotation",
+                " 2023.3",
+                "method .add_annotation will be discontinued for Sequence objects",
+            )
+
         annot = klass(self, *args, **kw)
         self.attach_annotations([annot])
         return annot
 
     def clear_annotations(self):
-        self.annotations = []
+        from cogent3.core.sequence import Sequence
 
-    def get_drawable(self, width=600, vertical=False):
-        """returns Drawable instance"""
-        from cogent3.draw.drawable import Drawable
+        if isinstance(self, Sequence):
+            from cogent3.util.warning import deprecated
 
-        drawables = self.get_drawables()
-        if not drawables:
-            return None
-        # we order by tracks
-        top = 0
-        space = 0.25
-        annotes = []
-        for feature_type in drawables:
-            new_bottom = top + space
-            for i, annott in enumerate(drawables[feature_type]):
-                annott.shift(y=new_bottom - annott.bottom)
-                if i > 0:
-                    annott._showlegend = False
-                annotes.append(annott)
+            deprecated(
+                "method",
+                "Sequence.clear_annotations",
+                "_Annotatable.clear_annotations",
+                " 2023.3",
+                "method .clear_annotations will be discontinued for Sequence objects",
+            )
 
-            top = annott.top
+        self.annotations = []
 
-        top += space
-        height = max((top / len(self)) * width, 300)
-        xaxis = dict(range=[0, len(self)], zeroline=False, showline=True)
-        yaxis = dict(range=[0, top], visible=False, zeroline=True, showline=True)
+    def attach_annotations(self, annots):
+        from cogent3.core.sequence import Sequence
 
-        if vertical:
-            all_traces = [t.T.as_trace() for t in annotes]
-            width, height = height, width
-            xaxis, yaxis = yaxis, xaxis
-        else:
-            all_traces = [t.as_trace() for t in annotes]
+        if isinstance(self, Sequence):
+            from cogent3.util.warning import deprecated
 
-        drawer = Drawable(
-            title=self.name, traces=all_traces, width=width, height=height
-        )
-        drawer.layout.update(xaxis=xaxis, yaxis=yaxis)
-        return drawer
+            deprecated(
+                "method",
+                "Sequence.attach_annotations",
+                "_Annotatable.attach_annotations",
+                " 2023.3",
+                "method .attach_annotations will be discontinued for Sequence objects",
+            )
 
-    def attach_annotations(self, annots):
         for annot in annots:
             if annot.parent is not self:
                 raise ValueError("doesn't belong here")
@@ -163,6 +253,18 @@ class _Annotatable:
             annot.attached = True
 
     def detach_annotations(self, annots):
+        from cogent3.core.sequence import Sequence
+
+        if isinstance(self, Sequence):
+            from cogent3.util.warning import deprecated
+
+            deprecated(
+                "method",
+                "Sequence.detach_annotations",
+                "_Annotatable.detach_annotations",
+                " 2023.3",
+                "method .detach_annotations will be discontinued for Sequence objects",
+            )
         for annot in annots:
             if annot.parent is not self:
                 raise ValueError("doesn't live here")
@@ -171,9 +273,6 @@ class _Annotatable:
                 self.annotations.remove(annot)
                 annot.attached = False
 
-    def add_feature(self, type, name, spans):
-        return self.add_annotation(Feature, type, name, spans)
-
     def get_annotations_matching(self, annotation_type, name=None, extend_query=False):
         """
 
@@ -189,6 +288,15 @@ class _Annotatable:
         -------
         list of AnnotatableFeatures
         """
+        from cogent3.util.warning import deprecated
+
+        deprecated(
+            "method",
+            "get_annotations_matching",
+            "get_features_matching",
+            " 2023.3",
+        )
+
         result = []
         if len(self.annotations) == 0:
             return result
@@ -208,6 +316,19 @@ class _Annotatable:
     def get_region_covering_all(
         self, annotations, feature_class=None, extend_query=False
     ):
+        from cogent3.core.sequence import Sequence
+
+        if isinstance(self, Sequence):
+            from cogent3.util.warning import deprecated
+
+            deprecated(
+                "method",
+                "Sequence.get_region_covering_all",
+                "_Annotatable.get_region_covering_all",
+                " 2023.3",
+                "method .get_region_covering_all will be discontinued for Sequence objects",
+            )
+
         if extend_query:
             annotations = [annot._projected_to_base(self) for annot in annotations]
         spans = []
@@ -236,7 +357,20 @@ class _Annotatable:
             current sequence are ignored.
 
         """
-        for annotation in self.get_annotations_matching(annotation_type, name):
+        from cogent3.core.sequence import Sequence
+
+        if isinstance(self, Sequence):
+
+            from cogent3.util.warning import deprecated
+
+            deprecated(
+                "method",
+                "Sequence.get_by_annotation",
+                "_Annotatable.get_by_annotation",
+                " 2023.3",
+                "method .get_by_annotation will be discontinued for Sequence objects",
+            )
+        for annotation in self.get_features_matching(annotation_type, name):
             try:
                 seq = self[annotation.map]
             except ValueError as msg:
@@ -289,6 +423,62 @@ class _Serialisable:
         return json.dumps(self.to_rich_dict())
 
 
+class FeatureNew(_AnnotationMixin, _Serialisable):
+    def __init__(self, parent, biotype, name, spans, reversed):
+        self.parent = parent
+        self.biotype = biotype
+        self.name = name
+        self.reversed = reversed
+        self.map = Map(locations=spans, parent_length=len(parent))
+
+    def get_slice(self):
+        # feature maps and spans are always oriented with respect to the + stand sequence, this means
+        # that when we really do need to be grabbing out the - strand sequence, that the features need to be
+        # be adjusted so that the features spans need to be made relative w.r.t. the end.
+
+        # todo: need to cope with cases where span is a lostSpan...
+
+        parent_reversed = self.parent._seq.reversed
+        feature_reversed = self.reversed
+
+        if not parent_reversed and not feature_reversed:
+            segments = [
+                str(self.parent[span.start : span.end]) for span in self.map.spans
+            ]
+
+        elif not parent_reversed and feature_reversed:
+            # because the underlying sequence is unchanged, just grab the segments out as strings
+            # join them, and then do the reverse complement
+            segments = [
+                str(self.parent[span.start : span.end]) for span in self.map.spans
+            ]
+            segments = self.parent.moltype.complement("".join(segments)[::-1])
+
+        elif parent_reversed and not feature_reversed:
+            parent = self.parent[::-1]
+            segments = [str(parent[span.start : span.end]) for span in self.map.spans]
+
+        elif parent_reversed and feature_reversed:
+            segments = [
+                str(self.parent[span.start : span.end])
+                for span in self.map.nucleic_reversed().reversed()
+            ]
+
+        return self.parent.moltype.make_seq(seq="".join(segments), name=self.name)
+
+    def get_children(self, biotype=None):
+        for child in self.parent.annotation_db.get_feature_children(
+            biotype=biotype, name=self.name
+        ):
+            yield FeatureNew(self.parent, **child)
+
+    def __len__(self):
+        return len(self.map)
+
+    def __repr__(self):
+        return f"{self.name} at {self.map}"
+
+
 # https://pythonspeed.com/products/filmemoryprofiler/
 
 
@@ -369,7 +559,7 @@ class _Feature(_Annotatable, _Serialisable):
 
     def get_shadow(self):
         return self.__class__(
-            self.parent, self.map.shadow(), type="region", name="not " + self.name
+            self.parent, self.map.shadow(), type="region", name=f"not {self.name}"
         )
 
     def __len__(self):
@@ -479,8 +669,8 @@ class _Variable(_Feature):
 def Variable(parent, type, name, xxy_list):
     """A variable that has 2 x-components (start, end) and a single y component.
     Currently used by Vestige - BMC Bioinformatics, 6:130, 2005."""
-    start = min([min(x1, x2) for ((x1, x2), y) in xxy_list])
-    end = max([max(x1, x2) for ((x1, x2), y) in xxy_list])
+    start = min(min(x1, x2) for ((x1, x2), y) in xxy_list)
+    end = max(max(x1, x2) for ((x1, x2), y) in xxy_list)
     if start != 0:
         xxy_list = [((x1 - start, x2 - start), y) for ((x1, x2), y) in xxy_list]
         end -= start
diff --git a/src/cogent3/core/annotation_db.py b/src/cogent3/core/annotation_db.py
new file mode 100644
index 000000000..d702627e9
--- /dev/null
+++ b/src/cogent3/core/annotation_db.py
@@ -0,0 +1,883 @@
+import collections
+import json
+import os
+import pathlib
+import re
+import sqlite3
+import typing
+
+import numpy
+
+from cogent3.util.deserialise import register_deserialiser
+from cogent3.util.misc import get_object_provenance
+from cogent3.util.table import Table
+
+
+__author__ = "Gavin Huttley"
+__copyright__ = "Copyright 2007-2022, The Cogent Project"
+__credits__ = ["Gavin Huttley"]
+__license__ = "BSD-3"
+__version__ = "2023.2.12a1"
+__maintainer__ = "Gavin Huttley"
+__email__ = "gavin.huttley@anu.edu.au"
+__status__ = "Production"
+
+
+OptionalInt = typing.Optional[int]
+OptionalStr = typing.Optional[str]
+OptionalBool = typing.Optional[bool]
+ReturnType = typing.Tuple[str, tuple]  # the sql statement and corresponding values
+
+
+# Define custom types for storage in sqlite
+# https://stackoverflow.com/questions/18621513/python-insert-numpy-array-into-sqlite3-database
+def array_to_sqlite(data):
+    return data.tobytes()
+
+
+def sqlite_to_array(data):
+    result = numpy.frombuffer(data, dtype=int)
+    dim = result.shape[0] // 2
+    return result.reshape((dim, 2))
+
+
+def dict_to_sqlite_as_json(data: dict) -> str:
+    return json.dumps(data)
+
+
+def sqlite_json_to_dict(data):
+    return json.loads(data)
+
+
+# registering the conversion functions with sqlite
+sqlite3.register_adapter(numpy.ndarray, array_to_sqlite)
+sqlite3.register_converter("array", sqlite_to_array)
+sqlite3.register_adapter(dict, dict_to_sqlite_as_json)
+sqlite3.register_converter("json", sqlite_json_to_dict)
+
+
+class FeatureDataType(typing.TypedDict):
+    seqid: str  # name of the sseq
+    biotype: str  # rename type attr of cogent3 Annotatables to match this?
+    name: str  # rename to name to match cogent3 Annotatable.name?
+    spans: list[tuple[int, int]]
+    reversed: bool  # True if feature on reverse strand
+
+
+@typing.runtime_checkable
+class SerialisableType(typing.Protocol):
+    def to_rich_dict(self) -> dict:
+        ...
+
+    def to_json(self) -> str:
+        ...
+
+    def from_dict(self, data):
+        ...
+
+
+@typing.runtime_checkable
+class SupportsQueryFeatures(typing.Protocol):  # should be defined centrally
+    def get_features_matching(
+        self,
+        biotype: OptionalStr = None,
+        name: OptionalStr = None,
+        start: OptionalInt = None,
+        end: OptionalInt = None,
+        strand: OptionalStr = None,
+        on_alignment: OptionalBool = None,
+    ) -> typing.Iterator[FeatureDataType]:
+        ...
+
+    def get_feature_children(
+        self, name: str, start: OptionalInt = None, end: OptionalInt = None
+    ) -> typing.List[FeatureDataType]:
+        ...
+
+    def get_feature_parent(
+        self, name: str, start: OptionalInt = None, end: OptionalInt = None
+    ) -> typing.List[FeatureDataType]:
+        ...
+
+
+@typing.runtime_checkable
+class SupportsWriteFeatures(typing.Protocol):  # should be defined centrally
+    def add_feature(
+        self,
+        seqid: str,
+        biotype: str,
+        name: str,
+        spans: typing.List[typing.Tuple[int, int]],
+        strand: str = None,
+        on_alignment: bool = False,
+    ) -> None:
+        ...
+
+    def add_records(
+        self,
+        records: typing.Sequence[dict],
+        # seqid required for genbank
+        seqid: typing.Optional[str] = None,
+    ) -> None:
+        ...
+
+
+@typing.runtime_checkable
+class SupportsFeatures(SupportsQueryFeatures, SupportsWriteFeatures, typing.Protocol):
+    ...
+
+
+def _make_table_sql(
+    table_name: str,
+    columns: dict,
+) -> str:
+    """makes the SQL for creating a table
+
+    Parameters
+    ----------
+    table_name : str
+        name of the table
+    columns : dict
+        {<column name>: <column SQL type>, ...}
+
+    Returns
+    -------
+    str
+    """
+    columns_types = ", ".join([f"{name} {ctype}" for name, ctype in columns.items()])
+    return f"CREATE TABLE IF NOT EXISTS {table_name} ({columns_types});"
+
+
+# for making sqlite regex queries see https://stackoverflow.com/questions/67209494/python3-sqlite-and-regex-query
+
+
+def _add_record_sql(
+    table_name: str,
+    data: dict,
+) -> ReturnType:
+    """creates SQL defining the table
+
+    Parameters
+    ----------
+    table_name : str
+        name of the table
+    data : dict
+        {<column name>: column type}
+
+    Returns
+    -------
+    str, tuple
+        the SQL statement and the tuple of values
+    """
+    cols = ", ".join(k.lower() for k in data)
+    pos = ", ".join("?" * len(data))
+    sql = f"INSERT INTO {table_name} ({cols}) VALUES ({pos});"
+    return sql, tuple(data.values())
+
+
+def _matching_conditions(
+    conditions: dict,
+    partial: bool = True,
+):
+    """creates WHERE clause
+
+    Parameters
+    ----------
+    conditions : dict
+        column name and values to be matched
+    partial : bool, optional
+        if False, only records within start, end are included. If True,
+        all records that overlap the segment defined by start, end are included.
+
+    Returns
+    -------
+    str, tuple
+        the SQL statement and the tuple of values
+    """
+    # todo this needs to support OR operation on some conditions, e.g. if the value of
+    # a condition is a tuple, do OR
+
+    start = conditions.pop("start", None)
+    end = conditions.pop("end", None)
+
+    sql = []
+    vals = ()
+    if conditions:
+        conds = []
+        vals = []
+        for col, val in conditions.items():
+            if val:
+                op = "LIKE" if "%" in val else "="
+                conds.append(f"{col} {op} ?")
+                vals.append(val)
+        sql.append(" AND ".join(conds))
+        vals = tuple(vals)
+
+    if isinstance(start, int) and isinstance(end, int):
+        # only matches within bounds
+        if partial:
+            cond = [
+                f"(start <= {start} AND end >= {start})",  # straddles beginning of segment
+                f"(start < {end} AND end >= {end})",  # straddles end of segment
+                f"(start <= {start} AND end >= {end})",  # includes segment
+            ]
+            cond = " OR ".join(cond)
+        else:
+            cond = f"start >= {start} AND end <= {end}"
+        sql.append(f"({cond})")
+    elif isinstance(start, int):
+        if partial:
+            cond = " OR ".join(
+                [f"start >= {start}", f"(start <= {start} AND end >= {start})"]
+            )
+        else:
+            cond = f"start >= {start}"
+        sql.append(f"({cond})")
+    elif isinstance(end, int):
+        if partial:
+            cond = " OR ".join([f"end <= {end}", f"(start <= {end} AND end >= {end})"])
+        else:
+            cond = f"end <= {end}"
+        sql.append(f"({cond})")
+
+    sql = f"{' AND '.join(sql)}"
+    return sql, vals
+
+
+def _del_records_sql(
+    table_name: str,
+    conditions: dict,
+    start: OptionalInt = None,
+    end: OptionalInt = None,
+    partial=True,
+) -> ReturnType:
+    """creates the SQL and values for identifying records to be deleted
+
+    Parameters
+    ----------
+    table_name : str
+        table to have records deleted from
+    conditions : dict
+        column name and values to be matched
+    start, end : OptionalInt
+        select records whose (start, end) values lie between start and end,
+        or overlap them if (partial is True)
+    partial : bool, optional
+        if False, only records within start, end are included. If True,
+        all records that overlap the segment defined by start, end are included.
+
+    Returns
+    -------
+    str, tuple
+        the SQL statement and the tuple of values
+    """
+    where, vals = _matching_conditions(
+        conditions=conditions, start=start, end=end, partial=partial
+    )
+    sql = f"DELETE FROM {table_name}"
+    if not where:
+        return sql
+
+    sql = f"{sql} WHERE {where};"
+    return sql, vals
+
+
+def _select_records_sql(
+    table_name: str,
+    conditions: dict,
+    columns: typing.Optional[typing.List[str]] = None,
+    start: OptionalInt = None,
+    end: OptionalInt = None,
+    partial=True,
+) -> ReturnType:
+    """create SQL select statement and values
+
+    Parameters
+    ----------
+    table_name : str
+        containing the data to be selected from
+    columns : Tuple[str]
+        values to select
+    conditions : dict
+        the WHERE conditions
+    start, end : OptionalInt
+        select records whose (start, end) values lie between start and end,
+        or overlap them if (partial is True)
+    partial : bool, optional
+        if False, only records within start, end are included. If True,
+        all records that overlap the segment defined by start, end are included.
+
+    Returns
+    -------
+    str, tuple
+        the SQL statement and the tuple of values
+    """
+
+    where, vals = _matching_conditions(conditions=conditions, partial=partial)
+    columns = f"{', '.join(columns)}" if columns else "*"
+    sql = f"SELECT {columns} FROM {table_name}"
+    if not where:
+        return sql
+
+    sql = f"{sql} WHERE {where};"
+    return sql, vals
+
+
+# todo add support for querying for text within the additional feature,
+# for example, add a "attrs_like" argument which users can provide text
+# that will be treated as surrounded by wild-cards
+class SqliteAnnotationDbMixin:
+    # table schema for user provided annotations
+    _user_schema = {
+        "biotype": "TEXT",
+        "seqid": "TEXT",
+        "name": "TEXT",
+        "parent_id": "TEXT",
+        "start": "INTEGER",
+        "end": "INTEGER",
+        "strand": "TEXT",
+        "spans": "array",
+        "on_alignment": "INT",
+    }
+
+    @property
+    def table_names(self) -> tuple[str]:
+        return self._table_names
+
+    def _init_tables(self):
+        # bit of magic
+        # assumes schema attributes named as `_<table name>_schema`
+        for table_name in self.table_names:
+            attr = getattr(self, f"_{table_name}_schema")
+            sql = _make_table_sql(table_name, attr)
+            self._execute_sql(sql)
+
+    @property
+    def db(self) -> sqlite3.Connection:
+        if self._db is None:
+            self._db = sqlite3.connect(
+                self.source, detect_types=sqlite3.PARSE_DECLTYPES
+            )
+            self._db.row_factory = sqlite3.Row
+
+        return self._db
+
+    def _execute_sql(self, cmnd, values=None):
+        with self.db:
+            # context manager ensures safe transactions
+            cursor = self.db.cursor()
+            cursor.execute(cmnd, values or [])
+            return cursor
+
+    def add_feature(
+        self,
+        seqid: str,
+        biotype: str,
+        name: str,
+        spans: typing.List[typing.Tuple[int, int]],
+        strand: str = None,
+        on_alignment: bool = False,
+    ) -> None:
+        """adds a record to user table
+
+        Parameters
+        ----------
+        seqid : str
+            name of the sequence feature resides on
+        biotype : str
+            biological type of the record
+        name : str
+            the name of a record, an identifier
+        spans : typing.List[typing.Tuple[int, int]], optional
+            this will be sorted
+        strand : str, optional
+            either +, -. Defaults to '+'
+        on_alignment : bool, optional
+            whether the annotation is an alignment annotation
+        """
+        spans = numpy.array(sorted(sorted(coords) for coords in spans), dtype=int)
+        # the start, end variables are added into the record via the loop
+        # on local variables
+        start = int(spans.min())
+        end = int(spans.max())
+        # we define record as all defined variables from local name space,
+        # excluding "self"
+        record = {k: v for k, v in locals().items() if k != "self" and v is not None}
+        sql, values = _add_record_sql("user", record)
+        self._execute_sql(sql, values=values)
+
+    def _get_records_matching(
+        self, table_name: str, **kwargs
+    ) -> typing.Iterator[FeatureDataType]:
+        """return all fields"""
+        columns = kwargs.pop("columns", None)
+        sql, vals = _select_records_sql(
+            table_name=table_name, conditions=kwargs, columns=columns
+        )
+        yield from self._execute_sql(sql, values=vals)
+
+    def get_records_matching(
+        self,
+        biotype: str = None,
+        seqid: str = None,
+        name: str = None,
+        start: int = None,
+        end: int = None,
+        strand: bool = None,
+        on_alignment: bool = None,
+    ) -> typing.Iterator[dict]:
+        """return all fields for matching records"""
+        # a record is Everything, a Feature is a subset
+        # we define query as all defined variables from local name space,
+        # excluding "self" and kwargs at default values
+        kwargs = {k: v for k, v in locals().items() if k != "self" and v is not None}
+        for table_name in self.table_names:
+            for result in self._get_records_matching(table_name, **kwargs):
+                yield {k: result[k] for k in result.keys()}
+
+    def get_features_matching(
+        self,
+        biotype: str = None,
+        seqid: str = None,
+        name: str = None,
+        start: int = None,
+        end: int = None,
+        strand: bool = None,
+        on_alignment: bool = None,
+    ) -> typing.Iterator[FeatureDataType]:
+        # returns essential values to create a Feature
+        # we define query as all defined variables from local name space,
+        # excluding "self" and kwargs at default values
+        kwargs = {k: v for k, v in locals().items() if k != "self" and v is not None}
+        columns = ("seqid", "biotype", "spans", "strand", "name")
+        for table_name in self.table_names:
+            for result in self._get_records_matching(
+                table_name=table_name, columns=columns, **kwargs
+            ):
+                yield FeatureDataType(
+                    seqid=result["seqid"],
+                    biotype=result["biotype"],
+                    name=result["name"],
+                    spans=[tuple(c) for c in result["spans"]],
+                    reversed=result["strand"] == "-",
+                )
+
+    @property
+    def describe(self) -> Table:
+        """top level description of the annotation db"""
+        sql_template = "SELECT DISTINCT {} FROM {};"
+        data = {}
+        for column in ("seqid", "biotype"):
+            data[column] = set()
+            for table in self.table_names:
+                sql = sql_template.format(column, table)
+                if result := self._execute_sql(sql).fetchall():
+                    data[column] |= {r[column] for r in result}
+
+        row_counts = []
+        for table in self.table_names:
+            result = self._execute_sql(f"SELECT COUNT(*) FROM {table}").fetchone()
+            row_counts.append(result["COUNT(*)"])
+
+        from cogent3 import make_table
+
+        table = make_table(
+            data={
+                "biotype": list(data.keys())
+                + [f"num_rows({t!r})" for t in self.table_names],
+                "count": [len(v) for v in data.values()] + row_counts,
+            }
+        )
+        return table
+
+    def biotype_counts(self) -> dict:
+        """return counts of biological types across all tables and seqids"""
+        sql_template = "SELECT biotype FROM {}"
+        counts = collections.Counter()
+        for table in self.table_names:
+            sql = sql_template.format(table)
+            if result := self._execute_sql(sql).fetchall():
+                counts.update(v for r in result if (v := r["biotype"]))
+        return counts
+
+    def to_rich_dict(self):
+        # todo drop columns from records with no value
+        # top level dict for each table_name
+        records = {}
+        for table in self.table_names:
+            records[table] = []
+            for record in self._execute_sql(f"SELECT * from {table};"):
+                record = {k: record[k] for k in record.keys()}
+                record["spans"] = record[
+                    "spans"
+                ].tolist()  # required for json serialisation
+                records[table].append(
+                    {k: v for k, v in record.items() if v is not None}
+                )
+
+        records["type"] = get_object_provenance(self)
+        return records
+
+    @classmethod
+    def from_dict(cls, data: dict):
+        data.pop("type")
+        # make an empty db
+        db = cls([])
+        for table_name, records in data.items():
+            for record in records:
+                record["spans"] = numpy.array(record["spans"], dtype=int)
+                sql, vals = _add_record_sql(
+                    table_name, {k: v for k, v in record.items() if v is not None}
+                )
+                db._execute_sql(sql, vals)
+        return db
+
+
+class GffAnnotationDb(SqliteAnnotationDbMixin):
+    _table_names = "gff", "user"
+    # am relying on name structured as _<table name>_schema
+    _gff_schema = {
+        "seqid": "TEXT",
+        "source": "TEXT",
+        "biotype": "TEXT",  # type in GFF
+        "start": "INTEGER",
+        "end": "INTEGER",
+        "score": "TEXT",  # check defn
+        "strand": "TEXT",
+        "phase": "TEXT",
+        "attributes": "TEXT",
+        "comments": "TEXT",
+        "spans": "array",  # aggregation of coords across records
+        "name": "TEXT",
+        "parent_id": "TEXT",
+    }
+
+    def __init__(self, data, db=None, source=":memory:"):
+        # note that data is destroyed
+        self._num_fakeids = 0
+        self.source = source
+        self._db = db
+        if db is None:
+            self._init_tables()
+        data = self._merged_data(data)
+        self.add_records(data)
+
+    def add_records(self, reduced):
+        # Can we really trust text case of "ID" and "Parent" to be consistent
+        # across sources of gff?? I doubt it, so more robust regex likely
+        # required
+        for record in reduced.values():
+            # our Feature code assumes start always < end,
+            # we record direction using Strand
+            spans = numpy.array(sorted(record["spans"]), dtype=int)  # sorts the rows
+            spans.sort(axis=1)
+            record["Start"] = int(spans.min())
+            record["End"] = int(spans.max())
+            record["spans"] = spans
+            cmnd, vals = _add_record_sql(
+                self.table_names[0],
+                {
+                    k: v
+                    for k, v in record.items()
+                    if isinstance(v, numpy.ndarray) or v not in (".", None)
+                },
+            )
+            self._execute_sql(cmnd=cmnd, values=vals)
+
+        del reduced
+
+    def _merged_data(self, records) -> dict:
+        field_template = r"(?<={}=)[^;\s]+"
+        name = re.compile(field_template.format("ID"))
+        parent_id = re.compile(field_template.format("Parent"))
+
+        reduced = collections.OrderedDict()
+        # collapse records with ID's occurring in multiple rows into one
+        # row, converting their coordinates
+        # extract the name from ID and add this into the table
+        # I am not convinced we can rely on gff files to be ordered,
+        # if we could, we could do this as one pass over the data
+        while records:
+            record = records.pop(0)
+            record["biotype"] = record.pop("Type")
+            attrs = record["Attributes"] or ""
+            if match := name.search(attrs):
+                record_id = match.group()
+            else:
+                record_id = f"unknown-{self._num_fakeids}"
+                self._num_fakeids += 1
+
+            record["name"] = record_id
+            if pid := parent_id.search(attrs):
+                record["parent_id"] = pid.group()
+
+            if record_id not in reduced:
+                reduced[record_id] = record
+                reduced[record_id]["spans"] = []
+
+            # should this just be an append?
+            reduced[record_id]["spans"].append((record["Start"], record["End"]))
+        return reduced
+
+    def _get_feature_by_id(
+        self, table_name: str, column: str, name: str, biotype: OptionalStr = None
+    ) -> typing.List[FeatureDataType]:
+        # we return the parent_id because `get_feature_parent()` requires it
+        sql, vals = _select_records_sql(
+            table_name=table_name,
+            conditions={column: name, "biotype": biotype},
+            columns=["biotype", "spans", "strand", "name", "parent_id"],
+        )
+        for result in self._execute_sql(sql, values=vals):
+            yield {
+                "biotype": result["biotype"],
+                "name": result["name"],
+                "spans": [tuple(c) for c in result["spans"]],
+                "reversed": result["strand"] == "-",
+                "parent_id": result["parent_id"],
+            }
+
+    def get_feature_children(
+        self, name: str, biotype: OptionalStr = None, **kwargs
+    ) -> typing.Iterator[FeatureDataType]:
+        """yields children of name"""
+        # kwargs is used because other classes need start / end
+        # just uses search matching
+        for table_name in self.table_names:
+            for result in self._get_feature_by_id(
+                table_name=table_name,
+                column="parent_id",
+                name=f"%{name}%",
+                biotype=biotype,
+            ):
+                result.pop("parent_id")  # remove invalid field for the FeatureDataType
+                yield result
+
+    def get_feature_parent(
+        self, name: str, **kwargs
+    ) -> typing.Iterator[FeatureDataType]:
+        """yields parents of name"""
+        for table_name in self.table_names:
+            for result in self._get_feature_by_id(
+                table_name=table_name, column="name", name=f"%{name}%"
+            ):
+                # multiple values for parent means this is better expressed
+                # as an OR clause
+                # todo modify the conditional SQL generation
+                for name in result["parent_id"].replace(" ", "").split(","):
+                    if parent := list(
+                        self._get_feature_by_id(
+                            table_name=table_name,
+                            column="name",
+                            name=name,
+                        )
+                    ):
+                        parent = parent[0]
+                        parent.pop(
+                            "parent_id"
+                        )  # remove invalid field for the FeatureDataType
+                        yield parent
+
+
+# The GenBank format is less clear on the relationship between identifiers,
+# e.g. it can be gene -> SRP_RNA -> exon, gene -> CDS, etc... Without having
+# the explicit type-hierarchy from NCBI (which I have not yet found) the only way
+# to establish the Parent of a feature is by knowing what other
+# features have the same ID and then which of those are candidates to contain
+# a feature based on the hierarchy of relationships.
+# In other words, we assume that children have the same ID as their parent BUT but
+# their span lies within the parents. Conversely, for a parent query, the parent
+# has the same ID as their child but their span contains the childs.
+
+
+class GenbankAnnotationDb(SqliteAnnotationDbMixin):
+    _table_names = "gb", "user"
+    # am relying on name structured as _<table name>_schema
+    _gb_schema = {
+        "seqid": "TEXT",
+        "source": "TEXT",
+        "biotype": "TEXT",  # type in GFF
+        "start": "INTEGER",
+        "end": "INTEGER",
+        "strand": "TEXT",
+        "comments": "TEXT",
+        "spans": "array",  # aggregation of coords across records
+        "name": "TEXT",
+        "symbol": "TEXT",
+        "parent_id": "TEXT",
+        "attributes": "json",
+    }
+
+    def __init__(self, data: typing.List[dict], seqid: str, source=":memory:", db=None):
+        # note that data is destroyed
+        self._db = db
+        self._num_fakeids = 0
+        self.source = source
+        self._init_tables()
+        self.add_records(data, seqid)
+
+    def add_records(self, records, seqid):
+        # Can we really trust case to be consistent across sources of gff??
+        # I doubt it, so more robust regex likely required
+        # Can we even rely on these concepts being represented by the text
+        # ID and Parent?
+
+        # need to capture genetic code from genbank, but what a
+
+        col_key_map = {"type": "biotype", "locus_tag": "name"}
+        exclude = {"translation", "location"}  # location is grabbed directly
+        for record in records:
+            # our Feature code assumes start always < end,
+            store = {"seqid": seqid}
+            # we create the location data directly
+            if location := record.get("location", None):
+                store["spans"] = numpy.array(
+                    sorted([sorted((s.first(), s.last())) for s in location]), dtype=int
+                )
+                if strand := location.strand:
+                    store["strand"] = strand
+                store["start"] = int(store["spans"].min())
+                store["end"] = int(store["spans"].max())
+
+            record_keys = col_key_map.keys() & record.keys()
+            attrs_keys = record.keys() - col_key_map.keys() - exclude
+            store |= {col_key_map[k]: record[k] for k in record_keys}
+            store["attributes"] = {k: record[k] for k in attrs_keys}
+            if "name" not in store:
+                store["name"] = [f"fakeid-{self._num_fakeids}"]
+                self._num_fakeids += 1
+
+            store["name"] = ",".join(store["name"])
+            cmnd, vals = _add_record_sql(
+                self.table_names[0],
+                store,
+            )
+            self._execute_sql(cmnd=cmnd, values=vals)
+
+    def _get_feature_by_id(
+        self,
+        table_name: str,
+        name: str,
+        start: int,
+        end: int,
+        biotype: OptionalStr = None,
+        partial: bool = False,
+    ) -> typing.List[FeatureDataType]:
+        # we return the parent_id because `get_feature_parent()` requires it
+        sql, vals = _select_records_sql(
+            table_name=table_name,
+            conditions={"name": name, "biotype": biotype},
+            start=start,
+            end=end,
+            columns=["biotype", "start", "end", "spans", "strand", "name", "parent_id"],
+            partial=partial,
+        )
+        for result in self._execute_sql(sql, values=vals):
+            yield {
+                "biotype": result["biotype"],
+                "name": name,
+                "start": result["start"],
+                "end": result["end"],
+                "spans": [tuple(c) for c in result["spans"]],
+                "reverse": result["strand"] == "-",
+            }
+
+    def get_feature_children(
+        self,
+        name: str,
+        biotype: OptionalStr = None,
+        exclude_biotype: OptionalStr = None,
+        start: OptionalInt = None,
+        end: OptionalInt = None,
+    ) -> typing.Iterator[FeatureDataType]:
+        """yields children of name"""
+        # children must lie within the parent coordinates
+        for table in self.table_names:
+            for feat in self._get_feature_by_id(
+                table_name=table,
+                name=name,
+                biotype=biotype,
+                start=int(start),
+                end=int(end),
+                partial=False,
+            ):
+                if feat["biotype"] == exclude_biotype:
+                    continue
+                cstart, cend = feat.pop("start"), feat.pop("end")
+                if not (start <= cstart < end and start < cend <= end):
+                    continue
+                yield feat
+
+    def get_feature_parent(
+        self,
+        name: str,
+        exclude_biotype: OptionalStr = None,
+        start: OptionalInt = None,
+        end: OptionalInt = None,
+    ) -> typing.Iterator[FeatureDataType]:
+        """yields parents of name"""
+        # parent must match or lie outside the coordinates
+        for table in self.table_names:
+            for feat in self._get_feature_by_id(
+                table_name=table,
+                name=name,
+                start=int(start),
+                end=int(end),
+                partial=False,
+            ):
+                # add support for != operation to SQL where clause generation
+                if feat["biotype"] == exclude_biotype:
+                    continue
+
+                cstart, cend = feat.pop("start"), feat.pop("end")
+                if cstart > start or end > cend:
+                    continue
+                yield feat
+
+
+@register_deserialiser(get_object_provenance(GffAnnotationDb))
+def deserialise_gff_db(data: dict):
+    return GffAnnotationDb.from_dict(data)
+
+
+@register_deserialiser(get_object_provenance(GenbankAnnotationDb))
+def deserialise_gb_db(data: dict):
+    return GenbankAnnotationDb.from_dict(data)
+
+
+def _db_from_genbank(path, db):
+    from cogent3 import open_
+    from cogent3.parse.genbank import MinimalGenbankParser
+
+    with open_(path) as infile:
+        data = list(MinimalGenbankParser(infile))
+
+    return GenbankAnnotationDb(data[0]["features"], data[0]["locus"], db=db)
+
+
+def _leave_attributes(*attrs):
+    return attrs[0]
+
+
+def _db_from_gff(path, seqids, db):
+    from cogent3.parse.gff import gff_parser
+
+    data = list(
+        gff_parser(
+            path,
+            seqids=seqids,
+            attribute_parser=_leave_attributes,
+        )
+    )
+    return GffAnnotationDb(data=data, db=db)
+
+
+def load_annotations(path: os.PathLike, seqids=None, db=None) -> SupportsFeatures:
+    if seqids is not None:
+        seqids = {seqids} if isinstance(seqids, str) else set(seqids)
+    path = pathlib.Path(path)
+    return (
+        _db_from_genbank(path, db=db)
+        if ".gb" in path.suffixes
+        else _db_from_gff(path, seqids=seqids, db=db)
+    )
diff --git a/src/cogent3/core/genetic_code.py b/src/cogent3/core/genetic_code.py
index dab09dd2e..a63ae0ca6 100644
--- a/src/cogent3/core/genetic_code.py
+++ b/src/cogent3/core/genetic_code.py
@@ -346,6 +346,13 @@ class GeneticCode:
 
         return "".join(mappings)
 
+    def get_alphabet(self, include_stop=False):
+        from cogent3.core.alphabet import Alphabet
+        from cogent3.core.moltype import DNA
+
+        motifs = list(self.codons) if include_stop else list(self.sense_codons)
+        return Alphabet(motifs, moltype=DNA)
+
 
 NcbiGeneticCodeData = [
     GeneticCode(*data)
diff --git a/src/cogent3/core/location.py b/src/cogent3/core/location.py
index ceca7e6ec..9ba70b2cb 100644
--- a/src/cogent3/core/location.py
+++ b/src/cogent3/core/location.py
@@ -575,11 +575,6 @@ class Map(object):
                     raise RuntimeError(
                         f"located outside sequence: {str((start, end, parent_length))}"
                     )
-                elif min(start, end) > parent_length:
-                    diff = max(start, end) - parent_length
-                    start = [start, parent_length][start > parent_length]
-                    end = [end, parent_length][end > parent_length]
-
                 span = Span(start, end, tidy, tidy, reverse=reverse)
                 if diff < 0:
                     spans += [LostSpan(-diff), span]
diff --git a/src/cogent3/core/moltype.py b/src/cogent3/core/moltype.py
index 99328fc46..ae93da4b0 100644
--- a/src/cogent3/core/moltype.py
+++ b/src/cogent3/core/moltype.py
@@ -29,8 +29,6 @@ from copy import deepcopy
 from random import choice
 from string import ascii_letters as letters
 
-import numpy
-
 from cogent3.core.alignment import (
     Alignment,
     ArrayAlignment,
@@ -66,7 +64,6 @@ from cogent3.util.misc import (
     FunctionWrapper,
     add_lowercase,
     get_object_provenance,
-    iterable,
 )
 from cogent3.util.transform import KeepChars, first_index_in_set
 from cogent3.util.warning import deprecated
@@ -528,6 +525,10 @@ AA_COLORS = _expand_colors(
 )
 
 
+def _do_nothing(x):
+    return x
+
+
 class MolType(object):
     """MolType: Handles operations that depend on the sequence type (e.g. DNA).
 
@@ -559,6 +560,7 @@ class MolType(object):
         make_alphabet_group=False,
         array_seq_constructor=None,
         colors=None,
+        coerce_string=None,
     ):
         """Returns a new MolType object. Note that the parameters are in flux.
 
@@ -701,6 +703,8 @@ class MolType(object):
 
         self._colors = colors or defaultdict(_DefaultValue("black"))
 
+        self._coerce_string = coerce_string or _do_nothing
+
     def __repr__(self):
         """String representation of MolType.
 
@@ -713,6 +717,9 @@ class MolType(object):
         data = self.to_rich_dict(for_pickle=True)
         return (), data
 
+    def coerce_str(self, data: str) -> str:
+        return self._coerce_string(data)
+
     def to_rich_dict(self, for_pickle=False):
         data = deepcopy(self._serialisable)
         if not for_pickle:  # we rely on reconstruction from label
@@ -748,7 +755,11 @@ class MolType(object):
 
     def make_seq(self, seq, name=None, **kwargs):
         """Returns sequence of correct type."""
-        return self._make_seq(seq, name, **kwargs)
+
+        name = name or getattr(seq, "name", None)
+        if isinstance(seq, ArraySequence):
+            seq = str(seq)
+        return self._make_seq(self.coerce_str(seq), name, **kwargs)
 
     def make_array_seq(self, seq, name=None, **kwargs):
         """
@@ -766,6 +777,7 @@ class MolType(object):
         -------
         ArraySequence
         """
+        name = name or getattr(seq, "name", None)
         alphabet = kwargs.pop("alphabet", None)
         if alphabet is None and hasattr(self, "alphabets"):
             alphabet = self.alphabets.degen_gapped
@@ -1295,6 +1307,14 @@ class MolType(object):
         return css, styles
 
 
+def _convert_to_rna(seq: str) -> str:
+    return seq.replace("t", "u").replace("T", "U")
+
+
+def _convert_to_dna(seq: str) -> str:
+    return seq.replace("u", "t").replace("U", "T")
+
+
 ASCII = MolType(
     # A default type for text read from a file etc. when we don't
     # want to prematurely assume DNA or Protein.
@@ -1316,6 +1336,7 @@ DNA = MolType(
     make_alphabet_group=True,
     array_seq_constructor=ArrayDnaSequence,
     colors=NT_COLORS,
+    coerce_string=_convert_to_dna,
 )
 
 RNA = MolType(
@@ -1329,6 +1350,7 @@ RNA = MolType(
     make_alphabet_group=True,
     array_seq_constructor=ArrayRnaSequence,
     colors=NT_COLORS,
+    coerce_string=_convert_to_rna,
 )
 
 PROTEIN = MolType(
@@ -1363,6 +1385,7 @@ BYTES = MolType(
     label="bytes",
 )
 
+
 # the None value catches cases where a moltype has no label attribute
 _style_defaults = {
     getattr(mt, "label", ""): defaultdict(
@@ -1380,6 +1403,10 @@ AB = MolType(
     label="ab",
 )
 
+# todo the _CodonAlphabet class should not exist,
+# the genetic code should have an alphabet, not
+# the alphabet has a genetic code
+
 
 class _CodonAlphabet(Alphabet):
     """Codon alphabets are DNA TupleAlphabets with a genetic code attribute and some codon-specific methods"""
@@ -1413,15 +1440,9 @@ def CodonAlphabet(gc=1, include_stop_codons=False):
     return a
 
 
-def _method_codon_alphabet(ignore, *args, **kwargs):
-    """If CodonAlphabet is set as a property, it gets self as extra 1st arg."""
-    return CodonAlphabet(*args, **kwargs)
-
-
 STANDARD_CODON = CodonAlphabet()
 
 # Modify NucleicAcidSequence to avoid circular import
-NucleicAcidSequence.codon_alphabet = _method_codon_alphabet
 NucleicAcidSequence.protein = PROTEIN
 ArrayRnaSequence.moltype = RNA
 ArrayRnaSequence.alphabet = RNA.alphabets.degen_gapped
diff --git a/src/cogent3/core/sequence.py b/src/cogent3/core/sequence.py
index d92c55784..60e89256c 100644
--- a/src/cogent3/core/sequence.py
+++ b/src/cogent3/core/sequence.py
@@ -10,12 +10,14 @@ Sequences are intended to be immutable. This is not enforced by the code for
 performance reasons, but don't alter the MolType or the sequence data after
 creation.
 """
+
+import contextlib
 import json
+import os
 import re
 import warnings
 
-from collections import defaultdict
-from functools import total_ordering
+from functools import singledispatch, total_ordering
 from operator import eq, ne
 from random import shuffle
 from typing import Generator, List, Tuple
@@ -35,13 +37,17 @@ from numpy import (
 from numpy.random import permutation
 
 from cogent3.core.alphabet import AlphabetError
-from cogent3.core.annotation import Map, _Annotatable
+from cogent3.core.annotation import FeatureNew, Map, _Annotatable
+from cogent3.core.annotation_db import (
+    GenbankAnnotationDb,
+    GffAnnotationDb,
+    load_annotations,
+)
 from cogent3.core.genetic_code import get_code
 from cogent3.core.info import Info as InfoClass
 from cogent3.format.fasta import alignment_to_fasta
 from cogent3.maths.stats.contingency import CategoryCounts
 from cogent3.maths.stats.number import CategoryCounter
-from cogent3.parse import gff
 from cogent3.util.dict_array import DictArrayTemplate
 from cogent3.util.misc import (
     DistanceFromMatrix,
@@ -774,8 +780,7 @@ class SequenceI(object):
         else:
             name = None
 
-        new_seq = self.__class__(str(self) + other_seq, name=name)
-        return new_seq
+        return self.__class__(str(self) + other_seq, name=name)
 
 
 @total_ordering
@@ -804,37 +809,19 @@ class Sequence(_Annotatable, SequenceI):
 
             check: if True (the default), validates against the MolType
         """
+
         if name is None and hasattr(seq, "name"):
             name = seq.name
         self.name = name
         orig_seq = seq
-        if isinstance(seq, Sequence):
-            seq = str(seq)
-        elif isinstance(seq, ArraySequence):
-            seq = str(seq)
-        elif isinstance(seq, bytes):
-            seq = seq.decode("utf-8")
-        elif not isinstance(seq, str):
-            try:
-                seq = "".join(seq)
-            except TypeError:
-                seq = "".join(map(str, seq))
 
-        seq = self._seq_filter(seq)
-
-        if isinstance(seq, SeqView):
-            if not preserve_case and not str(seq).isupper():
-                seq.seq = seq.seq.upper()
-            self._seq = seq
-
-        elif not preserve_case and not seq.isupper():
-            seq = seq.upper()
-
-        if isinstance(seq, str):
-            self._seq = SeqView(seq)
+        checker = (
+            (lambda x: self.moltype.verify_sequence(x, gaps_allowed, wildcards_allowed))
+            if check
+            else (lambda x: x)
+        )
 
-        if check:
-            self.moltype.verify_sequence(str(self), gaps_allowed, wildcards_allowed)
+        self._seq = _coerce_seq(seq, preserve_case, checker)
 
         if not isinstance(info, InfoClass):
             try:
@@ -842,10 +829,8 @@ class Sequence(_Annotatable, SequenceI):
             except TypeError:
                 info = InfoClass()
         if hasattr(orig_seq, "info"):
-            try:
+            with contextlib.suppress(Exception):
                 info.update(orig_seq.info)
-            except:
-                pass
         self.info = info
 
         if isinstance(orig_seq, _Annotatable):
@@ -854,6 +839,109 @@ class Sequence(_Annotatable, SequenceI):
 
         self._repr_policy = dict(num_pos=60)
 
+        self._annotation_db = None
+
+    @property
+    def annotation_offset(self):
+        """
+        The offset between annotation coordinates and sequence coordinates.
+
+        The offset can be used to adjust annotation coordinates to match the position
+        of the given Sequence within a larger genomic context. For example, if the
+        Annotations are with respect to a chromosome, and the sequence represents
+        a gene that is 100 bp from the start of a chromosome, the offset can be set to
+        100 to ensure that the gene's annotations are aligned with the appropriate
+        genomic positions.
+
+
+        Returns:
+            int: The offset between annotation coordinates and sequence coordinates.
+        """
+
+        return self._seq.offset
+
+    @annotation_offset.setter
+    def annotation_offset(self, value):
+        self._seq.offset = value
+
+    @property
+    def annotation_db(self):
+        return self._annotation_db
+
+    @annotation_db.setter
+    def annotation_db(self, value):
+        from cogent3.core.annotation_db import SupportsFeatures
+
+        if not isinstance(value, SupportsFeatures):
+            raise TypeError
+        self._annotation_db = value
+
+    def get_features_matching(
+        self,
+        feature_type=None,
+        name=None,
+        start=None,
+        stop=None,
+    ):
+        """yield features matching the given conditions"""
+
+        if self._annotation_db is None:
+            return None
+
+        query_start = self._seq.absolute_index(start) if start is not None else None
+        query_end = self._seq.absolute_index(stop) if stop is not None else None
+
+        for feature in self.annotation_db.get_features_matching(
+            seqid=self.name,
+            name=name,
+            biotype=feature_type,
+            start=query_start,
+            end=query_end,
+        ):
+            feature.pop("seqid")  # not required here as provided
+            yield FeatureNew(self, **feature)
+
+    def annotate_from_gff(self, f: os.PathLike, offset=None):
+        """copies annotations from a gff file to self,
+
+        Parameters
+        ----------
+        f : path to gff annotation file.
+        offset : Optional, the offset between annotation coordinates and sequence coordinates.
+
+        """
+        if self.annotation_db is None:
+            self.annotation_db = load_annotations(f, self.name)
+        elif isinstance(self.annotation_db, GenbankAnnotationDb):
+            raise ValueError("GenbankAnnotationDb already attached")
+        else:
+            self.annotation_db = load_annotations(
+                f, self.name, db=self.annotation_db._db
+            )
+
+        if offset:
+            self.annotation_offset = offset
+
+    def add_feature(
+        self,
+        biotype: str,
+        name: str,
+        spans,
+        strand: str = None,
+        on_alignment: bool = False,
+    ):
+        if self.annotation_db is None:
+            self.annotation_db = GffAnnotationDb([])
+
+        self.annotation_db.add_feature(
+            seqid=self.name,
+            biotype=biotype,
+            name=name,
+            spans=spans,
+            strand=strand,
+            on_alignment=on_alignment,
+        )
+
     def to_moltype(self, moltype):
         """returns copy of self with moltype seq
 
@@ -868,8 +956,11 @@ class Sequence(_Annotatable, SequenceI):
             raise ValueError(f"unknown moltype '{moltype}'")
 
         moltype = get_moltype(moltype)
-        make_seq = moltype.make_seq
-        new = make_seq(self, name=self.name)
+        s = moltype.coerce_str(self._seq.value)
+        moltype.verify_sequence(s, gaps_allowed=True, wildcards_allowed=True)
+        sv = SeqView(s)
+        new = moltype.make_seq(sv, name=self.name, info=self.info)
+
         new.clear_annotations()
         for ann in self.annotations:
             ann.copy_annotations_to(new)
@@ -901,58 +992,6 @@ class Sequence(_Annotatable, SequenceI):
 
         return offset + start
 
-    def annotate_from_gff(self, f, pre_parsed=False):
-        """annotates a Sequence from a gff file where each entry has the same SeqID"""
-        # only features with parent features included in the 'features' dict
-        gff_contents = f if pre_parsed else gff.gff_parser(f)
-        top_level = defaultdict(list)
-        grouped = defaultdict(list)
-        num_no_id = 0
-        for gff_dict in gff_contents:
-            if gff_dict["SeqID"] != self.name:
-                # we can only handle features for this sequence
-                continue
-
-            id_ = gff_dict["Attributes"]["ID"]
-            parents = gff_dict["Attributes"].get("Parent", None)
-            if parents is None and id_:
-                assert id_ not in top_level, f"non-unique id {id_}"
-                top_level[id_].append(
-                    self.add_feature(
-                        gff_dict["Type"], id_, [(gff_dict["Start"], gff_dict["End"])]
-                    )
-                )
-            elif parents is None:
-                id_ = f"no-id-{num_no_id}"
-                num_no_id += 1
-                self.add_feature(
-                    gff_dict["Type"], id_, [(gff_dict["Start"], gff_dict["End"])]
-                )
-            else:
-                for parent in parents:
-                    grouped[parent].append(gff_dict)
-
-        # we annotate the annotations
-        while grouped:
-            for key, features in top_level.items():
-                child_features = grouped.pop(key, [])
-                if child_features:
-                    break
-
-            for feature in features:
-                feature_start = self._get_feature_start(feature)
-                for gff_dict in child_features:
-                    id_ = gff_dict["Attributes"]["ID"]
-                    b = gff_dict["Start"]
-                    e = gff_dict["End"]
-                    type_ = gff_dict["Type"]
-                    sub_feat = feature.add_feature(
-                        type_,
-                        id_,
-                        [(b - feature_start, e - feature_start)],
-                    )
-                    top_level[gff_dict["Attributes"]["ID"]].append(sub_feat)
-
     def with_masked_annotations(
         self, annot_types, mask_char=None, shadow=False, extend_query=False
     ):
@@ -1046,6 +1085,30 @@ class Sequence(_Annotatable, SequenceI):
             seq = str(self)
         return f"{myclass}({seq})"
 
+    def __getitem__(self, index):
+
+        if hasattr(index, "get_slice"):
+            return index.get_slice()
+
+        if isinstance(index, Map):
+            new = self._mapped(index)
+
+        elif isinstance(index, (int, slice)):
+            new = self.__class__(
+                self._seq[index], name=self.name, check=False, info=self.info
+            )
+
+        if self.annotation_db is not None:
+            new.annotation_db = self.annotation_db
+
+        if hasattr(self, "_repr_policy"):
+            new._repr_policy.update(self._repr_policy)
+
+        return new
+
+    def __iter__(self):
+        yield from iter(str(self))
+
     def get_name(self):
         """Return the sequence name -- should just use name instead."""
         return self.name
@@ -1053,8 +1116,12 @@ class Sequence(_Annotatable, SequenceI):
     def __len__(self):
         return len(self._seq)
 
-    def __iter__(self):
-        return iter(self._seq)
+    def __str__(self):
+        result = str(self._seq)
+        if self._seq.reversed:
+            with contextlib.suppress(TypeError):
+                result = self.moltype.complement(result)
+        return result
 
     def gettype(self):  # pragma: no cover
         """Return the sequence type."""
@@ -1247,7 +1314,6 @@ class ProteinWithStopSequence(Sequence):
 class NucleicAcidSequence(Sequence):
     """Abstract base class for DNA and RNA sequences."""
 
-    PROTEIN = None  # will set in moltype
     codon_alphabet = None  # will set in moltype
 
     def reverse_complement(self):
@@ -1257,9 +1323,12 @@ class NucleicAcidSequence(Sequence):
 
     def rc(self):
         """Converts a nucleic acid sequence to its reverse complement."""
-        complement = self.moltype.rc(self)
-        rc = self.__class__(complement, name=self.name, info=self.info)
-        self._annotations_nucleic_reversed_on(rc)
+        rc = self.__class__(
+            self._seq[::-1], name=self.name, check=False, info=self.info
+        )
+        if self.annotation_db is not None:
+            rc.annotation_db = self.annotation_db
+        rc._seq = self._seq[::-1]
         return rc
 
     def has_terminal_stop(self, gc=None, allow_partial=False):
@@ -1310,7 +1379,7 @@ class NucleicAcidSequence(Sequence):
 
         return self.__class__(codons, name=self.name, info=self.info)
 
-    def get_translation(self, gc=None, incomplete_ok=False):
+    def get_translation(self, gc=None, incomplete_ok=False, include_stop=False):
         """translate to amino acid sequence
 
         Parameters
@@ -1326,7 +1395,7 @@ class NucleicAcidSequence(Sequence):
         sequence of PROTEIN moltype
         """
         gc = get_code(gc)
-        codon_alphabet = self.codon_alphabet(gc).with_gap_motif()
+        codon_alphabet = gc.get_alphabet(include_stop=include_stop).with_gap_motif()
         # translate the codons
         translation = []
         for posn in range(0, len(self._seq) - 2, 3):
@@ -1380,11 +1449,11 @@ class NucleicAcidSequence(Sequence):
 
     def to_rna(self):
         """Returns copy of self as RNA."""
-        return RnaSequence(self)
+        return self.to_moltype("rna")
 
     def to_dna(self):
         """Returns copy of self as DNA."""
-        return DnaSequence(self)
+        return self.to_moltype("dna")
 
     def strand_symmetry(self, motif_length=1):
         """returns G-test for strand symmetry"""
@@ -1455,7 +1524,7 @@ def _input_vals_neg_step(seqlen, start, stop, step):
 
 
 class SeqView:
-    __slots__ = ("seq", "start", "stop", "step")
+    __slots__ = ("seq", "start", "stop", "step", "offset")
 
     def __init__(self, seq, *, start: int = None, stop: int = None, step: int = None):
         if step == 0:
@@ -1467,8 +1536,22 @@ class SeqView:
         self.seq = seq
         self.start = start
         self.stop = stop
+        self.offset = None
         self.step = step
 
+    @property
+    def reversed(self):
+        return self.step < 0
+
+    def absolute_index(self, value):
+        # note: this is the positive/positive case...
+
+        if self.start + value > self.stop:
+            raise IndexError("Index out of bounds")
+
+        offset = self.start if self.offset is None else self.offset + self.start
+        return offset + (value * self.step)
+
     def _get_index(self, val):
         if len(self) == 0:
             raise IndexError(val)
@@ -2283,3 +2366,54 @@ class ArrayProteinSequence(ArraySequence):
 class ArrayProteinWithStopSequence(ArraySequence):
     moltype = None  # set to PROTEIN_WITH_STOP in moltype.py
     alphabet = None  # set to PROTEIN_WITH_STOP.alphabets.degen_gapped in moltype.py
+
+
+@singledispatch
+def _coerce_seq(data, preserve_case, checker):
+    from cogent3.core.alignment import Aligned
+
+    if isinstance(data, Aligned):
+        return _coerce_seq(str(data), preserve_case, checker)
+    raise NotImplementedError(f"{type(data)}")
+
+
+@_coerce_seq.register
+def _(data: SeqView, preserve_case, checker):
+    return data
+
+
+@_coerce_seq.register
+def _(data: Sequence, preserve_case, checker):
+    return _coerce_seq(str(data), preserve_case, checker)
+
+
+@_coerce_seq.register
+def _(data: ArraySequence, preserve_case, checker):
+    return _coerce_seq(str(data), preserve_case, checker)
+
+
+@_coerce_seq.register
+def _(data: str, preserve_case, checker):
+    if not preserve_case:
+        data = data.upper()
+    checker(data)
+    return SeqView(data)
+
+
+@_coerce_seq.register
+def _(data: bytes, preserve_case, checker):
+    if not preserve_case:
+        data = data.upper()
+    data = data.decode("utf8")
+    checker(data)
+    return SeqView(data)
+
+
+@_coerce_seq.register
+def _(data: tuple, preserve_case, checker):
+    return _coerce_seq("".join(data), preserve_case, checker)
+
+
+@_coerce_seq.register
+def _(data: list, preserve_case, checker):
+    return _coerce_seq("".join(data), preserve_case, checker)
diff --git a/src/cogent3/evolve/ns_substitution_model.py b/src/cogent3/evolve/ns_substitution_model.py
index d2089c323..c8f6c3712 100644
--- a/src/cogent3/evolve/ns_substitution_model.py
+++ b/src/cogent3/evolve/ns_substitution_model.py
@@ -1,6 +1,6 @@
 import numpy
 
-from cogent3.core import moltype
+from cogent3.core import genetic_code, moltype
 from cogent3.evolve.discrete_markov import PsubMatrixDefn
 from cogent3.evolve.predicate import MotifChange
 from cogent3.maths.optimisers import ParameterOutOfBoundsError
@@ -188,11 +188,11 @@ class NonReversibleTrinucleotide(NonReversibleNucleotide):
 class NonReversibleCodon(_Codon, NonReversibleNucleotide):
     """Base non-reversible codon substitution model."""
 
+    # todo deprecate alphabet argument
     @extend_docstring_from(Parametric.__init__)
     def __init__(self, alphabet=None, gc=None, **kw):
-        if gc is not None:
-            alphabet = moltype.CodonAlphabet(gc=gc)
-        kw["alphabet"] = alphabet or moltype.STANDARD_CODON
+        self.gc = genetic_code.get_code(gc)
+        kw["alphabet"] = self.gc.get_alphabet()
         NonReversibleNucleotide.__init__(self, **kw)
 
 
diff --git a/src/cogent3/evolve/substitution_model.py b/src/cogent3/evolve/substitution_model.py
index dc5af7d97..fbc02ca24 100644
--- a/src/cogent3/evolve/substitution_model.py
+++ b/src/cogent3/evolve/substitution_model.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 """
 substitution_model.py
 
@@ -40,7 +39,7 @@ import numpy
 
 from numpy.linalg import svd
 
-from cogent3.core import moltype
+from cogent3.core import genetic_code, moltype
 from cogent3.evolve import motif_prob_model, parameter_controller, predicate
 from cogent3.evolve.likelihood_tree import make_likelihood_tree_leaf
 from cogent3.evolve.substitution_calculation import (
@@ -1011,7 +1010,7 @@ class _Codon:
             return ndiffs == 1
 
     def get_predefined_predicates(self):
-        codon_preds = _CodonPredicates(self.get_alphabet().get_genetic_code())
+        codon_preds = _CodonPredicates(self.gc)
 
         preds = _TimeReversibleNucleotide.get_predefined_predicates(self)
         preds.update(
@@ -1028,9 +1027,9 @@ class _Codon:
 class TimeReversibleCodon(_Codon, _TimeReversibleNucleotide):
     """Core substitution model for codons"""
 
+    # todo deprecate alphabet argument
     @extend_docstring_from(_TimeReversibleNucleotide.__init__)
     def __init__(self, alphabet=None, gc=None, **kw):
-        if gc is not None:
-            alphabet = moltype.CodonAlphabet(gc=gc)
-        kw["alphabet"] = alphabet or moltype.STANDARD_CODON
+        self.gc = genetic_code.get_code(gc)
+        kw["alphabet"] = self.gc.get_alphabet()
         _TimeReversibleNucleotide.__init__(self, **kw)
diff --git a/src/cogent3/parse/gff.py b/src/cogent3/parse/gff.py
index 4a7ef79a1..5e96d2ca5 100644
--- a/src/cogent3/parse/gff.py
+++ b/src/cogent3/parse/gff.py
@@ -3,8 +3,6 @@ import functools
 import os
 import typing
 
-from pathlib import Path
-
 from cogent3.util.io import open_
 
 
@@ -22,10 +20,15 @@ __maintainer__ = "Peter Maxwell"
 __email__ = "pm67nz@gmail.com"
 __status__ = "Production"
 
+OptionalCallable = typing.Optional[typing.Callable]
+OptionalStrContainer = typing.Optional[typing.Union[str, typing.Sequence]]
+
 
 @functools.singledispatch
 def gff_parser(
-    f: typing.Union[str, os.PathLike, abc.Sequence, typing.IO], attribute_parser=None
+    f: typing.Union[str, os.PathLike, abc.Sequence, typing.IO],
+    attribute_parser: OptionalCallable = None,
+    seqids: OptionalStrContainer = None,
 ) -> typing.Iterable[dict]:
     """parses a gff file
 
@@ -43,23 +46,35 @@ def gff_parser(
     dict
         contains each of the 9 parameters specified by gff3, and comments.
     """
-    yield from _gff_parser(f, attribute_parser=attribute_parser)
+    yield from _gff_parser(f, attribute_parser=attribute_parser, seqids=seqids)
 
 
 @gff_parser.register
-def _(f: str, attribute_parser=None) -> typing.Iterable[dict]:
+def _(
+    f: str,
+    attribute_parser: OptionalCallable = None,
+    seqids: OptionalStrContainer = None,
+) -> typing.Iterable[dict]:
     with open_(f) as infile:
-        yield from _gff_parser(infile, attribute_parser=attribute_parser)
+        yield from gff_parser(infile, attribute_parser=attribute_parser, seqids=seqids)
 
 
 @gff_parser.register
-def _(f: os.PathLike, attribute_parser=None) -> typing.Iterable[dict]:
+def _(
+    f: os.PathLike,
+    attribute_parser: OptionalCallable = None,
+    seqids: OptionalStrContainer = None,
+) -> typing.Iterable[dict]:
     with open_(f) as infile:
-        yield from _gff_parser(infile, attribute_parser=attribute_parser)
+        yield from gff_parser(infile, attribute_parser=attribute_parser, seqids=seqids)
 
 
-def _gff_parser(f, attribute_parser=None):
+def _gff_parser(
+    f, attribute_parser: OptionalCallable = None, seqids: OptionalStrContainer = None
+) -> typing.Iterable[dict]:
     """parses a gff file"""
+    seqids = seqids or set()
+    seqids = {seqids} if isinstance(seqids, str) else set(seqids)
 
     gff3_header = "gff-version 3"
     if isinstance(f, list):
@@ -81,13 +96,16 @@ def _gff_parser(f, attribute_parser=None):
         if not line:
             continue
 
-        cols = line.split("\t")
+        cols = [c.strip() for c in line.split("\t")]
         # the final column (attributes) may be empty
         if len(cols) == 8:
             cols.append("")
         assert len(cols) == 9, len(line)
         seqid, source, type_, start, end, score, strand, phase, attributes = cols
 
+        if seqids and seqid not in seqids:
+            continue
+
         # adjust for 0-based indexing
         start, end = int(start) - 1, int(end)
         # start is always meant to be less than end in GFF
@@ -103,7 +121,7 @@ def _gff_parser(f, attribute_parser=None):
         # all attributes have an "ID" but this may not be unique
         attributes = attribute_parser(attributes, (start, end))
 
-        rtn = {
+        yield {
             "SeqID": seqid,
             "Source": source,
             "Type": type_,
@@ -115,10 +133,9 @@ def _gff_parser(f, attribute_parser=None):
             "Attributes": attributes,
             "Comments": comments,
         }
-        yield rtn
 
 
-def parse_attributes_gff2(attributes, span):
+def parse_attributes_gff2(attributes: str, span: typing.Tuple[int, int]) -> dict:
     """Returns a dict with name and info keys"""
     name = attributes[attributes.find('"') + 1 :]
     if '"' in name:
@@ -126,7 +143,7 @@ def parse_attributes_gff2(attributes, span):
     return {"ID": name, "Info": attributes}
 
 
-def parse_attributes_gff3(attributes, span):
+def parse_attributes_gff3(attributes: str, span: typing.Tuple[int, int]) -> dict:
     """Returns a dictionary containing all the attributes"""
     attributes = attributes.strip(";")
     attributes = attributes.split(";")
diff --git a/src/cogent3/util/deserialise.py b/src/cogent3/util/deserialise.py
index 3c3daddae..fc51d4f26 100644
--- a/src/cogent3/util/deserialise.py
+++ b/src/cogent3/util/deserialise.py
@@ -5,9 +5,7 @@ from importlib import import_module
 
 import cogent3
 
-from cogent3.core.alignment import Aligned
 from cogent3.core.genetic_code import get_code
-from cogent3.core.moltype import _CodonAlphabet, get_moltype
 from cogent3.util.io import open_, path_exists
 
 
@@ -161,6 +159,8 @@ def deserialise_result(data):
 @register_deserialiser("cogent3.core.moltype")
 def deserialise_moltype(data):
     """returns a cogent3 MolType instance, or a CodonAlphabet"""
+    from cogent3.core.moltype import _CodonAlphabet, get_moltype
+
     data.pop("version", None)
     label = data["moltype"]
     data["moltype"] = get_moltype(label)
@@ -178,6 +178,8 @@ def deserialise_moltype(data):
 @register_deserialiser("cogent3.core.alphabet")
 def deserialise_alphabet(data):
     """returns a cogent3 Alphabet instance"""
+    from cogent3.core.moltype import _CodonAlphabet, get_moltype
+
     data.pop("version", None)
     if _get_class(data.get("type")) == _CodonAlphabet:
         result = deserialise_moltype(data)
@@ -207,6 +209,7 @@ def deserialise_seq(data, aligned=False):
     -------
 
     """
+    from cogent3.core.alignment import Aligned
     from cogent3.core.moltype import get_moltype
 
     data.pop("version", None)
diff --git a/src/cogent3/util/misc.py b/src/cogent3/util/misc.py
index c66a7a69e..4688491c2 100644
--- a/src/cogent3/util/misc.py
+++ b/src/cogent3/util/misc.py
@@ -117,6 +117,10 @@ def bytes_to_string(data):
     """returns a string if data is bytes, otherwise returns original"""
     if isinstance(data, bytes):
         data = data.decode("utf_8")
+    elif isinstance(data, (list, tuple)):
+        data = "".join(str(d) for d in data)
+    else:
+        data = str(data)
     return data
 
 


Test patch:
diff --git a/tests/data/simple.gff b/tests/data/simple.gff
new file mode 100644
index 000000000..2748c1e8b
--- /dev/null
+++ b/tests/data/simple.gff
@@ -0,0 +1,6 @@
+##gff-version 3
+test_seq	.	exon	2	10	.	+	.	ID=exon1
+test_seq	.	CDS 	3	9	.	+	.	ID=cds1
+test_seq	.	exon	12	20	.	+	.	ID=exon2
+test_seq	.	CDS	    4	10	.	+	.	ID=cds2
+test_seq	.	CpG 	3	12	.	+	.	ID=CpG1
diff --git a/tests/data/simple2.gff b/tests/data/simple2.gff
new file mode 100644
index 000000000..641aa020d
--- /dev/null
+++ b/tests/data/simple2.gff
@@ -0,0 +1,6 @@
+##gff-version 3
+test_seq	.	exon	2	10	.	+	.	ID=exon3;parent=cds1
+test_seq	.	CDS 	3	9	.	+	.	ID=cds3
+test_seq	.	exon	12	20	.	+	.	ID=exon4;parent=cds1
+test_seq	.	CDS	    4	10	.	+	.	ID=cds4
+test_seq	.	CpG 	3	12	.	+	.	ID=CpG2
diff --git a/tests/test_core/test_alignment.py b/tests/test_core/test_alignment.py
index fc1f0fd22..39bdcb16c 100644
--- a/tests/test_core/test_alignment.py
+++ b/tests/test_core/test_alignment.py
@@ -45,7 +45,6 @@ from cogent3.core.alignment import (
     seqs_from_kv_pairs,
 )
 from cogent3.core.alphabet import AlphabetError
-from cogent3.core.annotation import Feature
 from cogent3.core.moltype import AB, ASCII, BYTES, DNA, PROTEIN, RNA
 from cogent3.core.sequence import ArraySequence, RnaSequence, Sequence
 from cogent3.maths.util import safe_p_log_p
@@ -922,7 +921,7 @@ class SequenceCollectionBaseTests(object):
 
         if self.Class == Alignment:
             aln_seq_3 = aln.get_seq("seq3")
-            matches = [m for m in aln_seq_3.get_annotations_matching("*")]
+            matches = [m for m in aln_seq_3.get_features_matching("*")]
             self.assertFalse("-" in matches[0].get_slice())
 
     def test_annotate_from_gff3(self):
@@ -942,14 +941,14 @@ class SequenceCollectionBaseTests(object):
         aln_seq = aln.named_seqs[name]
         if not hasattr(aln_seq, "annotations"):
             aln_seq = aln_seq.data
-        matches = [m for m in aln_seq.get_annotations_matching("*", extend_query=True)]
+        matches = [m for m in aln_seq.get_features_matching("*", extend_query=True)]
         # 13 features with one having 2 parents, so 14 instances should be found
         self.assertEqual(len(matches), 14)
-        matches = [m for m in aln_seq.get_annotations_matching("gene")]
+        matches = [m for m in aln_seq.get_features_matching("gene")]
         self.assertEqual(len(matches), 1)
-        matches = matches[0].get_annotations_matching("mRNA")
+        matches = matches[0].get_features_matching("mRNA")
         self.assertEqual(len(matches), 1)
-        matches = matches[0].get_annotations_matching("exon")
+        matches = matches[0].get_features_matching("exon")
         self.assertEqual(len(matches), 3)
 
     def test_add(self):
@@ -1430,7 +1429,7 @@ class SequenceCollectionBaseTests(object):
 
         if self.Class == SequenceCollection:
             # this class cannot slice
-            return
+            return True
 
         # should persist in slicing
         self.assertEqual(
@@ -1441,7 +1440,7 @@ class SequenceCollectionBaseTests(object):
         """the wrap argument affects the number of columns"""
         if self.Class == SequenceCollection:
             # this class does not have this method
-            return
+            return True
 
         # indirectly tested via counting number of occurrences of 'class="label"'
         seqs = self.Class({"a": "AAAAA", "b": "AAA--"})
@@ -2908,25 +2907,15 @@ class AlignmentTests(AlignmentBaseTests, TestCase):
         # generates an annotatable Alignment object
         aln = load_aligned_seqs(path, array_align=False, moltype="dna")
         # when the annotation is outside(before) boundary of the slice
-        aln.named_seqs["NineBande"].data.add_annotation(
-            Feature, "exon", "annot1", [(0, 10)]
-        )
+        aln.named_seqs["NineBande"].data.add_feature("exon", "annot1", [(0, 10)])
         # when the annotation is across boundary of the slice
-        aln.named_seqs["Mouse"].data.add_annotation(
-            Feature, "exon", "annot2", [(10, 21)]
-        )
+        aln.named_seqs["Mouse"].data.add_feature("exon", "annot2", [(10, 21)])
         # when the annotation is within boundary of the slice
-        aln.named_seqs["Human"].data.add_annotation(
-            Feature, "exon", "annot3", [(20, 25)]
-        )
+        aln.named_seqs["Human"].data.add_feature("exon", "annot3", [(20, 25)])
         # when the annotation is across boundary of the slice
-        aln.named_seqs["HowlerMon"].data.add_annotation(
-            Feature, "exon", "annot4", [(25, 32)]
-        )
+        aln.named_seqs["HowlerMon"].data.add_feature("exon", "annot4", [(25, 32)])
         # when the annotation is outside(after) boundary of the slice
-        aln.named_seqs["DogFaced"].data.add_annotation(
-            Feature, "exon", "annot5", [(40, 45)]
-        )
+        aln.named_seqs["DogFaced"].data.add_feature("exon", "annot5", [(40, 45)])
         aln = aln[20:30]
 
         # for these species, each has an annotation spanning slice boundary or within it
@@ -2955,9 +2944,7 @@ class AlignmentTests(AlignmentBaseTests, TestCase):
             self.assertEqual(len(new_seq.data.annotations), 1)
 
         # add another human annotation that is outside slice
-        aln.named_seqs["Human"].data.add_annotation(
-            Feature, "exon", "annot6", [(40, 45)]
-        )
+        aln.named_seqs["Human"].data.add_feature("exon", "annot6", [(40, 45)])
         # tests the case when sliced argument if False regarding the Human sequence
         new_seq = aln.named_seqs["Human"].deepcopy(sliced=False)
         self.assertEqual(len(new_seq.data), len(aln.named_seqs["Human"].data))
@@ -2995,9 +2982,9 @@ class AlignmentTests(AlignmentBaseTests, TestCase):
         s1 = Sequence("TTTTTTAAAA", name="test_seq1")
         s2 = Sequence("AAAATTTTTT", name="test_seq2")
         s3 = Sequence("AATTTTTAAA", name="test_seq3")
-        s1.add_annotation(Feature, "exon", "fred", [(0, 6)])
-        s2.add_annotation(Feature, "exon", "fred", [(4, 10)])
-        s3.add_annotation(Feature, "exon", "fred", [(2, 7)])
+        s1.add_feature("exon", "fred", [(0, 6)])
+        s2.add_feature("exon", "fred", [(4, 10)])
+        s3.add_feature("exon", "fred", [(2, 7)])
         data = {"seq1": s1, "seq2": s2, "seq3": s3}
         aln = self.Class(data=data)
         aln.add_feature("demo", "one", [(0, 1), (2, 4)])
@@ -3018,9 +3005,9 @@ class AlignmentTests(AlignmentBaseTests, TestCase):
         """get_annotations_from_any_seq returns correct annotations"""
         data = {"seq1": "ACGTACGTA", "seq2": "ACCGAA---", "seq3": "ACGTACGTT"}
         seqs = self.Class(data, moltype=DNA)
-        seqs.get_seq("seq1").add_annotation(Feature, "exon", "annotation1", [(3, 8)])
-        seqs.get_seq("seq2").add_annotation(Feature, "exon", "annotation2", [(1, 2)])
-        seqs.get_seq("seq3").add_annotation(Feature, "exon", "annotation3", [(3, 6)])
+        seqs.get_seq("seq1").add_feature("exon", "annotation1", [(3, 8)])
+        seqs.get_seq("seq2").add_feature("exon", "annotation2", [(1, 2)])
+        seqs.get_seq("seq3").add_feature("exon", "annotation3", [(3, 6)])
         got = seqs.get_annotations_from_any_seq()
         self.assertEqual(len(got), 3)
         self.assertEqual(str(got[0]), 'exon "annotation1" at [3:8]/9')
@@ -3047,7 +3034,7 @@ class AlignmentTests(AlignmentBaseTests, TestCase):
 
         data = {"seq1": "ACGTACGTA", "seq2": "ACCGAA---", "seq3": "ACGTACGTT"}
         seqs = self.Class(data, moltype=DNA)
-        x = seqs.get_seq("seq1").add_annotation(Feature, "exon", "fred", [(3, 8)])
+        x = seqs.get_seq("seq1").add_feature("exon", "fred", [(3, 8)])
         expect = str(x.get_slice())
         new = seqs.rename_seqs(lambda x: x.upper())
         got = list(new.get_annotations_from_any_seq("exon"))[0]
@@ -3396,6 +3383,21 @@ def test_upac_consensus_allow_gaps(moltype, array_align):
     assert iupac == "ACGG"
 
 
+def test_rc_iter():
+    dna = DNA.make_seq("ACG", name="seq1")
+    rc = dna.rc()
+
+    got = [x for x in rc]
+    expect = ["C", "G", "T"]
+    assert got == expect
+
+
+def test_to_dna_raises():
+    seq = make_seq("ETV", moltype="protein")
+    with pytest.raises(AlphabetError):
+        seq.to_moltype("dna")
+
+
 # run tests if invoked from command line
 if __name__ == "__main__":
     main()
diff --git a/tests/test_core/test_alphabet.py b/tests/test_core/test_alphabet.py
index ec334a3dd..8a1df9263 100644
--- a/tests/test_core/test_alphabet.py
+++ b/tests/test_core/test_alphabet.py
@@ -197,6 +197,30 @@ class EnumerationTests(TestCase):
         x = a * c
         self.assertEqual(x.gap, None)
 
+    def test_counts(self):
+        """Enumeration counts should count freqs in array"""
+        a = DnaBases
+        f = array([[0, 0, 1, 0, 0, 3]])
+        assert_equal(a.counts(f), array([4, 1, 0, 1]))
+        # check that it works with byte array
+        f = array([[0, 0, 1, 0, 0, 3]], "B")
+        assert_equal(a.counts(f), array([4, 1, 0, 1]))
+        # should ignore out-of-bounds items
+        g = [0, 4]
+        assert_equal(a.counts(g), array([1, 0, 0, 0]))
+        # make sure it works for long sequences, i.e. no wraparound at 255
+        h = [0, 3] * 70000
+        assert_equal(a.counts(h), array([70000, 0, 0, 70000]))
+        h2 = array(h).astype("B")
+        assert_equal(a.counts(h2), array([70000, 0, 0, 70000]))
+        i = array([0, 3] * 75000)
+        assert_equal(a.counts(i), array([75000, 0, 0, 75000]))
+        # make sure it works for long _binary_ sequences, e.g. the results
+        # of array comparisons.
+        a = array([0, 1, 2, 3] * 10000)
+        b = array([0, 0, 0, 0] * 10000)
+        same = a == b
+
 
 class CharAlphabetTests(TestCase):
     """Tests of CharAlphabets."""
diff --git a/tests/test_core/test_annotation.py b/tests/test_core/test_annotation.py
index b4df91056..7f590374f 100644
--- a/tests/test_core/test_annotation.py
+++ b/tests/test_core/test_annotation.py
@@ -23,10 +23,10 @@ def makeSampleSequence(with_gaps=False):
     cds = (15, 25)
     utr = (12, 15)
     if with_gaps:
-        raw_seq = raw_seq[:5] + "-----" + raw_seq[10:-2] + "--"
+        raw_seq = f"{raw_seq[:5]}-----{raw_seq[10:-2]}--"
     seq = DNA.make_seq(raw_seq)
-    seq.add_annotation(Feature, "CDS", "CDS", [cds])
-    seq.add_annotation(Feature, "5'UTR", "5' UTR", [utr])
+    seq.add_feature("CDS", "CDS", [cds])
+    seq.add_feature("5'UTR", "5' UTR", [utr])
     return seq
 
 
@@ -119,9 +119,8 @@ class TestAnnotations(unittest.TestCase):
             "LTR": {"FAKE01": "CCCTTTTT", "FAKE02": "CCCTTTTT"},
         }
         newaln = self.aln[:5] + self.aln[10:]
-        feature_list = newaln.get_annotations_matching("LTR")
         for annot_type in ["LTR", "misc_feature", "CDS", "5'UTR"]:
-            feature_list = newaln.get_annotations_matching(annot_type)
+            feature_list = newaln.get_features_matching(annot_type)
             new = newaln.get_region_covering_all(feature_list).get_slice().to_dict()
             expected = aln_expecteds[annot_type]
             assert expected == new, (annot_type, expected, new)
@@ -142,7 +141,7 @@ class TestAnnotations(unittest.TestCase):
 
     def test_feature_projection(self):
         expecteds = {"FAKE01": "CCCAAAATTTTTT", "FAKE02": "CCC-----TTTTT"}
-        aln_ltr = self.aln.get_annotations_matching("LTR")[0]
+        aln_ltr = self.aln.get_features_matching("LTR")[0]
         for seq_name in ["FAKE01", "FAKE02"]:
             expected = expecteds[seq_name]
             seq_ltr = self.aln.project_annotation(seq_name, aln_ltr)
diff --git a/tests/test_core/test_annotation_db.py b/tests/test_core/test_annotation_db.py
new file mode 100644
index 000000000..b259950de
--- /dev/null
+++ b/tests/test_core/test_annotation_db.py
@@ -0,0 +1,411 @@
+import pathlib
+
+import numpy
+import pytest
+
+from cogent3 import DNA, RNA, SequenceCollection, _Table, load_seq
+from cogent3.core.annotation import FeatureNew
+from cogent3.core.annotation_db import (
+    GffAnnotationDb,
+    SupportsFeatures,
+    load_annotations,
+)
+from cogent3.core.sequence import Sequence
+
+
+__author__ = "Gavin Huttley"
+__copyright__ = "Copyright 2007-2022, The Cogent Project"
+__credits__ = ["Gavin Huttley", "Katherine Caley"]
+__license__ = "BSD-3"
+__version__ = "2023.2.12a1"
+__maintainer__ = "Gavin Huttley"
+__email__ = "gavin.huttley@anu.edu.au"
+__status__ = "Production"
+
+DATA_DIR = pathlib.Path(__file__).parent.parent / "data"
+
+
+@pytest.fixture(scope="session")
+def gff_db():
+    paths = (
+        "sample_data/Homo_sapiens.GRCh38.109.chromosome.22.gff3.gz",
+        "data/c_elegans_WS199_shortened_gff.gff3",
+    )
+    return load_annotations(paths[1])
+
+
+@pytest.fixture()
+def seq_db():
+    seq = load_seq(DATA_DIR / "c_elegans_WS199_dna_shortened.fasta", moltype="dna")
+    db = load_annotations(DATA_DIR / "c_elegans_WS199_shortened_gff.gff3")
+
+    seq.annotation_db = db
+
+    return seq
+
+
+@pytest.fixture()
+def seq() -> Sequence:
+    return Sequence("ATTGTACGCCTTTTTTATTATT", name="test_seq")
+
+
+@pytest.fixture()
+def anno_db() -> GffAnnotationDb:
+    # an empty db that we can add to
+    return GffAnnotationDb([])
+
+
+@pytest.fixture()
+def simple_seq_gff_db() -> Sequence:
+    seq = Sequence("ATTGTACGCCTTTTTTATTATT", name="test_seq")
+    seq.annotate_from_gff(DATA_DIR / "simple.gff")
+    return seq
+
+
+def test_gff_describe(gff_db):
+    result = gff_db.describe
+    assert isinstance(result, _Table)
+
+
+def test_gff_features_matching(gff_db):
+    result = list(gff_db.get_features_matching(biotype="CDS"))
+    assert len(result)
+
+
+def test_gff_get_children(gff_db):
+    # an ID with 8 children
+    children = list(gff_db.get_feature_children(name="Transcript:B0019.1"))
+    assert len(children) == 8
+    # this is parent of the transcript, and due to the funky structure of
+    # this gff, also to one of the CDS entries
+    children = list(gff_db.get_feature_children(name="Gene:WBGene00000138"))
+    assert len(children) == 2
+
+
+@pytest.mark.parametrize(
+    "name,expected",
+    (
+        ("CDS:B0019.1", ("Transcript:B0019.1", "Gene:WBGene00000138")),
+        ("Transcript:B0019.1", ("Gene:WBGene00000138",)),
+    ),
+)
+def test_gff_get_parent(gff_db, name, expected):
+    # the first id has two parents, which is weird
+    got = list(gff_db.get_feature_parent(name=name))
+    assert len(got) == len(expected)
+    assert {g["name"] for g in got} == set(expected)
+
+
+def test_gff_counts(gff_db):
+    got = gff_db.biotype_counts()
+    assert len(got) > 0
+
+
+def test_gff_find_user_features(gff_db):
+    record = dict(
+        seqid="2", name="gene-01", biotype="gene", spans=[(23, 33)], strand="+"
+    )
+    gff_db.add_feature(**record)
+    # by biotype
+    found = any(
+        gene["name"] == "gene-01"
+        for gene in gff_db.get_features_matching(biotype="gene")
+    )
+    assert found
+    # by name
+    found = any(
+        gene["name"] == "gene-01"
+        for gene in gff_db.get_features_matching(name="gene-01")
+    )
+    assert found
+
+
+def test_to_rich_dict(gff_db):
+    # not a very good test! needs to be serialise / deserialize
+    rd = gff_db.to_rich_dict()
+    assert "type" in rd
+    _ = GffAnnotationDb.from_dict(rd)
+
+
+def test_empty_data():
+    _ = GffAnnotationDb([])
+
+
+# testing GenBank files
+@pytest.fixture(scope="session")
+def gb_db():
+    paths = ("data/annotated_seq.gb",)
+    return load_annotations(paths[0])
+
+
+@pytest.mark.parametrize("parent_biotype, name", (("gene", "CNA00110"),))
+def test_gb_get_children(gb_db, parent_biotype, name):
+    parent = list(gb_db.get_features_matching(biotype=parent_biotype, name=name))[0]
+    coords = numpy.array(parent["spans"])
+    child = list(
+        gb_db.get_feature_children(
+            name=name,
+            exclude_biotype=parent_biotype,
+            start=coords.min(),
+            end=coords.max(),
+        )
+    )[0]
+    assert child["biotype"] != parent["biotype"]
+    assert child["name"] == parent["name"]
+
+
+def test_gb_get_parent(gb_db):
+    cds_id = "CNA00110"
+    cds = list(gb_db.get_features_matching(biotype="CDS", name=cds_id))[0]
+    coords = numpy.array(cds["spans"])
+    parent = list(
+        gb_db.get_feature_parent(
+            name=cds_id,
+            exclude_biotype="CDS",
+            start=coords.min(),
+            end=coords.max(),
+        )
+    )[0]
+    assert parent["biotype"] != cds["biotype"]
+    assert parent["biotype"] == "gene"
+    assert parent["name"] == cds["name"]
+
+
+def test_protocol_adherence(gff_db, gb_db):
+    for db in (gff_db, gb_db):
+        assert isinstance(db, SupportsFeatures)
+
+
+def test_get_features_matching_no_annotation_db(seq):
+    """
+    Test that `get_features_matching` returns an empty list when no annotation database is attached to the sequence.
+    """
+    assert not list(seq.get_features_matching())
+
+
+def test_get_features_matching_no_matching_feature(seq, anno_db):
+    """
+    Test that `get_features_matching` returns an empty list when there is no matching feature in the annotation database.
+    """
+    seq.annotation_db = anno_db
+    anno_db.add_feature(
+        seqid=seq.name, biotype="exon", name="exon1", spans=[(1, 4)], strand="+"
+    )
+
+    assert not list(seq.get_features_matching(feature_type="exon", name="non_matching"))
+    assert not list(seq.get_features_matching(feature_type="CDS"))
+
+
+def test_get_features_matching_matching_feature(seq, anno_db):
+    """
+    Test that `get_features_matching` returns a list with one matching feature in the annotation database.
+    """
+    seq.annotation_db = anno_db
+    anno_db.add_feature(
+        seqid=seq.name, biotype="exon", name="exon1", spans=[(1, 4)], strand="+"
+    )
+    got = list(seq.get_features_matching(feature_type="exon"))
+
+    assert got[0].biotype == "exon"
+    assert got[0].name == "exon1"
+    assert len(got) == 1
+    assert got[0].reversed == False
+
+
+def test_get_features_matching_matching_features(anno_db: GffAnnotationDb, seq):
+    """
+    Test that `get_features_matching` returns a list with all matching features in the annotation database.
+    """
+    seq.annotation_db = anno_db
+
+    anno_db.add_feature(
+        seqid=seq.name, biotype="exon", name="exon1", spans=[(1, 4)], strand="+"
+    )
+    anno_db.add_feature(
+        seqid=seq.name, biotype="exon", name="exon2", spans=[(6, 10)], strand="+"
+    )
+    got = list(seq.get_features_matching(feature_type="exon"))
+
+    assert len(got) == 2
+
+
+def test_annotate_from_gff(seq):
+    seq.annotate_from_gff(DATA_DIR / "simple.gff")
+
+    got = list(seq.get_features_matching(feature_type="exon"))
+    assert len(got) == 2
+
+    got = list(seq.get_features_matching(feature_type="CDS"))
+    assert len(got) == 2
+
+    got = list(seq.get_features_matching(feature_type="CpG"))
+    assert len(got) == 1
+
+
+def test_get_features_matching_start_stop(seq):
+    # todo: need to implement LostSpans()
+    seq.annotate_from_gff(DATA_DIR / "simple.gff")
+    got = list(seq.get_features_matching(start=2, stop=10))
+    assert len(got) == 4
+
+@pytest.mark.xfail(reason="discuss with kath")
+def test_get_features_matching_start_stop_seqview(seq):
+    """testing that get_features_matching adjusts"""
+    seq.annotate_from_gff(DATA_DIR / "simple.gff")
+    seq_features = list(seq.get_features_matching(start=0, stop=3))
+    assert len(seq_features) == 3
+
+    # edge case, only 1 features that overlaps with index 12
+    # is actually returning [exon2 at [11:20]/13, CpG1 at [2:12]/13]
+    # possibly a bug in the SQL generating code
+    subseq = seq[9:]
+    seq_features_features = list(subseq.get_features_matching(start=3, stop=10))
+    assert len(seq_features_features) == 1
+
+
+def test_get_slice():
+    """get_slice should return the same as slicing the sequence directly"""
+    seq = Sequence("ATTGTACGCCCCTGA", name="test_seq")
+    feature_dict = {
+        "biotype": "CDS",
+        "name": "fake",
+        "spans": [
+            (5, 10),
+        ],
+        "reversed": False,
+    }
+
+    feature = FeatureNew(seq, **feature_dict)
+    got = feature.get_slice()
+    assert str(got) == str(seq[5:10])
+
+
+def test_feature_get_children(seq_db):
+    feat = list(seq_db.get_features_matching(name="Transcript:B0019.1"))[0]
+    new_feat_5pUTR = list(feat.get_children(biotype="five_prime_UTR"))
+    assert len(new_feat_5pUTR) == 1
+
+    new_feat_CDS = list(feat.get_children(biotype="CDS"))[0]
+    assert new_feat_CDS.name == "CDS:B0019.1"
+
+
+def test_db_persists_post_rc(seq_db):
+    """assert that the db persists after the .rc() method call"""
+    rc_seq = seq_db.rc()
+    assert rc_seq.annotation_db is not None
+
+
+def test_rc_get_slice_negative_feature(seq_db):
+    """given a feature on the - strand, the feature.get_slice() should return
+    the same sequence before and after the sequence is reverse complemented
+    """
+
+    feat = list(seq_db.get_features_matching(name="Transcript:B0019.1"))[0]
+    rc_seq = seq_db.rc()
+    r_feat = list(rc_seq.get_features_matching(name="Transcript:B0019.1"))[0]
+
+    assert feat.get_slice() == r_feat.get_slice()
+
+
+def test_rc_get_slice_positive_feature(anno_db):
+    """given a feature on the + strand, the feature.get_slice() should return
+    the same sequence before and after the sequence is reverse complemented
+    """
+
+    seq = DNA.make_seq("AAAAGGGG", name="seq1")
+
+    seq.annotation_db = anno_db
+    anno_db.add_feature(
+        seqid=seq.name, biotype="exon", name="exon1", spans=[(2, 6)], strand="+"
+    )
+
+    feat = list(seq.get_features_matching(name="exon1"))[0]
+    r_seq = seq.rc()
+    r_feat = list(r_seq.get_features_matching(name="exon1"))[0]
+
+    assert feat.get_slice() == r_feat.get_slice()
+
+
+def test_add_feature(seq):
+    """Sequence supports manual adding of features for a seq with no bound AnnotationDb"""
+    record = dict(name="gene-01", biotype="gene", spans=[(12, 16)], strand="+")
+    seq.add_feature(**record)
+    feats = list(seq.get_features_matching(feature_type="gene"))
+
+    assert seq.annotation_db is not None
+    assert len(feats) == 1
+    assert feats[0].biotype == "gene"
+
+
+def test_add_feature_existing_db(simple_seq_gff_db):
+    """Sequence supports manual adding of features for a seq with an existing AnnotationDb"""
+    record = dict(name="gene-01", biotype="gene", spans=[(12, 16)], strand="+")
+    simple_seq_gff_db.add_feature(**record)
+
+    # total features should be 5+1=6
+    all_feats = list(simple_seq_gff_db.get_features_matching())
+    assert len(all_feats) == 6
+
+
+def test__getitem__(simple_seq_gff_db):
+    """Sequence.__getitem__ should keep the underlying seq in the SeqView
+    and preserve any annotation_db"""
+
+    seq_sliced = simple_seq_gff_db[4:6]
+    assert seq_sliced == str(simple_seq_gff_db)[4:6]
+    # check the underlying seq is still the original sequence data
+    assert seq_sliced._seq.seq == str(simple_seq_gff_db)
+    # check the annotation_db is still attached and the same instance
+    assert (
+        seq_sliced.annotation_db
+        and seq_sliced.annotation_db is simple_seq_gff_db.annotation_db
+    )
+
+
+def test_to_moltype_dna():
+    """to_moltype("dna") ensures conversion from T to U"""
+    seq = DNA.make_seq("AAAAGGGGTTT", name="seq1")
+    rna = seq.to_moltype("rna")
+
+    assert "T" not in rna
+
+
+def test_to_moltype_rna():
+    """to_moltype("rna") ensures conversion from U to T"""
+    seq = RNA.make_seq("AAAAGGGGUUU", name="seq1")
+    rna = seq.to_moltype("dna")
+
+    assert "U" not in rna
+
+
+def test_annotate_from_gff_multiple_calls(seq):
+    """5 records in each gff file, total features on seq should be 10"""
+    seq.annotate_from_gff(DATA_DIR / "simple.gff")
+    seq.annotate_from_gff(DATA_DIR / "simple2.gff")
+    assert len(list(seq.get_features_matching())) == 10
+
+
+def test_sequence_collection_annotate_from_gff():
+    """providing a seqid to SequenceCollection.annotate_from_gff will
+    annotate the SequenceCollection, and the Sequence. Both of these will point
+    to the same AnnotationDb instance
+    """
+    seqs = {"test_seq": "ATCGATCGATCG", "test_seq2": "GATCGATCGATC"}
+    seq_collection = SequenceCollection(seqs)
+    seq_collection.annotate_from_gff(DATA_DIR / "simple.gff", seq_ids="test_seq")
+
+    # the seq for which the seqid was provided is annotated
+    seq = seq_collection.get_seq("test_seq")
+    assert seq_collection.get_seq("test_seq").annotation_db is not None
+    assert len(list(seq_collection.get_seq("test_seq").get_features_matching())) == 5
+    # the seq for which the seqid was NOT provided is NOT annotated
+    assert seq_collection.get_seq("test_seq2").annotation_db is None
+    # the annotation_db on the seq and the seq collection are the same object
+    assert (
+        seq_collection.get_seq("test_seq").annotation_db is seq_collection.annotation_db
+    )
+    got = list(seq.get_features_matching(feature_type="CDS"))
+    assert len(got) == 2
+
+    got = list(seq.get_features_matching(feature_type="CpG"))
+    assert len(got) == 1
diff --git a/tests/test_core/test_core_standalone.py b/tests/test_core/test_core_standalone.py
index 862fe1a4f..33f16017e 100644
--- a/tests/test_core/test_core_standalone.py
+++ b/tests/test_core/test_core_standalone.py
@@ -994,5 +994,37 @@ class SequenceTestMethods(unittest.TestCase):
         self.assertFalse(seq.has_terminal_stop(allow_partial=True))
 
 
+def test_load_seq_new():
+    """load single sequence"""
+    from cogent3 import Sequence
+
+    paths = (
+        "c_elegans_WS199_dna_shortened.fasta",
+        "annotated_seq.gb",
+        "brca1_5.250.paml",
+    )
+    seq_names = ("I", "AE017341", "NineBande")
+    data_dir = pathlib.Path(data_path)
+
+    for i, path in enumerate(paths):
+        got = load_seq(data_dir / path)
+        assert isinstance(got, Sequence)
+        assert got.info.source == str(data_dir / path)
+        assert got.name == seq_names[i]
+
+    # annotated with gb
+    got = load_seq(data_dir / "annotated_seq.gb")
+    assert isinstance(got, Sequence)
+    assert got.annotation_db is not None
+
+    # annotated with gff
+    got = load_seq(
+        data_dir / "annotated_seq.gb",
+        annotation_path=data_dir / "c_elegans_WS199_shortened_gff.gff3",
+    )
+    assert isinstance(got, Sequence)
+    assert got.annotation_db is not None
+
+
 if __name__ == "__main__":
     unittest.main()
diff --git a/tests/test_core/test_features.py b/tests/test_core/test_features.py
index 9bed39b37..44358b34d 100755
--- a/tests/test_core/test_features.py
+++ b/tests/test_core/test_features.py
@@ -24,8 +24,8 @@ class FeaturesTest(TestCase):
         self.s = DNA.make_seq(
             "AAGAAGAAGACCCCCAAAAAAAAAATTTTTTTTTTAAAAAAAAAAAAA", name="Orig"
         )
-        self.exon1 = self.s.add_annotation(Feature, "exon", "fred", [(10, 15)])
-        self.exon2 = self.s.add_annotation(Feature, "exon", "trev", [(30, 40)])
+        self.exon1 = self.s.add_feature("exon", "fred", [(10, 15)])
+        self.exon2 = self.s.add_feature("exon", "trev", [(30, 40)])
         self.nested_feature = self.exon1.add_feature("repeat", "bob", [(2, 5)])
 
     def test_exon_extraction(self):
@@ -37,14 +37,14 @@ class FeaturesTest(TestCase):
         self.assertEqual(str(self.s[self.exon1]), "CCCCC")
         self.assertEqual(str(self.exon1.get_slice()), "CCCCC")
 
-    def test_get_annotations_matching(self):
+    def test_get_features_matching(self):
         """correctly identifies all features of a given type"""
 
         # Usually the only way to get a Feature object like exon1
         # is to ask the sequence for it. There is one method for querying
         # annotations by type and optionally by name:
 
-        exons = self.s.get_annotations_matching("exon")
+        exons = self.s.get_features_matching("exon")
         self.assertEqual(
             str(exons), '[exon "fred" at [10:15]/48, exon "trev" at [30:40]/48]'
         )
@@ -55,30 +55,30 @@ class FeaturesTest(TestCase):
         seq = DNA.make_seq("AAAAAAAAA", name="x")
         exon = seq.add_annotation(Feature, "exon", "fred", [(3, 8)])
         nested_exon = exon.add_annotation(Feature, "exon", "fred", [(3, 7)])
-        exons = seq.get_annotations_matching("exon", extend_query=True)
+        exons = seq.get_features_matching("exon", extend_query=True)
         self.assertEqual(len(exons), 2)
         self.assertEqual(str(exons), '[exon "fred" at [3:8]/9, exon "fred" at [3:7]/5]')
         # tests multiple layers of nested annotations
         nested_exon.add_annotation(Feature, "exon", "fred", [(3, 6)])
-        exons = seq.get_annotations_matching("exon", extend_query=True)
+        exons = seq.get_features_matching("exon", extend_query=True)
         self.assertEqual(len(exons), 3)
         self.assertEqual(
             str(exons),
             '[exon "fred" at [3:8]/9, exon "fred" at [3:7]/5, exon "fred" at [3:6]/4]',
         )
         # tests extend_query=False, and only get back the base exon
-        exons = seq.get_annotations_matching("exon")
+        exons = seq.get_features_matching("exon")
         self.assertEqual(len(exons), 1)
         self.assertEqual(str(exons), '[exon "fred" at [3:8]/9]')
 
-    def test_get_annotations_matching2(self):
-        """get_annotations_matching returns empty feature if no matches"""
+    def test_get_features_matching2(self):
+        """get_features_matching returns empty feature if no matches"""
 
         # If the sequence does not have a matching feature
         # you get back an empty list, and slicing the sequence
         # with that returns a sequence of length 0.
 
-        dont_exist = self.s.get_annotations_matching("dont_exist")
+        dont_exist = self.s.get_features_matching("dont_exist")
         self.assertEqual(dont_exist, [])
         self.assertEqual(str(self.s[dont_exist]), "")
 
@@ -88,7 +88,7 @@ class FeaturesTest(TestCase):
         # To construct a pseudo-feature covering (or excluding)
         # multiple features, use get_region_covering_all:
 
-        exons = self.s.get_annotations_matching("exon")
+        exons = self.s.get_features_matching("exon")
         self.assertEqual(
             str(self.s.get_region_covering_all(exons)),
             'region "exon" at [10:15, 30:40]/48',
@@ -114,7 +114,7 @@ class FeaturesTest(TestCase):
         # Though .get_region_covering_all also guarantees
         # no overlaps within the result, slicing does not:
 
-        exons = self.s.get_annotations_matching("exon")
+        exons = self.s.get_features_matching("exon")
         self.assertEqual(
             str(self.s.get_region_covering_all(exons + exons)),
             'region "exon" at [10:15, 30:40]/48',
@@ -166,7 +166,7 @@ class FeaturesTest(TestCase):
         # but in those cases where a free-floating feature is created it can
         # ater be attached:
 
-        exons = self.s.get_annotations_matching("exon")
+        exons = self.s.get_features_matching("exon")
         self.assertEqual(len(self.s.annotations), 2)
         region = self.s.get_region_covering_all(exons)
         self.assertEqual(len(self.s.annotations), 2)
@@ -192,7 +192,7 @@ class FeaturesTest(TestCase):
         )
         self.assertEqual(str(plus_cds.get_slice()), "GGGGCCCCCTTTTTTTTTT")
         minus = plus.rc()
-        minus_cds = minus.get_annotations_matching("CDS")[0]
+        minus_cds = minus.get_features_matching("CDS")[0]
         self.assertEqual(str(minus_cds.get_slice()), "GGGGCCCCCTTTTTTTTTT")
 
     def test_feature_from_alignment(self):
@@ -586,15 +586,15 @@ class FeaturesTest(TestCase):
         seq.add_annotation(Feature, "exon", "myname", [(0, 5)])
         got = seq.to_json()
         new = deserialise_object(got)
-        feat = new.get_annotations_matching("exon")[0]
+        feat = new.get_features_matching("exon")[0]
         self.assertEqual(str(feat.get_slice()), "AAAAA")
 
         # now with a list span
         seq = seq[3:]
-        feat = seq.get_annotations_matching("exon")[0]
+        feat = seq.get_features_matching("exon")[0]
         got = seq.to_json()
         new = deserialise_object(got)
-        feat = new.get_annotations_matching("exon")[0]
+        feat = new.get_features_matching("exon")[0]
         self.assertEqual(str(feat.get_slice(complete=False)), "AA")
 
     def test_roundtripped_alignment(self):
@@ -620,7 +620,7 @@ class FeaturesTest(TestCase):
         expect = f.get_slice().to_dict()
         json = aln.to_json()
         new = deserialise_object(json)
-        got = list(new.get_annotations_matching("generic"))[0]
+        got = list(new.get_features_matching("generic"))[0]
         self.assertEqual(got.get_slice().to_dict(), expect)
 
         # annotations on both alignment and sequence
@@ -638,7 +638,7 @@ class FeaturesTest(TestCase):
         self.assertEqual(got_exons.get_slice().to_dict(), expect)
         ## get back the generic
         expect = f.get_slice().to_dict()
-        got = list(new.get_annotations_matching("generic"))[0]
+        got = list(new.get_features_matching("generic"))[0]
         self.assertEqual(got.get_slice().to_dict(), expect)
 
         # check masking of seq features still works
@@ -656,11 +656,10 @@ class FeaturesTest(TestCase):
         # at the alignment level
         sub_aln = aln[:-3]
         s = sub_aln.named_seqs["x"]
-        s.data.get_annotations_matching("exon", "E2")[0]
         d = s.data[:11]
         json = s.to_json()
         new = deserialise_object(json)
-        gf1, gf2 = list(new.data.get_annotations_matching("exon"))
+        gf1, gf2 = list(new.data.get_features_matching("exon"))
         self.assertEqual(str(gf1.get_slice()), "GGGGG")
         self.assertEqual(str(gf2.get_slice()), "C")
         # the sliced alignment
@@ -698,7 +697,7 @@ class FeaturesTest(TestCase):
         y_valued = seq.add_annotation(Variable, "SNP", "freq", xx_y)
         json = seq.to_json()
         new = deserialise_object(json)
-        got = list(new.get_annotations_matching("SNP"))[0]
+        got = list(new.get_features_matching("SNP"))[0]
         # annoyingly, comes back as list of lists
         self.assertEqual(got.xxy_list, [[list(xx), y] for xx, y in y_valued.xxy_list])
 
diff --git a/tests/test_core/test_genetic_code.py b/tests/test_core/test_genetic_code.py
index dfc5371f6..fe578e094 100644
--- a/tests/test_core/test_genetic_code.py
+++ b/tests/test_core/test_genetic_code.py
@@ -3,6 +3,8 @@
 """
 from unittest import TestCase, main
 
+import pytest
+
 from cogent3 import DNA, RNA
 from cogent3.core.genetic_code import (
     DEFAULT,
@@ -414,6 +416,15 @@ class GeneticCodeTests(TestCase):
         self.assertIn("Standard", got)
 
 
-# Run tests if called from command line
-if __name__ == "__main__":
-    main()
+# starting pytest versions
+@pytest.mark.parametrize("code", (1, 2))
+def test_get_alphabet(code):
+    from cogent3.core.alphabet import Alphabet
+
+    gc = get_code()
+    alpha_w_stop = gc.get_alphabet(include_stop=True)
+    assert tuple(alpha_w_stop) == tuple(gc.codons)
+    assert isinstance(alpha_w_stop, Alphabet)
+    alpha_wo_stop = gc.get_alphabet(include_stop=False)
+    assert tuple(alpha_wo_stop) == tuple(gc.sense_codons)
+    assert isinstance(alpha_wo_stop, Alphabet)
diff --git a/tests/test_core/test_moltype.py b/tests/test_core/test_moltype.py
index 52f734ab7..ae3d70cda 100644
--- a/tests/test_core/test_moltype.py
+++ b/tests/test_core/test_moltype.py
@@ -8,7 +8,6 @@ from cogent3.core.moltype import (
     DNA,
     PROTEIN,
     RNA,
-    STANDARD_CODON,
     AlphabetError,
     AlphabetGroup,
     CodonAlphabet,
@@ -975,6 +974,8 @@ class DinucAlphabet(_AlphabetTestCase):
 
 class TestCodonAlphabet(_AlphabetTestCase):
     def setUp(self):
+        from cogent3.core.moltype import STANDARD_CODON
+
         self.alpha = STANDARD_CODON
 
     def test_ambiguous_gaps(self):
diff --git a/tests/test_core/test_sequence.py b/tests/test_core/test_sequence.py
index 9eb6041f7..abed8aa69 100644
--- a/tests/test_core/test_sequence.py
+++ b/tests/test_core/test_sequence.py
@@ -93,15 +93,11 @@ class SequenceTests(TestCase):
     def test_copy(self):
         """correctly returns a copy version of self"""
         s = Sequence("TTTTTTTTTTAAAA", name="test_copy")
-        annot1 = s.add_annotation(Feature, "exon", "annot1", [(0, 10)])
-        annot2 = s.add_annotation(Feature, "exon", "annot2", [(10, 14)])
+        annot1 = s.add_feature("exon", "annot1", [(0, 10)])
+        annot2 = s.add_feature("exon", "annot2", [(10, 14)])
         got = s.copy()
-        got_annot1 = got.get_annotations_matching(
-            annotation_type="exon", name="annot1"
-        )[0]
-        got_annot2 = got.get_annotations_matching(
-            annotation_type="exon", name="annot2"
-        )[0]
+        got_annot1 = got.get_features_matching(feature_type="exon", name="annot1")[0]
+        got_annot2 = got.get_features_matching(feature_type="exon", name="annot2")[0]
         self.assertIsNot(got, s)
         self.assertIsNot(got_annot1, annot1)
         self.assertIsNot(got_annot2, annot2)
@@ -128,23 +124,15 @@ class SequenceTests(TestCase):
         self.assertEqual(r[-1], "G")
         self.assertEqual(r[1:3], "CA")
 
-    def test_conversion(self):
-        """Should convert t to u automatically"""
-        r = self.RNA("TCAtu")
-        self.assertEqual(str(r), "UCAUU")
-
-        d = self.DNA("UCAtu")
-        self.assertEqual(str(d), "TCATT")
-
     def test_to_dna(self):
         """Returns copy of self as DNA."""
-        r = self.RNA("TCA")
+        r = self.RNA("UCA")
         self.assertEqual(str(r), "UCA")
         self.assertEqual(str(r.to_dna()), "TCA")
 
     def test_to_rna(self):
         """Returns copy of self as RNA."""
-        r = self.DNA("UCA")
+        r = self.DNA("TCA")
         self.assertEqual(str(r), "TCA")
         self.assertEqual(str(r.to_rna()), "UCA")
 
@@ -247,10 +235,13 @@ class SequenceTests(TestCase):
 
         sequence = Sequence(seq)
         sequence.annotate_from_gff(gff3_path)
-        matches = [m for m in sequence.get_annotations_matching("*", extend_query=True)]
+        matches = [m for m in sequence.get_features_matching()]
         # 13 features with one having 2 parents, so 14 instances should be found
         self.assertEqual(len(matches), 14)
 
+    @pytest.mark.xfail(
+        reason="todo: annotate_from_gff needs to be using bound annotation_db"
+    )
     def test_annotate_gff_nested_features(self):
         """correctly annotate a sequence with nested features"""
         # the synthetic example
@@ -272,21 +263,21 @@ class SequenceTests(TestCase):
 
         # get the gene and check it has a single annotation and that
         # its slice is correct
-        ann = seq.get_annotations_matching("gene")
+        ann = seq.get_features_matching("gene")
         self.assertEqual(len(ann), 1)
         self.assertEqual(len(ann[0].annotations), 1)
         seq = ann[0].get_slice()
         self.assertEqual(str(seq), "GGAAAATTTTTTTTTAAGGGGGAAAAAAAAA")
 
         # the gene has 1 transcript
-        ann = seq.get_annotations_matching("mRNA", extend_query=True)
+        ann = seq.get_features_matching("mRNA", extend_query=True)
         self.assertEqual(len(ann), 1)
         self.assertEqual(len(ann[0].annotations), 2)  # 2 exons
         seq = ann[0].get_slice()
         self.assertEqual(str(seq), "AAAATTTTTTTTTAAGGGGGAAA")
 
         # the transcript has 2 exons
-        ann = seq.get_annotations_matching("exon", extend_query=True)
+        ann = seq.get_features_matching("exon", extend_query=True)
         self.assertEqual(len(ann), 2)
         exon_seqs = ("TTTTTTTTT", "GGGGG")
         for x in ann:
diff --git a/tests/test_parse/test_genbank.py b/tests/test_parse/test_genbank.py
index 2782ec1e9..74ce74196 100644
--- a/tests/test_parse/test_genbank.py
+++ b/tests/test_parse/test_genbank.py
@@ -446,7 +446,7 @@ ORIGIN
         parser = RichGenbankParser(infile, add_annotation=add_annotation)
 
         seq = [s for l, s in parser][0]
-        cds = dict([(f.name, f) for f in seq.get_annotations_matching("CDS")])
+        cds = dict([(f.name, f) for f in seq.get_features_matching("CDS")])
         expects = {
             "CNA00110": "MAGYDARYGNPLDPMSGGRPSPPETSQQDAYEYSKHGSSSGYLGQLPLGAD"
             "SAQAETASALRTLFGEGADVQALQEPPNQINTLAEGAAVAETGGVLGGDTTRSDNEALAIDPSL"
@@ -478,7 +478,7 @@ ORIGIN
             got_1 = [s for _, s in parser][0]
 
         # name formed from /product value
-        got = {f.name for f in got_1.get_annotations_matching("mRNA")}
+        got = {f.name for f in got_1.get_features_matching("mRNA")}
         self.assertEqual(got, {"conserved hypothetical protein", "chaperone, putative"})
 
         # the file defines itself as DNA
@@ -492,7 +492,7 @@ ORIGIN
 
             self.assertEqual(len(got_1.annotations), len(got_2.annotations))
             self.assertEqual(got_2.moltype.label, moltype)
-            got = {f.name for f in got_1.get_annotations_matching("mRNA")}
+            got = {f.name for f in got_1.get_features_matching("mRNA")}
             self.assertEqual(
                 got, {"conserved hypothetical protein", "chaperone, putative"}
             )
diff --git a/tests/test_parse/test_gff.py b/tests/test_parse/test_gff.py
index 6ac1dc6b7..dcced5ff3 100644
--- a/tests/test_parse/test_gff.py
+++ b/tests/test_parse/test_gff.py
@@ -1,12 +1,14 @@
-#!/usr/bin/env python
 """Unit tests for GFF and related parsers.
 """
 import os
+import re
 
 from io import StringIO
 from pathlib import Path
 from unittest import TestCase, main
 
+import pytest
+
 from cogent3.parse.gff import gff_parser, parse_attributes_gff2
 
 
@@ -19,6 +21,8 @@ __maintainer__ = "Matthew Wakefield"
 __email__ = "wakefield@wehi.edu.au"
 __status__ = "Production"
 
+DATA_DIR = Path(__file__).parent.parent / "data"
+
 headers = [
     """##gff-version 2 
 ##source-version <source> <version text> 
@@ -155,5 +159,50 @@ class GffTest(TestCase):
             self.assertIsInstance(result["Attributes"], str)
 
 
-if __name__ == "__main__":
-    main()
+@pytest.fixture
+def multi_seqid_path(tmp_path):
+    gff_path = DATA_DIR / "ensembl_sample.gff3"
+    data = gff_path.read_text().splitlines()
+    # add two new seqid's using the last 4 lines, twice
+    change = data[-4:]
+    data.extend(_modify_lines(change, "3"))
+    data.extend(_modify_lines(change, "4"))
+    outpath = tmp_path / "ensembl-edited.gff3"
+    outpath.write_text("\n".join(data))
+    return outpath
+
+
+def _modify_lines(change, seqid: str) -> list:
+    ident = re.compile(r"ENS[GT]\d+")
+    result = []
+    for line in change:
+        # change all identifiers so each seqid has unique identifiers
+        for match in ident.findall(line):
+            line = line.replace(match, f"{match}{seqid}")
+        line = line.split("\t")
+        line[0] = seqid
+        line = "\t".join(line)
+        result.append(line)
+    return result
+
+
+@pytest.mark.parametrize("seqids", ("22", ("3",), ("3", "4")))
+def test_seq_names(multi_seqid_path, seqids):
+    got = {r["SeqID"] for r in gff_parser(multi_seqid_path, seqids=seqids)}
+    expect = {seqids} if isinstance(seqids, str) else set(seqids)
+    assert got == expect
+
+
+def test_no_seq_names(multi_seqid_path):
+    got = {r["SeqID"] for r in gff_parser(multi_seqid_path, seqids=None)}
+    expect = {"22", "3", "4"}
+    assert got == expect
+
+
+def test_parse_field_spaces():
+    path = DATA_DIR / "simple.gff"
+    got = list(gff_parser(path))
+    for record in got:
+        for value in record.values():
+            if isinstance(value, str):
+                assert value.strip() == value, "should not have spaces!"