diff --git a/src/lab/experiments/cloning/_dna.py b/src/lab/experiments/cloning/_dna.py new file mode 100644 index 0000000..7885cf0 --- /dev/null +++ b/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) diff --git a/src/lab/experiments/cloning/domestication.py b/src/lab/experiments/cloning/domestication.py new file mode 100644 index 0000000..e051ebc --- /dev/null +++ b/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) diff --git a/src/lab/experiments/cloning/sequences.py b/src/lab/experiments/cloning/sequences.py new file mode 100644 index 0000000..4d733e9 --- /dev/null +++ b/src/lab/experiments/cloning/sequences.py @@ -0,0 +1,188 @@ +"""Calculate products from explicitly selected restriction digest fragments.""" + +from collections.abc import Mapping +from dataclasses import dataclass + +from lab._version import __version__ +from lab.artifacts import digest +from lab.experiments.cloning._dna import DnaSequence, Duplex +from lab.experiments.cloning.systems import AssemblyRecipe, FragmentSelection +from lab.provenance import ( + Activity, + Agent, + AgentKind, + Association, + Component, + DocumentSnapshot, + EvidenceState, + Plan, + Ref, + Sequence, + TopLevel, + Usage, +) +from lab.provenance.vocabulary import DNA, IUPAC_DNA, LAB + +CIRCULAR = "https://identifiers.org/SO:0000988" +LINEAR = "https://identifiers.org/SO:0000987" + + +@dataclass(frozen=True, kw_only=True) +class DigestFragment: + selection: FragmentSelection + watson: str + crick: str + overhang: int + + def record(self) -> Duplex: + return Duplex(self.watson, self.crick, self.overhang) + + +def sequence_record( + component: Ref[Component], + document: DocumentSnapshot, + resolved: Mapping[str, Ref[Component]] | None = None, +) -> DnaSequence: + actual = (resolved or {}).get(component.identity, component) + design = document.get(actual.identity, Component) + sequences = [document.get(ref.identity, Sequence) for ref in design.sequences] + dna = [sequence for sequence in sequences if sequence.encoding == IUPAC_DNA] + if DNA not in design.types or len(dna) != 1: + raise ValueError(f"{design.identity} needs exactly one explicit DNA sequence") + if not dna[0].elements or set(dna[0].elements.upper()) - set("ACGT"): + raise ValueError(f"{design.identity} needs a complete, unambiguous DNA sequence") + if (CIRCULAR in design.types) == (LINEAR in design.types): + raise ValueError( + f"{design.identity} needs exactly one explicit circular or linear topology" + ) + return DnaSequence(dna[0].elements.upper(), circular=CIRCULAR in design.types) + + +def digest_fragments( + component: Ref[Component], + *, + document: DocumentSnapshot, + enzyme: str, + resolved: Mapping[str, Ref[Component]] | None = None, +) -> tuple[DigestFragment, ...]: + record = sequence_record(component, document, resolved) + result: list[DigestFragment] = [] + for left, right, fragment in record.digest(enzyme): + result.append( + DigestFragment( + selection=FragmentSelection(component=component, left_cut=left, right_cut=right), + watson=fragment.watson, + crick=fragment.crick, + overhang=fragment.overhang, + ) + ) + return tuple(result) + + +@dataclass(frozen=True, kw_only=True) +class AssemblyDesign: + product: Component + sequence: Sequence + activity: Activity + agent: Agent + plan: Plan + + @property + def objects(self) -> tuple[TopLevel, ...]: + return self.product, self.sequence, self.activity, self.agent, self.plan + + +def calculate_assembly( + recipe: AssemblyRecipe, + *, + document: DocumentSnapshot, + resolved: Mapping[str, Ref[Component]] | None = None, +) -> AssemblyDesign: + selected: list[Duplex] = [] + inputs: dict[str, Ref[Component]] = {} + for selection in recipe.fragments: + fragments = digest_fragments( + selection.component, document=document, enzyme=recipe.enzyme, resolved=resolved + ) + matches = [ + fragment + for fragment in fragments + if fragment.selection.left_cut == selection.left_cut + and fragment.selection.right_cut == selection.right_cut + ] + if len(matches) != 1: + raise ValueError( + f"{selection.component.identity}: cut pair " + f"({selection.left_cut}, {selection.right_cut}) " + "does not select exactly one digest fragment" + ) + fragment = matches[0].record() + selected.append(fragment.reverse_complement() if selection.reverse_complement else fragment) + actual = (resolved or {}).get(selection.component.identity, selection.component) + inputs[actual.identity] = actual + product = selected[0] + try: + for fragment in selected[1:]: + product = product.ligate(fragment) + elements = product.close() if recipe.circular else product.linear_sequence() + except (ValueError, TypeError) as error: + raise ValueError( + f"{recipe.identity}: selected fragment ends are incompatible: {error}" + ) from error + try: + intended = document.get(recipe.product.identity, Component) + except KeyError: + intended = None + if intended is not None: + inputs[intended.identity] = intended.ref + if intended.sequences: + expected = sequence_record(intended.ref, document) + equal = len(expected.elements) == len(elements) and ( + elements in expected.elements * 2 + if recipe.circular + else elements == expected.elements + ) + if expected.circular != recipe.circular or not equal: + raise ValueError( + f"{recipe.identity}: calculated sequence does not match the requested design" + ) + calculation_id = ( + recipe.identity + "/calculation_" + digest((recipe, tuple(inputs), tuple(selected)))[:16] + ) + agent = Agent( + identity=calculation_id + "/calculator", + kind=AgentKind.SOFTWARE, + name="Lab sequence calculation", + software_version=f"lab-compiler {__version__}; Biopython 1.84", + ) + plan = Plan( + identity=calculation_id + "/method", + description=( + f"Ordered restriction-fragment ligation using {recipe.enzyme}; " + f"circular={recipe.circular}." + ), + ) + activity = Activity( + identity=calculation_id, + types=(LAB + "sequenceCalculation",), + evidence_state=EvidenceState.RECORDED, + usage=tuple(Usage(entity=ref) for ref in inputs.values()), + association=(Association(agent=agent.ref, plan=plan.ref),), + ) + sequence = Sequence( + identity=calculation_id + "/sequence", + elements=elements, + encoding=IUPAC_DNA, + generated_by=(activity.ref,), + ) + design = Component( + identity=calculation_id + "/product", + name=None if intended is None else intended.name, + types=(DNA, CIRCULAR if recipe.circular else LINEAR), + sequences=(sequence.ref,), + derived_from=tuple(inputs.values()), + generated_by=(activity.ref,), + ) + return AssemblyDesign( + product=design, sequence=sequence, activity=activity, agent=agent, plan=plan + ) diff --git a/tests/test_dna.py b/tests/test_dna.py new file mode 100644 index 0000000..c794568 --- /dev/null +++ b/tests/test_dna.py @@ -0,0 +1,98 @@ +"""Independent expected strands captured from pydna 5.5.0 / Biopython 1.84.""" + +import pytest + +from lab.experiments.cloning._dna import DnaSequence + + +@pytest.mark.parametrize( + "enzyme,elements,circular,expected", + [ + ( + "PstI", + "GGCTGCAGAAAAGCTGCAGTT", + False, + ( + (None, 7, "GGCTGCA", "GCC", 0), + (7, 18, "GAAAAGCTGCA", "GCTTTTCTGCA", 4), + (18, None, "GTT", "AACTGCA", 4), + ), + ), + ( + "SmaI", + "GGCCCGGGTTTCCCGGGAA", + False, + ( + (None, 5, "GGCCC", "GGGCC", 0), + (5, 14, "GGGTTTCCC", "GGGAAACCC", 0), + (14, None, "GGGAA", "TTCCC", 0), + ), + ), + ( + "EcoRI", + "AATTCCCCGGAATTCG", + True, + ( + (0, 10, "AATTCCCCGG", "AATTCCGGGG", -4), + (10, 0, "AATTCG", "AATTCG", -4), + ), + ), + ( + "PstI", + "CTGCAGAAAAGCTGCAGTT", + True, + ( + (5, 16, "GAAAAGCTGCA", "GCTTTTCTGCA", 4), + (16, 5, "GTTCTGCA", "GAACTGCA", 4), + ), + ), + ( + "BsaI", + "TTTGGTCTCAACGTTACGTTTGAGACCAAA", + False, + ( + (None, 10, "TTTGGTCTCA", "ACGTTGAGACCAAA", 0), + (10, 16, "ACGTTA", "AACGTA", -4), + (16, None, "CGTTTGAGACCAAA", "TTTGGTCTCA", -4), + ), + ), + ], +) +def test_digest_and_ligation_match_independent_strand_reference( + enzyme, elements, circular, expected +): + fragments = DnaSequence(elements, circular).digest(enzyme) + assert ( + tuple( + (left, right, fragment.watson, fragment.crick, fragment.overhang) + for left, right, fragment in fragments + ) + == expected + ) + joined = fragments[0][2] + for _, _, fragment in fragments[1:]: + joined = joined.ligate(fragment) + if circular: + assert len(joined.close()) == len(elements) and joined.close() in elements * 2 + else: + assert joined.linear_sequence() == elements + + +def test_reverse_complement_preserves_three_prime_ends(): + fragment = DnaSequence("GGCTGCAGAAAAGCTGCAGTT", False).digest("PstI")[0][2] + reversed_fragment = fragment.reverse_complement() + assert (reversed_fragment.watson, reversed_fragment.crick, reversed_fragment.overhang) == ( + "GCC", + "GGCTGCA", + 4, + ) + assert reversed_fragment.reverse_complement() == fragment + + +def test_incompatible_ends_and_unpaired_linear_products_fail_explicitly(): + eco = DnaSequence("GGGAATTCCCCGAATTCGG", False).digest("EcoRI")[1][2] + pst = DnaSequence("GGCTGCAGAAAAGCTGCAGTT", False).digest("PstI")[1][2] + with pytest.raises(ValueError, match="Incompatible ends"): + eco.ligate(pst) + with pytest.raises(ValueError, match="end repair"): + eco.linear_sequence()