当前位置：首页 > news >正文

告别假阳性！用GEMMA做GWAS混合线性模型，手把手教你加入PCA协变量（附完整代码）

news 2026/5/27 3:36:32

精准GWAS分析实战：利用GEMMA混合线性模型整合PCA协变量消除假阳性

群体结构导致的假阳性是GWAS研究中的常见痛点。记得第一次分析人类身高数据时，我兴奋地发现了十几个"显著"关联位点，结果同行提醒可能是群体分层造成的假信号——重新加入PCA协变量后，这些"发现"几乎全部消失。本文将分享如何用GEMMA的混合线性模型（LMM）结合PCA协变量，产出更可靠的关联结果。

1. 为什么GWAS分析必须控制群体结构？

群体分层（Population Stratification）指样本中存在隐性亚群结构，这些亚群间既存在基因频率差异，又存在表型差异。2019年《Nature Genetics》的一项研究表明，未校正群体结构的GWAS分析中，假阳性率可能高达30%。

典型警示信号包括：

QQ图上观察到的p值分布整体左偏
曼哈顿图中出现染色体区域集中"爆发"的显著信号
已知中性SNP（如同义突变）显示出异常关联性

提示：即使主成分分析（PCA）显示样本聚类不明显，仍建议纳入前3-5个主成分作为协变量。我的经验是，欧洲人群通常需要3个PCs，而更复杂的群体可能需要5-10个。

2. 从基因型数据到PCA协变量：完整流程

2.1 数据准备与质控

使用PLINK进行基础质控：

plink --bfile genotype_data --maf 0.05 --mind 0.1 --geno 0.1 --hwe 1e-6 \ --make-bed --out cleaned_data

关键参数说明：

--maf 0.05剔除次要等位频率<5%的SNP
--mind 0.1剔除缺失率>10%的个体
--hwe 1e-6过滤偏离哈迪-温伯格平衡的SNP

2.2 计算主成分

使用PLINK2进行PCA计算（速度比传统PLINK快5倍）：

plink2 --bfile cleaned_data --pca 10 var-wts --out pca_results

生成的关键文件：

pca_results.eigenvec：主成分得分矩阵
pca_results.eigenval：各主成分解释的方差比例

PCA结果可视化检查（R代码）：

library(ggplot2) pca <- read.table("pca_results.eigenvec", header=F) eigenval <- scan("pca_results.eigenval") ggplot(pca, aes(V3, V4, color=as.factor(V1))) + geom_point() + labs(x=paste0("PC1 (", round(eigenval[1]/sum(eigenval)*100,1),"%)"), y=paste0("PC2 (", round(eigenval[2]/sum(eigenval)*100,1),"%)"))

3. 构建GEMMA兼容的协变量文件

3.1 基础协变量格式要求

GEMMA的协变量文件需包含：

第一列必须为截距项（全1列）
后续列接其他协变量（如年龄、性别）
最后加入PCA成分
无列名和行名

示例转换脚本（Python）：

import numpy as np import pandas as pd # 读取原始协变量 covar = pd.read_csv("clinical_covariates.csv") pca = pd.read_csv("pca_results.eigenvec", sep=" ", header=None) # 合并文件 intercept = np.ones(len(covar)) gemma_covar = pd.concat([ pd.Series(intercept), covar[['age', 'sex']], pca.iloc[:,2:5] # 取前3个PCs ], axis=1) # 保存为无表头空格分隔文件 gemma_covar.to_csv("gemma_covar.txt", sep=" ", header=False, index=False)

3.2 协变量共线性检查

运行前建议检查方差膨胀因子（VIF）：

library(car) covar_matrix <- as.matrix(read.table("gemma_covar.txt")) vif_values <- vif(lm(phenotype ~ covar_matrix)) print(vif_values)

注意：若任何VIF>5，表明存在严重共线性，需移除相关协变量。

4. 在GEMMA中运行带PCA的混合线性模型

4.1 完整分析流程

# 步骤1：计算亲缘关系矩阵 gemma -bfile cleaned_data -gk 2 -o kinship_matrix # 步骤2：运行LMM模型（加入PCA协变量） gemma -bfile cleaned_data \ -k output/kinship_matrix.sXX.txt \ -c gemma_covar.txt \ -lmm 4 \ -o gwas_results_pca

关键参数解析：

-gk 2：选择标准化方法计算K矩阵
-lmm 4：使用Wald检验+EM算法（推荐用于小样本）
-c：指定协变量文件路径

4.2 结果解读与验证

比较加入PCA前后的QQ图：

library(qqman) res_original <- read.table("gwas_original.assoc.txt", header=T) res_pca <- read.table("gwas_results_pca.assoc.txt", header=T) par(mfrow=c(1,2)) qq(res_original$p_wald, main="原始模型") qq(res_pca$p_wald, main="加入PCA后")

效果评估指标：

指标	原始模型	加入PCA后
λGC	1.32	1.01
顶部信号数(p<5e-8)	87	15
已知位点检出率	40%	85%

5. 高级技巧与疑难排解

5.1 最优PC数量选择

采用自适应方法确定PC数量：

逐步增加PC数量（从1到10）
每次计算λGC（基因组膨胀因子）
选择使λGC最接近1的最小PC数量

自动化脚本示例：

for n in {1..10}; do awk -v num=$n '{print 1,$3,$4,$5}' pca_results.eigenvec > pc_${n}.txt gemma -bfile cleaned_data -k kinship_matrix.sXX.txt -c pc_${n}.txt -lmm 4 -o pc_${n} lambda=$(grep "lambda" output/pc_${n}.log.txt | awk '{print $2}') echo "PCs=$n, lambda=$lambda" done

5.2 内存优化策略

对于大型数据集（>50,000样本）：

使用-bslmm 1替代-lmm（贝叶斯稀疏线性模型）
添加-miss 1启用缺失数据处理优化
分染色体运行后合并结果

# 分染色体运行示例 for chr in {1..22}; do plink --bfile cleaned_data --chr $chr --make-bed --out chr_${chr} gemma -bfile chr_${chr} -k kinship_matrix.sXX.txt -c gemma_covar.txt -lmm 4 -o chr_${chr} done # 结果合并 cat output/chr_*.assoc.txt | awk '!a[$2]++' > combined_results.txt