Bio Population Genetics Linkage Disequilibrium

Version Compatibility

Reference examples tested with: PLINK 1.9 (1.90b7+), PLINK 2.0 (alpha 6+), scikit-allel 1.3+, numpy 1.26+.

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: PLINK 2.0 has NO bare --r2 (it was split into --r2-phased, the EM haplotype-frequency estimator, and --r2-unphased, the composite dosage estimator); PLINK 1.9 still uses bare --r2. PLINK 2.0 --indep-pairwise <window>['kb'] [step] <r2> requires the step to be 1 when the window is given in kb. allel.rogers_huff_r returns a CONDENSED upper-triangle array, not a square matrix - scipy.spatial.distance.squareform it before 2D indexing. The single source of truth for versions is this block, not headings.

Linkage Disequilibrium

"Measure LD between my variants" -> Quantify how predictably one variant's alleles co-occur with another's, choosing r2 or D' by the question being asked.

• CLI: plink2 --r2-unphased (composite dosage r2, no phasing) or --r2-phased (EM haplotype r2)
• Python: allel.rogers_huff_r(gn) ** 2 (composite r2, scikit-allel)

"Thin correlated variants before PCA or structure" -> Remove variants so no remaining pair is in high LD, leaving a near-independent marker set.

• CLI: plink2 --indep-pairwise 50 5 0.1 (genotype-blind, phenotype-agnostic)
• Python: allel.locate_unlinked(gn, size=50, step=5, threshold=0.1) (scikit-allel)

"Reduce my GWAS hits to one signal per locus" -> Group correlated associations around the most significant SNP per region.

• CLI: plink --clump sumstats --clump-p1 5e-8 --clump-r2 0.1 --clump-kb 250

Scope: LD measures, pruning, clumping, haplotype blocks, and LD decay. QC/conversion route to plink-basics; PCA/ADMIXTURE to population-structure; --glm GWAS to association-testing; resolving independent causal variants to causal-genomics/fine-mapping; phased haplotypes to phasing-imputation/haplotype-phasing.

The Single Most Important Insight -- D' and r2 answer different questions

1. D' measures whether recombination has been observed between two loci (a historical/structural question); r2 measures how well one variant predicts the other (a statistical/predictive question), so they are not interchangeable.
2. r2 owns tagging, imputation, pruning, and GWAS power because r2 = chi2/N exactly (N = number of haplotypes = 2 x sampled individuals, for the 1-df allelic 2x2 table; the composite estimator targets this gametic r2 under random mating): a proxy multiplies the association test's non-centrality by r2, so 1/r2 times the sample size is needed to recover the lost power (Pritchard & Przeworski 2001).
3. D'=1 routinely coexists with r2 ~ 0.05 because allele-frequency asymmetry caps r2: r2/r2max = D'^2 over most of frequency space (VanLiere & Rosenberg 2008), which is exactly why a common array SNP cannot tag a rare causal variant despite "complete LD".
4. Every operation that consumes an LD matrix (clumping, fine-mapping, LD-score regression) silently assumes that matrix matches the study sample's LD; a wrong-ancestry or underpowered reference does not error, it produces confidently wrong output.

Tool Taxonomy

Decision Tree by Scenario

Pairwise LD and the Phased/Unphased Choice

Goal: Compute pairwise r2 (and D' when recombination is the question) without an EM artifact under structure or missingness.

Approach: Default to the composite dosage estimator that needs no phasing; reach for EM haplotype r2 only when the sample is large and HWE-consistent, and use D' solely for block/recombination work.

# Composite dosage r2 (Rogers-Huff): robust default, no HWE assumption.
plink2 --bfile data --r2-unphased --ld-window-kb 1000 --ld-window-r2 0.2 --out ld_unphased

# EM haplotype-frequency r2: only when phase is trustworthy / sample is large and HWE-consistent.
plink2 --bfile data --r2-phased --ld-window-kb 1000 --out ld_phased

# D' for a single pair (recombination / block question), with observed haplotype frequencies.
plink2 --bfile data --ld rs123 rs456 --out pair

# PLINK 1.9 still has the bare --r2; add dprime to also report D'.
plink --bfile data --r2 dprime --ld-window-kb 500 --out ld_dprime

import allel
from scipy.spatial.distance import squareform

callset = allel.read_vcf('data.vcf.gz')
gn = allel.GenotypeArray(callset['calldata/GT']).to_n_alt()

# rogers_huff_r returns a CONDENSED upper-triangle vector; square it for r2, squareform for a matrix.
r2_matrix = squareform(allel.rogers_huff_r(gn[:200]) ** 2)

LD Pruning for Structure (genotype-blind)

Goal: Produce a near-independent marker set so PCA/ADMIXTURE/GRM are not dominated by a handful of LD blocks.

Approach: Exclude long-range-LD regions by coordinate FIRST (their internal r2 is high and real, so a threshold cannot remove them sensibly), then slide an r2 window over the survivors.

# 1. Drop long-range-LD regions by position (MHC chr6:25-35 Mb, 8p23.1, 17q21.31, LCT/2q21; coordinates are build-specific).
plink2 --bfile data --exclude range longrange_ld.txt --make-bed --out data_noLR

# 2. Windowed r2 prune (50 here is a variant count, so step 5 is fine; step MUST be 1 only with a kb window).
plink2 --bfile data_noLR --indep-pairwise 50 5 0.1 --out prune     # 50-variant window, step 5, r2 0.1
plink2 --bfile data_noLR --extract prune.prune.in --make-bed --out data_pruned

import allel

callset = allel.read_vcf('data.vcf.gz')
gn = allel.GenotypeArray(callset['calldata/GT']).to_n_alt()

loc_unlinked = allel.locate_unlinked(gn, size=50, step=5, threshold=0.1)  # boolean keep-mask
gn_pruned = gn.compress(loc_unlinked, axis=0)

Clumping GWAS Summary Statistics (p-value-aware)

Goal: Collapse a region of correlated associations to one lead SNP per independent locus.

Approach: Index on genome-wide-significant SNPs and absorb nearby LD partners, overriding PLINK's permissive defaults; use an ancestry-matched LD reference, ideally the study sample itself.

# Defaults (p1=1e-4, p2=1e-2, r2=0.5, kb=250) are neither genome-wide nor strict - set them explicitly.
plink --bfile ld_reference \
    --clump gwas_sumstats.txt \
    --clump-p1 5e-8 --clump-p2 1e-5 --clump-r2 0.1 --clump-kb 250 \
    --out clumped
# clumped.clumped lists one index SNP per locus. This is NOT conditional analysis or fine-mapping.

Haplotype Blocks and LD Decay

Goal: Map regions of little observed recombination, or characterize how LD decays with distance.

Approach: Use the Gabriel D'-CI block method for boundaries; for decay, bin composite r2 by physical distance after stratifying by population so admixture LD does not flatten the curve.

plink --bfile data --blocks no-pheno-req --out blocks   # blocks.blocks, blocks.blocks.det (Gabriel CI)

import allel, numpy as np

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

n = min(1000, gn.shape[0])
r2, dist = [], []
for i in range(n):
    for j in range(i + 1, min(i + 100, n)):
        r2.append(allel.rogers_huff_r(gn[[i, j]])[0] ** 2)  # condensed length-1 -> index [0]
        dist.append(pos[j] - pos[i])
r2, dist = np.array(r2), np.array(dist)

edges = np.arange(0, 100001, 1000)
decay = [np.mean(r2[(dist >= edges[k]) & (dist < edges[k + 1])]) if ((dist >= edges[k]) & (dist < edges[k + 1])).any() else np.nan for k in range(len(edges) - 1)]

Per-Method Failure Modes

Long-range-LD regions survive pruning

Trigger: --indep-pairwise run without excluding MHC/inversions by coordinate. Mechanism: their internal r2 is high and genuine, so a window prune keeps a dense cluster. Symptom: the top PCs capture the MHC (chr6:25-35Mb) or 17q21.31 inversion, not ancestry. Fix: --exclude range the long-range-LD list (MHC, 8p23.1, 17q21.31, LCT; Price 2008) before pruning, not a tighter r2.

D' read as "high LD" for tagging

Trigger: selecting tag SNPs or judging proxy adequacy from D'. Mechanism: D'=1 only says no recombinant was observed; frequency asymmetry caps r2 far below 1. Symptom: a "perfectly linked" common SNP retains almost no association power for a rare causal variant. Fix: use r2 for tagging/power and require r2 >= 0.8; reserve D' for blocks.

D' inflated at rare alleles / small N

Trigger: interpreting D' or --blocks where minor-allele count is below ~10-20. Mechanism: the fourth haplotype is unobserved by chance, pushing the D' MLE to ~1 even for independent loci. Symptom: spurious "perfect LD" blocks that vanish with more samples. Fix: do not interpret D' at low MAC; report r2, which is noisy but not systematically inflated.