Repository: cogent3/cogent3
Pull request #1400 (2023-04-28 03:50:16)

NEW: get_feature returns feature in relative coordinates


Patch:
diff --git a/src/cogent3/core/alignment.py b/src/cogent3/core/alignment.py
index 2f5522a01..7c91ce866 100644
--- a/src/cogent3/core/alignment.py
+++ b/src/cogent3/core/alignment.py
@@ -1712,7 +1712,6 @@ class _SequenceCollectionBase:
         )
         if isinstance(self, _Annotatable) and self.annotations:
             self._annotations_nucleic_reversed_on(rc)
-        rc.annotation_db = self.annotation_db
         return rc
 
     def reverse_complement(self):
diff --git a/src/cogent3/core/annotation.py b/src/cogent3/core/annotation.py
index 45f56785f..d39b94aea 100644
--- a/src/cogent3/core/annotation.py
+++ b/src/cogent3/core/annotation.py
@@ -7,7 +7,6 @@ from typing import Optional
 
 import numpy
 
-from cogent3.util import warning as c3warn
 from cogent3.util.misc import get_object_provenance
 
 from .location import Map, as_map
@@ -167,6 +166,39 @@ class _AnnotationCore:
     def _mapped(self, map):
         raise NotImplementedError
 
+    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 = []
+        annotation_types = []
+        for annot in annotations:
+            spans.extend(annot.map.spans)
+            if annot.type not in annotation_types:
+                annotation_types.append(annot.type)
+        map = Map(spans=spans, parent_length=len(self))
+        map = map.covered()  # No overlaps
+        name = ",".join(annotation_types)
+
+        if feature_class is None:
+            feature_class = _Feature
+
+        return feature_class(self, map, type="region", name=name)
+
     def _annotations_nucleic_reversed_on(self, new):
         """applies self.annotations to new with coordinates adjusted for
         reverse complement."""
@@ -424,7 +456,6 @@ class Annotation(_AnnotationCore, _Serialisable):
         self.biotype = biotype
         self.name = name
         self.parent = parent
-        self.seqid = seqid
         if isinstance(map, Map):
             assert map.parent_length == len(parent), (map, len(parent))
         else:
@@ -462,43 +493,18 @@ class Annotation(_AnnotationCore, _Serialisable):
             return self
         keep = self.map.nongap()
         return self.__class__(
-            self.parent,
-            self.map[keep],
-            biotype=self.biotype,
-            name=self.name,
-            seqid=self.seqid,
+            self.parent, self.map[keep], biotype=self.biotype, name=self.name
         )
 
     def as_one_span(self):
         new_map = self.map.get_covering_span()
         return self.__class__(
-            self.parent,
-            new_map,
-            biotype=self.biotype,
-            name=f"one-span {self.name}",
-            seqid=self.seqid,
+            self.parent, new_map, biotype=self.biotype, name=f"one-span {self.name}"
         )
 
-    def shadow(self):
-        """returns new instance corresponding to disjoint of self coordinates"""
-        return self.__class__(
-            parent=self.parent,
-            map=self.map.shadow(),
-            biotype=f"not {self.biotype}",
-            name=f"not {self.name}",
-            seqid=self.seqid,
-        )
-
-    @c3warn.deprecated_callable(
-        "2023.7", reason="new method", new="<instance>.shadow()"
-    )
     def get_shadow(self):
         return self.__class__(
-            self.parent,
-            self.map.shadow(),
-            type="region",
-            name=f"not {self.name}",
-            seqid=self.seqid,
+            self.parent, self.map.shadow(), type="region", name=f"not {self.name}"
         )
 
     def __len__(self):
@@ -506,7 +512,7 @@ class Annotation(_AnnotationCore, _Serialisable):
 
     def __repr__(self):
         name = f' "{self.name}"'
-        return f"{self.seqid!r} {self.biotype}{name} at {self.map}"
+        return f"{self.biotype}{name} at {self.map}"
 
     def _projected_to_base(self, base):
         if self.parent == base:
@@ -532,37 +538,6 @@ class Annotation(_AnnotationCore, _Serialisable):
         for child in db.get_feature_children(biotype=biotype, name=self.name):
             yield make_feature(child)
 
-    def union(self, features):
-        """return as a single Annotation
-
-        Notes
-        -----
-        Overlapping spans are merged
-        """
-        combined = self.map.spans[:]
-        feat_names = set()
-        biotypes = set()
-        seqids = set()
-        for feature in features:
-            if feature.parent is not self.parent:
-                raise ValueError(f"cannot merge annotations from different objects")
-
-            combined.extend(feature.map.spans)
-            if feature.name:
-                feat_names.add(feature.name)
-            if feature.seqid:
-                seqids.add(feature.seqid)
-            if feature.biotype:
-                biotypes.add(feature.biotype)
-        name = ", ".join(feat_names)
-        map = Map(spans=combined, parent_length=len(self.parent))
-        map = map.covered()  # No overlaps
-        seqid = ", ".join(seqids) if seqids else None
-        biotype = ", ".join(biotypes)
-        return self.__class__(
-            parent=self.parent, seqid=seqid, map=map, biotype=biotype, name=name
-        )
-
 
 # https://pythonspeed.com/products/filmemoryprofiler/
 
@@ -689,42 +664,6 @@ class _Feature(_Annotatable, _Serialisable):
         annotatable.attach_annotations([new])
         return annotatable
 
-    @c3warn.deprecated_callable(
-        "2023.7", reason="replaced by", new="<instance>.union()"
-    )
-    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 = []
-        annotation_types = []
-        for annot in annotations:
-            spans.extend(annot.map.spans)
-            if annot.type not in annotation_types:
-                annotation_types.append(annot.type)
-        map = Map(spans=spans, parent_length=len(self))
-        map = map.covered()  # No overlaps
-        name = ",".join(annotation_types)
-
-        if feature_class is None:
-            feature_class = _Feature
-
-        return feature_class(self, map, type="region", name=name)
-
 
 class AnnotatableFeature(_Feature):
     """These features can themselves be annotated."""
diff --git a/src/cogent3/core/location.py b/src/cogent3/core/location.py
index d0f334d2a..89733bba5 100644
--- a/src/cogent3/core/location.py
+++ b/src/cogent3/core/location.py
@@ -575,13 +575,35 @@ class Map(object):
                     raise RuntimeError(
                         f"located outside sequence: {str((start, end, parent_length))}"
                     )
-                span = Span(start, end, tidy, tidy, reverse=reverse)
-                if diff < 0:
-                    spans += [LostSpan(-diff), span]
-                elif diff > 0:
-                    spans += [span, LostSpan(diff)]
+                if max(start, end) > parent_length and min(start, end) < 0:
+                    l_diff = min(start, end)
+                    r_diff = max(start, end) - parent_length
+                    start, end = (
+                        (0, parent_length) if start < end else (parent_length, 0)
+                    )
+                    spans += [
+                        LostSpan(abs(l_diff)),
+                        Span(start, end, tidy, tidy, reverse=reverse),
+                        LostSpan(abs(r_diff)),
+                    ]
+                elif min(start, end) < 0:
+                    diff = min(start, end)
+                    start = 0 if start < 0 else start
+                    end = 0 if end < 0 else end
+                    spans += [
+                        LostSpan(abs(diff)),
+                        Span(start, end, tidy, tidy, reverse=reverse),
+                    ]
+                elif max(start, end) > parent_length:
+                    diff = max(start, end) - parent_length
+                    start = parent_length if start > parent_length else start
+                    end = parent_length if end > parent_length else end
+                    spans += [
+                        Span(start, end, tidy, tidy, reverse=reverse),
+                        LostSpan(abs(diff)),
+                    ]
                 else:
-                    spans += [span]
+                    spans += [Span(start, end, tidy, tidy, reverse=reverse)]
 
         self.offsets = []
         self.useful = False
@@ -690,7 +712,7 @@ class Map(object):
             last_x = x
             last_y = y
         assert y == 0
-        return Map(locations=result, parent_length=self.parent_length)
+        return Map(result, parent_length=self.parent_length)
 
     def reversed(self):
         """Reversed location on same parent"""
diff --git a/src/cogent3/core/sequence.py b/src/cogent3/core/sequence.py
index f389f1a08..60662eb6b 100644
--- a/src/cogent3/core/sequence.py
+++ b/src/cogent3/core/sequence.py
@@ -20,7 +20,7 @@ import warnings
 from functools import singledispatch, total_ordering
 from operator import eq, ne
 from random import shuffle
-from typing import Generator, List, Optional, Tuple, Union
+from typing import Generator, List, Optional, Tuple
 
 from numpy import (
     arange,
@@ -933,8 +933,14 @@ class Sequence(_Annotatable, SequenceI):
         # we set include_boundary=True because stop is exclusive indexing,
         # i,e., the stop can be equal to the length of the view
         query_end = self._seq.absolute_position(stop, include_boundary=True)
-        # todo gah check logic of handling of negative indices
-        query_start = max(query_start, 0)
+
+        # if the view is reversed, then (query_start > query_end)
+        query_start, query_end = (
+            (query_start, query_end)
+            if query_start < query_end
+            else (query_end, query_start)
+        )
+
         for feature in self.annotation_db.get_features_matching(
             seqid=self.name,
             name=name,
@@ -946,8 +952,34 @@ class Sequence(_Annotatable, SequenceI):
             # spans need to be converted from absolute to relative positions
             # DO NOT do adjustment in make_feature since that's user facing,
             # and we expect them to make a feature manually
+
+            feature["spans"] = self._relative_spans(feature["spans"])
+
             yield self.make_feature(feature)
 
+    def _relative_spans(self, spans):
+        r_spans = []
+
+        for (s, e) in spans:
+            # reverse feature determined by absolute position
+            reverse = s > e
+
+            # todo: kath, think about logic of stop=true for reverse locations
+            # I think, we need to know the whether it is reverse so that the adjustment
+            # for step (given stop=True) is in the correct direction.
+            # However, we need to keep s > e for a reverse span because this information
+            # is used in the Map?
+            if reverse:
+                start = self._seq.relative_position(s, stop=True)
+                end = self._seq.relative_position(e)
+            else:
+                start = self._seq.relative_position(s)
+                end = self._seq.relative_position(e, stop=True)
+
+            r_spans += [(start, end)]
+
+        return r_spans
+
     @deprecated_callable(
         "2023.7", reason="simpler name", new="<instance>.get_features()"
     )
@@ -962,7 +994,7 @@ class Sequence(_Annotatable, SequenceI):
         Parameters
         ----------
         feature
-            dict of key data to make a Annotation instance
+            dict of key data to make an Annotation instance
 
         Notes
         -----
@@ -1432,7 +1464,7 @@ class NucleicAcidSequence(Sequence):
         )
         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):


Test patch:
diff --git a/tests/test_core/test_annotation_db.py b/tests/test_core/test_annotation_db.py
index b10010752..a143ffd75 100644
--- a/tests/test_core/test_annotation_db.py
+++ b/tests/test_core/test_annotation_db.py
@@ -549,3 +549,29 @@ def test_gff_update_existing_specify_seqid(gff_db, gff_small_db, seqids):
     )
     gff_db.update(gff_small_db, seqids=seqids)
     assert gff_db.num_matches() == expect
+
+
+def test_relative_position_negative_feature(seq_db):
+    orig_feat_span = list(seq_db.get_features(name="Transcript:B0019.1"))[0].map
+
+    view = seq_db[5:]
+    view_feat_span = list(view.get_features(name="Transcript:B0019.1"))[0].map
+
+    assert orig_feat_span[0].start - 5 == view_feat_span[0].start
+    assert orig_feat_span[0].end - 5 == view_feat_span[0].end
+
+
+def test_relative_position_positive_feature(anno_db):
+    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="+"
+    )
+
+    orig_feat_span = list(seq.get_features_matching(name="exon1"))[0].map
+    view = seq[2:]
+    view_feat_span = list(view.get_features_matching(name="exon1"))[0].map
+
+    assert orig_feat_span[0].start - 2 == view_feat_span[0].start
+    assert orig_feat_span[0].end - 2 == view_feat_span[0].end
diff --git a/tests/test_core/test_features.py b/tests/test_core/test_features.py
index ea22bf34e..a0011ca6f 100755
--- a/tests/test_core/test_features.py
+++ b/tests/test_core/test_features.py
@@ -89,32 +89,85 @@ class FeaturesTest(TestCase):
         self.assertEqual(dont_exist, [])
         self.assertEqual(str(self.s[dont_exist]), "")
 
-    def test_union(self):
-        """combines multiple features into one"""
+    @pytest.mark.xfail(reason="todo gah duplicated in test_annotation?")
+    def test_get_region_covering_all(self):
+        """combines multiple features into one or their shadow"""
 
         # To construct a pseudo-feature covering (or excluding)
         # multiple features, use get_region_covering_all:
 
-        exons = list(self.s.get_features(biotype="exon"))
-        exon1 = exons.pop(0)
-        combined = exon1.union(exons)
-        self.assertEqual(str(combined.get_slice()), "CCCCCTTTTTAAAAA")
+        exons = self.s.get_features(biotype="exon")
+        self.assertEqual(
+            str(self.s.get_region_covering_all(exons)),
+            'region "exon" at [10:15, 30:40]/48',
+        )
+        self.assertEqual(
+            str(self.s.get_region_covering_all(exons).get_shadow()),
+            'region "not exon" at [0:10, 15:30, 40:48]/48',
+        )
 
-    def test_shadow(self):
-        """combines multiple features into shadow"""
+        # eg: all the exon sequence:
 
-        # To construct a pseudo-feature covering (or excluding)
-        # multiple features, use get_region_covering_all:
+        self.assertEqual(
+            str(self.s.get_region_covering_all(exons).get_slice()), "CCCCCTTTTTAAAAA"
+        )
 
-        exons = list(self.s.get_features(biotype="exon"))
-        expect = str(
-            self.s[: exons[0].map.start]
-            + self.s[exons[0].map.end : exons[1].map.start]
-            + self.s[exons[1].map.end :]
-        )
-        exon1 = exons.pop(0)
-        shadow = exon1.union(exons).shadow()
-        assert str(shadow.get_slice()) == expect
+        # or with slice notation:
+
+        self.assertEqual(str(self.s[self.exon1, self.exon2]), "CCCCCTTTTTAAAAA")
+
+    @pytest.mark.xfail(reason="todo gah update test to use latest API")
+    def test_slice_errors_from_merged(self):
+        """no overlap in merged features"""
+
+        # Though .get_region_covering_all also guarantees
+        # no overlaps within the result, slicing does not:
+
+        exons = self.s.get_features(biotype="exon")
+        self.assertEqual(
+            str(self.s.get_region_covering_all(exons + exons)),
+            'region "exon" at [10:15, 30:40]/48',
+        )
+        with self.assertRaises(ValueError):
+            self.s[self.exon1, self.exon1, self.exon1, self.exon1, self.exon1]
+
+        # You can use features, maps, slices or integers,
+        # but non-monotonic slices are not allowed:
+
+        with self.assertRaises(ValueError):
+            self.s[15:20, 5:16]
+
+    @pytest.mark.xfail(reason="todo gah delete test, not supporting this")
+    def test_feature_slicing(self):
+        """features can be sliced"""
+
+        # Features are themselves sliceable:
+
+        self.assertEqual(str(self.exon1[0:3].get_slice()), "CCC")
+
+        # When sequences are concatenated they keep their (non-overlapping) annotations:
+
+        c = self.s[self.exon1[4:]] + self.s
+        self.assertEqual(len(c), 49)
+
+        answers = [
+            'exon "fred" at [-4-, 0:1]/49',
+            'exon "fred" at [11:16]/49',
+            'exon "trev" at [31:41]/49',
+        ]
+
+        array = []
+
+        for feat in c.annotations:
+            array.append(str(feat))
+
+        self.assertEqual(answers, array)
+
+        # Since features know their parents you can't
+        # use a feature from one sequence to slice another:
+
+        with self.assertRaises(ValueError):
+            c[self.exon1]
 
     @pytest.mark.xfail(reason="todo gah delete test? or update to new API")
     def test_feature_attach_detach(self):