R绘制G矩阵聚类图(环形、横向、纵向)

环境: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
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。

推荐阅读更多精彩内容