R语言相关分析热图代码 r语言绘制相关性热图

上一篇文章推送的是怎样调整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 Roman
par(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)  # 恢复默认图形设置参数

图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)

图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|样本距离矩阵。