ICode9

精准搜索请尝试: 精确搜索
首页 > 其他分享> 文章详细

GWAS分析 一般线性模型GLM + 协变量中meta、se、T、p值的计算

2022-07-27 11:34:17  阅读:147  来源: 互联网

标签:GWAS GLM PC1 gwas meta test home root plink


 

001、plink

root@PC1:/home/test# ls
gwas_test.map  gwas_test.ped
root@PC1:/home/test# plink --file gwas_test --pca 3 1> /dev/null
root@PC1:/home/test# ls
gwas_test.map  gwas_test.ped  plink.eigenval  plink.eigenvec  plink.log
root@PC1:/home/test# head -n 3 plink.eigenvec
fam1 id1 -0.0102664 0.0442708 -0.0658835
fam2 id2 0.0421124 -0.0313769 -0.0705913
fam3 id3 -0.0467758 -0.00993805 -0.000487215
root@PC1:/home/test# plink --file gwas_test --covar plink.eigenvec --linear hide-covar 1> /dev/null
root@PC1:/home/test# ls
gwas_test.map  gwas_test.ped  plink.assoc.linear  plink.eigenval  plink.eigenvec  plink.log
root@PC1:/home/test# head -n 5 plink.assoc.linear
 CHR        SNP         BP   A1       TEST    NMISS       BETA         STAT            P
   1       snp1       2802    G        ADD      541     0.9887       0.1258       0.8999
   1       snp2       2823    T        ADD      541      6.096       0.6524       0.5144
   1       snp3       4512    G        ADD      541      9.957        1.091       0.2759
   1       snp4      16529    T        ADD      541      5.631       0.5394       0.5899
root@PC1:/home/test# head -n 5 plink.assoc.linear | awk '{if(NR == 1) {print "SE"} else {print $(NF -2 )/$(NF - 1)}}'
SE
7.8593
9.34396
9.12649
10.4394

 

 

002、R

root@PC1:/home/test# ls
gwas_test.map  gwas_test.ped
root@PC1:/home/test# plink --file gwas_test --recode A 1> /dev/null
root@PC1:/home/test# ls
gwas_test.map  gwas_test.ped  plink.log  plink.raw
root@PC1:/home/test# plink --file gwas_test --pca 3 header 1> /dev/null
root@PC1:/home/test# ls
gwas_test.map  gwas_test.ped  plink.eigenval  plink.eigenvec  plink.log  plink.raw
root@PC1:/home/test# awk '{print $2, $6}' plink.raw | paste -d " " - <(cut -d " " -f 3- plink.eigenvec ) | paste -d " " - <(cut -d " " -f 7- plink.raw ) > dat.txt
root@PC1:/home/test# head -n 3 dat.txt | cut -d " " -f 1-10 | column -t
IID  PHENOTYPE  PC1         PC2         PC3         snp1_G  snp2_T  snp3_G  snp4_T  snp5_G
id1  580        -0.0102664  0.0442708   -0.0658835  0       1       1       0       1
id2  690        0.0421124   -0.0313769  -0.0705913  0       0       0       0       0

 

library(data.table)
dat <- fread("dat.txt", data.table = F)
dat[1:3,1:8]

result <- data.frame()
for (i in 6:12) {
  reg = lm(PHENOTYPE ~  1 + PC1 + PC2 + PC3 + dat[,i], data=dat)
  result<- rbind(result, summary(reg)$coefficients[5,])
}
names(result) <- c("beta", "se", "t", "p")
result

 

标签:GWAS,GLM,PC1,gwas,meta,test,home,root,plink
来源: https://www.cnblogs.com/liujiaxin2018/p/16524215.html

本站声明: 1. iCode9 技术分享网(下文简称本站)提供的所有内容,仅供技术学习、探讨和分享;
2. 关于本站的所有留言、评论、转载及引用,纯属内容发起人的个人观点,与本站观点和立场无关;
3. 关于本站的所有言论和文字,纯属内容发起人的个人观点,与本站观点和立场无关;
4. 本站文章均是网友提供,不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属;如您发现该文章侵犯了您的权益,可联系我们第一时间进行删除;
5. 本站为非盈利性的个人网站,所有内容不会用来进行牟利,也不会利用任何形式的广告来间接获益,纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

专注分享技术,共同学习,共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有