全基因组关联分析GWAS实战¶

一句话概述¶

GWAS（Genome-Wide Association Study）通过在全基因组水平检测遗传变异与表型性状的关联，揭示复杂性状的遗传基础，本教程覆盖PLINK/GCTA工具链的完整分析流程。

目录¶

1. 核心知识点
2. 各步骤详解
3. 实战命令/代码
4. 面试常问点
5. 易错/易混淆点
6. 补充知识

1. 核心知识点¶

序号	步骤	工具	目的	输入	输出	关键参数
1	数据格式转换	PLINK2	统一基因型数据格式	VCF/BED/PED	PLINK binary (bed/bim/fam)	`--make-bed`
2	样本质控	PLINK2	去除低质量样本	bed/bim/fam	过滤后数据集	`--mind 0.05`
3	变异位点质控	PLINK2	去除低质量SNP	bed/bim/fam	过滤后数据集	`--geno 0.02 --maf 0.01 --hwe 1e-6`
4	LD修剪	PLINK2	去除高LD冗余位点	过滤后数据集	独立SNP集合	`--indep-pairwise 50 5 0.2`
5	群体结构校正	PLINK2/GCTA	计算PCA控制群体分层	独立SNP集合	PC协变量文件	`--pca 20`
6	亲缘关系检测	PLINK2/KING	检测并去除近亲	独立SNP集合	亲缘关系矩阵	`--king-cutoff 0.0884`
7	关联分析	PLINK2/GCTA	检测SNP-表型关联	基因型+表型+协变量	关联结果汇总统计	`--glm`/`--mlma`
8	多重检验校正	R/Python	控制假阳性率	P值列表	显著位点	Bonferroni/FDR
9	结果可视化	R (qqman)	Manhattan图/QQ图	关联结果	出版级图表	`manhattan()`/`qq()`
10	后续分析	FUMA/MAGMA	功能注释与通路分析	显著SNP列表	基因/通路富集	基因集/LD窗口

2. 各步骤详解¶

步骤1：数据格式转换¶

白话解释：GWAS数据可能来自不同的基因分型平台（如Illumina芯片、WGS变异检测），格式各异（VCF、PLINK PED、Oxford GEN等）。第一步是将所有数据统一转换为PLINK的二进制格式（.bed/.bim/.fam），这是后续所有质控和分析步骤的标准输入格式。

技术细节：

PLINK二进制格式包含三个文件： - .bed：二进制基因型数据（每个基因型用2bit编码） - .bim：变异位点信息（染色体、SNP ID、遗传距离、物理位置、等位基因） - .fam：样本信息（家系ID、个体ID、父母ID、性别、表型）

PLINK2（plink2）是PLINK1.9的升级版，支持多等位基因变异、剂量数据和更大的数据集。对于现代GWAS数据（动辄数十万样本、数百万SNP），PLINK2的内存效率和速度优势明显。

当输入为VCF格式时，需要注意：（1）多等位基因位点的处理（通常拆分为双等位基因）；（2）参考等位基因的一致性；（3）样本表型信息需要单独提供（VCF通常不包含表型）。

# VCF转PLINK格式
plink2 --vcf genotypes.vcf.gz \
  --make-bed \
  --out study_data \
  --allow-extra-chr \
  --max-alleles 2

# 更新表型信息
plink2 --bfile study_data \
  --pheno phenotype.txt \
  --make-bed \
  --out study_data_pheno

步骤2：样本质控¶

白话解释：有些样本可能因为DNA质量差、实验操作失误等原因导致大量位点缺失（call rate低），这些样本如果不去除会影响整体分析质量。同时还需要检查性别一致性、杂合度异常等问题。

技术细节：

样本质控主要包括以下几个方面：

样本缺失率（missingness per sample）：通常要求每个样本的缺失率<5%（--mind 0.05）。缺失率高的样本通常意味着DNA质量差或实验出了问题。
性别检查：通过X染色体的杂合度推断遗传性别，与报告性别比较。F值（近交系数）：女性通常<0.2，男性通常>0.8。不一致的样本可能存在样本混淆。
杂合度检查：全基因组杂合度异常（过高或过低）可能提示DNA污染（过高）或近交（过低）。通常排除偏离均值±3SD的样本。
种族/族群确认：与参考面板（如1000 Genomes）进行PCA，确认样本的遗传背景与自报族群一致。

# 计算样本缺失率
plink2 --bfile study_data_pheno \
  --missing \
  --out qc_stats

# 过滤高缺失率样本
plink2 --bfile study_data_pheno \
  --mind 0.05 \
  --make-bed \
  --out study_qc1

# 性别检查
plink --bfile study_qc1 \
  --check-sex \
  --out sex_check

# 杂合度检查
plink --bfile study_qc1 \
  --het \
  --out het_check

步骤3：变异位点质控¶

白话解释：不是所有SNP都适合做关联分析。缺失率太高的SNP信息不可靠、频率太低的SNP统计功效不足、严重偏离Hardy-Weinberg平衡的SNP可能是基因分型错误——这些都需要过滤掉。

技术细节：

SNP质控的标准阈值：

SNP缺失率（--geno）：通常<2%或<5%。这个阈值应该比样本缺失率更严格，因为一个坏的SNP会影响所有样本。
最小等位基因频率（MAF, --maf）：通常>1%或>5%。低频变异的统计功效很低（需要极大的样本量），而且基因分型错误率相对更高。
Hardy-Weinberg平衡检验（--hwe）：在对照组中检验，P<1e-6通常提示基因分型错误。注意：(a) 只在对照组中检测HWE，因为case组中真正的风险位点可能偏离HWE；(b) 对于X染色体，男性是半合子，需要分别处理。
差异缺失检验：case和control之间的缺失率差异可能产生虚假关联（--test-missing）。

# 综合SNP质控
plink2 --bfile study_qc1 \
  --geno 0.02 \
  --maf 0.01 \
  --hwe 1e-6 \
  --make-bed \
  --out study_qc2

# 差异缺失检验（PLINK1.9语法）
plink --bfile study_qc1 \
  --test-missing \
  --out diff_missing

# 排除差异缺失显著的SNP
awk '$5 < 1e-5 {print $2}' diff_missing.missing > fail_diffmiss.txt
plink2 --bfile study_qc2 \
  --exclude fail_diffmiss.txt \
  --make-bed \
  --out study_qc3

步骤4：LD修剪（LD Pruning）¶

白话解释：基因组上相邻的SNP往往是高度相关的（连锁不平衡，LD），在做PCA和亲缘关系估计时，如果不去除这种冗余，某些基因组区域会被过度代表。LD修剪就是保留一组相对独立的SNP子集。

技术细节：

LD修剪使用滑窗法：在每个窗口内，计算所有SNP对的r²值，如果某对SNP的r²超过阈值，则移除MAF较低的那个。

--indep-pairwise 50 5 0.2 含义： - 50：窗口大小（50个SNP） - 5：步长（每次滑动5个SNP） - 0.2：r²阈值（超过0.2的SNP对需要修剪）

注意：LD修剪后的SNP集合仅用于PCA和亲缘关系估计，关联分析本身使用所有通过质控的SNP。

特殊区域处理：MHC区域（chr6:25-35Mb）、染色体倒位区域等长程LD区域应额外排除，否则它们会主导PCA的前几个成分。

# LD修剪
plink2 --bfile study_qc3 \
  --indep-pairwise 50 5 0.2 \
  --out ld_pruned

# 排除长程LD区域
# long_range_ld.txt包含已知的长程LD区域坐标
plink2 --bfile study_qc3 \
  --extract ld_pruned.prune.in \
  --exclude range long_range_ld.txt \
  --make-bed \
  --out study_pruned

步骤5：群体结构校正（PCA）¶

白话解释：如果研究样本来自不同的遗传背景（如欧洲人和亚洲人混合），那么任何在两个群体间频率不同的SNP都可能显示出与表型的"虚假关联"。PCA能捕捉群体结构，将其作为协变量纳入关联模型来消除这种混杂。

技术细节：

主成分分析（PCA）将高维基因型数据投影到低维空间，前几个主成分（PC）通常反映群体结构。常见做法是取前10-20个PC作为协变量。

选择多少个PC纳入模型需要权衡： - 太少：群体分层未充分校正，膨胀因子（λ）偏高 - 太多：可能过度校正，降低统计功效

判断标准： - 观察scree plot（碎石图），选择拐点之前的PC数 - 检查QQ图的膨胀程度（λ应接近1.0） - Tracy-Widom检验判断哪些PC是显著的

GCTA的PCA实现（--pca）适用于大规模数据，FlashPCA2是另一个高效的选择。

# PLINK2计算PCA
plink2 --bfile study_pruned \
  --pca 20 \
  --out pca_results

# 或使用GCTA
gcta64 --bfile study_pruned \
  --make-grm \
  --out grm_data

gcta64 --grm grm_data \
  --pca 20 \
  --out pca_gcta

步骤6：亲缘关系检测¶

白话解释：如果研究中有近亲（如双胞胎、兄弟姐妹），他们之间共享的基因型比随机个体多得多，这会导致统计检验的假设被违反。我们需要找到并去除近亲关系的个体对中的一个。

技术细节：

KING算法是检测亲缘关系的标准方法，使用IBS（Identity-By-State）信息估计kinship系数。常用阈值： - 一级亲属（父母-子女、全同胞）：kinship > 0.177 - 二级亲属（祖孙、半同胞、叔侄）：0.0884 < kinship ≤ 0.177 - 三级亲属：0.0442 < kinship ≤ 0.0884

通常使用0.0884作为阈值（排除二级及以上亲属）。在需要去除一个个体时，通常保留缺失率更低或表型信息更完整的那个。

PLINK2的--king-cutoff可以自动执行这一步，选择要移除的最小样本集使得没有一对样本超过阈值。

# KING亲缘关系检测
plink2 --bfile study_pruned \
  --king-cutoff 0.0884 \
  --out king_filter

# 去除近亲样本
plink2 --bfile study_qc3 \
  --remove king_filter.king.cutoff.out.id \
  --make-bed \
  --out study_unrelated

步骤7：关联分析¶

白话解释：这是GWAS的核心步骤——对每个SNP逐一检验其基因型与表型之间是否存在统计学上的显著关联。对于连续性状用线性回归，对于二分类性状（如疾病有/无）用逻辑回归。

技术细节：

PLINK2的--glm命令执行广义线性模型（GLM）关联分析，支持： - 连续性状：线性回归 - 二分类性状：逻辑回归/Firth回归 - 加性模型（默认）、显性模型、隐性模型

模型：Y = β₀ + β₁*SNP + β₂*PC1 + β₃*PC2 + ... + βₖ*Covariates + ε

GCTA的混合线性模型（MLM/MLMA）可以更好地控制群体结构和隐性亲缘关系： - --mlma：将GRM作为随机效应纳入模型 - 比固定效应PCA方法更能控制虚假关联 - BOLT-LMM和SAIGE是MLM的更高效实现

对于大规模GWAS（>100k样本），建议使用BOLT-LMM（连续性状）或SAIGE（二分类性状，特别是case-control比例不平衡时）。

# PLINK2 GLM关联分析（连续性状）
plink2 --bfile study_unrelated \
  --pheno phenotype.txt \
  --covar pca_results.eigenvec \
  --covar-col-nums 3-12 \
  --glm hide-covar \
  --out gwas_results

# PLINK2 GLM关联分析（二分类性状，Firth回归）
plink2 --bfile study_unrelated \
  --pheno phenotype_binary.txt \
  --1 \
  --covar pca_results.eigenvec \
  --covar-col-nums 3-12 \
  --glm firth-fallback hide-covar \
  --out gwas_binary

# GCTA MLMA（混合线性模型）
gcta64 --bfile study_unrelated \
  --grm grm_data \
  --pheno phenotype.txt \
  --qcovar pca_results.eigenvec \
  --mlma \
  --out gwas_mlma \
  --thread-num 16

步骤8：多重检验校正¶

白话解释：GWAS检验了数百万个SNP，即使完全随机，也会有很多SNP看起来"显著"。我们需要严格的多重检验校正来控制假阳性率。

技术细节：

GWAS中最常用的显著性阈值是5×10⁻⁸，这来源于Bonferroni校正：假设全基因组有约100万个独立检验（考虑LD后），0.05/1,000,000 ≈ 5×10⁻⁸。

另外还有一个"建议性显著"阈值：1×10⁻⁵，用于标记值得在独立队列中验证的位点。

基因组通胀因子（λ, genomic inflation factor）： - 计算方式：观察到的中位数χ²统计量 / 期望中位数（0.456） - λ ≈ 1.0：群体结构控制良好 - λ > 1.1：可能存在未校正的群体分层或其他系统偏差 - LD Score Regression (LDSC) 可以区分真正的多基因信号和技术偏差导致的膨胀

# R语言中计算基因组通胀因子
results <- read.table("gwas_results.PHENO1.glm.linear", header=TRUE)
chisq <- qchisq(1 - results$P, 1)
lambda <- median(chisq) / qchisq(0.5, 1)
cat("Genomic inflation factor λ =", lambda, "\n")

# 如果λ > 1，进行基因组控制校正
results$P_GC <- pchisq(chisq / lambda, 1, lower.tail=FALSE)

步骤9：结果可视化¶

白话解释：Manhattan图展示全基因组每个SNP的关联强度（-log10 P值），就像纽约曼哈顿的天际线，突出的"高楼"就是显著关联的区域。QQ图则用于评估整体分析质量——点应沿对角线分布，末端的偏离代表真阳性信号。

技术细节：

Manhattan图的解读： - X轴：染色体位置（按染色体编号排列） - Y轴：-log10(P value) - 红线：全基因组显著阈值（5×10⁻⁸，-log10 = 7.3） - 蓝线：建议性显著阈值（1×10⁻⁵，-log10 = 5） - 显著的"峰"代表与表型关联的基因组区域（locus）

QQ图的解读： - X轴：期望-log10(P)（均匀分布） - Y轴：观察到的-log10(P) - 大部分点应沿对角线（y=x）分布 - 全面性的偏离提示群体分层或系统偏差 - 尾部的偏离（大P值区域）代表真正的关联信号

LocusZoom图：放大单个显著locus，展示lead SNP和周围SNP的P值及LD关系。

# 安装和加载包
# install.packages("qqman")
library(qqman)

# 读取关联结果
gwas <- read.table("gwas_results.PHENO1.glm.linear", header=TRUE)
gwas_plot <- data.frame(
  SNP = gwas$ID,
  CHR = gwas$`#CHROM`,
  BP = gwas$POS,
  P = gwas$P
)

# Manhattan图
png("manhattan_plot.png", width=12, height=6, units="in", res=300)
manhattan(gwas_plot,
  chr = "CHR", bp = "BP", snp = "SNP", p = "P",
  main = "GWAS Manhattan Plot",
  col = c("steelblue", "coral"),
  suggestiveline = -log10(1e-5),
  genomewideline = -log10(5e-8),
  cex = 0.6,
  cex.axis = 0.8
)
dev.off()

# QQ图
png("qq_plot.png", width=6, height=6, units="in", res=300)
qq(gwas_plot$P, main="GWAS QQ Plot")
dev.off()

步骤10：后续分析（Post-GWAS）¶

白话解释：找到显著的SNP只是第一步。接下来需要弄清楚这些SNP在哪些基因附近、影响什么生物学通路、在什么组织中有调控作用——这就是Post-GWAS分析。

技术细节：

主要的Post-GWAS分析包括：

基因映射（Gene Mapping）：将显著SNP映射到基因（位置映射、eQTL映射、染色质相互作用映射）
基因集/通路分析：MAGMA、PASCAL等工具将SNP级别的信号汇聚到基因和通路水平
组织/细胞类型富集：使用GTEx数据评估哪些组织的基因表达受关联信号调控
遗传力估计：LDSC估计SNP遗传力和遗传相关性
精细定位（Fine-mapping）：在每个locus中确定最可能的因果变异（FINEMAP、SuSiE）
多基因风险评分（PRS）：将GWAS结果应用于个体风险预测

FUMA（Functional Mapping and Annotation）是一个集成平台，自动化完成大部分Post-GWAS分析。

# MAGMA基因分析
# 1. 注释SNP到基因
magma --annotate \
  --snp-loc gwas_results.snploc \
  --gene-loc NCBI38.gene.loc \
  --out magma_annot

# 2. 基因分析
magma --bfile ref_panel \
  --pval gwas_results.pval ncol=1 \
  --gene-annot magma_annot.genes.annot \
  --out magma_genes

# 3. 基因集分析
magma --gene-results magma_genes.genes.raw \
  --set-annot MSigDB.sets.txt \
  --out magma_pathways

# LD Score Regression遗传力估计
ldsc.py --h2 gwas_sumstats.sumstats.gz \
  --ref-ld-chr eur_w_ld_chr/ \
  --w-ld-chr eur_w_ld_chr/ \
  --out h2_estimate

3. 实战命令/代码（完整流程）¶

#!/bin/bash
# =============================================
# GWAS完整分析流程
# 假设：二分类性状（case-control研究）
# 输入：VCF基因型数据 + 表型文件
# =============================================

set -euo pipefail

# ---- 全局参数 ----
PLINK2="plink2"
GCTA="gcta64"
THREADS=16
INPUT_VCF="raw_genotypes.vcf.gz"
PHENO="phenotype.txt"         # FID IID PHENO (1=control, 2=case)
COVAR="covariates.txt"        # FID IID AGE SEX ...
PREFIX="gwas_study"
OUTDIR="results"
mkdir -p ${OUTDIR}

# ---- Step 1: 格式转换 ----
echo "=== Step 1: 格式转换 ==="
${PLINK2} --vcf ${INPUT_VCF} \
  --make-bed \
  --max-alleles 2 \
  --pheno ${PHENO} \
  --out ${OUTDIR}/${PREFIX}_raw

# ---- Step 2: 样本质控 ----
echo "=== Step 2: 样本质控 ==="
# 计算缺失率
${PLINK2} --bfile ${OUTDIR}/${PREFIX}_raw \
  --missing \
  --out ${OUTDIR}/${PREFIX}_miss

# 过滤：缺失率>5%的样本
${PLINK2} --bfile ${OUTDIR}/${PREFIX}_raw \
  --mind 0.05 \
  --make-bed \
  --out ${OUTDIR}/${PREFIX}_s1

# 性别检查（PLINK1.9）
plink --bfile ${OUTDIR}/${PREFIX}_s1 \
  --check-sex \
  --out ${OUTDIR}/${PREFIX}_sexcheck

# 杂合度检查
plink --bfile ${OUTDIR}/${PREFIX}_s1 \
  --het \
  --out ${OUTDIR}/${PREFIX}_het

# ---- Step 3: SNP质控 ----
echo "=== Step 3: SNP质控 ==="
${PLINK2} --bfile ${OUTDIR}/${PREFIX}_s1 \
  --geno 0.02 \
  --maf 0.01 \
  --hwe 1e-6 \
  --make-bed \
  --out ${OUTDIR}/${PREFIX}_v1

echo "质控后SNP数: $(wc -l < ${OUTDIR}/${PREFIX}_v1.bim)"

# ---- Step 4: LD修剪 ----
echo "=== Step 4: LD修剪 ==="
${PLINK2} --bfile ${OUTDIR}/${PREFIX}_v1 \
  --indep-pairwise 50 5 0.2 \
  --out ${OUTDIR}/${PREFIX}_ldprune

${PLINK2} --bfile ${OUTDIR}/${PREFIX}_v1 \
  --extract ${OUTDIR}/${PREFIX}_ldprune.prune.in \
  --make-bed \
  --out ${OUTDIR}/${PREFIX}_pruned

echo "修剪后独立SNP数: $(wc -l < ${OUTDIR}/${PREFIX}_pruned.bim)"

# ---- Step 5: PCA ----
echo "=== Step 5: PCA ==="
${PLINK2} --bfile ${OUTDIR}/${PREFIX}_pruned \
  --pca 20 \
  --out ${OUTDIR}/${PREFIX}_pca

# ---- Step 6: 亲缘关系检测 ----
echo "=== Step 6: 亲缘关系检测 ==="
${PLINK2} --bfile ${OUTDIR}/${PREFIX}_pruned \
  --king-cutoff 0.0884 \
  --out ${OUTDIR}/${PREFIX}_king

# 去除近亲
if [ -f ${OUTDIR}/${PREFIX}_king.king.cutoff.out.id ]; then
  NRELATED=$(wc -l < ${OUTDIR}/${PREFIX}_king.king.cutoff.out.id)
  echo "去除 ${NRELATED} 个近亲样本"
  ${PLINK2} --bfile ${OUTDIR}/${PREFIX}_v1 \
    --remove ${OUTDIR}/${PREFIX}_king.king.cutoff.out.id \
    --make-bed \
    --out ${OUTDIR}/${PREFIX}_final
else
  echo "未检测到近亲"
  cp ${OUTDIR}/${PREFIX}_v1.bed ${OUTDIR}/${PREFIX}_final.bed
  cp ${OUTDIR}/${PREFIX}_v1.bim ${OUTDIR}/${PREFIX}_final.bim
  cp ${OUTDIR}/${PREFIX}_v1.fam ${OUTDIR}/${PREFIX}_final.fam
fi

# ---- Step 7: 关联分析 ----
echo "=== Step 7: 关联分析 ==="
${PLINK2} --bfile ${OUTDIR}/${PREFIX}_final \
  --pheno ${PHENO} \
  --covar ${OUTDIR}/${PREFIX}_pca.eigenvec \
  --covar-col-nums 3-12 \
  --glm firth-fallback hide-covar cols=chrom,pos,ref,alt,a1freq,firth,test,nobs,orbeta,se,ci,tz,p \
  --out ${OUTDIR}/${PREFIX}_gwas

echo "=== GWAS完成 ==="
echo "结果文件: ${OUTDIR}/${PREFIX}_gwas.PHENO1.glm.logistic.hybrid"

#!/usr/bin/env Rscript
# =============================================
# GWAS结果可视化与后处理
# =============================================

library(qqman)
library(dplyr)
library(ggplot2)

# ---- 读取结果 ----
gwas <- read.table("results/gwas_study_gwas.PHENO1.glm.logistic.hybrid",
                   header=TRUE, comment.char="")
colnames(gwas)[1] <- "CHR"

# 准备qqman格式
gwas_clean <- gwas %>%
  filter(!is.na(P)) %>%
  select(SNP=ID, CHR, BP=POS, P) %>%
  filter(CHR %in% 1:22)

gwas_clean$CHR <- as.numeric(gwas_clean$CHR)

# ---- 基因组通胀因子 ----
chisq <- qchisq(1 - gwas_clean$P, 1)
lambda <- median(chisq, na.rm=TRUE) / qchisq(0.5, 1)
cat(sprintf("Genomic inflation factor: λ = %.4f\n", lambda))

# ---- Manhattan图 ----
png("results/manhattan.png", width=14, height=6, units="in", res=300)
manhattan(gwas_clean,
  main = paste0("GWAS Manhattan Plot (λ = ", round(lambda, 3), ")"),
  col = c("#1B9E77", "#D95F02"),
  suggestiveline = -log10(1e-5),
  genomewideline = -log10(5e-8),
  cex = 0.5,
  ylim = c(0, max(-log10(gwas_clean$P), na.rm=TRUE) + 1)
)
dev.off()

# ---- QQ图 ----
png("results/qq_plot.png", width=7, height=7, units="in", res=300)
qq(gwas_clean$P, main = paste0("QQ Plot (λ = ", round(lambda, 3), ")"))
dev.off()

# ---- 显著位点汇总 ----
sig_snps <- gwas_clean %>%
  filter(P < 5e-8) %>%
  arrange(P)

cat(sprintf("\n全基因组显著位点数: %d\n", nrow(sig_snps)))
cat(sprintf("建议性显著位点数: %d\n", sum(gwas_clean$P < 1e-5)))

write.table(sig_snps, "results/significant_snps.txt",
            row.names=FALSE, quote=FALSE, sep="\t")

# ---- 独立位点识别（LD clumping） ----
# 需要回到PLINK执行clumping
system(paste0(
  "plink2 --bfile results/gwas_study_final ",
  "--clump results/gwas_study_gwas.PHENO1.glm.logistic.hybrid ",
  "--clump-p1 5e-8 --clump-p2 1e-5 --clump-r2 0.1 --clump-kb 1000 ",
  "--out results/gwas_clumped"
))

4. 面试常问点¶

Q1：什么是GWAS？其基本原理是什么？

GWAS是全基因组关联研究，通过在全基因组水平上系统性地检测数十万到数百万个遗传变异（通常是SNP）与表型性状之间的统计关联。基本原理基于"常见疾病-常见变异"假说：复杂性状受到多个常见遗传变异的共同影响，每个变异的效应较小。通过在大样本中检测每个SNP的等位基因频率在case和control之间是否有显著差异（或与连续性状的相关性），来定位与疾病/性状相关的基因组区域。

Q2：GWAS的全基因组显著性阈值为什么是5×10⁻⁸？

这个阈值来源于Bonferroni校正。虽然现代芯片可以检测数百万个SNP，但由于连锁不平衡（LD），独立统计检验的数目约为100万个。因此Bonferroni阈值为0.05/1,000,000 = 5×10⁻⁸。这是一个保守的阈值，在不同的基因组/人群中，真正独立检验的数目可能不同，但5×10⁻⁸已成为领域共识。

Q3：什么是群体分层（population stratification）？如何校正？

群体分层指研究样本中包含遗传背景不同的亚群体，这些亚群体之间的等位基因频率差异可能与表型差异混杂，导致虚假关联。校正方法：（1）PCA：计算前几个PC作为协变量；（2）混合线性模型（MLM）：将遗传关系矩阵作为随机效应；（3）基因组控制（GC）：用λ校正所有P值（最粗糙的方法）。

Q4：Manhattan图和QQ图分别告诉我们什么？

Manhattan图直观展示全基因组每个SNP的关联信号强度，横轴是染色体位置，纵轴是-log10(P)。突出的"峰"代表与表型关联的区域。QQ图将观察到的P值与理论期望进行比较：大部分点应沿对角线分布（代表无效位点），末端偏离代表真阳性信号。QQ图全面偏离对角线提示存在系统性偏差（如未校正的群体分层），通过λ值量化。

Q5：LD clumping和LD pruning的区别是什么？

LD pruning（修剪）：以滑窗方式迭代去除高LD的SNP，不考虑关联信号的强弱，目的是得到一组独立SNP集合（用于PCA等）。LD clumping（聚类）：以关联P值最小的SNP为seed，将其LD block内的其他SNP归入同一组，目的是识别独立的关联信号（用于报告GWAS结果和PRS构建）。

Q6：PLINK的GLM和GCTA的MLMA有什么区别？什么情况下选择哪种？

PLINK GLM使用固定效应模型（PC作为协变量），计算速度快，适合样本量不太大且群体结构简单的情况。GCTA MLMA使用混合线性模型（GRM作为随机效应），能更好地控制隐性群体结构和远亲关系，但计算量大。对于生物样本库规模的数据（>100k），BOLT-LMM和SAIGE是更高效的MLM实现。

Q7：什么是多基因风险评分（PRS）？如何构建？

PRS将多个SNP的效应累加为一个综合评分，预测个体的遗传风险。构建步骤：（1）从GWAS获取效应量（β或OR）；（2）选择纳入PRS的SNP（可用不同P值阈值或贝叶斯方法如PRScs）；（3）考虑LD（clumping或LDpred）；（4）在独立验证队列中评估预测性能。PRS = Σ(βᵢ × dosageᵢ)

Q8：Hardy-Weinberg平衡检验在GWAS质控中的作用是什么？为什么只在对照组中检验？

HWE检验用于识别基因分型错误：在无进化力作用的随机交配群体中，基因型频率应符合p² + 2pq + q² = 1。严重偏离HWE的SNP很可能存在技术问题。只在对照组中检验是因为：真正与疾病相关的位点在case组中可能偏离HWE（如显性效应），如果在全样本中检验可能会错误地去除真阳性信号。

5. 易错/易混淆点¶

混淆LD pruning和LD clumping的应用场景：Pruning用于PCA/亲缘关系估计（需要独立SNP集），clumping用于识别独立关联信号和PRS。不能将pruning后的SNP集直接用于报告GWAS hits。
HWE阈值设置错误：在case组中使用过严格的HWE过滤可能去除真正的风险位点。标准做法是仅在control组中以1e-6为阈值过滤。
忽略X/Y染色体和线粒体的特殊处理：男性X染色体是半合子（只有一个拷贝），HWE检验和关联分析的模型不同。Y染色体和线粒体仅在特定样本中有意义。
表型编码错误：PLINK中case/control表型：1=control, 2=case（不是0/1）。连续表型不需要特殊编码，但缺失值应标记为-9或NA。编码错误会完全颠倒OR方向。
PC数目选择不当：使用太少的PC无法充分控制群体分层（λ偏高），使用太多PC可能导致过度校正，降低真阳性的检测能力。应通过scree plot和λ值综合判断。
忽略批次效应：如果样本在不同批次分型，且批次与case/control状态相关，会产生大量假阳性。应将批次作为协变量或在同一批次内进行分析。
对罕见变异使用标准GWAS方法：MAF<1%的变异不适合单SNP关联分析（统计功效极低）。应使用burden test、SKAT等聚合方法（如RVTESTS、SAIGE-GENE）。
报告结果时忽略效应等位基因方向：不同参考面板的效应等位基因和非效应等位基因可能不同。进行meta-analysis或PRS构建时，必须对齐效应等位基因方向。

6. 补充知识¶

GWAS的历史与发展¶

2005年：第一个成功的GWAS（年龄相关性黄斑变性，CFH基因）
2007-2010年：GWAS黄金期，大量疾病和性状的关联位点被发现
2010年代：样本量从千人级扩展到百万级（UK Biobank等）
当前趋势：多族群GWAS、罕见变异关联、全基因组测序（WGS-based GWAS）

关键数据库与资源¶

GWAS Catalog (https://www.ebi.ac.uk/gwas/)：所有已发表GWAS结果的集合
UK Biobank：50万人的基因型和深度表型数据
1000 Genomes：全球人群参考面板
gnomAD：等位基因频率参考数据库
FUMA (https://fuma.ctglab.nl/)：GWAS后分析在线平台
LD Score Regression (https://github.com/bulik/ldsc)：遗传力和遗传相关性

样本量与统计功效¶

GWAS需要大样本量是因为单个常见变异的效应量很小（通常OR 1.1-1.5）。检测功效取决于： - 样本量（N）：最重要的因素 - 效应量（β/OR）：越大越容易检测 - 等位基因频率（MAF）：中等频率的变异更容易检测 - 多重检验负担：更多的检验需要更大的样本量

经验法则：要以80%的功效检测OR=1.2、MAF=0.2的位点，需要约5000 case + 5000 control。

新兴方法¶

SAIGE：适用于不平衡case-control比例的混合模型
REGENIE：适用于大规模生物样本库的两步关联方法
Multi-ancestry GWAS：联合多个族群增加发现能力
Mendelian Randomization：利用GWAS结果推断因果关系
TWAS/PWAS：整合表达/蛋白质组学数据的关联分析

全基因组关联分析GWAS实战¶

一句话概述¶

目录¶

1. 核心知识点¶

2. 各步骤详解¶

步骤1：数据格式转换¶

步骤2：样本质控¶

步骤3：变异位点质控¶

步骤4：LD修剪（LD Pruning）¶

步骤5：群体结构校正（PCA）¶

步骤6：亲缘关系检测¶

步骤7：关联分析¶

步骤8：多重检验校正¶

步骤9：结果可视化¶

步骤10：后续分析（Post-GWAS）¶

3. 实战命令/代码（完整流程）¶

4. 面试常问点¶

5. 易错/易混淆点¶

6. 补充知识¶

GWAS的历史与发展¶

关键数据库与资源¶

样本量与统计功效¶

新兴方法¶

📚 相关文章推荐