参考资料
- https://github.com/MareesAT/GWA_tutorial/
- 全基因组关联分析学习资料（GWAS tutorial）
- 论文 A tutorial on conducting genome-wide association studies: Quality control and statistical analysis
继续对数据集进行质控
检查个体杂合性，去除杂合性偏离平均值3个标准差的个体
这一步突然多出来一个inversion.txt文件，怎么来的还不太清楚
使用到的命令是
plink --bfile HapMap_3_r3_9 --exclude inversion.txt --range --indep-pairwise 50 5 0.2 --out indepSNP
plink --bfile HapMap_3_r3_9 --extract indepSNP.prune.in --het --out R_check
对结果进行可视化
het <- read.table("R_check.het", head=TRUE)
head(het)
het$HET_RATE = (het$"N.NM." - het$"O.HOM.")/het$"N.NM."
hist(het$HET_RATE, xlab="Heterozygosity Rate", ylab="Frequency", main= "Heterozygosity Rate")
ggplot(het,aes(x=HET_RATE,y=..count..))+
geom_histogram(bins=90,fill="orange")+
theme_bw()
image.png
选择杂合性超过平均值3个标准差的个体
het_fail <- subset(het, (het$HET_RATE < mean(het$HET_RATE)-3*sd(het$HET_RATE)) | (het$HET_RATE > mean(het$HET_RATE)+3*sd(het$HET_RATE)));
het_fail$HET_DST <- (het_fail$HET_RATE-mean(het$HET_RATE))/sd(het$HET_RATE);
head(het_fail)
write.table(het_fail, "fail-het-qc.txt", row.names=FALSE)
结果中筛选出两个个体，在数据集中去除这两个个体
sed 's/"// g' fail-het-qc.txt | awk '{print$1, $2}'> het_fail_ind.txt
plink --bfile HapMap_3_r3_9 --remove het_fail_ind.txt --make-bed --out HapMap_3_r3_10
检查个体间的亲缘关系
plink --bfile HapMap_3_r3_10 --extract indepSNP.prune.in --genome --min 0.2 --out pihat_min0.2
awk '{ if ($8 >0.9) print $0 }' pihat_min0.2.genome > zoom_pihat.genome
对结果可视化
relatedness <- read.table("pihat_min0.2.genome", header=T)
par(pch=16, cex=1)
with(relatedness,plot(Z0,Z1, xlim=c(0,1), ylim=c(0,1), type="n"))
with(subset(relatedness,RT=="PO") , points(Z0,Z1,col=4))
with(subset(relatedness,RT=="UN") , points(Z0,Z1,col=3))
legend(1,1, xjust=1, yjust=1, legend=levels(relatedness$RT), pch=16, col=c(4,3))
image.png
relatedness_zoom <- read.table("zoom_pihat.genome", header=T)
par(pch=16, cex=1)
with(relatedness_zoom,plot(Z0,Z1, xlim=c(0,0.02), ylim=c(0.98,1), type="n"))
with(subset(relatedness_zoom,RT=="PO") , points(Z0,Z1,col=4))
with(subset(relatedness_zoom,RT=="UN") , points(Z0,Z1,col=3))
legend(0.02,1, xjust=1, yjust=1, legend=levels(relatedness$RT), pch=16, col=c(4,3))
image.png
hist(relatedness[,10],main="Histogram relatedness", xlab= "Pihat")
image.png
plink --bfile HapMap_3_r3_10 --filter-founders --make-bed --out HapMap_3_r3_11
plink --bfile HapMap_3_r3_11 --extract indepSNP.prune.in --genome --min 0.2 --out pihat_min0.2_in_founders
plink --bfile HapMap_3_r3_11 --missing
plink --bfile HapMap_3_r3_11 --remove 0.2_low_call_rate_pihat.txt --make-bed --out HapMap_3_r3_12
至此教程的第一部分就结束了，但是关于亲缘关系这部分代码还完全看不懂！还是先重复完流程再说吧！
欢迎大家关注我的公众号
小明的数据分析笔记本
公众号二维码.jpg
网友评论