Bio Applied Population Genetics

Population Genetics and Molecular Evolution

When to Use

• Testing whether observed genotype counts (e.g. from a SNP panel) deviate from Hardy-Weinberg equilibrium, or computing an inbreeding coefficient F
• Simulating genetic drift (Wright-Fisher), directional/balancing/frequency-dependent selection, or bottleneck effects on allele frequency
• Estimating divergence time via the molecular clock, Jukes-Cantor distance, or dN/dS (Nei-Gojobori) on coding sequences
• Computing neutrality statistics (Tajima's D, McDonald-Kreitman/alpha) or population differentiation (Fst) and linkage disequilibrium (r²) from allele frequencies or haplotypes
• The request mentions "allele frequency," "selection coefficient," "effective population size," "dN/dS," "Tajima's D," "Fst," or "linkage disequilibrium" without a specific VCF/PLINK workflow already in hand

Version Compatibility

Python ≥3.10, NumPy ≥1.24, SciPy ≥1.10, pandas ≥2.0. Pure-Python/NumPy math — no genomics-specific package versions to track. For real VCF/PLINK-scale data, pair with scikit-allel or PLINK (see See Also).

Prerequisites

pip install numpy scipy pandas matplotlib. Assumes basic diploid genetics (alleles, genotypes) and familiarity with numpy.random.Generator (rng = np.random.default_rng(seed)).

Hardy-Weinberg Equilibrium and Inbreeding

Goal: decide whether a biallelic locus's genotype counts match random-mating expectations, and quantify departure with an inbreeding coefficient. Approach: derive allele frequencies from genotype counts, compute expected counts under HWE, run a chi-squared test (df=1 for a biallelic locus), then F = 1 − H_obs/H_exp.

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 (p_allele_freq, chi2_stat, p_value, F_inbreeding).
    F > 0 indicates excess homozygotes (inbreeding or Wahlund effect);
    F < 0 indicates excess heterozygotes.
    """
    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_stat = np.sum((obs - exp) ** 2 / exp)
    p_value = stats.chi2.sf(chi2_stat, df=1)
    F = 1 - (obs_Aa / n) / (2 * p * q) if p * q > 0 else 0.0
    return p, chi2_stat, p_value, F


p, chi2, pval, F = hwe_test(obs_AA=233, obs_Aa=200, obs_aa=67)  # MN blood group, n=500
print(f"p(M)={p:.4f}  chi2={chi2:.3f}  p-value={pval:.4f}  F={F:.4f}")

Wright-Fisher Drift and Selection

Goal: project how an allele frequency changes across generations under drift alone, or under drift + directional/balancing selection. Approach: each generation, draw the post-selection allele count from a Binomial(2N, p_selected); track until fixation (p=1) or loss (p=0), which are absorbing states.

import numpy as np


def wright_fisher_trajectory(N, p0, n_generations, rng, s=0.0, h=0.5):
    """Simulate one Wright-Fisher trajectory with optional additive selection.

    N: diploid population size. p0: initial allele frequency of A.
    s: selection coefficient (fitness of aa = 1-s, AA = 1, Aa = 1-h*s).
    Returns an array of allele frequencies of length n_generations+1.
    """
    freq = np.empty(n_generations + 1)
    freq[0] = p0
    p = p0
    for g in range(n_generations):
        q = 1 - p
        if s != 0.0:
            w_bar = p**2 + 2 * p * q * (1 - h * s) + q**2 * (1 - s)
            p = (p**2 + p * q * (1 - h * s)) / w_bar   # deterministic selection step
        k = rng.binomial(2 * N, p)                      # then binomial drift
        p = k / (2 * N)
        freq[g + 1] = p
        if p in (0.0, 1.0):
            freq[g + 2:] = p   # absorbing state: stays fixed/lost
            break
    return freq


rng = np.random.default_rng(42)
neutral = wright_fisher_trajectory(N=100, p0=0.3, n_generations=200, rng=rng)
selected = wright_fisher_trajectory(N=100, p0=0.01, n_generations=200, rng=rng, s=0.1, h=0.5)
# Neutral: P(fixation) = p0 (Kimura 1962). Selection dominates drift when 2*N*s >> 1.

dN/dS, Tajima's D, Fst, and LD

Goal: quantify selection pressure on coding sequences (dN/dS), detect non-neutral allele frequency spectra (Tajima's D), and measure population differentiation (Fst) or co-inheritance of variants (LD). Approach: JC69-correct raw p-distances before rate ratios; Tajima's D contrasts two theta estimators (pi vs Watterson); Fst contrasts within- vs total heterozygosity; r² is the standard LD measure.

import numpy as np
from itertools import combinations


def jukes_cantor_distance(p_diff):
    """JC69-corrected distance; multiple hits make this unreliable above p≈0.75."""
    arg = 1.0 - (4.0 / 3.0) * p_diff
    return np.inf if arg <= 0 else -(3.0 / 4.0) * np.log(arg)


def tajimas_d(seqs):
    """Tajima's D from a list of equal-length aligned sequence strings.
    D < 0: excess rare variants (sweep/expansion). D > 0: excess intermediate
    variants (balancing selection/bottleneck). Returns (pi, theta_W, D, S)."""
    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
    pi = np.mean([sum(a != b for a, b in zip(seqs[i], seqs[j]))
                  for i, j in combinations(range(n), 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, c2 = b1 - 1 / a1, b2 - (n + 2) / (a1 * n) + a2 / a1**2
    var_d = (c1 / a1) * S + (c2 / (a1**2 + a2)) * 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(p1, p2, n1=1, n2=1):
    """Fst between two populations at a biallelic locus (Wright's fixation index)."""
    p_bar = (n1 * p1 + n2 * p2) / (n1 + n2)
    h_t = 2 * p_bar * (1 - p_bar)
    h_s = (2 * p1 * (1 - p1) * n1 + 2 * p2 * (1 - p2) * n2) / (n1 + n2)
    return (h_t - h_s) / h_t if h_t > 0 else 0.0


def ld_r2(hap_a, hap_b):
    """r^2 LD between two 0/1 haplotype-indicator arrays of equal length."""
    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

Pitfalls

• HWE departure does not identify the cause: chi-squared only flags departure; genotyping errors, selection, inbreeding, and stratification all produce the same test statistic
• Wahlund effect: pooling samples from differentiated subpopulations creates apparent HWE departure (excess homozygotes) without any real selection or inbreeding — stratify by population first
• Wright-Fisher assumes non-overlapping generations: no age structure; use coalescent models for organisms with overlapping generations
• Jukes-Cantor breaks down above p ≈ 0.5–0.75: multiple substitutions per site saturate the correction; use HKY/GTR for highly diverged sequences
• dN/dS > 1 needs many sites to detect: gene-wide averaging masks episodic selection at specific codons — use PAML codeml or HyPhy for site-specific tests
• Tajima's D and MK test are both confounded by demography: population expansion mimics a sweep (D < 0); interpret alongside a demographic null model, not in isolation

See Also

• bio-population-genetics-selection-statistics — genome-scale selection scans (iHS, XP-EHH, nSL) beyond single-locus dN/dS
• bio-population-genetics-population-structure — PCA/ADMIXTURE-based structure inference on real genotype matrices
• bio-population-genetics-linkage-disequilibrium — LD decay, haplotype blocks, and pairwise D'/r² at genome scale
• bio-population-genetics-scikit-allel-analysis — running these same statistics directly on VCF-derived genotype arrays