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 倍高速です。出力ファイルは、後続の可視化（マンハッタンプロット／QQ プロット）やメタ解析ツールと互換性があります。

使用する場面

• ジェノタイピングアレイの QC 後に、症例対照コホートまたは量的形質コホートで GWAS を実施する
• サンプル QC を実施する：欠損率、ヘテロ接合度の外れ値、性別チェック、潜在的な血縁関係
• PCA または血縁関係の推定に向けて、ゲノムワイドな LD プルーニングを計算する
• 遺伝子型データに対して PCA を実行し、集団階層化を特定する
• PLINK バイナリ、VCF、BGEN 形式間で変換する
• VCF／インピュテーションデータのバリアントを MAF、HWE、欠損率、または INFO スコアでフィルタリングする
• 関連解析の結果からマンハッタンプロットと QQ プロットを作成するには、omics-plotting SKILL を使用する
• 集団構造を制御するための混合モデルによる関連解析が必要なバイオバンク規模の GWAS（サンプル数 >100k）では、代わりに regenie または SAIGE を使用する
• VCF 固有の集団遺伝学統計には、代替ツールとして VCFtools を使用する

前提条件

• ソフトウェア：PLINK2（コンパイル済みバイナリ。pip/conda パッケージはありません）
• 入力：ジェノタイピングアレイまたはインピュテーション由来の PLINK バイナリファイル（.bed/.bim/.fam）または VCF/BGEN

インストール前に確認してください：ツールはすでに利用可能な場合があります（例：pixi / conda 環境内）。まず必ず command -v plink2 を実行し、パスが返された場合はインストール用ブロックをスキップしてください。pixi プロジェクト内でツールを実行する場合は、通常の <tool> より pixi run <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：サンプル QC — 欠損率とヘテロ接合度

欠損率が高いサンプル、またはヘテロ接合度の外れ値を除外します。

# 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：バリアント QC — MAF、HWE、INFO スコアによるフィルタリング

マイナーアレル頻度、ハーディ・ワインベルグ平衡、インピュテーション品質に基づいてバリアントをフィルタリングします。

# 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：マンハッタンプロットと QQ プロットを作成する

関連解析の結果を CHR/BP/P テーブルとして準備し、エクスポートした CSV に対して skills/data-visualization/omics-plotting/SKILL.md を読み、「Manhattan」と「QQ plot」のレシピに従ってください（→ 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...）です。PCはここに追加します
--pca	—	整数	LDプルーニング済みデータから計算するPCの数です
--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を行い、血縁ペアごとに1人を除外する

# 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％信頼区間、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 に緩和し、入力データの品質を確認します
ロジスティック回帰が収束しません	まれなバリアント、症例数の少なさ、または共変量の共線性です	Firth回帰には --logistic 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統計量の外れ値を削除します
大規模データセットでメモリエラーが発生します	50万を超えるバリアントに対してRAMが不足しています	RAMの上限を設定するために --memory 32000 を使用し、染色体ごとに分割します

参考資料

• PLINK2ドキュメント — 公式のコマンドリファレンスとチュートリアル
• PLINK2 GitHub: chrchang/plink-ng — ソースコードとリリース
• Chang CC et al. (2015) "Second-generation PLINK: rising to the challenge of larger and richer datasets" — GigaScience 4:7. DOI:10.1186/s13742-015-0047-8
• GWASのベストプラクティス（Marees et al. 2018） — 包括的なQCと解析のガイド