Designing a combinatorial mutation library with GGAssembler

A combinatorial library is not a list of mutants. If you want every combination of two options at six positions, that is 64 proteins, and ordering 64 genes is absurd when the positions differ from the parent by one codon each. The economical answer is to order a few degenerate oligos — oligos synthesised with ambiguous bases, so one tube carries several sequences — and assemble them combinatorially.

Golden Gate does the assembly, and that is where the difficulty actually lives. The fragments are joined by four-base overhangs, every overhang in one reaction competes with every other, and mis-ligation between two similar overhangs silently produces the wrong protein. So the design problem is a search: cut the gene into the fewest segments such that each variable position sits inside a segment, every junction is an overhang that ligates efficiently, and no two junctions in the pot are confusable.

GGAssembler solves that search. You declare which positions vary and into what, and it returns a set of junctions and the oligos to order, scored against measured ligation data rather than a GC rule of thumb. The method is described in Hoch and colleagues, Protein Science 2024 (doi:10.1002/pro.5169), which reports one-pot camelid antibody libraries of more than 10^5 variants at under $0.007 of DNA per variant, with over 93% of the intended variants present in the product by deep sequencing.

This skill is about turning a design into an order. Getting a parent coding sequence that is clean enough to assemble — codon optimisation, GC windows, removing the assembly enzyme's own recognition sites — is the step before it, and the dna-chisel skill covers that.

What you need before the first run

• A Rust compiler. GGAssembler's shortest-path core is a Rust extension built at install time, so rustc and cargo must be on the path. There is no published wheel.
• Network access to a codon-usage service for one step only. The degenerate-table step scores candidate ambiguous codons by how often each of their codons is used in your host, and it fetches those counts from the Kazusa codon-usage database at request time. Nothing else in the pipeline touches the network. If that host is unreachable you can still build the table yourself — there is a block below that does.
• No account, no key, no GPU. The whole pipeline is CPU work and finishes in seconds for a gene-sized sequence.

GGAssembler is licensed AGPL-3.0, which permits commercial use. The network clause that distinguishes AGPL from GPL applies to offering the tool itself as a hosted service over a network; it does not reach the oligo sequences you design with it, and it does not constrain running it locally or inside your own pipeline. If you intend to expose GGAssembler to third parties as a web service, read the licence — that is the case it was written for.

# needs rustc on PATH; the extension is compiled during install
pip install "git+https://github.com/Fleishman-Lab/GGAssembler.git"

The import name is dawdlib, not ggassembler. Nothing is published on the package index under either name — a package called ggassembler there is an unrelated project by a different author, so pip install ggassembler installs the wrong software without erroring.

Four objects, in order

The pipeline is a chain of four files, and almost every confusion comes from mixing them up:

object	what it says
resfile	which amino-acid positions vary, and to which residues
degenerate codon table	which ambiguous codons realise that choice, per position
gate path	where the Golden Gate junctions go
oligo table	the DNA to order, with prefixes and suffixes attached

A resfile is the Rosetta format, one line per position:

nataa
start
65	A	PIKAA	ST
203	A	PIKAA	TY

PIKAA ST means "at position 65, allow only S or T". The parser reads the first whitespace field as the position and the last as the allowed residues; the chain column is ignored.

Step 1 — the degenerate codon table

dawdlib degenerate_table --res_file input.resfile --output_csv deg_table.csv

--organism_id is an NCBI taxonomy id and defaults to 37762 (Escherichia coli); pass 4932 for S. cerevisiae, or any other id whose codon usage the service carries. This is the step that needs the network, and when the codon-usage host cannot be reached it raises requests.exceptions.ConnectionError naming www.kazusa.or.jp rather than degrading quietly.

The output has one row per varied position:

column	meaning
AA_POS	amino-acid position, 1-based
DNA_POS	1-based position of the first base of that codon — AA_POS * 3 - 2
ENCODED_AAS	the residues the chosen codons actually encode
ENCODED_COUNT	how many codons encode each of those residues
AMBIGUOUS_CODONS1…4	the ambiguous codons for this position, blank-padded

Two of those columns carry the traps. DNA_POS counts bases from 1, so a table built by hand with 0-based offsets places every variable position one base early and the design is wrong without failing. And a position needs more than one ambiguous codon whenever the residues you want cannot be covered by a single IUPAC codon — {T, Y} differ at two of the three bases, so ACC and TAC are listed side by side and the assembly ends up with two oligos for that segment instead of one.

The library's size is the product of the per-position codon counts:

import ast
import pandas as pd

deg_table = pd.read_csv("deg_table.csv", na_filter=True, keep_default_na=False)
diversity = deg_table.ENCODED_COUNT.apply(ast.literal_eval).apply(sum).prod()
print(f"the table encodes {diversity} DNA variants")

Building the table without the codon-usage service

What the service buys is ranking — among the ambiguous codons that cover a residue set, prefer the one whose codons are common in the host. The coverage itself is arithmetic on the IUPAC code and the standard genetic code, so a table can be built offline. Use this when you want to pin the codons yourself, or when the codon-usage host is down:

import itertools
import pandas as pd
from Bio.Data.IUPACData import ambiguous_dna_values as AMBIG
from Bio.Seq import Seq

def codon_residues(codon):
    """Every residue an ambiguous codon encodes, with the number of codons for each."""
    counts = {}
    for bases in itertools.product(*[AMBIG[b] for b in codon]):
        aa = str(Seq("".join(bases)).translate())
        counts[aa] = counts.get(aa, 0) + 1
    return counts

def exact_codons(residues):
    """Single ambiguous codons encoding exactly this residue set and no stop."""
    want = frozenset(residues)
    hits = []
    for codon in map("".join, itertools.product("ACGTRYSWKMBDHVN", repeat=3)):
        counts = codon_residues(codon)
        if "*" not in counts and frozenset(counts) == want:
            hits.append(codon)
    return hits

print("S,T ->", exact_codons("ST")[:4])
print("T,Y ->", exact_codons("TY"))        # empty: needs two codons

def deg_table_from_design(design):
    """design maps a 1-based AA position to the list of ambiguous codons to use there."""
    rows = []
    for aa_pos, codons in sorted(design.items()):
        counts = {}
        for codon in codons:
            for aa, n in codon_residues(codon).items():
                counts[aa] = counts.get(aa, 0) + n
        if "*" in counts:
            raise ValueError(f"position {aa_pos} encodes a stop codon")
        residues = sorted(counts)
        row = {"AA_POS": aa_pos, "DNA_POS": aa_pos * 3 - 2,
               "ENCODED_AAS": str(residues),
               "ENCODED_COUNT": str([counts[aa] for aa in residues])}
        for i in range(4):
            row[f"AMBIGUOUS_CODONS{i + 1}"] = codons[i] if i < len(codons) else ""
        rows.append(row)
    return pd.DataFrame(rows)

exact_codons returning an empty list is the signal that the position needs a multi-codon entry, not that the design is impossible.

Step 2 — the parent gene must carry no site for the assembly enzyme

Golden Gate cuts with a type IIS enzyme, so a recognition site anywhere inside the parent gene is cut as well and the assembly loses a fragment. GGAssembler checks this for you and hands back which positions offend:

from dawdlib.golden_gate.utils import check_for_restriction_sites

clean, enzyme, positions = check_for_restriction_sites(dna, ["BsaI"])
if not clean:
    raise SystemExit(f"{enzyme} site(s) at {positions} — remove them before designing")

The check runs on both strands — BsaI is GGTCTC, and a gene with no GGTCTC anywhere can still carry GAGACC on the minus strand and be cut. Grepping the site as written is the mistake this check exists to catch. Clearing it is a synonymous codon substitution; dna-chisel does it as a declared constraint, and the ## Try it block below does it in a few lines for the one site in the example gene.

Note that a chosen degenerate codon can also create a site that the parent did not have. Nothing catches that until the reaction simulation in step 4, which is one reason not to skip it.

Step 3 — placing the junctions

Use the Python API, not the create_goldengates subcommand. At master 1bfee36 the dawdlib create_goldengates and dawdlib ggsim subcommands both raise

TypeError: Requirements.__new__() missing 2 required positional arguments: 'oligo_suffix' and 'const_cost'

because the Requirements tuple gained fields that those two call sites were never updated for. The degenerate_table, combine and embl subcommands are unaffected. The library path below is what the project's own example notebook uses and it works.

from dawdlib.golden_gate.utils import RequirementsFactory, expand_dna_var_poss
from dawdlib.golden_gate.gate_data import GGData
from dawdlib.golden_gate.graph_maker import GraphMaker, make_default_graph
from dawdlib.golden_gate.find_gg import deg_table_to_dict
from dawdlib.dijkstra.len_limit import all_shortest_paths
from dawdlib.degenerate_dna.deg_table import TableColNames

PREFIX = "GACATTGGTCTCA"                        # ends with the BsaI site
SUFFIX = "TGAGACCAACGACGCCGTACTCTTTGTCAAC"      # starts with the BsaI site

reqs = RequirementsFactory(
    min_oligo_length=20, max_oligo_length=200, min_const_oligo_length=20,
    min_efficiency=0.25, min_fidelity=0.1,
    oligo_prefix=PREFIX, oligo_suffix=SUFFIX,
    const_cost=40, filter_gc_overhangs=False,
)
ggdata = GGData(temperature=reqs.gg_temp, hours=reqs.gg_hours,
                min_efficiency=reqs.min_efficiency, min_fidelity=reqs.min_fidelity)

var_poss = expand_dna_var_poss(deg_table[TableColNames.DNA_POS.value].tolist())
graph, src, target = make_default_graph(
    GraphMaker(ggdata), dna, var_poss, deg_table_to_dict(deg_table), reqs, 4)
print(f"{graph.number_of_nodes()} candidate junctions, {graph.number_of_edges()} segments")

Four things in that call decide whether you get a usable answer.

• RequirementsFactory, not Requirements. The factory subtracts the prefix and suffix lengths from max_oligo_length, so the number you pass is the length of the oligo you will actually order. Construct the tuple directly and you will order oligos that are 44 bases longer than your synthesis limit.
• max_oligo_length bounds only the segments that carry a varied position. A segment with nothing variable in it has a minimum (min_const_oligo_length) and no maximum at all, so constant segments routinely come back longer than the synthesis limit. That is the design working as intended — you amplify those off the parent plasmid rather than buying them — but it means a blanket length check over the whole oligo table fails, and the check you want is over the non-constant rows. Both bounds are strict, so a segment exactly max_oligo_length long is rejected.
• min_fidelity is a mis-ligation budget, not a fidelity floor in the usual direction. It is compared against the reaction fidelity returned by the ligation model; 0.1 is permissive and what the project's example uses.
• The last argument is the overhang length, 4 for BsaI and BsmBI. SapI leaves 3-base overhangs, and switching to it means passing 3 and swapping the ligation table and the prefix and suffix for that enzyme's recognition sequence.
• PREFIX and SUFFIX must contain the enzyme's recognition site, because they are what the enzyme cuts to release each fragment. They are not primer tails.

Then enumerate shortest paths and keep the first whose overhang set ligates cleanly:

chosen = None
for i, (path, cost) in enumerate(all_shortest_paths(
        graph, src, target, weight="weight", len_cutoff=12)):
    junctions = [g for g in path if not g.src_or_target]
    try:
        fidelity = ggdata.reaction_fidelity(*[g.bps for g in junctions])[0]
    except ValueError:
        continue                      # overhang pair absent from the ligation table
    if fidelity > reqs.min_fidelity:
        chosen = (path, cost, fidelity)
        break
    if i > 2000:
        break

path, cost, fidelity = chosen
print(f"{len([g for g in path if not g.src_or_target])} junctions, "
      f"cost {cost}, reaction fidelity {fidelity:.4f}")

cost is the thing being minimised and it is in bases of synthesis, weighted by const_cost for each constant segment — that is what "economical" means here, and it is why the search prefers few long segments over many short ones.

Ligation data and reaction conditions

The fidelity model is a published measurement, not a heuristic — the T4 ligase four-base-overhang ligation profiling of Potapov and colleagues (doi:10.1021/acssynbio.8b00333), shipped with the package as CSV tables. That is what makes "this overhang set ligates cleanly" a claim with numbers behind it. Only three condition combinations are bundled, and GGData raises ValueError for anything else:

import dawdlib.golden_gate.resources as resources
from dawdlib.golden_gate.gate_data import GGData

print(sorted(resources.ligation_data))
print({k: v.rsplit("/", 1)[-1] for k, v in GGData.df_files.items()})   # (hours, °C)
try:
    GGData(temperature=16, hours=4)
except ValueError as exc:
    print("unsupported conditions ->", exc)

Watch the defaults: GGData defaults to 25 °C while Requirements defaults gg_temp to 37, so passing reqs.gg_temp through — as above — is not the same as taking either default. Set the pair to the conditions you will actually run, and say which in the design summary.

Step 4 — oligos, and simulating the assembly

from dawdlib.gg_dc_combine.gg_dc_combine import dc_df_codon_list, gate_cdn_oligos
from dawdlib.golden_gate.reaction_sim import ReactionSim

oligos = gate_cdn_oligos(path, dc_df_codon_list(deg_table), dna,
                         reqs.oligo_prefix, reqs.oligo_suffix, "GFP-SPECTRAL")
oligos.to_csv("oligos.csv", index=False)

sim = ReactionSim(ggdata, reqs, ["BsaI"])
problem = sim.create_reaction_graph("oligos.csv")
if problem is not None:
    message, entry = problem
    raise SystemExit(f"{message}\ncodons {entry.oligo_codons} in {entry.name} "
                     f"created a new restriction site — change them in the table")

wt, = sim.get_wt_dna()
assert wt.dna == dna, "the simulated parent product is not the input gene"
ok, products, offender = sim.verify_reaction(
    wt.end - wt.start, wt.no_segments, wt.fidelity_sum)
print(f"simulation {'passed' if ok else 'FAILED'}, {products} distinct products")
if not ok:
    print("first badly assembled product:", offender)

verify_reaction returns three values — the verdict, the product count, and the first product that failed, empty on success. Its own docstring describes a 2-tuple, so unpacking into two names raises ValueError: too many values to unpack rather than anything that points at the cause.

verify_reaction is the check that matters. It walks the assembly graph and counts the distinct products the reaction can make; that number must equal the diversity the degenerate table encodes. If it is lower, diversity was lost somewhere — a segment that cannot ligate, an overhang collision, a codon that got absorbed. If it is higher, something is assembling that should not. Either way the order is wrong and the design has to change, not the assertion.

create_reaction_graph returning None means no oligo created a new restriction site. When it returns a tuple, the fix is in the degenerate table — pick a different synonymous codon for the position it names and re-run from step 3. Do not relax the enzyme list.

The oligo table has one row per segment, and two boolean columns decide what you buy:

• wt marks the parent sequence for that segment, const marks a segment with no varied position in it.
• The rows that carry diversity are ~wt & ~const — those are the degenerate oligos.
• The wt & const rows are the invariant backbone pieces. You need them for the assembly, but only once, you may already have them from a previous order, and they are the rows that can exceed your synthesis limit — amplify them from the parent.

full_oligo_dna is what you paste into an order form — segment plus prefix and suffix, with the ambiguous bases in place. oligo_dna is the segment alone.

Step 5 — annotate the design so a human can read it

from Bio import SeqIO
from dawdlib.embl_utils.embl_maker import create_dc_features, create_path_features
from dawdlib.degenerate_dna.utils import parse_degenerate_codon_csv

record = SeqIO.read("parent.fasta", "fasta")
record.annotations["molecule_type"] = "DNA"
record.features.extend(create_dc_features(parse_degenerate_codon_csv("deg_table.csv")))
record.features.extend(create_path_features(path))
SeqIO.write(record, "design.embl", "embl")

That file opens in any sequence viewer with the varied codons and the junctions marked, which is the artefact to attach to a design review. molecule_type must be set or Biopython refuses to write EMBL.

Expect one BiopythonWarning per feature on the way out. GGAssembler names each feature after what it is — ST=1 for a degenerate codon, G.178.TTGT for a junction — and those are not valid INSDC feature keys. The file is written and reads correctly in a viewer; the warnings are not a failure and there is nothing to fix in the design.

What to hand back

The deliverable is the order table plus a short design summary — not a transcript of the run. Keep them separate, and never leave a placeholder in either.

• Proposed — the junctions chosen, the segment count, and the oligos to order with their lengths.
• Inferred — the diversity the table encodes and the product count the simulation found, stated as the same number or flagged as a mismatch.
• Decided — enzyme, overhang length, reaction temperature and time, synthesis length limit, and whether constant segments were included in the order.
• Outstanding — anything left to a human. A position whose residue set needed two codons and so doubled a segment's oligo count is the usual one.

Try it

Data. The native Aequorea victoria green fluorescent protein mRNA, ENA accession M62653.1, fetched from the ENA Browser API. INSDC records carry no access conditions and need no account; confirmed reachable 14 Oct 2026. The coding sequence is taken from the record's own CDS feature rather than assumed, and the design varies six positions that avGFP numbering makes recognisable — 64, 65, 145, 148, 167 and 203 — two options each, so 64 variants.

The block is routed through two traps on purpose. The native CDS carries no GGTCTC anywhere and is still cut by BsaI, on the minus strand, so it must be repaired before it can be assembled. And position 203 wants {T, Y}, which no single ambiguous codon encodes, so that segment ends up with two oligos while the others have one.

The degenerate table is built inline by the offline routine above rather than by dawdlib degenerate_table, so this runs with no network beyond the one ENA fetch. That is a deliberate substitution and it changes nothing that is checked here — the codons are pinned by the design, and what the codon-usage service would have supplied is a preference among codons that encode the same residues.

Run. Self-contained, from an empty directory. Needs rustc on the path.

pip install "git+https://github.com/Fleishman-Lab/GGAssembler.git"

import ast
import itertools
import urllib.request

import pandas as pd
from Bio import Restriction, SeqIO
from Bio.Data.IUPACData import ambiguous_dna_values as AMBIG
from Bio.Seq import Seq

from dawdlib.degenerate_dna.deg_table import TableColNames
from dawdlib.dijkstra.len_limit import all_shortest_paths
from dawdlib.gg_dc_combine.gg_dc_combine import dc_df_codon_list, gate_cdn_oligos
from dawdlib.golden_gate.find_gg import deg_table_to_dict
from dawdlib.golden_gate.gate_data import GGData
from dawdlib.golden_gate.graph_maker import GraphMaker, make_default_graph
from dawdlib.golden_gate.reaction_sim import ReactionSim
from dawdlib.golden_gate.utils import (RequirementsFactory,
                                       check_for_restriction_sites,
                                       expand_dna_var_poss)

# ---- parent gene: native avGFP from ENA
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(f"CDS {len(cds)} bp, protein {len(protein) - 1} aa")

# Trap 1: no GGTCTC in the text, still a BsaI site -- on the minus strand.
assert "GGTCTC" not in cds
print("BsaI cut positions in the native CDS:", Restriction.BsaI.search(Seq(cds)))
print("check_for_restriction_sites(native):", check_for_restriction_sites(cds, ["BsaI"]))

# Repair it with synonymous substitutions.
synonyms = {}
for codon in map("".join, itertools.product("ACGT", repeat=3)):
    synonyms.setdefault(str(Seq(codon).translate()), []).append(codon)

dna = cds
while Restriction.BsaI.search(Seq(dna)):
    site = next(p for p in (str(Restriction.BsaI.site),
                            str(Seq(str(Restriction.BsaI.site)).reverse_complement()))
                if p in dna)
    start = dna.find(site)
    for ci in range(start // 3, (start + len(site) - 1) // 3 + 1):
        here = dna[ci * 3:ci * 3 + 3]
        swap = next((alt for alt in synonyms[str(Seq(here).translate())]
                     if alt != here
                     and not Restriction.BsaI.search(
                         Seq(dna[:ci * 3] + alt + dna[ci * 3 + 3:]))), None)
        if swap:
            dna = dna[:ci * 3] + swap + dna[ci * 3 + 3:]
            break
    else:
        raise SystemExit("no synonymous substitution clears the site")

# Invariants: same protein, same length, no site on either strand.
assert str(Seq(dna).translate()) == protein, "repair changed the protein"
assert len(dna) == len(cds)
assert check_for_restriction_sites(dna, ["BsaI"])[0]
print("synonymous edits to clear BsaI:",
      sum(1 for a, b in zip(cds, dna) if a != b))
open("parent.fasta", "w").write(f">parent\n{dna}\n")

# ---- degenerate codon table, built offline
def codon_residues(codon):
    counts = {}
    for bases in itertools.product(*[AMBIG[b] for b in codon]):
        aa = str(Seq("".join(bases)).translate())
        counts[aa] = counts.get(aa, 0) + 1
    return counts

DESIGN = {64: ["TTM"], 65: ["AST"], 145: ["TWC"],
          148: ["YAC"], 167: ["AYT"], 203: ["ACC", "TAC"]}   # trap 2: 203 needs two

rows = []
for aa_pos, codons in sorted(DESIGN.items()):
    counts = {}
    for codon in codons:
        for aa, n in codon_residues(codon).items():
            counts[aa] = counts.get(aa, 0) + n
    assert "*" not in counts, aa_pos
    residues = sorted(counts)
    assert protein[aa_pos - 1] in residues, f"parent residue dropped at {aa_pos}"
    row = {"AA_POS": aa_pos, "DNA_POS": aa_pos * 3 - 2,
           "ENCODED_AAS": str(residues),
           "ENCODED_COUNT": str([counts[aa] for aa in residues])}
    for i in range(4):
        row[f"AMBIGUOUS_CODONS{i + 1}"] = codons[i] if i < len(codons) else ""
    rows.append(row)
deg_table = pd.DataFrame(rows)
deg_table.to_csv("deg_table.csv", index=False)

# Invariant: DNA_POS is the 1-based first base of the codon for that residue.
for row in rows:
    codon = dna[row["DNA_POS"] - 1:row["DNA_POS"] + 2]
    assert str(Seq(codon).translate()) == protein[row["AA_POS"] - 1], row

diversity = deg_table.ENCODED_COUNT.apply(ast.literal_eval).apply(sum).prod()
print("intended diversity:", diversity)

# ---- junctions
PREFIX = "GACATTGGTCTCA"
SUFFIX = "TGAGACCAACGACGCCGTACTCTTTGTCAAC"
reqs = RequirementsFactory(
    min_oligo_length=20, max_oligo_length=200, min_const_oligo_length=20,
    min_efficiency=0.25, min_fidelity=0.1,
    oligo_prefix=PREFIX, oligo_suffix=SUFFIX,
    const_cost=40, filter_gc_overhangs=False)
assert reqs.max_oligo_length == 200 - len(PREFIX) - len(SUFFIX)

ggdata = GGData(temperature=reqs.gg_temp, hours=reqs.gg_hours,
                min_efficiency=reqs.min_efficiency, min_fidelity=reqs.min_fidelity)
graph, src, target = make_default_graph(
    GraphMaker(ggdata), dna,
    expand_dna_var_poss(deg_table[TableColNames.DNA_POS.value].tolist()),
    deg_table_to_dict(deg_table), reqs, 4)
print(f"graph: {graph.number_of_nodes()} nodes, {graph.number_of_edges()} edges")

chosen = None
for i, (candidate, cost) in enumerate(all_shortest_paths(
        graph, src, target, weight="weight", len_cutoff=12)):
    junctions = [g for g in candidate if not g.src_or_target]
    try:
        fidelity = ggdata.reaction_fidelity(*[g.bps for g in junctions])[0]
    except ValueError:
        continue
    if fidelity > reqs.min_fidelity:
        chosen = (candidate, cost, fidelity, len(junctions))
        break
    if i > 2000:
        break
assert chosen, "no path cleared the fidelity threshold"
path, cost, fidelity, n_junctions = chosen

# ---- oligos and simulation
oligos = gate_cdn_oligos(path, dc_df_codon_list(deg_table), dna,
                         reqs.oligo_prefix, reqs.oligo_suffix, "GFP-SPECTRAL")
oligos.to_csv("oligos.csv", index=False)

sim = ReactionSim(ggdata, reqs, ["BsaI"])
problem = sim.create_reaction_graph("oligos.csv")
assert problem is None, f"an oligo created a new restriction site: {problem}"
wt, = sim.get_wt_dna()
assert wt.dna == dna, "simulated parent product is not the input gene"
ok, products, offender = sim.verify_reaction(
    wt.end - wt.start, wt.no_segments, wt.fidelity_sum)

# Invariants: the reaction makes exactly the library the table encodes.
assert ok, f"verify_reaction failed on {offender}"
assert products == diversity, (products, diversity)
to_order = oligos[~oligos.wt & ~oligos.const]
assert len(to_order) == 4, len(to_order)
# The length bound applies to segments carrying a varied position, not to constant ones.
assert all(len(o) <= 200 for o in oligos[~oligos.const].full_oligo_dna), \
    "a variable oligo exceeds the synthesis limit"
assert any(len(o) > 200 for o in oligos[oligos.const].full_oligo_dna), \
    "expected at least one constant segment past the limit -- those are amplified"

# Observed values -- drift to investigate, not failure.
print(f"junctions: {n_junctions}, cost: {cost}, fidelity: {fidelity:.4f}")
print(f"segments in the reaction: {wt.no_segments}")
print(f"oligo rows: {len(oligos)}, degenerate oligos to order: {len(to_order)}")
print("overhangs:", [g.bps for g in path if not g.src_or_target])
print(f"products: {products} (intended {diversity})")
print("OK")

Expect.

Invariants — a failure here means the skill or the design is wrong:

• The native CDS contains no GGTCTC and Restriction.BsaI.search still reports a cut.
• The repaired gene translates to the same 238-residue protein and is the same length.
• DNA_POS - 1 indexes the codon for AA_POS in the repaired gene.
• create_reaction_graph returns None.
• The simulated parent product equals the input gene, and verify_reaction reports exactly the diversity the table encodes.
• No oligo carrying a varied position is longer than the 200-base limit passed to RequirementsFactory, and at least one constant segment is longer — the bound does not apply to those.

Observed on 14 Oct 2026, against GGAssembler at master 1bfee36 — a mismatch here is drift to investigate, not a bug:

CDS 717 bp, protein 238 aa
BsaI cut positions in the native CDS: [639]
check_for_restriction_sites(native): (False, 'BsaI', [639])
synonymous edits to clear BsaI: 1
intended diversity: 64
graph: 609 nodes, 58326 edges
junctions: 5, cost: 460, fidelity: 0.9778
segments in the reaction: 8
oligo rows: 10, degenerate oligos to order: 4
overhangs: ['TTGT', 'TATG', 'CAAC', 'CTTT', 'CAAT']
products: 64 (intended 64)
OK

Sixty-four proteins from four degenerate oligos and six backbone segments is the point of the method, and the product count is the only evidence that the library you ordered is the library you designed.