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.

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.

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).

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