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.

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.

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).

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))

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).

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: 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

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

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.