bio-applied-genetic-engineering-in-silico

name	bio-applied-genetic-engineering-in-silico
description	Simulate restriction digests, overhang compatibility, and primer Tm (Wallace/SantaLucia NN) in Python; plot agarose gel bands. Use when planning cloning, enzyme compatibility, or PCR primer design for a target Tm.
tool_type	python
primary_tool	re

Genetic Engineering In Silico

When to Use

• Planning a subcloning strategy and need to pick restriction enzymes that cut a vector/insert at the right places
• Checking whether two enzymes (e.g. SalI/XhoI, XbaI/SpeI) produce compatible sticky ends for ligation
• Designing PCR/cloning primers with a target Tm and appending a restriction site + clamp
• Predicting fragment sizes from a single or double digest before running a real gel
• QC-ing a primer pair for 3' homopolymer runs, hairpins, or primer-dimer risk

Version Compatibility

Pure-Python + stdlib (re, math, itertools) — no version sensitivity. Gel simulation uses matplotlib ≥3.7. Python ≥3.9 (uses list[int] style hints). Tm models correspond to Wallace 1979 rule and SantaLucia 1998 unified nearest-neighbor parameters; for production primer design cross-check against Primer3 (bio-primer-design-primer-basics).

Prerequisites

• pip install matplotlib (only needed for gel plotting; digestion/Tm logic is stdlib-only)
• Familiarity with IUPAC ambiguity codes and 5'/3' sequence orientation
• Related skill: bio-sequence-manipulation-reverse-complement for strand math

Restriction Digestion & Compatibility

Goal: find cut sites for one or two enzymes, split a DNA sequence into fragments (linear or circular), and determine whether two enzymes leave ligatable ends. Approach: build a IUPAC-aware regex per enzyme, find all site starts, offset by the enzyme's top-strand cut position, then slice between consecutive cuts (wrapping around for circular DNA). Overhang compatibility compares the single-stranded sequence left by each enzyme, not the recognition site.

import re
import itertools

# name: (recognition_seq, cut_top, cut_bottom) — 0-indexed from site start,
# cut_bottom counted on the bottom strand read 5'->3'
RESTRICTION_ENZYMES = {
    'EcoRI':   ('GAATTC', 1, 5),   # 5'-AATT overhang
    'BamHI':   ('GGATCC', 1, 5),   # 5'-GATC
    'HindIII': ('AAGCTT', 1, 5),   # 5'-AGCT
    'SalI':    ('GTCGAC', 1, 5),   # 5'-TCGA (compatible with XhoI -> ligation scar)
    'XhoI':    ('CTCGAG', 1, 5),   # 5'-TCGA
    'XbaI':    ('TCTAGA', 1, 5),   # 5'-CTAG (compatible with SpeI -> scar)
    'SpeI':    ('ACTAGT', 1, 5),   # 5'-CTAG
    'NotI':    ('GCGGCCGC', 2, 6), # 5'-GGCC (8-cutter, rare site)
    'SmaI':    ('CCCGGG', 3, 3),   # blunt
    'XmaI':    ('CCCGGG', 1, ),   
    :    (, , ),   
    :    (, , ),   
    :    (, , ),   
}

IUPAC = {: , : , : , : , : ,
         : , : , : , : ,
         : , : , : , : , : , : }


 () -> :
    
     .join(IUPAC.get(b, b)  b  seq.upper())


 () -> :
    comp = {: , : , : , : , : ,
            : , : , : , : , : , : }
     .join(comp.get(b, )  b  (seq.upper()))


 () -> :
    
    site, cut_top, _ = RESTRICTION_ENZYMES[enzyme_name]
    pattern = iupac_to_regex(site)
     (m.start() + cut_top  m  re.finditer(, dna.upper()))


 () -> :
    
    cuts = find_cut_positions(dna, enzyme_name)
      cuts:
         [dna]
     circular:
        first = cuts[]
        rotated = dna[first:] + dna[:first]
        adjusted = [c - first  c  cuts[:]] + [(dna)]
        frags, prev = [], 
         c  adjusted:
            frags.append(rotated[prev:c])
            prev = c
         frags
    boundaries = [] + cuts + [(dna)]
     [dna[boundaries[i]:boundaries[i + ]]  i  ((boundaries) - )]


 () -> :
    
    all_cuts = ((find_cut_positions(dna, enzyme1) + find_cut_positions(dna, enzyme2)))
      all_cuts:
         [dna]
     circular:
         digest(dna, enzyme1, circular=)  enzyme1 == enzyme2  _slice_circular(dna, all_cuts)
    boundaries = [] + all_cuts + [(dna)]
     [dna[boundaries[i]:boundaries[i + ]]  i  ((boundaries) - )]


 () -> :
    first = cuts[]
    rotated = dna[first:] + dna[:first]
    adjusted = [c - first  c  cuts[:]] + [(dna)]
    frags, prev = [], 
     c  adjusted:
        frags.append(rotated[prev:c])
        prev = c
     frags


 () -> :
    
    site, cut_top, cut_bot = RESTRICTION_ENZYMES[enzyme_name]
     cut_top == cut_bot:
         (, )
     cut_top < cut_bot:
         (site[cut_top:cut_bot], )
     (reverse_complement(site)[cut_bot:cut_top], )


 () -> :
    
    oh1, type1 = compute_overhang(enzyme1)
    oh2, type2 = compute_overhang(enzyme2)
     type1 != type2:
         
       type1 ==   oh1 == oh2


 __name__ == :
    plasmid = 
    frags = digest(plasmid, , circular=)
     ((f)  f  frags) == (plasmid)
     are_compatible(, )     
     are_compatible(, )  
    ()

Primer Design and Tm Models

Goal: grow a primer from a template position to hit a target Tm, then QC it and append a restriction site for cloning. Approach: three Tm estimators of increasing accuracy — 4+2 rule (short/quick), salt-adjusted Wallace, and SantaLucia 1998 nearest-neighbor thermodynamics (most accurate for 18–30 bp). Use NN for real ordering decisions.

import math

NN_PARAMS = {  # SantaLucia 1998 (dH kcal/mol, dS cal/mol/K)
    'AA': (-7.9, -22.2), 'AT': (-7.2, -20.4), 'TA': (-7.2, -21.3), 'CA': (-8.5, -22.7),
    'GT': (-8.4, -22.4), 'CT': (-7.8, -21.0), 'GA': (-8.2, -22.2), 'CG': (-10.6, -27.2),
    'GC': (-9.8, -24.4), 'GG': (-8.0, -19.9), 'AC': (-7.8, -21.0), 'TC': (-8.2, -22.2),
    'AG': (-7.8, -21.0), 'TG': (-8.5, -22.7), 'TT': (-7.9, -22.2), 'CC': (-8.0, -19.9),
}
NN_INIT_GC = (0.1, -2.8)
NN_INIT_AT = (2.3, 4.1)


def gc_content() -> :
    s = seq.upper()
     (s.count() + s.count()) / (s)


 () -> :
    
    s = primer.upper()
      * (s.count() + s.count()) +  * (s.count() + s.count())


 () -> :
    
    n = (primer)
    gc = gc_content(primer)
      +  * math.log10(salt_mm / ) +  * gc -  / n


 () -> :
    
    seq, R = primer.upper(),   
    dH = dS = 
     end_base  (seq[], seq[-]):
        h, s = NN_INIT_GC  end_base    NN_INIT_AT
        dH += h
        dS += s
     i  ((seq) - ):
        h, s = NN_PARAMS.get(seq[i:i + ], (-, -))
        dH += h
        dS += s
    dS +=  * ((seq) - ) * math.log(salt_mm / )  
    ct = dna_conc_nm * 
     (dH * ) / (dS + R * math.log(ct / )) - 


 () -> :
    comp = {: , : , : , : , : }
     .join(comp.get(b, )  b  (seq.upper()))


 () -> :
    
    best = 
     length  (min_len, max_len + ):
         direction == :
            seq = template[start:start + length]
        :
            seq = reverse_complement(template[start - length + :start + ])
         (seq) < length:
            
        tm = tm_fn(seq)
        candidate = {: seq, : length, : tm, : gc_content(seq)}
         best    (tm - target_tm) < (best[] - target_tm):
            best = candidate
         tm >= target_tm:
            
     best


 () -> :
    
     ((primer[-run_len:].upper())) == 


 () -> :
    
    tail1 = primer1[-check_len:].upper()
    tail2 = reverse_complement(primer2[-check_len:]).upper()
     (a == b  a, b  (tail1, tail2))


 () -> :
    
     clamp + site + primer


 __name__ == :
    mcs = 
    fwd = design_primer(mcs, start=, direction=, target_tm=)
    rev = design_primer(mcs, start=(mcs) - , direction=, target_tm=)
      <= fwd[] <=   (fwd[] - ) < 
      check_3prime_run(fwd[], run_len=)
    cloning_fwd = add_re_site_to_primer(fwd[], )  
     cloning_fwd.startswith()
    ()

Gel Electrophoresis Simulation

Goal: visualize expected fragment sizes from digests as a synthetic agarose gel. Approach: migration distance is proportional to log10(size); plot each fragment as a horizontal band at y = log10(size) per lane, alongside a DNA ladder.

import math
import matplotlib.pyplot as plt


def plot_gel(lanes: dict, ladder_sizes=None):
    """Render a simulated agarose gel. `lanes` maps a lane label to a list of fragment sequences (or sizes)."""
    if ladder_sizes is None:
        ladder_sizes = [10000, 8000, 6000, 5000, 4000, 3000, 2000, 1500, 1000, 750, 500, 250, 100]
    all_lanes = {'Ladder': [str(s) for s in ladder_sizes]}
    all_lanes.update(lanes)

    fig, ax = plt.subplots(figsize=(2.5 * len(all_lanes), 5))
    ax.set_facecolor('#f5f0e8')
    ax.set_ylim(1.8, 4.1)
    ax.set_xlim(0, len(all_lanes))
    ax.set_xticks([])
    ax.set_yticks([])

    for lane_idx, (label, frags) in enumerate(all_lanes.items()):
        x = lane_idx + 0.5
        sizes = [int(f) if isinstance(f, str) and f.isdigit() else len(f) for f in frags]
         size  sizes:
             size < :
                
            y = math.log10(size)
            color =   label ==   
            ax.plot([x - , x + ], [y, y], color=color, linewidth=, solid_capstyle=)
        ax.text(x, , label, ha=, va=, fontsize=, fontweight=)
    ax.set_title()
    plt.tight_layout()
     fig

Pitfalls

• Compatible ends create scars: SalI/XhoI and XbaI/SpeI are ligation-compatible but leave a hybrid site — check the scar sequence if the junction falls in a coding region.
• Isoschizomers vs neoschizomers: SmaI and XmaI recognize the same site but cut differently (blunt vs 5'-CCGG overhang) — verify cut position, not just recognition sequence.
• Nearest-neighbor Tm is a solution-phase estimate: real PCR annealing temp is typically 3-5°C below calculated Tm; optimize with gradient PCR rather than trusting Tm alone.
• IUPAC degenerate bases: enzyme databases use IUPAC codes (e.g. AvaI = CYCGRG) — always translate via iupac_to_regex before pattern search, or matches will silently fail.
• Circular vs linear digest: for plasmids, n cuts → n fragments; for linear DNA, n cuts → n+1 fragments. Passing circular=False on a plasmid undercounts fragments by one.
• Primer 3' stability: design_primer optimizes for Tm only — always separately check check_3prime_run and count_3prime_complementarity before ordering.

See Also

• bio-restriction-analysis-restriction-mapping — building full restriction maps and multi-enzyme digest tables
• bio-restriction-analysis-enzyme-selection — choosing enzymes absent from a sequence (silent cloning sites)
• bio-primer-design-primer-validation — deeper primer QC (dimers, hairpins, specificity via BLAST)
• bio-primer-design-qpcr-primers — Tm/amplicon design tuned for qPCR assays