plink2-gwas-analysis

PLINK2 — GWAS 与群体遗传学

概述

PLINK2 是 PLINK 1.9 的高性能后继版本，专为大型队列的全基因组关联研究（GWAS）和群体遗传学分析而设计。它可处理 PLINK 二进制格式（.bed/.bim/.fam）、VCF 和 BGEN 格式的基因型数据，并执行样本和变异质控（QC）、亲缘关系估算、主成分分析（PCA）以及线性／逻辑回归关联检验。得益于多线程和优化的 I/O，PLINK2 在大多数任务上的速度比 PLINK 1.9 快 10–100 倍。其输出文件与下游可视化（Manhattan／QQ 图）和荟萃分析工具兼容。

使用场景

• 在完成基因分型芯片质控后，对病例-对照或定量性状队列运行 GWAS
• 执行样本质控：缺失率、杂合度异常值、性别检查、隐匿亲缘关系
• 计算全基因组 LD 剪枝，用于 PCA 或亲缘关系估算
• 对基因型数据运行 PCA，以识别群体分层
• 在 PLINK 二进制、VCF 和 BGEN 格式之间进行转换
• 在 VCF／填补数据中按 MAF、HWE、缺失率或 INFO 分数筛选变异
• 使用 omics-plotting SKILL 根据关联结果生成 Manhattan 图和 QQ 图
• 对于需要通过混合模型关联分析控制群体结构的生物样本库规模 GWAS（>100k 个样本），改用 regenie 或 SAIGE
• 使用 VCFtools 作为 VCF 专用群体遗传学统计分析的替代工具

先决条件

• 软件：PLINK2（预编译二进制文件；无 pip/conda 包）
• 输入：来自芯片基因分型或基因型填补的 PLINK 二进制文件（.bed/.bim/.fam）或 VCF/BGEN 文件

安装前检查：该工具可能已可用（例如，在 pixi / conda 环境中）。请务必先运行 command -v plink2；如果返回路径，则跳过安装部分。在 pixi 项目中执行工具时，优先使用 pixi run <tool>，而不是直接运行 <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

预先访谈

编写任何分析代码前，请先与用户确认以下事项。

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 取决于 D4，因为这些成分本身就是协变量——应纳入多少个成分，取决于协变量集合，而不是一个独立的参数。完全省略它们，是产生全基因组显著结果、但实际反映的是祖源而非性状的典型原因。

快速开始

# 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

工作流

步骤 1：将 VCF／填补数据转换为 PLINK 二进制格式

将输入的基因型数据转换为 PLINK 二进制格式，以便快速处理。

# 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)"

步骤 2：样本质控 — 缺失率和杂合度

移除缺失率较高或杂合度异常的样本。

# Compute per-sample and per-variant missingness
plink2 \
    --bfile cohort_plink \
    --missing \
    --out qc/sample_missingness \
    --threads 8

# Remove samples with > 2% missingness and variants with > 5% missing
plink2 \
    --bfile cohort_plink \
    --mind 0.02 \
    --geno 0.05 \
    --make-bed \
    --out cohort_sample_qc \
    --threads 8

echo "After sample QC:"
echo "Samples: $(wc -l < cohort_sample_qc.fam)"
echo "Variants: $(wc -l < cohort_sample_qc.bim)"

步骤 3：变异质控 — MAF、HWE 和 INFO 分数筛选

按次要等位基因频率、Hardy-Weinberg 平衡和基因型填补质量筛选变异。

# Variant QC: MAF, HWE, missingness
plink2 \
    --bfile cohort_sample_qc \
    --maf 0.01 \
    --hwe 1e-6 \
    --geno 0.05 \
    --make-bed \
    --out cohort_variantqc \
    --threads 8

echo "After variant QC:"
echo "Variants remaining: $(wc -l < cohort_variantqc.bim)"

# For imputed data: also filter by INFO score (using VCF INFO field)
# plink2 --bfile cohort_variantqc --var-min-qual 0.8 --make-bed --out cohort_info_qc

步骤 4：LD 剪枝和群体分层 PCA

根据经过 LD 剪枝的变异计算主成分。

# LD pruning: remove one variant from each pair with r² > 0.2
plink2 \
    --bfile cohort_variantqc \
    --indep-pairwise 50 5 0.2 \
    --out qc/ld_pruned \
    --threads 8

# Compute PCA on LD-pruned variants (top 20 PCs)
plink2 \
    --bfile cohort_variantqc \
    --extract qc/ld_pruned.prune.in \
    --pca 20 \
    --out pca/cohort_pca \
    --threads 8

echo "PCA files: pca/cohort_pca.eigenvec (PCs) and pca/cohort_pca.eigenval (variance)"
head pca/cohort_pca.eigenvec

步骤 5：运行 GWAS 关联分析

执行逻辑回归（病例-对照）或线性回归（定量性状）。

# Case-control GWAS: logistic regression with covariates
plink2 \
    --bfile cohort_variantqc \
    --pheno phenotypes.txt \
    --1 \
    --covar covariates.txt \
    --covar-variance-standardize \
    --logistic hide-covar \
    --ci 0.95 \
    --out results/gwas_cc \
    --threads 8

# Quantitative trait GWAS: linear regression
plink2 \
    --bfile cohort_variantqc \
    --pheno phenotypes.txt \
    --covar covariates.txt \
    --covar-variance-standardize \
    --linear hide-covar \
    --ci 0.95 \
    --out results/gwas_qt \
    --threads 8

echo "GWAS complete. Files: results/gwas_*.glm.*"
wc -l results/gwas_qt.*.glm.linear

步骤 6：绘制 Manhattan 图和 QQ 图

将关联结果整理为 CHR/BP/P 表，然后阅读 skills/data-visualization/omics-plotting/SKILL.md，并按照其中的 “Manhattan” 和 “QQ plot” 配方处理导出的 CSV（→ figures/manhattan.png、figures/qq.png）。

import pandas as pd

# Load GWAS results (PLINK2 linear output) as a CHR/BP/P table for plotting
df = pd.read_csv("results/gwas_qt.PHENO1.glm.linear", sep="\t",
                 usecols=["#CHROM", "POS", "ID", "P"])
df.columns = ["CHR", "BP", "SNP", "P"]
df = df.dropna(subset=["P"])
df["P"] = pd.to_numeric(df["P"], errors="coerce")
df = df[df["P"] > 0].copy()
df.to_csv("gwas.csv", index=False)

print(f"Variants: {len(df)}")
print(f"Genome-wide significant (p<5e-8): {(df['P'] < 5e-8).sum()}")
# Render with the omics-plotting SKILL (`skills/data-visualization/omics-plotting/SKILL.md`) "Manhattan" and "QQ plot" recipes (input columns CHR, BP, P).

关键参数

参数	默认值	范围/选项	作用
--maf	—	0–0.5	次要等位基因频率阈值；0.01 会移除单例变异
--hwe	—	0–1	Hardy-Weinberg 平衡检验 p 值阈值；通常使用 1e-6
--geno	—	0–1	变异缺失率上限；0.05 = 5%
--mind	—	0–1	样本缺失率上限；0.02 = 2%
--linear / --logistic	—	标志	回归模式：定量性状使用线性回归，病例-对照研究使用逻辑回归
--covar	—	文件路径	协变量文件（FID IID COV1 COV2...）；在此添加主成分
--pca	—	整数	从 LD 剪枝后的数据中计算的主成分数量
--indep-pairwise	—	窗口 步长 r²	LD 剪枝：窗口大小（变异数）、步长、r² 阈值
--threads	1	1–64	CPU 线程数；集群上通常使用 8–16
--ci	—	0–1	效应大小输出的置信区间（95% 置信区间使用 0.95）
--1	关闭	标志	重编码表型：1=对照，2=病例 → 0=对照，1=病例（逻辑回归）

常用操作示例

操作示例 1：亲缘关系 QC，并从每对亲缘样本中移除一个

# Compute kinship coefficients (King-robust estimator)
plink2 \
    --bfile cohort_variantqc \
    --extract qc/ld_pruned.prune.in \
    --king-cutoff 0.0625 \
    --out qc/kinship \
    --threads 8

# Remove related individuals (king-cutoff auto-generates exclusion list)
plink2 \
    --bfile cohort_variantqc \
    --remove qc/kinship.king.cutoff.out.id \
    --make-bed \
    --out cohort_unrelated \
    --threads 8

echo "After relatedness QC: $(wc -l < cohort_unrelated.fam) samples"

操作示例 2：解析 GWAS 结果并报告显著位点

import pandas as pd
import numpy as np

def load_gwas_results(filepath: str, p_threshold: float = 5e-8) -> pd.DataFrame:
    """Load PLINK2 GWAS linear/logistic output and return genome-wide significant hits."""
    df = pd.read_csv(filepath, sep="\t",
                     usecols=["#CHROM", "POS", "ID", "REF", "ALT", "A1", "BETA", "SE", "P"])
    df.columns = ["CHR", "POS", "SNP", "REF", "ALT", "A1", "BETA", "SE", "P"]
    df = df.dropna(subset=["P"])
    df["P"] = pd.to_numeric(df["P"], errors="coerce")
    sig = df[df["P"] < p_threshold].sort_values("P")
    return sig

# Load results for one phenotype
hits = load_gwas_results("results/gwas_qt.PHENO1.glm.linear")
print(f"Genome-wide significant loci: {len(hits)}")
print(hits.head(10)[["CHR", "POS", "SNP", "BETA", "SE", "P"]])
hits.to_csv("gwas_top_hits.tsv", sep="\t", index=False)

预期输出

输出	格式	描述
*.glm.linear	TSV	线性 GWAS 结果：每个变异包含 CHR、POS、ID、BETA、SE、P
*.glm.logistic.hybrid	TSV	逻辑回归 GWAS 结果：每个变异包含 OR、95% CI、P
*.eigenvec	TSV	PCA 载荷：FID IID PC1 PC2 ... PC20
*.eigenval	文本	用于解释 PCA 方差的特征值
*.king.cutoff.out.id	TSV	因亲缘关系 QC 而需移除的样本 ID
*.bed/.bim/.fam	二进制	QC 过滤后的 PLINK 二进制格式

故障排除

问题	原因	解决方案
Error: No samples left after --mind filter	缺失率阈值过严	放宽至 --mind 0.05；检查输入数据质量
逻辑回归未能收敛	变异罕见、病例数量少或协变量共线	使用 --logistic firth 进行 Firth 回归；移除相关协变量
PCA 显示明显的离群群组	混合祖源群体或样本污染	使用 PC 截止值移除离群值；使用 1000 Genomes 参考面板检查祖源
--hwe 移除了过多变异	群体结构导致 HWE 检验结果偏高	仅对对照组应用 HWE 过滤：--hwe 1e-6 ctrl-only
BGEN 导入失败	ref-first/ref-last 标志设置错误	查看填补服务文档；尝试 --bgen file.bgen ref-last
关联检验速度很慢	协变量过多或文件过大	预先筛选到 GWAS 显著区域；使用 --threads 16
性别不匹配警告	X 染色体杂合度超出性别特异性阈值	运行 --check-sex 并移除 F 统计量离群值
大型数据集出现内存错误	RAM 不足以处理 500k+ 个变异	使用 --memory 32000 限制 RAM；按染色体拆分

参考资料

• PLINK2 文档 — 官方命令参考与教程
• PLINK2 GitHub：chrchang/plink-ng — 源代码与发行版本
• Chang CC 等人（2015）“第二代 PLINK：迎接更大、更丰富数据集的挑战” — GigaScience 4:7。DOI:10.1186/s13742-015-0047-8
• GWAS 最佳实践（Marees 等人，2018） — 全面的 QC 与分析指南