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 --pca 3 1> /dev/null    ## 生成PCA协变量
root@PC1:/home/test# ls
gwas_case_cont.map  gwas_case_cont.ped  plink.eigenval  plink.eigenvec  plink.log
root@PC1:/home/test# head -n 3 plink.eigenvec | column -t
A1  A1  0.0161698   0.108569    -0.0182556
A2  A2  -0.0582408  0.00512212  0.0141237
A3  A3  -0.0363562  -0.0640895  0.0172107
root@PC1:/home/test# plink --file gwas_case_cont --covar plink.eigenvec --logistic hide-covar 1> /dev.null    ## 逻辑回归， 添加协变量
root@PC1:/home/test# ls
gwas_case_cont.map  gwas_case_cont.ped  plink.assoc.logistic  plink.eigenval  plink.eigenvec  plink.log
root@PC1:/home/test# head -n 3 plink.assoc.logistic       ## 查看结果
 CHR        SNP         BP   A1       TEST    NMISS         OR         STAT            P
   1       snp1       3046    A        ADD      288      1.069        0.237       0.8126
   1       snp2       3092    T        ADD      288      1.088       0.3019       0.7627

 

002、R

root@PC1:/home/test# ls
gwas_case_cont.map  gwas_case_cont.ped
root@PC1:/home/test# plink --file gwas_case_cont --pca 3 1> /dev/null
root@PC1:/home/test# ls
gwas_case_cont.map  gwas_case_cont.ped  plink.eigenval  plink.eigenvec  plink.log
root@PC1:/home/test# head -n 3 plink.eigenvec
A1 A1 0.0161698 0.108569 -0.0182556
A2 A2 -0.0582408 0.00512212 0.0141237
A3 A3 -0.0363562 -0.0640895 0.0172107
root@PC1:/home/test# plink --file gwas_case_cont --recode A 1> /dev/null
root@PC1:/home/test# ls
gwas_case_cont.map  gwas_case_cont.ped  plink.eigenval  plink.eigenvec  plink.log  plink.raw
root@PC1:/home/test# head -n 5 plink.raw | cut -d " " -f 1-8
FID IID PAT MAT SEX PHENOTYPE snp1_A snp2_T
A1 A1 0 0 1 1 1 1
A2 A2 0 0 1 1 0 0
A3 A3 0 0 1 1 0 0
A4 A4 0 0 1 1 0 0

 

dir()
pca <- read.table("plink.eigenvec", header = F)
dat <- fread("plink.raw", data.table = F)
dat <- cbind(dat[,c(2,6)], pca[,3:5], dat[,-c(1:6)])
names(dat)[3:5] <- c("pc1", "pc2", "pc3")
dat[,2] <- dat[,2] - 1
dat[1:3, 1:8]

result <- data.frame()
for (i in 6:10) {
  logis <- glm(dat[,2]~dat[,3] + dat[,4] + dat[,5] + dat[,i], family = "binomial",
               data = dat)                                                            ## 逻辑回归，增加协变量， 同时计算OR，OR = exp(beta); se = beta/stat
  result <- rbind(result, c(exp(logis$coefficients[5]),summary(logis)$coefficients[5,]))
}

names(result) <- c("OR", names(summary(logis)$coefficients[5,]))
result

 

 

 

plink验证：