Python reference¶

@dataclass(frozen=True, slots=True)
class SequenceRecord:
    """A DNA sequence with its topology, features, primers and notes.

    Parameters
    ----------
    sequence
        Stored upper-case.
    topology
        ``"linear"`` or ``"circular"``.
    name
        The name shown on a map.
    features, primers
        Every feature and primer binding site must fit the sequence under its topology. On a
        circular record a feature's first segment fits, the rest follow it in ascending order,
        and the whole spans at most one turn.
    notes
        Descriptive fields, such as ``"Description"``, keyed by name.
    extras
        Format-specific data a reader keeps for its writer. Ignored by ``==``.

    Raises
    ------
    ValueError
        If `sequence` holds a letter outside `IUPAC_DNA`, `topology` is unknown, or a span does
        not fit.
    """

    sequence: str
    _: KW_ONLY
    topology: Topology = "linear"
    name: str = ""
    features: tuple[Feature, ...] = ()
    primers: tuple[Primer, ...] = ()
    notes: Mapping[str, str] = field(default_factory=dict, hash=False)
    extras: Mapping[str, object] = field(default_factory=dict, hash=False, compare=False)

    def __post_init__(self) -> None:
        """Upper-case the sequence and check every span fits."""
        if self.topology not in ("linear", "circular"):
            raise ValueError(f"topology must be 'linear' or 'circular', got {self.topology!r}")
        object.__setattr__(self, "sequence", _dna(self.sequence))
        for feature in self.features:
            self._check_feature(feature)
        for primer in self.primers:
            for site in primer.binding_sites:
                self._check_fits(site.start, site.end, f"primer {primer.name!r}")

    def __len__(self) -> int:
        """Return the number of bases."""
        return len(self.sequence)

    def extract(self, span: Feature | Segment) -> str:
        """Return the bases under a feature or segment, reading across the origin.

        A reverse-strand feature reads as the reverse complement of its joined segments. A bare
        segment reads off the top strand. Each segment reads as `bases` reads it.

        Raises
        ------
        ValueError
            If the feature or segment does not fit this record, as an annotation must.

        Examples
        --------
        >>> SequenceRecord("AACCGGTTAC", topology="circular").extract(Segment(8, 12))
        'ACAA'
        """
        if isinstance(span, Feature):
            self._check_feature(span)
        else:
            self._check_fits(span.start, span.end, "extracted")
        segments = span.segments if isinstance(span, Feature) else (span,)
        bases = "".join(self.bases(segment.start, segment.end) for segment in segments)
        if isinstance(span, Feature) and span.strand == Strand.REVERSE:
            return reverse_complement(bases)
        return bases

    def bases(self, start: int, end: int) -> str:
        """Return the top-strand bases in ``[start, end)``, reading across the origin.

        On a circular record a position counts round the circle, so `start` may lie outside the
        record; the span is at most one turn long.

        Raises
        ------
        ValueError
            If `end` comes before `start`, or the span runs off a linear record or is longer than
            a circular one.

        Examples
        --------
        >>> SequenceRecord("AACCGGTTAC", topology="circular").bases(-2, 2)
        'ACAA'
        """
        n = len(self.sequence)
        span = f"span {start}-{end}"
        if end < start:
            raise ValueError(f"{span} ends before it starts")
        if self.topology == "linear":
            if start < 0 or end > n:
                raise ValueError(f"{span} runs past the end of a linear record of {n} bases")
            return self.sequence[start:end]
        if end - start > n:
            raise ValueError(f"{span} is longer than a circular record of {n} bases")
        if start == end:
            return ""
        first = start % n
        last = first + end - start
        if last <= n:
            return self.sequence[first:last]
        return self.sequence[first:] + self.sequence[: last - n]

    def covers(
        self, span: Feature | Segment | BindingSite, inner: int | Segment | BindingSite
    ) -> bool:
        """Whether `span` holds a position, or the whole of a shorter span, across the origin too.

        A feature holds what one of its segments holds. On a circular record a position counts
        round the circle, so one past the length names the base it reaches there.

        Examples
        --------
        >>> SequenceRecord("AACCGGTTAC", topology="circular").covers(Segment(8, 12), 1)
        True
        """
        segments = span.segments if isinstance(span, Feature) else (span,)
        first, last = (inner, inner + 1) if isinstance(inner, int) else (inner.start, inner.end)
        n = len(self.sequence)
        if self.topology == "linear":
            return any(one.start <= first and last <= one.end for one in segments)
        return any(
            (first - one.start) % n + last - first <= one.end - one.start for one in segments
        )

    def fits(self, start: int, end: int) -> bool:
        """Whether ``[start, end)`` lies on this record, as every span it annotates must.

        On a linear record it lies inside the bases. On a circular one it starts inside them and
        is at most one turn long, ending past the length across the origin. An empty span fits
        wherever a base could be inserted.

        Examples
        --------
        >>> SequenceRecord("AACCGGTTAC", topology="circular").fits(8, 12)
        True
        """
        n = len(self.sequence)
        if not 0 <= start <= end:
            return False
        if self.topology == "linear":
            return end <= n
        return start < n and end - start <= n

    def span(self, start: int, width: int) -> Segment:
        """Return the `width` bases at `start`, counted round the origin, clipped on a line.

        A circular record counts `start` round the circle, so a span reaching past the last
        base ends past the record's length. A linear record has nothing outside its bases, so
        a span hanging off either end is cut back to them. `bases` is handed the far end
        instead of a width.

        Raises
        ------
        ValueError
            If nothing is left: a width of no bases, or a linear span wholly off the record.

        Examples
        --------
        >>> SequenceRecord("AACCGGTTAC", topology="circular").span(18, 4)
        Segment(start=8, end=12, name='', color=None)
        >>> SequenceRecord("AACCGGTTAC").span(-2, 6)
        Segment(start=0, end=4, name='', color=None)
        """
        if self.topology == "circular":
            first = start % len(self)
            return Segment(first, first + width)
        return Segment(max(start, 0), min(start + width, len(self)))

    def _check_feature(self, feature: Feature) -> None:
        owner, segments = f"feature {feature.name!r}", feature.segments
        if self.topology == "linear":
            for segment in segments:
                self._check_fits(segment.start, segment.end, owner)
            return
        for before, after in pairwise(segments):
            if after.start < before.start:
                raise ValueError(
                    f"{owner} lists span {after.start}-{after.end} after {before.start}-"
                    f"{before.end}: past the origin of a circular record, a span starts past "
                    "its length"
                )
        self._check_fits(segments[0].start, max(one.end for one in segments), owner)

    def _check_fits(self, start: int, end: int, owner: str) -> None:
        if self.fits(start, end):
            return
        n = len(self.sequence)
        span = f"{owner} span {start}-{end}"
        if self.topology == "linear":
            raise ValueError(f"{span} runs past the end of a linear record of {n} bases")
        raise ValueError(f"{span} does not fit a circular record of {n} bases")

@dataclass(frozen=True, slots=True)
class Feature:
    """A named, typed annotation over one or more segments.

    Parameters
    ----------
    name
        The label shown on a map.
    type
        A GenBank feature key, such as ``"CDS"`` or ``"promoter"``.
    segments
        In the order the feature reads them along the top strand, so a reverse-strand feature
        reads them last to first. On a circular record they start in ascending order: a segment
        after the origin starts past the record's length, as one across it ends past it. A linear
        record has no origin to cross, so a feature cut apart at its two ends lists the higher
        start first.
    strand
        The strand the feature reads along.
    qualifiers
        GenBank-style qualifiers, each holding one or more values.
    color
        ``"#rrggbb"``, or ``None`` for the writer's default.

    Raises
    ------
    ValueError
        If `segments` is empty.
    """

    name: str
    type: str
    segments: tuple[Segment, ...]
    _: KW_ONLY
    strand: Strand = Strand.NONE
    qualifiers: Mapping[str, tuple[str | int, ...]] = field(default_factory=dict, hash=False)
    color: str | None = None

    def __post_init__(self) -> None:
        """Refuse a feature with no segment."""
        if not self.segments:
            raise ValueError(f"feature {self.name!r} has no segment")

    def counted_round(self, length: int) -> "Feature":
        """Return this feature as a circular record of `length` bases holds it.

        Each segment keeps its width and its place in the reading order, starts where it lies in
        the record, and is then counted round past the length if it comes after the origin.

        Examples
        --------
        >>> Feature("f", "CDS", (Segment(91, 98), Segment(2, 8))).counted_round(100).segments[1]
        Segment(start=102, end=108, name='', color=None)
        """
        starts = counted_round((one.start % length for one in self.segments), length)
        return dataclasses.replace(
            self,
            segments=tuple(
                dataclasses.replace(one, start=start, end=start + one.end - one.start)
                for start, one in zip(starts, self.segments, strict=True)
            ),
        )

@dataclass(frozen=True, slots=True)
class Segment:
    """One contiguous span of a feature.

    Parameters
    ----------
    start, end
        0-based, half-open. `end` passes the length of a circular record when the segment
        crosses the origin.
    name
        A label for this segment alone, such as a promoter's ``"-10"``. A ``.dna`` file carries
        it both ways; a GenBank file read here carries it in, and nothing here writes it out.
    color
        ``"#rrggbb"``, or ``None`` to take the feature's colour.

    Raises
    ------
    ValueError
        Unless ``0 <= start < end``.
    """

    start: int
    end: int
    _: KW_ONLY
    name: str = ""
    color: str | None = None

    def __post_init__(self) -> None:
        """Refuse an empty or negative span."""
        _check_span(self.start, self.end)

@dataclass(frozen=True, slots=True)
class Primer:
    """A named oligonucleotide.

    Parameters
    ----------
    name
        The name it is ordered under.
    sequence
        5' to 3', tail included. Stored upper-case.
    binding_sites
        Where it anneals to the record that carries it.
    description
        Free text.

    Raises
    ------
    ValueError
        If `sequence` holds a letter outside `IUPAC_DNA`.
    """

    name: str
    sequence: str
    _: KW_ONLY
    binding_sites: tuple[BindingSite, ...] = ()
    description: str = ""

    def __post_init__(self) -> None:
        """Upper-case and check the sequence."""
        object.__setattr__(self, "sequence", _dna(self.sequence))

class Strand(IntEnum):
    """The strand a feature or primer binding site lies on."""

    FORWARD = 1
    REVERSE = -1
    NONE = 0
    #: SnapGene's bidirectional.
    BOTH = 2

@dataclass(frozen=True, slots=True)
class Check:
    """One verdict on a primer, a pair or an assembled product, with the value it judged.

    Parameters
    ----------
    name
        What was measured, such as ``"gc_clamp"``.
    status
        One of `STATUSES`, or ``None`` where no sourced threshold judges it.
    value
        The measurement, in the unit the module that judged it documents.
    detail
        What a reader needs besides the number.
    """

    name: str
    status: Status | None
    value: float
    detail: str = ""

@dataclass(frozen=True, slots=True)
class Plan:
    """One planned Golden Gate experiment.

    Parameters
    ----------
    vector, inserts
        The records the plan was made from, the inserts in the order they go round the product
        and each on the strand that goes in.
    span
        The vector bases the assembly replaces, after any junction slide.
    choice
        The enzyme the plan uses, and what using it costs.
    ranking
        Every candidate enzyme, best first.
    overhangs
        The overhang every junction takes, with the fidelity of the whole set. An assembly of n
        inserts has n + 1 junctions.
    linearised_vector
        The part the vector is opened into.
    insert_parts
        The part each insert is amplified into, in insert order.
    assembly
        The simulated product, the parts and the junctions.
    colony
        The colony PCR that reads every junction and tells a correct clone from an empty vector.
        No insert can go in the other way round, so there is no lane for one.
    reads
        A sequencing primer reading in from outside the first junction and the last.
    amounts
        What to put in the assembly reaction, vector first.
    phenotype
        What the product says about itself.
    designed_oligos
        Every designed oligo and what it is for, in order: each part's PCR, the colony PCR, the
        reads.
    host, polymerase
        The choices the protocol names.
    cleanup_kit
        The spin-column kit the protocol names, as `mbio.bench.materials.kit` resolves it.
    thresholds
        What those oligos were designed and judged by, for each role, so a page prints the band
        beside the value.
    """

    vector: SequenceRecord
    inserts: tuple[SequenceRecord, ...]
    span: tuple[int, int]
    choice: EnzymeChoice
    ranking: tuple[EnzymeChoice, ...]
    overhangs: OverhangSet
    linearised_vector: Part
    insert_parts: tuple[Part, ...]
    assembly: Assembly
    colony: ColonyCheck
    reads: tuple[SangerRead, SangerRead]
    amounts: tuple[Amount, ...]
    phenotype: Phenotype
    designed_oligos: tuple[DesignedOligo, ...]
    host: str
    polymerase: Polymerase
    cleanup_kit: Material = field(default_factory=kit)
    thresholds: Mapping[PrimerRole, Thresholds] = THRESHOLDS_FOR

    @property
    def enzyme(self) -> Enzyme:
        """The Type IIS enzyme the assembly is cut with."""
        return self.assembly.enzyme

    @property
    def product(self) -> SequenceRecord:
        """The circular plasmid the assembly makes."""
        return self.assembly.product

    @property
    def parts(self) -> tuple[Part, ...]:
        """The parts that go into the reaction, the linearised vector first."""
        return (self.linearised_vector, *self.insert_parts)

    @property
    def oligos(self) -> tuple[Primer, ...]:
        """Every oligo the plan designs, in the order the sheet lists them."""
        return tuple(report.primer for report in self.reports)

    @property
    def reports(self) -> tuple[PrimerReport, ...]:
        """Every designed oligo's evaluation, in the order the sheet lists them."""
        return tuple(oligo.report for oligo in self.designed_oligos)

    @property
    def checks(self) -> tuple[Check, ...]:
        """The product's checks, with one more for the oligos."""
        return (*self.assembly.checks, primer_check(self.reports))

    @property
    def status(self) -> Status:
        """The worst status of any check."""
        return status(self.checks)

    def protocol(self) -> Protocol:
        """Return the bench protocol for this plan."""
        made = protocol_for(
            vector=self.vector,
            span=self.span,
            overhangs=self.overhangs,
            linearised_vector=self.linearised_vector,
            insert_parts=self.insert_parts,
            assembly=self.assembly,
            colony=self.colony,
            reads=self.reads,
            amounts=self.amounts,
            phenotype=self.phenotype,
            oligos=self.designed_oligos,
            checks=self.checks,
            host=self.host,
            polymerase=self.polymerase,
            cleanup_kit=self.cleanup_kit,
            thresholds=self.thresholds,
        )
        return ordered_from_sheet(made)

    def write(self, directory: str | os.PathLike[str]) -> Files:
        """Write the product, the amplicons, the primer sheet, the protocol, its page and maps.

        The directory is made when it is not there. The files are named by `PRODUCT_FILE`,
        `amplicon_files` and `PRIMER_FILE`, and by `mbio.protocol.render` for the protocol
        pair, and a second run over the same inputs writes the same bytes.
        """
        out = Path(directory)
        out.mkdir(parents=True, exist_ok=True)
        product = out / PRODUCT_FILE
        write_dna(self.product, product)
        amplicons = tuple(out / name for name in amplicon_files(one.name for one in self.parts))
        for part, path in zip(self.parts, amplicons, strict=True):
            write_dna(part.amplicon, path)
        sheet = out / PRIMER_FILE
        sheet.write_text(primer_sheet(self.reports), encoding="utf-8")
        written = write_run_files(as_project(self.protocol()), out)
        return Files(product, amplicons, sheet, written.data, written.page, written.maps)