企业官网建设流程全解析

在上一篇里，我们已经把PLINK的基础操作和数据格式过了一遍，从PED/MAP到BED/BIM/FAM，以及怎么用--recode、--make-bed做格式转换。不少朋友留言说按照那套流程已经能把数据跑起来了，但真正开始做关联分析的时候，又发现结果总是不对劲——要么曼哈顿图上一堆假阳性，要么QQ图的lambda值高得离谱。这其实不是关联分析的命令写错了，而是漏掉了最关键的一步：质控（QC）。

GWAS圈子里有句老话：“垃圾进，垃圾出”。如果你的样本和位点没有经过严格的筛选，那后面无论你用GEMMA、SAIGE还是BOLT-LMM，算出来的p值都不可信。这一篇，我们就专门讲PLINK做GWAS的质控环节。我会直接把我在实际项目中用的参数、脚本和判断逻辑摊开来说，尤其是那些文档里不会写、但踩过坑才知道的细节。

1. 为什么GWAS必须做质控：三个真实场景

1.1 基因分型误差如何变成假阳性

先聊个反直觉的事情。很多人觉得，现在芯片分型技术这么成熟，出来的数据还需要什么质控？但实际上，芯片数据里的错误远比你想的多。DNA样本在采集、运输、DNA提取、扩增、杂交、扫描的每一个环节都可能引入误差，而这些误差在统计上会直接表现为假阳性。

我给你举一个我实际遇到过的例子。某个项目里有2000个病例和2000个对照，做完关联分析以后，发现rs1234567这个位点的p值到了1e-20，效应量OR=2.3。这个结果看着非常漂亮，但仔细一查发现，病例组的DNA样本提取时间比对照组早了一年，两组样本在芯片上跑的批次也不一样。重新做质控，把分型成功率低于95%的样本剔除，再把缺失率偏高的位点过滤掉以后，这个位点的p值直接变成了0.3。

为什么会这样？因为基因芯片的荧光信号强度会受到样本DNA质量的影响。DNA质量差的样本，在某个基因型簇上的信号就会模糊，分型算法给出的结果就不稳定。如果这种不稳定性恰好和病例/对照的分组相关，就会产生系统性的分型误差，在统计上呈现出显著关联。

1.2 人群分层会让结果完全失真

第二个常见问题是人群分层。GWAS的前提假设是：除了我们关注的位点外，病例组和对照组的遗传背景应该是一致的。但现实里这个假设几乎不成立。

举个更生活化的例子。假设你在研究某个跟东亚人群相关的疾病，收了500个病例和500个对照。病例来自北方某医院，对照来自南方某社区。这时候你会发现，哪怕你随便拿10000个与疾病无关的SNP去做检验，也会有一大堆“显著”位点。原因很简单：北方和南方人群的等位基因频率本身就有差异，这些差异和疾病无关，但在统计检验里会被当成关联信号。

人群分层在统计上会拉高基因膨胀因子lambda。lambda接近1是理想状态，超过1.1就要警惕了。如果超过1.15，说明数据里有严重的人群结构问题，这个问题的根源不是关联分析本身，而是前面的QC没做到位。

1.3 样本关系混杂会高估显著性

第三个问题出在样本关系上。如果你项目里有 unknowingly 收录了亲属样本——比如一对父子、一对同卵双胞胎——那这些样本的基因型高度相似，在做关联分析时会人为地增加统计功效，导致假阳性率升高。

我见过一个极端的案例。某GWAS项目做完QC后，样本量从1500降到了1100，就是因为有相当一部分样本的PL_HAT（亲缘关系估计值）大于0.2。如果不做这个剔除，最终的关联信号会被严重夸大，而且复现性极差——你换一个独立队列验证，结果可能就消失了。

所以，QC这步绝对不能省。它看似是在“丢失”数据，实际上是在帮你提高后续分析的信噪比。下面我们就进入正题，用PLINK把整个过程跑一遍。

2. 个体层面的QC：把劣质样本挡在门外

个体层面的质控，目标是剔除那些DNA质量差、性别信息有误、样本之间存在亲缘关系的个体。这一步做扎实了，后面SNP层面的过滤才有意义。

2.1 检查缺失率与性别一致性

先说最基础的：样本缺失率。如果一个样本在大部分位点上都没有成功分型，说明这个样本的DNA质量不行，或者样本本身的浓度有问题。PLINK的命令很简单：

plink --bfile gwas_raw --missing --out qc_snp

运行后会生成plink.imiss文件，其中F_MISS列表示每个样本的缺失率。一般我们会以0.05为阈值，也就是缺失率超过5%的样本直接剔除。

但这里要插一句经验之谈。对于老样本、FFPE样本或者抽提很久的DNA，缺失率阈值可以放宽到0.1。因为这类样本本身质量就差，如果你卡5%，可能一半样本都没了。这时候更合理的做法是，把阈值设在0.1，同时结合其他指标综合判断。

性别检查也是必做项。基因芯片上会包含一些X染色体和Y染色体上的性别标记位点，PLINK可以根据这些位点的杂合度来推断样本的遗传性别。如果这个推断结果和样本记录的临床性别不一致，那几乎可以断定样本被搞混了或者污染了，直接剔除。

plink --bfile gwas_raw --check-sex --out qc_sex

生成的plink.sexcheck文件里，STATUS列显示OK或PROBLEM。碰到PROBLEM的样本，我的习惯是直接看PEDSEX和SNPSEX的差异。如果差异巨大，那就剔除；如果只是因为个别位点噪声导致的边界情况，可以结合其他指标再定。

2.2 用杂合率偏离程度筛查DNA污染

接下来是一个经常被忽略但特别有效的指标：常染色体杂合率。如果一个样本的DNA被另一种DNA污染了，杂合率会明显偏离群体平均水平。

这个概念我可以用一个很有意思的类比来说清楚。你把两种不同颜色的豆子混在一起，虽然颜色比例没变，但你随机抓一把时，抓到“一红一绿”这种异色组合的概率会变高。DNA污染就是这样，会人为地增加“看起来是杂合”的机会。

计算杂合率的命令：

plink --bfile gwas_raw --het --out qc_het

这里的F列是近交系数估算值，O(HOM)和E(HOM)分别代表观察到的纯合子计数和期望纯合子计数。杂合率异常高的样本，往往意味着DNA污染；杂合率异常低的样本，往往说明样本有近亲关系或者DNA质量极差。

那阈值怎么定呢？我的经验是，采用均值加减3倍标准差的策略。先用R快速统计一下F列的分布，然后把超出3倍标准差的样本剔除。这样做的好处是阈值是根据数据本身分布自动调整的，比硬编码一个固定值更稳健。

2.3 亲缘关系筛选：剔除duplicate和close relatives

第三步是检查样本间的亲缘关系。PLINK里最常用的是--genome命令，它会估算两两样本之间的PI_HAT值。PI_HAT大于0.1875相当于三阶亲属关系，大于0.5就是同卵双胞胎或者重复样本。

plink --bfile gwas_raw --genome --min 0.2 --out qc_related

这里我用的是--min 0.2，直接输出PI_HAT大于0.2的样本对。拿到结果后，需要手动决定保留哪个样本。通常的做法是保留缺失率更低的那个。

这里有个细节：做亲缘关系筛选前，最好先用一个独立的LD修剪后的SNP集，而不是全基因组所有位点。因为高LD区域会让亲缘估算产生偏差。你要先在QC的基础上用--indep-pairwise做一次LD pruning，然后用修剪后的SNP集跑--genome。

我在实际项目中一般会用：

plink --bfile gwas_qc_indiv --indep-pairwise 50 5 0.2 --out qc_prune plink --bfile gwas_qc_indiv --extract qc_prune.prune.in --genome --min 0.2 --out qc_related

3. SNP层面的QC：把不可靠的位点过滤掉

个体层面清干净以后，我们把目光转到每个SNP位点上。

3.1 SNP缺失率与差异缺失率

SNP层面的缺失率，跟个体层面的逻辑类似。一个位点在很多样本里都没分出来，那这个位点的芯片探针设计可能有问题，或者这个位点附近的序列比对有多义性。常规阈值是0.05，但也要结合具体芯片和样本情况来定。

差异缺失率这个指标更值得关注。它的意思是在病例和对照两组之间，某个位点的缺失率是否存在显著差异。如果一个位点在病例组缺失率5%，在对照组缺失率0.1%，即便单独看都不超标，但两组差异本身就预示着这个位点的分型结果可能受某种与表型相关的因素干扰，这在统计上会直接导致假阳性。

PLINK支持用--test-missing直接做这个检验：

plink --bfile gwas_qc_indiv --test-missing --out qc_diffmiss

输出的P列如果小于1e-4，就得把这个位点滤掉。实操里我会把差异缺失率的阈值卡得更严，直接用--missing输出两份统计，再用R计算Fisher精确检验，把P<0.0001的位点剔除。原因还是那句话：这种位点即便p值好看，进了下游分析就是隐患。

3.2 MAF过滤的逻辑与阈值选择

MAF（最小等位基因频率）过滤，是GWAS里最容易让人困惑的一个环节。很多人不理解：为什么低频位点不能被直接纳入分析？

核心原因有两个。第一，低频位点的基因型计数很少，统计检验的病态行为会很严重——比如某个位点只有两个样本是杂合，这俩样本恰好都在病例组，p值就可能极其显著。这在统计上叫小样本偏差。第二，低频位点在芯片上的可靠性本身就差，重复性低，验证成本高。

MAF阈值怎么选，要看你项目的研究目的。如果你的样本量在几千这个量级，建议卡0.05，也就是5%的MAF。如果你的样本量上万，可以考虑降到0.01。如果是为了做罕见变异分析，那就不应该用PLINK的常规设置来做，而是应该用专门的工具和专门的质控流程。

plink --bfile gwas_qc_indiv --maf 0.05 --make-bed --out gwas_qc_snp1

这里我想多提一句：很多人以为MAF过滤只是简单的“把低频的去掉”，其实不是。MAF过滤还承担着一个功能，就是减少后续关联分析的检验负担。GWAS动辄几十万、几百万个位点，多保留一些低频位点，多重检验校正的压力就更大。虽然现在计算能力不是瓶颈，但如果你用的是Bonferroni校正，那0.05/500000和0.05/800000的差别还是不可忽视的。

3.3 HWE检验：哪些位点应该被过滤

HWE过滤的逻辑有一点容易搞反，我特意把它单独拿出来说。

在做GWAS时，我们的做法是：在对照组里做HWE检验，然后把偏离HWE的位点剔除。阈值一般是1e-6甚至更严。

为什么是在对照组而不是全样本？因为在病例组里，如果一个位点真的与疾病相关，那病例组的等位基因频率本来就偏离群体期望，HWE检验自然也会偏离。如果你在全样本或者只在病例组做HWE过滤，很可能会把真正与疾病相关的位点给误删了。

在对照组里做HWE，就可以过滤掉那些因为分型错误导致的偏离HWE的位点，同时保证真实关联信号不受影响。

PLINK的命令：

plink --bfile gwas_qc_snp1 --filter-controls --hwe 1e-6 --make-bed --out gwas_qc_snp2

注意--filter-controls这个参数在PLINK 1.9里是默认的，它会只在对照样本中执行HWE检验并过滤，PLINK 2.0需要显式写--hwe 1e-6 include-nonctrl来做全样本过滤，但GWAS标准流程里我仍然建议只过滤对照组。

不过这里要区分一个场景。如果是一个纯病例的case-only研究，比如某些肿瘤GWAS没有正常对照组，那HWE过滤就可以直接用全样本，并把阈值放宽到1e-10甚至1e-12。为什么更严？因为在这种设计里，你没法通过健康对照来预估群体频率，任何HWE偏离都可能是分型错误引起的，宁可错杀也不放过。

4. 完整的QC流程串联：从原始数据到干净数据集

4.1 一套可直接复用的PLINK QC脚本

按前面的思路，我把我实际在项目里用的流程整理成一份可以直接跑的脚本。注意这里使用了--make-bed逐步覆盖中间文件，每一步都有对应的输出日志和统计文件，方便回溯问题。

# 第一步：个体缺失率过滤 plink --bfile gwas_raw --mind 0.05 --make-bed --out step1_mind # 第二步：SNP缺失率过滤 plink --bfile step1_mind --geno 0.05 --make-bed --out step2_geno # 第三步：性别检查（只输出报告，根据结果手动剔除异常样本） plink --bfile step2_geno --check-sex --out step3_sex # 手动查看 step3_sex.sexcheck，选出性别异常的样本ID，放入 file_remove_sex.txt plink --bfile step2_geno --remove file_remove_sex.txt --make-bed --out step4_nosex # 第四步：SNP差异缺失率检验 plink --bfile step4_nosex --test-missing --out step4_diffmiss # 根据输出结果过滤差异缺失显著的位点 plink --bfile step4_nosex --exclude snps_diffmiss.txt --make-bed --out step5_diffmiss # 第五步：MAF过滤 plink --bfile step5_diffmiss --maf 0.05 --make-bed --out step6_maf # 第六步：HWE过滤（只针对对照） plink --bfile step6_maf --hwe 1e-6 --make-bed --out step7_hwe # 第七步：杂合率异常样本过滤 plink --bfile step7_hwe --het --out step8_het # 根据F值均值±3SD计算阈值，手动剔除异常样本 plink --bfile step7_hwe --remove file_remove_het.txt --make-bed --out step9_het # 第八步：亲缘关系过滤 plink --bfile step9_het --indep-pairwise 50 5 0.2 --out step10_prune plink --bfile step9_het --extract step10_prune.prune.in --genome --min 0.2 --out step10_related # 根据PI_HAT结果保留缺失率更低的样本 plink --bfile step9_het --remove file_remove_rel.txt --make-bed --out final_qc

每一步都值得留意一下输出里的统计量——比如--mind跑完后的plink.imiss里F_MISS的最大值分布，--geno跑完后总共有多少个SNP被剔除。这些数字能帮助你判断数据质量到底怎么样，也为后面写methods section收集素材。

有人可能会问：为什么先做个体缺失率，再做SNP缺失率？这个顺序其实有点讲究。如果先做SNP过滤，再回过头看个体缺失率，很多样本的缺失率会明显下降，但你没法确定到底是样本本身好还是凑巧做了过多位点过滤。先清掉劣质样本，再做位点过滤，能让后面的SNP统计更准确。

4.2 质控报告怎么看：关键数值速查

每次跑完QC以后，我都会整理一份简单的质控报告，包含以下核心数字。不要小看这个习惯，它在你写论文methods的时候，以及被审稿人质疑数据质量的时候，能帮你省下大量时间。

当然，这是一套常规阈值，不是铁律。比如做的是超大样本（10万+）的全基因组测序数据，MAF可以放宽到0.001甚至更低，因为测序对低频变异的检出能力远超芯片。反过来，如果你的样本量只有几百，MAF卡0.05都会让你的有效位点数大打折扣，这时候可以考虑卡0.1，虽然损失信息，但统计上更可靠。

我个人的建议是：第一次跑项目时，先按标准阈值全流程跑一遍，再根据每一步的统计量和你的样本量、研究设计去微调。不要一上来就改阈值，因为你还没有掌握数据的整体分布情况。

4.3 lambda值和PCA：QC做完后必须检查的两个指标

质控做完之后，并不代表万事大吉。在跑关联分析之前，我会习惯性地做两次“体检”：

第一，用一组独立于关联分析的中性位点（或者全基因组LD修剪后的位点）去做一个简单的卡方检验，计算基因组膨胀因子lambda。

第二，做PCA主成分分析，看看样本在遗传空间上的分布是否均匀，是否有明显的人群分层。

# LD修剪 plink --bfile final_qc --indep-pairwise 50 5 0.2 --out final_prune # 计算PCA plink --bfile final_qc --extract final_prune.prune.in --pca 10 --out final_pca

跑完PCA以后，把前两个主成分画出来。如果病例和对照在图上明显分成两团，说明人群结构没有消除干净，后面的关联分析里必须把PC1、PC2甚至更多主成分作为协变量放进去。

lambda的计算也很简单。跑一个不加协变量的简单关联：

plink --bfile final_qc --assoc --out qc_check_assoc

然后在R里读取qc_check_assoc.assoc的P值列，计算卡方统计量的中位数除以0.4549（自由度为1的卡方分布中位数）。

p <- read.table("qc_check_assoc.assoc", header=TRUE)$P chi2 <- qchisq(1 - p, 1) lambda <- median(chi2) / qchisq(0.5, 1) print(lambda)

lambda在1到1.1之间，说明QC做得比较干净。如果在1.1到1.15之间，说明还有轻微的人群分层，需要把PCA协变量加入模型。如果大于1.15，我建议你回头检查QC步骤里有没有遗漏问题——比如样本重复、隐藏的亲缘关系、或者批次效应。

5. 常见问题与排查技巧实录

5.1 为什么做完QC后，显著SNP反而变少了

这其实是正常现象，但很多人第一次遇到时会慌张。QC之前看到的“显著”位点，很多都是分型错误或人群分层造成的假阳性。QC把这些位点滤掉以后，真正的关联信号反而更容易浮出水面——虽然p值的绝对值可能没有之前那么惊人了，但可信度大幅提高了。

我在一个代谢疾病项目中就遇到过类似情况。QC之前，曼哈顿图上一眼看过去最少有十几个“基因组显著”的位点。QC之后，只剩两个位点通过了显著性阈值。一开始合作方有点失望，觉得“信号变少了”。但后来这两个位点在独立队列中完美复现，而那些被滤掉的位点，一个都没能复现。这就是QC的价值。

5.2 性别检查报告里大量PROBLEM，是数据有问题还是参数问题

先别急着删样本。性别检查的--check-sex默认阈值在PLINK 1.9里对X染色体杂合率的标准是用0.8和0.2做分界，但不同芯片平台的数据分布会有差异。建议先画出X染色体杂合率的分布图，看看有没有明显的双峰，再决定threshold。

如果只有少量样本落在中间灰色地带，可以手动检查原始表型记录。如果大量样本都处在中间地带，更可能是芯片的性别标记有问题，而不是样本的问题。这时候我建议用--check-sex的--set-hh-missing参数重新跑一次，或者干脆以X染色体杂合率分布图的直观判断为准。

另外有一个容易忽略的点：在某些染色体异常情况下（比如XXY个体、X0个体），性别检查结果也会显示异常。这类样本不一定需要直接剔除，但对后续分析会有影响，建议标记出来单独处理。

5.3 亲缘关系过滤后样本量骤减怎么办

这确实是个现实问题。尤其在队列研究里，同一个家庭的多位成员可能都被纳入了。如果你的项目是做人群关联分析，剔除相关样本是正确做法，因为不剔除会高估检验统计量。

但如果样本量本来就紧张，可以退一步：不直接剔除，而是在混合线性模型（如GEMMA、BOLT-LMM）中通过kinship matrix把亲缘关系作为随机效应校正。这种方式比直接剔除更充分利用数据，前提是你要用对工具、用对模型。但也有个前提——亲缘关系不能过于复杂，如果存在大量一级亲属关系，混合模型也可能无法完全校正，稳妥做法还是剔掉高亲缘样本。

5.4 逻辑回归必须加入PC协变量吗

我的建议是：看PCA结果说话。这里分享一个简单的经验法则，如果你做完PCA后，PC1、PC2或者PC3在病例对照之间有显著差异（用t检验或Mann-Whitney检验看P值），那就必须加PC协变量。

有人喜欢直接用固定数量的PC，比如根据自己的经验取PC1-PC5。也有人会用Traces-Widom检验或者“肘部法则”来选PC数量。我更推荐的是在保证稳定性的前提下，优先看PC解释方差的比例。通常前5到10个PC已经能解释大部分人群结构变异，再加更多的PC对结果影响很小，反而可能吸收掉真正的关联信号。

顺带说一个细节：PLINK的--pca默认不输出特征值比例。如果你要看每个PC的方差解释比例，可以在跑PCA时加--pca 10 header，然后读final_pca.eigenval，用每个特征值除以所有特征值的总和就能得到对应PC的方差解释比例。

5.5 一个踩了很多次的坑：不同染色体数据没有合并就做QC

最后说一个特别基础的坑。有些公共数据是分染色体存放的，比如每一条染色体一个BED文件。如果你直接把所有染色体的BED文件拼在一起做QC，PLINK会默认它们是同一个个体的不同染色体——但前提是FID和IID必须完全一致。

实际操作里，我经常遇到不同染色体的文件里样本ID顺序不一样、甚至ID编码风格不一样的情况。如果直接合并，PLINK会以为某些样本有缺失，反而把染色体间的真实样本串位。

正确的做法是：先把每一条染色体单独过一遍基础QC（--mind、--geno）、确认样本ID格式统一，再用--merge合并。合并后如果遇到Multiple instances of a variant这类报错，说明不同染色体文件里可能有重复位点，先用--exclude去掉重复，再重新merge。

这些经验都是我在实际项目中一个一个踩出来的。现在做GWAS的流程虽然越来越标准化，但数据永远是最磨人的环节。你用的参考基因组版本是什么？芯片是哪个平台？样本有没有做过QC？这些问题每个都藏着坑。

好在PLINK这套QC流程本身足够成熟，只要把每一步的逻辑搞清楚、参数理解透，它就能帮你把数据里90%的隐患找出来。剩下的10%，就需要结合你对这个项目的理解和经验来做判断了。