Bio Population Genetics Association Testing

Association Testing

GWAS analysis using PLINK 2.0's unified --glm command for case-control and quantitative traits.

PLINK 2.0 Association Testing

Basic Case-Control (Binary Phenotype)

# Basic logistic regression
plink2 --bfile data --glm --out results

# With phenotype file
plink2 --bfile data --pheno pheno.txt --glm --out results

Quantitative Trait (Continuous Phenotype)

# Linear regression for quantitative traits
plink2 --bfile data --pheno pheno.txt --glm --out results

With Covariates

# Include covariates (sex, age, PCs)
plink2 --bfile data \
    --pheno pheno.txt \
    --covar covariates.txt \
    --glm --out results

# Specify which covariates to use
plink2 --bfile data \
    --pheno pheno.txt \
    --covar covariates.txt \
    --covar-name PC1,PC2,PC3,age,sex \
    --glm --out results

Covariate Files

Phenotype File Format

# pheno.txt: FID IID pheno
# For binary: 1=control, 2=case, -9=missing
# For quantitative: continuous values
FAM001 IND001 2
FAM002 IND002 1
FAM003 IND003 1.5

Covariate File Format

# covariates.txt: FID IID cov1 cov2 ...
FAM001 IND001 0.15 35 1
FAM002 IND002 -0.22 42 2
FAM003 IND003 0.08 28 1

GLM Options

Phenotype Handling

# Multiple phenotypes (test all)
plink2 --bfile data --pheno pheno_multi.txt --glm --out results

# Specific phenotype column
plink2 --bfile data --pheno pheno_multi.txt --pheno-name trait1 --glm --out results

# Missing phenotype handling
plink2 --bfile data --glm allow-no-covars --out results

Model Options

# Additive model (default)
plink2 --bfile data --glm --out results

# Dominant model
plink2 --bfile data --glm dominant --out results

# Recessive model
plink2 --bfile data --glm recessive --out results

# Genotypic (2df test)
plink2 --bfile data --glm genotypic --out results

# Hide covariates from output (cleaner output)
plink2 --bfile data --covar cov.txt --glm hide-covar --out results

Firth Regression (Rare Variants)

# Enable Firth fallback for case-control (default in PLINK 2.0)
plink2 --bfile data --glm firth-fallback --out results

# Force Firth regression
plink2 --bfile data --glm firth --out results

# Disable Firth
plink2 --bfile data --glm no-firth --out results

Output Format

Output Columns

# Default output: results.PHENO1.glm.logistic or results.PHENO1.glm.linear
# Columns: CHROM, POS, ID, REF, ALT, A1, FIRTH?, TEST, OBS_CT, OR/BETA, SE, Tstat, P

Custom Output Columns

# Add specific columns
plink2 --bfile data --glm cols=+a1freq,+machr2 --out results

# Available columns:
# +a1freq: A1 allele frequency
# +machr2: MaCH R-squared
# +ax: Reference allele dosage
# +err: Standard errors

Population Stratification Control

Include Principal Components

# 1. Run PCA
plink2 --bfile data --pca 10 --out pca_results

# 2. Use PCs as covariates
plink2 --bfile data \
    --pheno pheno.txt \
    --covar pca_results.eigenvec \
    --covar-name PC1,PC2,PC3,PC4,PC5 \
    --glm --out results

Combined Workflow

# QC, PCA, and GWAS in sequence
plink2 --bfile raw --maf 0.01 --geno 0.05 --hwe 1e-6 --make-bed --out qc
plink2 --bfile qc --pca 10 --out pca
plink2 --bfile qc \
    --pheno pheno.txt \
    --covar pca.eigenvec \
    --covar-name PC1-PC5 \
    --glm hide-covar --out gwas

Result Filtering

Command Line Filtering

# Filter significant results
awk 'NR==1 || $13 < 5e-8' results.PHENO1.glm.logistic > significant.txt

# Extract top hits
sort -k13 -g results.PHENO1.glm.logistic | head -100 > top_hits.txt

Python Analysis

import pandas as pd

results = pd.read_csv('results.PHENO1.glm.logistic', sep='\t')

significant = results[results['P'] < 5e-8]
print(f'Genome-wide significant hits: {len(significant)}')

suggestive = results[results['P'] < 1e-5]
print(f'Suggestive hits: {len(suggestive)}')

Visualization

Manhattan Plot (Python)

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np

results = pd.read_csv('results.PHENO1.glm.logistic', sep='\t')
results = results[results['TEST'] == 'ADD']
results['-log10P'] = -np.log10(results['P'])

chrom_colors = ['#1f77b4', '#ff7f0e']
results['color'] = results['#CHROM'].apply(lambda x: chrom_colors[x % 2])

cumulative_pos = []
offset = 0
for chrom in sorted(results['#CHROM'].unique()):
    chrom_data = results[results['#CHROM'] == chrom]
    cumulative_pos.extend(chrom_data['POS'] + offset)
    offset += chrom_data['POS'].max()

results['cumulative_pos'] = cumulative_pos

plt.figure(figsize=(14, 6))
plt.scatter(results['cumulative_pos'], results['-log10P'], c=results['color'], s=1)
plt.axhline(y=-np.log10(5e-8), color='red', linestyle='--', label='Genome-wide (5e-8)')
plt.axhline(y=-np.log10(1e-5), color='blue', linestyle='--', label='Suggestive (1e-5)')
plt.xlabel('Chromosome')
plt.ylabel('-log10(P)')
plt.legend()
plt.savefig('manhattan.png', dpi=150)

QQ Plot (Python)

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

results = pd.read_csv('results.PHENO1.glm.logistic', sep='\t')
observed_p = results[results['TEST'] == 'ADD']['P'].dropna().sort_values()

n = len(observed_p)
expected_p = np.arange(1, n + 1) / (n + 1)

plt.figure(figsize=(6, 6))
plt.scatter(-np.log10(expected_p), -np.log10(observed_p), s=1)
plt.plot([0, 8], [0, 8], 'r--')
plt.xlabel('Expected -log10(P)')
plt.ylabel('Observed -log10(P)')

lambda_gc = np.median(stats.chi2.ppf(1 - observed_p, 1)) / stats.chi2.ppf(0.5, 1)
plt.title(f'QQ Plot (λ = {lambda_gc:.3f})')
plt.savefig('qqplot.png', dpi=150)

Genomic Inflation

from scipy import stats
import numpy as np

results = pd.read_csv('results.PHENO1.glm.logistic', sep='\t')
pvalues = results[results['TEST'] == 'ADD']['P'].dropna()

chisq = stats.chi2.ppf(1 - pvalues, 1)
lambda_gc = np.median(chisq) / stats.chi2.ppf(0.5, 1)
print(f'Genomic inflation factor: {lambda_gc:.3f}')
# Good: 1.0-1.05, Acceptable: 1.05-1.1, Concerning: >1.1

Related Skills

• plink-basics - Data preparation and QC
• population-structure - PCA for stratification control
• linkage-disequilibrium - LD pruning before analysis

Association Testing - Usage Guide

Overview

GWAS identifies genetic variants associated with traits using regression models. PLINK 2.0's --glm command provides unified testing for binary (case-control) and quantitative traits with covariate support.

Prerequisites

conda install -c bioconda plink2

# For visualization
pip install pandas matplotlib scipy

Quick Start

Tell your AI agent what you want to do:

• "Run a GWAS for my case-control phenotype"
• "Test association with a quantitative trait"
• "Perform association testing with population covariates"
• "Find genome-wide significant hits"
• "Generate a Manhattan plot from my GWAS results"

Example Prompts

Basic Association

"Run a genome-wide association study on my case-control data"

"Test SNP associations with my quantitative phenotype"

"Perform logistic regression GWAS for disease status"

With Covariates

"Run GWAS including age, sex, and the first 5 PCs as covariates"

"Test associations while controlling for population stratification"

"Perform association testing with a custom covariate file"

Results Analysis

"Extract all genome-wide significant variants from my GWAS"

"Calculate genomic inflation factor for my results"

"Find the top 100 hits and annotate them with gene names"

"Create Manhattan and QQ plots from my association results"

What the Agent Will Do

1. Verify input data quality (post-QC PLINK files) 2. Check phenotype file format and distribution 3. Generate PCs for stratification control if needed 4. Run association testing with appropriate model 5. Filter results by significance threshold 6. Calculate genomic inflation (lambda) 7. Generate visualization if requested

Tips

• Always include population PCs as covariates to control stratification
• Start with --glm hide-covar to simplify output
• Genome-wide significance is 5e-8; suggestive is 1e-5
• Lambda > 1.1 suggests residual stratification or relatedness
• For family data, use mixed models (GCTA, BOLT-LMM) instead

Standard GWAS Workflow

1. Quality Control

plink2 --bfile raw \
    --maf 0.01 --geno 0.05 --mind 0.05 --hwe 1e-6 \
    --make-bed --out qc

2. Calculate PCs for Stratification

plink2 --bfile qc --pca 10 --out pca

3. Run Association

plink2 --bfile qc \
    --pheno phenotypes.txt \
    --covar pca.eigenvec \
    --covar-name PC1-PC5 \
    --glm hide-covar \
    --out gwas

4. Identify Significant Hits

awk '$13 < 5e-8' gwas.PHENO1.glm.* > significant.txt

Significance Thresholds

Common Issues

High Genomic Inflation (lambda > 1.1)

• Add more PCs as covariates
• Check for cryptic relatedness
• Consider mixed models (GCTA, BOLT-LMM)

No Significant Results

• Check phenotype file format
• Verify sample sizes
• May be underpowered

Separation Issues (Logistic)

• Firth regression automatically applied
• Check for very rare variants

Resources

• PLINK 2.0 GLM Documentation
• GWAS Tutorial