🔥 热搜 🔥
上海习近平新疆鄂州父女瓜乌鲁木齐疫情H工口小学生赛高习明泽芊川一笑图包印尼排华
分类
社会娱乐国际人权科技经济其它
🔥 热搜 🔥
百度今日热点微信公众平台贴吧opggdnf私服百度贴吧知乎dnf公益服百度傻逼
分类
社会娱乐国际人权科技经济其它
查看原文
其他

science 级别的组合图表函数升级了

文涛聊科研 微生信生物 2022-07-05


有朋友层问询是否可以升级这个图表，当然可以，晚饭小憩一会，看着日落，为大家带来这次升级。

这是之前写的函数


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

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

source("./plot_mantel_cor.R")

p = plot_mantel_cor(env = env.st,report = report,title = c("A") )
p

下面是升级后的组合图表

这是新的函数


source("./plot_mantel_cor_3_more.R")

当单个矩阵的时候和之前一样

report = read.csv("./report.csv",row.names = 1)
p = plot_mantel_cor_3_more(env = env.st,report = report,title = c("A") )
p

两个矩阵和因子做相关

这里我模拟了一个新的列，代表新的矩阵同环境因子做相关。


report = read.csv("./report.csv",row.names = 1)
report$aa = report$R
# report$bb = report$R
p = plot_mantel_cor_3_more(env = env.st,report = report,title = c("A","B") )
p

三个矩阵和因子做相关

无论是两个还是三个或者四个矩阵同因子做相关，只需要在report后面多添加几列就够了，名字随便命名。我在函数中都做了修改。

report = read.csv("./report.csv",row.names = 1)
report$aa = report$R
report$bb = report$R
p = plot_mantel_cor_3_more(env = env.st,report = report,title = c("A","B","C") )
p

四个矩阵和因子做相关

report = read.csv("./report.csv",row.names = 1)
report$aa = report$R
report$bb = report$R
report$cc = report$R
p = plot_mantel_cor_3_more(env = env.st,report = report,title = c("A","B","C","D") )
p

下面附函数

如果大家有时间，可以添加方向代码，让这个图形可以随意调换方向。图例可以精修一下。




######环境变量之间的关系··································································


# title = c("A","B","c")
plot_mantel_cor_3_more = function(env = env.st,report = report,title = "16S_microbiology"){
colnames(report) = c(colnames(report)[1],paste("R",1:(dim(report) [2] - 1),sep = ""))
report
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

# p = ggcorrplot(occor.r,p.mat = occor.p, method = "circle",lab = TRUE,insig = "blank",outline.color = "white",
# ggtheme = ggplot2::theme_bw())
#
# p
#
###############



# write.csv(occor.r,"./occor.r.csv",quote = F)
library(reshape2)

occor.r2 = occor.r[lower.tri(occor.r, diag = TRUE)]
# occor.r2 = occor.r[lower.tri(occor.r, diag = TRUE)]
length(colnames(occor.r))
# a = rep(0,length(occor.r2))
# i = 1
### 构造环境变量三角矩阵图

## 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


# cc = rep(1:length(colnames(occor.r)))
# for (i in 1:(length(colnames(occor.r))-1)) {
# ee = c(1:(length(colnames(occor.r))-i))
# rr = c(cc,ee)
#
# cc = rr
# }
#
# rr
##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
}

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 = cbind(data2,as.data.frame(report[2:dim(report)[2]]^2*500))
data2

axa = as.matrix(report[2:dim(report)[2]])
axa[axa>0] <- 1
axa[axa<0] <- -1
axa[axa==0] <- 0
colnames(axa) = paste(colnames(axa),"d",sep = "")
data2 = cbind(data2,axa)
data2


data2$label = colnames(env.st)
data2


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")
p
##如果三角之外只有一个点

if (dim(report) [2] == 2) {
p = p+
geom_curve(aes(x = max(wwt$x)*3/4, y = max(wwt$y)*3/4, xend = x, yend = y,group = as.factor(data2$R1d),color = as.factor(data2$R1d)),curvature = 0.2,data2,size = data2$R1) +
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")




}
## 有两个群落
if (dim(report) [2] == 3) {
p = p+
geom_curve(aes(x = max(wwt$x)*2/3, y = max(wwt$y), xend = x, yend = y,group = as.factor(data2$R1d),color = as.factor(data2$R1d)),curvature = 0.2,data2,size = data2$R1) +
geom_point(aes(x = x, y = y),pch = 21,size =4,data2,color = "black",fill = "#FFF5EB")+
geom_point(aes(x = max(wwt$x)*2/3, y = max(wwt$y)),pch = 21,size = 6,color = "black",fill = "#FEE6CE")

p = p+
geom_curve(aes(x = max(wwt$x)*4.2/5, y = max(wwt$y)*3/4, xend = x, yend = y,group = as.factor(data2$R2d),color = as.factor(data2$R2d)),curvature = 0.2,data2,size = data2$R2) +
geom_point(aes(x = x, y = y),pch = 21,size =4,data2,color = "black",fill = "#FFF5EB")+
geom_point(aes(x = max(wwt$x)*4.2/5, y = max(wwt$y)*3/4),pch = 21,size = 6,color = "black",fill = "#FEE6CE")

}
## 有三个群落
if (dim(report) [2] == 4) {
p = p+
geom_curve(aes(x = max(wwt$x)*1.5/3, y = max(wwt$y*1.5), xend = x, yend = y,group = as.factor(data2$R1d),color = as.factor(data2$R1d)),curvature = 0.2,data2,size = data2$R1) +
geom_point(aes(x = x, y = y),pch = 21,size =4,data2,color = "black",fill = "#FFF5EB")+
geom_point(aes(x = max(wwt$x)*1.5/3, y = max(wwt$y*1.5)),pch = 21,size = 6,color = "black",fill = "#FEE6CE")
p
p = p+
geom_curve(aes(x = max(wwt$x)*3.5/4, y = max(wwt$y), xend = x, yend = y,group = as.factor(data2$R2d),color = as.factor(data2$R2d)),curvature = 0.2,data2,size = data2$R2) +
geom_point(aes(x = x, y = y),pch = 21,size =4,data2,color = "black",fill = "#FFF5EB")+
geom_point(aes(x = max(wwt$x)*3.5/4, y = max(wwt$y)),pch = 21,size = 6,color = "black",fill = "#FEE6CE")


p = p+
geom_curve(aes(x = max(wwt$x)*4.8/5, y = max(wwt$y)*2/4, xend = x, yend = y,group = as.factor(data2$R3d),color = as.factor(data2$R3d)),curvature = 0.2,data2,size = data2$R3) +
geom_point(aes(x = x, y = y),pch = 21,size =4,data2,color = "black",fill = "#FFF5EB")+
geom_point(aes(x = max(wwt$x)*4.8/5, y = max(wwt$y)*2/4),pch = 21,size = 6,color = "black",fill = "#FEE6CE")


}
## 有4个群落
if (dim(report) [2] == 5) {
p = p+
geom_curve(aes(x = 2, y = max(wwt$y*1.5), xend = x, yend = y,group = as.factor(data2$R1d),color = as.factor(data2$R1d)),curvature = 0.2,data2,size = data2$R1) +
geom_point(aes(x = x, y = y),pch = 21,size =4,data2,color = "black",fill = "#FFF5EB")+
geom_point(aes(x = 2, y = max(wwt$y*1.5)),pch = 21,size = 6,color = "black",fill = "#FEE6CE")
p
p = p+
geom_curve(aes(x = 5, y = max(wwt$y*1.4), xend = x, yend = y,group = as.factor(data2$R2d),color = as.factor(data2$R2d)),curvature = 0.2,data2,size = data2$R2) +
geom_point(aes(x = x, y = y),pch = 21,size =4,data2,color = "black",fill = "#FFF5EB")+
geom_point(aes(x = 5, y = max(wwt$y*1.4)),pch = 21,size = 6,color = "black",fill = "#FEE6CE")
p

p = p+
geom_curve(aes(x = max(wwt$x)*3.5/4, y = max(wwt$y), xend = x, yend = y,group = as.factor(data2$R3d),color = as.factor(data2$R3d)),curvature = 0.2,data2,size = data2$R3) +
geom_point(aes(x = x, y = y),pch = 21,size =4,data2,color = "black",fill = "#FFF5EB")+
geom_point(aes(x = max(wwt$x)*3.5/4, y = max(wwt$y)),pch = 21,size = 6,color = "black",fill = "#FEE6CE")

p
p = p+
geom_curve(aes(x = max(wwt$x)*4.8/5, y = max(wwt$y)*2/4, xend = x, yend = y,group = as.factor(data2$R4d),color = as.factor(data2$R4d)),curvature = 0.2,data2,size = data2$R4) +
geom_point(aes(x = x, y = y),pch = 21,size =4,data2,color = "black",fill = "#FFF5EB")+
geom_point(aes(x = max(wwt$x)*4.8/5, y = max(wwt$y)*2/4),pch = 21,size = 6,color = "black",fill = "#FEE6CE")

p
}


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)
p
##添加群落矩阵标签
if (dim(report) [2] == 2) {
p = p +
geom_text(aes(x = max(wwt$x)*3/4+1, y = max(wwt$y)*3/4+1,label= title),size = 6)
p

}
# title = c("A","B","c")
if (dim(report) [2] == 3) {
p = p +
geom_text(aes(x = max(wwt$x)*2/3 +1, y = max(wwt$y),label= title[1]),size = 6) +
geom_text(aes(x = max(wwt$x)*4.2/5 +1, y = max(wwt$y)*3/4,label= title[2]),size = 6)

p

}
if (dim(report) [2] == 4) {
p = p +
geom_text(aes(x = max(wwt$x)*1.5/3+1, y = max(wwt$y*1.5)+1,label= title[1]),size = 6) +
geom_text(aes(x = max(wwt$x)*3.5/4+1, y = max(wwt$y)+1,label= title[2]),size = 6) +
geom_text(aes(x = max(wwt$x)*4.8/5+1, y = max(wwt$y)*2/4+1,label= title[3]),size = 6)

p

}
#四个群落
if (dim(report) [2] == 5) {
p = p +
geom_text(aes(x = 2+1, y = max(wwt$y*1.5)+1,label= title[1]),size = 6) +
geom_text(aes(x = 5+1, y = max(wwt$y*1.4)+1,label= title[2]),size = 6) +
geom_text(aes(x = max(wwt$x)*3.5/4+1, y = max(wwt$y)+1,label= title[3]),size = 6) +
geom_text(aes(x = max(wwt$x)*4.8/5+1, y = max(wwt$y)*2/4+1,label= title[4]),size = 6)

p

}





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



}



示例数据和代码可以在微生信生物群中得到。添加小编微信，带你入群。以后便随时得到一些有用的脚本。

或者后台留言，我会准备一个网盘链接，发送给你。

历史目录

R语言分析技术

  1. 《扩增子16s核心OTU挑选-基于otu_table的UpSet和韦恩图》

  2. 《分类堆叠柱状图顺序排列及其添加合适条块标签》

  3. 《R语言绘制带有显著性字母标记的柱状图》

  4. 《使用R实现批量方差分析(aov)和多重比较(LSD)》

  5. Rstudio切换挂载R版本及本地安装一些包

  6. 玩转R包

  7. 用ggplot画你们心中的女神

  8. ggplot版钢铁侠

  9. 想用ggplot做一张致谢ppt，但是碰到了520，画风就变了

扩增子专题

  1. 《16s分析之Qiime数据整理+基础多样性分析》

  2. 《16s分析之差异展示（热图）》

  3. 迅速提高序列拼接效率，得到后续分析友好型输入，依托qiime

  4. https://mp.weixin.qq.com/s/6zuB9JKYvDtlomtAlxSmGw》

  5. 16s分析之网络分析一(MENA)

  6. 16s分析之进化树+差异分析（一）

  7. 16s分析之进化树+差异分析（二）

  8. Qiime2学习笔记之Qiime2网站示例学习笔记

  9. PCA原理解读

  10. PCA实战

  11. 16s分析之LEfSe分析

  12. 16s分析之FishTaco分析

基于phyloseq的微生物群落分析

  1. 开年工作第一天phyloseq介绍

  2. phyloseq入门

  3. R语言微生物群落数据分析:从原始数据到结果（No1）

  4. R语言微生物群落数据分析:从原始数据到结果（No2））

  5. phyloseq进化树可视化

  6. 基于phyloseq的排序分析

代谢组专题

  1. 非靶向代谢组学数据分析连载（第零篇引子）

  2. 非靶向代谢组学分析连载（第一篇：缺失数据处理、归一化、标准化）

当科研遇见python

1.当科研遇见python

科学知识图谱

1.citespace 的基本术语与下载安装

杂谈

  1. 我的生物信息之路

  2. graphlan进一步学习笔记之进化树

  3. 如约 为大家带来了graphlan的物种分类树

欢迎加入微生信生物讨论群，扫描下方二维码添加小编微信，小编带你入伙啦，大牛如云，让交流变得简单。

哦对了，还记得之前我写了一个R包吗？好像只有一个功能

该整理代码封个R包了-来自复兴号的学习笔记

我刚才把这个函数和测试数据集放到easyMicrobiome包中了，额，还没有push·······

总之先放出代码，空了push到github，大家更方便安装指定测试。

library(easyMicrobiome)
data("env.st")
data("report")
#one
p = plot_mantel_cor_3_more(env = env.st,report = report,title = c("A") )
p


#two
report$aa = report$R
# report$bb = report$R
p = plot_mantel_cor_3_more(env = env.st,report = report,title = c("A","B") )
p
#three
report$aa = report$R
report$bb = report$R
p = plot_mantel_cor_3_more(env = env.st,report = report,title = c("A","B","C") )
p

#four
report$aa = report$R
report$bb = report$R
report$cc = report$R
p = plot_mantel_cor_3_more(env = env.st,report = report,title = c("A","B","C","D") )
p

或者提前添加小编微信，我把包发给你。

欢迎加入微生信生物讨论群，扫描下方二维码添加小编微信，小编带你入伙啦，大牛如云，让交流变得简单。


您可能也对以下帖子感兴趣

文章有问题？点此查看未经处理的缓存