参考资料
- 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
得把教程对应的文章多看几遍！ 今天就到这里了。 未完待续！
