环境：R 4.2.2，Rstudio 2022.07.2，win11
1、gcta计算G矩阵
gcta64 --bfile 20221115-393id --make-grm --out 20221115-393id
2、安装R包ggtree、amap
install.packages('BiocManager')
install.packages('amap')
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("ggtree")
需要先安装BiocManager，再安装amap和ggtree
3、读取G矩阵
#R读取GRM文件G矩阵
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提供代码，只能读取G矩阵对角线
aa = ReadGRMBin(prefix="F:/gcta/YY-4637ID/20221115-393id")
#此处路径需要填写路径+上面--out 后面的文件名
G_mat = matrix(0,length(aa$diag),length(aa$diag))
diag(G_mat) = aa$diag
lowerTriangle(G_mat,byrow = T) = aa$off
G_mat = G_mat+t(G_mat)-diag(diag(G_mat))
rownames(G_mat) = colnames(G_mat) = aa$id$V2
G_mat[1:10,1:10]
n = sample(1:393,100)
Gmat1 = G_mat[n,n]
library(ggtree)
library(amap)
clu <- hclusterpar(Gmat1)
4、画图
ggtree(clu,layout="roundrect")+geom_tiplab(size=3)+xlim(NA,2)
img_v2_65ed523a-46a4-4479-9880-304320a58d9g.jpg
ggtree(clu,layout="circular")+geom_tiplab2(offset=0.1)+xlim(NA,3)
img_v2_3b02e601-fa47-44a9-bdca-ca6ca85a591g.jpg
