这个章节,主要是人类性别的信息的质控,主要是根据性染色上SNP的比值,判断性别,然后把性别错误的个体去掉或者更改性别信息。对其它物种参考意义不大,因为在动物中一般把性别信息的SNP去掉,植物中一般都是雌雄同体的,不涉及到这个问题,之所以会有这一篇,是因为原文中有这个信息,而且plink 也有​​--check-sex​​的参数,所以操作一下,留下笔记。


原理:

检查性别差异。先验信息,女性的受试者的F值必须小于0.2,男性的受试者的F值必须大于0.8。这个F值是基于X染色体近交(纯合子)估计。不符合这些要求的受试者被PLINK标记为“PROBLEM”。

上一步,去掉缺失信息后,现在有文件是过滤缺失后的文件:

HapMap_3_r3_5.bed  HapMap_3_r3_5.fam  HapMap_3_r3_5.irem
HapMap_3_r3_5.bim HapMap_3_r3_5.hh HapMap_3_r3_5.log

1. 检查性别冲突

plink --bfile HapMap_3_r3_5 --check-sex

结果文件:​​plink.sexcheck​

第一列为家系ID,第二列为个体ID,第三列为系谱中的性别,第四列为SNP推断的性别,第五列是否正常,第六列为F值。

笔记 | GWAS 操作流程2-2:性别质控_r语言

使用R语言作图:

gender <- read.table("plink.sexcheck", header=T,as.is=T)

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()

笔记 | GWAS 操作流程2-2:性别质控_r语言_02

笔记 | GWAS 操作流程2-2:性别质控_r语言_03

图中可以看出,woman中,大部分都是小于0.2,有一个为1,应该是错误的ID。

2. 提取错误的ID

我们使用​​awk​​过滤一下:

根据STATUS列,如果有问题的话,为“PROBLEM”,我们可以根据这个关键词将有问题的行打印出来。

grep "PROBLEM" plink.sexcheck

​1349 NA10854 2 1 PROBLEM 0.99​

可以看出,个体NA10854是有问题的。

将相关错误的ID提取出来(家系ID,个体ID),之所以提取家系ID和个体ID,因为plink有参数​​remove​​可以根据ID进行筛选。

grep 'PROBLEM' plink.sexcheck | awk '{print $1,$2}' >sex_discrepancy.txt

我们将结果保存在​​sex_discrepancy.txt​​。

3. 使用​​remove​​去掉个体

plink --bfile HapMap_3_r3_5 --remove sex_discrepancy.txt --make-bed --out HapMap_3_r3_6

当然,你也可以对个体进行判定填充,这是用​​--impute-sex​​就可以实现,这样的话那个错误的个体会根据统计量更改性别信息。这里我们选择的是删掉这个个体。

4. 过滤的关键词

去掉个体或者SNP,关键词不一样,容易混淆,这里总结一下。

保留或去掉个体:

--keep <filename>
--remove <filename>

--keep-fam <filename>
--remove-fam <filename>

保留或去掉SNP:

--extract ['range'] <filename>
--exclude ['range'] <filename>