今天分享gwas教程中的性别质控。
这个章节,主要是人类性别的信息的质控,主要是根据性染色上SNP的比值,判断性别,然后把性别错误的个体去掉或者更改性别信息。对其它物种参考意义不大,因为在动物中一般把性别信息的SNP去掉,植物中一般都是雌雄同体的,不涉及到这个问题,之所以会有这一篇,是因为原文中有这个信息,而且plink 也有--check-sex的参数,所以操作一下,留下笔记。
原理:检查性别差异。先验信息,女性的受试者的F值必须小于0.2,男性的受试者的F值必须大于0.8。这个F值是基于X染色体近交(纯合子)估计。不符合这些要求的受试者被PLINK标记为“PROBLEM”。
上一步,去掉缺失信息后,现在有文件是过滤缺失后的文件:
ls HapMap_3_r3_3.*
HapMap_3_r3_3.bed HapMap_3_r3_3.fam
HapMap_3_r3_3.bim HapMap_3_r3_3.hh
HapMap_3_r3_3.log
检查性别冲突
plink --bfile HapMap_3_r3_3 --check-sex
结果文件:plink.sexcheck
第一列为家系ID
第二列为个体ID
第三列为系谱中的性别
第四列为SNP推断的性别
第五列是否正常
第六列为F值。
使用 R 语言作图
# 先要设置工作路径
gender = read.table("plink.sexcheck", header=T,as.is=T)
# 显示图片
hist(gender[,6],main="Gender", xlab="F")
male=subset(gender, gender$PEDSEX==1)
hist(male[,6],main="Men",xlab="F")
female=subset(gender, gender$PEDSEX==2)
hist(female[,6],main="Women",xlab="F")
# 保存图片
pdf("Gender_check.pdf")
hist(gender[,6],main="Gender", xlab="F")
dev.off()
pdf("Men_check.pdf")
male=subset(gender, gender$PEDSEX==1)
hist(male[,6],main="Men",xlab="F")
dev.off()
pdf("Women_check.pdf")
female=subset(gender, gender$PEDSEX==2)
hist(female[,6],main="Women",xlab="F")
dev.off()
结果:
图中可以看出,woman中,大部分都是小于0.2,有一个为1,应该是错误的ID。
提取错误的 ID
我们使用grep过滤一下:根据STATUS列,如果有问题的话,为“PROBLEM”,我们可以根据这个关键词将有问题的行打印出来。这里,可以通过在git终端下执行:
代码:
grep "PROBLEM" plink.sexcheck
将相关错误的ID提取出来(家系ID,个体ID),之所以提取家系ID和个体ID,因为plink有参数remove可以根据ID进行筛选。
使用 remove 去掉个体
plink --bfile HapMap_3_r3_3 --remove sex_discrepancy.txt --make-bed --
out HapMap_3_r3_6
当然,你也可以对个体进行判定填充,这是用就可以实现,
plink --bfile HapMap_3_r3_3 --impute-sex --recode --out test1
再对test1进行check-sex,就会发现性别信息没有问题了,即:test1的plink文件,已经将错误的性别信息,根据基因型的检测结果矫正过来了。
过滤的参数去掉个体或者SNP,关键词不一样,容易混淆,这里总结一下。
保留或去掉个体
--keep
--remove
--keep-fam
--remove-fam
保留或去掉SNP:
--extract ['range']
--exclude ['range']