Bio Population Genetics Selection Statistics

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

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

Selection Statistics - Usage Guide

Overview

Selection statistics detect signatures of natural selection in genomic data. Different methods detect different selection types and timescales, from recent sweeps (iHS) to ancient differentiation (Fst).

Prerequisites

pip install scikit-allel
conda install -c bioconda vcftools

Quick Start

Tell your AI agent what you want to do:

• "Calculate Fst between my two populations"
• "Scan for selection signatures using Tajima's D"
• "Compute iHS to detect ongoing selective sweeps"
• "Find regions under balancing selection"
• "Compare selection pressures between populations"

Example Prompts

Diversity Statistics

"Calculate Tajima's D in 50kb windows across the genome"

"Compute nucleotide diversity (pi) for each population"

"Find regions with unusually low diversity suggesting sweeps"

Population Differentiation

"Calculate Fst between European and African samples"

"Find highly differentiated SNPs between cases and controls"

"Generate a Manhattan plot of Fst values"

Haplotype-Based Tests

"Compute iHS scores to detect ongoing selection"

"Run XP-EHH between my populations to find completed sweeps"

"Identify haplotypes under positive selection"

Multi-Statistic Analysis

"Scan for selection using Tajima's D, Fst, and iHS together"

"Find regions significant in multiple selection tests"

"Compare selection signatures across chromosomes"

What the Agent Will Do

1. Assess data format and phase status 2. Calculate requested statistics genome-wide or in windows 3. Standardize/normalize scores where appropriate 4. Identify outlier regions exceeding thresholds 5. Generate visualizations (Manhattan plots, histograms) 6. Report candidate regions with coordinates

Tips

• Haplotype-based tests (iHS, XP-EHH) require phased data
• Demographic history can mimic selection signals
• Use multiple statistics to reduce false positives
• Always adjust for recombination rate variation
• Empirical outlier cutoffs (top 1%) are often more reliable than p-values

Selection Signatures Reference

Positive Selection (Hard Sweep)

Signs:

• Low Tajima's D (< -2)
• High |iHS| (> 2)
• High Fst between populations
• Reduced diversity (Pi)

Balancing Selection

Signs:

• High Tajima's D (> 2)
• Elevated heterozygosity
• Old alleles maintained

Recent Selection

Use haplotype-based methods:

• iHS for ongoing sweeps
• XP-EHH for completed sweeps

Ancient Selection

Use diversity-based methods:

• Fst for differentiation
• dN/dS for coding regions

Interpretation Caveats

• Demographic history mimics selection
• Recombination rate affects EHH statistics
• Multiple testing correction needed
• Functional validation recommended

Resources

• Selection Tutorial
• vcftools Manual