Implement ORF searcher with corrected terminal translation

docs/api/orf_searcher.rst

@@ -1,6 +1,49 @@
 ORF searcher
 ============
 
+The V2 package provides an adapter for the tested ``orffinder==1.8`` release.
+It receives nucleotide DNA and returns translated protein candidates; it does
+not return nucleotide ORFs. The pipeline remains responsible for sorting and
+deduplicating those candidates.
+
+The default minimum length is 75 **nucleotides**, matching the package and the
+legacy workflow. For a complete ORF, this includes the terminal stop codon: for
+example, ``ATGAAATAG`` has a 9-nt ORF and yields the 2-aa protein ``MK`` after
+the terminal stop is removed. This is not a protein-length threshold.
+
+By default, only ATG start codons are supported. In the installed 1.8 release,
+custom ``start_codons`` are silently ignored, so the adapter rejects custom
+codons rather than suggesting that they work. Translation uses Biopython's
+standard genetic code; the release accepts but ignores its ``translation_table``
+argument, so the adapter does not expose a translation-table option.
+
+The adapter does not use ``getORFProteins`` for translation: in 1.8 that
+function can truncate a full codon for candidates ending at the sequence
+boundary (and can return an empty protein for ``ATG``). Instead it uses
+``getORFs`` to locate candidates and extracts complete codons from the original
+input in the reported strand direction. It excludes the stop codon for complete
+ORFs, and trims only the 1 or 2 leftover nucleotides from incomplete terminal
+ORFs. No artificial bases are added.
+
+``remove_nested=False`` disables the package's nested-candidate filter. It does
+not guarantee that every overlapping ORF is enumerated; that depends on the
+finder's candidate-generation logic. ``trim_trailing=False`` (the default)
+allows ORFs without a terminal stop at a sequence boundary. The installed
+release internally reports incomplete terminal candidates one nucleotide
+shorter than their actual nucleotide span; the adapter compensates for this
+and enforces ``minimum_length`` against the actual span. For example,
+``ATGAAATG`` has an 8-nt terminal candidate and yields ``MK`` at minimums of 7
+or 8 nt, but not 9 nt. Setting ``trim_trailing=True`` excludes that candidate.
+
+Inputs must be nonempty and contain only case-insensitive IUPAC DNA symbols
+``ACGTRYSWKMBDHVN``. Ambiguous codons can translate to ``X``. Whitespace,
+gaps, RNA ``U``, and other symbols are rejected; input is not padded or
+silently cleaned. ORFFinder scans all six frames itself.
+
 .. autoclass:: lverage.orf_searcher.OrfSearcherTemplate
    :members:
    :show-inheritance:
+
+.. autoclass:: lverage.orf_searcher.OrffinderOrfSearcher
+   :members:
+   :show-inheritance:

docs/singularity_setup.md

@@ -6,6 +6,7 @@ This document tracks all software, external tools, databases, environment config
 
 - Python 3.11
 - Biopython
+- orffinder 1.8 (Python ORF finder)
 - Requests
 
 Developer/documentation dependencies:

pyproject.toml

@@ -11,6 +11,7 @@ requires-python = ">=3.11"
 license = {file = "LICENSE"}
 dependencies = [
     "biopython>=1.85",
+    "orffinder==1.8",
     "requests>=2.32",
 ]
 

src/lverage/orf_searcher.py

@@ -25,6 +25,11 @@
 """
 
 from abc import ABC, abstractmethod
+from collections.abc import Sequence
+
+from Bio.Seq import Seq
+from Bio.SeqRecord import SeqRecord
+from orffinder import orffinder
 
 class OrfSearcherTemplate(ABC):
     """
@@ -54,3 +59,153 @@ def get_orfs(self, sequence : str) -> list[str]:
         """
 
         raise NotImplementedError
+
+
+class OrffinderOrfSearcher(OrfSearcherTemplate):
+    """Find protein-producing ORFs using ``orffinder``.
+
+    Parameters
+    ----------
+    minimum_length : int, optional
+        Minimum ORF length in nucleotides. Complete ORFs include the stop
+        codon in this count. The default of 75 matches ``orffinder`` and the
+        legacy workflow.
+    start_codons : sequence of str, optional
+        Start codons to recognize. Installed ``orffinder`` 1.8 silently
+        ignores custom start codons, so only ATG is supported.
+    remove_nested : bool, optional
+        Whether to ask ``orffinder`` to remove fully nested candidates.
+        Defaults to False, which disables its nested-candidate filter.
+    trim_trailing : bool, optional
+        Whether to discard candidates that reach a sequence end without a
+        stop codon. Defaults to False, matching the legacy workflow.
+
+    Notes
+    -----
+    Input is a nonempty, case-insensitive DNA sequence using the IUPAC
+    alphabet ``ACGTRYSWKMBDHVN``. Ambiguous codons may translate to ``X``.
+    Whitespace, gaps, RNA ``U``, and other symbols are rejected rather than
+    removed or converted. Sequences are passed without padding; the finder
+    scans all six reading frames itself.
+
+    ORFs are translated with the standard genetic code, excluding a complete
+    ORF's terminal stop codon. No translation-table option is exposed because
+    ``orffinder`` 1.8 accepts but ignores that argument. Its protein-producing
+    function can truncate complete codons for incomplete terminal ORFs; this
+    adapter instead extracts codons from the original sequence using the
+    reported strand and start coordinate. Minimum length is checked against
+    the actual nucleotide span, including any incomplete trailing bases.
+    """
+
+    _DNA_ALPHABET = frozenset("ACGTRYSWKMBDHVN")
+
+    def __init__(self,
+                 minimum_length: int = 75,
+                 start_codons: Sequence[str] = ("ATG",),
+                 remove_nested: bool = False,
+                 trim_trailing: bool = False):
+        if isinstance(minimum_length, bool) or not isinstance(minimum_length, int):
+            raise TypeError("minimum_length must be an integer")
+        if minimum_length <= 0:
+            raise ValueError("minimum_length must be greater than zero")
+
+        if isinstance(start_codons, (str, bytes)) or not isinstance(start_codons, Sequence):
+            raise TypeError("start_codons must be a sequence of codon strings")
+        if not start_codons:
+            raise ValueError("start_codons must contain ATG")
+        if any(not isinstance(codon, str) for codon in start_codons):
+            raise TypeError("start_codons must contain only strings")
+        normalized_codons = tuple(codon.upper() for codon in start_codons)
+        if any(len(codon) != 3 or set(codon) - set("ACGT") for codon in normalized_codons):
+            raise ValueError("start codons must be three-base DNA codons")
+        if set(normalized_codons) != {"ATG"}:
+            raise ValueError("orffinder 1.8 supports only the ATG start codon")
+
+        if not isinstance(remove_nested, bool):
+            raise TypeError("remove_nested must be a boolean")
+        if not isinstance(trim_trailing, bool):
+            raise TypeError("trim_trailing must be a boolean")
+
+        self.minimum_length = minimum_length
+        self.start_codons = ("ATG",)
+        self.remove_nested = remove_nested
+        self.trim_trailing = trim_trailing
+
+    def get_orfs(self, sequence: str) -> list[str]:
+        """Return translated ORF candidates from a nucleotide sequence.
+
+        Parameters
+        ----------
+        sequence : str
+            Nonempty DNA sequence containing only IUPAC symbols
+            ``ACGTRYSWKMBDHVN`` (case-insensitive).
+
+        Returns
+        -------
+        list of str
+            Translated protein candidates without terminal stop symbols or
+            empty proteins. Returns an empty list when no usable candidate is
+            found.
+
+        Raises
+        ------
+        TypeError
+            If ``sequence`` is not a string.
+        ValueError
+            If ``sequence`` is empty or contains unsupported symbols.
+        """
+
+        if not isinstance(sequence, str):
+            raise TypeError("sequence must be a string")
+        if not sequence:
+            raise ValueError("sequence must not be empty")
+        invalid_symbols = sorted(set(sequence.upper()) - self._DNA_ALPHABET)
+        if invalid_symbols:
+            raise ValueError(
+                "sequence contains unsupported symbols: " + "".join(invalid_symbols)
+            )
+
+        nucleotide_sequence = Seq(sequence.upper())
+        sequence_length = len(nucleotide_sequence)
+        # orffinder 1.8 counts incomplete terminal candidates one nucleotide
+        # short, so scan one nt below our threshold and enforce it accurately
+        # against the extracted span below.
+        finder_minimum_length = max(1, self.minimum_length - 1)
+        loci = orffinder.getORFs(
+            SeqRecord(nucleotide_sequence),
+            minimum_length=finder_minimum_length,
+            remove_nested=self.remove_nested,
+            trim_trailing=self.trim_trailing,
+        )
+
+        forward = str(nucleotide_sequence)
+        reverse = str(nucleotide_sequence.reverse_complement())
+        proteins = []
+        for locus in loci:
+            is_trailing = locus["trailing"]
+            strand_sequence = forward if locus["sense"] == "+" else reverse
+            start_index = (
+                locus["start"] - 1
+                if locus["sense"] == "+"
+                else sequence_length - locus["start"] + 1
+            )
+
+            if is_trailing:
+                nucleotide_length = sequence_length - start_index
+                if nucleotide_length < self.minimum_length:
+                    continue
+                coding_length = nucleotide_length - nucleotide_length % 3
+            else:
+                nucleotide_length = locus["length"]
+                if nucleotide_length < self.minimum_length:
+                    continue
+                coding_length = nucleotide_length - 3  # Exclude terminal stop.
+
+            if coding_length <= 0:
+                continue
+            coding_sequence = strand_sequence[start_index:start_index + coding_length]
+            protein = str(Seq(coding_sequence).translate()).rstrip("*")
+            if protein:
+                proteins.append(protein)
+
+        return proteins

tests/v2/test_orf_searcher.py

@@ -1,6 +1,12 @@
 import unittest
 
-from lverage.orf_searcher import OrfSearcherTemplate
+from Bio.Seq import Seq
+
+from lverage.domain_scanner import DomainScannerTemplate
+from lverage.motif_database import MotifDBTemplate
+from lverage.orf_searcher import OrfSearcherTemplate, OrffinderOrfSearcher
+from lverage.ortholog_searcher import OrthologSearcherTemplate
+from lverage.pipeline import Lverage, LverageCode
 
 
 class IncompleteOrfSearcher(OrfSearcherTemplate):
@@ -13,6 +19,33 @@ def get_orfs(self, sequence : str) -> list[str]:
         return [sequence]
 
 
+class StubMotifDatabase(MotifDBTemplate):
+
+    name = "Stub Motif Database"
+
+    def search(self, request):
+        return []
+
+    def check_species_validity(self, species_tax_id: int):
+        return True
+
+
+class RecordingDomainScanner(DomainScannerTemplate):
+
+    def __init__(self):
+        self.scanned = []
+
+    def get_domains(self, sequence: str):
+        self.scanned.append(sequence)
+        return []
+
+
+class StubOrthologSearcher(OrthologSearcherTemplate):
+
+    def get_orthologs(self, sequence: str):
+        return []
+
+
 class OrfSearcherTemplateTests(unittest.TestCase):
 
     def test_template_cannot_be_instantiated(self):
@@ -29,5 +62,119 @@ def test_concrete_subclass_returns_orfs(self):
         self.assertEqual(searcher.get_orfs("ATGAAATAG"), ["ATGAAATAG"])
 
 
+class OrffinderOrfSearcherTests(unittest.TestCase):
+
+    def test_finds_forward_and_reverse_strand_proteins(self):
+        searcher = OrffinderOrfSearcher(minimum_length=9)
+
+        self.assertEqual(searcher.get_orfs("ATGAAATAG"), ["MK"])
+        self.assertEqual(searcher.get_orfs("CTATTTCAT"), ["MK"])
+
+    def test_minimum_length_is_inclusive_and_measured_in_nucleotides(self):
+        self.assertEqual(OrffinderOrfSearcher(minimum_length=9).get_orfs("ATGAAATAG"), ["MK"])
+        self.assertEqual(OrffinderOrfSearcher(minimum_length=10).get_orfs("ATGAAATAG"), [])
+
+    def test_incomplete_terminal_orfs_follow_trimming_and_length_settings(self):
+        self.assertEqual(OrffinderOrfSearcher(minimum_length=7).get_orfs("ATGAAATG"), ["MK"])
+        self.assertEqual(OrffinderOrfSearcher(minimum_length=8).get_orfs("ATGAAATG"), ["MK"])
+        self.assertEqual(OrffinderOrfSearcher(minimum_length=9).get_orfs("ATGAAATG"), [])
+        self.assertEqual(
+            OrffinderOrfSearcher(minimum_length=7, trim_trailing=True).get_orfs("ATGAAATG"),
+            [],
+        )
+
+    def test_incomplete_orfs_preserve_full_codons_for_each_suffix_remainder(self):
+        cases = (
+            ("ATGAAAGGG", 0, "MKG"),
+            ("ATGAAAGGGA", 1, "MKG"),
+            ("ATGAAAGGGAA", 2, "MKG"),
+            ("ATGAAAGGGAAA", 0, "MKGK"),
+        )
+        for nucleotide_sequence, expected_remainder, expected_protein in cases:
+            with self.subTest(nucleotide_sequence=nucleotide_sequence):
+                self.assertEqual(len(nucleotide_sequence) % 3, expected_remainder)
+                self.assertEqual(
+                    OrffinderOrfSearcher(minimum_length=9).get_orfs(nucleotide_sequence),
+                    [expected_protein],
+                )
+                self.assertEqual(
+                    OrffinderOrfSearcher(minimum_length=9).get_orfs(
+                        str(Seq(nucleotide_sequence).reverse_complement())
+                    ),
+                    [expected_protein],
+                )
+
+    def test_incomplete_orf_minimum_length_uses_actual_nucleotide_span(self):
+        self.assertEqual(OrffinderOrfSearcher(minimum_length=9).get_orfs("ATGAAAGGG"), ["MKG"])
+        self.assertEqual(OrffinderOrfSearcher(minimum_length=10).get_orfs("ATGAAAGGG"), [])
+        self.assertEqual(OrffinderOrfSearcher(minimum_length=10).get_orfs("ATGAAAGGGA"), ["MKG"])
+        self.assertEqual(OrffinderOrfSearcher(minimum_length=11).get_orfs("ATGAAAGGGA"), [])
+
+    def test_start_codon_alone_produces_one_residue_not_an_empty_candidate(self):
+        self.assertEqual(OrffinderOrfSearcher(minimum_length=1).get_orfs("ATG"), ["M"])
+
+    def test_removes_terminal_stop_but_preserves_internal_translation(self):
+        searcher = OrffinderOrfSearcher(minimum_length=9)
+
+        self.assertEqual(searcher.get_orfs("ATGAAATAG"), ["MK"])
+        self.assertEqual(searcher.get_orfs("ATGNNNTAA"), ["MX"])
+
+    def test_uses_the_installed_standard_translation_table(self):
+        self.assertEqual(OrffinderOrfSearcher(minimum_length=6).get_orfs("ATGTGATAG"), ["M"])
+
+    def test_returns_empty_list_when_no_orf_is_found(self):
+        searcher = OrffinderOrfSearcher(minimum_length=9)
+
+        self.assertEqual(searcher.get_orfs("CCCCCCCCCC"), [])
+
+    def test_remove_nested_false_disables_the_nested_candidate_filter(self):
+        sequence = "AATATCAATGGGTCAGCGGTACGCCGGGGGTACCGTAAGGCACATAGTATTAGAATGATA"
+
+        unfiltered = OrffinderOrfSearcher(minimum_length=9).get_orfs(sequence)
+        filtered = OrffinderOrfSearcher(minimum_length=9, remove_nested=True).get_orfs(sequence)
+
+        self.assertEqual(unfiltered, ["MGQRYAGGTVRHIVLE", "MCLTVPPAYR"])
+        self.assertEqual(filtered, ["MGQRYAGGTVRHIVLE"])
+
+    def test_rejects_invalid_sequences(self):
+        searcher = OrffinderOrfSearcher(minimum_length=9)
+
+        for invalid in ("", " ", "ATGUAA", "ATG-AA"):
+            with self.subTest(invalid=invalid):
+                with self.assertRaises(ValueError):
+                    searcher.get_orfs(invalid)
+        for invalid in (None, 123):
+            with self.subTest(invalid=invalid):
+                with self.assertRaises(TypeError):
+                    searcher.get_orfs(invalid)
+
+    def test_validates_configuration(self):
+        with self.assertRaises(TypeError):
+            OrffinderOrfSearcher(minimum_length=True)
+        with self.assertRaises(ValueError):
+            OrffinderOrfSearcher(minimum_length=0)
+        with self.assertRaises(ValueError):
+            OrffinderOrfSearcher(start_codons=["GTG"])
+        with self.assertRaises(ValueError):
+            OrffinderOrfSearcher(start_codons=[])
+        with self.assertRaises(TypeError):
+            OrffinderOrfSearcher(start_codons="ATG")
+        with self.assertRaises(TypeError):
+            OrffinderOrfSearcher(remove_nested=1)
+
+    def test_adapter_passes_protein_candidates_through_pipeline(self):
+        domain_scanner = RecordingDomainScanner()
+        pipeline = Lverage(
+            [StubMotifDatabase()],
+            OrffinderOrfSearcher(minimum_length=9),
+            domain_scanner,
+            StubOrthologSearcher(),
+        )
+
+        self.assertEqual(pipeline.run("ATGAAATAG"), [])
+        self.assertEqual(domain_scanner.scanned, ["MK"])
+        self.assertEqual(pipeline.lverage_code, LverageCode.NO_VALID_DOMAIN)
+
+
 if __name__ == "__main__":
     unittest.main()