上一篇写的是：重复一篇文献的GWAS（一）：基因型数据整理
这一篇主要讲GWAS的核心流程和画图。
软件使用参考：
使用GEMMA进行复杂性状全基因组关联分析(GWAS)
https://www.jianshu.com/p/d31404620c9b
github
https://github.com/genetics-statistics/GEMMA
https://github.com/genetics-statistics/GEMMA/blob/master/doc/manual.pdf
1. 将ped/map转换为fam/bed/bim
plink --file 172sample_maf_filter_snpID_LD_filter --make-bed --out clean_snp
#172sample_maf_filter_snpID_LD_filter是map，ped文件的前缀，会生成以clean_snp为前缀以bed, bim, fam结尾的三个文件
以上几个文件格式需要熟悉一下。
2. 修改表型文件
前面没有输入表型信息，172sample_maf_filter_snpID_LD_filter.ped文件和clean_snp.fam文件的第6列都是-9，默认值。所有样本的表型值在文献补充材料中都能找到，注意样本顺序替换表型值即可。
3. 生成kinship（亲缘关系）矩阵
nohup gemma -bfile clean_snp -gk -o kinship &
#gk参数可以指定kinship矩阵的类型，默认是1
$ tree ./output/
./output/
├── kinship.cXX.txt
└── kinship.log.txt
4. 单变量的混合线性模型
文献描述：A GWAS for seed area was conducted using a univariate mixed linear model method (LMM) in GEMMA software using the default parameters
混合线性模型既考虑了群体分层，也考虑了样本之间的关系。
nohup gemma -bfile clean_snp -lmm -k ./output/kinship.cXX.txt -o GWAS &
$ tree output/
output/
├── GWAS.assoc.txt
├── GWAS.log.txt
├── kinship.cXX.txt
└── kinship.log.txt
$ lsx GWAS.assoc.txt | head -n 5
chr rs ps n_miss allele1 allele0 af beta se logl_H1 l_remle p_wald
1 1_63 63 0 C T 0.064 -1.706603e-03 3.725438e-03 3.910054e+02 5.240987e+00 6.474695e-01
1 1_92 92 0 C A 0.401 4.789343e-04 1.974829e-03 3.909931e+02 6.030706e+00 8.086701e-01
1 1_138 138 0 C T 0.064 8.253728e-03 4.108313e-03 3.929843e+02 4.610281e+00 4.611623e-02
1 1_266 266 0 A G 0.064 4.685363e-03 4.019586e-03 3.916598e+02 5.506291e+00 2.453956e-01
上述GWAS.assoc.txt文件就是我们需要的结果，可以用来做判断以及画图。但是p值应该如何确定呢？
5. p值选择
Bonferroni-corrected：0.05除以标记数，以此作为显著水平线。原文是1.5E-07，我这里是1.49E-07。换句话说就是，原来的p值乘以标记数仍然能小于0.05的就符合要求。
到这儿我看了一下，和原文结果有一些出入，可能原因是：19个样本缺失导致；前面用LD过滤，可能在同一个LD block中原文留下的是A点，我留下的是B点，导致位置有偏差。
6. 换一个工具跑一下
使用rMVP包
> library(rMVP)
> MVP.Data(fileBed="clean_snp",filePhe="clean_value.txt",fileKin=TRUE,filePC=TRUE,out="mvp")
#clean_snp为plink二进制文件的前缀，clean_value.txt为样本与表型值的对应关系
会生成GWAS分析的所有数据文件，以及PCA作图的文件
$ ls mvp.*
mvp.geno.bin mvp.geno.desc mvp.geno.ind mvp.geno.map mvp.kin.bin mvp.kin.desc mvp.pc.bin mvp.pc.desc mvp.phe
读取这些文件
> pheno <- read.table("mvp.phe", header = TRUE)
> geno <- attach.big.matrix("mvp.geno.desc")
> map <- read.table("mvp.geno.map", header = TRUE)
> Kinship <- attach.big.matrix("mvp.kin.desc")
> Covariates <- attach.big.matrix("mvp.pc.desc")
运行程序，此处计算方差组分的方法是"EMMA"
> MVP(phe=pheno,geno=geno,map=map,K=Kinship,CV.MLM=Covariates,priority="speed",vc.method="EMMA",method=c("MLM"))
我将两次（一次是严格按照文献步骤--GEMMA软件，一次是用MVP包--EMMA算法）的结果取top 100 SNPs看交集，可以看到70%的重叠。我认为这个重叠还算合理，所以觉得我前面的分析问题不大，和原文献有差别的原因可能是：①19样本缺失；②根据LD过滤后余下的SNP有差异（因为对于原文找到的top SNPs来看，在我过滤得到的335.6K的SNP中基本找不到）
尽管文献的结果没有重复出来，形式还是要走的嘛！
接下来画图。
7. 画图
文献中，关于GWAS这一块的图有：①所有样本地理位置与种子大小的图；②样本表型值的概率密度图与直方图；③曼哈顿图；④top SNPs的累积效应图。
①、④还没有画过，有时间补一下；②之前写过——正态分布检验，QQ图的原理也在这一篇中写过；这里画一下曼哈顿图和QQ图。
曼哈顿图、QQ图
$ head -n 4 for_mvp.txt
rs chr pos seed_size
1_63 1 63 6.474695e-01
1_92 1 92 8.086701e-01
1_138 1 138 4.611623e-02
> library(rMVP)
> df <- read.table("for_mvp.txt",header=T,sep="\t")
> MVP.Report(df,plot.type="m",LOG10=TRUE,threshold=0.05/nrow(df),threshold.lty=2,threshold.lwd=1,threshold.col="grey",amplify=TRUE,col=c("dodgerblue4","deepskyblue"),signal.col="red",cex=0.7,chr.den.col=NULL,file="pdf",memo="",dpi=600)
> MVP.Report(df,plot.type="q",conf.int.col=NULL,box=TRUE,file="pdf",dpi=600)
网友评论