Repository: cogent3/cogent3
Pull request #1711 (2024-02-12 05:55:50)

Steps towards simplifying Map


Patch:
diff --git a/src/cogent3/core/alignment.py b/src/cogent3/core/alignment.py
index 4cabd43a8..0e8663081 100644
--- a/src/cogent3/core/alignment.py
+++ b/src/cogent3/core/alignment.py
@@ -2184,14 +2184,9 @@ class Aligned:
         -------
         a copy of self
         """
-        new_seq = self.data.copy(exclude_annotations=exclude_annotations)
-        db = new_seq.annotation_db
+        new_seq = self.data.copy(exclude_annotations=exclude_annotations, sliced=sliced)
         if sliced:
-            span = self.map.get_covering_span()
-            new_seq = type(new_seq)(
-                str(new_seq[span.start : span.end]), info=new_seq.info, name=self.name
-            )
-            new_seq.annotation_offset = self.map.start
+            db = new_seq.annotation_db
             if self.map.reverse or exclude_annotations:
                 new_seq.annotation_db = None
             else:
@@ -2252,13 +2247,30 @@ class Aligned:
 
     def __getitem__(self, slice):
         new_map = self.map[slice]
+        data = (
+            self.data[new_map.start : new_map.end] if new_map.useful else self.data[0:0]
+        )
         if new_map.reverse:
             # A reverse slice means we should have an empty sequence
             new_map = type(new_map)(locations=(), parent_length=len(self.data))
-        return Aligned(new_map, self.data)
+        elif new_map.useful:
+            new_map = new_map.zeroed()
+        return Aligned(new_map, data)
 
     def rc(self):
-        return Aligned(self.map.reversed(), self.data)
+        # the map now no longer stores history of reversing,
+        # as this done by Sequence & SeqView
+        new_map = self.map.nucleic_reversed()
+        # in this new approach, the map is always plus strand, so
+        # the following comprehension ensures the order of spans
+        # and the span start/end satisfy this
+        spans = [s if s.lost else s.reversed() for s in new_map.spans]
+        spans.reverse()
+        new_map = type(new_map)(
+            spans=spans,
+            parent_length=new_map.parent_length,
+        )
+        return Aligned(new_map, self.data.rc())
 
     def to_rna(self):
         return Aligned(self.map, self.data.to_rna())
@@ -5215,10 +5227,11 @@ class Alignment(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?
+        revd = feature.pop("strand", None) == "-"
+        feature["strand"] = "-" if revd else "+"
         fmap = Map(parent_length=len(self), locations=feature.pop("spans"))
-        if feature.pop("reversed", None):
+        if revd:
             fmap = fmap.nucleic_reversed()
-        feature.pop("strand", None)
         return Feature(parent=self, map=fmap, **feature)
 
     def _get_seq_features(
@@ -5260,12 +5273,12 @@ class Alignment(AlignmentI, SequenceCollection):
 
         for seqid in seqids:
             seq = self.named_seqs[seqid]
-
-            start, end = seq.map.start, seq.map.end
+            # we use parent seqid, stored on SeqView
+            parent_id, start, end, _ = seq.data.parent_coordinates()
             offset = seq.data.annotation_offset
 
             for feature in self.annotation_db.get_features_matching(
-                seqid=seqid,
+                seqid=parent_id,
                 biotype=biotype,
                 name=name,
                 on_alignment=False,
@@ -5341,12 +5354,13 @@ class Alignment(AlignmentI, SequenceCollection):
                 raise RuntimeError(f"{on_alignment=} {feature=}")
             if seq_map is None:
                 seq_map = self.seqs[0].map
+                *_, strand = self.seqs[0].data.parent_coordinates()
 
             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
+            feature["strand"] = "-" if strand == -1 else "+"
             yield self.make_feature(feature=feature, on_alignment=on_al)
 
 
diff --git a/src/cogent3/core/annotation.py b/src/cogent3/core/annotation.py
index 7c7e80b40..5f8a4270b 100644
--- a/src/cogent3/core/annotation.py
+++ b/src/cogent3/core/annotation.py
@@ -18,10 +18,13 @@ class Feature:
         "_name",
         "_serialisable",
         "_id",
+        "_strand",
     )
 
     # todo gah implement a __new__ to trap args for serialisation purposes?
-    def __init__(self, *, parent, seqid: str, map: Map, biotype: str, name: str):
+    def __init__(
+        self, *, parent, seqid: str, map: Map, biotype: str, name: str, strand: str
+    ):
         # _serialisable is used for creating derivative instances
         d = locals()
         exclude = ("self", "__class__", "kw")
@@ -35,6 +38,7 @@ class Feature:
         data = [id(self.parent), tuple(self.map.get_coordinates())]
         data.extend((self.seqid, self.biotype, self.name))
         self._id = hash(tuple(data))
+        self._strand = strand
 
     def __eq__(self, other):
         return self._id == other._id
@@ -89,8 +93,19 @@ class Feature:
         if not (complete or map.complete):
             map = map.without_gaps()
         if not allow_gaps:
-            return self.parent[map]
-        return self.parent[map.start : map.end]
+            if self.reversed:
+                map = map.reversed()
+            result = self.parent[map]
+            if self.reversed:
+                result = result.rc()
+            return result
+
+        # all slicing now requires start < end
+        start, end = min(map.start, map.end), max(map.start, map.end)
+        result = self.parent[start:end]
+        if self.reversed:
+            result = result.rc()
+        return result
 
     def without_lost_spans(self):
         """Keeps only the parts which are actually present in the underlying sequence"""
@@ -242,9 +257,13 @@ class Feature:
             **self._serialisable,
             **dict(
                 spans=self.map.get_coordinates(),
-                strand="-" if self.map.reverse else "+",
             ),
         }
         for key in ("map", "parent"):
             result.pop(key, None)
         return result
+
+    @property
+    def reversed(self):
+        """whether Feature is on the reverse strand relative to bound object"""
+        return self._strand == "-"
diff --git a/src/cogent3/core/annotation_db.py b/src/cogent3/core/annotation_db.py
index b9013f9b2..9e05b2e7f 100644
--- a/src/cogent3/core/annotation_db.py
+++ b/src/cogent3/core/annotation_db.py
@@ -66,7 +66,7 @@ class FeatureDataType(typing.TypedDict):
     biotype: str  # rename type attr of cogent3 Annotatables to match this?
     name: str  # rename to name to match cogent3 Annotatable.name?
     spans: list[tuple[int, int]]
-    reversed: bool  # True if feature on reverse strand
+    strand: str  # "-" if feature on reverse strand
     on_alignment: bool  # True if feature on an alignment
 
 
@@ -669,7 +669,6 @@ class SqliteAnnotationDbMixin:
             result = dict(zip(result.keys(), result))
             result["on_alignment"] = result.get("on_alignment")
             result["spans"] = [tuple(c) for c in result["spans"]]
-            result["reversed"] = result.pop("strand", None) == "-"
             yield result
 
     def get_feature_children(
@@ -786,7 +785,6 @@ class SqliteAnnotationDbMixin:
                 result = dict(zip(result.keys(), result))
                 result["on_alignment"] = result.get("on_alignment")
                 result["spans"] = [tuple(c) for c in result["spans"]]
-                result["reversed"] = result.pop("strand", None) == "-"
                 yield result
 
     def num_matches(
diff --git a/src/cogent3/core/location.py b/src/cogent3/core/location.py
index 8aeef14ec..db6c9d67c 100644
--- a/src/cogent3/core/location.py
+++ b/src/cogent3/core/location.py
@@ -823,12 +823,17 @@ class Map:
         data = self.to_rich_dict()
         zeroed = deserialise_map_spans(data)
         zeroed.parent_length = len(zeroed.get_covering_span())
-        min_val = min(zeroed.start, zeroed.end)
+        shift = min(zeroed.start, zeroed.end)
+        new_end = 0
         for span in zeroed.spans:
             if span.lost:
                 continue
-            span.start -= min_val
-            span.end -= min_val
+            span.start -= shift
+            span.end -= shift
+            new_end = max(new_end, span.end)
+
+        zeroed.start = 0
+        zeroed.end = new_end
 
         return zeroed
 
diff --git a/src/cogent3/core/sequence.py b/src/cogent3/core/sequence.py
index 4e0f7d1f7..d72cfa81f 100644
--- a/src/cogent3/core/sequence.py
+++ b/src/cogent3/core/sequence.py
@@ -1030,11 +1030,9 @@ class Sequence(SequenceI):
         """
         feature = dict(feature)
         seq_rced = self._seq.reverse
-        # 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?
-        revd = feature.pop("reversed", None) or feature.pop("strand", None) == "-"
         spans = feature.pop("spans", None)
+        revd = feature.pop("strand", None) == "-"
+        feature["strand"] = "+" if revd == seq_rced else "-"
 
         vals = array(spans)
         pre = abs(vals.min()) if vals.min() < 0 else 0
@@ -1264,6 +1262,7 @@ class Sequence(SequenceI):
                 name=None,
                 biotype=None,
                 map=Map(locations=[], parent_length=len(self)),
+                strand="+",
             )
         else:
             region = annotations[0].union(annotations[1:])
@@ -1630,21 +1629,22 @@ class Sequence(SequenceI):
         drawer.layout.update(xaxis=xaxis, yaxis=yaxis)
         return drawer
 
-    def parent_coordinates(self) -> Tuple[int, int, int]:
-        """returns start, stop, strand of this sequence on its parent
+    def parent_coordinates(self) -> Tuple[str, int, int, int]:
+        """returns seqid, start, stop, strand of this sequence on its parent
 
         Notes
         -----
-        Returned coordinates are with respect to the plus strand, irrespective
-        of whether the sequence has been reversed complemented or not.
+        seqid is the identifier of the parent. Returned coordinates are with
+        respect to the plus strand, irrespective of whether the sequence has
+        been reversed complemented or not.
 
         Returns
         -------
-        start, end, strand of this sequence on the parent. strand is either
+        seqid, start, end, strand of this sequence on the parent. strand is either
         -1 or 1.
         """
         strand = -1 if self._seq.reverse else 1
-        return self._seq.parent_start, self._seq.parent_stop, strand
+        return self._seq.seqid, self._seq.parent_start, self._seq.parent_stop, strand
 
 
 class ProteinSequence(Sequence):


Test patch:
diff --git a/tests/test_core/test_alignment.py b/tests/test_core/test_alignment.py
index 70a331630..d5614bdef 100644
--- a/tests/test_core/test_alignment.py
+++ b/tests/test_core/test_alignment.py
@@ -238,24 +238,6 @@ class SequenceCollectionBaseTests(object):
         self.assertEqual(list(self.ordered1.iter_selected()), ["A"] * 5 + ["B"] * 5)
         self.assertEqual(list(self.ordered2.iter_selected()), ["B"] * 5 + ["A"] * 5)
 
-    def test_iter_selected(self):
-        """SequenceCollection iter_selected() should iterate over items in correct order"""
-        # should work if one row
-        self.assertEqual(list(self.one_seq.iter_selected()), ["A"] * 5)
-        # should take order into account
-        self.assertEqual(list(self.ordered1.iter_selected()), ["A"] * 5 + ["B"] * 5)
-        self.assertEqual(list(self.ordered2.iter_selected()), ["B"] * 5 + ["A"] * 5)
-        # should allow row and/or col specification
-        r = self.ragged_padded
-        self.assertEqual(
-            list(r.iter_selected(seq_order=["c", "b"], pos_order=[5, 1, 3])),
-            list("-AA-A-"),
-        )
-        # should not interfere with superclass iteritems()
-        i = list(r.named_seqs.items())
-        i.sort()
-        self.assertEqual(i, [("a", "AAAAAA"), ("b", "AAA---"), ("c", "AAAA--")])
-
     def test_take_seqs(self):
         """SequenceCollection take_seqs should return new SequenceCollection with selected seqs."""
         a = self.ragged_padded.take_seqs(list("bc"))
@@ -1060,20 +1042,6 @@ class AlignmentBaseTests(SequenceCollectionBaseTests):
             list(map(list, ["AAA", "AAA", "AAA", "A-A", "A--", "A--"])),
         )
 
-    def test_take_positions(self):
-        """SequenceCollection take_positions should return new alignment w/ specified pos"""
-        self.assertEqual(
-            self.gaps.take_positions([5, 4, 0]), {"a": "AAA", "b": "A-A", "c": "--A"}
-        )
-        self.assertTrue(
-            isinstance(self.gaps.take_positions([0]), _SequenceCollectionBase)
-        )
-        # should be able to negate
-        self.assertEqual(
-            self.gaps.take_positions([5, 4, 0], negate=True),
-            {"a": "AAAA", "b": "--AA", "c": "A---"},
-        )
-
     def test_take_positions_info(self):
         aln = self.Class(
             {"a": "AAAAAAA", "b": "A--A-AA", "c": "AA-----"}, info={"key": "value"}
@@ -1924,49 +1892,6 @@ class ArrayAlignmentTests(AlignmentBaseTests, TestCase):
 class AlignmentTests(AlignmentBaseTests, TestCase):
     Class = Alignment
 
-    def test_sliced_deepcopy(self):
-        """correctly deep copy aligned objects in an alignment"""
-
-        def eval_data(data):
-            orig = self.Class(data)
-            aln = orig[2:5]
-
-            notsliced = aln.deepcopy(sliced=False)
-            sliced = aln.deepcopy(sliced=True)
-            for name in orig.names:
-                # if not sliced, underlying seq data should be same length
-                # as original
-                self.assertEqual(
-                    len(notsliced.named_seqs[name].data),
-                    len(orig.named_seqs[name].data),
-                )
-                # and the map.parent_length attributes should match
-                self.assertEqual(
-                    notsliced.named_seqs[name].map.parent_length,
-                    orig.named_seqs[name].map.parent_length,
-                )
-                # and map.parent_length and len(data) should match
-                self.assertEqual(
-                    sliced.named_seqs[name].map.parent_length,
-                    len(sliced.named_seqs[name].data),
-                )
-
-                # if sliced, seq data should be < orig
-                self.assertLess(
-                    len(sliced.named_seqs[name].data), len(orig.named_seqs[name].data)
-                )
-                # and map.parent_length and len(data) should match
-                self.assertEqual(
-                    sliced.named_seqs[name].map.parent_length,
-                    len(sliced.named_seqs[name].data),
-                )
-
-        eval_data({"seq1": "ACAACGACG", "seq2": "ACGACGACG"})
-        eval_data({"seq1": "-CAACGACG", "seq2": "ACGACGACG"})
-        eval_data({"seq1": "ACAACGAC-", "seq2": "ACGACGACG"})
-        eval_data({"seq1": "AC-ACGACG", "seq2": "ACGACGACG"})
-        eval_data({"seq1": "ACA-CGAC-", "seq2": "ACGACGACG"})
-
     def test_sliding_windows(self):
         """sliding_windows should return slices of alignments."""
         alignment = self.Class(
@@ -2585,6 +2510,13 @@ def test_copy_annotations_incompat_type_fails(seqcoll_db, gb_db):
         seqcoll_db.copy_annotations({"a": "ACGGT"})
 
 
+def test_aligned_deepcopy_sliced():
+    a = Aligned(*DNA.make_seq("GCAAGGCGCCAA").parse_out_gaps())
+    sliced = a[2:3]
+    sl_cp = sliced.deepcopy(sliced=True)
+    assert str(sl_cp) == str(sliced)
+
+
 def test_deepcopy_with_features(DATA_DIR):
     """correctly deepcopy Aligned objects in an alignment"""
     path = DATA_DIR / "brca1_5.paml"
@@ -3533,6 +3465,120 @@ def test_get_degapped_relative_to(cls):
         aln.get_degapped_relative_to("nameX")
 
 
+@pytest.mark.parametrize(
+    "data",
+    (
+        {"seq1": "ACAACGACG", "seq2": "ACGACGACG"},
+        {"seq1": "-CAACGACG", "seq2": "ACGACGACG"},
+        {"seq1": "ACAACGAC-", "seq2": "ACGACGACG"},
+        {"seq1": "AC-ACGACG", "seq2": "ACGACGACG"},
+        {"seq1": "ACA-CGAC-", "seq2": "ACGACGACG"},
+    ),
+)
+@pytest.mark.parametrize("name", ("seq1", "seq2"))
+@pytest.mark.parametrize("rev", (False, True))
+def test_sliced_deepcopy(data, name, rev):
+    """correctly deep copy aligned objects in an alignment"""
+
+    orig = Alignment(data, moltype="dna")
+    slice_start, slice_end = 3, 5
+    aln = orig[slice_start:slice_end]
+    if rev:
+        aln = aln.rc()
+
+    # the map just mirrors the slice effect, it no longer has "memory" of rev complement
+    assert aln.named_seqs[name].map.reverse == False
+
+    notsliced = aln.deepcopy(sliced=False)
+    # the annotation offsets should match original object
+    assert {s.data.annotation_offset for s in notsliced.seqs} == {
+        s.data.annotation_offset for s in aln.seqs
+    }
+    sliced = aln.deepcopy(sliced=True)
+    assert sliced.to_dict() == notsliced.to_dict()
+    sliced_seq = sliced.named_seqs[name]
+    notsliced_seq = notsliced.named_seqs[name]
+    assert str(sliced_seq) == str(notsliced_seq)
+
+    # the annotation offsets should match original object
+    assert {s.data.annotation_offset for s in sliced.seqs} == {
+        s.data.annotation_offset for s in aln.seqs
+    }
+    # if not sliced in copy, underlying seq data is identical to original
+    assert (
+        notsliced.named_seqs[name].data._seq.seq is orig.named_seqs[name].data._seq.seq
+    )
+    # but not the same for sliced
+    assert (
+        sliced.named_seqs[name].data._seq.seq is not orig.named_seqs[name].data._seq.seq
+    )
+
+    # the map just mirrors the slice effect, it has no  "memory" of rev complement
+    assert notsliced.named_seqs[name].map.reverse == False
+    assert sliced.named_seqs[name].map.reverse == False
+
+    assert sliced.named_seqs[name].map.parent_length == len(
+        str(sliced_seq).replace("-", "")
+    )
+    assert notsliced.named_seqs[name].map.parent_length == len(
+        str(notsliced_seq).replace("-", "")
+    )
+    # and map.parent_length and len(data) should match
+    assert sliced.named_seqs[name].map.parent_length == len(
+        sliced.named_seqs[name].data
+    )
+
+    # if sliced, seq data should be < orig
+    assert len(sliced.named_seqs[name].data) < len(orig.named_seqs[name].data)
+
+
+def test_aligned_deepcopy_sliced_map_matches_data():
+    m, seq = DNA.make_seq("ACAACGACG", name="seq1").parse_out_gaps()
+    aligned = Aligned(m, seq)
+    sliced = aligned[3:5]
+    sliced_copy = sliced.deepcopy(sliced=True)
+    assert sliced_copy.map.parent_length == len(sliced_copy.data)
+
+
+@pytest.mark.parametrize("cls", (SequenceCollection, Alignment, ArrayAlignment))
+def test_iter_selected(cls):
+    """SequenceCollection iter_selected() should iterate over items in correct order"""
+    # should work if one row
+    one_seq = cls({"a": "AAAAA"})
+    ragged_padded = cls({"a": "AAAAAA", "b": "AAA---", "c": "AAAA--"})
+    ordered1 = cls({"a": "AAAAA", "b": "BBBBB"}, names=["a", "b"])
+    ordered2 = cls({"a": "AAAAA", "b": "BBBBB"}, names=["b", "a"])
+
+    assert list(one_seq.iter_selected()) == (["A"] * 5)
+    # should take order into account
+    assert list(ordered1.iter_selected()) == (["A"] * 5 + ["B"] * 5)
+    assert list(ordered2.iter_selected()) == (["B"] * 5 + ["A"] * 5)
+    # should allow row and/or col specification
+    r = ragged_padded
+    assert list(r.iter_selected(seq_order=["c", "b"], pos_order=[5, 1, 3])) == list(
+        "-AA-A-"
+    )
+    # should not interfere with superclass iteritems()
+    i = list(r.named_seqs.items())
+    i.sort()
+    assert i == ([("a", "AAAAAA"), ("b", "AAA---"), ("c", "AAAA--")])
+
+
+@pytest.mark.parametrize("cls", (Alignment, ArrayAlignment))
+def test_take_positions(cls):
+    """SequenceCollection take_positions should return new alignment w/ specified pos"""
+    gaps = cls({"a": "AAAAAAA", "b": "A--A-AA", "c": "AA-----"})
+    assert gaps.take_positions([5, 4, 0]) == {"a": "AAA", "b": "A-A", "c": "--A"}
+    assert isinstance(gaps.take_positions([0]), _SequenceCollectionBase)
+
+    # should be able to negate
+    assert gaps.take_positions([5, 4, 0], negate=True) == {
+        "a": "AAAA",
+        "b": "--AA",
+        "c": "A---",
+    }
+
+
 @pytest.mark.parametrize("array_align", (True, False))
 @pytest.mark.parametrize(
     "seq_moltype", (("ACCT--GT", "dna"), ("ACCT--GU", "rna"), ("MTST--T", "protein"))
diff --git a/tests/test_core/test_annotation_db.py b/tests/test_core/test_annotation_db.py
index 2eac5427a..28e45da99 100644
--- a/tests/test_core/test_annotation_db.py
+++ b/tests/test_core/test_annotation_db.py
@@ -568,7 +568,7 @@ def test_get_slice():
         "spans": [
             (5, 10),
         ],
-        "reversed": False,
+        "strand": "+",
     }
     feature = seq.make_feature(feature_data)
     got = feature.get_slice()
diff --git a/tests/test_core/test_features.py b/tests/test_core/test_features.py
index 5330b2746..d51dea406 100755
--- a/tests/test_core/test_features.py
+++ b/tests/test_core/test_features.py
@@ -693,7 +693,12 @@ def test_feature_not_equal(ann_seq):
 def test_feature_not_equal_attr(ann_seq, attr):
     (f1,) = list(ann_seq.get_features(biotype="gene"))
     attrs = dict(
-        parent=f1.parent, seqid=f1.seqid, biotype=f1.biotype, map=f1.map, name=f1.name
+        parent=f1.parent,
+        seqid=f1.seqid,
+        biotype=f1.biotype,
+        map=f1.map,
+        name=f1.name,
+        strand="-" if f1.map.reverse else "+",
     )
     value = attrs["map"][:4] if attr == "map" else "different"
     attrs[attr] = value
diff --git a/tests/test_core/test_sequence.py b/tests/test_core/test_sequence.py
index 4abdb351c..90464ccfc 100644
--- a/tests/test_core/test_sequence.py
+++ b/tests/test_core/test_sequence.py
@@ -2598,19 +2598,24 @@ def test_gapped_by_map_segment_iter():
 @pytest.mark.parametrize("sliced", (False, True))
 @pytest.mark.parametrize("start_stop", ((None, None), (3, 7)))
 def test_copied_parent_coordinates(sliced, rev, start_stop):
-    seq = DNA.make_seq("ACGGTGGGAC")
+    orig_name = "orig"
+    seq = DNA.make_seq("ACGGTGGGAC", name=orig_name)
     start, stop = start_stop
     start = start or 0
     stop = stop or len(seq)
     sl = slice(start, stop)
     seq = seq[sl]
+    sliced_name = "sliced"
+    seq.name = sliced_name
+    assert seq.name == sliced_name
     if rev:
         seq = seq.rc()
     copied = seq.copy(sliced=sliced)
+    assert copied.name == sliced_name
     # matches original
     assert copied.parent_coordinates() == seq.parent_coordinates()
-    # and expected
-    assert copied.parent_coordinates() == (start, stop, -1 if rev else 1)
+    # and expected -- the coordinate name always reflects the underlying sequence
+    assert copied.parent_coordinates() == (orig_name, start, stop, -1 if rev else 1)
 
 
 @pytest.mark.parametrize("rev", (False, True))
@@ -2619,8 +2624,9 @@ def test_parent_coordinates(rev):
     seq = seq[1:1]
     if rev:
         seq = seq.rc()
-
-    assert seq.parent_coordinates() == (0, 0, 1)
+    seq.name = "sliced"  # this assignment does not affect the
+    # note that when a sequence has zero length, the parent seqid is None
+    assert seq.parent_coordinates() == (None, 0, 0, 1)
 
 
 def test_seqview_seqid():