今天度过了求知的一天,求知的快乐就是这么朴实无华且枯燥。 今天同事问了我一个问题,为什么plink计算的pca和GCTA计算得不一样?然后就引出的今天的查看说明文档,也证明了是介绍就怕认真二字。 我们的发现: 1,GCTA的说明文档中,有bug,公式没有写全:最后一个公式还要除以N。给出的2010 NG上有写,但是软件的说明文档中不完整。 2,GCTA计算PCA时,中间要构建G矩阵,G矩阵构建的方法有两种:- VanRaden #GS中GBLUP构建的G矩阵方法
两种方法计算PCA的代码: system("gcta --bfile aaa --autosome-num 29 --make-grm --make-grm-alg 1 --out kinship") system("gcta --grm kinship --pca 3 --out gcta")
system("gcta --bfile aaa --autosome-num 29 --make-grm --make-grm-alg 0 --out kinship0") system("gcta --grm kinship0 --pca 3 --out gcta0")
3,plink计算PCA时,用的是yang 的方法所以,如果如果plink的PCA和GCTA的VanRaden方法相遇时,结果就不一致了。 主要是因为计算PCA时G矩阵的方法不一致。 4,怎么证明呢?有一个包叫AGHmatrix 包,里面有个Gmatrix ,它构建矩阵时可以选择构建的方法: 结果证明了两者确实不一样。 5,自己构建G矩阵,手动计算PCA# 计算特征值和特征向量 re = eigen(Gmat)
# 计算解释百分比 por = re$values/sum(re$values)
# 整理格式 pca_re1 = re$vectors[,1:3] pca_re2 = data.frame(pca_re1,Ind = iid)
6,plink的--kinship 构建的G矩阵plink的--kinship 构建的G矩阵不是VanRaden的矩阵,而是yang 的矩阵,所以不能用于GBLUP的分析 7,pca用什么方法?推荐用Yang 的方法构建G矩阵,得到的PCA结果。也就是plink的--pca 的结果,同样也是gcta 默认的计算PCA的参数,--make-grm-alg 0 。 8,为什么要用GCTA计算PCA?因为GCTA给出了每个PCA的特征值,可以用于计算PCA解释的百分比。plink默认没有给出所有的(应该也可以指定PCA的个数,然后手动计算,待验证)。 也可以用plink的--kinship 计算yang 的G矩阵,然后手动计算PCA,这样就可以计算百分比了,计算代码: # 计算特征值和特征向量 re = eigen(Gmat)
# 计算解释百分比 por = re$values/sum(re$values)
特别致谢同事的一个问题,让我理清了很多东西,也让我知道了“世界上就怕认真二字”,也更让我坚信“说明文档是最香的”,多查资料,相互佐证,自己验证,然后豁然开朗。 工具人的自我修养,调包侠的核心技能!
基因型数据绘制PCA图和聚类分析图
GWAS分析新软件 | GMATs:解析复杂性状和复杂遗传机制的高效工具
|