GWAS logistic回归分析

 

001、plink

root@PC1:/home/test# ls
gwas_case_cont.map  gwas_case_cont.ped
root@PC1:/home/test# plink --file gwas_case_cont --logistic beta 1> /dev/null       ## plink 逻辑回归
root@PC1:/home/test# ls
gwas_case_cont.map  gwas_case_cont.ped  plink.assoc.logistic  plink.log
root@PC1:/home/test# head plink.assoc.logistic
 CHR        SNP         BP   A1       TEST    NMISS       BETA         STAT            P
   1       snp1       3046    A        ADD      288    0.06317       0.2228       0.8237
   1       snp2       3092    T        ADD      288    0.08104       0.2895       0.7722
   1       snp3       3174    T        ADD      288   -0.08676      -0.2947       0.7682
   1       snp4      32399    T        ADD      288   -0.03895      -0.1416       0.8874
   1       snp5      32402    G        ADD      288   -0.05129      -0.2792       0.7801
   1       snp6      32406    G        ADD      288   -0.07479       -0.411       0.6811
   1       snp7      32443    C        ADD      288    -0.1864      -0.8866       0.3753
   1       snp8      32548    T        ADD      288   -0.05275      -0.2156       0.8293
   1       snp9      45044    A        ADD      288     0.1644       0.6751       0.4996

 

002、R语言实现

root@PC1:/home/test# ls
gwas_case_cont.map  gwas_case_cont.ped
root@PC1:/home/test# plink --file gwas_case_cont --recode A 1> /dev/null     ## 利用plink将基因型数据转换为数值型
root@PC1:/home/test# ls
gwas_case_cont.map  gwas_case_cont.ped  plink.log  plink.raw
root@PC1:/home/test# head plink.raw | cut -d " " -f 1-10
FID IID PAT MAT SEX PHENOTYPE snp1_A snp2_T snp3_T snp4_T
A1 A1 0 0 1 1 1 1 0 0
A2 A2 0 0 1 1 0 0 0 0
A3 A3 0 0 1 1 0 0 0 0
A4 A4 0 0 1 1 0 0 0 0
A5 A5 0 0 1 1 0 0 0 0
A6 A6 0 0 1 1 0 0 0 0
A7 A7 0 0 1 1 0 0 0 0
A8 A8 0 0 1 1 0 0 0 0
A9 A9 0 0 1 1 0 0 0 0

 

dir()
library(data.table)
dat <- fread("plink.raw", data.table = F)
dat[,6] <- dat[,6] - 1              ## 表型由 1、2改为0和1
dat <- dat[,-c(1,3:5)]

result <- data.frame()
for (i in 3:10) {
  logis <- glm(dat[,2]~dat[,i], family = "binomial", data = dat)       ## 逻辑回归
  result <- rbind(result, c(exp(logis$coefficients[2]),summary(logis)$coefficients[2,]))   ## 同时计算OR值
}
names(result) <- c("OR",names(summary(logis)$coefficients[2,]))
result

 

参考：

001、https://mp.weixin.qq.com/s?__biz=MzIwODA1MzI4Mg%3D%3D&chksm=809f9dc7b7e814d1c5928e46974e94ce5a2c8e35f13163065bb46df37e7a6e9f1d52bfdb2908&idx=1&mid=2456012938&scene=21&sn=783a6589eb200c065c9c1993649050ae#wechat_redirect

002、https://cloud.tencent.com/developer/article/1708838