farshidk/codon-optimizer
1
1# test.py2# Simple C-wobbling optimizer: prefers wobble C (C > G > A > T) and3# avoids long runs of 'CCC' by locally switching to 'CCG' when needed.4 5from collections import defaultdict, Counter6from typing import Dict, List, Tuple7import numpy as np8import pandas as pd9 10# ---------------- Genetic code (DNA) ----------------11DNA_Codons: Dict[str, str] = {12 # 'M' - START, '_' - STOP13 "GCT": "A", "GCC": "A", "GCA": "A", "GCG": "A",14 "TGT": "C", "TGC": "C",15 "GAT": "D", "GAC": "D",16 "GAA": "E", "GAG": "E",17 "TTT": "F", "TTC": "F",18 "GGT": "G", "GGC": "G", "GGA": "G", "GGG": "G",19 "CAT": "H", "CAC": "H",20 "ATA": "I", "ATT": "I", "ATC": "I",21 "AAA": "K", "AAG": "K",22 "TTA": "L", "TTG": "L", "CTT": "L", "CTC": "L", "CTA": "L", "CTG": "L",23 "ATG": "M",24 "AAT": "N", "AAC": "N",25 "CCT": "P", "CCC": "P", "CCA": "P", "CCG": "P",26 "CAA": "Q", "CAG": "Q",27 "CGT": "R", "CGC": "R", "CGA": "R", "CGG": "R", "AGA": "R", "AGG": "R",28 "TCT": "S", "TCC": "S", "TCA": "S", "TCG": "S", "AGT": "S", "AGC": "S",29 "ACT": "T", "ACC": "T", "ACA": "T", "ACG": "T",30 "GTT": "V", "GTC": "V", "GTA": "V", "GTG": "V",31 "TGG": "W",32 "TAT": "Y", "TAC": "Y",33 "TAA": "_", "TAG": "_", "TGA": "_"34}35 36# Codons ordered by wobble priority C > G > A > T, uppercase37aa_to_codons: Dict[str, List[str]] = {38 'A': ['GCC', 'GCG', 'GCA', 'GCT'],39 'R': ['CGC', 'CGG', 'CGA', 'CGT', 'AGG', 'AGA'],40 'N': ['AAC', 'AAT'],41 'D': ['GAC', 'GAT'],42 'C': ['TGC', 'TGT'],43 'Q': ['CAG', 'CAA'],44 'E': ['GAG', 'GAA'],45 'G': ['GGC', 'GGG', 'GGA', 'GGT'],46 'H': ['CAC', 'CAT'],47 'I': ['ATC', 'ATA', 'ATT'],48 'L': ['CTC', 'CTG', 'CTA', 'CTT', 'TTG', 'TTA'],49 'K': ['AAG', 'AAA'],50 'M': ['ATG'],51 'F': ['TTC', 'TTT'],52 'P': ['CCC', 'CCG', 'CCA', 'CCT'],53 'S': ['TCC', 'TCG', 'TCA', 'TCT', 'AGC', 'AGT'],54 'T': ['ACC', 'ACG', 'ACA', 'ACT'],55 'W': ['TGG'],56 'Y': ['TAC', 'TAT'],57 'V': ['GTC', 'GTG', 'GTA', 'GTT'],58 '*': ['TGA', 'TAG', 'TAA']59}60 61# ---------------- Helpers ----------------62def aminoacid_percentage(codons: List[str]):63 """64 Build {AA: {codon: fraction}} and counts for the chosen codon list.65 Fractions are rounded to 2 decimals.66 """67 amino_dict_count = defaultdict(list)68 amino_dict_per = defaultdict(list)69 for v in codons:70 amino = DNA_Codons[v]71 amino_dict_per[amino].append(v)72 for k, v in amino_dict_per.items():73 c = Counter(v)74 sub_dict = {kk: np.round(vv / len(v), 2) for kk, vv in c.items()}75 amino_dict_per[k] = sub_dict76 amino_dict_count[k] = c77 return amino_dict_per, amino_dict_count78 79def translate_seq_to_aa(seq,remove_stop=True):80 if isinstance(seq, list):81 seq[0]=seq[0].upper()82 else:83 seq=[seq.upper()]84 85 if 'U' in seq:86 cds_ref=RNA_Codons87 else:88 cds_ref=DNA_Codons89 seq_list=seq_to_cds(seq,sep=' ',sos_eos=False,remove_stop=remove_stop)90 seq_list=seq_list[0].split()91 aa_list=""92 for i in seq_list:93 aa_list+=cds_ref[i]94 return aa_list95 96def gc_content(sequences: List[str]) -> pd.DataFrame:97 """98 GC% overall + per codon position for each nucleotide sequence provided.99 """100 gc_vals = [[(s.lower().count('g') + s.lower().count('c')) / len(s) * 100] for s in sequences]101 for idx, seq in enumerate(sequences):102 seq = seq.lower()103 for i in range(3):104 pos_nt = seq[i::3]105 gc_count = pos_nt.count('g') + pos_nt.count('c')106 pct = (gc_count / max(1, len(pos_nt))) * 100107 gc_vals[idx].append(pct)108 return pd.DataFrame(gc_vals,109 columns=['Original', 'Position One', 'Position Two', 'Position Three']110 ).round(1)111 112def c_content(sequences: List[str]) -> pd.DataFrame:113 """114 C% overall + per codon position for each nucleotide sequence provided.115 """116 c_vals = [[(s.lower().count('c') / max(1, len(s))) * 100] for s in sequences] # <-- fixed extra ')'117 for idx, seq in enumerate(sequences):118 seq = seq.lower()119 for i in range(3):120 pos_nt = seq[i::3]121 c_count = pos_nt.count('c')122 pct = (c_count / max(1, len(pos_nt))) * 100123 c_vals[idx].append(pct)124 return pd.DataFrame(125 c_vals,126 columns=['Original C%', 'Position One C%', 'Position Two C%', 'Position Three C%']127 ).round(1)128 129 130def seq_to_cds(seqs,sep=None,sos_eos=False,remove_stop=False):131 """seperate nucleotide sequence to codons (triple form).132 133 Parameters134 ----------135 seqs: the list of sequences for converting to codon format.136 sep: the separator you choose to separate codons based on it.137 if you don't specify any separator(default) each codons138 will separate by whitespace.139 sos_eos: if True, start of sequence character SOS and end of140 sequence character EOS will be added to the start 141 and the end of a sequence.142 143 Example:144 import structure as s145 cds=s.fasta_to_list("C:\\Users\\farsh\\Downloads\\cds.txt")146 codons=s.seq_to_cds(cds)147 """148 seqs_list=[]149 for seq in seqs:150 str_cds=[]151 for pos in range(0,len(seq)-len(seq)%3,3):152 str_cds.append(seq[pos:pos+3])153 if remove_stop:154 if str_cds[-1] in ("UAA","UAG","UGA","TAA","TAG","TGA"): #remove stop codons155 del str_cds[-1]156 if sos_eos:157 str_cds.insert(0,'[SOS]')158 str_cds.append('[EOS]')159 sep=' '160 seqs_list.append(sep.join(str_cds))161 return seqs_list162 163# ---------------- C-wobbling policy ----------------164def c_wobble(aa_seq: str) -> List[str]:165 """166 Choose the first codon in our wobble-ordered list (C>G>A>T),167 with a small local rule to avoid long 'CCC' runs.168 """169 aa_seq = aa_seq.upper()170 chosen: List[str] = []171 for i, amino in enumerate(aa_seq):172 cd = aa_to_codons[amino][0] # prefers wobble C by construction173 if i >= 1:174 # if current is 'CCC' and previous ends with 'CC', break the run:175 if cd == 'CCC' and chosen[i-1] in ('ACC', 'GCC', 'TCC'):176 cd = aa_to_codons[amino][1] # -> 'CCG'177 # still two 'CCC' in a row? flip previous to 'CCG'178 if cd == 'CCC' and chosen[i-1] == 'CCC':179 chosen[i-1] = 'CCG'180 if i > 1 and cd.startswith("C") and chosen[i-1] == "CCC" and chosen[i-2].endswith("C"):181 chosen[i-1] = "CCG"182 chosen.append(cd)183 return chosen184 185# ---------------- prevent restriction enzyme site ----------------186def codons_interval(seq,res_site):187 seq = seq.upper()188 res_site = res_site.upper()189 positions = []190 codonsIntervals = []191 start = 0192 while True:193 pos = seq.find(res_site, start)194 if pos == -1:195 break196 positions.append(pos)197 start = pos + 1198 for p in positions:199 codons = [p//3,(p+len(res_site)-1)//3]200 codonsInterval = list(range(codons[0], codons[-1] + 1))201 codonsIntervals.append(codonsInterval)202 return codonsIntervals,positions203 204def extract_specific_codons(aa,wobble):205 "extract specific codons that ends with a specific wobble"206 wobble = wobble.upper()207 specific_codons = []208 for c in aa_to_codons[aa]:209 if c[-1] == wobble:210 specific_codons.append(c)211 return specific_codons212 213# ---------------- Public API (compatible with your app) ----------------214def do_c_wobble(summary_path: str | None,aa_seq: str,**kwargs):215 """216 Returns: (designed_nt, aa_percent_dict, gc_percent_df, log_info)217 'summary_path' kept for signature compatibility; not used here.218 """219 codons = c_wobble(aa_seq)220 designed_nt = ''.join(codons)221 aa_percent_dict, _counts = aminoacid_percentage(codons)222 gc_df = gc_content([designed_nt])223 c_df = c_content([designed_nt])224 log_info = [{"mode": "c_wobbling_simple"}]225 return designed_nt, aa_percent_dict, gc_df, c_df, log_info226 227def change_codon(seqs, res_site):228 """229 Change codons to remove restriction sites while preserving amino acids.230 231 Returns232 -------233 new_seq : str234 Modified nucleotide sequence235 modified_positions : list[int]236 Nucleotide indices that were changed237 bold_intervals : list[tuple[int, int]]238 Nucleotide intervals of detected restriction sites, stored as239 Python-style half-open intervals: (start, end)240 meaning start <= i < end241 """242 243 seqs = seqs.upper()244 res_site = res_site.upper()245 246 modified_positions = []247 bold_intervals = []248 249 codons_intervals, nucleotides = codons_interval(seqs, res_site)250 seqs_list = seq_to_cds([seqs])[0].split()251 252 for interval, nuc_start in zip(codons_intervals, nucleotides):253 start_codon = interval[0]254 end_codon = interval[-1] + 1255 region_codons = seqs_list[start_codon:end_codon]256 257 changed = False258 259 for priority in ['C', 'G', 'A', 'T']:260 for i, old_codon in enumerate(region_codons):261 aa = DNA_Codons[old_codon]262 candidate_codons = extract_specific_codons(aa, priority)263 264 for new_codon in candidate_codons:265 if new_codon != old_codon:266 temp_codons = region_codons.copy()267 temp_codons[i] = new_codon268 269 # check if restriction site is removed from this local region270 if res_site not in ''.join(temp_codons):271 global_codon_index = start_codon + i272 old_global_codon = seqs_list[global_codon_index]273 274 # record which nucleotide positions changed275 for j in range(3):276 if old_global_codon[j] != new_codon[j]:277 modified_positions.append(global_codon_index * 3 + j)278 279 # apply change to the full sequence280 seqs_list[global_codon_index] = new_codon281 282 # store bold interval in nucleotide coordinates283 # half-open interval: [nuc_start, nuc_start + len(res_site))284 bold_intervals.append((nuc_start, nuc_start + len(res_site)))285 286 changed = True287 break288 289 if changed:290 break291 if changed:292 break293 294 return ''.join(seqs_list), modified_positions, bold_intervals