Repository: cogent3/cogent3
Pull request #2167 (2024-12-13 01:22:17)

MAINT: new collections and seqs auto-convert nucleic acid seqs
## Summary by Sourcery

Enhance sequence and alignment creation functions to automatically convert nucleic acid sequences to the appropriate molecular type, and add tests to verify this behavior.

Enhancements:
- Add automatic conversion of nucleic acid sequences to the appropriate molecular type in sequence and alignment creation functions.

Tests:
- Introduce tests to verify the automatic conversion of nucleic acid sequences to the correct molecular type for both aligned and unaligned sequences.


Patch:
diff --git a/src/cogent3/core/new_alignment.py b/src/cogent3/core/new_alignment.py
index dfbb3d348..b0ece4bca 100644
--- a/src/cogent3/core/new_alignment.py
+++ b/src/cogent3/core/new_alignment.py
@@ -2381,7 +2381,7 @@ def _(
 
 @coerce_to_raw_seq_data.register
 def _(seq: str, moltype: new_moltype.MolType, name: OptStr = None) -> raw_seq_data:
-    return raw_seq_data(seq=seq, name=name)
+    return coerce_to_raw_seq_data(seq.encode("utf8"), moltype, name)
 
 
 @coerce_to_raw_seq_data.register
@@ -2395,6 +2395,8 @@ def _(
 
 @coerce_to_raw_seq_data.register
 def _(seq: bytes, moltype: new_moltype.MolType, name: OptStr = None) -> raw_seq_data:
+    seq = seq.upper()
+    seq = moltype.coerce_to(seq) if moltype.coerce_to else seq
     return raw_seq_data(seq=seq, name=name)
 
 
@@ -3666,21 +3668,18 @@ class AlignedSeqsData(AlignedSeqsDataABC):
                 names=self.names,
             )
 
-        old = self.alphabet.as_bytes()
-        new = alphabet.as_bytes()
-        convert_old_to_bytes = new_alphabet.array_to_bytes(old)
-        convert_bytes_to_new = new_alphabet.bytes_to_array(
-            new,
-            dtype=new_alphabet.get_array_type(len(new)),
-        )
         gapped = numpy.empty(
             (len(self.names), self.align_len),
             dtype=self.alphabet.dtype,
         )
 
         for i in range(len(self.names)):
-            as_new_alpha = convert_bytes_to_new(convert_old_to_bytes(self._gapped[i]))
-
+            seq_data = self._gapped[i]
+            as_new_alpha = self.alphabet.convert_seq_array_to(
+                seq=seq_data,
+                alphabet=alphabet,
+                check_valid=False,
+            )
             if check_valid and not alphabet.is_valid(as_new_alpha):
                 raise new_moltype.MolTypeError(
                     f"Changing from old alphabet={self.alphabet} to new "
diff --git a/src/cogent3/core/new_alphabet.py b/src/cogent3/core/new_alphabet.py
index 39df257b5..9b310bd66 100644
--- a/src/cogent3/core/new_alphabet.py
+++ b/src/cogent3/core/new_alphabet.py
@@ -57,7 +57,13 @@ class AlphabetError(TypeError): ...
 class convert_alphabet:
     """convert one character set into another"""
 
-    def __init__(self, src: bytes, dest: bytes, delete: OptBytes = None):
+    def __init__(
+        self,
+        src: bytes,
+        dest: bytes,
+        delete: OptBytes = None,
+        allow_duplicates: bool = False,
+    ):
         """
         Parameters
         ----------
@@ -67,6 +73,8 @@ class convert_alphabet:
             characters in src will be mapped to these in order
         delete
             characters to be deleted from the result
+        allow_duplicates
+            whether to allow duplicates in src and dest
 
         Notes
         -----
@@ -84,7 +92,8 @@ class convert_alphabet:
         if len(src) != len(dest):
             raise ValueError(f"length of src={len(src)} != length of new {len(dest)}")
 
-        consistent_words(src, length=1)
+        if not allow_duplicates:
+            consistent_words(src, length=1)
 
         self._table = b"".maketrans(src, dest)
         self._delete = delete or b""
@@ -526,7 +535,12 @@ class CharAlphabet(tuple, AlphabetABC, MonomerAlphabetABC):
             motif_subset = [m for m in self if m not in motif_subset]
         gap = self.gap_char if self.gap_char in motif_subset else None
         missing = self.missing_char if self.missing_char in motif_subset else None
-        return self.__class__(tuple(motif_subset), gap=gap, missing=missing)
+        return make_alphabet(
+            chars=tuple(motif_subset),
+            gap=gap,
+            missing=missing,
+            moltype=self.moltype,
+        )
 
     def convert_seq_array_to(
         self,
@@ -591,16 +605,31 @@ def make_converter(
     common = b"".join(src_bytes & other_bytes)
 
     diff = b"".join(src_bytes - other_bytes)
-    src_indices = numpy.concatenate(
-        [src_alpha.to_indices(common), src_alpha.to_indices(diff)],
-    )
-    dest_indices = numpy.concatenate(
-        [
-            dest_alpha.to_indices(common),
-            numpy.array([len(dest_alpha)] * len(diff), dtype=dest_alpha.dtype),
-        ],
+    src_indices = [src_alpha.to_indices(common), src_alpha.to_indices(diff)]
+    dest_indices = [
+        dest_alpha.to_indices(common),
+        numpy.array([len(dest_alpha)] * len(diff), dtype=dest_alpha.dtype),
+    ]
+
+    # if we are going to a nucleic acid alphabet, we need to add an alias
+    # mapping for the character NOT present in that alphabet, i.e. if dest is
+    # RNA, we need to map T to U, if dest is DNA we need to map U to T
+    # we do this by adding these on at the end of the common characters
+    if dest_alpha.moltype.is_nucleic and "T" in dest_alpha:
+        src_indices.append(src_alpha.to_indices("U"))
+        dest_indices.append(dest_alpha.to_indices("T"))
+    elif dest_alpha.moltype.is_nucleic:
+        src_indices.append(src_alpha.to_indices("T"))
+        dest_indices.append(dest_alpha.to_indices("U"))
+
+    src_indices = numpy.concatenate(src_indices)
+    dest_indices = numpy.concatenate(dest_indices)
+
+    return convert_alphabet(
+        src_indices.tobytes(),
+        dest_indices.tobytes(),
+        allow_duplicates=True,
     )
-    return convert_alphabet(src_indices.tobytes(), dest_indices.tobytes())
 
 
 @register_deserialiser(get_object_provenance(CharAlphabet))


Test patch:
diff --git a/tests/test_core/test_new_alignment.py b/tests/test_core/test_new_alignment.py
index bca93ca43..7f0783453 100644
--- a/tests/test_core/test_new_alignment.py
+++ b/tests/test_core/test_new_alignment.py
@@ -5507,3 +5507,32 @@ def test_aligned_from_indel_map_and_aligned_seq_view():
     )
     new_al = new_alignment.Aligned.from_map_and_aligned_data_view(new_map, al.data)
     assert str(new_al) == "---AC--GTC"
+
+
+@pytest.mark.parametrize(
+    "mk_cls",
+    [new_alignment.make_unaligned_seqs, new_alignment.make_aligned_seqs],
+)
+@pytest.mark.parametrize(("moltype", "seq"), [("dna", "AUGC"), ("rna", "ATGC")])
+def test_coerce_moltype(mk_cls, moltype, seq):
+    coll = mk_cls({"s1": seq}, moltype=moltype)
+    assert coll.moltype.name == moltype
+    expect = "ATGC" if moltype == "dna" else "AUGC"
+    assert str(coll.get_seq("s1")) == expect
+
+
+@pytest.mark.parametrize(
+    "mk_cls",
+    [new_alignment.make_unaligned_seqs, new_alignment.make_aligned_seqs],
+)
+def test_coerce_moltype_obj(mk_cls):
+    data = dict(
+        [
+            ("s1", "ACGTA"),
+            ("s2", "ACGTC"),
+            ("s3", "ACGTT"),
+            ("s4", "ACGAT"),
+        ],
+    )
+    coll = mk_cls(data, moltype="text")
+    coll = coll.to_moltype("rna")
diff --git a/tests/test_core/test_new_alphabet.py b/tests/test_core/test_new_alphabet.py
index 7efd6be64..4ec4356c6 100644
--- a/tests/test_core/test_new_alphabet.py
+++ b/tests/test_core/test_new_alphabet.py
@@ -734,6 +734,7 @@ def test_get_subset():
     got = dna_alpha.get_subset("-", excluded=True)
     assert got.gap_char is None
     assert len(got) == 4
+    assert got.moltype is new_moltype.DNA
 
 
 def test_get_subset_invalid():
diff --git a/tests/test_core/test_new_sequence.py b/tests/test_core/test_new_sequence.py
index 20eaecc7f..04f32b740 100644
--- a/tests/test_core/test_new_sequence.py
+++ b/tests/test_core/test_new_sequence.py
@@ -154,10 +154,6 @@ def test_sequence_to_moltype():
     trev = list(got.get_features(name="trev"))[0]
     assert str(got[trev]) == "AAAA"
 
-    # should raise exception if moltype not compatible with sequence data
-    with pytest.raises(new_moltype.MolTypeError):
-        s.to_moltype("rna")
-
     # calling with a null object should raise an exception
     with pytest.raises(ValueError):
         s.to_moltype(None)
@@ -3015,3 +3011,11 @@ def test_make_seq_wrong_order_alpha():
 def test_make_seq_from_types(raw_seq):
     seq = new_moltype.DNA.make_seq(seq=raw_seq)
     assert str(seq) == "GGTAC"
+
+
+@pytest.mark.parametrize(("moltype", "seq"), [("dna", "AUGC"), ("rna", "ATGC")])
+def test_coerce_moltype(moltype, seq):
+    seq = cogent3.make_seq(seq=seq, moltype=moltype, new_type=True)
+    assert seq.moltype.name == moltype
+    expect = "ATGC" if moltype == "dna" else "AUGC"
+    assert str(seq) == expect