Bio Population Genetics Selection Statistics

Version Compatibility

Reference examples tested with: scikit-allel 1.3+, numpy 1.26+, selscan 2.0+, SweepFinder2 1.0+.

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures
• CLI: <tool> --version then <tool> --help to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Version traps that change results, not just syntax: allel.standardize_by_allele_count(score, aac, ...) is for iHS and nSL (it bins by derived-allele count), while XP-EHH and iHH12 use plain genome-wide allel.standardize(score); conflating them is a real bug. allel.nsl(h) takes no pos/map_pos (nSL is map-free by construction) whereas allel.ihs(h, pos, map_pos=...) distorts without a genetic map. selscan 2.0 adds --unphased for multilocus genotypes; selscan 1.x requires phased haplotypes. The single source of truth for versions is this block, not headings.

Selection Statistics

"Scan my population for signatures of natural selection" -> Contrast each locus against the genome-wide neutral expectation that already carries the demographic history, using statistics that read complementary features of a sweep.

• Python: allel.ihs(), allel.xpehh(), allel.nsl(), allel.garud_h(), allel.windowed_tajima_d(), allel.hudson_fst() (scikit-allel)
• CLI: selscan --ihs|--xpehh|--nsl then norm for standardization; SweepFinder2 -lrb for CLR with a background-selection map

Scope: selection-scan DESIGN - SFS neutrality tests, haplotype/EHH sweep statistics, differentiation outliers (FST, PBS), CLR scans, and polygenic Qx, with standardization and demography-aware outlier interpretation. PCA/ADMIXTURE population assignment routes to population-structure; LD/EHH mechanics to linkage-disequilibrium; the scikit-allel array-API mechanics (signatures, accessibility masks, data loading) to scikit-allel-analysis; phasing to phasing-imputation/haplotype-phasing; dN/dS and McDonald-Kreitman to comparative-genomics/positive-selection.

The Single Most Important Insight -- no single statistic separates selection from demography at one locus

1. A bottleneck removes variation genome-wide and skews the SFS toward rare variants exactly as a sweep does, so a single Tajima's D cannot tell them apart; recent expansion gives genome-wide negative D that mimics a post-sweep signal, and population structure gives positive D that mimics balancing selection.
2. The field's response is structural, not statistical: rank windows and flag empirical genome-wide OUTLIERS (the bulk absorbs the shared demography), intersect MULTIPLE ORTHOGONAL signals (a sweep distorts haplotype, SFS, diversity, and differentiation together while each confound acts differently), and calibrate against an explicit demographic-model SIMULATION (dadi/momi/fastsimcoal2 + msprime/SLiM).
3. iHS and nSL detect INCOMPLETE/ongoing sweeps and collapse to ~0 once the allele fixes; XP-EHH and Rsb catch FIXED/near-fixed sweeps by borrowing a second population, so they are COMPLEMENTARY not redundant (high XP-EHH with near-zero iHS is the textbook completed-sweep signature).
4. An absolute cutoff such as "Tajima's D < -2 means a sweep" is wrong by construction; -2 means nothing without the genome-wide context, which is why the entire skill is outlier-based.

Tool Taxonomy

Decision Tree by Scenario

FST and PBS - Differentiation

Goal: Quantify allele-frequency divergence between populations and flag local-adaptation candidates without being fooled by rare variants or unequal samples.

Approach: Count alleles per subpopulation, compute the Hudson estimator per SNP, and report mean FST as a ratio-of-averages over windows after an MAF filter (rare variants deflate FST, Bhatia 2013).

import allel
import numpy as np

callset = allel.read_vcf('data.vcf.gz')
gt = allel.GenotypeArray(callset['calldata/GT'])
pos = callset['variants/POS']

subpops = {'pop1': [0, 1, 2, 3, 4], 'pop2': [5, 6, 7, 8, 9]}
ac_subpops = gt.count_alleles_subpops(subpops)
ac1, ac2 = ac_subpops['pop1'], ac_subpops['pop2']

# MAF filter first: rare variants systematically deflate FST (Bhatia 2013).
maf = np.minimum(ac1.to_frequencies()[:, 1], ac2.to_frequencies()[:, 1])
keep = (maf > 0.05) | (1 - maf > 0.05)

# Hudson estimator: robust to unequal sample sizes. Mean FST = ratio-of-averages, never mean of per-SNP ratios.
num, den = allel.hudson_fst(ac1[keep], ac2[keep])
fst_mean = np.nansum(num) / np.nansum(den)

# Windowed scan for outlier localization.
fst_win, windows, n_snps = allel.windowed_hudson_fst(pos[keep], ac1[keep], ac2[keep], size=100000, step=50000)
outlier = fst_win > np.nanpercentile(fst_win, 99)

PBS (Population Branch Statistic, Yi 2010) needs three populations: convert three pairwise FST values to branch lengths via T = -log(1 - FST) and isolate the focal branch as PBS = (T_12 + T_13 - T_23) / 2. A long focal branch in a small/bottlenecked population can be drift, not selection.

iHS and nSL - Within-Population Incomplete Sweeps

Goal: Detect ongoing sweeps from extended haplotype homozygosity around derived core alleles.

Approach: Filter to segregating biallelic SNPs, compute the raw integrated-EHH score, then standardize WITHIN derived-allele-frequency bins (the raw score depends strongly on derived frequency and is uninterpretable unbinned).

import allel
import numpy as np

h = gt.to_haplotypes()
ac = h.count_alleles()
flt = (ac[:, 0] > 1) & (ac[:, 1] > 1)
h_flt, pos_flt, ac_flt = h.compress(flt, axis=0), pos[flt], ac.compress(flt, axis=0)

# iHS needs a genetic map (map_pos) where recombination rate varies; physical distance distorts iHH.
ihs_raw = allel.ihs(h_flt, pos_flt, min_maf=0.05, include_edges=True)
# Standardize WITHIN derived-allele-count bins (NOT genome-wide) - the most common iHS bug.
ihs_std, _ = allel.standardize_by_allele_count(ihs_raw, ac_flt[:, 1])
ihs_outlier = np.abs(ihs_std) > 2

# nSL: map-free (no pos/map_pos), robust to recombination-rate variation; still bin-standardized.
nsl_raw = allel.nsl(h_flt)
nsl_std, _ = allel.standardize_by_allele_count(nsl_raw, ac_flt[:, 1])

A null iHS does NOT mean no selection: once the derived allele fixes, the ancestral haplotype is gone and iHS collapses to zero. Switch to XP-EHH/Rsb for completed sweeps. Report the windowed PROPORTION of |iHS|>2 SNPs, not single noisy hits.

XP-EHH - Cross-Population Completed Sweeps

Goal: Catch fixed or near-fixed sweeps in one of two populations, where within-population haplotype tests are blind.

Approach: Compute the cross-population integrated-EHH ratio on shared positions, then apply a plain GENOME-WIDE z-score (XP-EHH is not strongly correlated with derived frequency, so frequency-bin standardization is the wrong rule).

import allel
import numpy as np

h1 = h_flt.take(pop1_hap_idx, axis=1)
h2 = h_flt.take(pop2_hap_idx, axis=1)

xpehh_raw = allel.xpehh(h1, h2, pos_flt, include_edges=True)
# Genome-wide z-score for XP-EHH/iHH12 - NOT standardize_by_allele_count (that is for iHS/nSL).
xpehh_std = allel.standardize(xpehh_raw)
completed_sweep = np.abs(xpehh_std) > 2

A shared ancestral sweep at the same locus in BOTH populations cancels in the contrast (false negative). Differential phasing quality between the two populations systematically biases the difference statistic.

H12 - Soft Sweeps

Goal: Recover sweeps on standing variation that iHS and CLR miss because no single long haplotype dominates.

Approach: Compute Garud's H over haplotype-frequency spectra; H12 collapses the two most common haplotypes into one, and H2/H1 tends (not classifies) toward hard-vs-soft.

h1, h12, h123, h2_h1 = allel.garud_h(h_flt)
h12_win = allel.moving_garud_h(h_flt, size=100)  # SNP-count windows, not bp

SFS-CLR with selscan and SweepFinder2

# selscan haplotype scans (phased VCF + genetic map), then norm for standardization.
selscan --ihs --vcf phased.vcf --map genetic.map --out scan        # within-pop incomplete sweeps
selscan --xpehh --vcf pop1.vcf --vcf-ref pop2.vcf --map genetic.map --out xp  # completed sweeps
selscan --nsl --vcf phased.vcf --out nsl_scan                       # no map; map-free
selscan --ihs --unphased --vcf unphased.vcf --map genetic.map --out scan_unphased  # selscan 2.0 multilocus-genotype

# norm bins iHS/nSL by derived-allele frequency (--bins default 100); XP-EHH/iHH12 get a genome-wide z-score.
norm --ihs --files scan.ihs.out --bins 100
norm --xpehh --files xp.xpehh.out

# SweepFinder2 CLR: build a genome-WIDE background SFS, then scan with recombination + B-value (BGS) map.
SweepFinder2 -f genome.freq genome.spect                            # empirical background spectrum
# -lrb takes N1 (current ingroup Ne), N2 (ancestral Ne), T (divergence time in generations) before OutFile.
SweepFinder2 -lrb 2000 region.freq genome.spect region.rec bvalue.map "$N1" "$N2" "$T" out.clr  # G=2000 grid; -lrb deconfounds BGS

Polarization and Polygenic Selection

Derived-allele tests (Fay & Wu H, Zeng E, unfolded SFS, the iHS sign) are POLARIZATION traps: one mispolarized CpG site masquerades as a high-frequency-derived variant and fakes a strongly negative H. Polarize with a probabilistic substitution model and two or more outgroups (Hernandez 2007), never a single chimp allele. Conservation scores (phyloP, GERP, phastCons) measure CONSTRAINT (purifying selection), the opposite sign of a recent positive sweep, and must not be read as "under selection." Polygenic Qx (Berg & Coop 2014) inherits every stratification bias in the GWAS effect sizes it uses: the celebrated European height-selection signal was shown to be a stratification artifact (Sohail 2019; Berg 2019), so use within-family/sibling effect estimates and treat standard-GWAS-beta Qx as provisional.

Per-Method Failure Modes

Absolute Tajima's D cutoff

Trigger: flagging windows on |D|>2 as selection. Mechanism: growth/bottleneck/structure produce the same sign as a sweep/balancing selection (non-identifiable at one locus). Symptom: genome-wide false positives that are pure demography. Fix: use the genome-wide empirical distribution and intersect orthogonal statistics; report D relative to the genomic background.

Unstandardized or mis-standardized iHS

Trigger: reporting raw iHS, or applying a genome-wide z-score to iHS. Mechanism: raw iHH ratio depends strongly on derived-allele frequency. Symptom: uninterpretable scores, false tails at particular frequencies. Fix: standardize_by_allele_count (DAF bins) for iHS/nSL; reserve standardize (genome-wide) for XP-EHH/iHH12.