PLINK实战：用--indep-pairwise和R脚本搞定GWAS杂合率质控（附完整代码）

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.2LD阈值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 家系IDIID 个体IDO(HOM) 观测纯合数E(HOM) 期望纯合数N(NM) 非缺失SNP数F 近交系数用R进行数据可视化和异常值筛选# 读取数据 het - read.table(het_results.het, headerTRUE) # 计算杂合率 het$HET_RATE - (het$N.NM. - het$O.HOM.)/het$N.NM. # 绘制分布直方图 pdf(heterozygosity_distribution.pdf) hist(het$HET_RATE, breaks30, collightblue, xlabHeterozygosity Rate, mainSample Heterozygosity Distribution) abline(vmean(het$HET_RATE) 3*sd(het$HET_RATE), colred, lty2) abline(vmean(het$HET_RATE) - 3*sd(het$HET_RATE), colred, lty2) 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.namesFALSE, col.namesFALSE, quoteFALSE)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 # 步骤1LD修剪 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 # 步骤3R分析 Rscript -e het - read.table(paste0(${OUTPUT}_het.het), headerTRUE) 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.namesFALSE, col.namesFALSE, quoteFALSE) # 步骤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在实际项目中我发现将杂合率与样本缺失率分析结合能更有效识别问题样本。例如同时偏离杂合率和缺失率阈值的样本很可能是实验操作问题导致的低质量样本应优先剔除。