终于你也可以做这张图了

其实最近无论是实验还是写作任务都很繁重，夜深人静的时候，总是在想明天会好一些吗？可能会的吧。

前一段时间，有个图形被大家讨论过很多次，也出现了许多的绘制方式。宏基因组也曾转载过一篇文章，讨论了这种图形的实现方式。Science组合图表解读。

就是下面这张图片

许多人对这种图形进行了探索

R语言之照猫画虎1

R语言之照猫画虎2

Science组合图表解读

总的来说这些推送对这张图形的实现由两种思路。1，使用基础绘图包，2，使用ggplot。似乎基础包更容易一些，使用corplot来绘制相关图形，然后添加其他内容，第三篇推送还做了一些图形和含义的解读。

为什么我还要自己探索呢？有以下几个原因：

1. 我不愿意使用基础包实现这种图形的绘制，而以上推送就ggplot的实现方式距离实际运用还有一段距离。需要简化出图。

2. 这张图用于微生物生态方面很适用，结合mantel检验和相关分析的部分我用得到。

很久没和大家见面了，答应大家的事情很多，完成的却总也不够。今天我把这幅图带给大家。虽然在之前有预告，但是却没有告诉大家数据输入格式，今天你们将看到封装好的函数，我将提供示例数据，并且可以一键出图，让这种组合图形也难不住我们努力奋进的人。

下面便是

两份数据格式很简单，只要整理成这个形式即可一键出图

env.st = read.csv("./env.st.csv",row.names = 1)
env.st


report = read.csv("./report.csv",row.names = 1)
report


运行下面一条代码即可

p = plot_mantel_cor(env = env.st,report = report,title = "好久不见-我的朋友" )
p

就是下面这个函数

其实有很多参数我没有放进去，比如组合图形的相对位置，标签位置，颜色设置，线条粗细，额！是不是太多参数了，反正都是ggplot出图，大家自己加上吧，如果有朋友完善了这个函数，记得发给我哦。


plot_mantel_cor = function(env = env.st,report = report,title = "16S_microbiology"){
library(ggcorrplot)
library(igraph)
library(psych)

occor = corr.test(env.st,use="pairwise",method="spearman",adjust="fdr",alpha=.05)
occor.r = occor$r # 取相关性矩阵R值

occor.p = occor$p # 取相关性矩阵p值
occor.r[occor.p>0.05] = 0
# write.csv(occor.r,"./occor.r.csv",quote = F)
library(reshape2)
occor.r2 = occor.r[lower.tri(occor.r, diag = TRUE)]
length(colnames(occor.r))

### 构造环境变量三角矩阵图
## dd：构建三交矩阵横坐标
cc = rep(length(colnames(occor.r)),(length(colnames(occor.r))))
for (i in 1:(length(colnames(occor.r))-1)) {
ee = rep(length(colnames(occor.r))-i,(length(colnames(occor.r))-i))
dd = c(cc,ee)

cc = dd
}
dd = (length(colnames(occor.r))+1) - dd
##ww:构建三角矩阵纵坐标
gg = rep(length(colnames(occor.r)):1)
for (i in 1:(length(colnames(occor.r))-1)) {
ee = c((length(colnames(occor.r))-i):1)
ww = c(gg,ee)

gg = ww
}

##构造变量对应关系
wwt =data.frame(x = dd,y = ww)
##提取相关值大小，这是下三角矩阵转化
# xa = melt(occor.r)
wwt$mantelR = occor.r2
wwt

### 下面添加群落和环境因子的相关结果
##构造每个环境因子的坐标
x = c(1:(length(colnames(occor.r))) )
y = c((length(colnames(occor.r))+1) :2)
data2 = data.frame(x = x,y = y)
data2
###添加线的粗细和颜色映射
data2$count = report$R
data2$count = as.character(data2$count)
data2$count = as.numeric(data2$count)
data2

as = rep("a",length(data2$count))
for (i in 1:length(data2$count)) {
if (data2$count[i] > 0) {
as[i] = "+"
}

if (data2$count[i] < 0) {
as[i] = "-"
}
if (data2$count[i] == 0) {
as[i] = "-"
}
}
as

data2$group = as
data2$label = colnames(env.st)
wwt

data2$count1 = data2$count^2*500
# data2$count1 = c(1:length(data2$count))
##这和图形整体绘制

###geom_curve添加size参数如果在aes内部的话默认线非常宽，但是size放到sea外就不会这样，这里我选择放到外面
p = ggplot() +
geom_tile(aes(x = x, y = y),wwt,fill=NA,color='gray',size=0.5)+
geom_point(aes(x = x, y = y,size=mantelR,fill=mantelR),wwt, shape=22,color='white')+
scale_size(range = c(1, 8))+
scale_fill_distiller(palette="RdYlBu")+
geom_curve(aes(x = max(wwt$x)*3/4, y = max(wwt$y)*3/4, xend = x, yend = y,group = group,color = group),curvature = 0.2,data2,size = data2$count1) +
geom_point(aes(x = x, y = y),pch = 21,size =4,data2,color = "black",fill = "#FFF5EB")+
geom_point(aes(x = max(wwt$x)*3/4, y = max(wwt$y)*3/4),pch = 21,size = 6,color = "black",fill = "#FEE6CE")

p

### 添加环境因子添加标签
p = p + geom_text(aes(x = x, y = 0,label= colnames(env.st)),size=4,data2) +
geom_text(aes(x = 0, y = y-1,label= colnames(env.st)),size=4,data2)
p = p + geom_text(aes(x = x+0.5, y = y+0.5,label=data2$label),size=4,data2)
##添加群落矩阵标签
wwt

# max(wwt$x)*3/4
# max(wwt$y)*3/4

p = p +
geom_text(aes(x = max(wwt$x)*3/4, y = max(wwt$y)*3/4,label= title),size = 6)
p

### 修改主题
p = p +theme_void()#横纵坐标群去掉

p

}