[7/11] Calculate cloning sequences and explicit domestication proposals

src/lab/experiments/cloning/_dna.py

@@ -0,0 +1,118 @@
+"""Small, explicit double-stranded restriction/ligation calculations.
+
+Both strands are written 5' to 3'. ``overhang`` is Watson-start minus
+Crick-start when both strands are aligned left to right. Enzyme recognition and
+Watson cleavage positions come from Biopython. No end repair is implicit.
+"""
+
+from dataclasses import dataclass
+
+from Bio.Restriction.Restriction import RestrictionBatch
+from Bio.Seq import Seq
+
+
+def reverse_complement(sequence: str) -> str:
+    return str(Seq(sequence).reverse_complement())
+
+
+@dataclass(frozen=True)
+class Duplex:
+    watson: str
+    crick: str
+    overhang: int = 0
+
+    @property
+    def left_end(self) -> tuple[str, str]:
+        if self.overhang < 0:
+            return "5", self.watson[: -self.overhang]
+        if self.overhang > 0:
+            return "3", self.crick[-self.overhang :]
+        return "blunt", ""
+
+    @property
+    def right_end(self) -> tuple[str, str]:
+        difference = len(self.watson) - len(self.crick) + self.overhang
+        if difference > 0:
+            return "3", self.watson[-difference:]
+        if difference < 0:
+            return "5", self.crick[:-difference]
+        return "blunt", ""
+
+    def reverse_complement(self) -> "Duplex":
+        return Duplex(self.crick, self.watson, self.overhang + len(self.watson) - len(self.crick))
+
+    def ligate(self, other: "Duplex") -> "Duplex":
+        if not compatible(self.right_end, other.left_end):
+            raise ValueError(f"Incompatible ends: {self.right_end} and {other.left_end}")
+        return Duplex(self.watson + other.watson, other.crick + self.crick, self.overhang)
+
+    def close(self) -> str:
+        if not compatible(self.right_end, self.left_end) or len(self.watson) != len(self.crick):
+            raise ValueError("Fragment ends cannot close into a circular duplex")
+        return self.watson
+
+    def linear_sequence(self) -> str:
+        if self.left_end[0] != "blunt" or self.right_end[0] != "blunt":
+            raise ValueError("Linear products with unpaired ends need explicit end repair")
+        return self.watson
+
+
+def compatible(right: tuple[str, str], left: tuple[str, str]) -> bool:
+    return right[0] == left[0] and reverse_complement(right[1]) == left[1]
+
+
+@dataclass(frozen=True)
+class DnaSequence:
+    elements: str
+    circular: bool
+
+    def cuts(self, enzyme: str) -> tuple[int, ...]:
+        restriction = next(iter(RestrictionBatch([enzyme])))
+        positions = {
+            position - 1
+            for position in restriction.search(Seq(self.elements), linear=not self.circular)
+        }
+        if self.circular:
+            return tuple(sorted(position % len(self.elements) for position in positions))
+        return tuple(
+            sorted(
+                position
+                for position in positions
+                if 0 <= position <= len(self.elements)
+                and 0 <= position - restriction.ovhg <= len(self.elements)
+            )
+        )
+
+    def digest(self, enzyme: str) -> tuple[tuple[int | None, int | None, Duplex], ...]:
+        restriction = next(iter(RestrictionBatch([enzyme])))
+        # Enzymes with two cleavage events or unknown cuts need a different model.
+        if restriction.ovhg is None or restriction.scd5 is not None or restriction.scd3 is not None:
+            raise ValueError(f"{enzyme} does not have one supported, known cleavage pair")
+        cuts = self.cuts(enzyme)
+        if not cuts:
+            return ()
+        bounds: tuple[int | None, ...] = (*cuts, cuts[0]) if self.circular else (None, *cuts, None)
+        result: list[tuple[int | None, int | None, Duplex]] = []
+        length = len(self.elements)
+
+        def region(start: int, end: int) -> str:
+            if self.circular:
+                return "".join(self.elements[index % length] for index in range(start, end))
+            return self.elements[start:end]
+
+        for left, right in zip(bounds, bounds[1:], strict=False):
+            watson_start = 0 if left is None else left
+            watson_end = length if right is None else right
+            if self.circular and watson_end <= watson_start:
+                watson_end += length
+            crick_start = 0 if left is None else watson_start - restriction.ovhg
+            crick_end = length if right is None else watson_end - restriction.ovhg
+            if min(watson_end, crick_end) <= max(watson_start, crick_start):
+                raise ValueError("Overlapping cuts do not leave a double-stranded fragment")
+            fragment = Duplex(
+                region(watson_start, watson_end),
+                reverse_complement(region(crick_start, crick_end)),
+                watson_start - crick_start,
+            )
+            result.append((left, right, fragment))
+        return tuple(result)

src/lab/experiments/cloning/domestication.py

@@ -0,0 +1,144 @@
+"""Sequence edit proposals that require an explicit caller selection."""
+
+from dataclasses import dataclass, replace
+
+from lab.experiments.cloning._dna import DnaSequence
+from lab.experiments.cloning.sequences import sequence_record
+from lab.provenance import (
+    Activity,
+    Component,
+    Document,
+    DocumentSnapshot,
+    EvidenceState,
+    Ref,
+    Sequence,
+    Usage,
+)
+from lab.provenance.types import require_iri
+from lab.provenance.vocabulary import IUPAC_DNA, LAB
+
+
+@dataclass(frozen=True, kw_only=True)
+class SequenceEdit:
+    position: int
+    before: str
+    after: str
+
+    def __post_init__(self) -> None:
+        if type(self.position) is not int or self.position < 0:
+            raise ValueError("Edit positions must be nonnegative integers")
+        if (
+            self.before not in tuple("ACGT")
+            or self.after not in tuple("ACGT")
+            or self.before == self.after
+        ):
+            raise ValueError("An edit replaces one unambiguous base with a different base")
+
+
+@dataclass(frozen=True, kw_only=True)
+class EditProposal:
+    """A single-base candidate; no claim is made about biological function.
+
+    Positions are zero-based. Further sites can remain after this edit. Applying
+    a proposal creates a new design; it does not alter inventory or a recipe.
+    """
+
+    component: Ref[Component]
+    sequence: Ref[Sequence]
+    enzyme: str
+    edit: SequenceEdit
+    remaining_sites: int
+
+    def apply(self, document: DocumentSnapshot, *, identity: str) -> DocumentSnapshot:
+        require_iri(identity)
+        original = document.get(self.component.identity, Component)
+        sequence = document.get(self.sequence.identity, Sequence)
+        if sequence.ref not in original.sequences:
+            raise ValueError("The proposal sequence does not belong to its component")
+        if (
+            sequence.elements[self.edit.position : self.edit.position + 1].upper()
+            != self.edit.before
+        ):
+            raise ValueError("Edit proposal no longer matches the input sequence")
+        if identity == original.identity:
+            raise ValueError("An edited design needs a new identity")
+        elements = (
+            sequence.elements[: self.edit.position]
+            + self.edit.after
+            + sequence.elements[self.edit.position + 1 :]
+        )
+        activity = Activity(
+            identity=identity + "/edit",
+            types=(LAB + "sequenceEdit",),
+            usage=(Usage(entity=original.ref), Usage(entity=sequence.ref)),
+            evidence_state=EvidenceState.RECORDED,
+        )
+        edited_sequence = replace(
+            sequence,
+            identity=identity + "/sequence",
+            namespace=None,
+            elements=elements,
+            derived_from=(sequence.ref,),
+            generated_by=(activity.ref,),
+        )
+        # Coordinate-bearing annotations require review after an edit. Keep the
+        # original linked rather than copying annotations onto changed bases.
+        edited = Component(
+            identity=identity,
+            types=original.types,
+            roles=original.roles,
+            name=original.name,
+            sequences=(edited_sequence.ref,),
+            derived_from=(original.ref,),
+            generated_by=(activity.ref,),
+        )
+        result = Document.from_snapshot(document)
+        result.add(activity, edited_sequence, edited)
+        return result.freeze()
+
+
+def propose_edits(
+    component: Ref[Component],
+    *,
+    document: DocumentSnapshot,
+    enzyme: str,
+    editable_positions: tuple[int, ...],
+) -> tuple[EditProposal, ...]:
+    """Enumerate substitutions only at positions explicitly declared editable.
+
+    Retain candidates that reduce this enzyme's cut count. The caller must review
+    coding/regulatory effects and other assembly-system constraints before use.
+    """
+    record = sequence_record(component, document)
+    design = document.get(component.identity, Component)
+    sequence = next(
+        sequence
+        for ref in design.sequences
+        if (sequence := document.get(ref.identity, Sequence)).encoding == IUPAC_DNA
+    )
+    baseline = len(record.cuts(enzyme))
+    if not isinstance(editable_positions, tuple) or len(set(editable_positions)) != len(
+        editable_positions
+    ):
+        raise ValueError("Editable positions must be a tuple of unique positions")
+    result: list[EditProposal] = []
+    for position in sorted(editable_positions):
+        if type(position) is not int or not 0 <= position < len(sequence.elements):
+            raise ValueError("Editable position lies outside the sequence")
+        before = sequence.elements[position].upper()
+        for after in "ACGT":
+            if after == before:
+                continue
+            candidate = record.elements[:position] + after + record.elements[position + 1 :]
+            remaining = len(DnaSequence(candidate, record.circular).cuts(enzyme))
+            if remaining < baseline:
+                result.append(
+                    EditProposal(
+                        component=component,
+                        sequence=sequence.ref,
+                        enzyme=enzyme,
+                        edit=SequenceEdit(position=position, before=before, after=after),
+                        remaining_sites=remaining,
+                    )
+                )
+    return tuple(result)