# Bio Core Computational Genetics > Translate DNA per-frame, score codon usage bias (RSCU/CAI), simulate restriction digests/ORFs, and test three-point-cross mapping and Hardy-Weinberg equilibrium. Use for CAI, virtual digests, crossover mapping, or HWE tests. - Skill: `pavel-kravchenko/bio-core-computational-genetics` (Agent Skill) - Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-core-computational-genetics` - Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-core-computational-genetics/raw - Safety review: pending - Works with: Claude Code, Claude.ai, OpenAI Codex - Category: Coding & Dev Tools - Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko) - Updated: 2026-09-17 - Page: https://skillmd.com/skills/pavel-kravchenko/bio-core-computational-genetics --- # Computational Genetics ## When to Use - Translating DNA in a given reading frame or finding ORFs across all 6 frames - Computing codon usage bias (RSCU) or Codon Adaptation Index (CAI) for codon optimization / heterologous expression - Simulating a restriction enzyme digest (fragment sizes, virtual gel) for cloning design - Mapping gene order and distance from two-point or three-point cross data - Testing observed genotype counts against Hardy-Weinberg equilibrium, or computing Ts/Tv ratio and CpG O/E from aligned sequences ## Version Compatibility Python >= 3.10, NumPy >= 1.24, SciPy >= 1.11 (for exact chi-squared p-values), Matplotlib >= 3.8. No external bioinformatics package is required for the core logic below; Biopython (`Bio.Seq`, `Bio.Data.CodonTable`) is a drop-in replacement for the hand-rolled genetic code if already in your environment. ## Prerequisites - `pip install numpy scipy matplotlib` - Comfortable with basic Mendelian genetics (dominant/recessive, linkage, crossover) and DNA/codon notation - Related skills: `bio-sequence-manipulation-codon-usage`, `bio-sequence-manipulation-transcription-translation`, `bio-restriction-analysis-restriction-mapping` ## Genetic Code, Translation, and Codon Usage **Goal:** Build the standard genetic code table, translate DNA in a chosen frame, and score a coding sequence's codon usage bias (RSCU, CAI) against a reference. **Approach:** Generate the 64-codon table programmatically (avoids typos), translate stopping at the first in-frame stop codon, then compute RSCU per synonymous family and CAI as the geometric mean of relative codon adaptiveness. ```python from collections import defaultdict, Counter import math bases = 'TCAG' amino_acids = 'FFLLSSSSYY**CC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG' codon_table = {} for i, first in enumerate(bases): for j, second in enumerate(bases): for k, third in enumerate(bases): codon_table[first + second + third] = amino_acids[i * 16 + j * 4 + k] stop_codons = {c for c, a in codon_table.items() if a == '*'} def translate(dna_sequence, codon_table, stop_at_stop=True): """Translate a DNA sequence codon-by-codon, stopping at the first stop codon.""" seq = dna_sequence.upper() protein = [] for i in range(0, len(seq) - 2, 3): codon = seq[i:i + 3] aa = codon_table.get(codon, 'X') if aa == '*' and stop_at_stop: break protein.append(aa) return ''.join(protein) def translate_frame(dna_sequence, codon_table, frame=0): """Translate starting from a given reading frame (0, 1, or 2).""" return translate(dna_sequence[frame:], codon_table) def count_codons(coding_sequence): """Count in-frame codons in a coding sequence (must start at position 0).""" seq = coding_sequence.upper() counts = Counter() for i in range(0, len(seq) - 2, 3): codon = seq[i:i + 3] if len(codon) == 3 and 'N' not in codon: counts[codon] += 1 return counts def compute_rscu(codon_counts, codon_table): """RSCU_ij = observed count / expected count under uniform synonymous usage.""" aa_groups = defaultdict(list) for codon, aa in codon_table.items(): if aa != '*': aa_groups[aa].append(codon) rscu = {} for aa, synonyms in aa_groups.items(): total = sum(codon_counts.get(c, 0) for c in synonyms) expected = total / len(synonyms) for codon in synonyms: rscu[codon] = (codon_counts.get(codon, 0) / expected) if expected > 0 else 0.0 return rscu def compute_cai(gene_sequence, rscu_reference, codon_table): """CAI = geometric mean of w_i = RSCU_i / max(RSCU) within each synonymous family.""" aa_groups = defaultdict(list) for codon, aa in codon_table.items(): if aa != '*': aa_groups[aa].append(codon) max_rscu = {aa: max(rscu_reference.get(c, 0) for c in syns) for aa, syns in aa_groups.items()} w = {c: (rscu_reference.get(c, 0) / max_rscu[aa] if max_rscu[aa] > 0 else 0) for c, aa in codon_table.items() if aa != '*'} log_w = [] seq = gene_sequence.upper() for i in range(0, len(seq) - 2, 3): codon = seq[i:i + 3] wi = w.get(codon, 0) if wi > 0: log_w.append(math.log(wi)) return math.exp(sum(log_w) / len(log_w)) if log_w else 0.0 test_seq = 'TAATGCCCGAATTTGCCTAAATGGGCAAATAG' for frame in range(3): print(f'Frame {frame}: {translate_frame(test_seq, codon_table, frame) or "(no ORF)"}') ``` ## Virtual Restriction Digest and ORF Finding **Goal:** Predict fragment sizes from a multi-enzyme restriction digest and locate ORFs in all 6 reading frames. **Approach:** Store each enzyme as `{site, cut}` (cut offset from the start of the recognition site on the top strand), scan for all occurrences, sort cut positions, and take consecutive differences as fragment sizes. ORF finding scans both strands in 3 frames for ATG...stop runs. ```python restriction_enzymes = { 'EcoRI': {'site': 'GAATTC', 'cut': 1}, 'BamHI': {'site': 'GGATCC', 'cut': 1}, 'HindIII': {'site': 'AAGCTT', 'cut': 1}, 'SmaI': {'site': 'CCCGGG', 'cut': 3}, 'NotI': {'site': 'GCGGCCGC', 'cut': 2}, } def reverse_complement(seq): """Reverse complement of a DNA sequence.""" return seq.upper().translate(str.maketrans('ACGT', 'TGCA'))[::-1] def find_cut_sites(sequence, enzyme_info): """Return absolute cut positions for one enzyme in a linear sequence.""" site, cut_offset = enzyme_info['site'], enzyme_info['cut'] seq_upper = sequence.upper() cuts, start = [], 0 while True: pos = seq_upper.find(site, start) if pos == -1: break cuts.append(pos + cut_offset) start = pos + 1 return cuts def virtual_digest(sequence, selected_enzymes, restriction_enzymes): """Multi-enzyme virtual digest of a linear sequence -> fragment sizes (largest first).""" all_cuts = sorted( (cut, name) for name in selected_enzymes for cut in find_cut_sites(sequence, restriction_enzymes[name]) ) boundaries = [0] + [c for c, _ in all_cuts] + [len(sequence)] fragments = [boundaries[i + 1] - boundaries[i] for i in range(len(boundaries) - 1) if boundaries[i + 1] > boundaries[i]] return sorted(fragments, reverse=True), all_cuts def find_orfs(sequence, codon_table, min_length=30): """Find ORFs (ATG...stop) across all 6 reading frames; returns list of dicts sorted by length.""" seq = sequence.upper() stops = {c for c, a in codon_table.items() if a == '*'} def scan(strand_seq, label): orfs = [] for frame in range(3): i, orf_start = frame, None while i <= len(strand_seq) - 3: codon = strand_seq[i:i + 3] if codon == 'ATG' and orf_start is None: orf_start = i elif codon in stops and orf_start is not None: if i - orf_start >= min_length: protein = translate(strand_seq[orf_start:i], codon_table, stop_at_stop=False) orfs.append({'start': orf_start, 'end': i + 3, 'frame': frame + 1, 'strand': label, 'length_nt': i - orf_start, 'protein': protein}) orf_start = None i += 3 return orfs return sorted(scan(seq, '+') + scan(reverse_complement(seq), '-'), key=lambda o: -o['length_nt']) plasmid = 'GAATTC' + 'A' * 500 + 'GGATCC' + 'A' * 800 + 'AAGCTT' + 'A' * 300 + 'GAATTC' + 'A' * 400 fragments, cuts = virtual_digest(plasmid, ['EcoRI', 'BamHI', 'HindIII'], restriction_enzymes) print(f'Fragments (bp): {fragments}, total={sum(fragments)} (== {len(plasmid)})') ``` ## Classical Genetics: Three-Point Cross and Hardy-Weinberg **Goal:** Order genes and compute map distances from three-point testcross progeny; test genotype counts for Hardy-Weinberg equilibrium. **Approach:** The rarest progeny classes are double crossovers — the gene that differs between the two double-CO classes and the parental classes is the middle gene. Map distance (cM) = (single COs in that region + double COs) / total x 100. For HWE, estimate allele frequencies from observed genotype counts, compute expected Hardy-Weinberg counts, and run a chi-squared goodness-of-fit test (1 df for a biallelic locus). ```python from scipy.stats import chi2 as chi2_dist three_point_data = {'+++': 580, 'abc': 592, 'ab+': 85, '++c': 83, 'a++': 218, '+bc': 214, 'a+c': 22, '+b+': 18} total = sum(three_point_data.values()) # 'a+c' and '+b+' are the smallest (double-CO) classes -> b is the middle gene double_co = three_point_data['a+c'] + three_point_data['+b+'] dist_a_b = (three_point_data['a++'] + three_point_data['+bc'] + double_co) / total * 100 dist_b_c = (three_point_data['ab+'] + three_point_data['++c'] + double_co) / total * 100 print(f'Order a-b-c: a-b = {dist_a_b:.1f} cM, b-c = {dist_b_c:.1f} cM') expected_double_co = (dist_a_b / 100) * (dist_b_c / 100) * total coc = double_co / expected_double_co print(f'Coefficient of coincidence = {coc:.3f}, interference = {1 - coc:.3f}') def hwe_chi_squared(n_AA, n_Aa, n_aa): """Estimate allele freqs, expected HWE genotype counts, and chi-squared test (df=1).""" n_total = n_AA + n_Aa + n_aa p = (2 * n_AA + n_Aa) / (2 * n_total) q = 1 - p expected = {'AA': p ** 2 * n_total, 'Aa': 2 * p * q * n_total, 'aa': q ** 2 * n_total} observed = {'AA': n_AA, 'Aa': n_Aa, 'aa': n_aa} chi2_stat = sum((observed[g] - expected[g]) ** 2 / expected[g] for g in observed) p_value = chi2_dist.sf(chi2_stat, df=1) return p, q, chi2_stat, p_value p, q, chi2_stat, p_value = hwe_chi_squared(320, 480, 200) print(f'p={p:.3f}, q={q:.3f}, chi2={chi2_stat:.3f}, p-value={p_value:.4f}') print('Deviates from HWE' if p_value < 0.05 else 'Consistent with HWE') ``` ## Pitfalls - **Undefined `translate()`**: naive `translate_frame` implementations call a `translate()` helper that was never defined — always define it explicitly (shown above) rather than assuming Biopython's `Seq.translate()` is in scope. - **Coordinate systems**: BED is 0-based half-open; VCF/GFF/genetic-map cM positions are 1-based inclusive — mixing them causes off-by-one errors in cut sites and ORF coordinates. - **CAI reference set**: CAI is only meaningful relative to a reference set of highly expressed genes (ribosomal proteins, glycolytic enzymes) from the *same organism* — using the wrong organism or a non-HEG reference gives misleading scores. - **Double crossovers**: the smallest progeny classes in a three-point cross are the double-CO classes; misidentifying them flips the inferred gene order. - **HWE p-value approximation**: without SciPy, chi-squared p-values need a proper CDF (Wilson-Hilferty or similar) — a crude approximation can flip significance calls near p = 0.05. ## See Also - `bio-sequence-manipulation-codon-usage` - `bio-restriction-analysis-restriction-mapping` - `bio-population-genetics-population-structure` - `bio-population-genetics-selection-statistics`