企业官网建设流程全解析

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.lmiss：SNP缺失率统计
• miss_report.hh：临时文件
• miss_report.log：日志文件

关键指标解读：

• F_MISS（在.imiss中）：样本的基因型缺失比例
• F_MISS（在.lmiss中）：SNP在所有样本中的缺失比例

2.2 可视化缺失率分布

使用R脚本可视化缺失率分布能更直观地识别异常值：

# 读取缺失率数据 ind_miss <- read.table("miss_report.imiss", header=TRUE) snp_miss <- read.table("miss_report.lmiss", header=TRUE) # 绘制样本缺失率直方图 pdf("sample_missingness.pdf") hist(ind_miss$F_MISS, breaks=50, col="lightblue", main="Sample Missingness", xlab="Missing Rate") abline(v=0.05, col="red", lty=2) # 常用阈值线 dev.off() # 绘制SNP缺失率直方图 pdf("snp_missingness.pdf") hist(snp_miss$F_MISS, breaks=50, col="lightgreen", main="SNP Missingness", xlab="Missing Rate") abline(v=0.02, col="red", lty=2) # 常用阈值线 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染色体杂合度是判断遗传性别的重要指标：

• 女性（XX）：X染色体杂合度较高（F值<0.2）
• 男性（XY）：X染色体杂合度较低（F值>0.8）

3.1 执行性别检查

plink --bfile data_clean_miss --check-sex --out sex_check

输出文件sex_check.sexcheck包含关键列：

• PEDSEX：.fam文件中报告的性别（1=男，2=女）
• SNPSEX：基于基因型推断的性别
• STATUS：一致性状态（OK/PROBLEM）
• F：X染色体近交系数

3.2 可视化性别检查结果

sex <- read.table("sex_check.sexcheck", header=TRUE) pdf("gender_validation.pdf", width=10, height=5) par(mfrow=c(1,2)) # 全体样本F值分布 hist(sex$F, breaks=50, col="gray", main="All Samples", xlab="F value") abline(v=c(0.2, 0.8), col=c("pink", "blue"), lty=2) # 按报告性别分组展示 male <- subset(sex, PEDSEX==1) female <- subset(sex, PEDSEX==2) boxplot(list(Male=male$F, Female=female$F), col=c("blue","pink"), main="F value by Reported Gender", ylab="F value") abline(h=0.5, col="red", lty=2) 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阈值选择需考虑样本量：

4.1 计算和可视化MAF分布

# 计算MAF plink --bfile data_clean_sex --freq --out maf_check

maf <- read.table("maf_check.frq", header=TRUE) pdf("maf_distribution.pdf") hist(maf$MAF, breaks=50, col="lightyellow", main="MAF Distribution", xlab="Minor Allele Frequency") abline(v=0.01, col="red", lty=2) 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($1>=1 && $1<=22) 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_maf

5. Hardy-Weinberg平衡检验

HWE检验识别偏离遗传平衡的SNP，这些偏离可能源于自然选择、基因分型错误或群体分层。

5.1 执行HWE检验

plink --bfile data_clean_maf --hardy --out hwe_check

5.2 分析HWE结果

hwe <- read.table("hwe_check.hwe", header=TRUE) pdf("hwe_analysis.pdf", width=10, height=5) par(mfrow=c(1,2)) # 全部SNP的p值分布 hist(hwe$P, breaks=50, col="lightgray", main="HWE p-value Distribution", xlab="p-value") # 显著偏离SNP的p值分布（p<1e-5） hwe_sig <- subset(hwe, P<1e-5) hist(hwe_sig$P, breaks=20, col="pink", main="Significant HWE Deviations (p<1e-5)", xlab="p-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_check

6.2 杂合率可视化与分析

het <- read.table("het_check.het", header=TRUE) het$HET_RATE <- (het$N.NM. - het$O.HOM.)/het$N.NM. pdf("heterozygosity_analysis.pdf", width=10, height=5) par(mfrow=c(1,2)) # 杂合率分布 hist(het$HET_RATE, breaks=30, col="lightblue", main="Sample Heterozygosity", xlab="Heterozygosity Rate") # 杂合率与缺失率的关系 miss <- read.table("miss_report.imiss", header=TRUE) merge_data <- merge(het, miss, by=c("FID","IID")) plot(merge_data$HET_RATE, merge_data$F_MISS, pch=16, col="blue", xlab="Heterozygosity Rate", ylab="Missing Rate", main="Heterozygosity 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.names=FALSE, quote=FALSE)

移除异常样本：

plink --bfile data_clean_hwe --remove het_outliers.txt --make-bed --out data_clean_het

7. 亲缘关系分析与样本去重

隐性亲缘关系会导致假阳性关联，需要通过IBD分析识别。

7.1 计算亲缘关系

plink --bfile data_clean_het --extract prune.prune.in --genome --min 0.2 --out ibd_check

7.2 分析亲缘关系结果

ibd <- read.table("ibd_check.genome", header=TRUE) pdf("ibd_analysis.pdf", width=8, height=6) plot(ibd$Z0, ibd$Z1, col=as.numeric(factor(ibd$RT)), pch=16, xlab="Z0 (IBD=0)", ylab="Z1 (IBD=1)", main="Identity by Descent (IBD) Analysis") legend("topright", legend=levels(factor(ibd$RT)), pch=16, col=1:nlevels(factor(ibd$RT))) dev.off()

7.3 处理相关个体

对于每对相关个体，通常保留数据质量较高者：

# 假设已生成需移除样本列表ibd_remove.txt plink --bfile data_clean_het --remove ibd_remove.txt --make-bed --out data_final

8. 质控后数据评估与归档

完成所有质控步骤后，应对最终数据集进行全面评估：

# 生成最终统计报告 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数量变化：

最终数据集已准备好用于后续的群体分层分析和关联研究。建议保存完整的质控脚本和中间文件，以确保研究可重复性。在实际项目中，可能需要根据数据特点调整各步骤的阈值参数，平衡数据质量与信息保留之间的权衡。