ADNI数据库SNP数据质控实战:用Plink一步步搞定GWAS前的数据清洗(附R脚本) ADNI数据库SNP数据质控全流程解析从原始数据到GWAS-ready数据集在基因组关联分析(GWAS)研究中数据质量直接决定了研究结果的可靠性。ADNI数据库作为阿尔茨海默病研究的重要资源其SNP数据需要经过严格的质控流程才能用于后续分析。本文将详细介绍如何使用Plink和R语言完成从原始数据到高质量数据集的完整质控流程。1. 质控流程概述与数据准备SNP数据质控是GWAS研究不可逾越的关键步骤。一个完整的质控流程通常包括六个核心环节缺失率检查、性别校验、MAF过滤、HWE检验、杂合率控制和亲缘关系排查。这套流程能够系统性地识别和剔除低质量数据确保后续分析结果的可靠性。ADNI数据库提供的SNP数据通常以PLINK二进制格式存储包含三个核心文件.bed文件存储基因型数据.bim文件存储SNP标记信息.fam文件存储样本信息在开始质控前建议先创建独立的工作目录并将原始数据文件统一命名如data.bed、data.bim、data.fam这能有效避免后续命令中的文件路径错误。mkdir gwas_qc cd gwas_qc cp /path/to/original_data.* . rename s/original_/data/ original_*2. 缺失率分析与过滤策略缺失率是评估数据质量的首要指标。高缺失率可能源于实验技术问题或样本/DNA质量不佳。我们分两步进行缺失率分析先检查原始数据的缺失情况再逐步过滤不合格的SNP和样本。2.1 初始缺失率检查使用Plink生成缺失率报告plink --bfile data --missing --out miss_report此命令会生成四个文件miss_report.imiss样本缺失率统计miss_report.lmissSNP缺失率统计miss_report.hh临时文件miss_report.log日志文件关键指标解读F_MISS在.imiss中样本的基因型缺失比例F_MISS在.lmiss中SNP在所有样本中的缺失比例2.2 可视化缺失率分布使用R脚本可视化缺失率分布能更直观地识别异常值# 读取缺失率数据 ind_miss - read.table(miss_report.imiss, headerTRUE) snp_miss - read.table(miss_report.lmiss, headerTRUE) # 绘制样本缺失率直方图 pdf(sample_missingness.pdf) hist(ind_miss$F_MISS, breaks50, collightblue, mainSample Missingness, xlabMissing Rate) abline(v0.05, colred, lty2) # 常用阈值线 dev.off() # 绘制SNP缺失率直方图 pdf(snp_missingness.pdf) hist(snp_miss$F_MISS, breaks50, collightgreen, mainSNP Missingness, xlabMissing Rate) abline(v0.02, colred, lty2) # 常用阈值线 dev.off()2.3 分阶段缺失率过滤推荐采用两阶段过滤策略先宽松后严格避免一次性过滤过多数据第一阶段过滤宽松阈值# 过滤缺失率20%的SNP plink --bfile data --geno 0.2 --make-bed --out data_step1 # 过滤缺失率20%的样本 plink --bfile data_step1 --mind 0.2 --make-bed --out data_step2第二阶段过滤严格阈值# 过滤缺失率2%的SNP plink --bfile data_step2 --geno 0.02 --make-bed --out data_step3 # 过滤缺失率2%的样本 plink --bfile data_step3 --mind 0.02 --make-bed --out data_clean_miss注意必须按顺序先过滤SNP再过滤样本。若顺序颠倒可能残留高缺失率SNP。3. 性别校验与不一致处理性别校验是通过比较遗传性别与报告性别的一致性来识别样本错误或污染。X染色体杂合度是判断遗传性别的重要指标女性XXX染色体杂合度较高F值0.2男性XYX染色体杂合度较低F值0.83.1 执行性别检查plink --bfile data_clean_miss --check-sex --out sex_check输出文件sex_check.sexcheck包含关键列PEDSEX.fam文件中报告的性别1男2女SNPSEX基于基因型推断的性别STATUS一致性状态OK/PROBLEMFX染色体近交系数3.2 可视化性别检查结果sex - read.table(sex_check.sexcheck, headerTRUE) pdf(gender_validation.pdf, width10, height5) par(mfrowc(1,2)) # 全体样本F值分布 hist(sex$F, breaks50, colgray, mainAll Samples, xlabF value) abline(vc(0.2, 0.8), colc(pink, blue), lty2) # 按报告性别分组展示 male - subset(sex, PEDSEX1) female - subset(sex, PEDSEX2) boxplot(list(Malemale$F, Femalefemale$F), colc(blue,pink), mainF value by Reported Gender, ylabF value) abline(h0.5, colred, lty2) dev.off()3.3 处理性别不一致样本发现性别不一致样本时有两种处理方式方案1直接删除不一致样本grep PROBLEM sex_check.sexcheck | awk {print $1,$2} sex_discrepancy.txt plink --bfile data_clean_miss --remove sex_discrepancy.txt --make-bed --out data_clean_sex方案2用遗传性别修正报告性别plink --bfile data_clean_miss --impute-sex --make-bed --out data_clean_sex提示对于阿尔茨海默病研究建议保留性别信息用于后续分析中的协变量调整。4. 最小等位基因频率(MAF)过滤MAF过滤能去除低频变异提高统计效力并减少多重检验负担。MAF阈值选择需考虑样本量样本量推荐MAF阈值1,0000.051,000-5,0000.015,0000.0054.1 计算和可视化MAF分布# 计算MAF plink --bfile data_clean_sex --freq --out maf_checkmaf - read.table(maf_check.frq, headerTRUE) pdf(maf_distribution.pdf) hist(maf$MAF, breaks50, collightyellow, mainMAF Distribution, xlabMinor Allele Frequency) abline(v0.01, colred, lty2) dev.off()4.2 执行MAF过滤plink --bfile data_clean_sex --maf 0.01 --make-bed --out data_clean_maf对于特定研究可先提取常染色体SNP再进行MAF过滤# 提取1-22号染色体SNP awk {if($11 $122) print $2} data_clean_sex.bim autosome_snps.txt plink --bfile data_clean_sex --extract autosome_snps.txt --make-bed --out data_autosome plink --bfile data_autosome --maf 0.01 --make-bed --out data_clean_maf5. Hardy-Weinberg平衡检验HWE检验识别偏离遗传平衡的SNP这些偏离可能源于自然选择、基因分型错误或群体分层。5.1 执行HWE检验plink --bfile data_clean_maf --hardy --out hwe_check5.2 分析HWE结果hwe - read.table(hwe_check.hwe, headerTRUE) pdf(hwe_analysis.pdf, width10, height5) par(mfrowc(1,2)) # 全部SNP的p值分布 hist(hwe$P, breaks50, collightgray, mainHWE p-value Distribution, xlabp-value) # 显著偏离SNP的p值分布p1e-5 hwe_sig - subset(hwe, P1e-5) hist(hwe_sig$P, breaks20, colpink, mainSignificant HWE Deviations (p1e-5), xlabp-value) dev.off()5.3 基于HWE的过滤针对不同分析采用不同阈值# 对对照组过滤宽松 plink --bfile data_clean_maf --hwe 1e-6 --make-bed --out data_hwe_filtered # 对全体样本过滤严格 plink --bfile data_hwe_filtered --hwe 1e-10 --hwe-all --make-bed --out data_clean_hwe注意疾病研究中病例组中疾病相关SNP可能自然偏离HWE因此需谨慎设置阈值。6. 杂合率分析与异常样本识别杂合率异常可能提示样本污染、近亲繁殖或基因分型问题。我们将通过以下步骤识别异常样本。6.1 计算杂合率# 先进行LD修剪 plink --bfile data_clean_hwe --indep-pairwise 50 5 0.2 --out prune # 基于修剪后SNP计算杂合率 plink --bfile data_clean_hwe --extract prune.prune.in --het --out het_check6.2 杂合率可视化与分析het - read.table(het_check.het, headerTRUE) het$HET_RATE - (het$N.NM. - het$O.HOM.)/het$N.NM. pdf(heterozygosity_analysis.pdf, width10, height5) par(mfrowc(1,2)) # 杂合率分布 hist(het$HET_RATE, breaks30, collightblue, mainSample Heterozygosity, xlabHeterozygosity Rate) # 杂合率与缺失率的关系 miss - read.table(miss_report.imiss, headerTRUE) merge_data - merge(het, miss, byc(FID,IID)) plot(merge_data$HET_RATE, merge_data$F_MISS, pch16, colblue, xlabHeterozygosity Rate, ylabMissing Rate, mainHeterozygosity vs Missingness) dev.off()6.3 识别和移除异常样本识别杂合率偏离均值±3SD的样本het$Z_SCORE - scale(het$HET_RATE) outliers - subset(het, abs(Z_SCORE) 3) write.table(outliers[,1:2], het_outliers.txt, row.namesFALSE, quoteFALSE)移除异常样本plink --bfile data_clean_hwe --remove het_outliers.txt --make-bed --out data_clean_het7. 亲缘关系分析与样本去重隐性亲缘关系会导致假阳性关联需要通过IBD分析识别。7.1 计算亲缘关系plink --bfile data_clean_het --extract prune.prune.in --genome --min 0.2 --out ibd_check7.2 分析亲缘关系结果ibd - read.table(ibd_check.genome, headerTRUE) pdf(ibd_analysis.pdf, width8, height6) plot(ibd$Z0, ibd$Z1, colas.numeric(factor(ibd$RT)), pch16, xlabZ0 (IBD0), ylabZ1 (IBD1), mainIdentity by Descent (IBD) Analysis) legend(topright, legendlevels(factor(ibd$RT)), pch16, col1:nlevels(factor(ibd$RT))) dev.off()7.3 处理相关个体对于每对相关个体通常保留数据质量较高者# 假设已生成需移除样本列表ibd_remove.txt plink --bfile data_clean_het --remove ibd_remove.txt --make-bed --out data_final8. 质控后数据评估与归档完成所有质控步骤后应对最终数据集进行全面评估# 生成最终统计报告 plink --bfile data_final --missing --out final_miss plink --bfile data_final --freq --out final_freq plink --bfile data_final --hardy --out final_hwe plink --bfile data_final --het --out final_het # 归档质控流程 tar -czvf gwas_qc_pipeline.tar.gz *.log *.pdf *.txt创建质控流程文档记录各步骤的样本和SNP数量变化质控步骤剩余样本数剩余SNP数过滤原因初始数据8122,379,855-缺失率过滤7921,241,966高缺失率性别校验7921,241,966性别不一致MAF过滤7921,073,372MAF 0.01HWE过滤7921,241,966HWE偏离(p1e-10)杂合率过滤7711,241,966杂合率异常亲缘关系过滤7651,241,966隐性亲缘关系(PI_HAT0.2)最终数据集已准备好用于后续的群体分层分析和关联研究。建议保存完整的质控脚本和中间文件以确保研究可重复性。在实际项目中可能需要根据数据特点调整各步骤的阈值参数平衡数据质量与信息保留之间的权衡。