您现在的位置是：首页 > 数据库

当前栏目

各种软件/包构建G矩阵结果比较

计算

2023-03-20 14:56:00 时间

飞哥注：这篇是我同事苏惠写的，内容更全面，代码更完整，我的上一篇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更接近群体真实的聚类情况，不过不知道能不能行得通，猜测很有可能对于聚类明确的群体，使用哪个方法都能得到正确聚类结果，而对于聚类不明确或遗传背景不均一的群体，两种方法得到的结果会差异较大。待验证。。。

阅读原文，查看作者知乎主页。

猜你喜欢

Python中的函数与方法以及Bound Method和Unbound Method
从本体论开始说起——运营商关系图谱的构建及应用
一篇运维老司机的大数据平台监控宝典（2）-联通大数据集群平台监控体系详解
一篇运维老司机的大数据平台监控宝典（1）-联通大数据集群平台监控体系进程详解
Flask中的请求上下文和应用上下文
深入探讨Java中的异常与错误处理
研究学习Kotlin的一些方法
如何成为一名数据科学家？
金融服务领域的大数据：即时分析
影响大数据、机器学习和人工智能未来发展的8个因素
从未见过的堂兄杀了人，你的DNA是关键证据
一文贯通python文件读取
数据显示Java热度持续下落，日子屈指可数？
从0开始构建一个属于你自己的PHP框架
如何将Hadoop集成到工作流程中？这6个优秀实践必看
2017年5月编程语言排行榜：Java与C语言优势正开始缩小
SEO公司使用大数据优化其模型的5种方法
Java多线程之内置锁与显示锁
关于Web Workers你需要了解的七件事
20个安全可靠的免费数据源，各领域数据任你挑

zl程序教程

当前栏目

各种软件/包构建G矩阵结果比较

相关文章