GWAS分析中一般线性模型GML + 协变量中snp的回归系数r2的计算

 

001、plink + 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

 

读入R，计算r2：

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

result <- vector()
for (i in 6:12) {
  reg = lm(PHENOTYPE ~  1 + PC1 + PC2 + PC3 + dat[,i], data=dat)
  result[i - 5] <- summary(reg)$r.squared
}
result