Bio Codon Usage

Version Compatibility

Reference examples tested with: BioPython 1.83+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Codon Usage

Analyze codon usage patterns, score adaptation to a host, and optimize coding sequences for expression.

"Analyze codon usage" -> Count codons in a coding sequence, compute frequencies and bias metrics.

• Python: Counter on in-frame triplets + RSCU/Nc helpers (BioPython + standard library)

"Score this gene against a host" -> Compute the Codon Adaptation Index from a reference set of highly expressed genes.

• Python: CodonAdaptationIndex(reference_seqs).calculate(query) (BioPython)

"Optimize codons for expression" -> Replace each codon with the host's single most-preferred synonymous codon.

• Python: CodonAdaptationIndex(reference_seqs).optimize(seq) (BioPython)

The governing principle

CAI is meaningless without an expression-biased reference. The relative-adaptiveness weights (w) must be built from the highly expressed genes of the TARGET organism (ribosomal proteins, elongation factors). A CAI computed against a whole-genome average, or against the wrong organism, is a number with no biological meaning. There is no bundled reference index in modern Biopython, so the reference set is always the caller's responsibility.

Two silent traps dominate this skill:

• Out-of-frame input is silently corrupted. calculate blindly steps range(0, len, 3) from position 0; it never checks the reading frame. A frame-shifted CDS returns a plausible CAI computed from garbage codons. Frame correctness is the caller's job (length divisible by 3, starts at the first base of codon 1).
• Naive max-CAI optimization can REDUCE expression. optimize() is the textbook max-CAI output (single most-frequent codon per amino acid). It is blind to the translation ramp, 5' mRNA structure, GC extremes, cryptic regulatory elements, and codon-pair bias. It is a starting point to screen, never a final design. Failure is silent: correct protein, high CAI, poor expression.

CRITICAL API migration (Biopython 1.82)

Bio.SeqUtils.CodonUsage and Bio.SeqUtils.CodonUsageIndices were removed in Biopython 1.82. Any code calling generate_index(), set_cai_index(), cai_for_gene(), print_index(), or importing SharpEcoliIndex raises ImportError on any modern install. The replacement is a redesigned class imported directly from Bio.SeqUtils:

from Bio.SeqUtils import CodonAdaptationIndex  # NOT Bio.SeqUtils.CodonUsage (removed)

Required Imports

from Bio import SeqIO
from Bio.Seq import Seq
from Bio.SeqUtils import CodonAdaptationIndex, GC123
from Bio.Data.CodonTable import standard_dna_table
from Bio.Data import CodonTable
from collections import Counter

Codon Adaptation Index (CAI)

Goal: Measure how closely a gene's codon usage matches the highly expressed genes of a host organism.

Approach: Build a CodonAdaptationIndex from a reference set of highly expressed CDS (the constructor computes per-codon relative adaptiveness w), then score query sequences with calculate (0-1, higher = better adapted).

from Bio import SeqIO
from Bio.Seq import Seq
from Bio.SeqUtils import CodonAdaptationIndex
from Bio.Data.CodonTable import standard_dna_table

# Reference = highly expressed genes of the TARGET host (ribosomal proteins, EFs).
# Parse a FASTA yourself; there is no bundled index. Pass str/Seq/SeqRecord.
reference_seqs = list(SeqIO.parse('highly_expressed_genes.fasta', 'fasta'))
cai = CodonAdaptationIndex(reference_seqs, table=standard_dna_table)

query = Seq('ATGAAACGTGCTGAAGCTAAATAA')
score = cai.calculate(query)   # 0-1; the query MUST be in-frame (see governing principle)
print(f'CAI: {score:.3f}')

CodonAdaptationIndex is a dict subclass -- the codon->w mapping is the object itself. Inspect or override weights directly:

print(cai['GCT'])         # relative adaptiveness of Ala codon GCT
cai.update({'GCT': 0.9})  # override a weight (replaces the old set_cai_index)

Verified behavior (Biopython >=1.82):

• ATG (Met) and TGG (Trp) are excluded from CAI -- single-codon families, w is always 1.
• Stop codons are excluded.
• Unobserved codons get w = 0.5 (Sharp & Li), softly down-weighted; no division-by-zero.
• Case-insensitive -- both the constructor and calculate uppercase internally.
• An illegal codon (non-ACGT) in a reference raises ValueError; an illegal or trailing-partial codon in a query raises TypeError. Out-of-frame input does NOT raise -- it is silently mis-scored.

RSCU = w is built from these ratios

CAI weights come from RSCU. w_ij = RSCU_ij / RSCU_jmax = (codon count) / (count of the most-used synonymous codon in that family); CAI = exp((1/L) * sum ln w). RSCU itself = observed count / expected-if-uniform within a synonymous family (=1 no bias, >1 over-used, <1 under-used). RSCU normalizes away amino-acid composition, which is why w is built from RSCU ratios rather than raw frequencies.

Goal: Quantify synonymous codon bias to detect translational selection or mutational pressure.

Approach: Group codons by amino acid via the codon table, then divide each codon's observed count by the family mean.

from Bio.Data import CodonTable
from collections import Counter

def count_codons(seq):
    s = str(seq).upper()
    return Counter(s[i:i+3] for i in range(0, len(s) - 2, 3))

def calculate_rscu(seq, table_id=1):
    '''RSCU per codon: observed / expected-if-uniform within its synonymous family'''
    table = CodonTable.unambiguous_dna_by_id[table_id]
    counts = count_codons(seq)
    back_table = {}
    for codon, aa in table.forward_table.items():
        back_table.setdefault(aa, []).append(codon)
    rscu = {}
    for aa, codons in back_table.items():
        total = sum(counts.get(c, 0) for c in codons)
        expected = total / len(codons) if codons else 0
        for codon in codons:
            rscu[codon] = counts.get(codon, 0) / expected if expected > 0 else 0
    return rscu

Codon optimization for expression

Goal: Generate a host-adapted CDS that preserves the protein.

Approach: optimize() swaps each amino acid for the host's single most-preferred synonymous codon (max-CAI). Always confirm the protein is unchanged, then screen the design against the tradeoffs below.

opt = cai.optimize(query, seq_type='DNA', strict=True)
assert opt.translate() == query.translate()   # protein must be identical

optimize(sequence, seq_type='DNA'|'RNA'|'protein', strict=True): strict=True raises ValueError on a tie (two equally-preferred codons, e.g. 'TTT and TTC are equally preferred.'); strict=False warns and picks one.

Why naive max-CAI can HURT expression

optimize() is blind to everything except single-codon frequency. Screen the output for:

• The translation ramp (Tuller et al. 2010): a conserved profile of slow (rare) codons over the first ~30-50 codons spaces ribosomes; flattening it can lower yield and increase misfolding.
• 5' mRNA secondary structure: strong folding near the start codon impedes initiation. Minimize 5' free energy, sometimes against CAI.
• GC extremes: swaps that push GC very high create stable hairpins; very low GC destabilizes.
• Cryptic elements created by swaps: splice sites, internal Shine-Dalgarno/RBS, polyadenylation signals, restriction sites, AU-rich destabilizing elements -- silent in protein, corrupting in expression.
• Codon-pair bias: decoding efficiency depends on adjacent codon pairs; CAI scores single codons only.