Repository: cogent3/cogent3
Pull request #2297 (2025-05-01 02:09:28)

ENH: added new duplicated_seqs() method
[NEW] on new type sequence collections. Returns list of list of sequence
    names whose hashes are equal.

## Summary by Sourcery

Add a new method `duplicated_seqs()` to SequenceCollection and Alignment classes to identify sequences with identical content

New Features:
- Added `duplicated_seqs()` method to identify sequences with identical content in SequenceCollection and Alignment classes

Enhancements:
- Implemented a hash-based approach to detect duplicate sequences
- Supports both aligned and unaligned sequence collections

Tests:
- Added comprehensive test cases for `duplicated_seqs()` method with various sequence scenarios


Patch:
diff --git a/changelog.d/20250406_125550_Gavin.Huttley.md b/changelog.d/20250406_125550_Gavin.Huttley.md
index c0369a52e..6cbec8c7a 100644
--- a/changelog.d/20250406_125550_Gavin.Huttley.md
+++ b/changelog.d/20250406_125550_Gavin.Huttley.md
@@ -30,6 +30,11 @@ Uncomment the section that is right (remove the HTML comment wrapper).
   with cogent3! The main brca1 alignment and associated Murphy tree are now
   included with the package. Use `cogent3.available_datasets()` and
   `cogent3.get_dataset()` to create the instances that use these.
+- Added `duplicated_seqs()` method to new type SequenceCollection and Alignments.
+  As the name implies, it identifies the names of sequences that are equal. On
+  Alignments, the gapped sequence is used. Added companion method
+  `drop_duplicated_seqs()` which does as the name implies. The first name in each
+  group of duplicates is retained.
 
 <!--
 ### BUG
diff --git a/src/cogent3/core/new_alignment.py b/src/cogent3/core/new_alignment.py
index 8e78deabb..2f8a39b09 100644
--- a/src/cogent3/core/new_alignment.py
+++ b/src/cogent3/core/new_alignment.py
@@ -1,7 +1,9 @@
 from __future__ import annotations
 
+import collections
 import contextlib
 import dataclasses
+import hashlib
 import json
 import os
 import re
@@ -93,6 +95,18 @@ MolTypes = str | new_moltype.MolType
 EPS = 1e-6
 
 
+def array_hash64(data: numpy.ndarray) -> int:
+    """returns 64-bit hash of numpy array.
+
+    Notes
+    -----
+    This function does not introduce randomisation and so
+    is reproducible between processes.
+    """
+    h = hashlib.md5(data.tobytes(), usedforsecurity=False)
+    return int.from_bytes(h.digest()[:8], byteorder="little")
+
+
 class _SeqNamer:
     def __init__(
         self,
@@ -2405,9 +2419,32 @@ class SequenceCollection:
                 raise TypeError(msg)
             self._repr_policy["wrap"] = wrap
 
+    def duplicated_seqs(self) -> list[list[str]]:
+        """returns the names of duplicated sequences"""
+        seq_hashes = collections.defaultdict(list)
+        for s in self.seqs:
+            seq_hashes[array_hash64(numpy.array(s))].append(s.name)
+        return [v for v in seq_hashes.values() if len(v) > 1]
+
+    def drop_duplicated_seqs(self) -> typing_extensions.Self:
+        """returns self without duplicated sequences
+
+        Notes
+        -----
+        Retains the first sequence of each duplicte group.
+        """
+        dupes = self.duplicated_seqs()
+        if not dupes:
+            return self
+
+        omit = []
+        for group in dupes:
+            omit.extend(group[1:])
+        return self.take_seqs(omit, negate=True)
+
 
 @register_deserialiser(get_object_provenance(SequenceCollection))
-def deserialise_sequence_collection(data) -> SequenceCollection:
+def deserialise_sequence_collection(data: dict) -> SequenceCollection:
     return SequenceCollection.from_rich_dict(data)
 
 
@@ -6892,6 +6929,19 @@ class Alignment(SequenceCollection):
             s.name: s.seq.strand_symmetry(motif_length=motif_length) for s in self.seqs
         }
 
+    def duplicated_seqs(self) -> list[list[str]]:
+        """returns the names of duplicated sequences
+
+        Notes
+        -----
+        The gapped sequence is used.
+        """
+        if not len(self):
+            # all have zero lengths
+            return [] if self.num_seqs < 2 else [list(self.names)]
+
+        return super().duplicated_seqs()
+
 
 @register_deserialiser(get_object_provenance(Alignment))
 def deserialise_alignment(data: dict[str, str | dict[str, str]]) -> Alignment:


Test patch:
diff --git a/tests/test_core/test_new_alignment.py b/tests/test_core/test_new_alignment.py
index 0148f76c9..fda5050a9 100644
--- a/tests/test_core/test_new_alignment.py
+++ b/tests/test_core/test_new_alignment.py
@@ -5900,3 +5900,94 @@ def test_coll_storage_degap_propagates_type():
     dg = seqcoll.degap(storage_backend="h5seqs_unaligned")
     assert not isinstance(dg.storage, builtin)
     assert isinstance(dg.storage, seqcoll.storage.__class__)
+
+
+@pytest.mark.parametrize(
+    "mk_cls", [new_alignment.make_aligned_seqs, new_alignment.make_unaligned_seqs]
+)
+def test_duplicated_seqs_duplicates(mk_cls):
+    data = {
+        "seq1": "ATC-G",
+        "seq2": "TAGCC",
+        "seq1-1": "ATC-G",
+    }
+    seqcoll = mk_cls(data, moltype="dna")
+    got = [set(v) for v in seqcoll.duplicated_seqs()]
+    assert got == [{"seq1", "seq1-1"}]
+
+    # all seqs length=0
+    data = {
+        "seq1": "",
+        "seq2": "",
+    }
+    seqcoll = mk_cls(data, moltype="dna")
+    got = [set(v) for v in seqcoll.duplicated_seqs()]
+    assert got == [{"seq1", "seq2"}]
+
+    data = {
+        "seq1": "ATC-G",
+        "seq2": "TAGCC",
+        "seq1-1": "ATC-G",
+        "seq2-1": "TAGCC",
+        "seq2-2": "TAGCC",
+    }
+    seqcoll = mk_cls(data, moltype="dna")
+    got = {frozenset(v) for v in seqcoll.duplicated_seqs()}
+    assert got == {
+        frozenset({"seq1", "seq1-1"}),
+        frozenset({"seq2", "seq2-1", "seq2-2"}),
+    }
+
+
+@pytest.mark.parametrize(
+    "mk_cls", [new_alignment.make_aligned_seqs, new_alignment.make_unaligned_seqs]
+)
+def test_duplicated_seqs_no_duplicates(mk_cls):
+    data = {
+        "seq1": "ATC-G",
+        "seq2": "TAGCC",
+    }
+    seqcoll = mk_cls(data, moltype="dna")
+    got = [set(v) for v in seqcoll.duplicated_seqs()]
+    assert got == []
+
+
+@pytest.mark.parametrize(
+    "mk_cls", [new_alignment.make_aligned_seqs, new_alignment.make_unaligned_seqs]
+)
+def test_drop_duplicated_seqs(mk_cls):
+    data = {
+        "seq1": "ATC-G",
+        "seq2": "TAGCC",
+        "seq1-1": "ATC-G",
+        "seq2-1": "TAGCC",
+        "seq2-2": "TAGCC",
+        "seq2-3": "TAGCC",
+    }
+    seqcoll = mk_cls(data, moltype="dna")
+    got = seqcoll.drop_duplicated_seqs()
+    names = set(got.names)
+    assert len(names & {"seq1", "seq1-1"}) == 1
+    assert len(names & {"seq2", "seq2-1", "seq2-2", "seq2-3"}) == 1
+
+    # on a collection with zero length seqs
+    data = {
+        "seq1": "",
+        "seq2": "",
+    }
+    seqcoll = mk_cls(data, moltype="dna")
+    got = seqcoll.drop_duplicated_seqs()
+
+
+@pytest.mark.parametrize(
+    "mk_cls", [new_alignment.make_aligned_seqs, new_alignment.make_unaligned_seqs]
+)
+def test_drop_duplicated_seqs_no_dupes(mk_cls):
+    data = {
+        "seq1": "ATC-G",
+        "seq2": "TAGCC",
+        "seq3": "T-GCC",
+    }
+    seqcoll = mk_cls(data, moltype="dna")
+    got = seqcoll.drop_duplicated_seqs()
+    assert got is seqcoll