# Population Structure Qc

> Detect and correct for population stratification and cryptic relatedness in genotype data using PCA, kinship/IBD estimation, and genomic inflation factor (lambda) diagnostics before running a GWAS. Use when doing ancestry PCA, checking sample relatedness, computing genomic inflation, or QC'ing genotype data for stratification confounding.

- Skill: `pavel-kravchenko/population-structure-qc` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/population-structure-qc`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/population-structure-qc/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/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.

```python
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.

```python
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.

```python
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