/plink2-gwas-analysis

name: "plink2-gwas-analysis"
description: "GWAS and population genetics tool. Processes PLINK (.bed/.bim/.fam), VCF, and BGEN; runs QC (MAF, HWE, missingness), IBD estimation, PCA, and linear/logistic regression GWAS. Outputs Manhattan-ready summary stats. Use regenie or SAIGE for biobanks (>100k samples) needing mixed models."
license: "GPL-3.0"

PLINK2 — GWAS and Population Genetics

Overview

PLINK2 is the high-performance successor to PLINK 1.9, designed for genome-wide association studies (GWAS) and population genetics analysis on large cohorts. It processes genotype data in PLINK binary format (.bed/.bim/.fam), VCF, and BGEN formats — performing sample and variant quality control (QC), kinship estimation, principal component analysis (PCA), and linear/logistic regression association testing. PLINK2 is 10–100× faster than PLINK 1.9 on most tasks due to multithreading and optimized I/O. Output files are compatible with downstream visualization (Manhattan/QQ plots) and meta-analysis tools.

When to Use

• Running GWAS on a case-control or quantitative trait cohort after genotyping array QC
• Performing sample QC: missingness, heterozygosity outliers, sex check, cryptic relatedness
• Computing genome-wide LD pruning for PCA or relatedness estimation
• Running PCA on genotype data to identify population stratification
• Converting between PLINK binary, VCF, and BGEN formats
• Filtering variants by MAF, HWE, missingness, or INFO score in VCF/imputed data
• Use **omics-plotting** SKILL to render Manhattan and QQ plots from the association results
• Use **regenie** or **SAIGE** instead for biobank-scale GWAS (>100k samples) requiring mixed model association to control for population structure
• Use **VCFtools** as an alternative for VCF-specific population genetics statistics

Prerequisites

• **Software**: PLINK2 (pre-compiled binary; no pip/conda package)
• **Input**: PLINK binary files (.bed/.bim/.fam) or VCF/BGEN from array genotyping or imputation

> **Check before installing**: The tool may already be available (e.g., inside a `pixi` / `conda` env). Always run `command -v plink2` first and skip the install block if it returns a path. When executing tools inside a pixi project, prefer `pixi run <tool>` over plain `<tool>`.

# Skip install if already present
if command -v plink2 >/dev/null 2>&1; then
    echo "plink2 already installed: $(plink2 --version)"
else
    # Download PLINK2 pre-compiled binary (Linux)
    wget https://s3.amazonaws.com/plink2-assets/alpha6/plink2_linux_avx2_20241112.zip
    unzip plink2_linux_avx2_20241112.zip
    chmod +x plink2
    export PATH="$PWD:$PATH"

    # macOS
    # wget https://s3.amazonaws.com/plink2-assets/alpha6/plink2_mac_20241112.zip
    # unzip plink2_mac_20241112.zip

    plink2 --version
    # PLINK v2.00a6LM
fi

# Python for downstream analysis
pip install pandas numpy matplotlib scipy

Pre-flight Interview

Settle these with the user before writing any analysis code.

decisions:
  - id: D1
    param: phenotype
    kind: required
    source: data
    ask: "Which column holds the trait being tested, and is it a measurement or a case/control status?"
    default: null

  - id: D2
    param: regressionMode
    kind: derived
    source: data
    depends_on: [D1]
    ask: "Should the association test be linear or logistic?"
    default: "linear for a quantitative trait, logistic for case/control"

  - id: D3
    param: qcThresholds
    kind: required
    source: user
    ask: "Where should variants and samples be cut for allele frequency, missingness, and Hardy-Weinberg departure?"
    default: "MAF 0.01, variant missingness 0.05, sample missingness 0.02, HWE 1e-6"

  - id: D4
    param: covariates
    kind: required
    source: data
    ask: "Which covariates belong in the model - age, sex, batch, study site?"
    default: null

  - id: D5
    param: populationStructure
    kind: required
    source: user
    ask: "How many principal components should be included to absorb ancestry differences?"
    default: "10, computed from LD-pruned variants"

  - id: D6
    param: ldPruning
    kind: optional
    source: user
    depends_on: [D5]
    ask: "Which window and correlation threshold should thin variants before computing the components?"
    default: "50 variant window, step 5, r-squared 0.2"

  - id: D7
    param: threads
    kind: never_ask
    source: data
    reason: "Affects runtime only, not the association statistics"
    default: "min(8, available_cores)"

D5 hangs on D4 because the components are themselves covariates - how many to include is a question about the covariate set, not a separate knob. Omitting them entirely is the classic way to produce a genome-wide significant result that reflects ancestry rather than the trait.

Quick Start

# Run GWAS: linear regression for quantitative trait
plink2 \
    --bfile cohort_qc \
    --pheno phenotypes.txt \
    --covar covariates.txt \
    --linear hide-covar \
    --out results/gwas_result \
    --threads 8

# View top hits
head results/gwas_result.*.glm.linear | sort -k12,12g | head -20

Workflow

Step 1: Convert VCF/Imputed Data to PLINK Binary Format

Convert input genotype data to PLINK binary format for fast processing.

# Convert VCF to PLINK binary
plink2 \
    --vcf cohort_genotyped.vcf.gz \
    --make-bed \
    --out cohort_plink \
    --threads 8

# Convert BGEN (imputed data from Michigan/TopMed imputation server)
plink2 \
    --bgen cohort_imputed.bgen ref-first \
    --sample cohort_imputed.sample \
    --make-bed \
    --out cohort_imputed_plink \
    --threads 8

echo "Files created:"
ls -lh cohort_plink.{bed,bim,fam}
echo "Samples: $(wc -l < cohort_plink.fam)"
echo "Variants: $(wc -l < cohort_plink.bim)"

Step 2: Sample QC — Missingness and Heterozygosity

Remove samples with high missingness or heterozygosity outliers.

# Compute per-sample and per-variant missingness
plink2 \