Bio applied genetic engineering in silico

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, 5),   # isoschizomer of SmaI, different cut -> 5'-CCGG
    'KpnI':    ('GGTACC', 5, 1),   # 3'-GTAC overhang
    'PstI':    ('CTGCAG', 5, 1),   # 3'-TGCA overhang
    'AvaI':    ('CYCGRG', 1, 5),   # degenerate site: C[CT]CG[AG]G
}

IUPAC = {'A': 'A', 'T': 'T', 'G': 'G', 'C': 'C', 'N': '[ATGC]',
         'R': '[AG]', 'Y': '[CT]', 'W': '[AT]', 'S': '[GC]',
         'M': '[AC]', 'K': '[GT]', 'B': '[CGT]', 'D': '[AGT]', 'H': '[ACT]', 'V': '[ACG]'}


def iupac_to_regex(seq: str) -> str:
    """Translate an IUPAC-degenerate recognition sequence into a regex pattern."""
    return ''.join(IUPAC.get(b, b) for b in seq.upper())


def reverse_complement(seq: str) -> str:
    comp = {'A': 'T', 'T': 'A', 'G': 'C', 'C': 'G', 'N': 'N',
            'R': 'Y', 'Y': 'R', 'W': 'W', 'S': 'S', 'M': 'K', 'K': 'M'}
    return ''.join(comp.get(b, 'N') for b in reversed(seq.upper()))


def find_cut_positions(dna: str, enzyme_name: str) -> list:
    """Return sorted top-strand cut positions (0-indexed) for all sites of `enzyme_name` in `dna`."""
    site, cut_top, _ = RESTRICTION_ENZYMES[enzyme_name]
    pattern = iupac_to_regex(site)
    return sorted(m.start() + cut_top for m in re.finditer(f'(?={pattern})', dna.upper()))


def digest(dna: str, enzyme_name: str, circular: bool = False) -> list:
    """Return fragment sequences after single-enzyme digestion.
    Linear DNA with n cuts -> n+1 fragments; circular DNA with n cuts -> n fragments.
    """
    cuts = find_cut_positions(dna, enzyme_name)
    if not cuts:
        return [dna]
    if circular:
        first = cuts[0]
        rotated = dna[first:] + dna[:first]
        adjusted = [c - first for c in cuts[1:]] + [len(dna)]
        frags, prev = [], 0
        for c in adjusted:
            frags.append(rotated[prev:c])
            prev = c
        return frags
    boundaries = [0] + cuts + [len(dna)]
    return [dna[boundaries[i]:boundaries[i + 1]] for i in range(len(boundaries) - 1)]


def double_digest(dna: str, enzyme1: str, enzyme2: str, circular: bool = False) -> list:
    """Digest with two enzymes simultaneously; merges cut positions before slicing."""
    all_cuts = sorted(set(find_cut_positions(dna, enzyme1) + find_cut_positions(dna, enzyme2)))
    if not all_cuts:
        return [dna]
    if circular:
        return digest(dna, enzyme1, circular=True) if enzyme1 == enzyme2 else _slice_circular(dna, all_cuts)
    boundaries = [0] + all_cuts + [len(dna)]
    return [dna[boundaries[i]:boundaries[i + 1]] for i in range(len(boundaries) - 1)]


def _slice_circular(dna: str, cuts: list) -> list:
    first = cuts[0]
    rotated = dna[first:] + dna[:first]
    adjusted = [c - first for c in cuts[1:]] + [len(dna)]
    frags, prev = [], 0
    for c in adjusted:
        frags.append(rotated[prev:c])
        prev = c
    return frags


def compute_overhang(enzyme_name: str) -> tuple:
    """Return (overhang_seq, overhang_type) with type in {'5prime', '3prime', 'blunt'}."""
    site, cut_top, cut_bot = RESTRICTION_ENZYMES[enzyme_name]
    if cut_top == cut_bot:
        return ('', 'blunt')
    if cut_top < cut_bot:
        return (site[cut_top:cut_bot], '5prime')
    return (reverse_complement(site)[cut_bot:cut_top], '3prime')


def are_compatible(enzyme1: str, enzyme2: str) -> bool:
    """True if the two enzymes' ends can be ligated (same overhang type + sequence, or both blunt)."""
    oh1, type1 = compute_overhang(enzyme1)
    oh2, type2 = compute_overhang(enzyme2)
    if type1 != type2:
        return False
    return True if type1 == 'blunt' else oh1 == oh2


if __name__ == '__main__':
    plasmid = 'AAGCTTGCATGCCTGCAGGTCGACGGATCCCCGGAATTCGAGCTCGGTACCCGGGGATCCTCTAGAGTCGAC'
    frags = digest(plasmid, 'EcoRI', circular=True)
    assert sum(len(f) for f in frags) == len(plasmid)
    assert are_compatible('SalI', 'XhoI') is True   # both leave 5'-TCGA
    assert are_compatible('EcoRI', 'BamHI') is False
    print('digestion + compatibility checks passed')

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(seq: str) -> float:
    s = seq.upper()
    return (s.count('G') + s.count('C')) / len(s)


def tm_basic(primer: str) -> float:
    """4+2 rule: Tm = 4*(G+C) + 2*(A+T). Only reliable for <20 bp primers."""
    s = primer.upper()
    return 4 * (s.count('G') + s.count('C')) + 2 * (s.count('A') + s.count('T'))


def tm_salt_adjusted(primer: str, salt_mm: float = 50.0) -> float:
    """Wallace rule variant: Tm = 81.5 + 16.6*log10([Na+]) + 41*GC - 675/N."""
    n = len(primer)
    gc = gc_content(primer)
    return 81.5 + 16.6 * math.log10(salt_mm / 1000) + 41 * gc - 675 / n


def tm_nearest_neighbor(primer: str, dna_conc_nm: float = 250.0, salt_mm: float = 50.0) -> float:
    """SantaLucia 1998 nearest-neighbor Tm; most accurate for 18-30 bp primers."""
    seq, R = primer.upper(), 1.987  # cal/mol/K
    dH = dS = 0.0
    for end_base in (seq[0], seq[-1]):
        h, s = NN_INIT_GC if end_base in 'GC' else NN_INIT_AT
        dH += h
        dS += s
    for i in range(len(seq) - 1):
        h, s = NN_PARAMS.get(seq[i:i + 2], (-8.0, -20.0))
        dH += h
        dS += s
    dS += 0.368 * (len(seq) - 1) * math.log(salt_mm / 1000)  # salt correction
    ct = dna_conc_nm * 1e-9
    return (dH * 1000) / (dS + R * math.log(ct / 4)) - 273.15


def reverse_complement(seq: str) -> str:
    comp = {'A': 'T', 'T': 'A', 'G': 'C', 'C': 'G', 'N': 'N'}
    return ''.join(comp.get(b, 'N') for b in reversed(seq.upper()))


def design_primer(template: str, start: int, direction: str, target_tm: float = 60.0,
                   min_len: int = 18, max_len: int = 28, tm_fn=tm_nearest_neighbor) -> dict:
    """Grow a primer from `start` until it reaches `target_tm` (or hits max_len).
    direction='fwd' reads left-to-right; 'rev' reads right-to-left and returns the reverse complement.
    """
    best = None
    for length in range(min_len, max_len + 1):
        if direction == 'fwd':
            seq = template[start:start + length]
        else:
            seq = reverse_complement(template[start - length + 1:start + 1])
        if len(seq) < length:
            break
        tm = tm_fn(seq)
        candidate = {'seq': seq, 'length': length, 'tm': tm, 'gc': gc_content(seq)}
        if best is None or abs(tm - target_tm) < abs(best['tm'] - target_tm):
            best = candidate
        if tm >= target_tm:
            break
    return best


def check_3prime_run(primer: str, run_len: int = 4) -> bool:
    """True if the 3' end is a homopolymer run of `run_len`+ identical bases (risks mispriming)."""
    return len(set(primer[-run_len:].upper())) == 1


def count_3prime_complementarity(primer1: str, primer2: str, check_len: int = 5) -> int:
    """Score (0-check_len) of 3'-end complementarity between two primers; >=3 flags dimer risk."""
    tail1 = primer1[-check_len:].upper()
    tail2 = reverse_complement(primer2[-check_len:]).upper()
    return sum(a == b for a, b in zip(tail1, tail2))


def add_re_site_to_primer(primer: str, site: str, clamp: str = 'GCGC') -> str:
    """Prepend a clamp + restriction recognition site to a primer for cloning (e.g. site='GGATCC')."""
    return clamp + site + primer


if __name__ == '__main__':
    mcs = 'GAATTCGAGCTCGGTACCCGGGGATCCTCTAGAGTCGACCTGCAGGCATGCAAGCTTGGCGTAATCATGGTCATAGCTG'
    fwd = design_primer(mcs, start=0, direction='fwd', target_tm=60)
    rev = design_primer(mcs, start=len(mcs) - 1, direction='rev', target_tm=60)
    assert 18 <= fwd['length'] <= 28 and abs(fwd['tm'] - 60) < 15
    assert not check_3prime_run(fwd['seq'], run_len=6)
    cloning_fwd = add_re_site_to_primer(fwd['seq'], 'GGATCC')  # BamHI
    assert cloning_fwd.startswith('GCGCGGATCC')
    print(f'fwd Tm={fwd["tm"]:.1f} rev Tm={rev["tm"]:.1f} diff={abs(fwd["tm"]-rev["tm"]):.1f}')

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]
        for size in sizes:
            if size < 50:
                continue
            y = math.log10(size)
            color = '#444444' if label == 'Ladder' else '#e05c1a'
            ax.plot([x - 0.35, x + 0.35], [y, y], color=color, linewidth=4, solid_capstyle='butt')
        ax.text(x, 1.85, label, ha='center', va='top', fontsize=9, fontweight='bold')
    ax.set_title('Simulated Agarose Gel')
    plt.tight_layout()
    return 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