1. 项目概述与核心需求解析做GWAS全基因组关联分析的人都知道基因型数据拿到手之后质控环节是决定后续分析成败的关键。这个项目标题“笔记 GWAS 操作流程2-4哈温平衡检验”一看就是系列笔记中的一环前面应该已经完成了样本层面的质控和位点层面的缺失率过滤现在轮到哈温平衡检验Hardy-Weinberg Equilibrium简称HWE检验了。哈温平衡检验在GWAS流程里扮演什么角色简单说它是一道“筛子”用来判断每个SNP位点的基因型频率是否符合群体遗传学预期。如果一个位点严重偏离HWE可能意味着几种情况基因分型错误、拷贝数变异干扰、群体分层、自然选择作用或者就是测序/芯片本身的系统误差。对于常见病关联研究通常会把严重偏离HWE的位点剔除避免假阳性关联信号。这个内容适合谁看如果你正在跑GWAS流程手头有PLINK格式的基因型数据已经完成了初步质控接下来需要做HWE过滤那这篇笔记就是为你准备的。如果你用的是其他工具链比如GCTA、bcftools、或R里的SNPRelate思路相通但具体命令需要调整。我会以PLINK为主来展开因为这是最通用、最不容易踩坑的工具。核心关键词哈温平衡检验、GWAS质控、PLINK、HWE过滤、基因型数据清洗。这几个词会贯穿全文你会在每个实操环节看到它们的具体落地方式。2. 哈温平衡检验的原理与实操前准备2.1 哈温平衡到底在检验什么哈温平衡描述的是在一个理想群体中随机交配、无突变、无选择、无迁移、无限大群体等位基因频率和基因型频率之间的关系。对于一个双等位基因位点假设A等位基因频率为pa等位基因频率为q1-p那么预期基因型频率为AA p²Aa 2pqaa q²。HWE检验就是拿实际观测到的基因型计数和这个预期做比较通常用卡方检验或精确检验。卡方检验的统计量公式是χ² Σ (观测值 - 预期值)² / 预期值自由度对于双等位基因位点通常是1因为p和q加起来等于1只有一个独立参数。算出来的p值如果小于某个阈值比如1e-6就认为该位点显著偏离HWE。但这里有个关键点在病例对照研究中HWE检验通常只在对照组中做。为什么因为如果疾病和位点真的关联病例组的基因型频率本来就不应该符合HWE——这正是关联信号的表现。如果你在病例组里做HWE过滤可能把真正的关联位点给筛掉了。所以标准做法是对照组做HWE检验病例组不做或者两组分别做但只根据对照组结果过滤。2.2 实操前的数据状态确认在跑HWE之前你需要确认几件事数据已经完成了样本缺失率过滤比如--mind 0.1数据已经完成了位点缺失率过滤比如--geno 0.05数据已经完成了MAF过滤比如--maf 0.01数据格式是PLINK的bed/bim/fam或者可以转换成这个格式为什么顺序很重要因为缺失率高的位点本身就会影响HWE检验的可靠性。如果一个位点有30%的样本缺失剩下的70%样本算出来的基因型频率可能完全不能代表真实群体HWE检验结果也就没有意义。所以先做缺失率过滤再做HWE是更稳妥的流程。另外MAF过滤也要在HWE之前做。极低频的位点比如MAF0.01在HWE检验中容易出现假阳性偏离因为稀有等位基因的预期基因型计数很小卡方检验的渐近假设不成立。先过滤掉这些位点可以减少后续HWE检验的噪音。2.3 工具选型为什么用PLINKPLINK做HWE检验有几个明显优势速度快百万级SNP可以在几分钟内跑完输出格式简单直接给出每个位点的HWE p值过滤命令一行搞定不需要写额外脚本支持病例对照设计可以指定只在对照组中计算其他工具比如bcftools也有HWE检验功能bcftools stats可以输出HWE相关统计但过滤逻辑不如PLINK直接。R里的HardyWeinberg包可以做更精细的检验比如精确检验、置换检验但需要把基因型数据读进R对于大规模数据来说效率是个问题。所以我的建议是用PLINK做全基因组的HWE扫描和过滤如果后续发现某些位点需要精细检验再单独用R处理。这样兼顾了效率和灵活性。3. PLINK实操哈温平衡检验完整流程3.1 第一步生成HWE检验结果文件假设你的数据已经完成了前序质控文件前缀是gwas_qc即gwas_qc.bed、gwas_qc.bim、gwas_qc.fam。现在要计算所有位点的HWE p值plink --bfile gwas_qc \ --hardy \ --out gwas_qc_hwe这条命令会生成一个gwas_qc_hwe.hwe文件内容格式如下CHR SNP TEST A1 A2 GENO O(HET) E(HET) P 1 rs123456 ALL A G 120/340/180 0.531 0.498 0.023 1 rs234567 AFF C T 80/200/90 0.541 0.502 0.041 1 rs234567 UNAFF C T 40/140/90 0.519 0.494 0.312关键列说明TESTALL表示所有样本AFF表示病例组UNAFF表示对照组GENO基因型计数格式是“纯合A1/杂合/纯合A2”O(HET)观测杂合率E(HET)期望杂合率PHWE检验p值如果你只想要对照组的HWE结果可以加--filter-controls参数plink --bfile gwas_qc \ --hardy \ --filter-controls \ --out gwas_qc_hwe_ctrl这样输出的hwe文件里只有UNAFF行的结果后续过滤也只用这些p值。3.2 第二步根据p值阈值过滤位点拿到hwe文件后下一步是剔除显著偏离HWE的位点。常用的阈值是1e-6对于全基因组关联研究来说这个阈值比较严格能有效去除分型错误的位点同时保留可能的真实关联信号。PLINK可以直接在过滤时指定HWE阈值plink --bfile gwas_qc \ --hwe 1e-6 \ --make-bed \ --out gwas_qc_hwe_filtered注意--hwe默认是在所有样本中计算HWE p值。如果你只想在对照组中过滤需要结合--filter-controlsplink --bfile gwas_qc \ --filter-controls \ --hwe 1e-6 \ --make-bed \ --out gwas_qc_hwe_filtered_ctrl但这里有个坑--filter-controls和--hwe一起用时PLINK会先筛选对照组样本然后只在这些样本中计算HWE再根据阈值过滤位点。最终输出的数据仍然包含所有样本病例对照只是位点被过滤了。这个逻辑是对的但如果你不确认很容易误以为输出数据只剩对照组了。3.3 第三步检查过滤前后的位点数量变化过滤完成后一定要检查一下删掉了多少位点wc -l gwas_qc.bim wc -l gwas_qc_hwe_filtered.bim假设过滤前有500,000个位点过滤后剩495,000个说明删掉了5,000个位点占比1%。这个比例是否合理一般来说在质量较好的芯片数据中HWE过滤删掉的位点应该在0.5%到2%之间。如果删掉超过5%可能意味着数据质量有问题或者阈值设得太宽松/太严格。你也可以从hwe文件中直接统计显著位点的数量awk NR1 $9 1e-6 gwas_qc_hwe.hwe | wc -l这个数字应该和过滤前后位点数量的差值大致吻合如果只用ALL行的p值过滤的话。3.4 第四步可视化HWE p值分布光看数字不够直观我习惯画一个HWE p值的直方图看看分布是否合理。如果数据质量好p值应该接近均匀分布但在0附近会有一个小峰代表偏离HWE的位点。用R快速画图hwe - read.table(gwas_qc_hwe.hwe, headerTRUE, stringsAsFactorsFALSE) hwe_all - subset(hwe, TESTALL) png(hwe_pvalue_hist.png, width800, height600) hist(hwe_all$P, breaks100, mainHWE P-value Distribution, xlabHWE P-value, colsteelblue, borderwhite) abline(v1e-6, colred, lty2, lwd2) dev.off()如果直方图在0附近有一个明显的尖峰说明有一部分位点严重偏离HWE需要关注。如果整个分布很平坦只在最左边有一点点翘起那是正常的。注意画图时不要用breaks1000这种太细的分组否则图形会很难看。100个分组对于百万级位点来说已经足够了。4. 常见问题与排查技巧实录4.1 HWE过滤删掉太多位点怎么办这是最常见的问题。如果你发现HWE过滤后位点数量骤降比如从50万降到40万那肯定有问题。排查思路如下第一检查阈值是否设得太宽松。1e-6是常用阈值但如果你用的是1e-4那删掉的位点会多很多。对于全基因组数据1e-6是比较稳妥的选择。如果你做的是候选基因研究位点数量少可以用1e-4甚至1e-3但GWAS不建议。第二检查是否在病例组中做了HWE检验。如果你用了--hwe但没有加--filter-controlsPLINK会在所有样本中计算HWE。如果疾病和位点有关联病例组的偏离会拉低整体p值导致真实关联位点被误删。解决办法就是加--filter-controls只在对照组中检验。第三检查前序质控是否到位。如果缺失率过滤没做好很多位点有大量缺失HWE检验结果会很不稳定。建议先确认--geno和--mind已经执行过并且阈值合理比如--geno 0.05、--mind 0.1。第四检查群体分层。如果你的样本来自不同 ancestry 的群体混合在一起做HWE检验会导致大量位点偏离。解决办法是先做群体分层分析比如PCA然后分别在每个群体内做HWE检验或者只保留一个相对同质的群体。4.2 HWE p值全是0或全是1这种情况通常意味着数据有问题。如果所有位点的HWE p值都是0可能是基因型编码出了问题比如所有样本都是纯合子没有杂合子。你可以检查一下bim文件里的等位基因编码或者用--freq看看等位基因频率是否正常。如果所有p值都是1那可能是样本量太小或者位点全是单态MAF0。先做MAF过滤再重新跑HWE。4.3 如何选择HWE阈值1e-6还是1e-4这个问题没有绝对答案取决于你的研究设计和数据质量。我的一般建议是研究类型推荐HWE阈值理由全基因组关联研究GWAS1e-6严格过滤避免假阳性候选基因研究1e-4位点少可以适当放宽拷贝数变异相关研究1e-4 ~ 1e-3CNV区域HWE容易偏离放宽阈值群体遗传学研究1e-3 ~ 1e-2关注的是群体结构不是关联信号另外如果你做的是meta分析各个队列的HWE阈值最好统一否则合并结果会有偏倚。4.4 HWE过滤后是否需要重新做缺失率过滤不需要。HWE过滤只是删除了部分位点不会引入新的缺失。但如果你在HWE过滤后做了其他操作比如合并数据集、转换格式那可能需要重新检查缺失率。不过有一个例外如果你在HWE过滤时用了--filter-controls输出数据仍然包含所有样本但位点是根据对照组筛选的。这种情况下病例组的缺失率可能和对照组不同但一般不需要重新过滤因为位点层面的缺失率在之前已经控制过了。4.5 常见问题速查表问题现象可能原因解决方法HWE过滤删掉超过5%位点阈值太宽松或病例组参与检验改用1e-6加--filter-controlsHWE p值全为0基因型编码错误或全纯合检查bim文件用--freq验证HWE p值全为1样本量太小或位点单态先做MAF过滤增加样本量过滤后位点数不变阈值太严格或命令写错检查--hwe参数是否生效病例组和对照组HWE结果差异大真实关联或群体分层检查PCA结果分群体检验5. 进阶技巧与个人经验分享5.1 分群体做HWE检验的策略如果你的样本包含多个群体比如不同地区、不同种族混合做HWE检验会引入很大偏差。我通常的做法是先用PCA或ADMIXTURE确定每个样本的群体归属按群体拆分数据--keep指定样本列表在每个群体内分别做HWE检验只保留在所有群体中都符合HWE的位点这样虽然会删掉更多位点但能有效避免群体分层导致的假阳性。代价是可能丢失一些群体特异性位点但对于关联研究来说稳妥比激进更重要。5.2 用精确检验替代卡方检验PLINK默认用的是卡方检验对于低频位点或小样本卡方检验的渐近假设可能不成立。这时候可以用精确检验。PLINK本身不直接支持精确检验但你可以把数据导出到R用HardyWeinberg包做library(HardyWeinberg) # 假设你有一个基因型计数矩阵 # 每行是一个位点三列分别是AA、Aa、aa的计数 hw - HWExact(genotype_counts)精确检验更准确但计算量大不适合全基因组扫描。我的建议是先用PLINK卡方检验做全基因组过滤对过滤后保留的位点如果p值在1e-6到1e-4之间可以用精确检验复核。这样既保证了效率又提高了准确性。5.3 HWE过滤与后续分析的衔接HWE过滤完成后你的数据就可以进入下一步了。通常接下来的步骤是连锁不平衡剪枝LD pruning群体分层校正PCA关联分析logistic回归或线性回归这里有一个经验HWE过滤后的位点列表最好保存下来方便后续追溯。你可以把过滤后的bim文件单独备份或者记录下被删掉的位点ID。如果后续分析发现某些重要位点被HWE过滤掉了你可以快速定位并决定是否重新纳入。5.4 一个容易忽略的细节X染色体PLINK的--hwe默认只对常染色体位点做检验。如果你需要检验X染色体位点需要加--chr X或者--hwe配合--nonfounders等参数。但X染色体的HWE检验逻辑和常染色体不同因为男性只有一条X基因型频率的计算方式不一样。对于GWAS来说通常只对常染色体做HWE过滤X染色体单独处理。如果你确实需要处理X染色体建议用--chr X --hwe 1e-6并且只对女性样本做检验因为女性有两条X符合HWE假设。男性样本的X染色体基因型只有纯合子不适用HWE检验。5.5 实操心得不要盲目追求“干净”的数据我见过一些人做GWAS恨不得把所有偏离HWE的位点都删掉结果最后剩下的位点全是MAF很高、HWE完美符合的“无聊”位点。这样虽然减少了假阳性但也可能丢掉了真正有生物学意义的信号。我的建议是HWE过滤是质控手段不是目的。1e-6的阈值已经足够严格不需要再收紧。如果你发现某个位点HWE p值在1e-5左右但文献报道它和你的研究表型有关那可以保留它在后续分析中重点关注。GWAS的本质是发现关联不是制造完美数据。另外HWE过滤的结果应该和缺失率过滤、MAF过滤的结果一起记录在质控日志里。这样你写论文的时候可以直接引用质控流程审稿人也会觉得你的分析很规范。5.6 批量处理多个数据集的脚本模板如果你有多个数据集需要做HWE过滤可以写一个简单的bash脚本批量处理#!/bin/bash for prefix in dataset1 dataset2 dataset3; do echo Processing $prefix... plink --bfile $prefix \ --filter-controls \ --hwe 1e-6 \ --make-bed \ --out ${prefix}_hwe_filtered echo Before: $(wc -l ${prefix}.bim) variants echo After: $(wc -l ${prefix}_hwe_filtered.bim) variants done这个脚本会依次处理每个数据集输出过滤前后的位点数量方便你快速检查。注意--filter-controls要求你的fam文件里病例对照编码正确通常是1对照2病例或者反过来取决于你的编码方式。6. 总结与后续扩展方向哈温平衡检验在GWAS质控流程中是一个承上启下的环节。它承接了缺失率过滤和MAF过滤又为后续的关联分析提供了更干净的数据。PLINK的--hwe命令让这个步骤变得非常简单但背后的逻辑和参数选择需要仔细斟酌。我个人在实际操作中的体会是HWE过滤的关键不在于命令本身而在于理解你的数据和研究设计。病例对照研究只在对照组做检验多群体样本要分群体做低频位点要谨慎处理。这些细节决定了你的质控是“有效过滤”还是“误伤好人”。这个内容后续还可以这样扩展如果你做的是全基因组测序数据HWE检验的逻辑和芯片数据略有不同因为测序数据的基因型 calling 误差模式不一样。另外如果你用的是imputation后的数据HWE检验的阈值可能需要调整因为imputation会引入额外的不确定性。这些话题都值得单独展开但核心思路是一致的用HWE检验筛掉明显有问题的位点但不要过度过滤。