Repository: cogent3/cogent3
Pull request #1409 (2023-05-03 03:36:58)

Annotation refactor


Patch:
diff --git a/src/cogent3/__init__.py b/src/cogent3/__init__.py
index 3de4e33b5..56e2acd9c 100644
--- a/src/cogent3/__init__.py
+++ b/src/cogent3/__init__.py
@@ -19,7 +19,6 @@ from cogent3.core.alignment import (
 )
 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 (
diff --git a/src/cogent3/core/alignment.py b/src/cogent3/core/alignment.py
index 112cdc651..defc74a95 100644
--- a/src/cogent3/core/alignment.py
+++ b/src/cogent3/core/alignment.py
@@ -65,7 +65,6 @@ from cogent3.core.genetic_code import get_code
 from cogent3.core.info import Info as InfoClass
 from cogent3.core.profile import PSSM, MotifCountsArray
 from cogent3.core.sequence import ArraySequence, Sequence, frac_same
-
 # which is a circular import otherwise.
 from cogent3.format.alignment import save_to_filename
 from cogent3.format.fasta import alignment_to_fasta
@@ -389,7 +388,7 @@ def merged_db_collection(seqs) -> SupportsFeatures:
                 f"annotation db types must be equal: {type(first)} != {type(db)}"
             )
 
-        if first is None or db is None:
+        if first is None or db is None or first is db:
             continue
         first.update(db)
 
@@ -597,6 +596,8 @@ class _SequenceCollectionBase:
         if value and not isinstance(value, SupportsFeatures):
             raise TypeError
         self._annotation_db = value
+        for seq in self.seqs:
+            seq.annotation_db = value
 
     def __str__(self):
         """Returns self in FASTA-format, respecting name order."""
@@ -2229,7 +2230,7 @@ class SequenceCollection(_SequenceCollectionBase):
         seqid: Optional[str] = None,
         biotype: Optional[str] = None,
         name: Optional[str] = None,
-        on_alignment: Optional[bool] = None,
+        allow_partial: bool = False,
     ) -> Iterator[Annotation]:
         """yields Annotation instances
 
@@ -2241,9 +2242,8 @@ class SequenceCollection(_SequenceCollectionBase):
             biotype of the feature, e.g. CDS, gene
         name
             name of the feature
-        on_alignment
-            limit query to features on Alignment, ignores sequences. Ignored on
-            SequenceCollection instances.
+        allow_partial
+            allow features partially overlaping self
 
         Notes
         -----
@@ -2251,8 +2251,10 @@ class SequenceCollection(_SequenceCollectionBase):
         yield a sequence segment that is consistently oriented irrespective
         of strand of the current instance.
         """
-        on_alignment = None if not isinstance(self, Alignment) else on_alignment
         if self.annotation_db is None:
+            # todo gah this should not be required
+            # the annotation db attribute should be defined on creation,
+            # or when annotation_db is assigned
             anno_db = merged_db_collection(self.seqs)
             self.annotation_db = anno_db
 
@@ -2260,41 +2262,27 @@ class SequenceCollection(_SequenceCollectionBase):
             return None
 
         seqids = [seqid] if isinstance(seqid, str) else seqid
-
-        if seqid and not set(seqid) & set(self.names):
-            raise ValueError(f"unknown {seqid=}")
-
-        if seqids and on_alignment:
-            raise ValueError(f"cannot provide seqids with on_alignment")
-
         if seqids is None:
-            seq_map = None
-            for feature in self.annotation_db.get_features_matching(
-                biotype=biotype, name=name, on_alignment=on_alignment
-            ):
-                if on_al := feature.pop("on_alignment", on_alignment):
-                    if seq_map is None:
-                        seq_map = self.seqs[0].map
-                    # gah todo use the Aligned map to transform the absolute
-                    # index from the db into a relative index
-                    spans = numpy.array(feature["spans"])
-                    spans = seq_map.relative_position(spans)
-                    feature["spans"] = spans.tolist()
-                    # and if i've been reversed...?
-                    feature["reversed"] = seq_map.reverse
-                    yield self.make_feature(feature=feature, on_alignment=on_al)
-                    continue
+            seqids = self.names
 
-                seqid = feature["seqid"]
-                seq = self.named_seqs[seqid]
-                # passing self only used when self is an Alignment
-                yield seq.make_feature(feature, self)
-            return
+        if seqid and not set(seqids) & set(self.names):
+            raise ValueError(f"unknown {seqid=}")
 
         for seqid in seqids:
             seq = self.named_seqs[seqid]
+            if isinstance(seq, Aligned):
+                start, end = seq.map.start, seq.map.end
+            else:
+                start, end = 0, len(seq)
+
             for feature in self.annotation_db.get_features_matching(
-                seqid=seqid, biotype=biotype, name=name, on_alignment=None
+                seqid=seqid,
+                biotype=biotype,
+                name=name,
+                on_alignment=False,
+                allow_partial=allow_partial,
+                start=start,
+                end=end,
             ):
                 # passing self only used when self is an Alignment
                 yield seq.make_feature(feature, self)
@@ -4747,6 +4735,20 @@ class Alignment(_Annotatable, AlignmentI, SequenceCollection):
 
         return f"{len(self.names)} x {self.seq_len} alignment: {seqs}"
 
+    @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 value and not isinstance(value, SupportsFeatures):
+            raise TypeError
+        self._annotation_db = value
+        for seq in self.seqs:
+            seq.data.annotation_db = value
+
     def _mapped(self, slicemap):
         align = []
         for name in self.names:
@@ -4769,9 +4771,9 @@ class Alignment(_Annotatable, AlignmentI, SequenceCollection):
             raise ValueError("Annotation does not belong to this alignment")
         return annot.remapped_to(target_aligned.data, target_aligned.map)
 
-    def get_projected_annotations(self, seq_name, *args):
-        aln_annots = self.get_features_matching(*args)
-        return [self.project_annotation(seq_name, a) for a in aln_annots]
+    def get_projected_annotations(self, *, seqid: str, **kwargs):
+        aln_annots = list(self.get_features(**kwargs))
+        return [self.project_annotation(seqid, a) for a in aln_annots]
 
     @c3warn.deprecated_callable(
         "2023.7",
@@ -5213,3 +5215,77 @@ class Alignment(_Annotatable, AlignmentI, SequenceCollection):
             fmap = fmap.nucleic_reversed()
         feature.pop("strand", None)
         return Annotation(parent=self, map=fmap, **feature)
+
+    def get_features(
+        self,
+        *,
+        seqid: Optional[str] = None,
+        biotype: Optional[str] = None,
+        name: Optional[str] = None,
+        on_alignment: Optional[bool] = None,
+        allow_partial: bool = False,
+    ) -> Iterator[Annotation]:
+        """yields Annotation instances
+
+        Parameters
+        ----------
+        seqid
+            limit search to features on this named sequence, defaults to search all
+        biotype
+            biotype of the feature, e.g. CDS, gene
+        name
+            name of the feature
+        on_alignment
+            limit query to features on Alignment, ignores sequences. Ignored on
+            SequenceCollection instances.
+        allow_partial
+            allow features partially overlaping self
+
+        Notes
+        -----
+        When dealing with a nucleic acid moltype, the returned features will
+        yield a sequence segment that is consistently oriented irrespective
+        of strand of the current instance.
+        """
+        # we only do on-alignment in here
+        if not on_alignment:
+            kwargs = {
+                k: v for k, v in locals().items() if k not in ("self", "__class__")
+            }
+            kwargs.pop("on_alignment")
+            yield from super().get_features(**kwargs)
+
+        if on_alignment == False:
+            return
+
+        if self.annotation_db is None:
+            anno_db = merged_db_collection(self.seqs)
+            self.annotation_db = anno_db
+
+        if self.annotation_db is None:
+            return None
+
+        seq_map = None
+        for feature in self.annotation_db.get_features_matching(
+            biotype=biotype,
+            name=name,
+            on_alignment=on_alignment,
+            allow_partial=allow_partial,
+        ):
+            if on_al := feature.pop("on_alignment", on_alignment):
+                if seq_map is None:
+                    seq_map = self.seqs[0].map
+                # gah todo use the Aligned map to transform the absolute
+                # index from the db into a relative index
+                spans = numpy.array(feature["spans"])
+                spans = seq_map.relative_position(spans)
+                feature["spans"] = spans.tolist()
+                # and if i've been reversed...?
+                feature["reversed"] = seq_map.reverse
+                yield self.make_feature(feature=feature, on_alignment=on_al)
+                continue
+
+                seqid = feature["seqid"]
+                seq = self.named_seqs[seqid]
+                # passing self only used when self is an Alignment
+                yield seq.make_feature(feature, self)
diff --git a/src/cogent3/core/annotation.py b/src/cogent3/core/annotation.py
index 42d3a5298..cd6d91f40 100644
--- a/src/cogent3/core/annotation.py
+++ b/src/cogent3/core/annotation.py
@@ -447,14 +447,34 @@ class Annotation(_AnnotationCore, _Serialisable):
         }
         return self.__class__(**kwargs)
 
-    def get_slice(self, complete=True):
-        """The corresponding sequence fragment.  If 'complete' is true
-        and the full length of this feature is not present in the sequence
-        then this method will fail."""
+    def get_slice(self, complete: bool = True, allow_gaps: bool = False):
+        """
+        The corresponding sequence fragment.
+
+        Parameters
+        ----------
+        complete
+            if feature not complete on parent,causes an exception to be
+            raised. If False, gaps are removed.
+        allow_gaps
+            if on an alignment, includes the gap positions
+
+        Returns
+        -------
+        a slice of self.parent
+
+        Notes
+        -----
+        If 'complete' is true and the full length of this feature is not
+        present in the sequence, then this method will fail.
+        """
+        # todo gah set allow_gaps=True as the default
         map = self.map
         if not (complete or map.complete):
             map = map.without_gaps()
-        return self.parent[map]
+        if not allow_gaps:
+            return self.parent[map]
+        return self.parent[map.start : map.end]
 
     def without_lost_spans(self):
         """Keeps only the parts which are actually present in the underlying sequence"""
@@ -509,11 +529,6 @@ class Annotation(_AnnotationCore, _Serialisable):
         name = f' "{self.name}"'
         return f'"{self.seqid}" {self.biotype}{name} at {self.map}'
 
-    def _projected_to_base(self, base):
-        if self.parent == base:
-            return self.__class__(base, self.map)
-        return self.remapped_to(base, self.parent._projected_to_base(base).map)
-
     def remapped_to(self, grandparent, gmap):
         kwargs = {
             **self._serialisable,
diff --git a/src/cogent3/core/annotation_db.py b/src/cogent3/core/annotation_db.py
index 087035101..158adeb3b 100644
--- a/src/cogent3/core/annotation_db.py
+++ b/src/cogent3/core/annotation_db.py
@@ -230,7 +230,7 @@ def _matching_conditions(
             # todo gah FIX this excludes a False value from being a condition!
             # conditions are filtered for None before here, so we should add
             # an else where the op is assigned !=
-            if val:
+            if val is not None:
                 op = "LIKE" if isinstance(val, str) and "%" in val else "="
                 conds.append(f"{col} {op} ?")
                 vals.append(val)
@@ -506,7 +506,7 @@ class SqliteAnnotationDbMixin:
         table_names = ["user"] if on_alignment else self.table_names
         for table_name in table_names:
             for result in self._get_records_matching(table_name, **kwargs):
-                yield {k: result[k] for k in result.keys()}
+                yield dict(zip(result.keys(), result))
 
     def get_features_matching(
         self,
@@ -527,23 +527,21 @@ class SqliteAnnotationDbMixin:
         table_names = ["user"] if on_alignment else self.table_names
         for table_name in table_names:
             columns = ("seqid", "biotype", "spans", "strand", "name")
+            query_args = {**kwargs}
+
             if table_name == "user":
                 columns += ("on_alignment",)
+            else:
+                query_args.pop("on_alignment", None)
+
             for result in self._get_records_matching(
-                table_name=table_name, columns=columns, **kwargs
+                table_name=table_name, columns=columns, **query_args
             ):
-                if "on_alignment" in result.keys():
-                    on_alignment = result["on_alignment"]
-                else:
-                    on_alignment = None
-                yield FeatureDataType(
-                    seqid=result["seqid"],
-                    biotype=result["biotype"],
-                    name=result["name"],
-                    spans=[tuple(c) for c in result["spans"]],
-                    reversed=result["strand"] == "-",
-                    on_alignment=on_alignment,
-                )
+                result = dict(zip(result.keys(), result))
+                result["on_alignment"] = result.get("on_alignment", None)
+                result["spans"] = [tuple(c) for c in result["spans"]]
+                result["reversed"] = result.pop("strand", None) == "-"
+                yield result
 
     def num_matches(
         self,
diff --git a/src/cogent3/core/sequence.py b/src/cogent3/core/sequence.py
index 541449e17..4e4970e54 100644
--- a/src/cogent3/core/sequence.py
+++ b/src/cogent3/core/sequence.py
@@ -890,6 +890,7 @@ class Sequence(_Annotatable, SequenceI):
     )
     def get_features(
         self,
+        *,
         biotype: Optional[str] = None,
         name: Optional[str] = None,
         start: Optional[int] = None,
@@ -1027,7 +1028,6 @@ class Sequence(_Annotatable, SequenceI):
         """
         feature = dict(feature)
         seq_rced = self._seq.reversed
-        feature.pop("seqid", None)  # not required here as provided
         # todo gah check consistency of relationship between reversed and strand
         # i.e. which object has responsibility for transforming the strand value
         # (a string) into a bool?
@@ -1081,6 +1081,7 @@ class Sequence(_Annotatable, SequenceI):
             fmap = fmap.nucleic_reversed().reversed()
 
         feature.pop("on_alignment", None)
+        feature.pop("seqid", None)
         return Annotation(parent=self, seqid=self.name, map=fmap, **feature)
 
     def annotate_from_gff(self, f: os.PathLike, offset=None):
@@ -1233,13 +1234,23 @@ class Sequence(_Annotatable, SequenceI):
         annotations = []
         annot_types = [annot_types, [annot_types]][isinstance(annot_types, str)]
         for annot_type in annot_types:
-            annotations += self.get_annotations_matching(
-                annot_type, extend_query=extend_query
+            annotations += list(
+                self.get_features(biotype=annot_type, allow_partial=True)
+            )
+
+        if not annotations:
+            region = Annotation(
+                parent=self,
+                seqid=self.name,
+                name=None,
+                biotype=None,
+                map=Map(locations=[], parent_length=len(self)),
             )
+        else:
+            region = annotations[0].union(annotations[1:])
 
-        region = self.get_region_covering_all(annotations, extend_query=extend_query)
         if shadow:
-            region = region.get_shadow()
+            region = region.shadow()
 
         i = 0
         segments = []


Test patch:
diff --git a/tests/test_core/test_annotation.py b/tests/test_core/test_annotation.py
index 6cc26fba3..6624d9864 100644
--- a/tests/test_core/test_annotation.py
+++ b/tests/test_core/test_annotation.py
@@ -105,60 +105,6 @@ class TestAnnotations(unittest.TestCase):
         with self.assertRaises(AssertionError):
             _ = annot.copy_annotations_to(seq[:-2])
 
-    @pytest.mark.xfail(reason="todo gah change to new API")
-    def test_reverse_complement(self):
-        """test correct translation of annotations on reverse complement."""
-        aln_expecteds = {
-            "misc_feature": {"FAKE01": "TTTGGGGGGGGGG", "FAKE02": "TTTGGGGGGGGGG"},
-            "CDS": {"FAKE01": "GGGGGGGGGG", "FAKE02": "GGGGGGGGGG"},
-            "5'UTR": {"FAKE01": "CC", "FAKE02": "CC"},
-            "LTR": {"FAKE01": "CCCAAAATTTTTT", "FAKE02": "CCC-----TTTTT"},
-        }
-
-        seq_expecteds = {
-            "CDS": {"FAKE01": "GGGGGGGGGG", "FAKE02": "GGGGGGGGGG"},
-            "5'UTR": {"FAKE01": "TTT", "FAKE02": "TTT"},
-        }
-
-        rc = self.aln.rc()
-        # rc'ing an Alignment or Sequence rc's their annotations too. This means
-        # slicing returns the same sequence as the non-rc'd alignment/seq
-        for annot_type in ["misc_feature", "CDS", "5'UTR", "LTR"]:
-            observed = list(self.aln.get_by_annotation(annot_type))[0].to_dict()
-            expected = aln_expecteds[annot_type]
-            assert observed == expected, ("+", annot_type, expected, observed)
-            observed = list(rc.get_by_annotation(annot_type))[0].to_dict()
-            expected = aln_expecteds[annot_type]
-            assert observed == expected, ("-", annot_type, expected, observed)
-
-            if annot_type in ["misc_feature", "LTR"]:
-                continue  # because seqs haven't been annotated with it
-            for name in self.aln.names:
-                observed = list(
-                    self.aln.named_seqs[name].data.get_by_annotation(annot_type)
-                )[0]
-                observed = str(observed)
-                expected = seq_expecteds[annot_type][name]
-                assert str(observed) == expected, (
-                    "+",
-                    annot_type,
-                    name,
-                    expected,
-                    observed,
-                )
-                observed = list(rc.named_seqs[name].data.get_by_annotation(annot_type))[
-                    0
-                ]
-                observed = str(observed)
-                expected = seq_expecteds[annot_type][name]
-                assert str(observed) == expected, (
-                    "-",
-                    annot_type,
-                    name,
-                    expected,
-                    observed,
-                )
-
 
 class TestMapSpans(unittest.TestCase):
     """Test attributes of Map & Spans classes critical to annotation
diff --git a/tests/test_core/test_annotation_db.py b/tests/test_core/test_annotation_db.py
index be52b9ac2..51285cc53 100644
--- a/tests/test_core/test_annotation_db.py
+++ b/tests/test_core/test_annotation_db.py
@@ -501,34 +501,35 @@ def test_annotate_from_gff_multiple_calls(seq):
     assert len(list(seq.get_features_matching())) == 10
 
 
-@pytest.mark.xfail(reason="todo gah implement lost spans")
 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")
+    seq_coll = SequenceCollection(seqs)
+    seq_coll.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
-    got = list(seq_collection.get_seq("test_seq").get_features(allow_partial=True))
-    assert len(got) == 5
-
-    # the seq for which the seqid was NOT provided is NOT annotated
-    assert seq_collection.get_seq("test_seq2").annotation_db is None
+    seq = seq_coll.get_seq("test_seq")
+    assert seq_coll.get_seq("test_seq").annotation_db is not 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
-    )
+    assert seq_coll.get_seq("test_seq").annotation_db is seq_coll.annotation_db
+    got = list(seq_coll.get_seq("test_seq").get_features(allow_partial=True))
+    assert len(got) == 5
     got = list(seq.get_features(feature_type="CDS"))
     assert len(got) == 2
 
     got = list(seq.get_features(feature_type="CpG"))
     assert len(got) == 1
 
+    # the seq for which the seqid was NOT provided also has a reference to the same db
+    seq2 = seq_coll.get_seq("test_seq2")
+    assert seq2.annotation_db is not None
+    # querying on that sequence returns []
+    got = list(seq2.get_features(feature_type="CDS"))
+    assert not got
+
 
 def test_seq_coll_query():
     """obtain same results when querying from collection as from seq"""
@@ -539,10 +540,17 @@ def test_seq_coll_query():
     seq = seq_coll.get_seq("test_seq")
     # the seq for which the seqid was provided is annotated
     assert seq.annotation_db is not None
-    expect = list(seq.get_features_matching())
-    got = seq_coll.get_features(seqid="test_seq")
-    # the seq for which the seqid was NOT provided is NOT annotated
-    assert seq_coll.get_seq("test_seq2").annotation_db is None
+    # todo gah this test fails when allow_partial=False because start / stop
+    # not used by seqcoll method
+    expect = set(
+        (f.seqid, f.biotype, f.name, str(f.map))
+        for f in seq.get_features(allow_partial=True)
+    )
+    got = set(
+        (f.seqid, f.biotype, f.name, str(f.map))
+        for f in seq_coll.get_features(seqid="test_seq", allow_partial=True)
+    )
+    assert got == expect
     # the annotation_db on the seq and the seq collection are the same object
     assert seq.annotation_db is seq_coll.annotation_db
     got = list(seq.get_features_matching(feature_type="CDS"))
diff --git a/tests/test_core/test_features.py b/tests/test_core/test_features.py
index f92f8ebe3..2481004d8 100755
--- a/tests/test_core/test_features.py
+++ b/tests/test_core/test_features.py
@@ -5,7 +5,6 @@ import pytest
 from cogent3 import ASCII, DNA, make_aligned_seqs
 from cogent3.core.annotation import Feature, Variable
 from cogent3.core.annotation_db import GffAnnotationDb
-
 # Complete version of manipulating sequence annotations
 from cogent3.util.deserialise import deserialise_object
 
@@ -104,57 +103,6 @@ class FeaturesTest(TestCase):
         minus_cds = minus.get_features(biotype="CDS")[0]
         self.assertEqual(str(minus_cds.get_slice()), "GGGGCCCCCTTTTTTTTTT")
 
-    @pytest.mark.xfail(reason="todo gah update test to use latest API")
-    def test_feature_from_alignment(self):
-        """seq features obtained from the alignment"""
-
-        # Sequence features can be accessed via a containing Alignment:
-
-        aln = make_aligned_seqs(
-            data=[["x", "-AAAAAAAAA"], ["y", "TTTT--TTTT"]], array_align=False
-        )
-        self.assertEqual(str(aln), ">x\n-AAAAAAAAA\n>y\nTTTT--TTTT\n")
-        exon = aln.get_seq("x").add_annotation(Feature, "exon", "fred", [(3, 8)])
-        aln_exons = aln.get_features(seqid="x", biotype="exon")
-        aln_exons = aln.get_features(biotype="exon")
-        # But these will be returned as **alignment**
-        # features with locations in alignment coordinates.
-
-        self.assertEqual(str(exon), 'exon "fred" at [3:8]/9')
-        self.assertEqual(str(aln_exons[0]), 'exon "fred" at [4:9]/10')
-        self.assertEqual(str(aln_exons[0].get_slice()), ">x\nAAAAA\n>y\n--TTT\n")
-        aln_exons[0].attach()
-        self.assertEqual(len(aln.annotations), 1)
-
-        # Similarly alignment features can be projected onto the aligned sequences,
-        # where they may end up falling across gaps:
-
-        exons = aln.get_projected_annotations("y", "exon")
-        self.assertEqual(str(exons), '[exon "fred" at [-2-, 4:7]/8]')
-        self.assertEqual(str(aln.get_seq("y")[exons[0].map.without_gaps()]), "TTT")
-
-        # We copy the annotations from another sequence,
-
-        aln = make_aligned_seqs(
-            data=[["x", "-AAAAAAAAA"], ["y", "TTTT--CCCC"]], array_align=False
-        )
-        self.s = DNA.make_seq("AAAAAAAAA", name="x")
-        exon = self.s.add_annotation(Feature, "exon", "fred", [(3, 8)])
-        exon = aln.get_seq("x").copy_annotations(self.s)
-        aln_exons = list(aln.get_features(seqid="x", biotype="exon"))
-        self.assertEqual(str(aln_exons), '[exon "fred" at [4:9]/10]')
-
-        # even if the name is different.
-
-        exon = aln.get_seq("y").copy_annotations(self.s)
-        aln_exons = list(aln.get_features(seqid="y", biotype="exon"))
-        self.assertEqual(str(aln_exons), '[exon "fred" at [3:4, 6:10]/10]')
-        self.assertEqual(str(aln[aln_exons]), ">x\nAAAAA\n>y\nTCCCC\n")
-
-        # default for get_annotations_from_any_seq is return all features
-        got = aln.get_features()
-        self.assertEqual(len(got), 2)
-
     @pytest.mark.xfail(reason="todo gah implement lost spans")
     def test_lost_spans(self):
         """features no longer included in an alignment represented by lost spans"""
@@ -172,204 +120,6 @@ class FeaturesTest(TestCase):
         self.assertEqual(str(copied), '[exon "A" at [5:5, -4-]/5]')
         self.assertEqual(str(copied[0].get_slice()), ">x\n----\n>y\n----\n")
 
-    @pytest.mark.xfail(reason="todo gah update test to use latest API")
-    def test_seq_shorter(self):
-        """lost spans on shorter sequences"""
-
-        # If the sequence is shorter, again you get a lost span.
-
-        aln = make_aligned_seqs(
-            data=[["x", "-AAAAAAAAA"], ["y", "TTTT--TTTT"]], array_align=False
-        )
-        diff_len_seq = DNA.make_seq("CCCCCCCCCCCCCCCCCCCCCCCCCCCC", "x")
-        diff_len_seq.add_feature("repeat", "A", [(12, 14)])
-        aln.get_seq("y").copy_annotations(diff_len_seq)
-        copied = list(aln.get_features(seqid="y", name="repeat"))
-        self.assertEqual(str(copied), '[repeat "A" at [10:10, -6-]/10]')
-
-    @pytest.mark.xfail(reason="todo gah update test to use latest API")
-    def test_terminal_gaps(self):
-        """features in cases of terminal gaps"""
-
-        # We consider cases where there are terminal gaps.
-
-        aln = make_aligned_seqs(
-            data=[["x", "-AAAAAAAAA"], ["y", "------TTTT"]], array_align=False
-        )
-        feat = dict(seqid="x", biotype="exon", name="fred", spans=[(3, 8)])
-        aln.add_feature(**feat)
-        aln_exons = list(aln.get_features(seqid="x", name="exon"))
-        self.assertEqual(str(aln_exons), '[exon "fred" at [4:9]/10]')
-        self.assertEqual(str(aln_exons[0].get_slice()), ">x\nAAAAA\n>y\n--TTT\n")
-        aln = make_aligned_seqs(
-            data=[["x", "-AAAAAAAAA"], ["y", "TTTT--T---"]], array_align=False
-        )
-        aln.add_feature(**feat)
-        aln_exons = list(aln.get_features(seqid="x", name="exon"))
-        self.assertEqual(str(aln_exons[0].get_slice()), ">x\nAAAAA\n>y\n--T--\n")
-
-    @pytest.mark.xfail(reason="todo gah change test to new api")
-    def test_annotated_region_masks(self):
-        """masking a sequence with specific features"""
-
-        # Annotated regions can be masked (observed sequence characters
-        # replaced by another), either through the sequence on which they
-        # reside or by projection from the alignment. Note that mask_char must
-        # be a valid character for the sequence MolType. Either the features
-        # (multiple can be named), or their shadow, can be masked.
-
-        # We create an alignment with a sequence that has two different annotation types.
-
-        aln = make_aligned_seqs(
-            data=[["x", "C-CCCAAAAAGGGAA"], ["y", "-T----TTTTG-GTT"]], array_align=False
-        )
-        self.assertEqual(str(aln), ">x\nC-CCCAAAAAGGGAA\n>y\n-T----TTTTG-GTT\n")
-        exon = aln.get_seq("x").add_feature(
-            biotype="exon", name="norwegian", spans=[(0, 4)]
-        )
-        self.assertEqual(str(exon.get_slice()), "CCCC")
-        repeat = aln.get_seq("x").add_feature(
-            biotype="repeat", name="blue", spans=[(9, 12)]
-        )
-        self.assertEqual(str(repeat.get_slice()), "GGG")
-        repeat = aln.get_seq("y").add_feature(
-            biotype="repeat", name="frog", spans=[(5, 7)]
-        )
-        self.assertEqual(str(repeat.get_slice()), "GG")
-
-        # Each sequence should correctly mask either the single feature,
-        # it's shadow, or the multiple features, or shadow.
-
-        self.assertEqual(
-            str(aln.get_seq("x").with_masked_annotations("exon", mask_char="?")),
-            "????AAAAAGGGAA",
-        )
-        self.assertEqual(
-            str(
-                aln.get_seq("x").with_masked_annotations(
-                    "exon", mask_char="?", shadow=True
-                )
-            ),
-            "CCCC??????????",
-        )
-        self.assertEqual(
-            str(
-                aln.get_seq("x").with_masked_annotations(
-                    ["exon", "repeat"], mask_char="?"
-                )
-            ),
-            "????AAAAA???AA",
-        )
-        self.assertEqual(
-            str(
-                aln.get_seq("x").with_masked_annotations(
-                    ["exon", "repeat"], mask_char="?", shadow=True
-                )
-            ),
-            "CCCC?????GGG??",
-        )
-        self.assertEqual(
-            str(aln.get_seq("y").with_masked_annotations("exon", mask_char="?")),
-            "TTTTTGGTT",
-        )
-        self.assertEqual(
-            str(aln.get_seq("y").with_masked_annotations("repeat", mask_char="?")),
-            "TTTTT??TT",
-        )
-        self.assertEqual(
-            str(
-                aln.get_seq("y").with_masked_annotations(
-                    "repeat", mask_char="?", shadow=True
-                )
-            ),
-            "?????GG??",
-        )
-
-        # The same methods can be applied to annotated Alignment's.
-
-        self.assertEqual(
-            str(aln.with_masked_annotations("exon", mask_char="?")),
-            ">x\n?-???AAAAAGGGAA\n>y\n-T----TTTTG-GTT\n",
-        )
-        self.assertEqual(
-            str(aln.with_masked_annotations("exon", mask_char="?", shadow=True)),
-            ">x\nC-CCC??????????\n>y\n-?----?????-???\n",
-        )
-        self.assertEqual(
-            str(aln.with_masked_annotations("repeat", mask_char="?")),
-            ">x\nC-CCCAAAAA???AA\n>y\n-T----TTTT?-?TT\n",
-        )
-        self.assertEqual(
-            str(aln.with_masked_annotations("repeat", mask_char="?", shadow=True)),
-            ">x\n?-????????GGG??\n>y\n-?----????G-G??\n",
-        )
-        self.assertEqual(
-            str(aln.with_masked_annotations(["repeat", "exon"], mask_char="?")),
-            ">x\n?-???AAAAA???AA\n>y\n-T----TTTT?-?TT\n",
-        )
-        self.assertEqual(
-            str(aln.with_masked_annotations(["repeat", "exon"], shadow=True)),
-            ">x\nC-CCC?????GGG??\n>y\n-?----????G-G??\n",
-        )
-
-    @pytest.mark.xfail(reason="todo gah implement support for projection")
-    def test_projected_to_base(self):
-        """tests a given annotation is correctly projected on the base sequence"""
-
-        seq = DNA.make_seq("AAAAAAAAATTTTTTTTT", name="x")
-        layer_one = seq.add_feature("repeat", "frog", [(1, 17)])
-        layer_two = layer_one.add_feature("repeat", "frog", [(2, 16)])
-        got = layer_two._projected_to_base(seq)
-        self.assertEqual(got.map.start, 3)
-        self.assertEqual(got.map.end, 17)
-        self.assertEqual(got.map.parent_length, len(seq))
-
-        layer_three = layer_two.add_feature("repeat", "frog", [(5, 10)])
-        got = layer_three._projected_to_base(seq)
-        self.assertEqual(got.map.start, 8)
-        self.assertEqual(got.map.end, 13)
-        self.assertEqual(got.map.parent_length, len(seq))
-
-        layer_four = layer_three.add_feature("repeat", "frog", [(0, 4)])
-        layer_five = layer_four.add_feature("repeat", "frog", [(1, 2)])
-        got = layer_five._projected_to_base(seq)
-        self.assertEqual(got.map.start, 9)
-        self.assertEqual(got.map.end, 10)
-        self.assertEqual(got.map.parent_length, len(seq))
-
-    @pytest.mark.xfail(reason="todo gah update test to use latest API")
-    def test_nested_annotated_region_masks(self):
-        """masking a sequence with specific features when nested annotations"""
-
-        aln = make_aligned_seqs(
-            data=[["x", "C-GGCAAAAATTTAA"], ["y", "-T----TTTTG-GTT"]], array_align=False
-        )
-        gene = aln.get_seq("x").add_feature("gene", "norwegian", [(0, 4)])
-        self.assertEqual(str(gene.get_slice()), "CGGC")
-        gene.add_feature("repeat", "blue", [(1, 3)])
-        # evaluate the sequence directly
-        masked = str(
-            aln.get_seq("x").with_masked_annotations(
-                "repeat", mask_char="?", extend_query=True
-            )
-        )
-        self.assertEqual(masked, "C??CAAAAATTTAA")
-
-        exon = aln.get_seq("y").add_feature("repeat", "frog", [(1, 4)])
-        self.assertEqual(str(exon.get_slice()), "TTT")
-        # evaluate the sequence directly
-        masked = str(
-            aln.get_seq("y").with_masked_annotations(
-                "repeat", mask_char="?", extend_query=True
-            )
-        )
-        self.assertEqual(masked, "T???TGGTT")
-
-        masked = aln.with_masked_annotations("gene", mask_char="?")
-        got = masked.to_dict()
-        self.assertEqual(got["x"], "?-???AAAAATTTAA")
-        self.assertEqual(got["y"], "-T----TTTTG-GTT")
-
     @pytest.mark.xfail(reason="todo gah update test to use latest API")
     def test_annotated_separately_equivalence(self):
         """allow defining features as a series or individually"""
@@ -576,11 +326,6 @@ class FeaturesTest(TestCase):
         # annoyingly, comes back as list of lists
         self.assertEqual(got.xxy_list, [[list(xx), y] for xx, y in y_valued.xxy_list])
 
-    @pytest.mark.xfail(reason="todo gah change test to use feature to get nested")
-    def test_nested_get_slice(self):
-        """check the get_slice method works on nested annotations"""
-        self.assertEqual(self.nested_feature.get_slice(), "CCC")
-
     @pytest.mark.xfail(reason="todo gah change test to use latest API")
     def test_nested_to_rich_dict(self):
         """check the to_rich_dict method works with nested annotations"""
@@ -760,3 +505,211 @@ def test_feature_query_parent_aln():
     parent = parent[0]
     assert parent.name == "GG"
     assert parent.get_slice().to_dict() == aln[0:10].to_dict()
+
+
+def test_aln_feature_lost_spans():
+    """features outside the sequence should not be returned"""
+    db = GffAnnotationDb(data=[])
+    db.add_feature(seqid="y", biotype="repeat", name="A", spans=[(12, 14)])
+    # If the sequence is shorter, again you get a lost span.
+    aln = make_aligned_seqs(
+        data={"x": "-AAAAAAAAA", "y": "TTTT--TTTT"}, array_align=False
+    )
+    aln.annotation_db = db
+    copied = list(aln.get_features(seqid="y", biotype="repeat"))
+    assert not copied
+
+
+def test_terminal_gaps():
+    """features in cases of terminal gaps"""
+
+    # We consider cases where there are terminal gaps.
+    db = GffAnnotationDb()
+    feat = dict(seqid="x", biotype="exon", name="fred", spans=[(3, 8)])
+    db.add_feature(**feat)
+    aln = make_aligned_seqs(
+        data=[["x", "-AAAAAAAAA"], ["y", "------TTTT"]], array_align=False
+    )
+    aln.annotation_db = db
+    aln_exons = list(aln.get_features(seqid="x", biotype="exon"))
+    assert 'exon "fred" at [4:9]/10' in str(aln_exons)
+    assert aln_exons[0].get_slice().to_dict() == dict(x="AAAAA", y="--TTT")
+    aln = make_aligned_seqs(
+        data=[["x", "-AAAAAAAAA"], ["y", "TTTT--T---"]], array_align=False
+    )
+    aln.annotation_db = db
+    aln_exons = list(aln.get_features(seqid="x", biotype="exon"))
+    assert aln_exons[0].get_slice().to_dict() == dict(x="AAAAA", y="--T--")
+
+
+def test_annotated_region_masks():
+    """masking a sequence with specific features"""
+
+    # Annotated regions can be masked (observed sequence characters
+    # replaced by another), either through the sequence on which they
+    # reside or by projection from the alignment. Note that mask_char must
+    # be a valid character for the sequence MolType. Either the features
+    # (multiple can be named), or their shadow, can be masked.
+
+    # We create an alignment with a sequence that has two different annotation types.
+    orig_data = {"x": "C-CCCAAAAAGGGAA", "y": "-T----TTTTG-GTT"}
+    db = GffAnnotationDb()
+    db.add_feature(seqid="x", biotype="exon", name="norwegian", spans=[(0, 4)])
+    db.add_feature(
+        biotype="repeat",
+        name="blue",
+        spans=[(9, 12)],
+        seqid="x",
+    )
+    db.add_feature(seqid="y", biotype="repeat", name="frog", spans=[(5, 7)])
+    aln = make_aligned_seqs(data=orig_data, array_align=False, moltype="dna")
+    aln.annotation_db = db
+
+    assert aln.to_dict() == {"x": "C-CCCAAAAAGGGAA", "y": "-T----TTTTG-GTT"}
+    x = aln.get_seq("x")
+    y = aln.get_seq("y")
+    exon = list(x.get_features(biotype="exon"))[0]
+    assert str(exon.get_slice()) == "CCCC"
+    repeat_x = list(x.get_features(biotype="repeat"))[0]
+    assert str(repeat_x.get_slice()) == "GGG"
+    repeat_y = list(y.get_features(biotype="repeat"))[0]
+    assert str(repeat_y.get_slice()) == "GG"
+
+    # Each sequence should correctly mask either the single feature,
+    # it's shadow, or the multiple features, or shadow.
+
+    assert (
+        str(aln.get_seq("x").with_masked_annotations("exon", mask_char="?"))
+        == "????AAAAAGGGAA"
+    )
+    assert (
+        str(
+            aln.get_seq("x").with_masked_annotations("exon", mask_char="?", shadow=True)
+        )
+        == "CCCC??????????"
+    )
+    assert (
+        str(aln.get_seq("x").with_masked_annotations(["exon", "repeat"], mask_char="?"))
+        == "????AAAAA???AA"
+    )
+    assert (
+        str(
+            aln.get_seq("x").with_masked_annotations(
+                ["exon", "repeat"], mask_char="?", shadow=True
+            )
+        )
+        == "CCCC?????GGG??"
+    )
+    assert (
+        str(aln.get_seq("y").with_masked_annotations("exon", mask_char="?"))
+        == "TTTTTGGTT"
+    )
+    assert (
+        str(aln.get_seq("y").with_masked_annotations("repeat", mask_char="?"))
+        == "TTTTT??TT"
+    )
+    assert (
+        str(
+            aln.get_seq("y").with_masked_annotations(
+                "repeat", mask_char="?", shadow=True
+            )
+        )
+        == "?????GG??"
+    )
+
+    # The same methods can be applied to annotated Alignment's.
+
+    assert aln.with_masked_annotations("exon", mask_char="?").to_dict() == {
+        "x": "?-???AAAAAGGGAA",
+        "y": "-T----TTTTG-GTT",
+    }
+    assert aln.with_masked_annotations(
+        "exon", mask_char="?", shadow=True
+    ).to_dict() == {"x": "C-CCC??????????", "y": "-?----?????-???"}
+    assert aln.with_masked_annotations("repeat", mask_char="?").to_dict() == {
+        "x": "C-CCCAAAAA???AA",
+        "y": "-T----TTTT?-?TT",
+    }
+    assert aln.with_masked_annotations(
+        "repeat", mask_char="?", shadow=True
+    ).to_dict() == {"x": "?-????????GGG??", "y": "-?----????G-G??"}
+    assert aln.with_masked_annotations(["repeat", "exon"], mask_char="?").to_dict() == {
+        "x": "?-???AAAAA???AA",
+        "y": "-T----TTTT?-?TT",
+    }
+    assert aln.with_masked_annotations(["repeat", "exon"], shadow=True).to_dict() == {
+        "x": "C-CCC?????GGG??",
+        "y": "-?----????G-G??",
+    }
+
+
+def test_nested_annotated_region_masks():
+    """masking a sequence with specific features when nested annotations"""
+    # gene = aln.get_seq("x").add_feature("gene", "norwegian", [(0, 4)])
+    # gene.add_feature("repeat", "blue", [(1, 3)])
+    # exon = aln.get_seq("y").add_feature("repeat", "frog", [(1, 4)])
+    db = GffAnnotationDb()
+    db.add_feature(seqid="x", biotype="gene", name="norwegian", spans=[(0, 4)])
+    db.add_feature(seqid="x", biotype="repeat", name="blue", spans=[(1, 3)])
+    db.add_feature(seqid="y", biotype="repeat", name="frog", spans=[(1, 4)])
+    aln = make_aligned_seqs(
+        data=[["x", "C-GGCAAAAATTTAA"], ["y", "-T----TTTTG-GTT"]], array_align=False
+    )
+    aln.annotation_db = db
+    gene = list(aln.get_seq("x").get_features(biotype="gene"))[0]
+    assert str(gene.get_slice()) == "CGGC"
+
+    # evaluate the sequence directly
+    masked = str(aln.get_seq("x").with_masked_annotations("repeat", mask_char="?"))
+    assert masked == "C??CAAAAATTTAA"
+
+    exon = list(aln.get_seq("y").get_features(biotype="repeat", name="frog"))[0]
+    assert str(exon.get_slice()) == "TTT"
+    # evaluate the sequence directly
+    masked = str(aln.get_seq("y").with_masked_annotations("repeat", mask_char="?"))
+    assert masked == "T???TGGTT"
+    masked = aln.with_masked_annotations("gene", mask_char="?")
+    got = masked.to_dict()
+    assert got["x"] == "?-???AAAAATTTAA"
+    assert got["y"] == "-T----TTTTG-GTT"
+
+
+def test_feature_from_alignment():
+    """seq features obtained from the alignment"""
+
+    # we no longer support copying annotations individually
+    # nor do we provide a mechanism for copying annotations from one
+    # sequence to another
+
+    # Sequence features can be accessed via a containing Alignment:
+    db = GffAnnotationDb()
+    aln = make_aligned_seqs(
+        data={"x": "-AAAAAAAAA", "y": "TTTT--TTTT"}, array_align=False
+    )
+    db.add_feature(seqid="x", biotype="exon", name="fred", spans=[(3, 8)])
+    aln.annotation_db = db
+    aln_exons = list(aln.get_features(seqid="x", biotype="exon"))
+    assert len(aln_exons) == 1
+    aln_exons = aln_exons[0]
+
+    # But these will be returned as **alignment**
+    # features with locations in alignment coordinates.
+    assert aln_exons.get_slice().to_dict() == {"x": "AAAAA", "y": "--TTT"}
+
+    # Similarly alignment features can be projected onto the aligned sequences,
+    # where they may end up falling across gaps:
+
+    exons = aln.get_projected_annotations(seqid="y", biotype="exon")
+    assert len(exons) == 1
+    assert str(aln.get_seq("y")[exons[0].map.without_gaps()]), "TTT"
+    assert 'exon "fred" at [-2-, 4:7]/8' in str(exons[0])
+
+
+def test_nested_get_slice():
+    """check the get_slice method works on nested annotations"""
+    s = DNA.make_seq("AAGAAGAAGACCCCCAAAAAAAAAATTTTTTTTTTAAAAAAAAAAAAA", name="Orig")
+    ex = s.add_feature(biotype="exon", name="fred", spans=[(10, 20)])
+    s.add_feature(biotype="exon", name="trev", spans=[(30, 40)])
+    s.add_feature(biotype="repeat", name="bob", spans=[(12, 17)], parent_id="fred")
+    f = list(ex.get_children())[0]
+    assert str(f.get_slice()) == str(s[12:17])
diff --git a/tests/test_core/test_maps.py b/tests/test_core/test_maps.py
index 444565f99..525ef1234 100644
--- a/tests/test_core/test_maps.py
+++ b/tests/test_core/test_maps.py
@@ -1,9 +1,5 @@
-#!/usr/bin/env python
-
 import unittest
 
-import pytest
-
 from cogent3 import DNA, make_aligned_seqs
 from cogent3.core.location import Map, Span
 
@@ -59,26 +55,27 @@ class MapTest(unittest.TestCase):
         self.assertEqual(len(answer), 1)
         self.assertEqual(str(seq[answer[0]]), "TCGAT")
 
-    @pytest.mark.xfail(reason="todo: see comments in test")
-    def test_get_by_seq_annotation(self):
-        aln = make_aligned_seqs(
-            data={"a": "ATCGAAATCGAT", "b": "ATCGA--TCGAT"}, array_align=False
-        )
-        # original version was putting annotation directly on seq
-        b = aln.get_seq("b")
-        f = aln.add_feature(
-            seqid="b", biotype="test_type", name="test_label", spans=[(4, 6)]
-        )
-        print(f)  # f should be a feature, but is currently None
 
-        # in the original get_by_seq_annotations(), inside the method it
-        # used self[feature.map.start : feature.map.end],
-        # so what was returned was an alignment slice compared to below,
-        # which is only non-gap positions of "b"
-        # todo gah add argument get_slice(strict=True), which means if parent is
-        # an alignment, only the columns from non-gap positions are returned.
-        # Argument has non effect on Sequence, OR add a new get_alignment_slice()
-        # method?
-        answer = list(aln.get_features(seqid="b", biotype="test_type"))
-        answer = answer[0].get_slice().to_dict()
-        self.assertEqual(answer, {"b": "A--T", "a": "AAAT"})
+def test_get_by_seq_annotation_allow_gaps():
+    aln = make_aligned_seqs(
+        data={"a": "ATCGAAATCGAT", "b": "ATCGA--TCGAT"}, array_align=False
+    )
+    # original version was putting annotation directly on seq
+    f = aln.add_feature(
+        seqid="b",
+        biotype="test_type",
+        name="test_label",
+        spans=[(4, 6)],
+        on_alignment=False,
+    )
+    # in the original get_by_seq_annotations(), inside the method it
+    # used self[feature.map.start : feature.map.end],
+    # so what was returned was an alignment slice compared to below,
+    # which is only non-gap positions of "b"
+    # when we set allow_gaps, we get gaps in an alignment
+    f = list(aln.get_features(seqid="b", biotype="test_type"))[0]
+    got = f.get_slice(allow_gaps=True).to_dict()
+    assert got == {"b": "A--T", "a": "AAAT"}
+    # otherwise we only get the exact positions represented
+    got = f.get_slice(allow_gaps=False).to_dict()
+    assert got == {"b": "AT", "a": "AT"}
diff --git a/tests/test_core/test_sequence.py b/tests/test_core/test_sequence.py
index fd439bafb..01888a02f 100644
--- a/tests/test_core/test_sequence.py
+++ b/tests/test_core/test_sequence.py
@@ -230,51 +230,6 @@ class SequenceTests(TestCase):
         with self.assertRaises(ValueError):
             s.to_moltype("")
 
-    @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
-        #          1111111111222222222333333333334
-        # 1234567890123456789012345678901234567890
-        #  **** biological_region
-        #                                     ** biological_region
-        #                                       * biological_region
-        #      *******************************  gene
-        #         *********************   mRNA
-        #            *********            exon
-        #                       *****     exon
-        # ACCCCGGAAAATTTTTTTTTAAGGGGGAAAAAAAAACCCCCCC...
-        seq = DNA.make_seq("ACCCCGGAAAATTTTTTTTTAAGGGGGAAAAAAAAACCCCCCC", name="22")
-        gff3_path = DATADIR / "ensembl_sample.gff3"
-        seq.annotate_from_gff(gff3_path)
-        # we have 1 "full chromosome" annotation, 3 generic regions and 1 gene
-        self.assertEqual(len(seq.annotations), 5)
-
-        # get the gene and check it has a single annotation and that
-        # its slice is correct
-        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_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_features_matching("exon", extend_query=True)
-        self.assertEqual(len(ann), 2)
-        exon_seqs = ("TTTTTTTTT", "GGGGG")
-        for x in ann:
-            self.assertEqual(len(x.annotations), 0)
-            self.assertTrue(str(x.get_slice()) in exon_seqs, msg=x.get_slice())
-
     def test_strip_degenerate(self):
         """Sequence strip_degenerate should remove any degenerate bases"""
         self.assertEqual(self.RNA("UCAG-").strip_degenerate(), "UCAG-")
@@ -2124,3 +2079,46 @@ def test_absolute_relative_roundtrip_reverse(
     rel_val = view.relative_position(abs_val)
     assert view.offset == offset
     assert (view[rel_val]).value == view[value].value
+
+
+def test_annotate_gff_nested_features():
+    """correctly annotate a sequence with nested features"""
+    # the synthetic example
+    #          1111111111222222222333333333334
+    # 1234567890123456789012345678901234567890
+    #  **** biological_region
+    #                                     ** biological_region
+    #                                       * biological_region
+    #      *******************************  gene
+    #         *********************   mRNA
+    #            *********            exon
+    #                       *****     exon
+    # ACCCCGGAAAATTTTTTTTTAAGGGGGAAAAAAAAACCCCCCC...
+    seq = DNA.make_seq("ACCCCGGAAAATTTTTTTTTAAGGGGGAAAAAAAAACCCCCCC", name="22")
+    gff3_path = DATADIR / "ensembl_sample.gff3"
+    seq.annotate_from_gff(gff3_path)
+    # we have 8 records in the gff file
+    assert seq.annotation_db.num_matches() == 8
+
+    # get the gene and check it has a single annotation and that
+    # its slice is correct
+    ann = list(seq.get_features(biotype="gene"))
+    assert len(ann) == 1
+    ann_seq = ann[0].get_slice()
+    assert str(ann_seq) == "GGAAAATTTTTTTTTAAGGGGGAAAAAAAAA"
+    # the gene has 1 transcript
+    gene = ann[0]
+    mrna = list(gene.get_children(biotype="mRNA"))
+    assert len(mrna) == 1
+    mrna = mrna[0]
+    ann_seq = mrna.get_slice()
+    assert str(ann_seq) == "AAAATTTTTTTTTAAGGGGGAAA"
+
+    # the transcript has 2 exons, from the parent feature
+    exons = list(mrna.get_children(biotype="exon"))
+    assert len(exons) == 2
+    # or the sequence
+    ann = list(seq.get_features(biotype="exon"))
+    assert len(ann) == 2
+    exon_seqs = ("TTTTTTTTT", "GGGGG")
+    assert tuple(str(ex.get_slice()) for ex in exons) == exon_seqs