企业官网建设流程全解析

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 × 2pq
• E(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.fam

bim文件有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虽然性能更好，但某些质控功能的参数设计不太一样，容易搞混。

除了PLINK，R语言里的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 '$3==p {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：染色体
• SNP：SNP ID
• TEST：检验类型，包括ALL、AFF、UNAFF三种
• A1：次等位基因
• A2：主等位基因
• GENO：基因型计数，格式是A1A1/A1A2/A2A2
• O(HET)：观察到的杂合子频率
• E(HET)：期望的杂合子频率
• P：HWE检验的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检验本身就不太可靠，卡方近似误差大，精确检验的计算也容易不稳定。

我自己的习惯是：

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.txt

5. 病例对照研究中的特殊处理与常见误区

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 常见误区速查表

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", header=TRUE) hwe$het_dev <- (hwe$O.HET. - hwe$E.HET.) / hwe$E.HET. # 看杂合子偏离的分布 hist(hwe$het_dev, breaks=100, main="Heterozygote 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), xlab="Expected -log10(P)", ylab="Observed -log10(P)", main="HWE QQ Plot", pch=20, cex=0.5) abline(0, 1, col="red")

如果点明显偏离对角线，尤其是在低端翘起，说明有大量位点偏离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偏低、杂合度异常，那它更可能是分型错误。这种综合判断的能力，需要在实践中慢慢积累。