# Bio Applied Gwas

> Run GWAS: SNP QC (MAF/HWE), PCA ancestry covariates, per-SNP logistic/linear regression, Manhattan/QQ plots, genomic inflation, LD clumping. Use when doing GWAS, SNP-array association, or stratification correction.

- Skill: `pavel-kravchenko/bio-applied-gwas` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-applied-gwas`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-applied-gwas/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/bio-applied-gwas

---


# Genome-Wide Association Studies (GWAS)

## When to Use

- Testing hundreds of thousands to millions of SNPs for association with a binary (case/control) or quantitative trait.
- Correcting for population stratification / cryptic relatedness before interpreting association signal.
- Producing Manhattan and QQ plots, and calling genome-wide significant hits (p < 5×10⁻⁸).
- Reducing a list of significant SNPs to independent signals via LD clumping, or prioritizing loci with GWAS Catalog lookups.
- Sanity-checking a PLINK/REGENIE association run (genomic inflation factor λ, expected vs. observed p-value distribution).

## Version Compatibility

- Python ≥ 3.10, NumPy ≥ 1.26, pandas ≥ 2.0, SciPy ≥ 1.11, scikit-learn ≥ 1.3, Matplotlib ≥ 3.8.
- Production-scale GWAS: PLINK2 ≥ 2.00a (`--glm`, `--clump`), REGENIE ≥ 3.x, or SAIGE for mixed models. The from-scratch code below is for understanding the statistics / prototyping on small cohorts; at biobank scale (>100K samples, >1M SNPs) use PLINK2/REGENIE, not per-SNP Python loops.

## Prerequisites

- `pip install numpy pandas scipy scikit-learn matplotlib`
- Concepts: Hardy-Weinberg equilibrium, minor allele frequency (MAF), linkage disequilibrium (LD), principal component analysis for ancestry, multiple-testing correction. See `bio-applied-population-genetics` and `bio-applied-statistics-for-bioinformatics`.
- Input: a genotype dosage matrix (samples × SNPs, values in {0,1,2}) plus a phenotype vector — e.g. loaded from PLINK `.bed/.bim/.fam` via `pandas-plink` or simulated as below.

## Study Design

Sample size: n > 5,000 for common variants, n > 50,000 for small effects. Include age, sex, ancestry, and batch as covariates. **Why p < 5×10⁻⁸:** ~1M approximately independent SNPs after LD pruning → Bonferroni α = 0.05 / 10⁶. **Genomic inflation factor λ** = median(χ²_obs) / median(χ²_expected_null): λ ≈ 1.0 is ideal, λ > 1.1 suggests stratification, cryptic relatedness, or polygenicity (LDSC intercept disentangles these at scale).

**Goal:** simulate a realistic genotype matrix and case/control phenotype with known causal SNPs, so the QC/association/plotting code below can be tested end-to-end.
**Approach:** draw per-SNP allele frequencies, sample genotypes under Hardy-Weinberg proportions, then inject a logistic effect from a handful of causal SNPs.

```python
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
from scipy.stats import chi2 as chi2_dist
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression

rng = np.random.default_rng(42)
N_SAMPLES, N_SNPS = 1000, 3000

mafs = rng.uniform(0.05, 0.45, N_SNPS)
p = mafs
freqs = np.column_stack([(1 - p) ** 2, 2 * p * (1 - p), p ** 2])  # HWE genotype proportions
G = np.vstack([rng.choice([0, 1, 2], size=N_SAMPLES, p=freqs[j]) for j in range(N_SNPS)]).T

# Inject a logistic-model phenotype driven by 3 causal SNPs
causal_idx = [300, 1200, 2500]
betas = [0.6, -0.5, 0.4]
logit = sum(b * G[:, c] for b, c in zip(betas, causal_idx)) + rng.normal(0, 0.5, N_SAMPLES)
prob = 1 / (1 + np.exp(-logit))
phenotype = (rng.random(N_SAMPLES) < prob).astype(int)

chrom_ids = np.repeat(["chr1", "chr2", "chr3"], [1000, 1000, 1000])
positions = np.concatenate([
    rng.integers(1e6, 250e6, 1000),
    rng.integers(1e6, 240e6, 1000),
    rng.integers(1e6, 198e6, 1000),
])
snp_df = pd.DataFrame({"chrom": chrom_ids, "pos": positions, "maf": mafs})
```

## SNP QC and Ancestry Correction

**Goal:** drop low-quality/low-frequency SNPs and derive ancestry covariates so association tests aren't confounded by population structure.
**Approach:** filter on MAF, call rate, and Hardy-Weinberg p-value (computed in controls only, since cases can show real HWE deviation at causal loci); then run PCA on standardized genotypes and keep the top PCs as covariates.

```python
def qc_snps(G, snp_df, phenotype, min_maf=0.01, min_call_rate=0.95, hwe_p_thresh=1e-6):
    """Filter SNPs on MAF, call rate, and HWE (controls only).

    Returns the filtered genotype matrix, filtered SNP metadata, and the boolean mask.
    """
    maf = np.minimum(G.mean(0) / 2, 1 - G.mean(0) / 2)
    call_rate = np.ones(G.shape[1])  # replace with real missingness fraction if data has NaNs

    controls = phenotype == 0
    G_ctrl = G[controls]
    obs_hets = (G_ctrl == 1).mean(0)
    p_ctrl = G_ctrl.mean(0) / 2
    exp_hets = 2 * p_ctrl * (1 - p_ctrl)
    chi2_hwe = (obs_hets - exp_hets) ** 2 / (exp_hets + 1e-9) * controls.sum()
    hwe_p = stats.chi2.sf(chi2_hwe, df=1)

    keep = (maf >= min_maf) & (call_rate >= min_call_rate) & (hwe_p >= hwe_p_thresh)
    return G[:, keep], snp_df[keep].reset_index(drop=True), keep


G_qc, snp_qc, keep_mask = qc_snps(G, snp_df, phenotype)

G_std = StandardScaler().fit_transform(G_qc)
pca = PCA(n_components=10)
PCs = pca.fit_transform(G_std)  # top PCs = ancestry covariates
```

## Association Testing, Manhattan/QQ Plots

**Association tests:** binary trait → logistic regression; quantitative trait → linear regression. Both include top PCs as ancestry covariates.

**Goal:** test every SNP for association with the phenotype while adjusting for ancestry PCs, then visualize genome-wide signal and calibration.
**Approach:** per-SNP Wald test from a logistic-regression fit (swap in `statsmodels.OLS`/`GLM` for quantitative traits or when you need standard errors from a proper solver rather than a manual information-matrix inversion).

```python
def gwas_logistic(G, phenotype, covariates):
    """Per-SNP logistic regression (additive model) with covariate adjustment.

    Returns a p-value array, one per column of G. Uses the Wald test on the
    SNP coefficient; covariates (e.g. ancestry PCs) are included but not tested.
    """
    n = len(phenotype)
    pvals = np.ones(G.shape[1])
    for j in range(G.shape[1]):
        snp = G[:, j].astype(float)
        X = np.column_stack([np.ones(n), snp, covariates])
        try:
            lr = LogisticRegression(solver="lbfgs", max_iter=300, C=1e9,
                                     fit_intercept=False, random_state=0)
            lr.fit(X, phenotype)
            mu = lr.predict_proba(X)[:, 1]
            w = mu * (1 - mu)                       # Fisher information weights
            XtWX = (X.T * w) @ X                     # observed information matrix
            cov_beta = np.linalg.inv(XtWX)
            se = np.sqrt(max(cov_beta[1, 1], 1e-15))
            z = lr.coef_[0][1] / se
            pvals[j] = chi2_dist.sf(z ** 2, df=1)     # Wald chi^2, 1 df
        except (np.linalg.LinAlgError, ValueError):
            pvals[j] = 1.0
    return pvals


pvals = gwas_logistic(G_qc, phenotype, covariates=PCs[:, :5])
snp_qc["pval"] = pvals
snp_qc["-log10p"] = -np.log10(snp_qc["pval"].clip(1e-300))
```

```python
fig, axes = plt.subplots(2, 1, figsize=(14, 8))

# Manhattan plot
chroms = ["chr1", "chr2", "chr3"]
colors = ["#1f77b4", "#ff7f0e", "#2ca02c"]
offset, xticks, xticklabels = 0, [], []
for i, c in enumerate(chroms):
    sub = snp_qc[snp_qc["chrom"] == c]
    axes[0].scatter(sub["pos"] + offset, sub["-log10p"], c=colors[i], s=2, alpha=0.5)
    xticks.append(offset + sub["pos"].median())
    xticklabels.append(c)
    offset += sub["pos"].max() + 10e6
axes[0].axhline(-np.log10(5e-8), color="red", lw=1.2, ls="--", label="5×10⁻⁸")
axes[0].axhline(-np.log10(1e-5), color="orange", lw=0.8, ls="--", label="1×10⁻⁵ (suggestive)")
axes[0].set_xticks(xticks); axes[0].set_xticklabels(xticklabels)
axes[0].set_ylabel("-log10(p)"); axes[0].set_title("Manhattan plot"); axes[0].legend(frameon=False)

# QQ plot + genomic inflation factor
observed = np.sort(pvals)
expected = np.arange(1, len(observed) + 1) / (len(observed) + 1)
lam = np.median(stats.chi2.ppf(1 - observed, 1)) / stats.chi2.ppf(0.5, 1)
axes[1].scatter(-np.log10(expected), -np.log10(observed), s=3, alpha=0.5, color="steelblue")
axes[1].plot([0, -np.log10(expected[-1])], [0, -np.log10(expected[-1])], "r--", lw=1)
axes[1].set_title(f"QQ plot (lambda = {lam:.3f})")
plt.tight_layout(); plt.show()
```

## LD Clumping and Fine-Mapping

Significant SNPs are rarely independent — nearby LD-linked SNPs tag the same signal. Clumping keeps only the lead (most significant) SNP per LD block:
1. Sort candidate SNPs by p-value.
2. For each lead SNP, drop all remaining SNPs within ±250 kb that have r² > 0.1 with it (PLINK: `plink2 --clump --clump-p1 5e-8 --clump-r2 0.1 --clump-kb 250`).

```python
def clump(snp_qc, r2_matrix, snp_indices, p_thresh=5e-8, r2_thresh=0.1):
    """Greedy LD clumping: keep top SNP per LD block, drop correlated neighbors.

    snp_indices maps rows/cols of r2_matrix to positions in snp_qc.
    """
    candidates = snp_qc.loc[snp_indices]
    candidates = candidates[candidates["pval"] < p_thresh].sort_values("pval")
    kept, removed = [], set()
    for idx in candidates.index:
        if idx in removed:
            continue
        kept.append(idx)
        local = snp_indices.index(idx)
        correlated = [snp_indices[k] for k, r2 in enumerate(r2_matrix[local]) if r2 > r2_thresh]
        removed.update(correlated)
    return kept


chr1_idx = snp_qc.index[snp_qc["chrom"] == "chr1"].tolist()[:50]
chr1_G = G_qc[:, chr1_idx].T
r2_matrix = np.corrcoef(chr1_G) ** 2
lead_snps = clump(snp_qc, r2_matrix, chr1_idx)
```

Statistical fine-mapping (SuSiE, FINEMAP) models a locus as L causal signals and computes a 95% credible set; combine with functional annotations (ENCODE, GTEx eQTLs) to prioritize credible-set members. See `bio-causal-genomics-fine-mapping`-style tools if available, or run `susieR::susie_rss` in R on summary statistics.

```r
# R: SuSiE fine-mapping from GWAS summary statistics (z-scores) and an LD matrix
library(susieR)
fit <- susie_rss(z = snp_qc$zscore, R = r2_matrix, n = nrow(G_qc), L = 10)
credible_sets <- susie_get_cs(fit, coverage = 0.95)
print(credible_sets)
```

## GWAS Catalog Lookup

```python
import json
import urllib.request


def query_gwas_catalog(rsid):
    """Look up known associations for an rsID via the GWAS Catalog REST API.

    Returns a list of (pvalue, [traits]) tuples, or an offline placeholder on failure.
    """
    url = (f"https://www.ebi.ac.uk/gwas/rest/api/singleNucleotidePolymorphisms/"
           f"{rsid}/associations?projection=associationBySnp")
    try:
        with urllib.request.urlopen(url, timeout=5) as r:
            data = json.load(r)
        assocs = data.get("_embedded", {}).get("associations", [])
        return [(a.get("pvalue", ""), [t["trait"] for t in a.get("efoTraits", [])])
                 for a in assocs[:3]]
    except Exception:
        return [("(offline)", ["GWAS Catalog requires internet access"])]


# e.g. rs429358 = APOE epsilon-4 allele, a well-known Alzheimer's risk variant
results = query_gwas_catalog("rs429358")
```

## GWAS Checklist

| Step | Method | Key Threshold |
|---|---|---|
| SNP QC | MAF, call rate, HWE | MAF > 0.01, HWE p > 10⁻⁶ |
| Sample QC | Missingness, ancestry outliers | missingness < 5% |
| Stratification | PCA → include PCs as covariates | top 10 PCs standard |
| Association | Logistic / linear regression | adjusted for covariates |
| Multiple testing | Bonferroni / GW threshold | p < 5×10⁻⁸ |
| LD / Clumping | PLINK `--clump` | r² < 0.1, window 250 kb |
| Fine-mapping | SuSiE / FINEMAP | 95% credible set |
| Lookup | GWAS Catalog REST API | validate prior associations |

## Pitfalls

- **Uncorrected stratification**: skipping PCA covariates inflates λ and produces false-positive hits that track ancestry, not phenotype — always plot PC1 vs PC2 colored by case/control before trusting results.
- **HWE test on all samples**: computing HWE using cases pools in real disease-associated deviation at causal loci; always test HWE in controls only.
- **Per-SNP Python loops at scale**: the logistic-regression loop above is O(SNPs × samples) and fine for <10K SNPs; for genome-wide data (>500K SNPs) use PLINK2 `--glm` or REGENIE, which use much faster approximations (e.g. saddlepoint/firth corrections for rare variants and case-control imbalance).
- **Ignoring LD when counting hits**: reporting every significant SNP without clumping overstates the number of independent discoveries.
- **Multiple testing**: use the standard p < 5×10⁻⁸ genome-wide threshold for SNP-level tests, or Benjamini-Hochberg FDR for other multi-feature settings (e.g. gene-based tests).

## See Also

- `bio-applied-population-genetics` — Hardy-Weinberg, LD, Fst, and population structure fundamentals.
- `bio-applied-statistics-for-bioinformatics` — hypothesis testing and multiple-testing correction background.
- `bio-applied-variant-calling-and-snp-analysis` — producing the genotype/VCF data that feeds a GWAS.
- `bio-applied-clinical-genomics` — interpreting and reporting GWAS-derived variants clinically.