利用基因型数据构建G矩阵

 

001、plink + sommer实现

root@PC1:/home/test# ls
outcome.map  outcome.ped
root@PC1:/home/test# plink --file outcome --recode A 1> /dev/null   ## 将A\T\C\G转换为0\1\2形式
root@PC1:/home/test# ls
outcome.map  outcome.ped  plink.log  plink.nosex  plink.raw

 

library(sommer)
library(data.table)
dir()

dat <- fread("plink.raw")
g = as.matrix(dat[,!c(1:6)])
rownames(g) = dat$IID
gmat = A.mat(g-1)

gmat_cor <- gmat

for (i in 1:nrow(gmat_cor)) {
  for (j in 1:ncol(gmat_cor)) {
    gmat_cor[i,j] <- gmat[i,j]/sqrt(gmat[i,i] * gmat[j,j])
  }
}
dim(gmat_cor)
gmat_cor[1:5,1:5]
write.table(gmat_cor, "sommer_g.txt",row.names = T, col.names = T, sep = "\t", quote = F)

 

 

002、gemma实现

root@PC1:/home/test# ls
outcome.map  outcome.ped
root@PC1:/home/test# plink --file outcome --make-bed 1> /dev/null   ## 生成二进制文件
root@PC1:/home/test# ls
outcome.map  outcome.ped  plink.bed  plink.bim  plink.fam  plink.log
root@PC1:/home/test# /home/software/gemma-0.98.5-linux-static-AMD64 -bfile plink -gk 2 -o gemma -outdir ./  ## 构建g矩阵
GEMMA 0.98.5 (2021-08-25) by Xiang Zhou, Pjotr Prins and team (C) 2012-2021
Reading Files ...
## number of total individuals = 60
## number of analyzed individuals = 60
## number of covariates = 1
## number of phenotypes = 1
## number of total SNPs/var        =    42629
## number of analyzed SNPs         =    42629
Calculating Relatedness Matrix ...
================================================== 100%
**** INFO: Done.
root@PC1:/home/test# ls
gemma.log.txt  gemma.sXX.txt  outcome.map  outcome.ped  plink.bed  plink.bim  plink.fam  plink.log
root@PC1:/home/test# head -n 5 gemma.sXX.txt | cut -f 1-5
0.9890465761    0.1073153377    0.1182919093    0.1263857559    0.1589786895
0.1073153377    0.9722035427    0.1886718834    0.2216020996    0.1799985087
0.1182919093    0.1886718834    0.9714687294    0.1556666191    0.1414708648
0.1263857559    0.2216020996    0.1556666191    0.9272709102    0.141052078
0.1589786895    0.1799985087    0.1414708648    0.141052078     0.9287239597

 

dir()
sommer <- fread("sommer_g.txt")
sommer[1:5, 1:5]
gemma <- fread("gemma.sXX.txt")
gemma[1:5, 1:5]
cor(sommer$DOR1, gemma$V1)   ## 计算相关系数

 

 

003、plink实现

root@PC1:/home/test# ls
outcome.map  outcome.ped
root@PC1:/home/test# plink --file outcome --make-rel square 1> /dev/null    ## plink生成G矩阵
root@PC1:/home/test# ls
outcome.map  outcome.ped  plink.log  plink.rel  plink.rel.id
root@PC1:/home/test# head -n 6 plink.rel | cut -f 1-6
1.04625 0.122221        0.132869        0.142839        0.178275        0.131377
0.122221        1.03097 0.209148        0.246329        0.201715        0.235886
0.132869        0.209148        1.02611 0.174476        0.160126        0.206312
0.142839        0.246329        0.174476        0.979224        0.16131 0.2332
0.178275        0.201715        0.160126        0.16131 0.98128 0.135685
0.131377        0.235886        0.206312        0.2332  0.135685        0.984521

dir()
sommer <- fread("sommer_g.txt")
sommer[1:3, 1:3]
gemma <- fread("gemma.sXX.txt")
gemma[1:3, 1:3]
plink <- fread("plink.rel")
plink[1:3, 1:3]
cor(sommer$DOR1, plink$V1)
cor(gemma$V1, plink$V1)

 

004、gcta实现

root@PC1:/home/test# ls
outcome.map  outcome.ped
root@PC1:/home/test# plink --file outcome --make-bed 1> /dev/null      ## 转二进制
root@PC1:/home/test# ls
outcome.map  outcome.ped  plink.bed  plink.bim  plink.fam  plink.log                                                        ## 构建矩阵
root@PC1:/home/test# /home/software/gcta_v1.94.0Beta_linux_kernel_3_x86_64/gcta_v1.94.0Beta_linux_kernel_3_x86_64_static --bfile plink --make-grm-gz --make-grm-alg 1 --out gcta 1> /dev/null
root@PC1:/home/test# ls
gcta.grm.gz  gcta.grm.id  gcta.log  outcome.map  outcome.ped  plink.bed  plink.bim  plink.fam  plink.log
root@PC1:/home/test# gunzip gcta.grm.
gzip: gcta.grm..gz: No such file or directory
root@PC1:/home/test# gunzip gcta.grm.gz
root@PC1:/home/test# ls
gcta.grm  gcta.grm.id  gcta.log  outcome.map  outcome.ped  plink.bed  plink.bim  plink.fam  plink.log
root@PC1:/home/test# head -n 3 gcta.grm    ## 结果文件
1       1       1.593672e+04    1.068706e+00
2       1       1.593672e+04    1.294735e-01
2       2       1.593672e+04    1.078836e+00
root@PC1:/home/test# cut -f 1,2,4 gcta.grm > gcta.grm2   ## 修改格式
root@PC1:/home/test# ls
gcta.grm  gcta.grm2  gcta.grm.id  gcta.log  outcome.map  outcome.ped  plink.bed  plink.bim  plink.fam  plink.log

 

dir()
sommer <- fread("sommer_g.txt")
gemma <- fread("gemma.sXX.txt")
plink <- fread("plink.rel")

dat <- read.table("gcta.grm2")
head(dat,3)
gcta <- matrix(0, max(dat[,1]), max(dat[,2]))

gcta[as.matrix(dat[,1:2])] <- dat[,3]       ## 三元组转换为矩阵形式
gcta[1:4,1:4]
gcta <- gcta + t(gcta) - diag(gcta)
gcta[1:4,1:4]
cor(sommer$DOR1, gcta[,1])
cor(gemma$V1, gcta[,1])
cor(plink$V1, gcta[,1])