"""
CodonAdaptPy Constrained Sequence Optimization
This module designs multiple synonymous coding sequences while preserving the
translated amino-acid sequence. Candidate ranking can increase or decrease
reference CAI, approach a target GC fraction, control CpG/UpA, avoid or require
motifs and restriction sites, limit homopolymers, and account for codon-pair
scores. Results expose each trade-off rather than hiding it in one sequence.
Classes:
- OptimizationConfig: Objectives, constraints, and reproducibility options.
- OptimizationCandidate: One designed sequence and its diagnostics.
- OptimizationResult: Ranked candidates and source metadata.
- CodonOptimizer: Reproducible stochastic constrained optimizer.
:Created: July 20, 2026
:Updated: July 20, 2026
:Author: Naveen Duhan
:Version: 1.0.2
"""
from __future__ import annotations
import random
import re
from dataclasses import asdict, dataclass, field
from typing import Any, Literal
from .exceptions import AnalysisError
from .genetic_code import GeneticCode
from .metrics.adaptation import AdaptationMetrics
from .metrics.composition import CompositionMetrics
from .metrics.pairs import PairMetrics
[docs]
@dataclass(slots=True)
class OptimizationConfig:
"""Configure candidate search objectives and hard/soft constraints."""
direction: Literal["optimize", "deoptimize", "match"] = "optimize"
candidates: int = 5
iterations: int = 5000
seed: int = 1
target_gc: float | None = None
gc_tolerance: float = 0.05
avoid_motifs: tuple[str, ...] = ()
preserve_motifs: tuple[str, ...] = ()
remove_restriction_sites: tuple[str, ...] = ()
add_restriction_sites: tuple[str, ...] = ()
max_homopolymer: int = 6
target_cpg: float | None = None
target_upa: float | None = None
cai_weight: float = 1.0
gc_weight: float = 1.0
motif_penalty: float = 10.0
pair_weight: float = 0.0
@dataclass(slots=True)
class OptimizationCandidate:
"""Store one synonymous candidate and transparent objective diagnostics."""
rank: int
sequence: str
score: float
cai: float
gc: float
cpg_representation: float
upa_representation: float
constraint_violations: list[str] = field(default_factory=list)
@dataclass(slots=True)
class OptimizationResult:
"""Collect source sequence, configuration, and ranked candidate designs."""
source_sequence: str
amino_acid_sequence: str
candidates: list[OptimizationCandidate]
config: OptimizationConfig
genetic_code: int
def to_dict(self) -> dict[str, Any]:
"""Return nested optimization output in JSON-compatible form."""
return asdict(self)
[docs]
class CodonOptimizer:
"""Generate and rank synonymous sequences using seeded stochastic search."""
def __init__(
self, reference_counts: dict[str, float], *, genetic_code: int = 1, pair_scores: dict[str, float] | None = None
) -> None:
"""Initialize optimizer with a target reference profile."""
self.code = GeneticCode.from_ncbi(genetic_code)
self.adaptation = AdaptationMetrics(self.code)
self.composition = CompositionMetrics()
self.pairs = PairMetrics(self.code)
self.weights = self.adaptation.weights_from_counts(reference_counts)
self.reference_counts = reference_counts
self.pair_scores = pair_scores or {}
[docs]
def optimize(self, sequence: str, config: OptimizationConfig | None = None) -> OptimizationResult:
"""Generate several ranked synonymous candidate coding sequences."""
config = config or OptimizationConfig()
normalized = "".join(sequence.upper().replace("U", "T").split())
if len(normalized) % 3:
raise AnalysisError("Optimization requires a sequence divisible by three.")
codons = [normalized[index : index + 3] for index in range(0, len(normalized), 3)]
amino_acids = [self.code.translate(codon) for codon in codons]
if any(amino_acid is None for amino_acid in amino_acids):
raise AnalysisError("Optimization does not accept ambiguous or invalid codons.")
terminal_stop = bool(amino_acids and amino_acids[-1] == "*")
if "*" in amino_acids[:-1]:
raise AnalysisError("Optimization does not accept internal stop codons.")
design_amino_acids = [aa for aa in amino_acids if aa is not None and aa != "*"]
stop = codons[-1] if terminal_stop else ""
rng = random.Random(config.seed)
designs: dict[str, tuple[float, dict[str, float], list[str]]] = {}
initial = "".join(codons[:-1] if terminal_stop else codons) + stop
for _ in range(max(config.iterations, config.candidates)):
candidate_codons = [self._choose_codon(aa, config.direction, rng) for aa in design_amino_acids]
candidate = "".join(candidate_codons) + stop
score, diagnostics, violations = self._score(candidate, candidate_codons, config)
previous = designs.get(candidate)
if previous is None or score > previous[0]:
designs[candidate] = (score, diagnostics, violations)
if initial not in designs:
initial_codons = codons[:-1] if terminal_stop else codons
designs[initial] = self._score(initial, initial_codons, config)
ordered = sorted(designs.items(), key=lambda item: item[1][0], reverse=True)[: config.candidates]
candidates = [
OptimizationCandidate(
rank,
candidate,
score,
diagnostics["cai"],
diagnostics["gc"],
diagnostics["cpg"],
diagnostics["upa"],
violations,
)
for rank, (candidate, (score, diagnostics, violations)) in enumerate(ordered, start=1)
]
return OptimizationResult(normalized, "".join(design_amino_acids), candidates, config, self.code.table_id)
def _choose_codon(self, amino_acid: str, direction: str, rng: random.Random) -> str:
"""Sample a synonymous codon according to the requested direction."""
family = self.code.synonymous_codons[amino_acid]
weights = [self.weights[codon] for codon in family]
if direction == "deoptimize":
weights = [1 / max(weight, 1e-9) for weight in weights]
elif direction == "match":
weights = [max(float(self.reference_counts.get(codon, 0.0)), 1e-9) for codon in family]
return rng.choices(family, weights=weights, k=1)[0]
def _score(
self, sequence: str, codons: list[str], config: OptimizationConfig
) -> tuple[float, dict[str, float], list[str]]:
"""Evaluate objectives and report every constraint violation."""
cai = self.adaptation.cai(codons, self.weights)
composition = self.composition.gc_metrics(sequence)
cpg = self.composition.representation(sequence, "CG")
upa = self.composition.representation(sequence, "TA")
cai_term = -cai if config.direction == "deoptimize" else cai
score = config.cai_weight * cai_term
if config.target_gc is not None:
score -= config.gc_weight * abs(composition["gc"] - config.target_gc)
if config.target_cpg is not None:
score -= abs(cpg - config.target_cpg)
if config.target_upa is not None:
score -= abs(upa - config.target_upa)
if self.pair_scores and config.pair_weight:
score += config.pair_weight * self.pairs.bias(codons, self.pair_scores)
violations: list[str] = []
forbidden = tuple(
item.upper().replace("U", "T") for item in config.avoid_motifs + config.remove_restriction_sites
)
required = tuple(
item.upper().replace("U", "T") for item in config.preserve_motifs + config.add_restriction_sites
)
violations.extend(f"forbidden motif present: {motif}" for motif in forbidden if motif and motif in sequence)
violations.extend(f"required motif absent: {motif}" for motif in required if motif and motif not in sequence)
if config.target_gc is not None and abs(composition["gc"] - config.target_gc) > config.gc_tolerance:
violations.append("GC content outside tolerance")
if config.max_homopolymer > 0 and re.search(rf"([ACGT])\1{{{config.max_homopolymer},}}", sequence):
violations.append(f"homopolymer longer than {config.max_homopolymer}")
score -= config.motif_penalty * len(violations)
return score, {"cai": cai, "gc": composition["gc"], "cpg": cpg, "upa": upa}, violations