Selection Statistics

Detect natural selection signatures using diversity statistics and extended haplotype homozygosity.

Fst - Population Differentiation

scikit-allel

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)

num, den = allel.hudson_fst(ac_subpops['pop1'], ac_subpops['pop2'])
fst_per_snp = num / den
print(f'Mean Fst: {np.nanmean(fst_per_snp):.4f}')

Windowed Fst

fst_windowed, windows, n_snps = allel.windowed_hudson_fst(
    pos, ac_subpops['pop1'], ac_subpops['pop2'],
    size=100000, step=50000)

import matplotlib.pyplot as plt
plt.figure(figsize=(14, 4))
plt.plot(windows[:, 0], fst_windowed)
plt.xlabel('Position')
plt.ylabel('Fst')
plt.savefig('fst_windows.png')

vcftools

# Calculate Fst between populations
vcftools --vcf data.vcf --weir-fst-pop pop1.txt --weir-fst-pop pop2.txt --out fst_result

# With window
vcftools --vcf data.vcf --weir-fst-pop pop1.txt --weir-fst-pop pop2.txt \
         --fst-window-size 100000 --fst-window-step 50000 --out fst_windowed

Tajima's D - Departures from Neutrality

scikit-allel

import allel
import numpy as np

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

D, windows, counts = allel.windowed_tajima_d(pos, ac, size=100000, step=50000)

plt.figure(figsize=(14, 4))
plt.plot(windows[:, 0], D)
plt.axhline(y=0, color='r', linestyle='--')
plt.xlabel('Position')
plt.ylabel("Tajima's D")
plt.savefig('tajima_d.png')

Interpretation

D Value	Interpretation
D < -2	Recent selective sweep or population expansion
D ≈ 0	Neutral evolution
D > 2	Balancing selection or population bottleneck

vcftools

vcftools --vcf data.vcf --TajimaD 100000 --out tajima
# Output: tajima.Tajima.D (CHROM, BIN_START, N_SNPS, TajimaD)

iHS - Integrated Haplotype Score

Detects ongoing selective sweeps.

import allel
import numpy as np

callset = allel.read_vcf('data.vcf.gz')
gt = allel.GenotypeArray(callset['calldata/GT'])
pos = callset['variants/POS']
h = gt.to_haplotypes()
ac = h.count_alleles()
flt = (ac[:, 0] > 1) & (ac[:, 1] > 1)
h_flt = h.compress(flt, axis=0)
pos_flt = pos[flt]
ac_flt = ac.compress(flt, axis=0)

ihs = allel.ihs(h_flt, pos_flt, include_edges=True)
ihs_std = allel.standardize_by_allele_count(ihs, ac_flt[:, 1])

significant_ihs = np.abs(ihs_std[0]) > 2
print(f'Significant iHS hits: {significant_ihs.sum()}')

Plot iHS

import matplotlib.pyplot as plt

plt.figure(figsize=(14, 4))
plt.scatter(pos_flt, ihs_std[0], s=1)
plt.axhline(y=2, color='r', linestyle='--')
plt.axhline(y=-2, color='r', linestyle='--')
plt.xlabel('Position')
plt.ylabel('Standardized iHS')
plt.savefig('ihs.png')

XP-EHH - Cross-Population Extended Haplotype Homozygosity

Detects completed sweeps by comparing populations.

import allel
import numpy as np

h = gt.to_haplotypes()
h_pop1 = h.take(pop1_hap_idx, axis=1)
h_pop2 = h.take(pop2_hap_idx, axis=1)

xpehh = allel.xpehh(h_pop1, h_pop2, pos, include_edges=True)

significant = np.abs(xpehh) > 2
print(f'Significant XP-EHH hits: {significant.sum()}')

NSL - Number of Segregating Sites by Length

Alternative to iHS, less sensitive to recombination rate variation.

nsl = allel.nsl(h_flt)
nsl_std = allel.standardize_by_allele_count(nsl, ac_flt[:, 1])

Garud's H Statistics

Detect soft sweeps.

h1, h12, h123, h2_h1 = allel.garud_h(h)

h12_windowed = allel.moving_garud_h(h, size=100)

Composite Selection Score

Combine multiple statistics:

import numpy as np
from scipy import stats

def composite_score(fst, tajD, ihs_abs):
    fst_rank = stats.rankdata(fst) / len(fst)
    tajD_rank = stats.rankdata(-tajD) / len(tajD)  # Low Tajima's D
    ihs_rank = stats.rankdata(ihs_abs) / len(ihs_abs)
    return (fst_rank + tajD_rank + ihs_rank) / 3

css = composite_score(fst_per_snp, tajD_values, np.abs(ihs_values))

Complete Selection Scan

import allel
import numpy as np
import matplotlib.pyplot as plt

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

flt = ac.is_segregating() & (ac.max_allele() == 1)
gt = gt.compress(flt, axis=0)
pos = pos[flt]
ac = ac.compress(flt, axis=0)

window_size = 100000
window_step = 50000

tajD, tajD_windows, _ = allel.windowed_tajima_d(pos, ac, size=window_size, step=window_step)

pi, pi_windows, _, _ = allel.windowed_diversity(pos, ac, size=window_size, step=window_step)

fig, axes = plt.subplots(2, 1, figsize=(14, 8), sharex=True)

axes[0].plot(tajD_windows[:, 0], tajD)
axes[0].axhline(0, color='r', linestyle='--')
axes[0].set_ylabel("Tajima's D")

axes[1].plot(pi_windows[:, 0], pi)
axes[1].set_ylabel('Pi')
axes[1].set_xlabel('Position')

plt.tight_layout()
plt.savefig('selection_scan.png', dpi=150)

Related Skills

• scikit-allel-analysis - Data loading and basic statistics
• population-structure - Population assignment for Fst
• linkage-disequilibrium - EHH depends on LD patterns