Repository: cogent3/cogent3
Pull request #1714 (2024-02-12 21:14:21)

Feature querying now works when sequences renamed


Patch:
diff --git a/changelog.d/20240213_075007_Gavin.Huttley.md b/changelog.d/20240213_075007_Gavin.Huttley.md
new file mode 100644
index 000000000..822293351
--- /dev/null
+++ b/changelog.d/20240213_075007_Gavin.Huttley.md
@@ -0,0 +1,50 @@
+<!--
+A new scriv changelog fragment.
+
+Uncomment the section that is right (remove the HTML comment wrapper).
+-->
+
+
+### Contributors
+
+- Katherine Caley for reviewing PRs and code contributions
+  related to SeqView taking on a more prominent role in recording
+  history of operations on sequences.
+
+
+
+### ENH
+
+- Annotation databases are now preserved after renaming sequences.
+  This is made possible because SeqView has the seqid attribute which
+  is preserved and is independent of the name of the enclosing Sequence.
+- Sequence.parent_coordinates() now returns information related to its
+  parent. For example, we have a sequence with name "1" which we then
+  sliced and rename the slice to "gene". Calling gene.parent_coordinates()
+  would return "1", slice start, slice end, strand.
+
+
+<!--
+### BUG
+
+- A bullet item for the BUG category.
+
+-->
+<!--
+### DOC
+
+- A bullet item for the DOC category.
+
+-->
+<!--
+### Deprecations
+
+- A bullet item for the Deprecations category.
+
+-->
+<!--
+### Discontinued
+
+- A bullet item for the Discontinued category.
+
+-->
diff --git a/src/cogent3/__init__.py b/src/cogent3/__init__.py
index ac4d511a7..7cd610989 100644
--- a/src/cogent3/__init__.py
+++ b/src/cogent3/__init__.py
@@ -95,9 +95,7 @@ def make_seq(seq, name=None, moltype=None):
     """
     moltype = moltype or "text"
     moltype = get_moltype(moltype)
-    seq = moltype.make_seq(seq)
-    if name is not None:
-        seq.name = name
+    seq = moltype.make_seq(seq, name=name)
     return seq
 
 
diff --git a/src/cogent3/core/alignment.py b/src/cogent3/core/alignment.py
index 2277b97f0..77a317b1a 100644
--- a/src/cogent3/core/alignment.py
+++ b/src/cogent3/core/alignment.py
@@ -1732,7 +1732,7 @@ class _SequenceCollectionBase:
             new[new_name] = new_seq
         result = self.__class__(data=new, info=self.info, moltype=self.moltype)
         result.info.name_map = name_map
-        result.annotation_db = DEFAULT_ANNOTATION_DB()
+        result.annotation_db = self.annotation_db
         return result
 
     @UI.display_wrap
@@ -2072,11 +2072,10 @@ class SequenceCollection(_SequenceCollectionBase):
         if not self.annotation_db:
             return None
 
-        seqids = [seqid] if isinstance(seqid, str) else seqid
-        if seqids is None:
-            seqids = self.names
-
-        if seqid and not set(seqids) & set(self.names):
+        seqid_to_seqname = {seq.parent_coordinates()[0]: seq.name for seq in self.seqs}
+        if seqid and (
+            seqid not in seqid_to_seqname and seqid not in seqid_to_seqname.values()
+        ):
             raise ValueError(f"unknown {seqid=}")
 
         for feature in self.annotation_db.get_features_matching(
@@ -2086,7 +2085,8 @@ class SequenceCollection(_SequenceCollectionBase):
             on_alignment=False,
             allow_partial=allow_partial,
         ):
-            seq = self.named_seqs[feature["seqid"]]
+            seqname = seqid_to_seqname[feature["seqid"]]
+            seq = self.named_seqs[seqname]
             if offset := seq.annotation_offset:
                 feature["spans"] = (array(feature["spans"]) - offset).tolist()
             yield seq.make_feature(feature, self)
@@ -4668,6 +4668,7 @@ class Alignment(AlignmentI, SequenceCollection):
             new = self.__class__(
                 seqs, name=self.name, info=self.info, moltype=self.moltype
             )
+            new.annotation_db = self.annotation_db
 
         if isinstance(index, (list, tuple)):
             raise TypeError(f"cannot slice using {type(index)}")
@@ -5266,15 +5267,27 @@ class Alignment(AlignmentI, SequenceCollection):
         if self.annotation_db is None:
             return None
 
+        seqid_to_seqname = {
+            seq.data.parent_coordinates()[0]: seq.name for seq in self.seqs
+        }
+
         seqids = [seqid] if isinstance(seqid, str) else seqid
         if seqids is None:
-            seqids = self.names
-
-        if seqid and not set(seqids) & set(self.names):
+            seqids = tuple(seqid_to_seqname)
+        elif set(seqids) & set(self.names):
+            # we've been given seq names, convert to parent names
+            seqids = [
+                self.named_seqs[seqid].data.parent_coordinates()[0] for seqid in seqids
+            ]
+        elif seqids and set(seqids) <= seqid_to_seqname.keys():
+            # already correct
+            pass
+        else:
             raise ValueError(f"unknown {seqid=}")
 
         for seqid in seqids:
-            seq = self.named_seqs[seqid]
+            seqname = seqid_to_seqname[seqid]
+            seq = self.named_seqs[seqname]
             # we use parent seqid, stored on SeqView
             parent_id, start, end, _ = seq.data.parent_coordinates()
             offset = seq.data.annotation_offset
diff --git a/src/cogent3/core/sequence.py b/src/cogent3/core/sequence.py
index d72cfa81f..390fc62e8 100644
--- a/src/cogent3/core/sequence.py
+++ b/src/cogent3/core/sequence.py
@@ -940,19 +940,25 @@ class Sequence(SequenceI):
         start = start + len(self) if start < 0 else start
         stop = stop + len(self) if stop < 0 else stop
 
-        # sort start / stop
         start, stop = (start, stop) if start < stop else (stop, start)
 
         # note: offset is handled by absolute_position
         # we set include_boundary=False because start is inclusive indexing,
         # i,e., the start cannot be equal to the length of the view
+        (
+            # parent_id can differ from self.name if there was a
+            # rename operation
+            parent_id,
+            *_,
+            strand,
+        ) = self.parent_coordinates()
         query_start = self._seq.absolute_position(start, include_boundary=False)
         # 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)
 
-        reversed = query_end < query_start
-        if reversed:
+        rev_strand = strand == -1
+        if rev_strand:
             query_start, query_end = query_end, query_start
 
         query_start = max(query_start, 0)
@@ -973,7 +979,7 @@ class Sequence(SequenceI):
         # flag which comes back from the db
 
         for feature in self.annotation_db.get_features_matching(
-            seqid=self.name,
+            seqid=parent_id,
             name=name,
             biotype=biotype,
             start=query_start,
@@ -982,13 +988,14 @@ class Sequence(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
+            # and we expect them to make a feature manually wrt to their
+            # current view
             spans = array(feature.pop("spans"), dtype=int)
             for i, v in enumerate(spans.ravel()):
                 rel_pos = self._seq.relative_position(v)
                 spans.ravel()[i] = rel_pos
 
-            if reversed:
+            if rev_strand:
                 # see above comment
                 spans = len(self) - spans
 


Test patch:
diff --git a/tests/test_core/test_alignment.py b/tests/test_core/test_alignment.py
index c37630f27..5dc949871 100644
--- a/tests/test_core/test_alignment.py
+++ b/tests/test_core/test_alignment.py
@@ -2549,14 +2549,17 @@ def test_deepcopy_aligned(cls):
 
 
 @pytest.mark.parametrize("cls", (SequenceCollection, Alignment))
-def test_seq_rename_drops_annotations(cls):
+def test_seq_rename_preserves_annotations(cls):
     """rename seqs discards all annotations"""
     data = {"seq1": "ACGTACGTA", "seq2": "ACCGAA---", "seq3": "ACGTACGTT"}
     seqs = cls(data, moltype=DNA)
     seqs.add_feature(seqid="seq1", biotype="exon", name="fred", spans=[(3, 8)])
     assert seqs.annotation_db is not None
     new = seqs.rename_seqs(lambda x: x.upper())
-    assert not len(new.annotation_db)
+    assert len(new.annotation_db) == 1
+    assert len(list(new.get_features(biotype="exon")))
+    # using original seq name should also work
+    assert len(list(new.get_features(seqid="seq1")))
 
 
 @pytest.mark.parametrize(
diff --git a/tests/test_core/test_annotation.py b/tests/test_core/test_annotation.py
index 0fcfcbc1a..e8f384cb8 100644
--- a/tests/test_core/test_annotation.py
+++ b/tests/test_core/test_annotation.py
@@ -3,6 +3,7 @@ import unittest
 import pytest
 
 from cogent3 import DNA, load_seq, make_aligned_seqs, make_unaligned_seqs
+from cogent3.core.alignment import Alignment, SequenceCollection
 from cogent3.core.location import Map, Span
 
 
@@ -299,3 +300,56 @@ def test_gbdb_get_children_get_parent(DATA_DIR):
     (child,) = list(orig.get_children("CDS"))
     parent, *_ = list(child.get_parent())
     assert parent == orig
+
+
+@pytest.mark.parametrize("rev", (False, True))
+def test_features_survives_seq_rename(rev):
+    segments = ["A" * 10, "C" * 10, "T" * 5, "C" * 5, "A" * 5]
+
+    seq = DNA.make_seq("".join(segments), name="original")
+    gene = seq.add_feature(biotype="gene", name="gene1", spans=[(10, 20), (25, 30)])
+    domain = seq.add_feature(
+        biotype="domain", name="domain1", spans=[(20, 25)], strand="-"
+    )
+    sliced = seq[5:-3]
+    sliced.name = "sliced"
+    sliced = sliced.rc() if rev else sliced
+    got = list(sliced.get_features(name="gene1"))[0]
+    assert str(got.get_slice()) == str(gene.get_slice())
+    got = list(sliced.get_features(name="domain1"))[0]
+    assert str(got.get_slice()) == str(domain.get_slice())
+
+
+@pytest.mark.parametrize("rev", (False, True))
+@pytest.mark.parametrize("cls", (SequenceCollection, Alignment))
+def test_features_survives_aligned_seq_rename(rev, cls):
+    segments = ["A" * 10, "C" * 10, "T" * 5, "C" * 5, "A" * 5]
+
+    seqs = cls({"original": "".join(segments)}, moltype="dna")
+    seqs.annotation_db.add_feature(
+        seqid="original", biotype="gene", name="gene1", spans=[(10, 20), (25, 30)]
+    )
+    seqs.annotation_db.add_feature(
+        seqid="original", biotype="domain", name="domain1", spans=[(20, 25)], strand="-"
+    )
+    seqs = seqs.rename_seqs(lambda x: "newname")
+    assert seqs.names == ["newname"]
+    seqs = seqs.rc() if rev else seqs
+    # quite different behaviour from Alignment and SequenceCollection
+    # so we convert to string to make comparison simpler
+    got = list(seqs.get_features(name="gene1"))[0]
+    sliced = str(got.get_slice()).splitlines()[-1]
+    assert sliced == "C" * 15
+
+
+@pytest.mark.parametrize("cls", (SequenceCollection, Alignment))
+def test_features_invalid_seqid(cls):
+    segments = ["A" * 10, "C" * 10, "T" * 5, "C" * 5, "A" * 5]
+
+    seqs = cls({"original": "".join(segments)}, moltype="dna")
+    seqs.annotation_db.add_feature(
+        seqid="original", biotype="domain", name="domain1", spans=[(20, 25)], strand="-"
+    )
+    with pytest.raises(ValueError):
+        # seqid does not exist
+        list(seqs.get_features(name="gene1", seqid="blah"))
diff --git a/tests/test_core/test_sequence.py b/tests/test_core/test_sequence.py
index 90464ccfc..52664911e 100644
--- a/tests/test_core/test_sequence.py
+++ b/tests/test_core/test_sequence.py
@@ -2709,3 +2709,8 @@ def test_sequences_propogates_seqid():
     seq = Sequence(SeqView("ACGGTGGGAC"), name="seq_name")
     assert seq.name == "seq_name"
     assert seq._seq.seqid == None
+
+
+def test_make_seq_assigns_to_seqview():
+    seq = cogent3.make_seq("ACGT", name="s1")
+    assert seq.name == seq._seq.seqid == "s1"