PLINK实战指南：GWAS杂合率质控全流程解析与代码实现
在基因组关联分析(GWAS)中，数据质量直接影响研究结果的可靠性。杂合率异常可能暗示样本污染或近亲繁殖等问题，而PLINK作为GWAS分析的瑞士军刀，配合R语言的数据处理能力，可以高效完成这一关键质控步骤。本文将手把手带你完成从SNP筛选到异常样本剔除的全流程，并提供可直接复用的代码模板。
1. 理解杂合率质控的核心逻辑
杂合率指个体中杂合基因型的比例。健康人群的杂合率通常处于特定范围，过高可能表明DNA样本污染（如不同个体混合），过低则可能反映近亲繁殖。GWAS分析中，我们通常剔除杂合率偏离群体平均值±3个标准差的个体。
但直接计算全基因组SNP的杂合率存在一个问题：连锁不平衡(LD)会导致结果偏差。因此需要先筛选出一组近似独立的SNP作为计算基础。这就是--indep-pairwise参数的核心作用——通过滑动窗口法去除高LD区域的SNP，得到相对独立的SNP集合。
关键概念：连锁不平衡(LD)指相邻SNP间非随机关联的现象，会导致统计检验效力下降。质控阶段需要控制LD的影响。
2. 环境准备与数据检查
开始前确保已安装以下工具：
- PLINK 1.9或更高版本
- R语言环境（建议4.0+）
- 基础Linux工具（awk/sed等）
检查数据格式是否符合PLINK二进制文件要求：
ls -l # 应看到.bed/.bim/.fam三组文件 # 如只有ped/map文件需先转换： plink --file yourdata --make-bed --out newname3. 独立SNP筛选：--indep-pairwise详解
执行LD修剪获取独立SNP集合：
plink --bfile yourdata \ --indep-pairwise 50 5 0.2 \ --out pruned参数解析：
50：滑动窗口大小（单位：SNP个数）5：每次滑动的步长（SNP个数）0.2：LD阈值（r² > 0.2的SNP将被剔除）
典型问题排查：
- 窗口太小（如20）可能导致保留过多相关SNP
- 步长过大（如10）可能遗漏局部LD区域
- r²阈值过严（如0.1）会大幅减少可用SNP数量
执行后会生成两个文件：
pruned.prune.in：保留的SNP列表pruned.prune.out：剔除的SNP列表
4. 杂合率计算与可视化
基于独立SNP计算杂合率：
plink --bfile yourdata \ --extract pruned.prune.in \ --het \ --out het_results生成的het_results.het文件包含6列：
- FID 家系ID
- IID 个体ID
- O(HOM) 观测纯合数
- E(HOM) 期望纯合数
- N(NM) 非缺失SNP数
- F 近交系数
用R进行数据可视化和异常值筛选：
# 读取数据 het <- read.table("het_results.het", header=TRUE) # 计算杂合率 het$HET_RATE <- (het$N.NM. - het$O.HOM.)/het$N.NM. # 绘制分布直方图 pdf("heterozygosity_distribution.pdf") hist(het$HET_RATE, breaks=30, col="lightblue", xlab="Heterozygosity Rate", main="Sample Heterozygosity Distribution") abline(v=mean(het$HET_RATE) + 3*sd(het$HET_RATE), col="red", lty=2) abline(v=mean(het$HET_RATE) - 3*sd(het$HET_RATE), col="red", lty=2) dev.off() # 筛选异常样本 het_fail <- subset(het, HET_RATE < mean(HET_RATE)-3*sd(HET_RATE) | HET_RATE > mean(HET_RATE)+3*sd(HET_RATE)) # 保存结果 write.table(het_fail[,1:2], "het_fail_ind.txt", row.names=FALSE, col.names=FALSE, quote=FALSE)5. 剔除异常样本与数据清理
使用PLINK剔除杂合率异常的样本：
plink --bfile yourdata \ --remove het_fail_ind.txt \ --make-bed \ --out cleaned_data验证处理效果：
# 计算处理后数据的杂合率 plink --bfile cleaned_data \ --extract pruned.prune.in \ --het \ --out postqc_het # 比较处理前后异常样本数 wc -l het_fail_ind.txt6. 高级技巧与问题排查
场景1：群体分层影响当研究群体存在明显亚结构时，建议分群体计算杂合率阈值。可以先进行PCA分析确定亚群，再分别计算各亚群的杂合率阈值。
场景2：数据量特别大时的优化对于超大规模数据，可以增加LD修剪的严格程度：
plink --bfile bigdata \ --indep-pairwise 100 10 0.1 \ --out strict_pruned常见错误处理：
错误："No valid SNPs left after filtering"
- 原因：LD阈值过严或初始数据质量差
- 解决：放宽r²阈值（如0.2→0.5）或检查原始数据
警告："Fewer than 100 SNPs used for heterozygosity calculation"
- 原因：保留的独立SNP过少
- 解决：调整
--indep-pairwise参数或检查原始数据密度
7. 完整流程自动化脚本
将上述步骤整合为可重复使用的bash脚本：
#!/bin/bash # GWAS杂合率质控自动化脚本 # 用法：./het_qc.sh input_prefix output_prefix INPUT=$1 OUTPUT=$2 # 步骤1：LD修剪 plink --bfile ${INPUT} \ --indep-pairwise 50 5 0.2 \ --out ${OUTPUT}_pruned # 步骤2：杂合率计算 plink --bfile ${INPUT} \ --extract ${OUTPUT}_pruned.prune.in \ --het \ --out ${OUTPUT}_het # 步骤3：R分析 Rscript -e ' het <- read.table(paste0("${OUTPUT}_het.het"), header=TRUE) het$HET_RATE <- (het$N.NM. - het$O.HOM.)/het$N.NM. fail_samples <- subset(het, HET_RATE < mean(HET_RATE)-3*sd(HET_RATE) | HET_RATE > mean(HET_RATE)+3*sd(HET_RATE)) write.table(fail_samples[,1:2], "${OUTPUT}_het_fail_ind.txt", row.names=FALSE, col.names=FALSE, quote=FALSE)' # 步骤4：剔除异常样本 plink --bfile ${INPUT} \ --remove ${OUTPUT}_het_fail_ind.txt \ --make-bed \ --out ${OUTPUT}_clean echo "质控完成！清洁数据保存在 ${OUTPUT}_clean.*"将此脚本保存为het_qc.sh后，运行：
chmod +x het_qc.sh ./het_qc.sh mydata cleaned_data在实际项目中，我发现将杂合率与样本缺失率分析结合能更有效识别问题样本。例如同时偏离杂合率和缺失率阈值的样本，很可能是实验操作问题导致的低质量样本，应优先剔除。
