bio-core-computational-genetics

name	bio-core-computational-genetics
description	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.
tool_type	python
primary_tool	SciPy

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)


 ():
    
    seq = coding_sequence.upper()
    counts = Counter()
     i  (, (seq) - , ):
        codon = seq[i:i + ]
         (codon) ==      codon:
            counts[codon] += 
     counts


 ():
    
    aa_groups = defaultdict()
     codon, aa  codon_table.items():
         aa != :
            aa_groups[aa].append(codon)
    rscu = {}
     aa, synonyms  aa_groups.items():
        total = (codon_counts.get(c, )  c  synonyms)
        expected = total / (synonyms)
         codon  synonyms:
            rscu[codon] = (codon_counts.get(codon, ) / expected)  expected >   
     rscu


 ():
    
    aa_groups = defaultdict()
     codon, aa  codon_table.items():
         aa != :
            aa_groups[aa].append(codon)
    max_rscu = {aa: (rscu_reference.get(c, )  c  syns)  aa, syns  aa_groups.items()}
    w = {c: (rscu_reference.get(c, ) / max_rscu[aa]  max_rscu[aa] >   )
          c, aa  codon_table.items()  aa != }

    log_w = []
    seq = gene_sequence.upper()
     i  (, (seq) - , ):
        codon = seq[i:i + ]
        wi = w.get(codon, )
         wi > :
            log_w.append(math.log(wi))
     math.exp((log_w) / (log_w))  log_w  


test_seq = 
 frame  ():
    ()

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 ():
    
    all_cuts = (
        (cut, name)
         name  selected_enzymes
         cut  find_cut_sites(sequence, restriction_enzymes[name])
    )
    boundaries = [] + [c  c, _  all_cuts] + [(sequence)]
    fragments = [boundaries[i + ] - boundaries[i]  i  ((boundaries) - )  boundaries[i + ] > boundaries[i]]
     (fragments, reverse=), all_cuts


 ():
    
    seq = sequence.upper()
    stops = {c  c, a  codon_table.items()  a == }

     ():
        orfs = []
         frame  ():
            i, orf_start = frame, 
             i <= (strand_seq) - :
                codon = strand_seq[i:i + ]
                 codon ==   orf_start  :
                    orf_start = i
                 codon  stops  orf_start   :
                     i - orf_start >= min_length:
                        protein = translate(strand_seq[orf_start:i], codon_table, stop_at_stop=)
                        orfs.append({: orf_start, : i + , : frame + ,
                                     : label, : i - orf_start, : protein})
                    orf_start = 
                i += 
         orfs

     (scan(seq, ) + scan(reverse_complement(seq), ), key= o: -o[])


plasmid =  +  *  +  +  *  +  +  *  +  +  * 
fragments, cuts = virtual_digest(plasmid, [, , ], restriction_enzymes)
()

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 = ( * n_AA + n_Aa) / ( * n_total)
    q =  - p
    expected = {: p **  * n_total, :  * p * q * n_total, : q **  * n_total}
    observed = {: n_AA, : n_Aa, : n_aa}
    chi2_stat = ((observed[g] - expected[g]) **  / expected[g]  g  observed)
    p_value = chi2_dist.sf(chi2_stat, df=)
     p, q, chi2_stat, p_value


p, q, chi2_stat, p_value = hwe_chi_squared(, , )
()
(  p_value <   )

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