Population Structure Qc

Population Structure QC

When to Use

• Checking a genotype matrix for hidden ancestry structure before running association tests
• Computing principal components (PCs) to use as covariates in a GWAS model
• Detecting cryptic relatedness (unreported relatives/duplicates) via kinship or IBD estimation
• Diagnosing test-statistic inflation (genomic inflation factor λ) after a GWAS to decide whether stratification correction is needed
• Simulating population structure to validate that a QC/correction pipeline actually removes confounding

Version Compatibility

• Python ≥ 3.10, NumPy ≥ 1.26, scikit-learn ≥ 1.4, SciPy ≥ 1.11, pandas ≥ 2.0
• Concepts apply identically to PLINK2 --pca / --king-cuteoff output and to in-memory genotype matrices

Prerequisites

• pip install numpy scipy scikit-learn pandas matplotlib
• Genotype matrix encoded as allele dosage (0/1/2 copies of the minor allele), samples × SNPs
• Prior skill: bio-population-genetics-plink-basics (genotype QC filters: call rate, HWE, MAF) — run that QC before stratification analysis

Goal: Detect population stratification in a genotype matrix and derive PCs to use as GWAS covariates.

Approach: Standardize the QC'd genotype matrix, run PCA, and inspect the top PCs for cluster structure (visually and via a simple silhouette/variance check) before feeding them into the association model.

import numpy as np
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler


def pca_stratification(G, n_components=10):
    """Compute PCs of a QC'd genotype matrix for use as stratification covariates.

    G: (n_samples, n_snps) allele-dosage matrix (0/1/2), already filtered
       for call rate / HWE / MAF.
    Returns (pcs, explained_variance_ratio).
    """
    G_std = StandardScaler().fit_transform(G)
    pca = PCA(n_components=n_components)
    pcs = pca.fit_transform(G_std)
    return pcs, pca.explained_variance_ratio_


def simulate_structured_genotypes(n_samples=1000, n_snps=5000, fst=0.05, seed=42):
    """Simulate two subpopulations with allele-frequency divergence ~Fst.

    Useful for validating that a stratification-correction pipeline works:
    run GWAS with and without PCs and compare genomic inflation (lambda).
    """
    rng = np.random.default_rng(seed)
    pop = rng.integers(0, 2, n_samples)          # 0 = pop A, 1 = pop B
    maf_root = rng.uniform(0.05, 0.5, n_snps)
    # Balding-Nichols-style divergence: shift MAF per population by ~fst
    delta = rng.normal(0, fst, n_snps)
    maf_a = np.clip(maf_root + delta, 0.01, 0.99)
    maf_b = np.clip(maf_root - delta, 0.01, 0.99)
    G = np.zeros((n_samples, n_snps), dtype=int)
    for j in range(n_snps):
        p = np.where(pop == 0, maf_a[j], maf_b[j])
        G[:, j] = rng.binomial(2, p)
    return G, pop

Goal: Quantify whether GWAS test statistics are inflated by stratification or relatedness.

Approach: Compute the genomic inflation factor λ from the observed vs. expected χ² distribution of association p-values; λ ≈ 1.0 is clean, λ > 1.05–1.1 signals residual confounding that more PCs (or a mixed model) should absorb.

import numpy as np
from scipy import stats


def genomic_inflation_factor(pvals):
    """Compute lambda_GC = median(observed chi2) / median(expected chi2 under null, df=1).

    pvals: array of GWAS p-values (one per SNP tested).
    lambda ~ 1.0 => no inflation; > 1.05-1.1 => likely stratification/relatedness.
    """
    pvals = np.clip(np.asarray(pvals), 1e-300, 1.0)
    chi2_obs = stats.chi2.ppf(1 - pvals, df=1)
    lam = np.median(chi2_obs) / stats.chi2.ppf(0.5, df=1)
    return lam


def qq_plot_data(pvals):
    """Return (expected -log10 p, observed -log10 p) sorted for a QQ plot."""
    pvals = np.clip(np.sort(np.asarray(pvals)), 1e-300, 1.0)
    n = len(pvals)
    expected = (np.arange(1, n + 1)) / (n + 1)
    return -np.log10(expected), -np.log10(pvals)

Goal: Flag cryptic relatedness (duplicates, unreported siblings/cousins) that biases both PCA and association tests.

Approach: Estimate pairwise kinship from genotype dosage (KING-robust style estimator) and flag pairs above standard relatedness thresholds; remove one member of each related pair before PCA/GWAS.

import numpy as np


def kinship_king_robust(G):
    """Pairwise KING-robust kinship estimator from allele-dosage genotypes.

    G: (n_samples, n_snps) dosage in {0,1,2}. Returns (n_samples, n_samples)
    kinship matrix. Thresholds (approx.): >0.354 duplicate/MZ twin,
    >0.177 1st-degree, >0.0884 2nd-degree, >0.0442 3rd-degree.
    """
    n = G.shape[0]
    het = (G == 1)
    kinship = np.zeros((n, n))
    for i in range(n):
        for k in range(i + 1, n):
            n_het_i = het[i].sum()
            n_het_k = het[k].sum()
            ibs0 = np.sum(np.abs(G[i] - G[k]) == 2)   # opposite homozygotes
            het_shared = np.sum(het[i] & het[k])
            denom = min(n_het_i, n_het_k)
            if denom == 0:
                continue
            phi = 0.5 + (het_shared - 2 * ibs0) / (2 * denom) - (n_het_i + n_het_k) / (4 * denom)
            kinship[i, k] = kinship[k, i] = phi
    return kinship


def flag_related_pairs(kinship, threshold=0.0884):
    """Return (i, j) index pairs with kinship above threshold (2nd-degree by default)."""
    n = kinship.shape[0]
    iu = np.triu_indices(n, k=1)
    mask = kinship[iu] > threshold
    return list(zip(iu[0][mask], iu[1][mask]))

Pitfalls

• PCA run before relatedness filtering: duplicates/close relatives distort PC1–PC2; drop related samples (flag_related_pairs) before computing stratification PCs.
• Too few PCs: using only 2–3 PCs under-corrects fine-scale structure in admixed cohorts; 10 PCs is the common default, but check the scree plot / elbow.
• Too many PCs: including PCs that capture batch effects or family structure (not ancestry) can absorb real signal — validate each PC against known ancestry labels when available.
• λ misinterpretation: λ > 1.1 is not proof of stratification alone — polygenicity also inflates λ; use LD Score Regression to separate polygenic signal from confounding if λ stays high after PC correction.
• HWE filtering on cases+controls: run Hardy-Weinberg checks on controls only — true causal SNPs deviate from HWE in cases, and mixing populations inflates apparent HWE violations.
• KING-robust on small SNP sets: kinship estimates are noisy below ~10k independent SNPs; LD-prune SNPs first for a stable kinship matrix.

See Also

• bio-population-genetics-plink-basics — genotype QC (call rate, HWE, MAF) upstream of stratification analysis
• bio-population-genetics-population-structure — full ADMIXTURE/STRUCTURE-style ancestry inference
• bio-population-genetics-association-testing — per-SNP regression using PCs as covariates
• bio-workflows-gwas-pipeline — end-to-end FASTQ/VCF → GWAS pipeline including this QC step