Repository: cogent3/cogent3
Pull request #1401 (2023-04-28 04:32:11)

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 2f5522a01..36ad54bf2 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
@@ -2246,11 +2245,21 @@ class SequenceCollection(_SequenceCollectionBase):
             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_alignment:
-                    yield self.make_feature(feature=feature, 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
 
                 seqid = feature["seqid"]
@@ -5169,7 +5178,8 @@ class Alignment(_Annotatable, AlignmentI, SequenceCollection):
         # there's no sequence to bind to, the feature is directly on self
         # todo gah check handling of strand etc..., maybe reuse code
         # in Sequence?
-        feature["map"] = Map(parent_length=len(self), locations=feature.pop("spans"))
-        feature.pop("reversed", None)
+        fmap = Map(parent_length=len(self), locations=feature.pop("spans"))
+        if feature.pop("reversed", None):
+            fmap = fmap.nucleic_reversed()
         feature.pop("strand", None)
-        return Annotation(parent=self, **feature)
+        return Annotation(parent=self, map=fmap, **feature)
diff --git a/src/cogent3/core/annotation.py b/src/cogent3/core/annotation.py
index 45f56785f..46134562e 100644
--- a/src/cogent3/core/annotation.py
+++ b/src/cogent3/core/annotation.py
@@ -506,7 +506,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.seqid}" {self.biotype}{name} at {self.map}'
 
     def _projected_to_base(self, base):
         if self.parent == base:
@@ -557,6 +557,11 @@ class Annotation(_AnnotationCore, _Serialisable):
         name = ", ".join(feat_names)
         map = Map(spans=combined, parent_length=len(self.parent))
         map = map.covered()  # No overlaps
+        # the covered method drops reversed status so we need to
+        # resurrect that, but noting we've not checked consistency
+        # across the features
+        if self.map.reverse != map.reverse:
+            map = map.reversed()
         seqid = ", ".join(seqids) if seqids else None
         biotype = ", ".join(biotypes)
         return self.__class__(
diff --git a/src/cogent3/core/location.py b/src/cogent3/core/location.py
index 928bf6be3..6621ddcc1 100644
--- a/src/cogent3/core/location.py
+++ b/src/cogent3/core/location.py
@@ -1,4 +1,3 @@
-#!/usr/bin/env python
 """Alignments and Sequences are _Annotatables
 _Annotatables hold a list of Maps.
 Maps can be Features, Variables or AlignedSequences.
@@ -42,6 +41,9 @@ import copy
 from bisect import bisect_left, bisect_right
 from functools import total_ordering
 from itertools import chain
+from typing import Union
+
+from numpy import array, ndarray
 
 from cogent3.util.misc import (
     ClassChecker,
@@ -853,6 +855,51 @@ class Map(object):
 
         return zeroed
 
+    T = Union[ndarray, int]
+
+    def absolute_position(self, rel_pos: T) -> T:
+        """converts rel_pos into an absolute position
+
+        Raises
+        ------
+        raises ValueError if rel_pos < 0
+        """
+        check = array([rel_pos], dtype=int) if isinstance(rel_pos, int) else rel_pos
+        if check.min() < 0:
+            raise ValueError(f"must positive, not {rel_pos=}")
+
+        if len(self) == self.parent_length:
+            # handle case of reversed here?
+            return rel_pos
+
+        return self.start + rel_pos
+
+    def relative_position(self, abs_pos: T) -> T:
+        """converts abs_pos into an relative position
+
+        Raises
+        ------
+        raises ValueError if abs_pos < 0
+        """
+        check = array([abs_pos], dtype=int) if isinstance(abs_pos, int) else abs_pos
+        if check.min() < 0:
+            raise ValueError(f"must positive, not {abs_pos=}")
+        return abs_pos - self.start
+
+    def convert_position_for_reversed(self, rel_pos: T) -> T:
+        """adjusts a relative position for changed orientation
+
+        Raises
+        ------
+        raises ValueError if rel_pos < 0 or rel_pos > len(self)
+        """
+        check = array([rel_pos], dtype=int) if isinstance(rel_pos, int) else rel_pos
+        if check.min() < 0 or check.max() > len(self):
+            raise ValueError(
+                f"must be a relative position within [0,{len(self)}], not {rel_pos=}"
+            )
+        return len(self) - rel_pos - 1
+
 
 class SpansOnly(ConstrainedList):
     """List that converts elements to Spans on addition."""
diff --git a/src/cogent3/core/sequence.py b/src/cogent3/core/sequence.py
index df7db23a7..535cd3e84 100644
--- a/src/cogent3/core/sequence.py
+++ b/src/cogent3/core/sequence.py
@@ -934,12 +934,26 @@ class Sequence(_Annotatable, SequenceI):
         # i,e., the stop can be equal to the length of the view
         query_end = self._seq.absolute_position(stop, include_boundary=True)
 
-        # 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)
-        )
+        reversed = query_end < query_start
+        if reversed:
+            query_start, query_end = query_end, query_start
+
+        query_start = max(query_start, 0)
+        # in the underlying db, features are always plus strand oriented
+        # (need to check that for user defined features)
+        # a '-' strand feature in the db may be [(2, 4), (6, 8)]
+        # so we would take 7-6, 3-2 and complement it
+        # if the current sequence is reverse complemented
+        # then a '+' feature from the db needs to be nucleic reversed
+        # and a '-' feature from the db needs to be nucleic reversed too
+
+        # making this more complicated is self.make_feature() assumes
+        # all coordinates are with respect to the current orientation
+        # so self.make_feature(data(spans=[(2,4)])) has different meaning if
+        # self is rc'ed, it would correspond to len(self)-2, etc...
+        # To piggy-back on that method we need to convert our feature spans
+        # into the current orientation. HOWEVER, we also have the reversed
+        # flag which comes back from the db
 
         for feature in self.annotation_db.get_features_matching(
             seqid=self.name,
@@ -952,9 +966,18 @@ 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"])
-
+            spans = array(feature.pop("spans"), dtype=int)
+            for i, v in enumerate(spans.ravel()):
+                rel_pos = self._seq.relative_position(v)
+                if rel_pos < 0 or rel_pos > len(self):
+                    raise NotImplementedError
+                spans.ravel()[i] = rel_pos
+
+            if reversed:
+                # see above comment
+                spans = len(self) - spans
+
+            feature["spans"] = spans.tolist()
             yield self.make_feature(feature)
 
     def _relative_spans(self, spans):


Test patch:
diff --git a/tests/test_core/test_alignment.py b/tests/test_core/test_alignment.py
index 207f3bb4f..7a2ec4ff5 100644
--- a/tests/test_core/test_alignment.py
+++ b/tests/test_core/test_alignment.py
@@ -3312,21 +3312,21 @@ def test_get_annotations_from_any_seq(cls):
     seqs.get_seq("seq3").add_feature(biotype="exon", name="annotation3", spans=[(3, 6)])
     got = list(seqs.get_features())
     assert len(got) == 3
-    assert str(got[0]) == 'exon "annotation1" at [3:8]/9'
-    assert str(got[1]) == 'exon "annotation2" at [1:2]/9'
-    assert str(got[2]) == 'exon "annotation3" at [3:6]/9'
+    assert 'exon "annotation1" at [3:8]/9' in str(got[0])
+    assert 'exon "annotation2" at [1:2]/9' in str(got[1])
+    assert 'exon "annotation3" at [3:6]/9' in str(got[2])
 
     got = list(seqs.get_features(name="annotation1"))
     assert len(got) == 1
-    assert str(got[0]) == 'exon "annotation1" at [3:8]/9'
+    assert 'exon "annotation1" at [3:8]/9' in str(got[0])
 
     got = list(seqs.get_features(biotype="exon", name="annotation2"))
     assert len(got) == 1
-    assert str(got[0]) == 'exon "annotation2" at [1:2]/9'
+    assert 'exon "annotation2" at [1:2]/9' in str(got[0])
 
     got = list(seqs.get_features(name="annotation3"))
     assert len(got) == 1
-    assert str(got[0]) == 'exon "annotation3" at [3:6]/9'
+    assert 'exon "annotation3" at [3:6]/9' in str(got[0])
 
 
 def test_annotate_matches_to():
diff --git a/tests/test_core/test_annotation.py b/tests/test_core/test_annotation.py
index 2dc10a479..1a0278de8 100644
--- a/tests/test_core/test_annotation.py
+++ b/tests/test_core/test_annotation.py
@@ -105,37 +105,7 @@ class TestAnnotations(unittest.TestCase):
             expected = seq_expecteds[annot_type]
             assert observed == expected
 
-    @pytest.mark.xfail(reason="todo gah get_region_covering_all not working")
-    def test_slice_aln_with_annotations(self):
-        """test that annotations of sequences and alignments survive alignment
-        slicing."""
-        aln_expecteds = {
-            "misc_feature": {"FAKE01": "TTTGGGGGGGGGG", "FAKE02": "TTTGGGGGGGGGG"},
-            "CDS": {"FAKE01": "GGGGGGGGGG", "FAKE02": "GGGGGGGGGG"},
-            "5'UTR": {"FAKE01": "CC", "FAKE02": "CC"},
-            "LTR": {"FAKE01": "CCCTTTTT", "FAKE02": "CCCTTTTT"},
-        }
-        newaln = self.aln[10:]
-        for annot_type in ["LTR", "misc_feature", "CDS", "5'UTR"]:
-            feature_list = newaln.get_features(biotype=annot_type)
-            new = newaln.get_region_covering_all(feature_list).get_slice().to_dict()
-            expected = aln_expecteds[annot_type]
-            assert expected == new, (annot_type, expected, new)
-            if annot_type in ["misc_feature", "LTR"]:
-                continue  # because seqs haven't been annotated with it
-            for name in self.aln.names:
-                orig = str(
-                    list(self.aln.get_features(seqid=name, biotype=annot_type))[
-                        0
-                    ].get_slice()
-                )
-                new = str(
-                    list(newaln.get_features(seqid=name, biotype=annot_type))[
-                        0
-                    ].get_slice()
-                )
-                assert orig == new, (name, annot_type, orig, new)
-
+    @pytest.mark.xfail(reason="todo gah update to latest API")
     def test_feature_projection(self):
         expecteds = {"FAKE01": "CCCAAAATTTTTT", "FAKE02": "CCC-----TTTTT"}
         aln_ltr = self.aln.get_features_matching("LTR")[0]
@@ -162,6 +132,7 @@ 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 = {
@@ -252,20 +223,28 @@ def test_constructing_collections(alignment):
     assert coll.annotation_db.num_matches() == expect
 
 
-@pytest.mark.xfail(reason="todo gah get_region_covering_all not working")
-def test_region_covering_all():
+@pytest.mark.parametrize("reversed", (False, True))
+@pytest.mark.parametrize("annot_type", ("LTR", "misc_feature", "CDS", "5'UTR"))
+def test_region_union_on_alignment(annot_type, reversed):
+    # >FAKE01
+    # AACCCAAAATTTTTTGGGGGGGGGGCCCC
+    # >FAKE02
+    # AACCC-----TTTTTGGGGGGGGGGCC--
+    #   **  5'UTR
+    #   ************* LTR
+    #                **********  CDS
+    #             *************  misc_feature
+
     aln = makeSampleAlignment()
     aln_expecteds = {
         "misc_feature": {"FAKE01": "TTTGGGGGGGGGG", "FAKE02": "TTTGGGGGGGGGG"},
         "CDS": {"FAKE01": "GGGGGGGGGG", "FAKE02": "GGGGGGGGGG"},
         "5'UTR": {"FAKE01": "CC", "FAKE02": "CC"},
-        "LTR": {"FAKE01": "CCCTTTTT", "FAKE02": "CCCTTTTT"},
+        "LTR": {"FAKE01": "CCCAAAATTTTTT", "FAKE02": "CCC-----TTTTT"},
     }
-    newaln = aln[10:]
-    for annot_type in ["LTR", "misc_feature", "CDS", "5'UTR"]:
-        feature_list = list(newaln.get_features(biotype=annot_type))
-        new = newaln.get_region_covering_all(feature_list).get_slice().to_dict()
-        expected = aln_expecteds[annot_type]
-        assert expected == new, (annot_type, expected, new)
-        if annot_type in ["misc_feature", "LTR"]:
-            continue  # because seqs haven't been annotated with it
+    newaln = aln.rc() if reversed else aln
+    feature_list = list(newaln.get_features(biotype=annot_type, on_alignment=True))
+    new = feature_list[0].union(feature_list[1:])
+    new = new.get_slice().to_dict()
+    expected = aln_expecteds[annot_type]
+    assert expected == new, (annot_type, expected, new)
diff --git a/tests/test_core/test_annotation_db.py b/tests/test_core/test_annotation_db.py
index a143ffd75..99aa4e696 100644
--- a/tests/test_core/test_annotation_db.py
+++ b/tests/test_core/test_annotation_db.py
@@ -262,9 +262,9 @@ def test_gav():
         on_alignment=False,
     )
     seq.annotation_db = db
-    plus = list(seq.get_features_matching(name="plus"))[0]
+    plus = list(seq.get_features(name="plus"))[0]
     assert str(plus.get_slice()) == plus_seq
-    minus = list(seq.get_features_matching(name="minus"))[0]
+    minus = list(seq.get_features(name="minus"))[0]
     assert str(minus.get_slice()) == minus_seq
 
     f = Feature(seq, "cds", "minus", reversed([reversed(s) for s in minus_spans]))
@@ -278,9 +278,9 @@ def test_gav():
     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]
+    plus = list(rced.get_features(name="plus"))[0]
     assert str(plus.get_slice()) == plus_seq
-    minus = list(rced.get_features_matching(name="minus"))[0]
+    minus = list(rced.get_features(name="minus"))[0]
     print(rf, minus)
     assert str(minus.get_slice()) == minus_seq
 
@@ -394,9 +394,9 @@ def test_rc_get_slice_negative_feature(seq_db):
     the same sequence before and after the sequence is reverse complemented
     """
 
-    feat = list(seq_db.get_features_matching(name="Transcript:B0019.1"))[0]
+    feat = list(seq_db.get_features(name="Transcript:B0019.1"))[0]
     rc_seq = seq_db.rc()
-    r_feat = list(rc_seq.get_features_matching(name="Transcript:B0019.1"))[0]
+    r_feat = list(rc_seq.get_features(name="Transcript:B0019.1"))[0]
 
     assert feat.get_slice() == r_feat.get_slice()
 
@@ -479,6 +479,7 @@ 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
@@ -491,26 +492,19 @@ def test_sequence_collection_annotate_from_gff():
     # the seq for which the seqid was provided is annotated
     seq = seq_collection.get_seq("test_seq")
     assert seq_collection.get_seq("test_seq").annotation_db is not None
-    assert (
-        len(
-            list(
-                seq_collection.get_seq("test_seq").get_features_matching(
-                    allow_partial=True
-                )
-            )
-        )
-        == 5
-    )
+    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
     # the annotation_db on the seq and the seq collection are the same object
     assert (
         seq_collection.get_seq("test_seq").annotation_db is seq_collection.annotation_db
     )
-    got = list(seq.get_features_matching(feature_type="CDS"))
+    got = list(seq.get_features(feature_type="CDS"))
     assert len(got) == 2
 
-    got = list(seq.get_features_matching(feature_type="CpG"))
+    got = list(seq.get_features(feature_type="CpG"))
     assert len(got) == 1
 
 
diff --git a/tests/test_core/test_features.py b/tests/test_core/test_features.py
index ea22bf34e..4f6543466 100755
--- a/tests/test_core/test_features.py
+++ b/tests/test_core/test_features.py
@@ -4,7 +4,6 @@ import pytest
 
 from cogent3 import ASCII, DNA, make_aligned_seqs
 from cogent3.core.annotation import Feature, Variable
-
 # Complete version of manipulating sequence annotations
 from cogent3.util.deserialise import deserialise_object
 
@@ -51,7 +50,8 @@ class FeaturesTest(TestCase):
 
         exons = list(self.s.get_features(biotype="exon"))
         self.assertEqual(
-            str(exons), '[exon "fred" at [10:15]/48, exon "trev" at [30:40]/48]'
+            str(exons),
+            '["Orig" exon "fred" at [10:15]/48, "Orig" exon "trev" at [30:40]/48]',
         )
 
     @pytest.mark.xfail(reason="todo gah update test to use latest API")
diff --git a/tests/test_core/test_location.py b/tests/test_core/test_location.py
index 08416e30f..ae5fdb846 100644
--- a/tests/test_core/test_location.py
+++ b/tests/test_core/test_location.py
@@ -1,5 +1,3 @@
-#!/usr/bin/env python
-
 """Unit tests for Range, Span and Point classes.
 """
 from unittest import TestCase, main
@@ -566,6 +564,34 @@ class MapTests(TestCase):
             got = gap_coords_to_map({20: 1}, len(seq))
 
 
-# run the following if invoked from command-line
-if __name__ == "__main__":
-    main()
+def test_map_plus_position():
+    # seq is 9 long
+    # plus coords  012345678
+    # +slice         **
+    # plus seq     AAACCCTGG
+
+    # orig = Aligned(*DNA.make_seq("AAACCCTGG", name="a").parse_out_gaps())
+    orig = Map([(0, 9)], parent_length=9)
+    assert orig.absolute_position(2) == 2
+    assert orig.absolute_position(6) == 6
+
+    assert orig.relative_position(2) == 2
+    assert orig.relative_position(6) == 6
+
+    subseq = orig[2:4]
+    assert subseq.absolute_position(0) == 2
+    assert subseq.absolute_position(4) == 6
+
+    assert subseq.relative_position(2) == 0
+    assert subseq.relative_position(6) == 4
+
+    # minus coords 012345678
+    # rel coord      01234
+    # -slice         *****
+    # minus seq    CCAGGGTTT
+    # plus coords  876543210
+    rc = orig.nucleic_reversed()
+    rcsubseq = rc[2:7]
+    abs_coord = rcsubseq.absolute_position(0)
+    abs_rc_coord = orig.convert_position_for_reversed(abs_coord)
+    assert abs_rc_coord == 6, (abs_coord, abs_rc_coord)