上一篇文章推送的是怎样调整corrplot热图的可视化参数,以修改字符和图例位置,数据可视化形式和字符小大和颜色等这篇是一个补充部分,记录怎样修改参数以变量排序方式和突出部分数据。本流程还是使用R统计绘图-环境因子相关性热图中的不同土壤环境因子数据进行相关性分析、绘制热图并进行图细节更改。流程开始按下图整理环境因子数据,行为样品名称,列为环境因子名称和分组信息,共有11个环境变量,3个分组信息。
图1|环境因子及分组信息表,env.csv。
1 设置工作路径并调用R包
#设置工作路径
#knitr::opts_knit$set(root.dir="D:\\EnvStat")#使用Rmarkdown进行程序运行
Sys.setlocale('LC_ALL','C') # Rmarkdown全局设置
#调用R包
library(Rmisc)
library(corrplot)
library(ggcorrplot)
library(RColorBrewer)
library(grDevices)
library(vegan)
2数据读入
options(stringsAsFactors=F)#R中环境变量设置,防止字符型变量转换为因子
#读入环境因子数据表
ENV=read.csv("env.csv",header = T,row.names = 1,sep = ",",comment.char = "",stringsAsFactors = F,colClasses = c(rep("character",4),rep("numeric",11)))head(ENV)#查看数据前几行#dim(ENV)#查看数据行、列数#str(ENV)#查看数据表每列的数据形式#解决数据中存在ties的警告,选用"pearson"fan方法计算相关性系数,不用运行此命令ENV[4:14]=rank(ENV[4:14],ties.method="random")#使用"kendall"或"spearman"方法计算相关性系数,可能会遇到此报错
3 绘制相关性热图
3.1 计算相关性系数、置信区间和显著性检验
# 3.1.1 环境因子相关性分析,根据自己的数据选择相关性计算方法:"pearson", "kendall", or "spearman"#env.cor <- round(cor(ENV[4:14], method = "spearman"),3) # round(),对输出结果取小数点前三位,计算spearman相关性系数#e?cor.mtestnv.cor <- round(cor(ENV[4:14], method = "kendall"),3) # round(),对输出结果取小数点前三位,计算kendall相关性系数env.cor <- round(cor(ENV[4:14], method = "pearson"),3) # round(),对输出结果取小数点前三位,计算pearson相关性系数,要求数据满足正态分布,正态分布检测方法见之前的推文。env.cor
# 3.1.2进行显著性检验
#env.p <-round(cor_pmat(ENV[4:14],method = "spearman"),3) # spearman系数
#env.p <-round(cor_pmat(ENV[4:14],method = "kendall"),3) # kendall系数
env.p <-round(cor_pmat(ENV[4:14],method = "pearson"),3) # pearson系数
env.p
# 3.1.3 使用corrplot的cor.mtest可以同时计算显著性检验和置信区间。返回是是一个列表
#cor.mtest(ENV[4:14], conf.level = 0.95)
env.uppCI = round(cor.mtest(ENV[4:14], conf.level = 0.95)$uppCI,3)
env.uppCI # 置信区间上界
env.lowCI = round(cor.mtest(ENV[4:14], conf.level = 0.95)$lowCI,3)
env.lowCI # 置信区间下界
图2|pearson相关性系数,env.cor 。
图3|pearson相关性系数的p值,env.p 。
图4|pearson相关性系数置信区间上界,env.uppCI 。
3.2 相关性系数热图绘制
corrplot绘制热图,变量顺序默认跟输入顺序一致,但这样的热图看起来不一定很直观。所以参数order可用于对相关性数据矩阵进行重排序。如果想要突出图中某个区域,可以使用addrect参数给指定区域添加矩形框。非相关性数据矩阵、NA和数学表达式也可以使用corrplot可视化。
3.2.1 更改变量可视化顺序
order参数可选:'AOE', 'FPC', 'hclust'和'alphabet'。
#### 3.2.1 更改变量顺序
?corrplot.mixed # 查看函数帮助信息
m=par(no.readonly = TRUE) # 保存默认图形设置参数
pdf("order.pdf",width = 8,height = 12,family = "Times") # 保存图片到本地,并设置字体为Time New Romanpar(mfrow=c(3,2))# 形成3行,2列的图形矩阵
# a图 使用corrplot.mixed()绘制upper和lower以不同形式展示相关性系数热图,默认参数变量排序。
corrplot.mixed(env.cor,env.p,
upper = "number",lower = "circle",
tl.pos = "lt",tl.col="black", diag="n")
# b图 order=“AOE”,AOE是 angular order of the eigenvectors。
corrplot.mixed(env.cor, env.p,
upper = "number",lower = "circle",
tl.pos = "lt",tl.col="black", diag="n",
order="AOE")
# c图 order="FPC",FPC表示first principal component order。
corrplot.mixed(env.cor, env.p,
upper = "number",lower = "circle",
tl.pos = "lt",tl.col="black", diag="n",
order="FPC")
# d图 order="hclust",hclust是hierarchical clustering order,设置层次聚类的同时,还需要选择聚类方法。## 设置'hclust.method参数,可选:'ward', 'single', 'complete', 'average', 'mcquitty', 'median'或'centroid'。corrplot.mixed(env.cor, env.p,upper = "number",lower = "circle",tl.pos = "lt",tl.col="black", diag="n",order="hclust",hclust="median")
# e图 order="hclust",hclust是hierarchical clustering order,设置层次聚类的同时,还需要选择聚类方法。## 设置'hclust.method参数,可选:'ward', 'single', 'complete', 'average', 'mcquitty', 'median'或'centroid'。corrplot.mixed(env.cor, env.p,upper = "number",lower = "circle",tl.pos = "lt",tl.col="black", diag="n",order="hclust",hclust="complete")
# f图 order="alphabet" alphabetical order表示按照字母顺序。
corrplot.mixed(env.cor, env.p,
upper = "number",lower = "circle",
tl.pos = "lt",tl.col="black", diag="n",
order="alphabet")
dev.off()
par(m)# 恢复默认图形设置参数
注:使用add=TRUE对两个相关性热图进行叠加,upper和lower需要设置相同的排序方式。R包seriation可以通过几种seriation/sequencing/ordination方法对矩阵、相异矩阵和树状图进行重新排序。感兴趣的同学可以查看seriation帮助文件进行学习。
图5|导入外部软件计算的相关性系数和p值绘图可能会报错。使用as.matrix()可以解决。
图6|corrplot修改order参数调整变量顺序。a: 默认参数;b:order=“AOE”;c:order="FPC";d:order="hclust",hclust="median";e:order="hclust",hclust="complete";f:order="alphabet"。
3.2.2 矩形框突出显示特定数据
corrplot()的addrect参数可以根据层次聚类发结果添加矩形框,以达到突出显示的效果。rect*参数可以修改矩形框的颜色和线宽度。需要注意的是,只有order="hclust"时此参数才起效。也可通过外部函数自定义矩形位置。corrRect.hclust()也可以根据层次聚类结果绘制矩形框,需要与corrplot()配合使用,可以可视化p值。corrRect()则可以自定义矩形框位置,需要corrplot()输出的列表数据和矩形位置信息作为输入,也需要与corrplot()配合使用进行可视化。
#### 3.2.2 矩形框突出显示特定数据
m=par(no.readonly = TRUE) # # 保存默认图形设置参
pdf("addrect.pdf",width = 8,height = 12,family = "Times") # 保存图片到本地,并设置字体为Time New Roman。par(mfrow=c(3,2))# 形成3行,2列的图形矩阵
# a图 addrect参数,此参数只能在type="full"以及没有p值矩阵的时候,才能正常运行。
corrplot(env.cor,
tl.col="black",
order = "hclust",addrect =2,
rect.col = "black",rect.lwd = 2) # 由图所示,pH与其它环境因子分开。
# b图 corrRect.hclust()根据层次聚类的结果获取每个cluster的成员,需要与corrplot()配合使用,可以可视化p值。使用1-corr作为距离进行层次聚类。?corrRect.hclust # 查看函数参数帮助信息corrplot(corr = env.cor,p.mat = env.p,order = "hclust",hclust.method = "ward.D",insig = "label_sig",sig.level = c(.05,0.01),tl.col = "black",pch.cex = 0.8)corrRect.hclust(env.cor, k = 3,method = 'ward.D') # method参数默认与corrplot设置的一致,k=3相当于addrect=3。帮助文件写可不在corrplot设置order参数,直接使用corrRect()中的method参数设置。但是我在实际操作中发现,必须在corrplot()中设置order和hclust.method参数,否则绘图行列变量名会与输入的一致,但是矩形框却是按照自定义的聚类方法形成的,导致绘图错误。所以建议直接在corrplot中设置聚类方法#corrRect()则可以自定义矩形框位置,需要corrplot()输出的列表数据和矩形位置信息(index,name或namesMat)作为输入。?corrRect # 查看函数参数帮助信息## c图,index参数使用序号数值指定矩形框起始终止位置,适用于行列名相同的矩阵。此处绘制与图b聚类结果相同的矩形框。p <- corrplot(corr = env.cor,p.mat = env.p,order = "hclust",hclust.method = "ward.D",insig = "label_sig",sig.level = c(.05,0.01),tl.col = "black",pch.cex = 0.8)corrRect(p,index = c(1,2, 8, 11),col = "black",lwd = 2) # index数值依次表示为第一个矩形的起始位置为变量1,第二个矩形的起始位置为变量2,第三个矩形的起始位置为变量8,第三个矩形的终止位置为变量11。
## d图 name参数使用变量名称指定矩形框起始终止位置,适用于行列名相同的矩阵。此处绘制与图b聚类结果相同的矩形框。name <- c("pH","AHN","TP","Ammonia")p <- corrplot(corr = env.cor,p.mat = env.p,order = "hclust",hclust.method = "ward.D",insig = "label_sig",sig.level = c(.05,0.01),tl.col = "black",pch.cex = 0.8)corrRect(p,name = name,col = "blue",lwd = 3) # name的参数依次表示为第一个矩形的起始位置为变量TN,第二个矩形的起始位置为变量TP,第三个矩形的起始位置为变量AK,第三个矩形的终止位置为变量OC。
## e图 namesMat参数需要4个字符向量或者四列字符矩阵,用于指定绘制矩形框的左、下、右和上侧边界。可用于行列变量名不一致的情况。nameMat <- rbind(c("pH","pH","pH","pH"),c("AHN","OC","OC","AHN"),c("TP","Ammonia","Ammonia","TP"))p <- corrplot(corr = env.cor,p.mat = env.p,order = "hclust",hclust.method = "ward.D",insig = "label_sig",sig.level = c(.05,0.01),tl.col = "black",pch.cex = 0.8)corrRect(p,namesMat = nameMat,col = "blue",lwd = 3)
# f图 绘制一个矩形,"|>"可以将corrplot()输出结果传递给corrRect(),要求R版本大于4.1.0。
nameMat <-c("pH","TP","TP","pH") # 设置矩形框的左、下、右和上侧的变量名称,即从左侧逆时针设置位置参数。if(getRversion() >= "4.1.0"){corrplot(corr = env.cor,p.mat = env.p,order = "hclust",hclust.method = "complete",insig = "label_sig",sig.level = c(.05,0.01),tl.col = "black",pch.cex = 0.8) |> corrRect(namesMat = nameMat,col = "red",lwd = 3) }
dev.off()
par(m)
注:corrRect.hclust()使用时method参数默认与corrplot设置的一致,k=3相当于addrect=3。帮助文件写可不在corrplot()设置order参数,直接使用corrRect()中的method参数设置。但是在实际操作中发现,必须在corrplot()中设置order和hclust.method参数,否则绘图行列变量名会与输入的一致,但是矩形框却是按照自定义的聚类方法形成的,导致绘图错误。所以建议直接在corrplot中设置聚类方法。R版本低于4.1.0可以使用magrittr包的“%>%”替代”|>“。
图7|添加矩形框热图。a: addrect参添加矩形框,此参数只能在type="full"以及没有p值矩阵的时候;b:corrRect.hclust()添加矩形框,需要与corrplot()配合使用,可以可视化p值;c:corrRec()使用index添加矩形框;d:corrRec()使用name添加矩形框;e:corrRec()使用namesMat添加矩形框;f:corrRec()使用namesMat添加矩形框,使用“|>”将corrplot()输出结果传递给corrRect(),要求R版本大于4.1.0。
3.2.3 可视化非相关性矩阵、NA和数学表达式标签
使用is.corr=FALSE和col.lim参数绘制非相关性矩阵,col.lim是按col设置颜色的范围。如果矩阵的行列变量数量不一致,绘图时可以使用win.asp调整图形纵横比,使图形成正方形形式。除非method为'circle'或'square',否则win.asp的值一般不大于1。corrplot()默认用?表示NA,调整na.label参数,可以用最多2个字符表示NA。使用"$"激活plotmath,可以在绘图时设置变量名为数学表达式。
m = par(no.readonly = TRUE) # 保存默认图形设置参数
pdf("dist.pdf",width=10,height = 10,family="Times") # 保存图片到本地
par(mfrow=c(2,2)) # 设置2行,2列的图形矩阵
# 可视化非相关性数据矩阵
## 基于环境因子数据计算样品间euclidean距离
env.dist=as.matrix(vegdist(ENV[1:10,4:14],method = "euclidean")) # 计算前10个样品的euclidean距离,并形成矩阵格式env.dist
## a图,绘制样本环境因子距离热图
?hcl.colors
hcl.pals() # 查看存在的颜色板名称
col3 = hcl.colors(12, "YlOrRd", rev = TRUE) # 选择了YlOrRD颜色板的12个颜色,并反转了颜色顺序。设置颜色适合[-N, 0]或[0, N]数据格式。corrplot(env.dist, is.corr = FALSE, col.lim=c(0,60),col = col3, win.asp = 1) #col.lim是按col设置颜色的范围,默认是[min(env.dist), max(env.dist)],自定义的范围需要超过数据矩阵的范围。
# 可视化NA
## 图b,na.label默认参数
diag(env.dist) = NA # 将env.dist数据矩阵的对角线值设置为NA。
corrplot(env.dist, is.corr = FALSE, col.lim=c(0,60),
col = col3, win.asp = 1) # 行列变量数不一样,可以设置win.asp值,调整图的纵横比。
## 图c,设置na.label="NA"
diag(env.dist) = NA # 将env.dist数据矩阵的对角线值设置为NA。
corrplot(env.dist, is.corr = FALSE, col.lim=c(0,60),
col = col3, win.asp = 1,na.label="NA")
# 可视化数学公式标签
## "$"激活plotmath,设置变量名为数学表达式。
colnames(env.dist) = rep(c('$alpha+beta', '$alpha[0]', '$alpha[beta]'),
c(3, 3, 4))
rownames(env.dist) = rep(c('$Sigma[i]^n', '$sigma', '$alpha[0]^100', '$alpha[beta]'), c(2, 3, 2, 3))## 图d,绘制变量名为数学表达式热图corrplot(env.dist, is.corr = FALSE, col.lim=c(0,60),col = col3, win.asp = 1,na.label="NA")dev.off() # 保存图片到本地par(m) # 恢复默认图形设置参数
图8|样本距离矩阵。
图9|corrplot的type参数修改热图展示方式。a:upper; b: lower。
扫描上方二维码,关注EcoEvoPhylo数据分析小店,购买数据分析书籍和服务。corrplot包推文就写到这里,其实还有很多函数可以设置,感兴趣的话,可以查看帮助文件进行学习。公众号后台发送“corrplot heatmap”获取数据和全部代码。
参考文献
http://www.datavis.ca/papers/corrgram.pdf
http://127.0.0.1:24605/library/corrplot/doc/corrplot-intro.html#reorder-a-correlation-matrix
推荐阅读
R中进行单因素方差分析并绘图
R绘图-相关性分析及绘图
R绘图-相关性系数图
R统计绘图-环境因子相关性热图
R统计绘图-corrplot绘制热图及颜色、字体等细节修改
R绘图-RDA排序分析
R统计绘图-RDA分析、Mantel检验及绘图
R统计-VPA分析(RDA/CCA)
R统计-PCA/PCoA/db-RDA/NMDS/CA/CCA/DCA等排序分析教程
开放转载
公众号原创文章开放转载,在文末留言告知即可
EcoEvoPhylo:主要分享微生物生态和phylogenomics的数据分析教程。扫描上方二维码,即可关注EcoEvoPhylo。让我们大家一起学习,互相交流,共同进步。