gemma 软件对数量性状 混合线性模型T值及P值的计算
001、
root@DESKTOP-1N42TVH:/home/test4# ls gwas_test.bed gwas_test.bim gwas_test.fam root@DESKTOP-1N42TVH:/home/test4# /home/software/gemma-0.98.5-linux-static-AMD64 -bfile gwas_test -gk -o kin 1> /dev/null **** INFO: Done. root@DESKTOP-1N42TVH:/home/test4# /home/software/gemma-0.98.5-linux-static-AMD64 -bfile gwas_test -k output/kin.cXX.txt -lmm -o result 1> /dev/null **** INFO: Done. root@DESKTOP-1N42TVH:/home/test4# ls gwas_test.bed gwas_test.bim gwas_test.fam output root@DESKTOP-1N42TVH:/home/test4# cd output/ root@DESKTOP-1N42TVH:/home/test4/output# ls kin.cXX.txt kin.log.txt result.assoc.txt result.log.txt root@DESKTOP-1N42TVH:/home/test4/output# head -n 3 result.assoc.txt chr rs ps n_miss allele1 allele0 af beta se logl_H1 l_remle p_wald 1 snp1 2802 0 G T 0.099 9.393582e-01 7.917546e+00 -3.077033e+03 2.285023e+01 9.056027e-01 1 snp2 2823 0 T C 0.063 5.449975e+00 9.371469e+00 -3.076887e+03 2.260596e+01 5.611132e-01 root@DESKTOP-1N42TVH:/home/test4/output# awk 'NR != 1 {print $(NF - 4)/$(NF - 3)}' result.assoc.txt > t.txt root@DESKTOP-1N42TVH:/home/test4/output# head -n 3 t.txt 0.118643 0.58155 0.985854 root@DESKTOP-1N42TVH:/home/test4/output# paste result.assoc.txt <(sed '1i t' t.txt ) > data.txt root@DESKTOP-1N42TVH:/home/test4/output# head data.txt chr rs ps n_miss allele1 allele0 af beta se logl_H1 l_remle p_wald t 1 snp1 2802 0 G T 0.099 9.393582e-01 7.917546e+00 -3.077033e+03 2.285023e+01 9.056027e-010.118643 1 snp2 2823 0 T C 0.063 5.449975e+00 9.371469e+00 -3.076887e+03 2.260596e+01 5.611132e-010.58155 1 snp3 4512 0 G A 0.067 9.014939e+00 9.144296e+00 -3.076572e+03 2.237338e+01 3.246472e-010.985854 1 snp4 16529 0 T C 0.055 -2.639543e+00 1.041215e+01 -3.076994e+03 2.293325e+01 7.999739e-01-0.253506 1 snp5 16578 0 G C 0.065 -2.507825e+00 9.000474e+00 -3.077046e+03 2.274653e+01 7.806337e-01-0.278633 1 snp6 16579 0 C A 0.275 1.898949e+00 5.196884e+00 -3.077006e+03 2.320171e+01 7.149551e-010.365401 1 snp7 16635 0 G C 0.088 4.344375e+00 8.115947e+00 -3.076915e+03 2.244365e+01 5.926710e-010.535289 1 snp8 20879 0 T G 0.130 3.405991e-01 6.905055e+00 -3.077041e+03 2.264882e+01 9.606777e-010.0493261 1 snp9 20908 0 C T 0.167 -1.752315e-01 6.246462e+00 -3.077057e+03 2.267830e+01 9.776303e-01-
002、R语言中验证
dir() dat <- read.table("data.txt", header = T) head(dat, 3) p_verify <- vector() for (i in 1:nrow(dat)) { p_verify[i] <- 2 * pt(-abs(dat[i,13]), df = 541 - 1) } dat$p_verify <- p_verify head(dat,5) cor(dat$p_wald, dat$p_verify)