Bio Workflows Gwas Pipeline

Version Compatibility

Reference examples tested with: PLINK 2.0 (alpha 5+), regenie 3.4+, SAIGE 1.3+, Eagle 2.4+ / SHAPEIT5, Minimac4 / Beagle 5.4, bcftools 1.19+, LDSC 1.0, qqman 0.1.9+, ggplot2 3.5+

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

• R: packageVersion('<pkg>') then ?function_name to verify parameters
• 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.

Note: regenie/SAIGE run a two-step design (step 1 whole-genome ridge / null GLMM on LD-pruned common markers; step 2 tests the imputed set with LOCO) — step-1 and step-2 sample IDs + variance-ratio file MUST match or calibration silently fails. Binary PLINK2 --glm defaults to firth-fallback (output .glm.logistic.hybrid). Sequencing/WGS and non-EUR LD shift the genome-wide threshold off the 5e-8 EUR-array folklore. Confirm in-tool before quoting.

GWAS Pipeline

"Run a GWAS from my genotype data" -> Orchestrate sample/variant QC (PLINK2), population stratification (PCA), association testing (linear/logistic regression), multiple testing correction, and Manhattan/QQ plot visualization.

This is a workflow skill: it owns the chaining decisions and hand-offs, not the internals of any one step. Every step below cross-references the component skill that teaches its mechanism.

The governing principle

A GWAS result is decided at four seams, not inside the association test.

1. The genome build + ancestry-matched imputation panel is a made-once commitment. The panel fixes what can be imputed (HRC MAF floor ~5e-4; TOPMed imputes far rarer) AND the build (HRC/1000G-P3 are GRCh37; TOPMed/1000G-NYGC are GRCh38). Ancestry match > panel size: a bigger ancestry-mismatched panel imputes WORSE because there are no matching haplotypes to copy.
2. Strand/allele harmonization is a made-once GATE, and the silent corruptor. Align every study variant's alleles to the panel REF/ALT, fix strand, and DROP unresolvable palindromes (A/T, C/G at MAF>0.4). A flipped palindrome or a build mismatch flips BETA and cancels/manufactures signal WITHOUT erroring — caught only on the AF-concordance plot (points on the y=1-x anti-diagonal), not by a crash.
3. QC runs before association and in a defensible internal order; the engine is chosen by structure, not convenience. Variant missingness before sample missingness; HWE in CONTROLS only (a true non-additive risk variant deviates in cases); impute cases+controls TOGETHER; exclude long-range-LD regions before PCA. PC-covariate GLM is valid only for unrelated, continuous-ancestry cohorts — relatedness/structure needs an LMM (regenie/SAIGE/BOLT) with LOCO, because PCs remove a mean shift, not a covariance.
4. Refuse the bare lambda. lambda_GC alone cannot separate polygenicity from confounding — read it WITH the LDSC intercept (~1 = polygenic, fine; >1 = confounding). Do not reflexively genomic-control (it over-corrects true signal). Filter INFO/R2 MAF-stratified (a flat cutoff is a hidden rare-variant filter).

Made-once commitments

Workflow Overview

VCF/PLINK files
    |
    v
[1. QC Filtering] ------> Sample and variant QC
    |
    v
[2. Phase + Impute] ----> Align to panel, phase, impute to dosages, filter by R2 (-> phasing-imputation)
    |
    v
[3. LD Pruning] --------> Independent variants for PCA
    |
    v
[4. Population Structure] --> PCA for covariates
    |
    v
[5. Association Testing] --> Logistic/linear regression on dosages
    |
    v
[6. Results] -----------> Manhattan plot, QQ plot
    |
    v
Significant associations

Step 1: Data Import and QC

Convert VCF to PLINK

# VCF to PLINK binary format
plink2 --vcf input.vcf.gz \
    --pheno phenotypes.txt \
    --make-bed \
    --out study   # load --pheno at conversion so PHENO1 is in the .fam for the controls-only HWE gate below

QC order is critical (-> population-genetics/plink-basics): variant missingness runs FIRST, in its own invocation, so a sample is not dropped for missingness driven by variants slated for removal.

Variant call-rate (first)

plink2 --bfile study --missing --out study_stats
plink2 --bfile study --geno 0.05 --make-bed --out study_var_qc   # variant missingness BEFORE sample missingness

Sample QC

# Sample missingness AFTER variant missingness.
plink2 --bfile study_var_qc --mind 0.05 --make-bed --out study_sample_qc

# Sex check: split the pseudoautosomal region first or male PAR het reads as a sex error (see plink-basics).
plink2 --bfile study_sample_qc --split-par hg38 --check-sex --out study_sex_check

# KING-robust relatedness - structure-robust and IBD-free (this is the point of KING; no --genome needed).
plink2 --bfile study_sample_qc --king-cutoff 0.0884 --make-bed --out study_unrelated

Variant QC: MAF and controls-only HWE

# Controls-only HWE with mid-p. plink2 is NOT controls-only by default, so gate to controls explicitly:
# a true risk variant depletes heterozygotes in cases and would fail a case-inclusive HWE test (see plink-basics).
plink2 --bfile study_unrelated --keep-if "PHENO1 == control" --hwe 1e-6 midp --write-snplist --out hwe_pass
plink2 --bfile study_unrelated --maf 0.01 --extract hwe_pass.snplist --make-bed --out study_qc
plink2 --bfile study_qc --freq --out study_qc

QC Checkpoint:

• Variant call rate >95% (applied before sample missingness)
• Sample call rate >95%
• MAF >1%
• HWE applied to controls only with mid-p (p>1e-6)

Step 2: Phasing and Imputation

Most GWAS impute the QC'd array genotypes up to a dense reference panel before association, to increase variant density and harmonize across platforms and studies. This stage is owned by the phasing-imputation skills; the workflow only orchestrates the handoff. The decisions that matter here, in order:

1. Select and prepare the panel (-> phasing-imputation/reference-panels). Match the panel ancestry to the cohort (TOPMed for diverse/admixed, HRC or 1000G for European, HGDP+1kGP for diverse-and-downloadable), reconcile genome build, and run the strand/allele harmonization check; a flipped palindromic SNP or build mismatch corrupts results without erroring.
2. Phase, then impute (-> phasing-imputation/haplotype-phasing, phasing-imputation/genotype-imputation). Pre-phase the QC'd VCF (Eagle2/SHAPEIT5) and impute against the panel per chromosome (Beagle/Minimac4/IMPUTE5), or upload to the Michigan/TOPMed server for controlled-access panels. Impute cases and controls TOGETHER; separate imputation manufactures false associations. The output carries dosages (DS), not hard calls.
3. Filter by quality (-> phasing-imputation/imputation-qc). Drop variants below an INFO/R2/DR2 cutoff paired with a MAF floor, MAF-stratified, because a flat cutoff is a hidden rare-variant filter. Carry dosages forward.

Goal: Increase variant density and harmonize across platforms by imputing the QC'd genotypes against a dense ancestry-matched panel, carrying dosages (not hard calls) into association.

Approach: Export the QC'd genotypes to VCF, hand off to the phasing-imputation skills for strand/build harmonization, phasing, and per-chromosome imputation (cases and controls together), then filter on the engine's quality field plus a MAF floor.

# Convert QC'd PLINK back to VCF, align to the panel, phase + impute (see phasing-imputation skills for the full commands)
plink2 --bfile study_qc --export vcf bgz --out study_qc
# ... reference-panels: strand/build harmonization; haplotype-phasing: phase; genotype-imputation: impute to dosages ...
# Post-imputation quality filter on the engine's field (DR2 Beagle / R2 Minimac), with a MAF floor
bcftools view -e 'INFO/DR2<0.3 || INFO/AF<0.01 || INFO/AF>0.99' imputed.vcf.gz -Oz -o imputed.qc.vcf.gz

QC Checkpoint: cases and controls imputed together; INFO/R2 filtered (MAF-stratified) with a MAF floor; dosages (DS), not hard calls, carried into association.

Step 3: LD Pruning for PCA

# Exclude long-range-LD regions and inversions FIRST (MHC, 8p23.1, 17q21.31, LCT) - their internal r2 is
# high and real, so a window prune keeps them and a PC then tracks the inversion (-> population-structure).
plink2 --bfile study_qc --exclude range longrange_ld.txt --make-bed --out study_noLR

# Identify independent variants (r2 0.1 matches the linkage-disequilibrium / population-structure default).
plink2 --bfile study_noLR --indep-pairwise 50 5 0.1 --out pruned
plink2 --bfile study_noLR --extract pruned.prune.in --make-bed --out study_pruned

Step 4: Population Structure (PCA)

# Calculate principal components
plink2 --bfile study_pruned \
    --pca 10 \
    --out study_pca

# The eigenvec file contains PCs for use as covariates

Visualize PCA

library(ggplot2)

# Load PCA results
pca <- read.table('study_pca.eigenvec', header = FALSE)
colnames(pca) <- c('FID', 'IID', paste0('PC', 1:10))

# Load phenotype for coloring
pheno <- read.table('phenotypes.txt', header = TRUE)
pca <- merge(pca, pheno, by = c('FID', 'IID'))

# Plot
ggplot(pca, aes(x = PC1, y = PC2, color = as.factor(PHENO))) +
    geom_point(alpha = 0.5) +
    labs(title = 'PCA of Study Samples', color = 'Phenotype') +
    theme_minimal()
ggsave('pca_plot.pdf', width = 8, height = 6)