← All skills

dna-chisel

utility

Turn a designed protein into orderable DNA with DNA Chisel — codon-optimize for a host, hold GC inside a window, strip Golden Gate and other restriction sites from both strands, avoid repeats and hairpins, export an annotated GenBank of the edits.

Designing a coding sequence with DNA Chisel

A protein design is not an order. Between the two sits a sequence that has to satisfy a dozen unrelated requirements at once — the host's codon preferences, the synthesis vendor's GC and repeat limits, the assembly method's forbidden restriction sites, and whatever else the downstream vector demands. Hand-editing codons until all of them hold is where errors enter, because fixing a GC window reintroduces a site you removed two edits ago.

DNA Chisel inverts that. You declare what must be true and what you would prefer, and a constraint solver finds a sequence inside the intersection. The declaration is the useful artefact: it is reviewable, it is diffable, and it is what you re-run when the vendor changes a rule.

Nothing here needs an account, an API key, or a GPU. DNA Chisel is MIT-licensed, installs from the package index, and runs on CPU in seconds for a gene-sized sequence. The one optional network call is noted where it appears.

Constraints and objectives are different things

This distinction decides whether the run succeeds, and it is the first thing to get right.

  • A constraint is a hard requirement. The solver must satisfy every one of them or it raises. Translation must be preserved; the BsaI site must be gone; GC must stay between 0.35 and 0.65 in every 50 bp window.
  • An objective is a preference, scored and maximized within whatever room the constraints leave. Codon optimization belongs here, almost always. Make it a constraint and you are asserting that a particular codon in a particular position is non-negotiable, which it is not.

The two run in order, and both calls are needed:

problem.resolve_constraints()   # find any sequence satisfying the hard requirements
problem.optimize()              # improve the objectives without breaking them

Calling optimize() alone on a sequence that violates a constraint optimizes a sequence that was never valid. Calling resolve_constraints() alone gives a legal sequence with no codon optimization at all.

Install

pip install dnachisel biopython

python_codon_tables and numpy come with it. Everything below was executed against dnachisel 3.2.16.

From protein to DNA

reverse_translate gives you a starting sequence. By default it picks one fixed codon per amino acid, which produces a low-complexity sequence with long repeats — fine as a solver starting point, bad as a final answer. Pass randomize_codons=True for a more realistic start, then let the constraints do the work.

from Bio.Seq import Seq
import dnachisel as dc

protein = "MKVLATGIVALQDLEYRFKAWCNHPSTM"

naive = dc.reverse_translate(protein)
varied = dc.reverse_translate(protein, randomize_codons=True)

for label, seq in [("one-codon-per-aa", naive), ("randomized", varied)]:
    print(label, len(seq), "bp  GC %.3f" % dc.biotools.gc_content(seq))
    print("   ", seq)
    print("    round-trips:", str(Seq(seq).translate()) == protein)

Neither output is a design. It is the input to one.

Audit the sequence before you optimize it

Run the checks before you change anything, so you know what the solver is actually fixing. The sequence here is the native Aequorea victoria GFP coding sequence, pulled from its ENA record.

import re, urllib.request
from Bio import SeqIO
import dnachisel as dc

urllib.request.urlretrieve(
    "https://www.ebi.ac.uk/ena/browser/api/embl/M62653.1", "gfp.embl")
record = SeqIO.read("gfp.embl", "embl")
cds = str(next(f for f in record.features if f.type == "CDS").extract(record.seq)).upper()

print("CDS", len(cds), "bp  overall GC %.3f" % dc.biotools.gc_content(cds))

# Restriction sites, counted on BOTH strands -- this is the part done wrong by hand.
for enzyme in ["BsaI", "BsmBI", "EcoRI"]:
    site = str(dc.EnzymeSitePattern(enzyme)).split("(")[1].rstrip(")")
    fwd = [m.start() for m in re.finditer(site, cds)]
    rev = [m.start() for m in re.finditer(dc.reverse_complement(site), cds)]
    print(f"{enzyme} ({site})  plus strand {fwd}  minus strand {rev}")

# Local GC, which is what a synthesis vendor actually rejects on.
windows = [dc.biotools.gc_content(cds[i:i + 50]) for i in range(0, len(cds) - 50, 10)]
print("50 bp GC window range %.3f - %.3f" % (min(windows), max(windows)))

Native avGFP is the right teaching case because its problems are invisible to a forward reading. It carries no GGTCTC anywhere — and it is still unusable for Golden Gate with BsaI, because GAGACC sits at position 643. Restriction sites are palindrome-free and strand-agnostic; a grep for the site as written finds half of them. AvoidPattern scans both strands by default, which is the whole reason to use it.

The design run

import urllib.request
from Bio import SeqIO
from Bio.Seq import Seq
import numpy
import dnachisel as dc

urllib.request.urlretrieve(
    "https://www.ebi.ac.uk/ena/browser/api/embl/M62653.1", "gfp.embl")
record = SeqIO.read("gfp.embl", "embl")
cds = str(next(f for f in record.features if f.type == "CDS").extract(record.seq)).upper()
numpy.random.seed(0)

problem = dc.DnaOptimizationProblem(
    sequence=cds,
    constraints=[
        dc.EnforceTranslation(location=(0, len(cds), 1)),   # the protein is fixed
        dc.EnforceGCContent(mini=0.35, maxi=0.65, window=50),
        dc.AvoidPattern("BsaI_site"),                        # Golden Gate
        dc.AvoidPattern("BsmBI_site"),
        dc.AvoidPattern("9xA"),                              # homopolymer run
        dc.UniquifyAllKmers(12),                             # no repeated 12-mers
    ],
    objectives=[
        dc.CodonOptimize(species="e_coli", method="use_best_codon",
                         location=(0, len(cds), 1)),
    ],
    logger=None,
)

print(problem.constraints_text_summary())      # what is wrong with the input
problem.resolve_constraints()
problem.optimize()
print(problem.constraints_text_summary())      # what the solver fixed
print(problem.objectives_text_summary())

designed = problem.sequence
print("edits:", sum(1 for a, b in zip(cds, designed) if a != b), "of", len(cds))
print("final GC %.3f" % dc.biotools.gc_content(designed))
print("protein unchanged:", str(Seq(designed).translate()) == str(Seq(cds).translate()))

Two details in that block carry weight.

location=(0, len(cds), 1) is not optional for EnforceTranslation. The third element is the strand, and the span tells the solver where the reading frame starts. Give it the wrong span and it will faithfully preserve the translation of the wrong frame.

CodonOptimize always pairs with EnforceTranslation over the same location. On its own it is free to change amino acids, because nothing told it not to.

Choosing a codon-optimization method

Three are available, and they answer different questions. All three were run against the same avGFP CDS with the same constraints and numpy.random.seed(0).

import urllib.request
from Bio import SeqIO
from Bio.Seq import Seq
import numpy
import dnachisel as dc

urllib.request.urlretrieve(
    "https://www.ebi.ac.uk/ena/browser/api/embl/M62653.1", "gfp.embl")
record = SeqIO.read("gfp.embl", "embl")
cds = str(next(f for f in record.features if f.type == "CDS").extract(record.seq)).upper()

def run(**codon_kwargs):
    numpy.random.seed(0)
    problem = dc.DnaOptimizationProblem(
        sequence=cds,
        constraints=[dc.EnforceTranslation(location=(0, len(cds), 1)),
                     dc.EnforceGCContent(mini=0.35, maxi=0.65, window=50),
                     dc.AvoidPattern("BsaI_site"), dc.AvoidPattern("BsmBI_site")],
        objectives=[dc.CodonOptimize(location=(0, len(cds), 1), **codon_kwargs)],
        logger=None)
    problem.resolve_constraints()
    problem.optimize()
    assert str(Seq(problem.sequence).translate()) == str(Seq(cds).translate())
    return problem

for kwargs in [dict(species="e_coli", method="use_best_codon"),
               dict(species="e_coli", method="match_codon_usage"),
               dict(species="e_coli", original_species="h_sapiens",
                    method="harmonize_rca")]:
    p = run(**kwargs)
    print("%-18s edits=%3d  GC=%.3f" % (
        kwargs["method"],
        sum(1 for a, b in zip(cds, p.sequence) if a != b),
        dc.biotools.gc_content(p.sequence)))
  • use_best_codon replaces every codon with the host's most frequent synonym. Watch the objectives summary and you will see it resolve to an objective reported as MaximizeCAI, which is what it is. It moves the most bases of the three and produces the most repetitive sequence, so it is the method that fights the repeat and GC constraints hardest — 166 edits here against 104 for the distribution-matching run.
  • match_codon_usage reproduces the host's codon distribution rather than its argmax, so rare codons survive at roughly their natural frequency. Preferred when translational pausing matters — folding-sensitive proteins, membrane proteins.
  • harmonize_rca maps each codon's usage in its original host onto a codon of comparable usage in the target, preserving the native rhythm. It requires both species and original_species; omitting the second raises ValueError: Provide either an species name or a codon usage table, which names the missing table rather than the missing argument.

Nine codon tables ship with the package and need no network — e_coli, h_sapiens, m_musculus, s_cerevisiae, b_subtilis, d_melanogaster, c_elegans, g_gallus and m_musculus_domesticus. Passing a numeric TaxID instead downloads the table at run time, so a pipeline keyed on a TaxID fails on an offline host while the same pipeline keyed on a name does not.

The solver is stochastic — seed it

Repeated runs of the same problem give different sequences. Across three unseeded runs of the design above the edit count moved between 160 and 164 and the final GC between 0.494 and 0.497, all of them valid. Nothing is wrong; the search starts from a random point.

For a design you intend to order, synthesise, or put in a paper, seed numpy's global generator before constructing the problem. Seeding it alone is sufficient — the random module is not involved.

import urllib.request
from Bio import SeqIO
import numpy
import dnachisel as dc

urllib.request.urlretrieve(
    "https://www.ebi.ac.uk/ena/browser/api/embl/M62653.1", "gfp.embl")
record = SeqIO.read("gfp.embl", "embl")
cds = str(next(f for f in record.features if f.type == "CDS").extract(record.seq)).upper()

def design(seed):
    numpy.random.seed(seed)
    problem = dc.DnaOptimizationProblem(
        sequence=cds,
        constraints=[dc.EnforceTranslation(location=(0, len(cds), 1)),
                     dc.EnforceGCContent(mini=0.35, maxi=0.65, window=50),
                     dc.AvoidPattern("BsaI_site"), dc.AvoidPattern("BsmBI_site"),
                     dc.AvoidPattern("9xA"), dc.UniquifyAllKmers(12)],
        objectives=[dc.CodonOptimize(species="e_coli", method="use_best_codon",
                                     location=(0, len(cds), 1))],
        logger=None)
    problem.resolve_constraints()
    problem.optimize()
    return problem.sequence

print("seed 0 twice identical:", design(0) == design(0))
print("seed 0 vs seed 1 differ:", design(0) != design(1))

Record the seed next to the sequence. A design nobody can regenerate cannot be audited.

The trap — some constraints rewrite the sequence on construction

Most constraints are checked and then fixed by the solver. A few are enforced by restricting the mutation space instead, and DNA Chisel applies that restriction when the problem object is built. AvoidRareCodons is one of them.

import urllib.request
from Bio import SeqIO
import numpy
import dnachisel as dc

urllib.request.urlretrieve(
    "https://www.ebi.ac.uk/ena/browser/api/embl/M62653.1", "gfp.embl")
record = SeqIO.read("gfp.embl", "embl")
cds = str(next(f for f in record.features if f.type == "CDS").extract(record.seq)).upper()

def build(extra):
    numpy.random.seed(0)
    return dc.DnaOptimizationProblem(
        sequence=cds,
        constraints=[dc.EnforceTranslation(location=(0, len(cds), 1)),
                     dc.AvoidPattern("BsaI_site")] + extra,
        logger=None)

plain = build([])
restricted = build([dc.AvoidRareCodons(0.1, species="e_coli")])

for label, problem in [("without AvoidRareCodons", plain),
                       ("with AvoidRareCodons", restricted)]:
    changed = sum(1 for a, b in zip(cds, problem.sequence) if a != b)
    print("%-24s untouched=%-5s codons already rewritten=%d  BsaI site at %d"
          % (label, problem.sequence == cds, changed, problem.sequence.find("GAGACC")))

The consequence is specific. problem.constraints_text_summary() called straight after construction reports on the already-adjusted sequence, not on your input — in the run above the minus-strand BsaI site has disappeared before any solver call, making the "before" report look cleaner than the sequence you handed in. Audit the input sequence itself, as in the section above, if you want to know what was wrong with it.

When there is no solution

Over-constrain the problem and the solver raises rather than returning something plausible. Read the message: it names the specification that failed and the span it failed on.

import numpy
import dnachisel as dc

numpy.random.seed(0)
# Tryptophan has exactly one codon (TGG), so a low-GC ceiling here is unsatisfiable.
sequence = dc.reverse_translate("MWWWWWWWWWWK")

try:
    problem = dc.DnaOptimizationProblem(
        sequence=sequence,
        constraints=[dc.EnforceTranslation(location=(0, len(sequence), 1)),
                     dc.EnforceGCContent(mini=0.0, maxi=0.30, window=15)],
        logger=None)
    problem.resolve_constraints()
    print("resolved -- not expected")
except dc.NoSolutionError as error:
    print(type(error).__name__)
    print(str(error).split("\n")[0])

The span in the message is where to look first. A failure localized to one short segment is usually a genuine conflict — a run of single-codon residues, or two constraints disagreeing over the same bases — and the fix is to relax one constraint over that span rather than loosen it globally.

Export what changed, not just the answer

A bare FASTA of the optimized sequence loses the reasoning. to_record writes a GenBank in which every constraint, every objective and every edited stretch is a feature, so the design opens in any sequence viewer with its own justification attached. Adjacent edited bases merge into one feature, so the edit-segment count runs below the number of changed bases.

import urllib.request
from Bio import SeqIO
import numpy
import dnachisel as dc

urllib.request.urlretrieve(
    "https://www.ebi.ac.uk/ena/browser/api/embl/M62653.1", "gfp.embl")
record = SeqIO.read("gfp.embl", "embl")
cds = str(next(f for f in record.features if f.type == "CDS").extract(record.seq)).upper()

numpy.random.seed(0)
problem = dc.DnaOptimizationProblem(
    sequence=cds,
    constraints=[dc.EnforceTranslation(location=(0, len(cds), 1)),
                 dc.EnforceGCContent(mini=0.35, maxi=0.65, window=50),
                 dc.AvoidPattern("BsaI_site"), dc.AvoidPattern("BsmBI_site")],
    objectives=[dc.CodonOptimize(species="e_coli", method="use_best_codon",
                                 location=(0, len(cds), 1))],
    logger=None)
problem.resolve_constraints()
problem.optimize()

annotated = problem.to_record(with_sequence_edits=True)
annotated.id = "avGFP_ecoli_seed0"
annotated.annotations["molecule_type"] = "DNA"
SeqIO.write(annotated, "avGFP_ecoli.gb", "genbank")

roles = [f.qualifiers.get("role", ["edit"])[0] for f in annotated.features]
print("features written:", len(annotated.features))
print("constraints:", roles.count("c"), " objectives:", roles.count("o"),
      " edit segments:", roles.count("edit"))

Keep that GenBank, the seed, and the constraint list together. They are the three things a reviewer needs to reproduce the design, and the sequence alone is none of them.

What this does not do

  • It does not simulate the assembly. Removing BsaI sites makes a fragment Golden-Gate-compatible; it does not design the overhangs or check the junctions.
  • It does not check synthesis complexity against a specific vendor. The constraints here approximate the common rules; vendors run their own screens and the authoritative answer is their quote.
  • It does not predict expression. Codon adaptation correlates with yield and does not determine it.
  • It does not do biosecurity screening. Sequences of concern are screened by the synthesis provider, and that step is separate from and not replaced by anything here.

The Edinburgh Genome Foundry also hosts a browser interface to these methods. It still runs, but its codebase has been untouched since 2022 — reasonable for a one-off check, not something to build a pipeline against.

Try it

Data. The native Aequorea victoria green fluorescent protein mRNA, ENA accession M62653.1, fetched from the ENA Browser API. INSDC records are free to use without an account or a token; confirmed reachable 2 Oct 2026. The coding sequence is extracted from the record's own CDS feature rather than assumed, so the example survives the record's untranslated regions.

The example is routed through the trap deliberately: this CDS has no GGTCTC anywhere and still contains a BsaI site, on the minus strand at position 643. A check that greps the site as written passes it.

Run. Self-contained, from an empty directory:

pip install dnachisel biopython

import re, urllib.request
from Bio import SeqIO
from Bio.Seq import Seq
import numpy
import dnachisel as dc

urllib.request.urlretrieve(
    "https://www.ebi.ac.uk/ena/browser/api/embl/M62653.1", "gfp.embl")
record = SeqIO.read("gfp.embl", "embl")
cds = str(next(f for f in record.features if f.type == "CDS").extract(record.seq)).upper()
protein = str(Seq(cds).translate())

print("CDS %d bp, protein %d aa, GC %.3f"
      % (len(cds), len(protein) - 1, dc.biotools.gc_content(cds)))
print("GGTCTC in input:", cds.find("GGTCTC"), "| GAGACC in input:", cds.find("GAGACC"))

numpy.random.seed(0)
problem = dc.DnaOptimizationProblem(
    sequence=cds,
    constraints=[dc.EnforceTranslation(location=(0, len(cds), 1)),
                 dc.EnforceGCContent(mini=0.35, maxi=0.65, window=50),
                 dc.AvoidPattern("BsaI_site"), dc.AvoidPattern("BsmBI_site"),
                 dc.AvoidPattern("9xA"), dc.UniquifyAllKmers(12)],
    objectives=[dc.CodonOptimize(species="e_coli", method="use_best_codon",
                                 location=(0, len(cds), 1))],
    logger=None)
problem.resolve_constraints()
problem.optimize()
designed = problem.sequence

# Invariants.
assert str(Seq(designed).translate()) == protein, "protein changed"
assert len(designed) == len(cds)
for site in ["GGTCTC", "CGTCTC"]:
    assert site not in designed and dc.reverse_complement(site) not in designed, site
assert "A" * 9 not in designed
windows = [dc.biotools.gc_content(designed[i:i + 50])
           for i in range(0, len(designed) - 50 + 1)]
assert 0.35 <= min(windows) and max(windows) <= 0.65, (min(windows), max(windows))
assert len(re.findall("GAGACC", cds)) == 1, "input should carry one minus-strand BsaI site"

# Observed values -- drift to investigate, not failure.
print("edits: %d" % sum(1 for a, b in zip(cds, designed) if a != b))
print("final GC: %.4f" % dc.biotools.gc_content(designed))
print("GC window range: %.3f - %.3f" % (min(windows), max(windows)))
print("objective score: %.2f" % problem.objective_scores_sum())
print("designed tail:", designed[-24:])
print("OK")

Expect.

Invariants — a failure here means the skill is wrong, not that upstream moved:

  • The translation is byte-identical before and after, and the length is unchanged. Codon optimization that alters the protein is the one unrecoverable failure, which is why it is asserted rather than printed.
  • Neither GGTCTC/GAGACC (BsaI) nor CGTCTC/GAGACG (BsmBI) appears in the output on either strand, and no run of nine adenines survives.
  • Every 50 bp window of the output has GC between 0.35 and 0.65 inclusive.
  • The input contains exactly one GAGACC and zero GGTCTC. If that stops holding, ENA's M62653.1 has changed and the minus-strand lesson needs a new example rather than a looser assertion.

Observed on 2 Oct 2026 with dnachisel 3.2.16, python_codon_tables 0.1.18, biopython 1.88, numpy 2.4.6, Python 3.11.15, numpy.random.seed(0):

  • CDS 717 bp, protein 238 aa, GC 0.385
  • GGTCTC in input: -1 | GAGACC in input: 643
  • edits: 166, final GC: 0.4993, GC window range: 0.360 - 0.640
  • objective score: -14.64
  • designed tail: GGCATGGATGAACTGTATAAATAA

The edit count, the score and the tail are seed-dependent and version-dependent. With the seed removed they move by a few edits per run (160–166 observed); a change in them under seed 0 after a dnachisel release is drift to note, not a break.