Interactive genome-wide association study — from genotyping arrays to Manhattan plots, fine-mapping, and functional follow-up of disease loci
A genome-wide association study begins with careful assembly of a case-control cohort — individuals with a disease of interest and matched, unaffected controls — followed by dense genome-wide genotyping. Modern GWAS routinely genotype hundreds of thousands to millions of people, relying on cost-effective SNP microarrays and statistical imputation to reach near-complete coverage of common genetic variation.
Case-control GWAS design compares allele frequencies between disease cases and unaffected controls:
Cohort assembly: • Cases: clinically confirmed diagnosis (ICD-10 codes, physician adjudication, biomarker thresholds) • Controls: population-representative, disease-free, ideally same ancestry and recruitment source as cases • Family-based alternative: trio design (proband + parents) controls for population stratification by construction, used in TDT (transmission disequilibrium test)
Minimizing confounding: • Batch effects: cases and controls genotyped on different plates/arrays create spurious associations — mitigate by randomizing case/control status across plates • Population stratification: allele frequency differences between subpopulations (e.g. ancestry) that correlate with disease prevalence can create false associations — corrected downstream with principal components • Winner's curse: effect sizes from discovery GWAS are systematically inflated; replication in independent cohorts is mandatory
Sample size and power: • Power to detect OR=1.1 at MAF=0.3 with p<5×10⁻⁸: needs ~50,000–100,000 cases for well-powered discovery • Power to detect OR=1.5 at MAF=0.05: achievable with ~5,000–10,000 cases • Rule of thumb: doubling sample size roughly doubles the number of discovered loci for polygenic traits
Direct genotyping arrays (Illumina Global Screening Array, Affymetrix UK Biobank Axiom) assay a curated set of ~500k–700k tag SNPs chosen to capture common variation via linkage disequilibrium. Statistical imputation then infers untyped genotypes using a dense reference haplotype panel:
Imputation pipeline: 1. Phasing: array genotypes phased into haplotypes (Eagle2, SHAPEIT4) 2. Reference panel matching: phased haplotypes matched against TOPMed (~97k diverse genomes) or 1000 Genomes Project (2,504 genomes) 3. Imputation: Minimac4 / IMPUTE5 use a hidden Markov model over reference haplotypes to infer dosage of untyped SNPs 4. Quality filter: imputation INFO/R² score >0.3–0.8 required to retain a variant
Sample-level QC exclusions: • Call rate <98% (excess missingness — DNA degradation, array failure) • Heterozygosity outliers (>3 SD from mean — contamination or inbreeding) • Sex-check mismatch (genotypic sex vs. reported sex — sample mix-up) • Cryptic relatedness (pairwise IBD >0.185 — remove one of each related pair, typically via KING) • Ancestry outliers (PCA projection far from cohort centroid)
Variant-level QC exclusions: • Minor allele frequency (MAF) <0.01 (underpowered, imputation-unreliable) • Hardy-Weinberg equilibrium p<1×10⁻⁶ in controls (genotyping artifact signature) • Call rate <98% across samples • Differential missingness between cases and controls (p<1×10⁻⁵, batch artifact)
The UK Biobank genotyped 500,000 participants on a custom Affymetrix array (~800k markers) then imputed to ~96 million variants using the Haplotype Reference Consortium panel — a single resource that has powered thousands of downstream GWAS across nearly every human disease.
With a curated, imputed genotype matrix in hand, the core statistical engine of GWAS runs a separate regression test at every one of tens of millions of SNPs, asking: does carrying the risk allele change the odds of disease, after adjusting for age, sex, and genetic ancestry? Modern GWAS pipelines (PLINK2, REGENIE, SAIGE) make this computationally tractable even for biobank-scale cohorts.
For a binary disease trait, each SNP is tested independently via logistic regression:
logit(P(disease)) = β₀ + β₁·dosage + β₂·age + β₃·sex + Σ βₖ·PCₖ
Where dosage ∈ [0,2] is the imputed expected count of the risk allele (0, 1, or 2 copies under the additive model, or a fractional dosage reflecting imputation uncertainty). The coefficient β₁ is exponentiated to give the odds ratio (OR = e^β₁): the multiplicative change in disease odds per additional risk allele copy.
For quantitative traits (height, LDL cholesterol, BMI), linear regression replaces logistic regression, with β₁ interpreted directly as the trait-unit change per allele copy.
Alternative test statistics: • Chi-square / allelic test: simpler 2×2 contingency table test, doesn't adjust for covariates — mostly superseded • Score test (SAIGE, REGENIE): computationally efficient approximation to full logistic regression, essential at biobank scale with millions of variants × hundreds of thousands of samples • Mixed models (BOLT-LMM, fastGWA): include a random effect capturing cryptic relatedness and fine-scale population structure via a genetic relationship matrix (GRM)
Multiple testing burden: testing ~1 million independent LD blocks across the genome at α=0.05 would yield ~50,000 false positives by chance alone — this motivates the stringent genome-wide significance threshold applied at the plotting stage.
Population stratification is the single most important confounder in GWAS: if cases and controls differ systematically in ancestry, and ancestry correlates with allele frequency at many SNPs, thousands of spurious associations can appear genome-wide.
Principal Component Analysis (PCA) on genome-wide genotype data captures the major axes of ancestry variation: 1. Prune SNPs to near-independence (LD r²<0.1, remove long-range LD regions like HLA and 8p23 inversion) 2. Compute genetic relationship matrix (GRM) or use randomized PCA (FlashPCA2, fastPCA) on ~100k–500k pruned SNPs 3. Top PCs correspond to continental ancestry (PC1–PC2 often separate European/African/East Asian clusters), then finer sub-structure (PC3+) 4. Include top 10–20 PCs as covariates in the association model — this "adjusts away" the ancestry-driven allele frequency differences
Genomic control (λ_GC): λ_GC = median(χ² test statistics) / 0.456
λ_GC=1.0 indicates no inflation; λ_GC>1.05 suggests residual stratification or polygenicity. For highly polygenic traits, LD Score Regression (LDSC) distinguishes true polygenic signal from confounding by regressing χ² statistics against LD scores — genuine polygenicity produces an LDSC intercept near 1.0 even when λ_GC is elevated by many true small effects.
The Manhattan plot is the iconic output of a GWAS: −log10(p-value) plotted against chromosomal position, with chromosomes displayed as alternating color bands left to right. Towers of points rising above the genome-wide significance line (p<5×10⁻⁸) mark loci where genetic variation is robustly associated with disease risk — named for its resemblance to the Manhattan skyline.
Construction of the plot:
X-axis: cumulative genomic position, chromosome 1 through 22 (then X), each chromosome plotted as a contiguous block, alternating shading (e.g. blue/indigo) purely to visually separate adjacent chromosomes
Y-axis: −log10(p-value) for each tested SNP — smaller p-values produce taller points; a p=5×10⁻⁸ sits at y=7.3, p=5×10⁻²⁰ sits at y=19.3
Why 5×10⁻⁸: this threshold approximates a Bonferroni correction for ~1 million independent tests — the effective number of independent LD blocks in the genome for a European-ancestry population (the true number of independent SNPs is far higher, but strong local LD means nearby SNPs are highly correlated, not independent tests). The threshold was formalized by the International HapMap Consortium in 2005 and has remained the field standard since.
Interpreting tower shapes: • Sharp single-SNP spike: rare, possibly low imputation quality — treat cautiously • Broad tower spanning tens of kb: classic LD-driven signal — many correlated SNPs all show association because they are inherited together, not because they are all independently causal • Multiple independent towers on the same chromosome: distinct causal loci, or one broad locus with multiple independent signals revealed by conditional analysis
The QQ-plot companion: observed −log10(p) plotted against the expected distribution under the null. Early deviation from the diagonal (many small p-values) indicates true signal or unaccounted-for stratification; if the entire distribution shifts uniformly upward, elevated λ suggests confounding rather than genuine massively polygenic effects.
A landmark human height GWAS (Yengo et al. 2022, Nature) analyzed 5.4 million individuals and identified over 12,000 independent genome-wide significant SNPs collectively explaining ~40% of height heritability — illustrating that for highly polygenic traits, Manhattan plots become saturated with hundreds of towers rather than a handful of dominant peaks.
A genome-wide significant locus typically contains dozens to hundreds of correlated SNPs spanning the same LD block — the statistical signal alone cannot tell you which variant, if any, is causal. Fine-mapping methods use the fine-grained correlation structure of local haplotypes to narrow a broad association peak down to a small credible set of variants most likely to be functionally responsible.
Linkage disequilibrium (LD) — the non-random co-inheritance of nearby alleles — is both the reason GWAS works with sparse arrays (tag SNPs capture untyped neighbors) and the reason a raw Manhattan peak cannot pinpoint causality.
LD is quantified as r² between two SNPs: the squared correlation of their allele dosages across the population. SNPs with r²>0.8 to the lead SNP are considered to be "tagging" the same underlying signal and are difficult to statistically distinguish.
Bayesian fine-mapping (SuSiE, FINEMAP, PAINTOR): 1. Compute the LD matrix (r between all SNP pairs) within the locus from the study cohort or a matched reference panel 2. Model observed association statistics as a linear combination of a small number of true causal effects plus LD-induced correlation 3. Assign each SNP a Posterior Inclusion Probability (PIP) — the model's estimated probability that this specific SNP is causal 4. Construct a 95% credible set: the smallest set of SNPs whose summed PIP reaches 0.95
SuSiE (Sum of Single Effects, Wang et al. 2020) explicitly models multiple independent causal signals within one locus via iterative Bayesian variable selection, outputting one credible set per detected independent signal.
Functional annotation weighting (PolyFun, PAINTOR) incorporates external biological priors — chromatin accessibility (ATAC-seq), transcription factor binding, conservation scores — to sharpen posterior probabilities toward functionally plausible variants.
Once a credible set is defined, translating statistical signal into biological hypothesis requires gene annotation:
Nearest-gene heuristic: simplest approach, assign locus to physically closest gene — correct roughly 60–70% of the time but frequently wrong, since regulatory variants can act on genes hundreds of kb away, skipping the nearest gene entirely
Coding variant priority: if the credible set contains a missense, nonsense, or splice-site variant, that gene is prioritized as the likely causal mechanism (protein-altering effect is directly interpretable)
Regulatory annotation: most GWAS hits (~90%) fall in non-coding regions — enhancers, promoters, and other regulatory elements. Overlap with ENCODE/Roadmap chromatin state maps, ATAC-seq peaks, and Hi-C/promoter-capture chromatin interaction data helps link a non-coding variant to its regulated gene, which may not be the nearest one
eQTL/pQTL overlap: if the credible set SNP is also an expression or protein quantitative trait locus (see functional follow-up stage) in a disease-relevant tissue, this substantially strengthens gene assignment
Example: the FTO obesity locus was long assumed to act through the FTO gene itself (nearest gene), until functional studies revealed the causal enhancer variant actually regulates IRX3/IRX5 expression in adipocyte precursors, hundreds of kilobases away — a cautionary tale for naive nearest-gene assignment.
A genome-wide significant, fine-mapped locus is still just a statistical statement until functional and replication evidence establishes a plausible causal mechanism. Colocalization with gene expression, cross-referencing the GWAS Catalog, and independent replication cohorts together transform a Manhattan plot peak into an actionable biological or therapeutic hypothesis.
Expression quantitative trait loci (eQTLs) are genetic variants associated with the expression level of a nearby (cis-eQTL) or distant (trans-eQTL) gene. The GTEx (Genotype-Tissue Expression) project catalogs eQTLs across 54 human tissues from ~950 donors, providing the reference resource for interpreting GWAS loci mechanistically.
Colocalization analysis (coloc, Giambartolomei et al. 2014) tests whether a GWAS signal and an eQTL signal at the same locus share a single underlying causal variant, versus arising from two distinct nearby variants in LD:
Five competing hypotheses per locus: • H0: no association with either trait • H1: association with GWAS trait only • H2: association with expression only • H3: both traits associated, but distinct causal variants (LD confounding) • H4: both traits share a single causal variant (true colocalization)
Posterior probability PP4 (support for H4) >0.8 is generally considered strong evidence that the GWAS association acts through modulating expression of that specific gene in that specific tissue — a powerful mechanistic anchor connecting a statistical peak to a testable hypothesis (e.g. "this variant increases disease risk by decreasing expression of Gene X in pancreatic islet cells").
Similar colocalization frameworks apply to protein QTLs (pQTLs, from Olink/SomaScan proteomic panels) and splicing QTLs (sQTLs), each adding independent mechanistic evidence.
The NHGRI-EBI GWAS Catalog aggregates published genome-wide significant associations across thousands of studies and traits, enabling rapid cross-referencing: is this locus already known for a related phenotype? Pleiotropy — one locus affecting multiple traits — is extremely common and often mechanistically informative (e.g. a lipid locus also associated with coronary artery disease strengthens a causal lipid-to-CAD hypothesis).
Replication requirements: • Independent cohort, ideally different ancestry/recruitment source, tests the same lead SNP • Consistent direction of effect and p<0.05 (or Bonferroni-adjusted for number of loci tested) considered replicated • Meta-analysis across discovery + replication cohorts (fixed or random-effects models, inverse-variance weighted) produces the final consensus effect size
Mendelian randomization (MR): uses genetic variants as instrumental variables to test causal relationships between an exposure (e.g. LDL cholesterol) and an outcome (e.g. heart disease), exploiting the fact that genotype is randomly assigned at conception and thus immune to reverse causation and much confounding — a genetic analog of a randomized controlled trial.
Translational impact: drug targets with human genetic support (i.e. a GWAS or rare-variant association pointing to the target gene) succeed in clinical trials roughly twice as often as targets without such support (Nelson et al. 2015; King et al. 2019) — making GWAS-derived target validation a cornerstone of modern drug discovery pipelines, exemplified by PCSK9 inhibitors for cardiovascular disease, originally motivated by human loss-of-function genetic variants.
PCSK9 is the textbook translational success story: rare loss-of-function variants discovered through human genetics were associated with lifelong low LDL cholesterol and dramatically reduced coronary artery disease risk with no adverse phenotype — directly motivating development of PCSK9 inhibitor antibodies (evolocumab, alirocumab), now approved cholesterol-lowering therapeutics.