Enforce design requirements and fix optimizer correctness

.github/workflows/ci.yml

@@ -39,8 +39,8 @@ jobs:
       - uses: actions/setup-python@v6
         with:
           python-version: "3.12"
-      # The default splice backend prefers MaxEntScan when it is installed, so
-      # without this job that path is never exercised in CI.
+      # Explicit/automatic human splice screening uses MaxEntScan when it is
+      # installed; exercise that optional backend as well as dependency-free runs.
       - run: pip install -e ".[dev]" "git+https://github.com/kepbod/maxentpy.git"
       - name: Confirm MaxEntScan really is in use
         run: python -c "from codonopt.splice import maxent_available; assert maxent_available()"

README.md

@@ -23,15 +23,17 @@
 </p>
 
 Give it a protein, or an existing coding sequence, and it returns a synonymous
-sequence that balances GC content, repeats, hairpins, homopolymers, forbidden
-motifs, cryptic splice sites and host codon usage in a single score. It is a
+sequence that resolves mandatory requirements before balancing GC content,
+repeats, hairpins, homopolymers, candidate splice signals and host codon usage.
+Avoided motifs are mandatory by default; other criteria are soft preferences
+unless explicitly required. It is a
 pure-Python package with no runtime dependencies, a command line, a Python API
 and a local web interface.
 
 ## Highlights
 
-- **One objective, not four passes.** Every concern is a weighted penalty in one score, so fixing one problem cannot silently undo another.
-- **Fast.** Scoring is incremental: a codon swap re-scores only what sits near it. An 800-residue protein takes 19 s, against 234 s for the scripts this replaces.
+- **Requirements before preferences.** Mandatory violations take precedence over the weighted score. Every result is independently checked and marked valid or unresolved; an unresolved result is not proof that no solution exists.
+- **Incremental scoring.** A codon swap re-scores nearby terms and updates repeat/hairpin indexes. Each edit still copies the DNA string; search cost is not strictly independent of gene length.
 - **The protein cannot change.** The only move is swapping a codon for a synonymous one, so correctness holds by construction, not by checking afterwards.
 - **Three ways in.** `codonopt` on the command line, a small Python API, and a browser or desktop interface.
 - **Checked against brute force.** The tests compare the incremental score with a from-scratch recomputation and with simple reference implementations.
@@ -54,24 +56,26 @@ codonopt optimize --add-stop --avoid golden_gate,EcoRI --seed 5 --splice-backend
 The sequence goes to stdout as FASTA (or `--out FILE`), the summary to stderr:
 
 ```
->protein codonopt table=human score=1.17
+>protein codonopt table=human score=1.08 valid=true
 ATGAAATGGGTGACATTCATCAGCTTACTGCTGCTGTTCTCTTCCGCTTACTCCCGGGGC
-GTGTTTAGACGGGACACTCATAAATCTGAGATCGCTCACAGGTTTAAGGACCTGGGCGAG
-GAGAACTTCAAGGCCTTGGTTTTAATCGCCTTCGCACAGTATCTGCAGCAGTGCCCTTTC
+GTGTTTAGACGGGACACACACAAATCTGAGATCGCTCACAGGTTTAAGGACCTGGGCGAG
+GAGAACTTCAAGGCCTTGGTTTTAATCGCCTTTGCACAGTATCTGCAGCAGTGCCCTTTC
 TGA
 
-Objective           54.44 -> 1.17  (improved by 53.28)
+Requirements        valid
+Objective           53.57 -> 1.08  (improved by 52.49)
 
 Property                      before       after
 ------------------------------------------------
 GC fraction                    0.481       0.503
-Codon adaptation index         0.774       0.837
+Codon adaptation index         0.774       0.842
 Forbidden motif hits               1           0
-Splice donor sites                 2           0
+Splice donor sites                 1           0
 Splice acceptor sites              1           0
 ```
 
 `--json FILE` writes the same numbers, the per-term score and the sequence.
+Reports also include validity and the effective settings.
 
 ## Web interface
 
@@ -84,6 +88,10 @@ table and what to keep out, and press **Optimize**. You get:
 - a before/after table, a GC plot along the sequence and the score breakdown;
 - copy and download for the DNA, FASTA and a report; runs can be cancelled and keep the best result so far.
 
+If mandatory requirements remain, the result is labelled **Unresolved candidate**.
+Downloading it does not certify it as a valid design. Splice screening is off
+by default; select it explicitly when relevant to the expression host.
+
 It is local only: the server binds to loopback, checks Host and Origin, requires
 a per-run token, sends a strict content security policy and never touches your
 files. Add `--no-window --no-browser` to just serve it.
@@ -92,13 +100,13 @@ files. Add `--no-window --no-browser` to just serve it.
 
 | Criterion | Why | Default |
 | --- | --- | --- |
-| Overall GC | Extreme GC complicates synthesis and PCR | target 0.50 |
+| Overall GC | Extreme GC complicates synthesis and PCR | target 0.50; length-scaled soft cost |
 | Local GC | Stretches matter, not just the average | 30 bp windows within 0.30–0.70 |
 | Direct repeats | Mis-assembly in synthesis, recombination in vivo | no repeated 12-mers |
 | Hairpins | Stall synthesis and reverse transcription | no 10 bp inverted repeat, 3–60 nt loop |
 | Homopolymers | Slippage in synthesis and sequencing | at most 6 identical bases |
-| Forbidden motifs | Cloning sites, poly(A) signals, anything you name | none |
-| Cryptic splice sites | Unwanted splicing of a transgene | MaxEntScan or a consensus screen |
+| Forbidden motifs | Cloning sites, poly(A) signals, anything you name | none supplied; mandatory when supplied |
+| Candidate splice signals | Screen for potential unwanted splicing | off; opt into MaxEntScan or a consensus screen |
 | Codon usage | Rare codons slow translation | host codon adaptation index |
 
 ## Command line
@@ -116,27 +124,67 @@ codonopt motifs                                             # known motifs and g
 | `--avoid SPEC` | Exclude motifs: a group (`golden_gate`, `common_cloning`, `regulatory`), an enzyme (`EcoRI`), an IUPAC pattern (`GGWCC`) or `NAME=PATTERN`. Both strands are checked. |
 | `--table NAME` / `--table-file F` | Host codon usage: `human`, `ecoli`, `uniform`, or your own JSON/CSV. |
 | `--weight NAME=VALUE` | Re-balance a term: `gc_global=200` holds the GC target harder, `usage=0` ignores codon preference. |
-| `--splice-backend` | `auto` (MaxEntScan if installed, else a consensus screen), `maxent`, `motif`, `none`. |
+| `--splice-backend` | Default `none`. `auto` screens only the bundled human host, using MaxEntScan if installed, otherwise a consensus screen. `maxent` and `motif` are explicit overrides for other hosts; neither is automatically validated for them. |
+| `--splice-strand` | Default `forward`, the expressed transcript. Select `reverse` or `both` only when appropriate. Restriction motifs still use both DNA strands. |
+| `--require NAME` | Make `gc_window`, `repeat`, `hairpin`, `homopolymer` or `splice` mandatory. Repeat the option for several criteria. Required splice screening needs an explicit backend. |
+| `--soft-motifs` | Make motif avoidance a preference instead of a requirement. |
+| `--left-flank SEQ` / `--right-flank SEQ` | Score immutable construct context, including features across CDS junctions. Output FASTA remains the editable CDS. |
+| `--protect-codon INDEX` | Preserve a zero-based codon of supplied starting DNA, including across restarts. Repeat for several positions. |
+| `--exclude-codon CODON` | Forbid a codon explicitly; a zero frequency in a usage table is not an exclusion. |
 | `--seed N` | Reproducible runs. |
 | `--restarts N`, `--sweeps N`, `--anneal N` | More search effort. |
 | `--seconds N` | Wall-clock cap (default 120). |
+| `--repair-candidates N` / `--repair-max-codons N` | Bound coordinated repair per restart (defaults 2000 complete candidates and up to 3 codons; maximum 4). |
+
+The CLI returns **0** when configured requirements are met, **2** when they
+remain unresolved, and **1** for invalid input. It still emits the best
+candidate and report on exit 2, including in quiet mode. A valid result can
+have a higher soft score than its input because validity takes precedence.
+Cancellation and time limits are recorded separately from validity.
 
 `codonopt optimize --help` lists the rest. MaxEntScan is optional and not on
 PyPI: `pip install "git+https://github.com/kepbod/maxentpy.git"`.
 
 ## Python API
 
 ```python
-from codonopt import DesignConfig, OptimizerSettings, get_table, optimize, resolve_motifs
+from codonopt import DesignConfig, OptimizerSettings, get_table, optimize, optimize_dna, resolve_motifs
 
 config = DesignConfig(table=get_table("human"), motifs=tuple(resolve_motifs(["golden_gate"])))
 result = optimize("MKWVTFISLLLLFSSAYS", config, OptimizerSettings(seed=1))
 result.sequence, result.score, result.breakdown.terms
+result.valid, result.unresolved_constraints, result.stop_reason
 ```
 
 `measure(sequence, config)` returns what `report` prints; `optimize_dna()`
 improves an existing sequence.
 
+For an existing construct with fixed context and mandatory local GC:
+
+```python
+config = DesignConfig(
+    motifs=tuple(resolve_motifs(["golden_gate"])),
+    required=("motif", "gc_window"),
+    gc_step=1,  # every base-offset window, when the synthesis policy needs it
+    left_flank="ACGT", right_flank="TGCA",
+    protected_codons=(0,),
+)
+result = optimize_dna("ATGGGTCTCAAA", config, OptimizerSettings(seed=5))
+```
+
+All context-dependent criteria scan the **complete supplied construct**,
+including fixed-only regions. A violation confined to a fixed flank remains
+unresolved. Codon usage and overall GC target apply to the CDS. Diagnostic
+positions and GC profiles use construct coordinates; `cds_start` records the
+offset. Flanks need not have lengths divisible by three.
+
+JSON reports preserve original input/first seed, winning restart seed, full
+effective configuration, actual codon-table data and hash, generated or supplied
+seed, search settings, software/backend versions, validity and stop reason.
+Custom scorer implementations and uninstalled backend sources additionally
+need their source revision recorded externally. A seed and deterministic
+proposal budgets replay a run; wall-clock truncation can differ by machine.
+
 ## How it works
 
 Each penalty term is tied to a position (a GC window, a homopolymer, a motif
@@ -145,12 +193,21 @@ it. Repeats and hairpins are not local, so an incrementally updated k-mer index
 carries them. A candidate swap costs time proportional to the window sizes, not
 the sequence length.
 
-The search sweeps every codon greedily until nothing improves, optionally
-followed by simulated annealing and several restarts.
+The search compares candidates by `(mandatory violation count, soft score)`.
+It sweeps codons greedily, uses bounded coordinated repair around unresolved
+regions, and optionally follows with simulated annealing and several restarts.
+Annealing cannot worsen mandatory violation count. Intermediate edits in a
+multi-codon repair are temporary; the complete candidate decides acceptance.
 `SequenceState.audit()` recomputes a score from scratch; the tests check the
 incremental path against it and against brute-force implementations in
 `codonopt.metrics`.
 
+Overall GC cost is `weight * abs(GC bases - target * CDS bases) / 300`.
+This retains the old scale at 100 codons while making each GC-base change
+carry the same tradeoff at different gene lengths. Local GC, repeats and the
+other occurrence-based costs remain sums. The default target is a design
+preference, not a universal expression optimum.
+
 ## Development
 
 The Python package needs nothing but Python. The interface is written in
@@ -191,6 +248,16 @@ benchmark are in [`legacy/README.md`](legacy/README.md).
   for production work. The consensus splice screen is a filter, not a model.
 - The search is local. A given `--seed` is reproducible, but another seed may do
   slightly better; use `--restarts` or `--anneal` when it counts.
+- Bounded coordinated repair can still leave feasible requirements unresolved.
+  The tool does not claim mathematical infeasibility from a search failure.
+- A valid result meets the **configured checks**, not an unprovided synthesis
+  vendor policy or a guarantee of improved expression. CAI and exact inverted
+  repeats are proxies, not validated predictions of expression or RNA folding.
+- The time limit covers search checks. Initial scoring and final independent
+  validation/reporting must still complete and can add time beyond that budget.
+- Version 0.3 changes motif validity, splice defaults/strand policy and GC cost
+  scaling. Previous seeds can produce different designs. Before/after DNA
+  reports now always use the original DNA, even when a later restart wins.
 
 ## License
 

benchmarks/compare_legacy.py

@@ -66,6 +66,8 @@ def main() -> int:
         help="skip the end-to-end comparison, which is the slow part",
     )
     args = parser.parse_args()
+    # The legacy functions use the module-level generator, not the local RNG.
+    random.seed(args.seed)
 
     splice_legacy = load_legacy(LEGACY_SPLICE)
     hairpin_legacy = load_legacy(LEGACY_HAIRPIN)

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "codonopt"
-version = "0.2.0"
+version = "0.3.0"
 description = "Codon-optimize an open reading frame for synthesis success and expression"
 readme = "README.md"
 requires-python = ">=3.10"

src/codonopt/__init__.py

@@ -2,9 +2,8 @@
 
 A coding sequence can be written in many ways without changing the protein it
 encodes, and some of those ways are much easier to synthesise and express than
-others. This package searches that space under one weighted objective covering
-GC content, direct repeats, hairpins, homopolymer runs, forbidden motifs,
-cryptic splice sites and host codon preference.
+others. This package resolves configured mandatory requirements before
+improving soft sequence-quality and host-codon-preference objectives.
 
 Typical use::
 
@@ -23,9 +22,10 @@
 
 from __future__ import annotations
 
-__version__ = "0.2.0"
+__version__ = "0.3.0"
 
 from .config import DesignConfig, OptimizerSettings, Weights
+from .constraints import constraint_violations
 from .index import KmerIndex
 from .metrics import (
     count_hairpin_pairs,
@@ -65,6 +65,7 @@
     "SequenceState",
     "Weights",
     "available_tables",
+    "constraint_violations",
     "count_hairpin_pairs",
     "find_hairpins",
     "find_splice_sites",

src/codonopt/cli.py

@@ -195,9 +195,9 @@ def _add_constraint_arguments(
     )
     group.add_argument(
         "--splice-backend",
-        default="auto",
+        default="none",
         choices=("auto", "maxent", "motif", "none"),
-        help="cryptic splice site scoring (default: auto)",
+        help="candidate splice signal screening (default: none; auto is human-only)",
     )
     group.add_argument(
         "--donor-cutoff",
@@ -211,6 +211,17 @@ def _add_constraint_arguments(
         metavar="F",
         help="highest acceptable acceptor score (default: the backend's)",
     )
+    group.add_argument("--splice-strand", choices=("forward", "reverse", "both"), default="forward")
+    group.add_argument("--require", action="append", default=[],
+                       choices=("motif", "gc_window", "repeat", "hairpin", "homopolymer", "splice"),
+                       help="make a criterion mandatory; repeat to require several")
+    group.add_argument("--soft-motifs", action="store_true", help="treat avoided motifs as preferences instead of requirements")
+    group.add_argument("--left-flank", default="", help="immutable upstream DNA used when scoring the construct")
+    group.add_argument("--right-flank", default="", help="immutable downstream DNA used when scoring the construct")
+    group.add_argument("--protect-codon", type=int, action="append", default=[], metavar="INDEX",
+                       help="protect a zero-based codon index in starting DNA; repeat as needed")
+    group.add_argument("--exclude-codon", action="append", default=[], metavar="CODON",
+                       help="forbid a codon independently of its frequency")
     if weights:
         group.add_argument(
             "--weight",
@@ -250,6 +261,8 @@ def _add_search_arguments(parser: argparse.ArgumentParser) -> None:
     group.add_argument(
         "--seed-mode", default="weighted", choices=("weighted", "best", "random")
     )
+    group.add_argument("--repair-candidates", type=int, default=2000)
+    group.add_argument("--repair-max-codons", type=int, default=3, choices=(2, 3, 4))
 
 
 class _UsageError(Exception):
@@ -261,7 +274,7 @@ def _build_config(args: argparse.Namespace) -> DesignConfig:
         load_table(args.table_file) if args.table_file else get_table(args.table)
     )
     try:
-        scorer = get_scorer(args.splice_backend)
+        scorer = get_scorer(args.splice_backend, host=table.name)
     except RuntimeError as exc:
         raise _UsageError(str(exc)) from None
 
@@ -286,6 +299,12 @@ def _build_config(args: argparse.Namespace) -> DesignConfig:
         donor_cutoff=args.donor_cutoff,
         acceptor_cutoff=args.acceptor_cutoff,
         weights=_build_weights(getattr(args, "weight", [])),
+        required=tuple(args.require) + (() if args.soft_motifs else ("motif",)),
+        splice_strand=args.splice_strand,
+        left_flank=args.left_flank,
+        right_flank=args.right_flank,
+        protected_codons=tuple(args.protect_codon),
+        excluded_codons=tuple(args.exclude_codon),
     )
     return config
 
@@ -352,7 +371,11 @@ def _run_optimize(args: argparse.Namespace) -> int:
         protein = validate_protein(raw)
     if args.add_stop and not protein.endswith("*"):
         protein += "*"
-        start = None  # the starting sequence no longer matches the protein
+        if start is not None:
+            options = [c for c in config.table.codons("*") if c not in config.excluded_codons]
+            if not options:
+                raise _UsageError("all stop codons are excluded")
+            start += max(options, key=config.table.frequency)
 
     log = (lambda message: None) if args.quiet else _stderr
     settings = OptimizerSettings(
@@ -361,10 +384,12 @@ def _run_optimize(args: argparse.Namespace) -> int:
         anneal_steps=args.anneal,
         start_temperature=args.temp_start,
         end_temperature=args.temp_end,
-        max_seconds=None if args.seconds <= 0 else args.seconds,
+        max_seconds=None if args.seconds == 0 else args.seconds,
         seed=args.seed,
         seed_mode=args.seed_mode,
         progress=log,
+        repair_candidates=args.repair_candidates,
+        repair_max_codons=args.repair_max_codons,
     )
     if config.splice_scorer is not None and not args.quiet:
         log(f"splice backend: {config.splice_scorer.name}")
@@ -374,7 +399,7 @@ def _run_optimize(args: argparse.Namespace) -> int:
     record = Record(
         name=args.name or name,
         description=(args.name or name)
-        + f" codonopt table={config.table.name} score={result.score:.2f}",
+            + f" codonopt table={config.table.name} score={result.score:.2f} valid={str(result.valid).lower()}",
         sequence=result.sequence,
     )
     text = format_fasta(record)
@@ -385,12 +410,15 @@ def _run_optimize(args: argparse.Namespace) -> int:
         sys.stdout.write(text)
     if args.json:
         Path(args.json).write_text(
-            json.dumps(report_dict(result, config), indent=2) + "\n"
+            json.dumps(report_dict(result, config), indent=2, allow_nan=False) + "\n"
         )
         log(f"wrote {args.json}")
     if not args.quiet:
         log("")
         log(format_report(result, config))
+    if not result.valid:
+        _stderr("codonopt: unresolved mandatory constraints; output is a candidate, not a valid design")
+        return 2
     return 0