farshidk/codon-optimizer
1
1from typing import Dict, List, Tuple, Optional, Any2import os, json, math, numpy as np, pandas as pd3from collections import defaultdict, Counter4 5# ---------------- Genetic code (DNA) ----------------6AA2CODONS = {7 'A':['GCT','GCC','GCA','GCG'],8 'R':['CGT','CGC','CGA','CGG','AGA','AGG'],9 'N':['AAT','AAC'],10 'D':['GAT','GAC'],11 'C':['TGT','TGC'],12 'Q':['CAA','CAG'],13 'E':['GAA','GAG'],14 'G':['GGT','GGC','GGA','GGG'],15 'H':['CAT','CAC'],16 'I':['ATT','ATC','ATA'],17 'L':['TTA','TTG','CTT','CTC','CTA','CTG'],18 'K':['AAA','AAG'],19 'M':['ATG'],20 'F':['TTT','TTC'],21 'P':['CCT','CCC','CCA','CCG'],22 'S':['TCT','TCC','TCA','TCG','AGT','AGC'],23 'T':['ACT','ACC','ACA','ACG'],24 'W':['TGG'],25 'Y':['TAT','TAC'],26 'V':['GTT','GTC','GTA','GTG'],27 '*':['TAA','TAG','TGA']28}29DNA_Codons = {30 # 'M' - START, '_' - STOP31 "GCT": "A", "GCC": "A", "GCA": "A", "GCG": "A",32 "TGT": "C", "TGC": "C",33 "GAT": "D", "GAC": "D",34 "GAA": "E", "GAG": "E",35 "TTT": "F", "TTC": "F",36 "GGT": "G", "GGC": "G", "GGA": "G", "GGG": "G",37 "CAT": "H", "CAC": "H",38 "ATA": "I", "ATT": "I", "ATC": "I",39 "AAA": "K", "AAG": "K",40 "TTA": "L", "TTG": "L", "CTT": "L", "CTC": "L", "CTA": "L", "CTG": "L",41 "ATG": "M",42 "AAT": "N", "AAC": "N",43 "CCT": "P", "CCC": "P", "CCA": "P", "CCG": "P",44 "CAA": "Q", "CAG": "Q",45 "CGT": "R", "CGC": "R", "CGA": "R", "CGG": "R", "AGA": "R", "AGG": "R",46 "TCT": "S", "TCC": "S", "TCA": "S", "TCG": "S", "AGT": "S", "AGC": "S",47 "ACT": "T", "ACC": "T", "ACA": "T", "ACG": "T",48 "GTT": "V", "GTC": "V", "GTA": "V", "GTG": "V",49 "TGG": "W",50 "TAT": "Y", "TAC": "Y",51 "TAA": "_", "TAG": "_", "TGA": "_"52}53 54# ---------------- Helpers ----------------55def aminoacid_percentage(codons):56 """% and count of codons per amino acid for a chosen codon list."""57 amino_dict_count = defaultdict(list)58 amino_dict_per = defaultdict(list)59 for v in codons:60 amino = DNA_Codons[v]61 amino_dict_per[amino].append(v)62 for k,v in amino_dict_per.items():63 c = Counter(v)64 sub_dict = {kk:np.round(vv/len(v),2) for kk,vv in c.items()}65 amino_dict_per[k] = sub_dict66 amino_dict_count[k] = c67 return amino_dict_per,amino_dict_count68 69def gc_content(sequence=None,fasta_dir=None):70 """71 GC% overall + per codon position for each sequence provided.72 Use with: gc_content([nt_string])73 """74 if sequence is not None:75 sequences=sequence76 else:77 raise ValueError("Provide sequence=[nt_string].")78 gc_content=[[(seq.lower().count('g')+seq.lower().count('c'))/len(seq)*100] for seq in sequences]79 for index,seq in enumerate(sequences):80 seq=seq.lower()81 for i in range(3):82 position_nucleotides=seq[i::3]83 gc_count = position_nucleotides.count('g') + position_nucleotides.count('c')84 total_position_nucleotides = len(position_nucleotides)85 gc_content_percentage = (gc_count / max(1,total_position_nucleotides)) * 10086 gc_content[index].append(gc_content_percentage)87 return pd.DataFrame(gc_content,columns=['Original','Position One','Position Two','Position Three']).round(1)88 89 90def c_content(sequences: List[str]) -> pd.DataFrame:91 """92 C% overall + per codon position for each nucleotide sequence provided.93 """94 c_vals = [[(s.lower().count('c') / max(1, len(s))) * 100] for s in sequences] # <-- fixed extra ')'95 for idx, seq in enumerate(sequences):96 seq = seq.lower()97 for i in range(3):98 pos_nt = seq[i::3]99 c_count = pos_nt.count('c')100 pct = (c_count / max(1, len(pos_nt))) * 100101 c_vals[idx].append(pct)102 return pd.DataFrame(103 c_vals,104 columns=['Original C%', 'Position One C%', 'Position Two C%', 'Position Three C%']105 ).round(1)106 107def parse_kmer_list(x) -> List[str]:108 """Parse semicolon-separated kmers 'AAA;TTT;...' into a list (uppercased)."""109 if x is None:110 return []111 s = str(x).strip()112 if not s:113 return []114 return [k.strip().upper() for k in s.split(";") if k.strip()]115 116def load_summary(summary_path: str) -> pd.DataFrame:117 """Read CSV or XLSX summary to DataFrame."""118 sp = summary_path.lower()119 if sp.endswith(".xlsx") or sp.endswith(".xls"):120 return pd.read_excel(summary_path)121 return pd.read_csv(summary_path)122 123def scale_interval_codon(a_codon: int, b_codon: int,124 L_train_cds: int, L_target_cds: int) -> Tuple[int, int]:125 """Percentage-map [a,b] (1-based codons) from training length to target length."""126 a2 = 1 + int(((a_codon - 1) / max(1, L_train_cds)) * max(1, L_target_cds - 1))127 b2 = int(math.ceil((b_codon / max(1, L_train_cds)) * L_target_cds))128 a2 = max(1, min(a2, L_target_cds))129 b2 = max(1, min(b2, L_target_cds))130 return a2, b2131 132def codon_region_to_nt_span(a_codon: int, b_codon: int) -> Tuple[int, int]:133 """Convert a codon region to 1-based nucleotide span [nt_start, nt_end]."""134 return 3*(a_codon-1)+1, 3*b_codon135 136def feasible_codons_with_pattern(aa: str, pattern: str) -> List[str]:137 """138 pattern like '.C.' where '.' = free nt, returns codons for this AA matching it.139 """140 outs = []141 for c in AA2CODONS[aa]:142 ok = True143 for i, ch in enumerate(pattern):144 if ch != '.' and c[i] != ch:145 ok = False146 break147 if ok:148 outs.append(c)149 return outs150 151# ---------------- Wobble preference helpers ----------------152def wobble_bonus(base: str) -> int:153 """3rd-base preference: C=+2, G=0, A/T=-1; '.' or others -> 0."""154 if base == "C":155 return 2156 if base == "G":157 return 0158 if base in ("A", "T"):159 return -1160 return 0161 162def expected_wobble_for_pattern(aa: str, pattern: str) -> float:163 """Average wobble bonus if wobble free, or fixed wobble bonus if specified."""164 wob = pattern[2]165 if wob in "ACGT":166 return float(wobble_bonus(wob))167 feas = feasible_codons_with_pattern(aa, pattern)168 if not feas:169 return 0.0170 return sum(wobble_bonus(c[2]) for c in feas) / len(feas)171 172def placement_wobble_score(constraints: Dict[int, str], aa_seq: str) -> float:173 """Sum wobble preferences across codons touched by a seed placement."""174 total = 0.0175 for ci, patt in constraints.items():176 if 0 <= ci < len(aa_seq):177 total += expected_wobble_for_pattern(aa_seq[ci], patt)178 return total179 180# ---------------- Seeding best_kmers (enumerate placements) ----------------181def place_kmer_seed_in_region_codon(aa_seq: str,182 a_codon: int, b_codon: int,183 kmer: str,184 fixed_nt: Dict[int, str]) -> List[Dict]:185 """186 Enumerate feasible placements of a k-mer inside [a_codon,b_codon] (1-based, inclusive).187 Returns list of dicts with {"start_nt","end_nt","constraints"} where nt indices are 0-based.188 """189 nt_start, nt_end = codon_region_to_nt_span(a_codon, b_codon)190 region_start_nt0, region_end_nt0 = nt_start-1, nt_end-1191 k = len(kmer)192 placements = []193 194 for u in range(region_start_nt0, region_end_nt0 - k + 2):195 v = u + k - 1196 197 # conflict with fixed nts?198 if any((u+t) in fixed_nt and fixed_nt[u+t] != kmer[t] for t in range(k)):199 continue200 201 # build codon-level constraints202 constraints: Dict[int, str] = {}203 codon_i0 = u // 3204 codon_i1 = v // 3205 ok = True206 for ci in range(codon_i0, codon_i1 + 1):207 if ci >= len(aa_seq):208 ok = False; break209 patt = list("...")210 for ofs in range(3):211 nt_idx = ci*3 + ofs212 if u <= nt_idx <= v:213 patt[ofs] = kmer[nt_idx - u]214 patt_s = "".join(patt)215 if not feasible_codons_with_pattern(aa_seq[ci], patt_s):216 ok = False; break217 constraints[ci] = patt_s218 219 if not ok:220 continue221 222 placements.append({"start_nt": u, "end_nt": v, "constraints": constraints})223 224 return placements225 226# ---------------- Scoring while filling ----------------227def _parse_klist(s):228 if pd.isna(s) or not str(s).strip():229 return []230 return [t.strip().upper() for t in str(s).split(";") if t.strip()]231 232def _parse_kwmap(s):233 """Parse 'kmer:weight;kmer:weight;...' -> dict. If weights absent, return {}."""234 if pd.isna(s) or not str(s).strip():235 return {}236 out = {}237 for tok in str(s).split(";"):238 tok = tok.strip()239 if not tok:240 continue241 if ":" in tok:242 k, w = tok.split(":", 1)243 try:244 out[k.strip().upper()] = float(w)245 except:246 pass247 return out248 249def len_weight(K: int, is_pos: bool) -> float:250 """251 Length-aware default weights.252 Pos: +0.5*K; Neg: -1.0*K253 """254 return (0.5 * K) if is_pos else (-1.0 * K)255 256def collect_right_known_nt(ci: int,257 fixed_nt: dict,258 L_target_cds: int,259 limit_nt: int) -> str:260 """261 Collect up to 'limit_nt' contiguous known nts to the RIGHT of codon index 'ci'.262 """263 out = []264 start_nt = (ci + 1) * 3265 nt = start_nt266 while len(out) < limit_nt:267 if nt in fixed_nt:268 out.append(fixed_nt[nt])269 nt += 1270 continue271 break272 return "".join(out)273 274def enumerate_local_windows_overlapping_new(block: str,275 new_start: int, # index in block (0-based)276 new_len: int,277 K: int):278 """Yield all K-length windows inside 'block' that overlap [new_start, new_start+new_len)."""279 L = len(block)280 if K > L:281 return282 new_end = new_start + new_len # exclusive283 s_min = max(0, new_end - K)284 s_max = min(new_start, L - K)285 s_lo = max(0, new_start - (K - 1))286 s_hi = min(L - K, new_end - 1)287 s_from = min(s_min, s_lo)288 s_to = max(s_max, s_hi)289 for s in range(s_from, s_to + 1):290 if s < new_end and (s + K) > new_start:291 yield s, block[s:s+K]292 293def score_increment_multiKs(tail_nt: str,294 new_codon: str,295 best_sets: dict, # {K: set(kmer)} (we'll pass {} during fill)296 avoid_sets: dict, # {K: set(kmer)}297 pos_w: dict = None, # {(K,kmer): weight}298 neg_w: dict = None, # {(K,kmer): weight}299 wobble: bool = True,300 wobble_scale: float = 1.0,301 scoring_mode: str = "local", # "local" or "suffix"302 right_known_nt: str = "",303 hard_forbid: set = None,304 require_full_windows: bool = True):305 """306 Incremental score for appending `new_codon`.307 - Only FULL K-mer windows that are fully known are scored.308 - We build a local block = left_tail + new_codon + right_known_nt (all known nts).309 """310 HARD_KILL = -1e9311 312 # Build local block313 if scoring_mode == "local":314 left = tail_nt or ""315 mid = new_codon316 right = right_known_nt or ""317 block = left + mid + right318 new_start = len(left)319 new_len = 3320 L_block = len(block)321 else:322 block = (tail_nt + new_codon) if tail_nt else new_codon323 new_start = len(block) - len(new_codon)324 new_len = 3325 L_block = len(block)326 327 gain = 0.0328 329 Ks = sorted(best_sets.keys() | avoid_sets.keys())330 for K in Ks:331 if K > L_block:332 continue333 334 for _, km in enumerate_local_windows_overlapping_new(block, new_start, new_len, K):335 if len(km) != K:336 continue337 338 if hard_forbid and km in hard_forbid:339 return HARD_KILL340 341 # negatives (penalty)342 if km in avoid_sets.get(K, set()):343 if neg_w and (K, km) in neg_w:344 gain += -abs(neg_w[(K, km)])345 else:346 gain += len_weight(K, is_pos=False)347 348 # positives (reward) — left as option; we pass {} during fill349 if km in best_sets.get(K, set()):350 if pos_w and (K, km) in pos_w:351 gain += abs(pos_w[(K, km)])352 else:353 gain += len_weight(K, is_pos=True)354 355 if wobble:356 gain += wobble_scale * wobble_bonus(new_codon[2])357 358 return gain359 360# ---------------- Row parsing helpers ----------------361def _parse_allowed_Ks_from_row(row) -> List[int] | None:362 s = row.get("K_allowed", None)363 if isinstance(s, str) and s.strip():364 out = []365 for tok in s.split(";"):366 tok = tok.strip()367 if tok:368 try:369 out.append(int(tok))370 except:371 pass372 return out or None373 return None374 375def _row_has_best(row) -> bool:376 """True if this region has any best_kmers (positives) to seed."""377 return bool(_parse_klist(row.get("best_kmers", "")))378 379def _row_has_avoid(row) -> bool:380 """True if this region provides any avoid motifs."""381 if _parse_klist(row.get("avoid_kmers", "")):382 return True383 for K in range(2, 10):384 if _parse_klist(row.get(f"K{K}_neg", "")):385 return True386 if str(row.get(f"K{K}_neg_w", "")).strip():387 return True388 if str(row.get(f"K{K}_neg_norm", "")).strip():389 return True390 return False391 392def _build_avoid_sets_from_row(row, allowed_Ks=None):393 """394 Return NEGATIVE scoring sets + weights:395 avoid_sets: {K: set(kmer)}396 neg_w: {(K,kmer): weight}397 """398 avoid_sets, neg_w = {}, {}399 400 # From avoid_kmers (variable K)401 for km in parse_kmer_list(row.get("avoid_kmers", "")):402 K = len(km)403 if allowed_Ks is not None and K not in allowed_Ks:404 continue405 avoid_sets.setdefault(K, set()).add(km.upper())406 407 # From per-K columns408 for K in range(2, 10):409 if allowed_Ks is not None and K not in allowed_Ks:410 continue411 neg_col = f"K{K}_neg"412 if neg_col in row and pd.notna(row[neg_col]):413 for km in _parse_klist(row.get(neg_col, "")):414 if len(km) == K:415 avoid_sets.setdefault(K, set()).add(km.upper())416 417 negw_col = f"K{K}_neg_w"418 if negw_col in row and pd.notna(row[negw_col]):419 for km, w in _parse_kwmap(row.get(negw_col, "")).items():420 if len(km) == K:421 neg_w[(K, km.upper())] = float(w)422 423 return avoid_sets, neg_w424 425def _compute_target_intervals(df: pd.DataFrame,426 L_train_cds: int,427 L_target_cds: int,428 use_percent_intervals: bool = True):429 """430 Use explicit 'tgt_start','tgt_end' if present; else scale if use_percent_intervals=True;431 else clip original [start,end] to target length.432 Returns list of tuples: (tgt_a, tgt_b, row).433 """434 has_tgt = ("tgt_start" in df.columns and "tgt_end" in df.columns)435 out = []436 for _, r in df.iterrows():437 a_c, b_c = int(r["start"]), int(r["end"])438 if has_tgt and pd.notna(r.get("tgt_start")) and pd.notna(r.get("tgt_end")):439 ta, tb = int(r["tgt_start"]), int(r["tgt_end"])440 elif use_percent_intervals:441 ta, tb = scale_interval_codon(a_c, b_c, L_train_cds, L_target_cds)442 else:443 ta, tb = max(1, a_c), min(L_target_cds, b_c)444 if ta <= tb:445 out.append((ta, tb, r))446 return out447 448# ---------------- Seeding (positives) ----------------449def build_promote_queue(row,450 allowed_Ks: List[int] | None = None,451 Ks_in_order: List[int] | None = None,452 exclude_winner_K_from_extras: bool = True) -> List[str]:453 """454 Build a seeding queue. (We won't use this if only_best=True;455 kept here for completeness.)456 """457 base = [km.upper() for km in _parse_klist(row.get("best_kmers",""))]458 if allowed_Ks:459 S = set(allowed_Ks)460 base = [km for km in base if len(km) in S]461 winner = base[0] if base else None462 queue = base[:]463 K_order = (Ks_in_order or allowed_Ks or list(range(2,11))).copy()464 if exclude_winner_K_from_extras and winner:465 K_order = [K for K in K_order if K != len(winner)]466 seen_in_queue = set(queue)467 extras_by_K: dict[int, dict[str, float]] = {K: {} for K in K_order}468 for K in K_order:469 for col in (f"K{K}_pos_w", f"K{K}_pos_norm", f"K{K}_pos"):470 cell = row.get(col, "")471 if not isinstance(cell, str) or not cell.strip():472 continue473 if col.endswith(("_pos_w","_pos_norm")):474 for km, sc in _parse_kwmap(cell).items():475 km = km.upper()476 if len(km) != K or km in seen_in_queue:477 continue478 prev = extras_by_K[K].get(km, float("-inf"))479 scf = float(sc)480 if scf > prev:481 extras_by_K[K][km] = scf482 else:483 for km in _parse_klist(cell):484 km = km.upper()485 if len(km) != K or km in seen_in_queue:486 continue487 extras_by_K[K].setdefault(km, 0.0)488 for K in K_order:489 items = sorted(extras_by_K[K].items(), key=lambda kv: (-kv[1], kv[0]))490 queue.extend([km for km, _ in items])491 out, seen = [], set()492 for km in queue:493 if km not in seen:494 out.append(km); seen.add(km)495 return out496 497def seed_promote_motifs(aa_seq: str,498 a_c: int, b_c: int,499 row: pd.Series,500 fixed_nt: Dict[int, str],501 use_wobble: bool = True,502 wobble_scale: float = 1.0,503 seed_cap: int = 3,504 min_gap_codons: int = 2,505 allowed_Ks: List[int] | None = None,506 only_best: bool = False) -> tuple[List[str], List[str]]:507 """508 Greedy seeding of Promote motifs (positives).509 If only_best=True, seed strictly from 'best_kmers' and ignore Kx_pos* extras.510 """511 if only_best:512 queue = [km.upper() for km in _parse_klist(row.get("best_kmers", ""))]513 if allowed_Ks:514 S = set(allowed_Ks)515 queue = [km for km in queue if len(km) in S]516 else:517 queue = build_promote_queue(row, allowed_Ks=allowed_Ks)518 if not queue:519 return [], []520 521 placed, skipped = [], []522 placed_starts_nt: List[int] = []523 min_gap_nt = 3 * max(1, int(min_gap_codons))524 525 for km in queue:526 if seed_cap is not None and len(placed) >= seed_cap:527 break528 places = place_kmer_seed_in_region_codon(aa_seq, a_c, b_c, km, fixed_nt)529 if not places:530 skipped.append(km); continue531 532 def _ok_start_gap(p):533 s_nt = p["start_nt"]534 return all(abs(s_nt - ps) >= min_gap_nt for ps in placed_starts_nt)535 536 feas = [p for p in places if _ok_start_gap(p)]537 if not feas:538 skipped.append(km); continue539 540 best_p, best_s = None, float("-inf")541 for p in feas:542 s = 0.0543 if use_wobble:544 s += wobble_scale * placement_wobble_score(p["constraints"], aa_seq)545 s += -0.001 * p["start_nt"]546 if s > best_s:547 best_s, best_p = s, p548 549 if best_p is None:550 skipped.append(km); continue551 552 u, v = best_p["start_nt"], best_p["end_nt"]553 if any((t in fixed_nt and fixed_nt[t] != km[t - u]) for t in range(u, v + 1)):554 skipped.append(km); continue555 556 for t in range(u, v + 1):557 fixed_nt[t] = km[t - u]558 placed_starts_nt.append(u)559 placed.append(km)560 561 return placed, skipped562 563# ---------------- Main: seed positives; avoid-only scoring on fill ----------------564def optimization(summary_path: str,aa_seq: str,565 use_wobble: bool = True,566 wobble_scale: float = 1.0,567 use_percent_intervals: bool = True,568 seed_cap: int | None = 3,569 min_gap_codons: int = 2,570 beam_size: int = 5):571 """572 - Seed ONLY rows that have `best_kmers` (positives).573 - Global left->right beam fill.574 - Inside mapped regions that have avoiders, apply NEGATIVE scoring (avoid sets).575 No positive rewards during fill. Outside all regions: wobble only.576 Returns: (designed_nt, aa_stats_stub, gc_percent_df, log_info)577 """578 df = load_summary(summary_path)579 L_target_cds = len(aa_seq)580 L_train_cds = int(df["end"].max())581 582 # 1) map every summary row to a target interval583 intervals = _compute_target_intervals(df, L_train_cds, L_target_cds, use_percent_intervals)584 585 # 2) seed positives (best_kmers) only586 fixed_nt: Dict[int, str] = {}587 seed_log: Dict[Tuple[int,int], Dict[str, Any]] = {}588 for (ta, tb, row) in intervals:589 if not _row_has_best(row):590 continue591 allowed_Ks = _parse_allowed_Ks_from_row(row)592 placed, skipped = seed_promote_motifs(593 aa_seq=aa_seq, a_c=ta, b_c=tb, row=row, fixed_nt=fixed_nt,594 use_wobble=use_wobble, wobble_scale=wobble_scale,595 seed_cap=seed_cap, min_gap_codons=min_gap_codons,596 allowed_Ks=allowed_Ks, only_best=True597 )598 seed_log[(ta, tb)] = {"placed": placed, "skipped": skipped}599 600 # 3) build per-codon avoid-only config601 per_ci_cfg: List[Optional[Tuple[Dict[int,set], Dict[Tuple[int,str],float], int]]] = [None] * L_target_cds602 Kmax_all = 2603 for (ta, tb, row) in intervals:604 if not _row_has_avoid(row):605 continue606 allowed_Ks = _parse_allowed_Ks_from_row(row)607 avoid_sets, neg_w = _build_avoid_sets_from_row(row, allowed_Ks)608 if not avoid_sets:609 continue610 Kmax_here = max(avoid_sets.keys())611 Kmax_all = max(Kmax_all, Kmax_here)612 for ci in range(ta - 1, tb):613 if per_ci_cfg[ci] is None:614 per_ci_cfg[ci] = ({K:set(v) for K,v in avoid_sets.items()}, dict(neg_w), Kmax_here)615 else:616 old_avoid, old_negw, old_kmax = per_ci_cfg[ci]617 for K, ss in avoid_sets.items():618 old_avoid.setdefault(K, set()).update(ss)619 old_negw.update(neg_w)620 per_ci_cfg[ci] = (old_avoid, old_negw, max(old_kmax, Kmax_here))621 622 # 4) one global left->right beam fill (avoid-only scoring, wobble everywhere)623 k_len_ref = Kmax_all624 tail_nt = ""625 beam: List[Tuple[float, str, List[str], Dict[int,str]]] = [(0.0, tail_nt, [], dict())]626 627 for ci in range(L_target_cds):628 aa = aa_seq[ci]629 # enforce seeded nts630 patt = list("...")631 for ofs in range(3):632 nt_idx = ci * 3 + ofs633 if nt_idx in fixed_nt:634 patt[ofs] = fixed_nt[nt_idx]635 patt_s = "".join(patt)636 cand_codons = feasible_codons_with_pattern(aa, patt_s) or AA2CODONS[aa]637 638 cfg = per_ci_cfg[ci]639 if cfg is None:640 avoid_sets, neg_w, Kmax_here = {}, {}, 2641 else:642 avoid_sets, neg_w, Kmax_here = cfg643 644 new_beam = []645 for score, tail, local_codons, local_fix in beam:646 right_known = collect_right_known_nt(ci, fixed_nt, L_target_cds, limit_nt=Kmax_here-1)647 for c in cand_codons:648 gain = score_increment_multiKs(649 tail_nt=tail,650 new_codon=c,651 best_sets={}, # << NO positive rewards during fill652 avoid_sets=avoid_sets,653 pos_w=None,654 neg_w=neg_w if neg_w else None,655 wobble=use_wobble,656 wobble_scale=wobble_scale,657 scoring_mode="local",658 right_known_nt=right_known,659 hard_forbid=None660 )661 if gain <= -1e5:662 continue663 tail2 = (tail + c)[-(k_len_ref - 1):] if k_len_ref > 1 else ""664 local2 = local_codons + [c]665 fix2 = dict(local_fix)666 for ofs, ch in enumerate(c):667 fix2[ci * 3 + ofs] = ch668 new_beam.append((score + gain, tail2, local2, fix2))669 670 if not new_beam:671 # safety fallback672 c = cand_codons[0]673 tail2 = (beam[0][1] + c)[-(k_len_ref - 1):] if k_len_ref > 1 else ""674 local2 = beam[0][2] + [c]675 fix2 = dict(beam[0][3])676 for ofs, ch in enumerate(c):677 fix2[ci * 3 + ofs] = ch678 new_beam = [(beam[0][0], tail2, local2, fix2)]679 680 new_beam.sort(key=lambda t: t[0], reverse=True)681 beam = new_beam[:beam_size]682 683 best_score, _, best_local, best_fix = max(beam, key=lambda t: t[0])684 685 # 5) finalize686 chosen_codons = best_local687 fixed_nt.update(best_fix)688 designed_nt = "".join(chosen_codons)689 gc_percent = gc_content([designed_nt])690 c_percent = c_content([designed_nt])691 # Build amino-acid → {codon: fraction} dict for UI692 aa_percent_dict, _aa_counts = aminoacid_percentage(chosen_codons)693 694 log_info = [{695 "mode": "global_fill_avoid_only",696 "beam_best_score": best_score,697 "beam_kept": beam_size,698 "seed_summary": seed_log699 }]700 701 # aa_percent placeholder kept for API compatibility702 return designed_nt, aa_percent_dict, gc_percent, c_percent, log_info703 704 705# ---------------- Example usage ----------------706if __name__ == "__main__":707 # Example:708 # summary_path = "region_sweep_summary.csv"709 # aa_seq = "M" + "A"*514 # 515 aa example710 # nt_seq, aa_stats, gc_df, log = optimization(711 # summary_path=summary_path,712 # aa_seq=aa_seq,713 # use_wobble=True,714 # wobble_scale=1.0,715 # use_percent_intervals=False, # True to scale; False if you have tgt_start/tgt_end716 # seed_cap=3,717 # min_gap_codons=2,718 # beam_size=5719 # )720 # print(nt_seq)721 # print(gc_df)722 pass723