Repository: cogent3/cogent3
Pull request #1377 (2023-04-21 06:49:09)

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 5aa546167..2751a2e0d 100644
--- a/src/cogent3/core/alignment.py
+++ b/src/cogent3/core/alignment.py
@@ -2010,7 +2010,7 @@ class SequenceCollection(_SequenceCollectionBase):
             if name in self.named_seqs:
                 self.named_seqs[name].copy_annotations(seq)
 
-    def annotate_from_gff(self, f: os.PathLike, seq_ids: [str]):
+    def annotate_from_gff(self, f: os.PathLike, seq_ids: list[str]):
         """copies annotations from a gff file to a sequence in self
 
         Parameters
diff --git a/src/cogent3/core/annotation.py b/src/cogent3/core/annotation.py
index b5abac1cf..2298709b6 100644
--- a/src/cogent3/core/annotation.py
+++ b/src/cogent3/core/annotation.py
@@ -3,6 +3,7 @@ import json
 
 from collections import defaultdict
 from fnmatch import fnmatch
+from typing import Optional
 
 import numpy
 
@@ -118,12 +119,9 @@ class _AnnotationMixin:
         return drawer
 
 
-class _Annotatable(_AnnotationMixin):
-    # default
-    annotations = ()
-
-    # Subclasses should provide __init__, getOwnTracks, and a _mapped for use by
-    # __getitem__
+class _AnnotationCore:
+    """core methods for annotation handling, used by both new and old
+    style features"""
 
     def _sliced_annotations(self, new, slice):
         result = []
@@ -147,15 +145,6 @@ class _Annotatable(_AnnotationMixin):
                             result.append(annot)
         return result
 
-    def _shifted_annotations(self, new, shift):
-        result = []
-        if self.annotations:
-            newmap = Map([(shift, shift + len(self))], parent_length=len(new))
-            for annot in self.annotations:
-                annot = annot.remapped_to(new, newmap)
-                result.append(annot)
-        return result
-
     def _as_map(self, index):
         """Can take a slice, integer, map or feature, or even a list/tuple of those"""
         if type(index) in [list, tuple]:
@@ -193,6 +182,66 @@ class _Annotatable(_AnnotationMixin):
     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."""
+        assert len(new) == len(self)
+        annotations = []
+        for annot in self.annotations:
+            new_map = annot.map.nucleic_reversed()
+            annotations.append(annot.__class__(new, new_map, annot))
+        new.attach_annotations(annotations)
+
+
+class _Annotatable(_AnnotationCore, _AnnotationMixin):
+    # default
+    annotations = ()
+
+    # Subclasses should provide __init__, getOwnTracks, and a _mapped for use by
+    # __getitem__
+
+    def _shifted_annotations(self, new, shift):
+        result = []
+        if self.annotations:
+            newmap = Map([(shift, shift + len(self))], parent_length=len(new))
+            for annot in self.annotations:
+                annot = annot.remapped_to(new, newmap)
+                result.append(annot)
+        return result
+
     def add_annotation(self, klass, *args, **kw):
         from cogent3.core.sequence import Sequence
 
@@ -313,39 +362,6 @@ class _Annotatable(_AnnotationMixin):
                 )
         return result
 
-    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 get_by_annotation(self, annotation_type, name=None, ignore_partial=False):
         """yields the sequence segments corresponding to the specified
         annotation_type and name one at a time.
@@ -380,16 +396,6 @@ class _Annotatable(_AnnotationMixin):
             seq.info["name"] = annotation.name
             yield seq
 
-    def _annotations_nucleic_reversed_on(self, new):
-        """applies self.annotations to new with coordinates adjusted for
-        reverse complement."""
-        assert len(new) == len(self)
-        annotations = []
-        for annot in self.annotations:
-            new_map = annot.map.nucleic_reversed()
-            annotations.append(annot.__class__(new, new_map, annot))
-        new.attach_annotations(annotations)
-
     def _projected_to_base(self, base):
         raise NotImplementedError
 
@@ -423,6 +429,96 @@ class _Serialisable:
         return json.dumps(self.to_rich_dict())
 
 
+class Annotation(_AnnotationCore, _Serialisable):
+    """new style annotation, created on demand"""
+
+    __slots__ = "parent", "map", "biotype", "name", "_serialisable", "base", "base_map"
+
+    # todo implement a __new__ to trap args for serialisation purposes
+    def __init__(self, parent, map: Map, biotype: str, name: str):
+        d = locals()
+        exclude = ("self", "__class__", "kw")
+        self._serialisable = {k: v for k, v in d.items() if k not in exclude}
+        self.biotype = biotype
+        self.name = name
+        self.parent = parent
+        if isinstance(map, Map):
+            assert map.parent_length == len(parent), (map, len(parent))
+        else:
+            map = Map(locations=map, parent_length=len(parent))
+
+        self.map = map
+
+    def get_drawable(self):
+        """returns plotly trace"""
+        from cogent3.draw.drawable import make_shape
+
+        return make_shape(type_=self)
+
+    def _mapped(self, slicemap):
+        name = f"{repr(slicemap)} of {self.name}"
+        return self.__class__(
+            parent=self.parent, map=slicemap, biotype=self.biotype, name=name
+        )
+
+    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."""
+        map = self.map
+        if not (complete or map.complete):
+            map = map.without_gaps()
+        return self.parent[map]
+
+    def without_lost_spans(self):
+        """Keeps only the parts which are actually present in the underlying sequence"""
+        if self.map.complete:
+            return self
+        keep = self.map.nongap()
+        return self.__class__(
+            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}"
+        )
+
+    def get_shadow(self):
+        return self.__class__(
+            self.parent, self.map.shadow(), type="region", name=f"not {self.name}"
+        )
+
+    def __len__(self):
+        return len(self.map)
+
+    def __repr__(self):
+        name = f' "{self.name}"'
+        return f"{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):
+        map = gmap[self.map]
+        return self.__class__(grandparent, map, biotype=self.biotype, name=self.name)
+
+    def get_coordinates(self):
+        """returns sequence coordinates of this Feature as
+        [(start1, end1), ...]"""
+        return self.map.get_coordinates()
+
+    def get_children(self, biotype: Optional[str] = None):
+        """generator returns sub-features of self optionally matching biotype"""
+        make_feature = self.parent.make_feature
+        db = self.parent.annotation_db
+        for child in db.get_feature_children(biotype=biotype, name=self.name):
+            yield make_feature(child)
+
+
 class FeatureNew(_AnnotationMixin, _Serialisable):
     def __init__(self, parent, biotype, name, spans, reversed):
         self.parent = parent
@@ -431,6 +527,28 @@ class FeatureNew(_AnnotationMixin, _Serialisable):
         self.reversed = reversed
         self.map = Map(locations=spans, parent_length=len(parent))
 
+    # remapped_to takes an alignment and the "inverted" Aligned map and returns a new instance of self
+    def from_seq_to_align(self, aln, aligned):
+        """transfer instance into aln through alignment coords
+
+        Parameters
+        ----------
+        aln : Alignment
+            self should be bound to an Aligned member of obj
+        """
+        assert len(self.parent) == len(aln)
+        # the Map.inverse() method transforms the seq coords into alignment
+        # coords, returning a new Map which can itself be sliced
+        mapp = aligned.map.inverse()
+        new_map = mapp[self.map]
+        return self.__class__(
+            aln,
+            biotype=self.biotype,
+            name=self.name,
+            spans=new_map.get_coordinates(),
+            reversed=self.reversed,
+        )
+
     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
@@ -438,6 +556,7 @@ class FeatureNew(_AnnotationMixin, _Serialisable):
 
         # todo: need to cope with cases where span is a lostSpan...
 
+        # separation of concerns issue, this prevents slicing using an alignment
         parent_reversed = self.parent._seq.reversed
         feature_reversed = self.reversed
 
@@ -470,7 +589,7 @@ class FeatureNew(_AnnotationMixin, _Serialisable):
         for child in self.parent.annotation_db.get_feature_children(
             biotype=biotype, name=self.name
         ):
-            yield FeatureNew(self.parent, **child)
+            yield Annotation(self.parent, **child)
 
     def __len__(self):
         return len(self.map)
@@ -697,3 +816,7 @@ def SimpleVariable(parent, type, name, data):
     assert len(data) == len(parent), (len(data), len(parent))
     map = Map([(0, len(data))], parent_length=len(parent))
     return _SimpleVariable(parent, map, type=type, name=name, data=data)
+
+
+def make_annotation_for_seq(seq, feature):
+    ...
diff --git a/src/cogent3/core/sequence.py b/src/cogent3/core/sequence.py
index 60e89256c..e1f204889 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, Tuple
+from typing import Generator, List, Optional, Tuple, Union
 
 from numpy import (
     arange,
@@ -37,8 +37,9 @@ from numpy import (
 from numpy.random import permutation
 
 from cogent3.core.alphabet import AlphabetError
-from cogent3.core.annotation import FeatureNew, Map, _Annotatable
+from cogent3.core.annotation import Annotation, Map, _Annotatable
 from cogent3.core.annotation_db import (
+    FeatureDataType,
     GenbankAnnotationDb,
     GffAnnotationDb,
     load_annotations,
@@ -56,7 +57,7 @@ from cogent3.util.misc import (
     get_setting_from_environ,
 )
 from cogent3.util.transform import for_seq, per_shortest
-from cogent3.util.warning import deprecated
+from cogent3.util.warning import deprecated, deprecated_args
 
 
 __author__ = "Rob Knight, Gavin Huttley, and Peter Maxwell"
@@ -876,14 +877,32 @@ class Sequence(_Annotatable, SequenceI):
             raise TypeError
         self._annotation_db = value
 
+    # todo deprecate feature_type for biotype
     def get_features_matching(
         self,
-        feature_type=None,
-        name=None,
-        start=None,
-        stop=None,
+        feature_type: Optional[str] = None,
+        name: Optional[str] = None,
+        start: Optional[int] = None,
+        stop: Optional[int] = None,
     ):
-        """yield features matching the given conditions"""
+        """yields Annotation instances
+
+        Parameters
+        ----------
+        feature_type
+            biotype of the feature
+        name
+            name of the feature
+        start, end
+            start, end positions to search between, relative to offset
+            of this sequence. If not provided, entire span of sequence is used.
+
+        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.
+        """
 
         if self._annotation_db is None:
             return None
@@ -891,6 +910,8 @@ class Sequence(_Annotatable, SequenceI):
         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
 
+        seq_rced = self._seq.reversed
+
         for feature in self.annotation_db.get_features_matching(
             seqid=self.name,
             name=name,
@@ -898,8 +919,45 @@ class Sequence(_Annotatable, SequenceI):
             start=query_start,
             end=query_end,
         ):
-            feature.pop("seqid")  # not required here as provided
-            yield FeatureNew(self, **feature)
+            yield self.make_feature(feature)
+
+    def make_feature(self, feature: FeatureDataType) -> Annotation:
+        """
+        return an Annotation instance from feature data
+
+        Parameters
+        ----------
+        feature
+            dict of key data to make a Annotation instance
+
+        Notes
+        -----
+        Unlike add_feature(), this method does not add the feature to the
+        database.
+        """
+        feature = dict(feature)
+        seq_rced = self._seq.reversed
+        feature.pop("seqid", None)  # not required here as provided
+        revd = feature.pop("reversed", None)
+        spans = feature.pop("spans", None)
+        fmap = Map(locations=spans, parent_length=len(self))
+        if revd and not seq_rced:
+            # the sequence is on the plus strand, and the
+            # feature coordinates are also for the plus strand
+            # but their order needs to be changed to indicate
+            # reverse complement is required
+            fmap = fmap.reversed()
+        elif seq_rced and not revd:
+            # plus strand feature, but the sequence reverse complemented
+            # so we need to nucleic-reverse the map
+            fmap = fmap.nucleic_reversed()
+        elif seq_rced:
+            # sequence is rc'ed, the feature was minus strand of
+            # original, so needs to be both nucleic reversed and
+            # then reversed
+            fmap = fmap.nucleic_reversed().reversed()
+
+        return Annotation(self, fmap, **feature)
 
     def annotate_from_gff(self, f: os.PathLike, offset=None):
         """copies annotations from a gff file to self,
@@ -1051,7 +1109,7 @@ class Sequence(_Annotatable, SequenceI):
                 else:
                     raise ValueError(f"gap(s) in map {map}")
             else:
-                seg = str(self._seq[span.start : span.end])
+                seg = str(self[span.start : span.end])
                 if span.reverse:
                     complement = self.moltype.complement
                     seg = [complement(base) for base in seg[::-1]]
@@ -1080,15 +1138,15 @@ class Sequence(_Annotatable, SequenceI):
         myclass = f"{self.__class__.__name__}"
         myclass = myclass.split(".")[-1]
         if len(self) > 10:
-            seq = f"{str(self._seq[:7])}... {len(self)}"
+            seq = f"{str(self[:7])}... {len(self)}"
         else:
             seq = str(self)
         return f"{myclass}({seq})"
 
     def __getitem__(self, index):
-
-        if hasattr(index, "get_slice"):
-            return index.get_slice()
+        # todo gah check this omission
+        # if hasattr(index, "get_slice"):
+        #     return index.get_slice()
 
         if isinstance(index, Map):
             new = self._mapped(index)
@@ -1545,7 +1603,7 @@ class SeqView:
 
     def absolute_index(self, value):
         # note: this is the positive/positive case...
-
+        # todo gah handle reversed
         if self.start + value > self.stop:
             raise IndexError("Index out of bounds")
 


Test patch:
diff --git a/tests/test_core/test_annotation_db.py b/tests/test_core/test_annotation_db.py
index b1809bc42..d3fb600ee 100644
--- a/tests/test_core/test_annotation_db.py
+++ b/tests/test_core/test_annotation_db.py
@@ -4,12 +4,12 @@ import numpy
 import pytest
 
 from cogent3 import DNA, RNA, SequenceCollection, _Table, load_seq
-from cogent3.core.annotation import FeatureNew
+from cogent3.core.annotation import Annotation
 from cogent3.core.annotation_db import (
     GffAnnotationDb,
     SupportsFeatures,
-    load_annotations,
     _matching_conditions,
+    load_annotations,
 )
 from cogent3.core.sequence import Sequence
 
@@ -210,7 +210,77 @@ def test_get_features_matching_matching_feature(seq, anno_db):
     assert got[0].biotype == "exon"
     assert got[0].name == "exon1"
     assert len(got) == 1
-    assert got[0].reversed == False
+
+
+def test_gav():
+    # gah todo also need to check a slice with stride
+    from cogent3 import make_seq
+    from cogent3.core.annotation import Feature
+
+    #            ++   ++++
+    #              ---    --
+    raw_seq = "AACCTTTGGGGAATTT"
+
+    plus_spans = [(2, 4), (7, 11)]
+    plus_seq = "".join(raw_seq[s:e] for s, e in plus_spans)
+
+    minus_spans = [(4, 7), (11, 13)]
+    minus_seq = "".join(raw_seq[s:e] for s, e in minus_spans)
+    minus_seq = "".join([{"T": "A", "A": "T"}[b] for b in minus_seq[::-1]])
+    seq = make_seq(raw_seq, name="s1", moltype="dna")
+    db = GffAnnotationDb(data=[])
+    db.add_feature(
+        seqid="s1",
+        biotype="cds",
+        name="plus",
+        spans=plus_spans,
+        strand="+",
+        on_alignment=False,
+    )
+    db.add_feature(
+        seqid="s1",
+        biotype="cds",
+        name="minus",
+        spans=minus_spans,
+        strand="-",
+        on_alignment=False,
+    )
+    seq.annotation_db = db
+    plus = list(seq.get_features_matching(name="plus"))[0]
+    assert str(plus.get_slice()) == plus_seq
+    minus = list(seq.get_features_matching(name="minus"))[0]
+    assert str(minus.get_slice()) == minus_seq
+
+    f = Feature(seq, "cds", "minus", reversed([reversed(s) for s in minus_spans]))
+    s = f.get_slice()
+    assert str(s) == minus_seq
+
+    # now reverse complement the sequence
+    rced = seq.rc()
+
+    rced_spans = [(5, 9), (12, 14)]
+    rf = Feature(rced, "cds", "minus", [reversed(s) for s in rced_spans])
+    got = rf.get_slice()
+
+    plus = list(rced.get_features_matching(name="plus"))[0]
+    assert str(plus.get_slice()) == plus_seq
+    minus = list(rced.get_features_matching(name="minus"))[0]
+    print(rf, minus)
+    assert str(minus.get_slice()) == minus_seq
+
+
+def test_gav2():
+    from cogent3 import make_seq
+    from cogent3.core import location as loc
+
+    seq = make_seq("AACCTTTGGGGAATTT", moltype="dna")
+    mmap = loc.Map(locations=[(4, 7), (11, 13)], parent_length=16)
+    expect = seq[mmap.reversed()]
+
+    rcseq = seq.rc()
+    rmap = mmap.nucleic_reversed().reversed()
+    got = rcseq[rmap]
+    assert str(got) == str(expect)
 
 
 def test_get_features_matching_matching_features(anno_db: GffAnnotationDb, seq):
@@ -251,7 +321,7 @@ def test_get_features_matching_start_stop(seq):
 
 
 def test_matching_conditions():
-    got, _ = _matching_conditions({'start': 1, 'end': 5}, partial=True)
+    got, _ = _matching_conditions({"start": 1, "end": 5}, partial=True)
     expect = "((start <= 1 AND end > 1) OR (start < 5 AND end >= 5) OR (start <= 1 AND end >= 5))"
     assert got == expect
 
@@ -273,7 +343,7 @@ def test_get_features_matching_start_stop_seqview(seq):
 def test_get_slice():
     """get_slice should return the same as slicing the sequence directly"""
     seq = Sequence("ATTGTACGCCCCTGA", name="test_seq")
-    feature_dict = {
+    feature_data = {
         "biotype": "CDS",
         "name": "fake",
         "spans": [
@@ -281,8 +351,7 @@ def test_get_slice():
         ],
         "reversed": False,
     }
-
-    feature = FeatureNew(seq, **feature_dict)
+    feature = seq.make_feature(feature_data)
     got = feature.get_slice()
     assert str(got) == str(seq[5:10])
 
diff --git a/tests/test_core/test_sequence.py b/tests/test_core/test_sequence.py
index abed8aa69..c6eb9f9fb 100644
--- a/tests/test_core/test_sequence.py
+++ b/tests/test_core/test_sequence.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 """Unit tests for Sequence class and its subclasses.
 """
 
@@ -1949,6 +1948,29 @@ def test_seqview_repr():
     assert repr(view) == expected
 
 
-# run if called from command-line
-if __name__ == "__main__":
-    main()
+@pytest.fixture(scope="session")
+def one_seq():
+    from cogent3 import make_seq
+
+    return make_seq("AACCTGGAACC", moltype="dna")
+
+
+@pytest.mark.parametrize("rc", (False, True))
+def test_seq_repr(one_seq, rc):
+    pat = re.compile("[ACGT]+")
+    dna = one_seq.moltype
+    if rc:
+        expect = dna.rc(str(one_seq))
+        seq = one_seq.rc()
+    else:
+        expect = str(one_seq)
+        seq = one_seq
+
+    got = pat.findall(repr(seq))[0]
+    assert expect.startswith(got), (expect, got)
+
+
+def test_gav3():
+    raise NotImplementedError(
+        "check annotations get correct sequence when there's a stride"
+    )