Bio Applied Molecular Evolution

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
    return np.array(freqs[:n_gen + 1])

rng = np.random.default_rng(42)
neutral = wright_fisher_trajectory(N=100, p0=0.3, n_gen=200, rng=rng)
selected = wf_selection_drift(N=100, p0=0.01, s=0.05, h=0.5, n_gen=200, 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':'Val','GTC':'Val','GTA':'Val','GTG':'Val',
    'TCT':'Ser','TCC':'Ser','TCA':'Ser','TCG':'Ser','CCT':'Pro','CCC':'Pro','CCA':'Pro','CCG':'Pro',
    'ACT':'Thr','ACC':'Thr','ACA':'Thr','ACG':'Thr','GCT':'Ala','GCC':'Ala','GCA':'Ala','GCG':'Ala',
    'TAT':'Tyr','TAC':'Tyr','TAA':'Stop','TAG':'Stop','CAT':'His','CAC':'His','CAA':'Gln','CAG':'Gln',
    'AAT':'Asn','AAC':'Asn','AAA':'Lys','AAG':'Lys','GAT':'Asp','GAC':'Asp','GAA':'Glu','GAG':'Glu',
    'TGT':'Cys','TGC':'Cys','TGA':'Stop','TGG':'Trp','CGT':'Arg','CGC':'Arg','CGA':'Arg','CGG':'Arg',
    'AGT':'Ser','AGC':'Ser','AGA':'Arg','AGG':'Arg','GGT':'Gly','GGC':'Gly','GGA':'Gly','GGG':'Gly',
}

def count_syn_sites(codon):
    """Fraction of a codon's 3 sites that are synonymous (Nei-Gojobori averaging)."""
    if codon not in GENETIC_CODE or GENETIC_CODE[codon] == 'Stop':
        return 0.0, 3.0
    aa_orig = GENETIC_CODE[codon]
    syn_sites = 0.0
    for pos in range(3):
        syn_count = total_nonstop = 0
        for base in _BASES:
            if base == codon[pos]:
                continue
            mutant = codon[:pos] + base + codon[pos + 1:]
            if mutant not in GENETIC_CODE or GENETIC_CODE[mutant] == 'Stop':
                continue
            total_nonstop += 1
            if GENETIC_CODE[mutant] == aa_orig:
                syn_count += 1
        if total_nonstop > 0:
            syn_sites += syn_count / total_nonstop
    return syn_sites, 3.0 - syn_sites

def pairwise_dn_ds(seq1, seq2):
    """dN, dS, omega between two aligned in-frame coding sequences (Nei-Gojobori 1986).
    omega < 1: purifying selection; omega ~ 1: neutral; omega > 1: positive selection."""
    assert len(seq1) == len(seq2) and len(seq1) % 3 == 0
    S_total = N_total = Sd = Nd = 0.0
    for i in range(0, len(seq1), 3):
        c1, c2 = seq1[i:i + 3], seq2[i:i + 3]
        if c1 not in GENETIC_CODE or c2 not in GENETIC_CODE:
            continue
        if GENETIC_CODE[c1] == 'Stop' or GENETIC_CODE[c2] == 'Stop':
            continue
        s1, n1 = count_syn_sites(c1)
        s2, n2 = count_syn_sites(c2)
        S_total += (s1 + s2) / 2
        N_total += (n1 + n2) / 2
        n_diffs = sum(c1[j] != c2[j] for j in range(3))
        if n_diffs == 0:
            continue
        if GENETIC_CODE[c1] == GENETIC_CODE[c2]:
            Sd += n_diffs
        else:
            Nd += n_diffs   # simplified: multi-hit codons assign all diffs by final aa change
    pS = Sd / S_total if S_total > 0 else 0.0
    pN = Nd / N_total if N_total > 0 else 0.0
    dS = jukes_cantor_distance(pS)
    dN = jukes_cantor_distance(pN)
    omega = dN / dS if dS > 0 else np.inf
    return {'S': S_total, 'N': N_total, 'dS': dS, 'dN': dN, 'omega': 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 + 1) / (3 * (n - 1))
    b2 = 2 * (n**2 + n + 3) / (9 * n * (n - 1))
    c1 = b1 - 1 / a1
    c2 = b2 - (n + 2) / (a1 * n) + a2 / a1**2
    e1, e2 = c1 / a1, c2 / (a1**2 + a2)
    var_d = e1 * S + e2 * S * (S - 1)

    D = (pi - theta_W) / np.sqrt(var_d) if var_d > 0 else 0.0
    return pi, theta_W, D, S

def fst(allele_freqs):
    """Fst for a biallelic locus across N populations (equal sample sizes).
    Wright's scale: <0.05 little, 0.05-0.15 moderate, 0.15-0.25 great, >0.25 very great."""
    freqs = np.asarray(allele_freqs, dtype=float)
    H_S = np.mean(2 * freqs * (1 - freqs))
    p_bar = freqs.mean()
    H_T = 2 * p_bar * (1 - p_bar)
    return (H_T - H_S) / H_T if H_T > 0 else 0.0

def ld_r2(hap_a, hap_b):
    """r^2 linkage disequilibrium between two 0/1-coded haplotype vectors."""
    p_A, p_B = hap_a.mean(), hap_b.mean()
    D = np.mean(hap_a * hap_b) - p_A * p_B
    denom = p_A * (1 - p_A) * p_B * (1 - p_B)
    return D**2 / denom if denom > 0 else 0.0

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