前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >全基因组关联分析(GWAS)学习笔记——3.1

全基因组关联分析(GWAS)学习笔记——3.1

作者头像
用户7010445
发布2020-03-03 15:05:04
1.8K0
发布2020-03-03 15:05:04
举报
参考资料
  • 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格式
代码语言:javascript
复制
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文件的行和列分别进行过滤 首先是看一下数据缺失率

代码语言:javascript
复制
plink --bfile HapMap_3_r3_1 --missing

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

代码语言:javascript
复制
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

image.png 首先是根据snp缺失和个体缺失,阈值设置为0.2

代码语言:javascript
复制
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

代码语言:javascript
复制
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 (这个是什么意思还不太明白)
代码语言:javascript
复制
plink --bfile HapMap_3_r3_5 --check-sex 

对结果进行展示

代码语言:javascript
复制
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

image.png

image.png 然后是删除sex discrepancy 的个体

代码语言:javascript
复制
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
代码语言:javascript
复制
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
接下来是统计最小等位基因频率
代码语言:javascript
复制
plink --bfile HapMap_3_r3_7 --freq --out MAF_check

对统计结果进行可视化

代码语言:javascript
复制
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对数据进行过滤

代码语言:javascript
复制
plink --bfile HapMap_3_r3_7 --maf 0.05 --make-bed --out HapMap_3_r3_8
接下来检测不符合哈迪温伯格定律的snp
代码语言:javascript
复制
plink --bfile HapMap_3_r3_8 --hardy

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

代码语言:javascript
复制
awk '{ if ($9 <0.00001) print $0 }' plink.hwe > plinkzoomhwe.hwe

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

代码语言:javascript
复制
Rscript --no-save hwe.R 

image.png

image.png

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

代码语言:javascript
复制
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

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

本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2019-12-27,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 小明的数据分析笔记本 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 第一步根据vcf文件得到了三个文件,分别是.bed;.fam;.bim
  • 第二步是对数据进行质控
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档