在上一篇里我们已经把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效应量OR2.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_related3. 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精确检验把P0.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的时候以及被审稿人质疑数据质量的时候能帮你省下大量时间。指标阈值判断标准样本缺失率F_MISS 0.05超过则剔除样本SNP缺失率GENO 0.05超过则剔除位点性别不一致STATUCPROBLEM直接剔除杂合率偏离均值±3SD超出则剔除亲缘关系PI_HAT 0.2超过则二选一保留MAF 0.05低于则剔除HWE对照P 1e-6低于则剔除差异缺失率Fisher P 1e-4低于则剔除当然这是一套常规阈值不是铁律。比如做的是超大样本10万的全基因组测序数据MAF可以放宽到0.001甚至更低因为测序对低频变异的检出能力远超芯片。反过来如果你的样本量只有几百MAF卡0.05都会让你的有效位点数大打折扣这时候可以考虑卡0.1虽然损失信息但统计上更可靠。我个人的建议是第一次跑项目时先按标准阈值全流程跑一遍再根据每一步的统计量和你的样本量、研究设计去微调。不要一上来就改阈值因为你还没有掌握数据的整体分布情况。4.3 lambda值和PCAQC做完后必须检查的两个指标质控做完之后并不代表万事大吉。在跑关联分析之前我会习惯性地做两次“体检”第一用一组独立于关联分析的中性位点或者全基因组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, headerTRUE)$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文件拼在一起做QCPLINK会默认它们是同一个个体的不同染色体——但前提是FID和IID必须完全一致。实际操作里我经常遇到不同染色体的文件里样本ID顺序不一样、甚至ID编码风格不一样的情况。如果直接合并PLINK会以为某些样本有缺失反而把染色体间的真实样本串位。正确的做法是先把每一条染色体单独过一遍基础QC--mind、--geno、确认样本ID格式统一再用--merge合并。合并后如果遇到Multiple instances of a variant这类报错说明不同染色体文件里可能有重复位点先用--exclude去掉重复再重新merge。这些经验都是我在实际项目中一个一个踩出来的。现在做GWAS的流程虽然越来越标准化但数据永远是最磨人的环节。你用的参考基因组版本是什么芯片是哪个平台样本有没有做过QC这些问题每个都藏着坑。好在PLINK这套QC流程本身足够成熟只要把每一步的逻辑搞清楚、参数理解透它就能帮你把数据里90%的隐患找出来。剩下的10%就需要结合你对这个项目的理解和经验来做判断了。