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