全基因组关联分析（GWAS）学习笔记——3.1

参考资料

https://github.com/MareesAT/GWA_tutorial/
全基因组关联分析学习资料（GWAS tutorial）
论文 A tutorial on conducting genome-wide association studies: Quality control and statistical analysis

数据集

HapMap_3_r3_1.bed
HapMap_3_r3_1.fam
HapMap_3_r3_1.bim

数据集的具体信息还没太搞懂

将上面三个数据转化为vcf格式

plink --bfile HapMap_3_r3_1 --recode vcf bgz --out gwasPra

参考
Create VCF from .bim, .bed and .fam files

第一步根据vcf文件得到了三个文件，分别是.bed；.fam；.bim

这个教程里是直接有这三个文件

第二步是对数据进行质控

根据snp的缺失率和个体标记的缺失率进行过滤

我目前理解为根据vcf文件的行和列分别进行过滤
首先是看一下数据缺失率

plink --bfile HapMap_3_r3_1 --missing

输出结果中plink.imiss文件是个体标记的缺失率；plink.lmiss是每个标记个体的缺失率
原教程中提供了R脚本对这两个文件使用直方图进行可视化，我这里选择ggplot2对结果进行可视化

indmiss<-read.table(file="plink.imiss",header=T)
snpmiss<-read.table(file="plink.lmiss",header=T)
head(indmiss)
head(snpmiss)
library(ggplot2)
ggplot(indmiss,aes(x=F_MISS,y=..count..))+
  geom_histogram(bins=90,fill="darkgreen")+
  theme_bw()
options(scipen = 999)
ggplot(snpmiss,aes(x=F_MISS,y=..count..))+
  geom_histogram(bins=30,fill="darkgreen")+
  theme_bw()

image.png

首先是根据snp缺失和个体缺失，阈值设置为0.2

plink --bfile HapMap_3_r3_1 --geno 0.2 --make-bed --out HapMap_3_r3_2
plink --bfile HapMap_3_r3_2 --mind 0.2 --make-bed --out HapMap_3_r3_3

再把阈值设置为0.02

plink --bfile HapMap_3_r3_3 --geno 0.02 --make-bed --out HapMap_3_r3_4
plink --bfile HapMap_3_r3_4 --mind 0.02 --make-bed --out HapMap_3_r3_5

不太明白为什么这块分两步进行，直接将阈值设置为0.02可以吗？

接下来是检查 sex discrepancy (这个是什么意思还不太明白)

plink --bfile HapMap_3_r3_5 --check-sex

对结果进行展示

gender <- read.table("plink.sexcheck", header=T,as.is=T)
head(gender)
colnames(gender)<-c(colnames(gender)[1:5],"F6")
ggplot(gender,aes(x=F6,y=..count..))+
  geom_histogram(bins=30,fill="darkgreen")+
  theme_bw()
male<-subset(gender, gender$PEDSEX==1)
colnames(male)<-c(colnames(gender)[1:5],"F6")
ggplot(male,aes(x=F6,y=..count..))+
  geom_histogram(bins=30,fill="darkgreen")+
  theme_bw()
female<-subset(gender, gender$PEDSEX==2)
colnames(female)<-c(colnames(gender)[1:5],"F6")
ggplot(female,aes(x=F6,y=..count..))+
  geom_histogram(bins=30,fill="darkgreen")+
  theme_bw()

image.png

然后是删除sex discrepancy 的个体

grep "PROBLEM" plink.sexcheck | awk '{print$1,$2}' >  sex_discrepancy.txt
plink --bfile HapMap_3_r3_5 --remove sex_discrepancy.txt --make-bed --out HapMap_3_r3_6

接下来是只保留常染色体SNP

awk '{ if ($1 >= 1 && $1 <= 22) print $2 }' HapMap_3_r3_6.bim > snp_1_22.txt
plink --bfile HapMap_3_r3_6 --extract snp_1_22.txt --make-bed --out HapMap_3_r3_7

接下来是统计最小等位基因频率

plink --bfile HapMap_3_r3_7 --freq --out MAF_check

对统计结果进行可视化

maf_freq<-read.table("MAF_check.frq",header=T,as.is=T)
head(maf_freq)
ggplot(maf_freq,aes(x=MAF,y=..count..))+
  geom_histogram(bins=90,fill="darkgreen")+
  theme_bw()

image.png

最小等位基因频率阈值设置为0.05对数据进行过滤

plink --bfile HapMap_3_r3_7 --maf 0.05 --make-bed --out HapMap_3_r3_8

接下来检测不符合哈迪温伯格定律的snp

plink --bfile HapMap_3_r3_8 --hardy

接下来的步骤操作明白，但是不太明白背后的意义

awk '{ if ($9 <0.00001) print $0 }' plink.hwe > plinkzoomhwe.hwe

这一步的数据文件有点大，直接按照教程中的命令在服务器上运行教程提供的脚本

Rscript --no-save hwe.R

image.png

对数据进行过滤，但是背后的意义不太明白

plink --bfile HapMap_3_r3_8 --hwe 1e-6 --make-bed --out HapMap_hwe_filter_step1
plink --bfile HapMap_hwe_filter_step1 --hwe 1e-10 --hwe-all --make-bed --out HapMap_3_r3_9

得把教程对应的文章多看几遍！
今天就到这里了。
未完待续！

欢迎大家关注我的公众号
小明的数据分析笔记本

公众号二维码.jpg

最后编辑于：2019.12.27 22:18:31