各种软件/包构建G矩阵结果比较
各种软件/包构建G矩阵结果比较
邓飞
发布于 2021-11-29 14:25:04
发布于 2021-11-29 14:25:04
代码可运行
运行总次数：0
代码可运行
飞哥注：这篇是我同事苏惠写的，内容更全面，代码更完整，我的上一篇plink计算的PCA为什么和GCTA计算的不一样？是一个引子，而且这一篇给出了plink --pca 样本数时，
，可以计算每个PCA的百分比，和我的设想一样。另外，发现GCTA真是一个宝库，
它也可以计算遗传力、遗传相关、GBLUP等参数。下一步学习起来！
下面是正文
故事是从PCA开始的，要算前3个PC解释的方差比例，然鹅Plink并不给输出所有PC的eigen value，于是转而用GCTA计算PCA以及前3个PC的eigen value。GCTA计算PCA首先需要构建kinship矩阵，也就是G矩阵，然后使用kinship矩阵计算PCA。算完PCA发现GCTA算的PCA结果居然和Plink不一样，然后就很好想知道为啥不一样，然后就开始研究各种软件/包构建G矩阵基于的算法和结果的异同。
GCTA构建GRM（Genetic Relationship Matrix，遗传关系矩阵，G矩阵）有2种方法可选，--make-grm-alg 0 基于Yang et al. 2010提出的方法，--make-grm-alg 1 基于Van Raden 2008提出的方法，默认0。
# GRM 0 Yang sum{[(xij- 2pi)*(xik- 2pi)] / [2pi(1-pi)]}（注：这是GCTA官网给出的公式，但这个应该是打错了，公式没写全，后面还有一个*1/N，GCTA的paper里的公式是对的） # GRM 1 Van Raden sum[(xij- 2pi)(xik- 2pi)] / sum[2pi(1-pi)]
用Van Raden方法构建GRM然后计算得到的PCA和Plink结果不一致，用Yang的方法得到的PCA结果和Plink的一致。这个很合理，因为虽然在--pca那里没提，但--make-rel那里Plink官方文档里说了是使用的Yang的方法。所以，虽然Plink表面上是直接输入基因型数据就输出PCA结果，但中间应该也是先构建了G阵，并且构建G阵的方法是使用Yang的方法，然后再基于G阵计算了PCA，只不过这个过程Plink直接帮我做了。
# GCTA GRM 0
gcta --bfile plink9996loci --autosome-num 29 --make-grm --make-grm-alg 0 --out kinship0
gcta --grm kinship --pca 3 --out gcta0
head -n 5 gcta0.eigenvec
0 sample1 0.0650311 0.0628848 -0.0760136
0 sample2 0.0804442 -0.0591764 0.0106865
0 sample3 -0.0243089 -0.00532782 0.0797794
0 sample4 0.0904922 0.0129799 -0.0417386
0 sample5 0.0458302 -0.0731895 0.0243814
# GCTA GRM 1
gcta --bfile plink9996loci --autosome-num 29 --make-grm --make-grm-alg 1 --out kinship
gcta --grm kinship --pca 3 --out gcta
head -n 5 gcta.eigenvec
0 sample1 -0.0349996 0.132768 0.158889
0 sample2 -0.0876233 -0.0556155 0.0798619
0 sample3 0.055747 -0.104418 0.00153495
0 sample4 -0.0998956 0.0625678 0.107681
0 sample5 -0.0652175 -0.0814918 -0.0759696
# Plink
plink --bfile plink9996loci --chr-set 29 --pca 3 header
suhui@suhui-MS-7B23:pca_testGCTA$ head -n 5 plink.eigenvec
FID IID PC1 PC2 PC3
0 sample1 -0.0650311 -0.0628848 -0.0760136
0 sample2 -0.0804442 0.0591764 0.0106865
0 sample3 0.0243089 0.00532782 0.0797794
0 sample4 -0.0904922 -0.0129799 -0.0417386前面说的不严谨，其实用Plink也不是不能算解释的方差比例，就有略有一丢丢麻烦，只要指定输出所有PC的eigen value其实就可以了，得到的结果和GCTA的一样：
system("plink --bfile plink9996loci --chr-set 29 --pca 150 header --out plink_pc150")
val <- scan("plink_pc150.eigenval")
var1 = round(val[1]/sum(val),4)
var2 = round(val[2]/sum(val),4)
var3 = round(val[3]/sum(val),4)
> head(val)
[1] 14.51180 10.94780 9.39109 8.68224 7.26250 7.01008
> sum(val)
[1] 179.1556
> var1
[1] 0.081
> var2
[1] 0.0611
> var3
[1] 0.0524
val <- scan("gcta0.eigenval")
var1 = round(val[1]/sum(val),4)
var2 = round(val[2]/sum(val),4)
var3 = round(val[3]/sum(val),4)
> head(val)
[1] 14.51180 10.94780 9.39109 8.68224 7.26250 7.01008
> sum(val)
[1] 179.1556
> var1
[1] 0.081
> var2
[1] 0.0611
> var3
[1] 0.0524以上是PCA对比结果，结束了。接下来是G矩阵：
### Plink构建G阵结果 ###
system("plink --bfile plink9996loci --chr-set 29 --make-rel square")
gmat_plink=read.table(file="plink.rel")
id=read.table(file="plink.rel.id")
colnames(gmat_plink)=id$V2
rownames(gmat_plink)=id$V2
> gmat_plink[1:5,1:5]
sample1 sample2 sample3 sample4 sample5
sample1 0.8703640 0.0272451 -0.1577940 0.1667670 -0.1117380
sample2 0.0272451 0.4955400 0.0446270 0.0691964 0.1130490
sample3 -0.1577940 0.0446270 1.0501800 -0.0994707 0.0471838
sample4 0.1667670 0.0691964 -0.0994707 0.6174680 -0.0740188
sample5 -0.1117380 0.1130490 0.0471838 -0.0740188 1.0315400
### 查看GCTA构建G阵结果 ###
# R script to read the GRM binary file | 来自GCTA官网
ReadGRMBin=function(prefix, AllN=F, size=4){
sum_i=function(i){
return(sum(1:i))
}
BinFileName=paste(prefix,".grm.bin",sep="")
NFileName=paste(prefix,".grm.N.bin",sep="")
IDFileName=paste(prefix,".grm.id",sep="")
id = read.table(IDFileName)
n=dim(id)[1]
BinFile=file(BinFileName, "rb");
grm=readBin(BinFile, n=n*(n+1)/2, what=numeric(0), size=size)
NFile=file(NFileName, "rb");
if(AllN==T){
N=readBin(NFile, n=n*(n+1)/2, what=numeric(0), size=size)
}
else N=readBin(NFile, n=1, what=numeric(0), size=size)
i=sapply(1:n, sum_i)
return(list(diag=grm[i], off=grm[-i], id=id, N=N))
}
GCTA GRM 0
> gmat_gcta0=ReadGRMBin("kinship0")
> head(gmat_gcta0$off)
[1] 0.02724512 -0.15779422 0.04462700 0.16676699 0.06919638 -0.09947072
# 和Plink结果一致
GCTA GRM 1
> gmat_gcta1=ReadGRMBin("kinship")
> head(gmat_gcta1$off)
[1] 0.01600686 -0.32093960 0.08792991 0.29738194 0.11218226 -0.22020614
# 和plink结果不一样
### 基于两种G阵构建方法写代码手算 ###
system("plink --bfile plink9996loci --chr-set 29 --recodeA --out plink9996loci")
# 读取数据
m012 = fread("plink9996loci.raw")
# 保留FID，IID和基因型数据
g012 = m012[,-c(3:6)]
dim(g012)
fid = g012$FID
iid = g012$IID
# 整理格式
setDF(g012)
rownames(g012) = g012$IID
g012$IID = NULL
g012$FID = NULL
g012=as.matrix(g012)
> g012[1:10,1:10]
1_5802_G 1_5823_C 1_5952_T 1_6121_T 1_6357_C 1_6597_T 1_7182_G 1_7934_A 1_8061_A 1_8282_T
sample1 0 0 0 0 0 0 0 0 0 0
sample2 0 0 0 0 0 0 0 0 1 1
sample3 0 0 0 0 0 0 0 0 0 0
sample4 0 0 0 0 0 0 0 0 0 0
sample5 0 2 0 0 0 0 1 0 1 1
sample6 0 2 1 0 0 1 1 0 1 1
sample7 2 2 0 0 0 0 0 0 1 0
sample8 1 1 0 0 0 0 0 0 1 1
sample9 0 2 0 0 0 1 0 0 1 1
sample10 0 0 0 0 1 0 0 0 0 0
# 计算MAF
system("plink --bfile plink9996loci --chr-set 29 --freq --out maf")
maf=read.table(file="maf.frq",header = T)
# 以样本1和2的G阵值为例
# grm 0
> (sum(((g012[1,]-2*maf$MAF)*(g012[2,]-2*maf$MAF))/(2*maf$MAF*(1-maf$MAF))))/nrow(maf)
[1] 0.02724698 # 和GCTA GRM 0结果一致
# grm 1
> sum((g012[1,]-2*maf$MAF)*(g012[2,]-2*maf$MAF))/sum(2*maf$MAF*(1-maf$MAF))
[1] 0.01600781 # 和GCTA GRM 1结果一致
### 使用DMU的Gmatrix模块 ###
# par文件不加任何质控和scale：
$MINMAF
0
$FREQMETHOD
1
$SCALEMETHOD
1
$DIAG_ADD
0
$G_ADD
0
$DIAG_ONE
0
$AGSCALE
0
$PROP_A_to_G
0
$CAL_DET
0
$OUT_GMATRIX
1
$OUT_IGMATRIX
0
$OUT_AMATRIX
0
gmat_dum=read.table(file="test.gmat")
> head(gmat_dum)
V1 V2 V3
1 1 1 1.23386414
2 1 2 0.01600686
3 1 3 -0.32093966
4 1 4 0.29738199
5 1 5 -0.22977475
6 1 6 0.10775250
# 使用的Van Raden 2008方法，和GCTA GRM 1结果一致
### 使用sommer包 ###
system("plink --bfile plink9996loci --chr-set 29 --recodeA --out plink9996loci")
# 读取数据
m012 = fread("plink9996loci.raw")
# 保留FID，IID和基因型数据
g012 = m012[,-c(3:6)]
dim(g012)
fid = g012$FID
iid = g012$IID
library(sommer)
# 整理格式，计算G矩阵
setDF(g012)
rownames(g012) = g012$IID
g012$IID = NULL
g012$FID = NULL
g012=as.matrix(g012)
Gmat = A.mat(g012-1)
> Gmat[1:5,1:5]
sample1 sample2 sample3 sample4 sample5
sample1 1.23386414 0.01600686 -0.32093966 0.2973820 -0.2297748
sample2 0.01600686 0.79154839 0.08792993 0.1121823 0.2459466
sample3 -0.32093966 0.08792993 1.14676773 -0.2202062 0.1277878
sample4 0.29738199 0.11218228 -0.22020617 0.9891461 -0.1771543
sample5 -0.22977475 0.24594663 0.12778778 -0.1771543 1.2713619
# 使用的Van Raden方法，和DMU-Gmatrix以及GCTA GRM 1结果一致
### 使用AGHmatrix包 ###
library(AGHmatrix)
gmat_van=Gmatrix(g012,method = "VanRaden")
gmat_yang=Gmatrix(g012,method = "Yang")
> gmat_van[1:5,1:5]
sample1 sample2 sample3 sample4 sample5
sample1 1.23386414 0.01600686 -0.32093966 0.2973820 -0.2297748
sample2 0.01600686 0.79154839 0.08792993 0.1121823 0.2459466
sample3 -0.32093966 0.08792993 1.14676773 -0.2202062 0.1277878
sample4 0.29738199 0.11218228 -0.22020617 0.9891461 -0.1771543
sample5 -0.22977475 0.24594663 0.12778778 -0.1771543 1.2713619
# 使用的Van Raden方法，和DMU-Gmatrix、GCTA GRM 1、sommer::A.mat结果一致
> gmat_yang[1:5,1:5]
sample1 sample2 sample3 sample4 sample5
sample1 1.19944230 0.02724512 -0.15779421 0.16676698 -0.11173845
sample2 0.02724512 1.04046826 0.04462700 0.06919638 0.11304946
sample3 -0.15779421 0.04462700 1.10626955 -0.09947072 0.04718378
sample4 0.16676698 0.06919638 -0.09947072 1.15327070 -0.07401880
sample5 -0.11173845 0.11304946 0.04718378 -0.07401880 1.26881311
# 使用的Yang方法，和Plink、GCTA GRM 0结果一致所以2种方法结果差异还蛮大的，那这2种方法的相关系数高么？
# 整个矩阵
> cor(as.vector(gmat_van),as.vector(gmat_yang))
[1] 0.9099529
# 矩阵对角线
> cor(gmat_gcta0$diag,gmat_gcta1$diag)
[1] 0.8409531
# 矩阵非对角线
> cor(gmat_gcta0$off,gmat_gcta1$off)
[1] 0.8960708相关系数也不是很高。所以现在就很困惑了，这2种方法到底哪个构建的G阵是对的呢？似乎GWAS系的软件大都用Yang的方法，GS系的软件/包大都用Van Raden的方法。。。why？？两个作者都是大神，两篇文献引用都很高，所以到底哪个是对的？？？
有一个设想来验证，通过PCA结果，用一个聚类信息明确的群体当做真值，分别使用2种方法构建的矩阵做PCA，看看哪个方法得到的PCA更接近群体真实的聚类情况，不过不知道能不能行得通，猜测很有可能对于聚类明确的群体，使用哪个方法都能得到正确聚类结果，而对于聚类不明确或遗传背景不均一的群体，两种方法得到的结果会差异较大。待验证。。。
阅读原文，查看作者知乎主页。
本文参与 腾讯云自媒体同步曝光计划，分享自微信公众号。
原始发表：2021-11-26，如有侵权请联系 cloudcommunity@tencent.com 删除
评论
登录后参与评论
推荐阅读
腾讯云开发者
