GWAS质控中哈温平衡检验的实操指南与避坑心得

1. 为什么哈温平衡检验是GWAS质控里最容易被忽视的一环做全基因组关联分析的人十有八九把注意力放在关联模型的显著性阈值、群体分层校正、曼哈顿图好不好看这些事情上而哈温平衡检验往往被当成一个顺手跑一下的步骤。我刚开始做GWAS那会儿也是这个心态觉得这个检验不就是算个卡方看看p值嘛能有什么门道。后来踩了几次坑才发现这个环节处理不好轻则损失一批真实的关联信号重则让整个分析结果出现系统性偏差审稿人一眼就能看出来。哈温平衡检验英文叫Hardy-Weinberg Equilibrium test简称HWE检验本质上是在检验一个群体里的等位基因频率和基因型频率是否符合理论预期。在一个理想的大规模随机交配群体中如果没有选择、突变、迁移和遗传漂变的影响基因型频率应该满足p²、2pq、q²这个经典分布。GWAS里做这个检验核心目的不是去研究群体遗传学本身而是把它当作一个质控过滤器——把那些明显偏离平衡的位点揪出来因为它们大概率是基因分型错误、拷贝数变异干扰或者样本污染导致的假信号。这里有个关键的分寸问题HWE偏离的位点到底该不该删删多少阈值设多少这些问题没有标准答案但有一套成熟的判断逻辑。我见过有人一刀切把所有p小于0.001的位点全删了结果把一个真实的疾病易感位点也删掉了因为那个位点本身就和疾病相关在病例组里天然就不符合哈温平衡。这就是典型的把洗澡水和孩子一起倒掉。所以这篇文章我想把HWE检验这件事从头到尾讲透包括它在整个GWAS流程里的位置、底层原理、实操代码、阈值选择的逻辑、病例对照研究里的特殊处理方式以及我自己踩过的那些坑。适合已经跑过基础GWAS流程、正在做质控优化的朋友也适合刚接触遗传统计、想搞清楚每一步到底在干什么的新手。2. HWE检验在GWAS质控体系中的定位与核心原理2.1 它在整个流程里到底排在第几步一个完整的GWAS质控流程大致会经过样本层面和位点层面两轮过滤。样本层面包括call rate过滤、杂合度检查、性别一致性验证、亲缘关系剔除、群体分层评估等位点层面则包括call rate过滤、次等位基因频率过滤、HWE检验、分型质量打分过滤等。HWE检验属于位点层面的质控通常排在MAF过滤之后、关联分析之前。为什么是这个位置因为MAF太低的位点本身统计功效就不足先删掉可以减轻后续计算负担而HWE检验放在关联分析前是为了避免那些分型质量差的位点进入统计模型污染结果。但这里有个细节很多人没注意HWE检验应该分群体做而不是把所有样本混在一起做。如果你手头有东亚、欧洲、非洲多个群体的样本混在一起算HWE由于群体分层效应几乎每个位点都会显著偏离平衡。正确做法是先按遗传背景分组在每个群体内部单独检验或者至少在关联分析时把群体主成分作为协变量纳入模型后再看HWE。2.2 卡方检验的数学逻辑到底是什么HWE检验最常用的方法是Pearson卡方检验。假设某个SNP有两个等位基因A和a在N个样本中我们观察到三种基因型AA、Aa、aa的实际数量分别是n_AA、n_Aa、n_aa。等位基因频率可以这样算p(A) (2×n_AA n_Aa) / (2N)q(a) (2×n_aa n_Aa) / (2N)理论上如果群体处于哈温平衡三种基因型的期望数量应该是E(AA) N × p²E(Aa) N × 2pqE(aa) N × q²然后构造卡方统计量χ² Σ [(观察值 - 期望值)² / 期望值]自由度是1因为三个基因型频率之和为1减去一个估计的参数p所以df 3-1-1 1。算出卡方值后查表或者用精确检验得到p值。但这里有个坑当样本量较小或者次等位基因频率很低时卡方近似不准确这时候应该用精确检验exact test。PLINK默认用的是精确检验这一点比很多自己手写卡方的人要靠谱。精确检验基于的是给定等位基因总数条件下观察到的基因型计数的条件概率分布计算量更大但结果更准确。还有一个容易被忽略的点对于病例对照研究HWE检验应该在对照组中做而不是在全体样本或病例组中做。原因很简单——如果这个位点真的和疾病相关那么病例组的基因型频率本来就会偏离平衡这是真实的生物学信号不是错误。在对照组中检验对照组是未患病的人群理论上应该符合哈温平衡如果偏离了那更可能是分型错误。2.3 偏离HWE到底意味着什么一个位点偏离HWE可能的原因有很多种我按常见程度排个序偏离原因典型表现处理方式基因分型错误杂合子缺失或过量删除位点或重新分型拷贝数变异杂合子过量簇图异常删除或单独处理样本污染多个位点同时偏离检查样本层面质控群体分层大量位点轻微偏离校正群体结构真实关联仅病例组偏离保留位点选择效应特定群体中偏离结合生物学背景判断近亲交配杂合子系统性缺失检查样本亲缘关系杂合子缺失heterozygote deficiency是最常见的偏离模式通常指向分型错误或者等位基因丢失。杂合子过量heterozygote excess则更可能和拷贝数变异、样本污染或者测序错误有关。这两种模式的排查思路完全不同后面我会详细讲。3. 实操前的准备工作与工具选型3.1 数据格式要求与前置检查做HWE检验之前你的数据至少应该是经过基础质控的。我一般要求满足这几个条件样本call rate 95%有些严格的项目要求98%位点call rate 95%已经完成性别一致性检查已经剔除亲缘关系过近的样本如果是多群体数据已经完成群体分层评估数据格式方面PLINK的bed/bim/fam三件套是最通用的。如果你用的是VCF可以先用PLINK转换成bed格式。转换命令plink --vcf input.vcf.gz --make-bed --out gwas_data --double-id --allow-extra-chr这里--double-id是为了处理VCF里样本ID和群体ID可能重复的问题--allow-extra-chr是允许非标准染色体名称具体要不要加看你的数据情况。转换完之后先看一眼bim文件和fam文件的基本情况wc -l gwas_data.bim wc -l gwas_data.fam head gwas_data.bim head gwas_data.fambim文件有6列染色体、SNP ID、遗传距离、物理位置、次等位基因、主等位基因。fam文件也有6列家系ID、个体ID、父本ID、母本ID、性别、表型。表型列里1通常代表对照2代表病例-9代表缺失。这些基础信息确认无误后再往下走。3.2 工具选型PLINK还是其他做HWE检验PLINK是最主流的选择没有之一。原因有几个一是速度快C语言写的底层几百万个位点几分钟就能跑完二是功能全--hardy参数直接输出HWE检验结果三是社区成熟遇到问题容易找到解决方案。PLINK有两个主要版本1.9和2.0。1.9的--hardy功能更成熟输出格式也更友好我一般推荐用1.9做HWE检验。2.0虽然性能更好但某些质控功能的参数设计不太一样容易搞混。除了PLINKR语言里的HardyWeinberg包也值得了解一下它提供了更丰富的检验方法包括精确检验、置换检验等适合做小规模数据的深入分析。但对于全基因组级别的数据还是PLINK更实用。还有一个工具叫bcftools它的hardy插件也能做HWE检验适合直接在VCF文件上操作省去格式转换的步骤。命令大概是bcftools hardy input.vcf.gz -o hwe_results.txt不过bcftools的HWE功能相对简单没有PLINK那么灵活的过滤选项。3.3 分群体检验的实操方案如果你的数据包含多个群体分群体检验是必须的。假设你的fam文件里有一列或者你有一个单独的群体注释文件可以这样操作# 假设群体信息在pheno.txt里第一列是FID第二列是IID第三列是群体标签 # 先按群体拆分数据 for pop in EUR EAS AFR; do awk -v p$pop $3p {print $1, $2} pheno.txt ${pop}_samples.txt plink --bfile gwas_data --keep ${pop}_samples.txt --make-bed --out gwas_${pop} done # 对每个群体单独做HWE检验 for pop in EUR EAS AFR; do plink --bfile gwas_${pop} --hardy --out hwe_${pop} done这样每个群体得到一个.hwe文件里面包含了每个位点的HWE检验结果。后续可以分别过滤也可以取交集或并集具体策略取决于你的研究目的。注意分群体检验时每个群体的样本量不能太小。如果某个群体只有几十个样本HWE检验的统计功效很低结果参考价值有限。一般建议每个群体至少100个样本以上再做分群体检验。4. HWE检验的完整实操流程与参数详解4.1 PLINK --hardy 命令的完整用法最基本的HWE检验命令plink --bfile gwas_data --hardy --out hwe_results这条命令会生成一个hwe_results.hwe文件包含以下列CHR染色体SNPSNP IDTEST检验类型包括ALL、AFF、UNAFF三种A1次等位基因A2主等位基因GENO基因型计数格式是A1A1/A1A2/A2A2O(HET)观察到的杂合子频率E(HET)期望的杂合子频率PHWE检验的p值这里TEST列很关键。ALL表示用所有样本做的检验AFF表示只用病例组UNAFF表示只用对照组。对于病例对照研究你应该重点关注UNAFF这一行。如果你只想要对照组的HWE结果可以用plink --bfile gwas_data --hardy --filter-controls --out hwe_controls--filter-controls会只保留表型为1的样本对照。但要注意这个参数在不同版本的PLINK里行为可能略有差异建议跑完之后检查一下样本数量对不对。4.2 阈值选择0.001还是1e-6这是被问得最多的问题。我的经验是没有万能阈值要看研究类型和数据质量。对于传统的候选基因研究或者样本量较小的GWAS几千人以下常用的阈值是p 0.001。对于大规模GWAS几万人以上由于统计功效很高即使是轻微的偏离也会显著这时候阈值要更严格常用p 1e-6甚至p 1e-10。但这里有个反直觉的点在大规模GWAS中HWE p值显著并不一定意味着位点有问题。样本量大了之后微小的分型偏差或者群体分层残留都会被放大。所以在大规模研究中我一般会结合其他质控指标一起判断比如分型质量打分info score、杂合子频率是否异常等。还有一个策略是分层过滤对常见位点MAF 5%用较严格的阈值如1e-6对罕见位点MAF 1%用较宽松的阈值如1e-4。因为罕见位点的HWE检验本身就不太可靠卡方近似误差大精确检验的计算也容易不稳定。我自己的习惯是研究类型样本量推荐HWE阈值备注候选基因 50000.001宽松避免丢失真实信号中等GWAS5000-500001e-4平衡假阳性和假阴性大规模GWAS 500001e-6严格控制假阳性罕见变异任意1e-4检验本身不可靠放宽病例对照任意仅看对照组病例组偏离可能是真实信号4.3 过滤命令与结果验证确定阈值后用--hwe参数过滤plink --bfile gwas_data --hwe 1e-6 --make-bed --out gwas_hwe_filtered如果是病例对照研究只想在对照组中过滤plink --bfile gwas_data --hwe 1e-6 --hwe-filter-controls --make-bed --out gwas_hwe_filtered跑完之后一定要验证一下过滤效果# 过滤前位点数 wc -l gwas_data.bim # 过滤后位点数 wc -l gwas_hwe_filtered.bim # 看删了多少我一般会记录过滤前后的位点数量如果删得太多比如超过5%就要回头检查是不是阈值设得太松或者数据质量有问题。正常情况下HWE过滤应该只删掉很小一部分位点通常小于1%。还有一个验证步骤是看被删位点的特征。把被删的位点单独拿出来看看它们的MAF分布、call rate分布、在染色体上的位置分布。如果被删的位点集中在某些染色体区域可能提示这些区域有拷贝数变异或者分型问题。# 提取被删的位点 comm -23 (cut -f2 gwas_data.bim | sort) (cut -f2 gwas_hwe_filtered.bim | sort) removed_snps.txt # 统计被删位点的数量 wc -l removed_snps.txt5. 病例对照研究中的特殊处理与常见误区5.1 为什么不能在全样本中做HWE这是我最想强调的一点。很多教程和自动化流程默认在全样本中做HWE检验这在病例对照研究中是错误的。原因用一句话说清楚如果这个位点真的和疾病相关那么病例组的基因型频率分布和对照组不同合并在一起算HWE自然会偏离平衡。但这个偏离是真实的生物学信号不是错误。如果你在全样本中过滤就会把这个真实的关联位点删掉。我见过一个真实的案例当然细节做了模糊处理某研究团队做精神分裂症的GWAS在全样本中做HWE过滤阈值设了0.001结果把一个已知的强关联位点删掉了。后来审稿人要求补充分析他们才发现这个问题重新在对照组中做HWE后那个位点不仅保留了下来还是全基因组最显著的信号之一。所以正确的做法是HWE检验只在对照组中做过滤也只在对照组中过滤。病例组的HWE结果可以参考但不用来过滤。5.2 对照组HWE偏离的排查思路如果对照组中也有大量位点偏离HWE那就要认真排查了。我一般按这个顺序查第一步检查样本污染。看样本的杂合度分布如果有样本的杂合度明显偏离群体均值比如超过3个标准差可能是样本污染或混样。可以用PLINK的--het参数算杂合度plink --bfile gwas_data --het --out het_check然后看het_check.het文件里的F值近交系数F值显著偏离0的样本要警惕。第二步检查群体分层。如果对照组本身包含多个遗传背景差异较大的亚群体HWE偏离是必然的。可以用PCA或者ADMIXTURE评估群体结构然后在回归模型中加入主成分作为协变量。第三步检查分型质量。看偏离位点的分型簇图cluster plot如果三个基因型的簇分不开或者杂合子簇偏移那就是分型算法的问题。不同平台如Affymetrix和Illumina的分型特征不同需要针对性排查。第四步检查拷贝数变异。杂合子过量的位点往往和CNV有关。可以用CNV检测工具如PennCNV扫描一遍看看偏离位点是否落在CNV区域。5.3 常见误区速查表误区后果正确做法全样本做HWE删除真实关联位点仅在对照组中做阈值一刀切要么删太多要么删太少按MAF和研究类型分层设阈值忽略群体分层大量位点假性偏离分群体检验或校正主成分只看p值不看方向漏掉杂合子缺失/过量信息同时看O(HET)和E(HET)过滤后不验证可能删了不该删的检查被删位点的特征分布罕见位点用严格阈值误删大量真实罕见变异罕见位点放宽阈值6. 结果解读、可视化与下游衔接6.1 怎么读懂.hwe文件拿到.hwe文件后不要只看P列。我一般会做几个衍生指标杂合子偏离度 (O(HET) - E(HET)) / E(HET)正值表示杂合子过量负值表示杂合子缺失等位基因频率 从GENO列反推p (2×n_A1A1 n_A1A2) / (2N)检验类型 确认是ALL、AFF还是UNAFF用R或者Python快速算一下hwe - read.table(hwe_results.hwe, headerTRUE) hwe$het_dev - (hwe$O.HET. - hwe$E.HET.) / hwe$E.HET. # 看杂合子偏离的分布 hist(hwe$het_dev, breaks100, mainHeterozygote Deviation) # 看显著位点的偏离方向 sig - hwe[hwe$P 1e-6 hwe$TEST UNAFF, ] table(sign(sig$het_dev))如果显著位点里杂合子缺失占多数那更可能是分型错误如果杂合子过量占多数那要怀疑CNV或污染。6.2 可视化QQ图和曼哈顿图上的HWE虽然HWE检验本身不直接出QQ图但你可以把HWE的p值做成QQ图看看整体分布是否偏离均匀分布。如果大量位点的HWE p值集中在低值区说明数据存在系统性问题。# HWE p值的QQ图 observed - sort(hwe$P[hwe$TEST UNAFF]) expected - ppoints(length(observed)) plot(-log10(expected), -log10(observed), xlabExpected -log10(P), ylabObserved -log10(P), mainHWE QQ Plot, pch20, cex0.5) abline(0, 1, colred)如果点明显偏离对角线尤其是在低端翘起说明有大量位点偏离HWE需要排查原因。另一个有用的可视化是杂合子偏离度沿染色体的分布图。如果某些染色体区域集中出现杂合子缺失或过量可能提示这些区域有结构变异或分型问题。6.3 与下游分析的衔接HWE过滤完成后数据就可以进入关联分析了。但这里有个衔接细节HWE过滤后的位点集要和后续的关联分析、meta分析保持一致。如果你做了多个批次的基因分型每个批次单独做HWE过滤最后合并时要注意位点集的一致性。另外如果你打算做条件分析或者精细定位HWE过滤的阈值可以适当放宽因为这些分析对位点的完整性要求更高删掉太多位点会影响精细定位的分辨率。还有一个实践中的建议保留一份未经过HWE过滤的数据备份。因为HWE过滤是有争议的步骤有些审稿人可能会要求你提供不过滤的结果作为敏感性分析。保留原始数据可以让你在需要时快速重新分析。7. 我踩过的坑与实操心得7.1 那些年我删错的位点刚开始做GWAS的时候我按照某个教程的推荐把HWE阈值设成了0.001在全样本中过滤。结果跑完关联分析发现一个之前文献里报道过的强关联位点不见了。查了半天才发现那个位点在病例组中显著偏离HWE因为它是真实关联位点在全样本检验中p值很小被我的过滤步骤删掉了。从那以后我养成了一个习惯每次HWE过滤后都把被删的位点和已知的关联位点数据库如GWAS Catalog做个比对。如果被删的位点里有已知的关联信号就要格外警惕重新评估过滤策略。7.2 群体分层带来的假阳性偏离有一次处理一个多群体数据我偷懒没有分群体做HWE直接在全样本中跑。结果发现将近10%的位点HWE p值小于0.001这个比例远远超出预期。后来分群体重新做每个群体内部的偏离比例都降到了正常水平小于1%。这就是典型的群体分层导致的假性偏离。这个经历让我明白HWE检验的结果必须结合群体结构来解读。如果你的数据有明显的群体分层分群体检验不是可选项而是必选项。7.3 罕见位点的HWE检验不可靠罕见变异MAF 1%的HWE检验有个根本性的问题当次等位基因频率很低时三种基因型的期望计数可能小于5卡方近似完全失效。虽然PLINK默认用精确检验但精确检验在极低频率下也不稳定。我的做法是对MAF 1%的位点不做HWE过滤或者用非常宽松的阈值如1e-4。因为这些位点本身在关联分析中的统计功效就很低HWE检验的假阳性率又高严格过滤只会误伤。7.4 过滤后的验证不能省我见过太多人跑完过滤命令就直接进入下一步不做任何验证。结果要么是过滤没生效参数写错了要么是过滤过头了阈值设错了等到关联分析结果出来才发现问题那时候排查成本就高了。我的标准操作是每次过滤后检查三个数——过滤前位点数、过滤后位点数、被删位点数。如果被删比例超过2%就要停下来看看是不是哪里出了问题。另外随机抽几个被删的位点手动算一下它们的HWE p值验证PLINK的输出是否正确。7.5 一个实用的小技巧如果你不确定阈值设多少合适可以做一个敏感性分析用几个不同的阈值如0.001、1e-4、1e-6分别过滤然后比较关联分析的结果。如果不同阈值下结果基本一致说明你的数据质量很好阈值选择不敏感如果结果差异很大说明数据质量有问题需要先解决质控问题再谈阈值。这个敏感性分析虽然多花一点时间但能帮你避免很多后续的麻烦。尤其是在投稿前审稿人如果质疑你的HWE过滤策略你可以直接拿出敏感性分析的结果来回应。最后再分享一个我在实际操作中的体会HWE检验不是一个孤立的步骤它和样本call rate、杂合度、群体分层、分型质量这些质控指标是相互关联的。单独看HWE的p值容易误判把多个指标放在一起看才能做出准确的判断。比如一个位点HWE p值很小但它的call rate很高、分型簇图很清晰、MAF也正常那它可能只是一个真实的关联信号反过来如果一个位点HWE p值只是轻微偏离但call rate偏低、杂合度异常那它更可能是分型错误。这种综合判断的能力需要在实践中慢慢积累。