笔记 GWAS 操作流程5-2：驾驭GEMMA混合模型：从G矩阵构建到群体结构校正

1. GEMMA混合模型实战：从G矩阵到群体结构校正

如果你已经用GEMMA跑过基础线性模型（LM），现在遇到更复杂的亲缘关系或群体结构数据，混合线性模型（LMM）就是你的下一站。我刚开始用LMM时，最大的困惑就是G矩阵怎么构建、PCA结果怎么塞进模型。今天咱们就手把手走通这个流程，顺便聊聊怎么避免群体结构带来的假阳性结果。

先说说为什么需要LMM。当样本间存在亲缘关系（比如家系数据）或群体分层（比如不同地理来源的群体），传统LM会把遗传背景差异误判为性状关联信号。去年我分析一组水稻数据时，LM模型跑出20多个显著位点，加入亲缘关系矩阵后只剩3个是真的——这就是为什么说LMM是GWAS的"防伪标签"。

2. 构建G矩阵：亲缘关系的量化

2.1 标准化方法选择

GEMMA提供四种G矩阵计算方法（-gk参数），实测下来最常用的是下面两种：

# 方法1：中心化矩阵（适合近交系） gemma -bfile genotype -gk 1 -o centered_matrix # 方法2：标准化矩阵（推荐用于异交群体） gemma -bfile genotype -gk 2 -o standardized_matrix

我习惯先用方法2，因为它的数值范围更稳定。曾经用小鼠数据对比过，当群体中存在强近交系时，方法1会过度放大某些亲缘关系。生成的G矩阵保存在output/result.sXX.txt，是个对称矩阵，可以用R可视化：

library(gplots) g_matrix <- as.matrix(read.table("result.sXX.txt")) heatmap.2(g_matrix, trace="none", col=bluered(100))

2.2 矩阵质量检查

三个必看指标：

1. 对角线值：正常范围在0.8-1.2之间，如果出现>2的值可能是样本污染
2. 非对角线分布：健康数据应该呈现右偏分布，我一般用以下代码快速检查：

hist(g_matrix[upper.tri(g_matrix)], breaks=50, main="Off-diagonal Values Distribution")

3. 特征值分解：用eigen()函数计算特征值，第一个特征值占比不应超过30%

3. 群体结构校正：PCA实战技巧

3.1 用plink生成PCA

虽然GEMMA可以直接读G矩阵，但群体结构校正还需要PCA结果。这里推荐用plink2更高效：

plink2 --bfile genotype --pca 20 var-wts --out pca_result

关键参数：

• --pca 20：保留前20个主成分（通常10-20足够）
• var-wts：输出载荷矩阵，方便后续解释

3.2 协变量文件制作

把PCA结果转为GEMMA能读的协变量格式，注意要包含截距项（第一列全1）：

# 添加截距列并保留前5个PC awk 'BEGIN{OFS="\t";print "ID","PC1","PC2","PC3","PC4","PC5"}{print $1,$2,$3,$4,$5,$6}' pca_result.eigenvec > cov.txt # 用R添加截距列更稳妥 pcs <- read.table("pca_result.eigenvec") cov_file <- cbind(rep(1,nrow(pcs)), pcs[,1:5]) write.table(cov_file, "cov.txt", col.names=F, row.names=F)

4. 运行LMM模型：参数详解

4.1 完整命令示例

gemma -bfile genotype \ -k output/result.sXX.txt \ -c cov.txt \ -lmm 4 \ -o lmm_results

参数解析：

• -lmm 4：使用似然比检验（LRT），比默认的Wald检验更稳定
• -miss 0.1：可选项，允许10%的缺失率（默认是5%）

4.2 结果文件解读

主要看output/lmm_results.assoc.txt里的几列：

• beta：效应值大小
• se：标准误
• p_lrt：LRT检验的p值（如果用-lmm 1就是p_wald）
• p_wald：Wald检验的p值

建议用QQ图验证模型校准情况：

res <- read.table("lmm_results.assoc.txt", header=T) lambda <- median(qchisq(1-res$p_lrt,1))/0.454 qqplot(res$p_lrt, main=paste("Lambda=",round(lambda,3)))

5. 校正效果验证：LM vs LMM对比

5.1 假阳性控制

拿我之前的大豆数据为例：

• LM模型：λGC=1.82（严重膨胀）
• LMM模型：λGC=1.02（理想范围）

可以用下面代码快速计算膨胀因子：

calc_lambda <- function(pvals) { chisq <- qchisq(1-pvals,1) median(chisq)/qchisq(0.5,1) }

5.2 信号一致性检查

健康数据中，LM和LMM的top信号应该有一定重叠：

lm_res <- read.table("lm.assoc.txt", header=T) lmm_res <- read.table("lmm.assoc.txt", header=T) merge_res <- merge(lm_res, lmm_res, by="rs") plot(-log10(merge_res$p_wald.x), -log10(merge_res$p_lrt.y), xlab="-log10(LM p)", ylab="-log10(LMM p)") abline(0,1,col="red")

如果发现大量LM显著但LMM不显著的位点，可能需要检查：

1. 群体结构是否校正充分（增加PCA数量）
2. G矩阵计算是否合适（尝试-gk 1或2）
3. 是否存在异常样本（检查G矩阵对角线）

6. 进阶技巧：当标准流程失效时

6.1 稀疏矩阵优化

遇到超大样本（>10万）时，可以启用稀疏矩阵计算：

gemma -bfile bigdata -gk 2 -snp -miss 0.05 -o sparse_matrix

6.2 多性状联合分析

GEMMA支持多性状LMM，只需准备多列表型文件：

gemma -bfile genotype -k matrix -lmm 1 \ -p multi_pheno.txt -n 1 2 3 \ -o multi_trait

其中-n 1 2 3指定使用第1、2、3列表型

6.3 内存控制技巧

对于大型分析，可以通过这些参数避免爆内存：

• -maf 0.01：过滤低频SNP
• -r2 0.1：LD剪枝
• -bsize 1000：分块处理大小

7. 常见报错解决方案

1. Error: eigen decomposition failed

   • 检查G矩阵是否包含NA/Inf
   • 尝试-gk 1或-gk 2重新生成矩阵

2. WARNING: Using only 1 thread

   • 编译时加上-fopenmp选项
   • 设置环境变量export OMP_NUM_THREADS=8

3. P值全部为NA

   • 检查表型数据是否有缺失
   • 确认协变量文件与表型样本顺序一致

最后提醒一个小坑：GEMMA默认会过滤掉MAF<0.01的位点，如果需要保留，务必加上-maf 0参数。我在分析稀有变异时就栽过这个跟头，花了三天才发现是默认过滤导致的信号丢失。