Bio Population Genetics Association Testing

Version Compatibility

Reference examples tested with: PLINK 2.0 (alpha 6+), SAIGE 1.3+, regenie 3.4+, BOLT-LMM 2.4+, GEMMA 0.98+, numpy 1.26+, pandas 2.2+, scipy 1.12+.

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures
• CLI: <tool> --version then <tool> --help to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Version traps that change results, not just syntax: PLINK 2.0 --glm firth-fallback is the DEFAULT for binary traits and writes .glm.logistic.hybrid (mixed logistic and Firth rows), not .glm.logistic. SAIGE --LOCO=TRUE is the recommended default and the step-1 and step-2 sample IDs plus variance-ratio file must match exactly. regenie applies Firth/SPA only when --firth/--spa are explicitly set in step 2 (not automatic for every variant). BOLT-LMM is calibrated only for quantitative traits at large N (case fraction >= 10%, MAF > 0.1% for binary coding). The single source of truth for versions is this block, not headings.

Association Testing

"Run a GWAS on my genotypes" -> Fit one regression per variant (additive dosage as predictor) under a confounder model chosen to make genotype independent of unmodeled phenotype drivers, then read effect sizes and p-values that are honest only insofar as that model holds.

• CLI: plink2 --glm firth-fallback hide-covar --covar pcs.eigenvec for an unrelated sample whose structure is captured by PCs
• CLI: SAIGE (SPA) or regenie --step 1/2 --firth for relatedness and/or case:control imbalance; GEMMA -lmm or BOLT-LMM --lmm for related quantitative traits

Scope: single-variant common-variant GWAS (linear/logistic GLM, linear mixed models, SPA, Firth) and the choice between them. Rare-variant GENE-BASED AGGREGATION (burden, SKAT, SKAT-O, ACAT, STAAR) routes to rare-variant-association. Fine-mapping, Mendelian randomization, and TWAS route to causal-genomics/fine-mapping and causal-genomics/mendelian-randomization. Polygenic risk scores route to clinical-databases/polygenic-risk. Imputed dosages enter from phasing-imputation/genotype-imputation.

The Single Most Important Insight -- a GWAS statistic is valid only if genotype is independent of unmodeled phenotype drivers AFTER the chosen covariates and random effects

1. Every classic GWAS pathology is one violation of that assumption with a specific signature: structure or cryptic relatedness lifts the whole QQ plot; a heritable covariate adjusted as a collider manufactures direction-flipped hits; case:control imbalance with low MAC under a score test makes the significant tail anti-conservative.
2. The corollary that organizes the toolchain: under a real polygenic trait genomic inflation (lambda) is EXPECTED to exceed 1 and is mostly true signal, so genomic control over-corrects and erases discoveries; the LDSC intercept (not lambda) separates confounding from polygenicity (Bulik-Sullivan 2015).
3. The fix differs per distortion, and the wrong fix makes it worse: lambda-correcting a polygenic trait, or adding PCs to a family sample, both silently degrade the result.
4. Choosing the engine is choosing the assumption: a PC-adjusted GLM assumes structure is a low-rank mean shift, an LMM assumes a polygenic covariance, SPA/Firth assume the score/Wald tail needs the cumulant-generating-function correction.

Tool Taxonomy

Decision Tree by Scenario

plink2 GLM (unrelated sample, PC-adjusted)

# Logistic for binary, linear for quantitative is auto-detected from the phenotype coding.
# firth-fallback is the binary-trait default: ordinary logistic, falling back to Firth only on
# non-convergence (separation). Output is .glm.logistic.hybrid (mixed logistic and Firth rows).
plink2 --bfile qc --pheno pheno.txt --glm firth-fallback hide-covar \
    --covar pcs.eigenvec --covar-name PC1-PC10 --out gwas

# Regress on imputed DOSAGES, not hard calls, so imputation uncertainty enters the SE. Reporting
# columns add A1 frequency and the imputation R2 (machr2 is meaningful only on dosage input) for
# downstream QC and effect-allele bookkeeping.
plink2 --pfile imputed --pheno pheno.txt --glm firth-fallback hide-covar cols=+a1freq,+machr2 \
    --covar covars.txt --covar-name PC1-PC10,age,sex --out gwas_dosage

PCs must be computed on LD-pruned, MAF-filtered genotypes with long-range-LD regions excluded (see population-structure); too few PCs leave residual stratification, too many absorb real signal. --glm sex adds sex as a covariate on chrX (no-x-sex suppresses it); split the pseudoautosomal region before any chrX test (see plink-basics).

Linear mixed model with LOCO (relatedness / structure)

# GEMMA: -gk builds the GRM (1=centered, 2=standardized), then -lmm fits one variance component.
# The covariate file passed to -c MUST contain an explicit intercept column of 1s (GEMMA does not add one).
gemma -bfile qc -gk 1 -o grm
gemma -bfile qc -k output/grm.cXX.txt -c covars_with_intercept.txt -lmm 4 -o lmm   # -lmm 4 = Wald+LRT+score

# BOLT-LMM: --lmm decides by cross-validation whether the non-infinitesimal model adds power.
# --lmmInfOnly forces the standard infinitesimal model; --lmmForceNonInf forces the mixture.
bolt --bfile=qc --phenoFile=pheno.txt --phenoCol=trait \
    --covarFile=covars.txt --qCovarCol=PC{1:10} --lmm --LDscoresFile=LDSCORE.1000G_EUR.tab.gz \
    --statsFile=bolt.stats

# SAIGE: step 1 fits the null GLMM (sparse GRM, variance ratio); step 2 runs the SPA score test per variant.
# --LOCO=TRUE (default, recommended) excludes the tested chromosome from the polygenic predictor.
step1_fitNULLGLMM.R --plinkFile=qc --phenoFile=pheno.txt --phenoCol=trait --traitType=binary \
    --covarColList=PC1,PC2,age,sex --sampleIDColinphenoFile=IID --outputPrefix=step1 --LOCO=TRUE
step2_SPAtests.R --vcfFile=chr1.vcf.gz --GMMATmodelFile=step1.rda --varianceRatioFile=step1.varianceRatio.txt \
    --minMAC=20 --is_Firth_beta=TRUE --LOCO=TRUE --SAIGEOutputFile=chr1.saige

# regenie: step 1 builds the whole-genome LOCO ridge predictor; step 2 tests with Firth (--approx for speed) or SPA.
regenie --step 1 --bed qc --phenoFile pheno.txt --covarFile covars.txt --bt --bsize 1000 --out step1
regenie --step 2 --bed qc --phenoFile pheno.txt --covarFile covars.txt --bt \
    --firth --approx --pThresh 0.01 --pred step1_pred.list --bsize 400 --out step2

LOCO (leave-one-chromosome-out) is not optional: if the tested variant's chromosome is in the GRM/predictor, the random effect explains part of the variant's own signal and deflates power (proximal contamination). A hand-rolled "GRM from all SNPs" silently throws away power at every true locus.

Diagnose inflation: lambda vs the LDSC intercept

import numpy as np
from scipy import stats

def lambda_gc(pvalues):
    chisq = stats.chi2.ppf(1 - pvalues, 1)
    return np.median(chisq) / stats.chi2.ppf(0.5, 1)

# lambda > 1 under a polygenic trait is EXPECTED and mostly true signal. Do NOT divide statistics by it.
# Rescale to lambda_1000 (per 1000 cases/1000 controls) before comparing studies of different N.
def lambda_1000(lam, n_cases, n_controls):
    return 1 + (lam - 1) * (1 / n_cases + 1 / n_controls) / (1 / 1000 + 1 / 1000)

# The confounding diagnostic is the LDSC INTERCEPT (run ldsc on the sumstats), not lambda:
# intercept ~ 1 with high lambda = polygenicity (clean); intercept materially > 1 = confounding.
# Prefer the attenuation ratio = (intercept - 1) / (mean(chi2) - 1): the fraction of inflation NOT due
# to polygenicity (~0-0.2 acceptable). The intercept is also inflated by sample overlap in bivariate LDSC.

Per-Method Failure Modes

Genomic control on a polygenic trait

Trigger: dividing every chi-square by lambda_GC because lambda > 1.1. Mechanism: lambda rises with N and h2 under true polygenicity, so it is mostly signal. Symptom: discoveries vanish; power destroyed. Fix: never lambda-correct on lambda alone; use the LDSC intercept and only deflate if intercept-minus-1 is materially > 0.

PCs for a related sample

Trigger: adding "10 PCs" to a family or cryptically-related cohort. Mechanism: relatedness is a pairwise covariance, not a low-rank mean shift, so no finite PC set removes it. Symptom: inflated, miscalibrated tail statistics despite the PCs. Fix: switch to an LMM (GEMMA/BOLT/SAIGE/regenie); more PCs is the wrong lever.

LMM without LOCO

Trigger: building the GRM from all chromosomes including the candidate. Mechanism: the variant partly explains itself as a random effect. Symptom: deflated statistics, lost power at true loci. Fix: use the LOCO variant (--mlma-loco, SAIGE --LOCO=TRUE, regenie step-1 predictor).

Wald collapse / score anti-conservatism at imbalance

Trigger: plain logistic for a rare or near-monomorphic-in-one-arm variant, or a score test at extreme case:control. Mechanism: the MLE diverges (Wald BETA/SE -> 0) or the score null is wrong in the tail. Symptom: real rare-variant signal looks non-significant, or the significant tail fills with artifacts. Fix: Firth (plink2 firth-fallback, regenie --firth) or SPA (SAIGE, regenie --spa).

HWE filtered in cases

Trigger: HWE filter on the case-only or combined sample as a discovery QC step. Mechanism: a true non-additive disease variant legitimately deviates from HWE in cases. Symptom: real associations removed before testing. Fix: compute HWE in CONTROLS (or founders) only (see plink-basics).