bio-applied-molecular-evolution

name	bio-applied-molecular-evolution
description	Test Hardy-Weinberg equilibrium, simulate Wright-Fisher drift/selection, and compute dN/dS, Tajima's D, and Fst with NumPy/SciPy. Use for neutral theory, molecular clock divergence time, selection scans, or effective population size (Ne) questions.
tool_type	python
primary_tool	numpy/scipy

Population Genetics and Molecular Evolution

When to Use

• Testing whether observed genotype counts (e.g., a SNP panel) depart from Hardy-Weinberg equilibrium
• Simulating genetic drift, bottlenecks, or selection trajectories under the Wright-Fisher model
• Estimating divergence time from sequence identity via the molecular clock (Jukes-Cantor)
• Computing dN/dS (omega) on aligned coding sequences to detect purifying/positive selection
• Scanning for selection footprints with Tajima's D, McDonald-Kreitman, or Fst/LD across populations

Version Compatibility

Python ≥3.10, NumPy ≥1.24, SciPy ≥1.11, pandas ≥2.0. No specialized package required — all methods below are pure NumPy/SciPy. For production-grade dN/dS (site models, branch-site tests) use PAML ≥4.10 or HyPhy ≥2.5 instead of the simplified Nei-Gojobori code here.

Prerequisites

• pip install numpy scipy pandas matplotlib
• Familiarity with allele frequencies, diploid genotypes, and codon tables
• For real data: aligned coding sequences (FASTA, in-frame, gap-free) or VCF-derived allele counts per population

Hardy-Weinberg Equilibrium

Goal: decide whether a biallelic locus's genotype counts are consistent with random mating (HWE), and quantify inbreeding if not. Approach: estimate p, q from counts, compute expected genotype counts under p²+2pq+q², run a 1-df chi-squared test, and compute the inbreeding coefficient F.

import numpy as np
from scipy import stats

def hwe_test(obs_AA, obs_Aa, obs_aa):
    """Chi-squared test for Hardy-Weinberg equilibrium at a biallelic locus.

    Returns (chi2_stat, p_value, F) where F is the inbreeding coefficient
    (F ~ 0: no inbreeding; F > 0: heterozygote deficit; F < 0: heterozygote excess).
    """
    n = obs_AA + obs_Aa + obs_aa
    p = (2 * obs_AA + obs_Aa) / (2 * n)
    q = 1 - p
    exp_AA, exp_Aa, exp_aa = p**2 * n, 2 * p * q * n, q**2 * n
    obs = np.array([obs_AA, obs_Aa, obs_aa], dtype=float)
    exp = np.array([exp_AA, exp_Aa, exp_aa], dtype=float)
    chi2 = np.sum((obs - exp) ** 2 / exp)
    p_value = stats.chi2.sf(chi2, df=1)   # df = classes(3) - 1 - allele_freqs_estimated(1)
    F = 1 - (obs_Aa / n) / (2 * p * q)
    return chi2, p_value, F

# Example: MN blood group, n=500 (233 MM, 200 MN, 67 NN)
chi2, p_value, F = hwe_test(233, 200, 67)
print(f"chi2={chi2:.3f} p={p_value:.4f} F={F:.4f}")

Wright-Fisher Drift and Selection

Goal: simulate how an allele frequency evolves under finite population size (drift), with optional directional selection and population-size bottlenecks. Approach: each generation draws the next allele count from Binomial(2N, p); selection reweights p before sampling; absorbing states (p=0 or p=1) short-circuit the loop.

import numpy as np

def wright_fisher_trajectory(N, p0, n_gen, rng):
    """Neutral Wright-Fisher trajectory. Returns array of length n_gen+1.
    Absorbing states: p=0 (loss) and p=1 (fixation)."""
    freq = np.empty(n_gen + 1)
    freq[0] = p = p0
    for g in range(n_gen):
        p = rng.binomial(2 * N, p) / (2 * N)
        freq[g + 1] = p
        if p == 0.0 or p == 1.0:
            freq[g + 2:] = p
            break
    return freq

def wf_selection_drift(N, p0, s, h, n_gen, rng):
    """Wright-Fisher with selection then drift. aa fitness=1-s, Aa fitness=1-h*s, AA fitness=1.
    Selection dominates drift when 2*N*s >> 1; drift dominates when 2*N*s << 1."""
    p = p0
    freqs = [p]
    for _ in range(n_gen):
        q = 1 - p
        w_bar = p**2 + 2 * p * q * (1 - h * s) + q**2 * (1 - s)
        p_sel = (p**2 + p * q * (1 - h * s)) / w_bar
        p = rng.binomial(2 * N, p_sel) / (2 * N)
        freqs.append(p)
        if p in (0.0, 1.0):
            freqs.extend([p] * (n_gen - len(freqs) + 1))
            break
     np.array(freqs[:n_gen + ])

rng = np.random.default_rng()
neutral = wright_fisher_trajectory(N=, p0=, n_gen=, rng=rng)
selected = wf_selection_drift(N=, p0=, s=, h=, n_gen=, rng=rng)

Effective population size Ne (always ≤ census N):

• Unequal sex ratio: Ne = 4*Nm*Nf / (Nm + Nf)
• Fluctuating size (bottlenecks): Ne = len(N_per_gen) / sum(1/N_per_gen) (harmonic mean) — a few generations of crash dominate long-run drift more than the arithmetic mean would suggest
• Human Ne ≈ 10,000 despite ~8 billion census size (bottlenecks + recent growth)

Molecular Clock, dN/dS, and Neutrality Statistics

Goal: estimate divergence time from raw sequence identity, quantify selection pressure on coding sequences (dN/dS), and detect selection footprints from a sample alignment (Tajima's D, Fst). Approach: correct raw p-distance for multiple hits (Jukes-Cantor), count synonymous/non-synonymous sites and differences per codon (Nei-Gojobori) for omega, and use segregating-site/pairwise-difference estimators of theta for Tajima's D.

import numpy as np

def jukes_cantor_distance(p_distance):
    """JC69-corrected distance from raw proportion of differing sites.
    Diverges (returns inf) at p_distance >= 0.75 (saturation)."""
    arg = 1.0 - (4.0 / 3.0) * p_distance
    if arg <= 0:
        return np.inf
    return -(3.0 / 4.0) * np.log(arg)

# Divergence time: T = d_JC / (2*mu), mu = substitution rate per site per generation
def divergence_time(p_distance, mu):
    return jukes_cantor_distance(p_distance) / (2 * mu)

# --- Nei-Gojobori dN/dS ---
_BASES = ['A', 'T', 'C', 'G']
GENETIC_CODE = {
    'TTT':'Phe','TTC':'Phe','TTA':'Leu','TTG':'Leu','CTT':'Leu','CTC':'Leu','CTA':'Leu','CTG':'Leu',
    'ATT':'Ile','ATC':'Ile','ATA':'Ile','ATG':'Met','GTT':,:,:,:,
    :,:,:,:,:,:,:,:,
    :,:,:,:,:,:,:,:,
    :,:,:,:,:,:,:,:,
    :,:,:,:,:,:,:,:,
    :,:,:,:,:,:,:,:,
    :,:,:,:,:,:,:,:,
}

 ():
    
     codon   GENETIC_CODE  GENETIC_CODE[codon] == :
         , 
    aa_orig = GENETIC_CODE[codon]
    syn_sites = 
     pos  ():
        syn_count = total_nonstop = 
         base  _BASES:
             base == codon[pos]:
                
            mutant = codon[:pos] + base + codon[pos + :]
             mutant   GENETIC_CODE  GENETIC_CODE[mutant] == :
                
            total_nonstop += 
             GENETIC_CODE[mutant] == aa_orig:
                syn_count += 
         total_nonstop > :
            syn_sites += syn_count / total_nonstop
     syn_sites,  - syn_sites

 ():
    
     (seq1) == (seq2)  (seq1) %  == 
    S_total = N_total = Sd = Nd = 
     i  (, (seq1), ):
        c1, c2 = seq1[i:i + ], seq2[i:i + ]
         c1   GENETIC_CODE  c2   GENETIC_CODE:
            
         GENETIC_CODE[c1] ==   GENETIC_CODE[c2] == :
            
        s1, n1 = count_syn_sites(c1)
        s2, n2 = count_syn_sites(c2)
        S_total += (s1 + s2) / 
        N_total += (n1 + n2) / 
        n_diffs = (c1[j] != c2[j]  j  ())
         n_diffs == :
            
         GENETIC_CODE[c1] == GENETIC_CODE[c2]:
            Sd += n_diffs
        :
            Nd += n_diffs   
    pS = Sd / S_total  S_total >   
    pN = Nd / N_total  N_total >   
    dS = jukes_cantor_distance(pS)
    dN = jukes_cantor_distance(pN)
    omega = dN / dS  dS >   np.inf
     {: S_total, : N_total, : dS, : dN, : omega}

Goal: compute Tajima's D and Fst to flag selection/demographic signals across a sample or across populations. Approach: Tajima's D compares two theta estimators (pairwise pi vs. Watterson) using Tajima 1989's exact variance; Fst compares within- vs. total-population heterozygosity.

def tajimas_d(seqs):
    """Tajima's D (1989) for a list of equal-length aligned DNA strings.
    Returns (pi, theta_W, D, S). D<0: excess rare variants (sweep/expansion);
    D>0: excess intermediate-freq variants (balancing selection/bottleneck)."""
    n, L = len(seqs), len(seqs[0])
    S = sum(1 for j in range(L) if len({s[j] for s in seqs}) > 1)
    if S == 0:
        return 0.0, 0.0, 0.0, 0

    total_diffs = sum(sum(seqs[i][k] != seqs[j][k] for k in range(L))
                       for i in range(n) for j in range(i + 1, n))
    pi = total_diffs / (n * (n - 1) / 2)

    a1 = sum(1 / k for k in range(1, n))
    a2 = sum(1 / k**2 for k in range(1, n))
    theta_W = S / a1

    b1 = (n + ) / ( * (n - ))
    b2 =  * (n** + n + ) / ( * n * (n - ))
    c1 = b1 -  / a1
    c2 = b2 - (n + ) / (a1 * n) + a2 / a1**
    e1, e2 = c1 / a1, c2 / (a1** + a2)
    var_d = e1 * S + e2 * S * (S - )

    D = (pi - theta_W) / np.sqrt(var_d)  var_d >   
     pi, theta_W, D, S

 ():
    
    freqs = np.asarray(allele_freqs, dtype=)
    H_S = np.mean( * freqs * ( - freqs))
    p_bar = freqs.mean()
    H_T =  * p_bar * ( - p_bar)
     (H_T - H_S) / H_T  H_T >   

 ():
    
    p_A, p_B = hap_a.mean(), hap_b.mean()
    D = np.mean(hap_a * hap_b) - p_A * p_B
    denom = p_A * ( - p_A) * p_B * ( - p_B)
     D** / denom  denom >   

McDonald-Kreitman neutrality index (real data — Drosophila Adh, McDonald & Kreitman 1991):

from scipy import stats

Ps, Pn, Ds, Dn = 43, 2, 17, 7   # synonymous/non-syn polymorphic and fixed-difference counts
NI = (Pn / Ps) / (Dn / Ds)       # NI < 1: positive selection; NI > 1: slightly deleterious polymorphism
alpha = 1 - (Ds * Pn) / (Dn * Ps)
_, p_fisher = stats.fisher_exact([[Pn, Ps], [Dn, Ds]])

Pitfalls

• HWE departures from genotyping error: batch-specific or low-quality genotype calls produce systematic HWE failures — check per-batch before biological interpretation
• Ne ≠ census N: using census size in drift calculations severely underestimates drift; human Ne ≈ 10,000
• Tajima's D is confounded: population expansion mimics a sweep (D<0) and a bottleneck mimics balancing selection (D>0) — always interpret alongside a demographic null model
• JC correction breaks down: assumes equal base frequencies and independent single substitutions; returns inf for p ≥ 0.75 and is unreliable above p ≈ 0.6
• dN/dS gene-averaging hides episodic selection: a handful of positively selected codons can be masked by genome-wide omega < 1 — use PAML/HyPhy site or branch-site models, not this simplified per-codon count, for publication-grade inference
• Multiple testing: apply Benjamini-Hochberg when scanning many loci for HWE, Tajima's D, or Fst outliers

See Also

• bio-applied-population-genetics — allele-frequency workflows, PCA/admixture, and structure inference at scale
• bio-core-phylogenetics — tree inference and distance methods that consume Jukes-Cantor-style corrected distances
• bio-core-comparative-genomics — ortholog identification and synteny for setting up dN/dS comparisons
• bio-applied-gwas — association testing once population structure (Fst) is accounted for