| name | molecular-cloning |
| description | Molecular cloning simulation and design. PCR amplicon prediction, restriction enzyme digestion, Golden Gate and Gibson assembly simulation, primer design, CRISPR sgRNA design, and plasmid annotation. For protein-level sequence analysis use biopython or esm; for database lookups use gene-database or ensembl-database. |
| category | biology |
| license | MIT license |
| metadata | {"skill-author":"InkVell Inc."} |
Molecular Cloning: Sequence Engineering & Cloning Design
Overview
Molecular Cloning provides computational tools for simulating and designing molecular cloning workflows. This skill covers PCR amplicon prediction with primer binding analysis, restriction enzyme digestion simulation, Golden Gate and Gibson assembly design and verification, primer design with thermodynamic calculations, CRISPR sgRNA design with off-target scoring, and plasmid feature annotation. All simulations use Biopython's Bio.Restriction and Bio.SeqUtils for accurate enzyme and sequence handling.
When to Use This Skill
- Predicting PCR amplicons from primer sequences and templates
- Simulating restriction enzyme digestions and predicting fragment sizes
- Designing Golden Gate assembly with 4bp overhang compatibility
- Planning Gibson assembly with overlap design
- Designing PCR primers with Tm and specificity constraints
- Designing CRISPR sgRNAs and scoring off-target potential
- Annotating plasmid features (promoters, CDS, terminators, origins)
- Generating plasmid maps in GenBank format
Related Skills: For protein-level sequence analysis use biopython or esm. For gene/transcript lookups use gene-database or ensembl-database. For synthetic biology circuit design use synthetic-biology.
Installation
uv pip install biopython primer3-py numpy
Quick Start
from Bio.Seq import Seq
from Bio.Restriction import BamHI, EcoRI
from Bio.SeqUtils import MeltingTemp as mt
sequence = Seq("ATCGATCGGGATCCATCGATCGAATTCATCGATCG")
print(f"BamHI cuts at: {BamHI.search(sequence)}")
print(f"EcoRI cuts at: {EcoRI.search(sequence)}")
primer = Seq("ATCGATCGGATCCATCGATCG")
tm = mt.Tm_NN(primer)
print(f"Primer Tm: {tm:.1f} C")
Core Capabilities
1. PCR Simulation
Predict amplicon from primer binding on template.
from Bio.Seq import Seq
from Bio.SeqUtils import MeltingTemp as mt
import re
def find_primer_binding(template, primer, max_mismatches=2):
"""Find primer binding sites on template (both strands).
Returns list of (position, strand, mismatches) tuples.
"""
template_str = str(template).upper()
primer_str = str(primer).upper()
rc_template = str(template.reverse_complement()).upper()
sites = []
for i in range(len(template_str) - len(primer_str) + 1):
region = template_str[i:i+len(primer_str)]
mismatches = sum(a != b for a, b in zip(primer_str, region))
if mismatches <= max_mismatches:
sites.append((i, '+', mismatches))
for i in range(len(rc_template) - len(primer_str) + 1):
region = rc_template[i:i+len(primer_str)]
mismatches = sum(a != b for a, b in zip(primer_str, region))
if mismatches <= max_mismatches:
pos = len(template_str) - i - len(primer_str)
sites.append((pos, '-', mismatches))
return sites
def simulate_pcr(template, fwd_primer, rev_primer, max_mismatches=2):
"""Simulate PCR and predict amplicon.
Args:
template: Bio.Seq template sequence (can be circular)
fwd_primer: forward primer sequence
rev_primer: reverse primer sequence (as ordered, 5'->3')
"""
fwd_sites = find_primer_binding(template, fwd_primer, max_mismatches)
rev_rc = Seq(str(rev_primer)).reverse_complement()
rev_sites = find_primer_binding(template, rev_rc, max_mismatches)
amplicons = []
for f_pos, f_strand, f_mm in fwd_sites:
if f_strand != '+':
continue
for r_pos, r_strand, r_mm in rev_sites:
if r_strand != '-':
continue
if r_pos > f_pos:
amp_len = r_pos + len(str(rev_primer)) - f_pos
if 50 < amp_len < 10000:
amplicon = template[f_pos:r_pos + len(str(rev_primer))]
amplicons.append({
'start': f_pos,
'end': r_pos + len(str(rev_primer)),
'length': amp_len,
'fwd_mismatches': f_mm,
'rev_mismatches': r_mm,
'sequence': str(amplicon)
})
fwd_tm = mt.Tm_NN(fwd_primer)
rev_tm = mt.Tm_NN(rev_primer)
print(f"Forward primer Tm: {fwd_tm:.1f} C")
print(f"Reverse primer Tm: {rev_tm:.1f} C")
print(f"Tm difference: {abs(fwd_tm - rev_tm):.1f} C")
print(f"Predicted amplicons: {len(amplicons)}")
for i, amp in enumerate(amplicons):
print(f" Amplicon {i+1}: {amp['length']} bp "
f"(pos {amp['start']}-{amp['end']}, "
f"mismatches: fwd={amp['fwd_mismatches']}, rev={amp['rev_mismatches']})")
return amplicons
2. Restriction Digestion
Simulate enzyme digestion and predict fragments.
from Bio.Seq import Seq
from Bio.Restriction import *
from Bio.Restriction import RestrictionBatch, Analysis
def restriction_digest(sequence, enzymes, is_linear=True):
"""Simulate restriction enzyme digestion.
Args:
sequence: Bio.Seq DNA sequence
enzymes: list of enzyme names (e.g., ['EcoRI', 'BamHI'])
is_linear: True for linear DNA, False for circular
"""
rb = RestrictionBatch()
for enz_name in enzymes:
rb.add(eval(enz_name))
analysis = Analysis(rb, sequence, linear=is_linear)
results = analysis.full()
all_cut_sites = []
for enzyme, sites in results.items():
if sites:
print(f"{enzyme}: cuts at positions {sites}")
all_cut_sites.extend(sites)
else:
print(f"{enzyme}: no cut sites")
if not all_cut_sites:
print(f"No cuts — single fragment: {len(sequence)} bp")
return [len(sequence)]
all_cut_sites = sorted(set(all_cut_sites))
if is_linear:
positions = [0] + all_cut_sites + [len(sequence)]
else:
positions = all_cut_sites
fragments = []
if is_linear:
for i in range(len(positions) - 1):
fragments.append(positions[i+1] - positions[i])
else:
for i in range(len(positions)):
next_i = (i + 1) % len(positions)
if next_i == 0:
frag = len(sequence) - positions[i] + positions[0]
else:
frag = positions[next_i] - positions[i]
fragments.append(frag)
fragments.sort(reverse=True)
print(f"\nFragments ({len(fragments)}): {fragments}")
print(f"Total: {sum(fragments)} bp")
return fragments
seq = Seq("ATCG" * 100 + "GGATCC" + "ATCG" * 50 + "GAATTC" + "ATCG" * 75)
frags = restriction_digest(seq, ['BamHI', 'EcoRI'], is_linear=True)
3. Golden Gate Assembly
Design and simulate Golden Gate cloning.
from Bio.Seq import Seq
def design_golden_gate(parts, enzyme='BsaI'):
"""Design Golden Gate assembly with verified overhang compatibility.
Args:
parts: list of dicts with 'name', 'sequence' keys
enzyme: Type IIS enzyme ('BsaI' or 'BpiI')
"""
standard_overhangs = [
'AATG',
'AGGT',
'TTCG',
'GCTT',
'CGCT',
]
if len(parts) + 1 > len(standard_overhangs):
raise ValueError(f"Too many parts ({len(parts)}) for standard overhang set")
overhangs_used = standard_overhangs[:len(parts) + 1]
for i, oh in enumerate(overhangs_used):
rc = str(Seq(oh).reverse_complement())
if oh == rc:
print(f"WARNING: Overhang {oh} is palindromic — may self-ligate")
for j, oh2 in enumerate(overhangs_used):
if i != j and oh == oh2:
raise ValueError(f"Duplicate overhangs: position {i} and {j}")
assembly = []
for i, part in enumerate(parts):
left_oh = overhangs_used[i]
right_oh = overhangs_used[i + 1]
if enzyme == 'BsaI':
recognition = 'GGTCTC'
spacer = 'N'
else:
recognition = 'GAAGAC'
spacer = 'NN'
assembled_part = f"{recognition}{spacer}{left_oh}{part['sequence']}{right_oh}"
assembly.append({
'name': part['name'],
'left_overhang': left_oh,
'right_overhang': right_oh,
'part_length': len(part['sequence']),
'total_length': len(assembled_part)
})
print(f"Part {i+1} ({part['name']}): [{left_oh}]--{len(part['sequence'])}bp--[{right_oh}]")
total_insert = sum(len(p['sequence']) for p in parts) + \
len(overhangs_used) * 4
print(f"\nAssembly: {len(parts)} parts")
print(f"Total insert: ~{total_insert} bp (excluding vector)")
print(f"Overhang set: {' -> '.join(overhangs_used)}")
return assembly
parts = [
{'name': 'Promoter', 'sequence': 'TTGACAATTAATCATCGGCTCG' * 5},
{'name': 'RBS', 'sequence': 'AAGGAGATATACAT'},
{'name': 'GFP_CDS', 'sequence': 'ATGGTGAGCAAGGGCGAG' * 40},
{'name': 'Terminator', 'sequence': 'TTTTTTTTTTT' * 4},
]
assembly = design_golden_gate(parts)
4. Gibson Assembly
Design overlapping fragments for Gibson assembly.
from Bio.Seq import Seq
from Bio.SeqUtils import MeltingTemp as mt
def design_gibson_assembly(fragments, overlap_length=30):
"""Design Gibson assembly with overlap primers.
Args:
fragments: list of dicts with 'name', 'sequence' keys (in assembly order)
overlap_length: overlap length in bp (20-40 recommended)
"""
assembly_plan = []
for i in range(len(fragments)):
current = fragments[i]
next_frag = fragments[(i + 1) % len(fragments)]
if i == 0:
fwd_primer = current['sequence'][:20]
else:
prev = fragments[i - 1]
overlap_5 = prev['sequence'][-overlap_length:]
fwd_primer = overlap_5 + current['sequence'][:20]
overlap_3 = next_frag['sequence'][:overlap_length]
rev_binding = str(Seq(current['sequence'][-20:]).reverse_complement())
rev_primer = str(Seq(overlap_3).reverse_complement()) + rev_binding
fwd_tm = mt.Tm_NN(Seq(current['sequence'][:20]))
rev_tm = mt.Tm_NN(Seq(current['sequence'][-20:]))
assembly_plan.append({
'fragment': current['name'],
'fragment_length': len(current['sequence']),
'fwd_primer': fwd_primer,
'rev_primer': rev_primer,
'fwd_tm': fwd_tm,
'rev_tm': rev_tm,
'overlap_length': overlap_length
})
print(f"Fragment {i+1} ({current['name']}): {len(current['sequence'])} bp")
print(f" Fwd primer: {fwd_primer[:40]}... ({len(fwd_primer)} nt, Tm={fwd_tm:.1f} C)")
print(f" Rev primer: {rev_primer[:40]}... ({len(rev_primer)} nt, Tm={rev_tm:.1f} C)")
total = sum(len(f['sequence']) for f in fragments)
print(f"\nTotal assembled length: ~{total} bp")
return assembly_plan
5. Primer Design
Design primers with thermodynamic constraints.
import primer3
def design_primers(template_seq, target_start, target_length,
product_size_range=(200, 500), tm_target=60):
"""Design PCR primers using primer3.
Args:
template_seq: template DNA sequence (string)
target_start: start position of target region
target_length: length of target region
product_size_range: (min, max) product size
tm_target: target melting temperature
"""
result = primer3.bindings.design_primers(
seq_args={
'SEQUENCE_TEMPLATE': template_seq,
'SEQUENCE_TARGET': [target_start, target_length],
},
global_args={
'PRIMER_NUM_RETURN': 5,
'PRIMER_OPT_SIZE': 20,
'PRIMER_MIN_SIZE': 18,
'PRIMER_MAX_SIZE': 25,
'PRIMER_OPT_TM': tm_target,
'PRIMER_MIN_TM': tm_target - 5,
'PRIMER_MAX_TM': tm_target + 5,
'PRIMER_MIN_GC': 40,
'PRIMER_MAX_GC': 60,
'PRIMER_PRODUCT_SIZE_RANGE': [list(product_size_range)],
'PRIMER_MAX_SELF_COMPLEMENT': 6,
'PRIMER_MAX_SELF_END': 3,
'PRIMER_MAX_HAIRPIN_TH': 47,
}
)
primers = []
for i in range(result.get('PRIMER_PAIR_NUM_RETURNED', 0)):
pair = {
'left_seq': result[f'PRIMER_LEFT_{i}_SEQUENCE'],
'right_seq': result[f'PRIMER_RIGHT_{i}_SEQUENCE'],
'left_tm': result[f'PRIMER_LEFT_{i}_TM'],
'right_tm': result[f'PRIMER_RIGHT_{i}_TM'],
'left_gc': result[f'PRIMER_LEFT_{i}_GC_PERCENT'],
'right_gc': result[f'PRIMER_RIGHT_{i}_GC_PERCENT'],
'product_size': result[f'PRIMER_PAIR_{i}_PRODUCT_SIZE'],
'penalty': result[f'PRIMER_PAIR_{i}_PENALTY'],
}
primers.append(pair)
print(f"Pair {i+1}: product={pair['product_size']}bp, penalty={pair['penalty']:.2f}")
print(f" Fwd: {pair['left_seq']} (Tm={pair['left_tm']:.1f}, GC={pair['left_gc']:.0f}%)")
print(f" Rev: {pair['right_seq']} (Tm={pair['right_tm']:.1f}, GC={pair['right_gc']:.0f}%)")
return primers
6. CRISPR sgRNA Design
Design guide RNAs for CRISPR-Cas9 knockout.
from Bio.Seq import Seq
import re
def design_crispr_guides(target_seq, pam='NGG', guide_length=20, top_n=10):
"""Design CRISPR sgRNAs for SpCas9.
Args:
target_seq: target gene/region sequence (string)
pam: PAM sequence ('NGG' for SpCas9)
guide_length: guide RNA length (typically 20)
top_n: number of top guides to return
"""
seq = str(target_seq).upper()
rc_seq = str(Seq(seq).reverse_complement())
guides = []
for i in range(len(seq) - guide_length - len(pam) + 1):
pam_site = seq[i + guide_length:i + guide_length + 3]
if re.match(pam.replace('N', '.'), pam_site):
guide_seq = seq[i:i + guide_length]
score = score_guide(guide_seq)
guides.append({
'sequence': guide_seq,
'pam': pam_site,
'position': i,
'strand': '+',
'score': score,
'full_target': guide_seq + pam_site
})
for i in range(len(rc_seq) - guide_length - len(pam) + 1):
pam_site = rc_seq[i + guide_length:i + guide_length + 3]
if re.match(pam.replace('N', '.'), pam_site):
guide_seq = rc_seq[i:i + guide_length]
score = score_guide(guide_seq)
guides.append({
'sequence': guide_seq,
'pam': pam_site,
'position': len(seq) - i - guide_length,
'strand': '-',
'score': score,
'full_target': guide_seq + pam_site
})
guides.sort(key=lambda x: x['score'], reverse=True)
print(f"Found {len(guides)} potential guide RNAs")
print(f"\nTop {min(top_n, len(guides))} guides:")
for i, g in enumerate(guides[:top_n]):
print(f" {i+1}. {g['sequence']} {g['pam']} "
f"(strand={g['strand']}, pos={g['position']}, score={g['score']:.2f})")
return guides[:top_n]
def score_guide(guide_seq):
"""Heuristic guide scoring based on known design rules."""
score = 50.0
gc = (guide_seq.count('G') + guide_seq.count('C')) / len(guide_seq)
if 0.4 <= gc <= 0.7:
score += 10
elif gc < 0.3 or gc > 0.8:
score -= 20
if 'TTTT' in guide_seq:
score -= 30
if guide_seq[-1] == 'G':
score += 5
if guide_seq[-2:] == 'GG':
score -= 5
rc = str(Seq(guide_seq).reverse_complement())
matches = sum(a == b for a, b in zip(guide_seq, rc))
if matches > 12:
score -= 15
return max(score, 0)
7. Plasmid Annotation
Identify and annotate features in plasmid sequences.
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.SeqFeature import SeqFeature, FeatureLocation
from Bio import SeqIO
from Bio.Restriction import RestrictionBatch, Analysis, CommOnly
def annotate_plasmid(sequence, name='plasmid'):
"""Annotate plasmid features and restriction sites.
Args:
sequence: plasmid DNA sequence (string)
name: plasmid name
"""
seq = Seq(sequence)
record = SeqRecord(seq, id=name, name=name,
description=f'{name} annotated plasmid',
annotations={'molecule_type': 'DNA', 'topology': 'circular'})
feature_patterns = {
'T7_promoter': 'TAATACGACTCACTATAG',
'lac_operator': 'AATTGTGAGCGGATAACAATT',
'RBS_consensus': 'AAGGAG',
'T7_terminator': 'CTAGCATAACCCCTTGGGGCCTCTAAACGGGTCTTGAGG',
'ColE1_origin': 'CCTGTTTTGGCGGATGAGAGAAG',
}
for feat_name, pattern in feature_patterns.items():
pos = str(seq).find(pattern)
if pos >= 0:
record.features.append(SeqFeature(
FeatureLocation(pos, pos + len(pattern)),
type='misc_feature',
qualifiers={'label': feat_name}
))
print(f"Found {feat_name} at position {pos}")
for strand, nuc in [(1, seq), (-1, seq.reverse_complement())]:
for frame in range(3):
trans = str(nuc[frame:].translate())
start = 0
while True:
start = trans.find('M', start)
if start == -1:
break
stop = trans.find('*', start)
if stop == -1:
stop = len(trans)
orf_len = (stop - start) * 3
if orf_len >= 300:
if strand == 1:
dna_start = frame + start * 3
dna_end = frame + stop * 3 + 3
else:
dna_end = len(seq) - frame - start * 3
dna_start = len(seq) - frame - stop * 3 - 3
record.features.append(SeqFeature(
FeatureLocation(min(dna_start, dna_end),
max(dna_start, dna_end), strand),
type='CDS',
qualifiers={'label': f'ORF_{orf_len}bp'}
))
print(f"ORF: {orf_len}bp at {dna_start}-{dna_end} (strand {'+' if strand==1 else '-'})")
start = stop + 1
rb = CommOnly
analysis = Analysis(rb, seq, linear=False)
unique_sites = {str(enz): sites for enz, sites in analysis.full().items()
if len(sites) == 1}
print(f"\nUnique restriction sites: {len(unique_sites)}")
for enz, sites in sorted(unique_sites.items()):
print(f" {enz}: {sites[0]}")
return record
Typical Workflows
Workflow 1: Design PCR Primers and Predict Amplicon
from Bio.Seq import Seq
template = Seq("ATCGATCGATCGATCGATCGATCGATCG" * 50)
fwd = Seq("ATCGATCGATCGATCGATCG")
rev = Seq("CGATCGATCGATCGATCGAT")
amplicons = simulate_pcr(template, fwd, rev)
if amplicons:
print(f"Amplicon size: {amplicons[0]['length']} bp")
Workflow 2: Plan Golden Gate Assembly with 4 Parts
parts = [
{'name': 'J23100_promoter', 'sequence': 'TTGACGGCTAGCTCAGTCCTAGGTACAGTGCTAGC'},
{'name': 'B0034_RBS', 'sequence': 'AAAGAGGAGAAA'},
{'name': 'GFP', 'sequence': 'ATGGTGAGCAAGGGCGAGGAG' + 'NNN' * 230 + 'TAA'},
{'name': 'B0015_terminator', 'sequence': 'CCAGGCATCAAATAAAACGAAAGGCTCAGTCGAAAG'},
]
assembly = design_golden_gate(parts)
Workflow 3: Design CRISPR Knockout sgRNAs
target_exon = "ATGCGATCGATCGATCGATCGATCGAGGCGATCGATCGATCGATCGATCGATCGATCGATCG"
guides = design_crispr_guides(target_exon, pam='NGG', top_n=5)
Best Practices
- Primer design — aim for Tm within 2 C between forward and reverse primers; 40-60% GC content; avoid 3' complementarity
- Golden Gate — verify 4bp overhangs are unique and non-palindromic; use standardized overhang sets from published studies
- Gibson assembly — use 20-40bp overlaps; ensure fragments have similar Tm at overlap junctions; minimum 2 fragments
- CRISPR guides — avoid poly-T (>4 T) which terminates U6/H1 transcription; prefer guides in early exons; validate against off-target databases
- Restriction digestion — always check both orientations; verify compatible cohesive ends for ligation; methylation sensitivity can block cutting
- Plasmid annotation — verify ORFs are in correct reading frame; check for cryptic promoters; map all unique restriction sites for future cloning
Troubleshooting
Problem: PCR simulation finds no amplicons
Solution: Check primer orientation (forward should match sense strand 5'->3'). Increase max_mismatches. Verify template contains the target region.
Problem: Restriction enzyme doesn't cut expected site
Solution: Check for CpG methylation sensitivity (dam/dcm methylation). Verify site isn't in a modified context. Use isoschizomers() to find alternatives.
Problem: Golden Gate assembly produces wrong product
Solution: Verify all overhangs are unique. Check enzyme is BsaI (not BsmBI) — different cut distances. Ensure parts are in correct orientation.
Problem: primer3 returns no primers
Solution: Relax constraints: widen Tm range, increase max size, lower GC requirement. Ensure target region has sufficient flanking sequence.
Problem: CRISPR guide design returns few candidates
Solution: Expand search region (use full gene, not just one exon). Try alternative PAM sequences (NAG for relaxed SpCas9). Consider Cas12a (TTTV PAM) for AT-rich regions.
Resources
