Bio Population Genetics Scikit Allel Analysis

Version Compatibility

Reference examples tested with: scikit-allel 1.3.13+, numpy 1.26+, zarr 2.18+.

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

• Python: pip show <package> then help(module.function) to check signatures

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: scikit-allel is in MAINTENANCE mode (latest line v1.3.x, e.g. v1.3.13 Sep 2024); the README names sgkit (xarray+dask) as the successor but states it is "not yet at feature parity", so scikit-allel remains the pragmatic choice for established stat workflows. average_patterson_f3/average_patterson_d are the current names; pre-2020 code used blockwise_patterson_* (gone). to_n_alt() default is fill=0, not fill=-1. read_vcf is eager and loads the whole file into RAM. The single source of truth for versions is this block, not headings.

scikit-allel Analysis

"Analyze population genetics in Python" -> Read a VCF into array structures, then compute frequency-, diversity-, differentiation-, and haplotype-based statistics with correct denominators.

• Python: allel.GenotypeArray, allel.AlleleCountsArray, allel.windowed_diversity(..., is_accessible=), allel.average_hudson_fst, allel.pca

Scope: in-memory and dask/zarr scikit-allel analysis - the array-API mechanics of the data model, diversity/SFS, FST and f/D admixture statistics, LD pruning, PCA, and selection-scan computation. PLINK-format QC routes to plink-basics; PCA/ADMIXTURE via CLI tools to population-structure; selection-scan DESIGN (standardization, outlier calling, demographic confounding) to selection-statistics; phased input to phasing-imputation/haplotype-phasing; VCF generation to variant-calling/vcf-basics.

The Single Most Important Insight -- one denominator bug, two faces

1. Almost every statistic scikit-allel computes is a RATIO or a DENSITY, and the two ways the denominator is silently wrong both return a finite, plausible, unflagged number.
2. Face (a) - the per-base denominator: pi, Watterson's theta, and Dxy divide a segregating-site sum by the number of bases; omit is_accessible= and they divide by total span (stop-start+1) instead of callable bp, deflating values 2-5x AND distorting the genome-wide landscape because callability varies per window.
3. Face (b) - FST is a ratio of variance components: the correct multi-locus estimate is sum(a)/(sum(a)+sum(b)+sum(c)), NOT mean(per_snp_fst) (which is rare-variant-dominated and biased, Bhatia 2013); scikit-allel returns the (a,b,c) components from weir_cockerham_fst and (num,den) from hudson_fst/patterson_fst precisely to force ratio-of-sums.
4. Both bugs run silently with no exception or warning, so correctness is a design decision the analyst makes, not something the library enforces.

Tool Taxonomy -- the data model and the scaling path

Decision Tree by Scenario

Reading Genotypes and Counting Alleles

Goal: Load a VCF region into a GenotypeArray and derive the allele-count currency, missing-aware.

Approach: Read only the needed fields (read_vcf is eager), wrap GT, count alleles per site (missing ignored), and get per-population counts in one pass with count_alleles_subpops.

import allel
import numpy as np

callset = allel.read_vcf('data.vcf.gz', fields=['samples', 'calldata/GT', 'variants/POS', 'variants/CHROM'], region='2L:1-5000000')
gt = allel.GenotypeArray(callset['calldata/GT'])   # (n_variants, n_samples, 2); -1 = missing
pos = callset['variants/POS']

ac = gt.count_alleles()                            # ignores -1, so per-site allele number varies
subpops = {'pop1': [0, 1, 2, 3, 4], 'pop2': [5, 6, 7, 8, 9]}
ac_subpops = gt.count_alleles_subpops(subpops)     # one pass, consistent variant axis
ac1, ac2 = ac_subpops['pop1'], ac_subpops['pop2']

Per-base Diversity with the Accessibility Mask

Goal: Compute pi and Watterson's theta as honest per-base quantities, not span-deflated ones.

Approach: Pass a boolean callability mask (one entry per base, from coverage/mappability, NOT from variant positions) as is_accessible; inspect the returned n_bases per window to confirm the denominator.

# is_accessible: bool array over genomic positions (a callable-loci mask), NOT the VCF variant sites.
pi = allel.sequence_diversity(pos, ac, is_accessible=is_accessible)
theta_w = allel.watterson_theta(pos, ac, is_accessible=is_accessible)

# windowed_diversity returns 4 values; n_bases is the accessible-bp denominator PER window.
pi_w, windows, n_bases, counts = allel.windowed_diversity(pos, ac, size=100000, is_accessible=is_accessible)
# windowed_tajima_d returns 3 values (no n_bases): Tajima's D is dimensionless, no is_accessible.
D, td_windows, td_counts = allel.windowed_tajima_d(pos, ac, size=100000)

FST as a Ratio of Sums

Goal: Get a genome-wide FST point estimate with a jackknife SE, and a per-window landscape, without the mean-of-ratios bias.

Approach: Let average_hudson_fst do the ratio-of-sums plus block-jackknife; if hand-aggregating, sum the components THEN divide; size blocks (blen) to exceed the LD scale.

# Genome-wide estimate + standard error (ratio-of-sums + delete-one-block jackknife):
fst, se, vb, vj = allel.average_hudson_fst(ac1, ac2, blen=2000)   # blen must exceed the LD decay length

# Hand-aggregating Hudson correctly (NEVER mean of per-SNP fst):
num, den = allel.hudson_fst(ac1, ac2)
fst_manual = np.sum(num) / np.sum(den)

# Weir-Cockerham returns per-allele components (a, b, c); aggregate over BOTH axes:
a, b, c = allel.weir_cockerham_fst(gt, subpops=[[0, 1, 2, 3, 4], [5, 6, 7, 8, 9]])
fst_wc = np.sum(a) / (np.sum(a) + np.sum(b) + np.sum(c))

# Landscape: per-window FST is already ratio-aggregated within each window.
fst_windows = allel.moving_hudson_fst(ac1, ac2, size=1000)

Admixture: f3 and D (block-jackknifed)

Goal: Formally test whether a population is admixed (f3) or whether gene flow violates a tree (D / ABBA-BABA).

Approach: Call the average_* form for the jackknife z-score; put the test population FIRST in f3; treat a significantly negative f3 (z < ~-3) as admixture and |z| > ~3 for D as treeness violation.

# f3(C; A, B): TEST population C is the FIRST argument. Returns (f3, se, z, vb, vj).
f3, se3, z3, vb3, vj3 = allel.average_patterson_f3(acc, aca, acb, blen=2000)
# Significantly negative f3 (z3 < ~-3) => C is admixed between A and B.

# D-statistic (ABBA-BABA): returns (d, se, z, vb, vj). |z| > ~3 flags gene flow.
d, sed, zd, vbd, vjd = allel.average_patterson_d(aca, acb, acc, acd, blen=2000)

LD Pruning and PCA

Goal: Project samples onto ancestry axes that reflect drift, not LD blocks or inversions.

Approach: Convert to a missing-free 012 matrix, LD-prune iteratively with locate_unlinked, mask known inversions by position, then run Patterson-scaled PCA (randomized at scale).

gn = gt.to_n_alt(fill=-1)                          # default fill=0 imputes missing to REF; use -1 then handle
gn = np.where(gn < 0, 0, gn)                        # impute-to-reference is a deliberate choice here

# Iterate LD pruning ~3 rounds; one pass leaves residual LD. Returns a KEEP mask (True = unlinked).
for _ in range(3):
    keep = allel.locate_unlinked(gn, size=100, step=20, threshold=0.1)
    gn = gn[keep]

coords, model = allel.randomized_pca(gn, n_components=10, scaler='patterson', random_state=0)
explained = model.explained_variance_ratio_        # scree; coords is (n_samples, n_components)

Selection Scans on Phased Haplotypes

Goal: Score the genome for recent selection from haplotype structure.

Approach: Reshape phased genotypes to a HaplotypeArray, compute the raw scan, then standardize - iHS/nSL binned by derived-allele frequency (standardize_by_allele_count), XP-EHH genome-wide (standardize); the raw scores are not directly interpretable.

h = gt.to_haplotypes()                              # VALID only if data are PHASED
ihs_raw = allel.ihs(h, pos, min_maf=0.05)           # unstandardized
ihs_std, bins = allel.standardize_by_allele_count(ihs_raw, ac[:, 1])   # bin by DERIVED count; ac[:,1] is derived only if REF is ancestral (polarize first); |z| > 2 flags candidates
h1, h12, h123, h2_h1 = allel.garud_h(h)             # soft-vs-hard-sweep haplotype-homozygosity stats

Per-Function Failure Modes

sequence_diversity / watterson_theta without is_accessible

Trigger: calling per-base diversity with no callability mask. Mechanism: divides the numerator by stop-start+1 (total span) instead of callable bp. Symptom: pi/theta deflated 2-5x and the genome-wide landscape distorted because callability varies per window. Fix: pass is_accessible= from a coverage/mappability callable-loci mask and inspect the returned n_bases.

mean(per_snp_fst) instead of ratio-of-sums

Trigger: averaging per-SNP a/(a+b+c) for the genome-wide FST. Mechanism: mean-of-ratios is dominated by low-frequency SNPs with tiny noisy denominators. Symptom: a biased FST that differs from published estimates of the same comparison (Bhatia 2013). Fix: sum(a)/(sum(a)+sum(b)+sum(c)), or average_hudson_fst/average_weir_cockerham_fst which do it plus a jackknife SE.

to_n_alt default fill=0

Trigger: gt.to_n_alt() with no fill. Mechanism: missing calls become 0 alt alleles = homozygous reference. Symptom: PCA/LD silently biased toward the reference allele. Fix: to_n_alt(fill=-1) then impute deliberately, or pre-filter for high call rate; state the imputation choice.

sfs() fed alt counts and run unfolded

Trigger: allel.sfs(ac[:, 1]) without a confident ancestral allele. Mechanism: sfs() is unfolded and treats the ALT count as the DERIVED count; ALT != DERIVED. Symptom: mis-polarized spectrum biasing demographic/DFE inference. Fix: use sfs_folded(ac) when polarization is uncertain; use sfs(dac) only with a verified ancestral allele.