Bio.Restriction

Core Principle: Bio.Restriction incorporates the REBASE restriction enzyme database to identify cleavage sites, simulate in silico digestions, and evaluate restriction fragment distributions.

A small module, but a good example of Biopython doing something genuinely hard to reimplement: it encodes recognition sites, cut positions, star activity, methylation sensitivity and isoschizomer relationships for thousands of enzymes.

Basics

4096.0

search() returns 1-based cut positions, not the start of the recognition site. Mixing these up shifts your fragment boundaries. See Coordinate Systems.

Multiple enzymes

BamHI      :  8.
EcoRI      :  2, 20.
HindIII    :  14.

{BamHI: [8], HindIII: [14]}

linear=True vs linear=False matters. A plasmid is circular, so a site spanning the origin is real and a linear analysis will miss it. Get this wrong on a vector and your predicted digest is wrong.

Enzyme collections

CommOnly is the practical default. There is no point designing a cloning strategy around an enzyme nobody sells.

Cloning workflow

# 1. which unique cutters are available in my insert's MCS?
ana = Analysis(CommOnly, insert, linear=True)
unique = ana.with_N_sites(1)

# 2. which of those do NOT cut inside my gene of interest?
gene_ana = Analysis(RestrictionBatch(list(unique)), gene, linear=True)
safe = gene_ana.without_site()

# 3. and are also present in the vector's MCS?
vec_ana = Analysis(RestrictionBatch(safe), vector, linear=False)   # circular!
usable = vec_ana.with_N_sites(1)

print(sorted(str(e) for e in usable))

This three-step intersection: cuts the vector once, cuts the insert flanks, does not cut the gene body: is the actual logic of restriction cloning, and it is much less error-prone done programmatically than by eye.

Compatible ends and isoschizomers

'GATC'

Compatible ends let you ligate fragments cut with different enzymes: BamHI and BglII both leave a GATC 5' overhang, so their products ligate (destroying both sites, which can be exactly what you want).

Methylation sensitivity

EcoRI.is_methylable()
DpnI.site                    # 'GATC': cuts ONLY Dam-methylated DNA

DpnI is the workhorse of site-directed mutagenesis: it digests the methylated template plasmid (grown in dam+ E. coli) while leaving the unmethylated PCR product intact.

Conversely, many enzymes are blocked by CpG methylation. A digest that works on PCR product and fails on genomic DNA is usually a methylation problem, and is_methylable() will tell you it was foreseeable.

Digest simulation and gel prediction

fragments = EcoRI.catalyse(plasmid, linear=False)
sizes = sorted((len(f) for f in fragments), reverse=True)
print(sizes)

rb = RestrictionBatch([EcoRI, BamHI])
frags = rb.search(plasmid, linear=False)

# predict a gel
fig, ax = plt.subplots(figsize=(2.5, 5))
for s in sizes:
    ax.axhline(s, xmin=0.2, xmax=0.8, lw=3, c="black")
ax.set_yscale("log")
ax.set_ylabel("Fragment size (bp)")
ax.set_xticks([])
ax.invert_yaxis()

Comparing a predicted gel to an observed one is the standard construct-verification step, and getting the prediction from code rather than by hand removes a whole class of arithmetic errors.

Bioinformatics examples

# RFLP: does a SNP create or destroy a restriction site?
ref_ctx = Seq(genome[chrom].seq[pos-10:pos+10])
alt_ctx = Seq(str(ref_ctx)[:10] + alt + str(ref_ctx)[11:])

ref_sites = Analysis(CommOnly, ref_ctx, linear=True).with_sites()
alt_sites = Analysis(CommOnly, alt_ctx, linear=True).with_sites()

created   = set(alt_sites) - set(ref_sites)
destroyed = set(ref_sites) - set(alt_sites)
print(f"created: {[str(e) for e in created]}")
print(f"destroyed: {[str(e) for e in destroyed]}")

That is a genuine assay-design task: an RFLP genotyping assay needs a variant that changes a site, and this finds candidates in one pass.

# how many fragments would a RAD-seq enzyme produce genome-wide?
for enzyme in [EcoRI, PstI, SbfI]:
    n = sum(len(enzyme.search(rec.seq)) for rec in SeqIO.parse("genome.fa", "fasta"))
    print(f"{enzyme}: {n:,} sites → ~{n * 2:,} tags")

# check a synthetic construct for unwanted sites
banned = RestrictionBatch([EcoRI, BamHI, XhoI, NotI])
hits = banned.search(construct, linear=True)
for enz, positions in hits.items():
    if positions:
        print(f"WARNING: {enz} cuts at {positions}")

# Golden Gate: check for internal BsaI sites
from Bio.Restriction import BsaI
if BsaI.search(part_seq, linear=True):
    print("Domesticate this part: internal BsaI site")

The Golden Gate check is a real requirement: type IIS assembly fails if any part contains an internal recognition site, and "domestication" (silently mutating it out) is a standard step.

Common mistakes

• search() positions read as site starts rather than cut positions.
• linear=True on a plasmid. Misses sites spanning the origin.
• Using AllEnzymes and choosing something unpurchasable.
• Ignoring methylation sensitivity, then debugging a failed genomic digest.
• Ignoring star activity: many enzymes cut degenerate sites under suboptimal conditions.
• Forgetting isoschizomers. The enzyme you want may be sold under a different name.
• Not checking the insert body, only the flanks.
• Assuming a single cut site without verifying with with_N_sites(1).

See also

Seq · SeqIO · Coordinate Systems · Bio.SeqUtils · Applied - Sequence Composition Analysis

Exercise

Bio.Restriction is not available here, so find the sites with str.find and work out the fragments yourself. Each enzyme is given as a recognition site with a caret marking the cut, for example G^AATTC: removing the caret gives the site, and the caret's index gives how many bases of the site lie before the cut, so a match beginning at index i cuts at i + offset. Report cut positions 0-based, unlike the library's 1-based search().

Write digest(topology, seq, enzymes), where topology is linear or circular and enzymes is a list of (name, spec) pairs. Print one line per enzyme in ascending name order as <name> cuts=<n> at=<comma-joined positions> with the positions sorted ascending, then a final fragments=<comma-joined sizes> line for the combined digest with the sizes sorted descending. On a circular sequence a site may span the origin, so search seq + seq[:len(site) - 1] and reduce each cut position modulo the length; a circular sequence with k cuts yields k fragments, and with exactly one cut that single fragment is the whole molecule rather than a fragment of length zero. On a linear sequence the cuts split it into k + 1 pieces; discard any piece of length zero. An enzyme that does not cut prints cuts=0 at= and contributes nothing. The topology is the first line of input, the sequence is the second, and the remaining lines are name and spec separated by a tab.