||
作者:林元震(https://yzhlinscau2018.netlify.com)
相关矩阵的热图为何带概率或显著性
使得热图更有统计学意义;
ggplot2能否实现这种图形。
ggplot2中热图的基础绘图方法如下:
ggplot(data, aes(x=,y=, fill=))+ geom_tile()
其中data代表数据集,aes代表数据映射,x代表变量名,y代表变量名或样本号,fill参数为数值型变量,geom_tile代表几何对象为瓦片图。
ggplot2不能直接绘制热图,需要做数据变换,最后的数据格式保留为3列,分别是变量名、样本号(或变量名)和数值,再分别映射给x、y和fill参数。此外,不同变量直接的数值差异较大,所以通过scale()对所有变量的数值标准化。这样就可降低值差,给图例一个更合理的色域,也利于图形美观。
程序包'RSTAT2D'可在http://yzhlin-asreml.ys168.com/免费下载。
data(d8.1.6,package='RSTAT2D')
mydata <- d8.1.6[,-1]
tdd <- as.data.frame(scale((mydata)))
head(tdd)
## h dbh v str wd tw
## 1 -0.7771651 0.1508848 -0.611934185 0.7352031 -0.5350928 -0.57061982
## 2 -0.4637808 0.1508848 -0.346070351 0.7352031 -0.5350928 -0.57061982
## 3 0.4260068 0.4495434 -0.918700148 -1.1221521 -0.7830405 -0.99018209
## 4 0.8904872 0.2172534 0.001597739 -1.1221521 -0.5350928 0.22009369
## 5 -0.4637808 -0.2307345 0.226559445 -0.5030337 0.4129422 1.04308123
## 6 1.3269868 -0.3302874 0.247010509 -1.1221521 -0.6080186 -0.04616698
## larea dis lt
## 1 0.5675137 -0.8680278 1.1899014
## 2 0.5675137 -0.8680278 1.1899014
## 3 0.4739996 1.1160357 1.1899014
## 4 -0.9661175 1.1160357 -0.8141431
## 5 -0.8351978 -0.8680278 -0.8141431
## 6 -1.5646078 1.1160357 -0.8141431
tdd$Prov=d8.1.6$Prov
head(tdd)
## h dbh v str wd tw
## 1 -0.7771651 0.1508848 -0.611934185 0.7352031 -0.5350928 -0.57061982
## 2 -0.4637808 0.1508848 -0.346070351 0.7352031 -0.5350928 -0.57061982
## 3 0.4260068 0.4495434 -0.918700148 -1.1221521 -0.7830405 -0.99018209
## 4 0.8904872 0.2172534 0.001597739 -1.1221521 -0.5350928 0.22009369
## 5 -0.4637808 -0.2307345 0.226559445 -0.5030337 0.4129422 1.04308123
## 6 1.3269868 -0.3302874 0.247010509 -1.1221521 -0.6080186 -0.04616698
## larea dis lt Prov
## 1 0.5675137 -0.8680278 1.1899014 Prov001
## 2 0.5675137 -0.8680278 1.1899014 Prov002
## 3 0.4739996 1.1160357 1.1899014 Prov003
## 4 -0.9661175 1.1160357 -0.8141431 Prov004
## 5 -0.8351978 -0.8680278 -0.8141431 Prov005
## 6 -1.5646078 1.1160357 -0.8141431 Prov006
# 使用reshape2的melt()进行数据融合
mdd <- reshape2::melt(tdd)
## Using Prov as id variableshead(mdd)
## Prov variable value
## 1 Prov001 h -0.7771651
## 2 Prov002 h -0.4637808
## 3 Prov003 h 0.4260068
## 4 Prov004 h 0.8904872
## 5 Prov005 h -0.4637808
## 6 Prov006 h 1.3269868
默认的热图见Figure 1:
(pheat=ggplot(data = mdd, aes(x=variable, y=Prov, fill=value))+geom_tile())
由于采用默认的填充色,色差较小,热图效果不好,这时可设置渐变色来加以调整。
pheat+ scale_fill_gradientn(colours = terrain.colors(20))
ggplot2中相关图的绘图方法如下:
ggplot(data, aes(x=,y=, fill=))+ geom_tile()
其中data代表数据集,aes代表数据映射,x代表变量名,y代表变量名,fill参数为数值型变量,geom_tile代表几何对象为瓦片图。
mcor <- round(cor(mydata),2)
cormat <- melt(mcor)
(pcorr=ggplot(data = cormat, aes(x=Var1, y=Var2, fill=value)) +
geom_tile()+scale_fill_gradient2(low="darkred", high="darkgreen"))
从生成的相关图看出,相关值正负排列比较混乱,因此可以对相关矩阵进行重排序(自定义函数),然后再绘制相关图。
reorder_cormat <- function(cormat){
dd <- as.dist((1-cormat)/2)
hc <- hclust(dd)
cormat <-cormat[hc$order, hc$order]
}
mcor1=reorder_cormat(mcor)
cormat1 <- melt(mcor1)
pcorr %+% cormat1
ggplot2不能直接绘制热图与相关图的结合图,这时,可以自行编程heatmap1()来实现。具体代码见本文的最后。
heatmap1(mydata,type='data',Sig=FALSE,Nbreaks = 8)
## Loading required package: agricolae
注:Nbreaks控制图例的显示数量。
将p值和显著性添加到图形中。
heatmap1(mydata,type='data',Sig=TRUE,Nbreaks = 8)
现在对相关矩阵进行聚类排序。
heatmap1(mydata,type='data',Sig=TRUE,order=TRUE,Nbreaks = 6)
事实上,除了数据集和相关矩阵,heatmap1()函数可以绘制任意矩阵和代标签的另一个矩阵,下述举一个简单的例子,当然读者应当明白下述例子没有什么科学意义。
set.seed(2018)
# 生成11×11的随机值矩阵
labellabel<-matrix(runif(121),nrow=11)
label<-round(label,2)
heatmap1(label)
# 两个矩阵的粘贴
# matrix a and b should have the same dim
matrix.ps<-function(a,b){
N=nrow(a)#=nrow(b)
a=round(a,2) # matrix a must be numeric
# matrix b should be character
c=matrix(0,nrow=N,ncol=N)
for(i in 1:N){
for(j in 1:N) c[i,j]=paste(a[i,j],b[i,j],sep='\n ')
} return(c)
}
sig<-c('*','','.','**','***')
b<-matrix(sig,nrow=11,ncol=11)
# 生成标签矩阵
label2<-matrix.ps(label,b)
heatmap1(label,type='matrix',df.label=label2,
gtitle='corr',Nbreaks = 8)
当然,也可以对之前的示例mydata上添加随意的标签矩阵,而不是p值和显著性。
# 生成9×9的随意标签矩阵
label3=label2[1:9,1:9]
heatmap1(mydata,type='data',df.label=label3,Nbreaks = 6)
通过上述示例可知,ggplot2包确实功能很强大,我自编的函数heatmap1()可以用于绘制遗传评估方面生成的遗传相关矩阵及其显著性,下回找个时间演示一下。
heatmap1()的具体代码如下:
#' Create a Heatmap
#' @param df A data.frame or matrix containing only numeric data.
#' @param type Identify df to be 'matrix'(default) or 'data'.
#' @param df.label A matrix for heatmap labels.
#' @param gtitle guide or legend title.
#' @param Nbreaks A number controls legend breaks.
#' @param Sig if TRUE put pvalue and sig level for heatmap labels.
#' @param data.only Logical, if TRUE returns correlation and pvalue.
heatmap1 <-
function(df,type='matrix', df.label=NULL,gtitle=NULL,
Nbreaks=NULL,Sig= FALSE,order=FALSE,
data.only = FALSE) {
require(ggplot2) # ggplot2
require(reshape2) # melt data
require(agricolae) # count corr and p.value
reorder_cormat <- function(cormat){
dd <- as.dist((1-cormat)/2)
hc <- hclust(dd)
cormat <-cormat[hc$order, hc$order]
}
if(is.null(gtitle)) gtitle<-'legend'
if(type=='data') {
df.cor.p<-correlation(df)
df1<-df.cor.p$correlation
df1.p<-df.cor.p$pvalue
if(order==TRUE){
df1<-reorder_cormat(df1)
nm1<-rownames(df1)
df1.p<-df1.p[nm1,nm1]
}
}
if(type=='matrix') df1<-df
if(type=='corr.matrix') {
df1<-df
if(order==TRUE){
df1<-reorder_cormat(df)
nm1<-rownames(df1)
if(!is.null(df.label)) df.label<-df.label[nm1,nm1]
}
}
test <- melt(df1)
if(type=='data'){
test$p.value<-round(melt(df1.p)$value,2)
ra<-abs(test$p.value)
NN<-nrow(test)
prefix<-rep('',NN)
for(i in 1:NN){
if(ra[i]<=0.1) prefix[i] <- '.'
if(ra[i]<=0.05) prefix[i] <- '*'
if(ra[i]<=0.01) prefix[i] <- '**'
if(ra[i]<=0.001) prefix[i] <- '***'
}
ra1<-paste(test$p.value,prefix,sep='\n')
test$rp<-ra1
}
if(is.null(df.label)) {
if(Sig==FALSE) test.label=round(test$value,2)
else test.label<-test$rp
}
if(!is.null(df.label)) {
test1<-melt(df.label)
if(is.numeric(test1$value)) test.label<-round(test1$value,2)
else test.label<-test1$value
}
p1<-ggplot(test,aes(x=Var1,y=Var2,fill=value,label=test.label))+
geom_tile() +
geom_text() +
labs(x="",y="",fill=gtitle)
if(is.null(Nbreaks))
p1<-p1+scale_fill_distiller(palette="Spectral",
trans = "reverse",
guide = "legend")
if(!is.null(Nbreaks)) {
bv<-unique(test$value)
bv<-bv[order(bv)]
Nbv<-length(bv)
nn<-seq(2,Nbv-1,by=Nbreaks)
bv1<-bv[c(1,nn,Nbv)]
bv2<-rev(bv1)
p1<-p1+scale_fill_distiller(palette = "Spectral",
trans = "reverse",
breaks = bv1,
guide = "legend")
}
if(data.only) {
return(test)
}
print(p1)
}
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-9-27 07:09
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社