From b37440467f6e264cd826a44bdd69f65c7fabb85b Mon Sep 17 00:00:00 2001 From: Thomas Shin Date: Sun, 4 Oct 2026 09:08:19 -0400 Subject: [PATCH] Implement ORF searcher with corrected terminal translation --- docs/api/orf_searcher.rst | 43 ++++++++++ docs/singularity_setup.md | 1 + pyproject.toml | 1 + src/lverage/orf_searcher.py | 155 ++++++++++++++++++++++++++++++++++ tests/v2/test_orf_searcher.py | 149 +++++++++++++++++++++++++++++++- 5 files changed, 348 insertions(+), 1 deletion(-) diff --git a/docs/api/orf_searcher.rst b/docs/api/orf_searcher.rst index 2a2014d..72ec8f4 100644 --- a/docs/api/orf_searcher.rst +++ b/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: diff --git a/docs/singularity_setup.md b/docs/singularity_setup.md index e909b01..497deb7 100644 --- a/docs/singularity_setup.md +++ b/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: diff --git a/pyproject.toml b/pyproject.toml index 95da89d..10ca740 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -11,6 +11,7 @@ requires-python = ">=3.11" license = {file = "LICENSE"} dependencies = [ "biopython>=1.85", + "orffinder==1.8", "requests>=2.32", ] diff --git a/src/lverage/orf_searcher.py b/src/lverage/orf_searcher.py index e5ae3ed..2536908 100644 --- a/src/lverage/orf_searcher.py +++ b/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 diff --git a/tests/v2/test_orf_searcher.py b/tests/v2/test_orf_searcher.py index f6a07ed..d180ab0 100644 --- a/tests/v2/test_orf_searcher.py +++ b/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()