Science：组合图表绘制

写在前面

有些事情做了很痛苦，但是不做更难受。比如今天谈到的这件事情，我花了一天加一个晚上时间才写了一个这个功能：Science组合图表的实现。这部分一共写了三个函数，一个门特尔检验，一个ggplot版本的相关热图绘制，还有连线图形绘制。为什么我花费了这么多时间呢？一方面这个功能的实现要想做的灵活点，就必须有一个ggplot版本的相关热图轮子，公开的轮子不好用，所以我自己造了一个，功能很基础哈，只是一个雏形，但是可以作为后面组合图表的轮子。这几个独立的函数也很有用，但是不细化了。

早在之前我们就推送过science组合图表的推文了，我们做的也是比较早的，ggplot版本的，例如；

最后最好用的出现了，虽然在分析微生物数据上还有点小瑕疵，但是已经非常舒服了，但是这个包已经被作者删除，厚哥表示不会在公开这个包。半年来大家的需求似乎不断，本来从第四个版本以来我放弃了这个图形的开发，因为开发的没有厚哥的好，所以就搁置了一年，到如今，我只能捡起来继续做一些工作，尽量让这个功能好用点。还有就是任何问题的解决方案都不应该只有一个，太容易就断掉了，所以增加一种问题的解决方法，虽然不一定是最好的，也是有必要的。我也没有参考厚哥之前的代码，自己实现，达到一个功能也是挺有意思的，就是浪费了好多时间，代码也不够成熟，慢慢提升吧，实用性后再优化。

MatCorPlot函数参数解释

• env.dat = env1 其他指标表格，参见要求

• tabOTU = tabOTU1 otu表格，list文件，最多三个。

• distance = “bray”群落距离计算方法，这里就两种 bray和jcd。

• method = “metal” 门特尔检验方法， 这里两种，另外一种是偏门特尔检验。

• x = F 是否设置x轴标签。

• y = F 是否设置y轴标签

• diag = T 是否设置对角线标签。

• sig = TRUE 是否仅仅展示显著的色块

• siglabel = FALSE 显著色标是否标记星号

• shownum = TRUE 是否展示相关值

• numpoint = NULL 色块point类型，默认22号方块，直接使用pch编号即可更换。注意必须更换具有fill属性的编号。

• numsymbol = 27 这是symbol，扩展了ggplot中point类型，如果要使用，则需要另外安装R包ggsymbol。需要github安装

• numpoint2 这部分只可以使用ggplot中的pch点形状，默认21号点，控制link部分点形状。

• lacx = “left” 相关热图展示方向，x轴，可选“right”

• lacy = “bottom” 相关热图展示方向，y轴,可选“top”

• range = 0.5 连线整体粗细变化，值越大变化越大。

• p.thur = 0.3 显著性阈值

• onlysig = F 是否只展示显著性的连线。

实战

注意这个功能为2020年11月29号加入，之前安装的没有这个功能。

devtools:: install_github("taowenmicro/ggClusterNet")

需要R包和数据（内置）

library(phyloseq)
library(ggClusterNet)
library(tidyverse)
data(ps)
ps1 = filter_taxa(ps, function(x) sum(x ) > 200 , TRUE);ps1

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 442 taxa and 18 samples ]
## sample_data() Sample Data:       [ 18 samples by 10 sample variables ]
## tax_table()   Taxonomy Table:    [ 442 taxa by 7 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 442 tips and 441 internal nodes ]

otu = as.data.frame(t(vegan_otu(ps1)))
mapping = as.data.frame( sample_data(ps1))
data(env1)

关于微生物组数据的要求：要求很简单，将otu表格作为list对象即可，就像下这样，多个optu表格自然都添加进去list对象即可。

其他指标数据要求：一张数据框，行名是样本，列名字是指标。

tabOTU1 = list(otu1 = otu)
data(env1)

默认参数会会给出一套解决方案，当然可以换

# ?MatCorPlot
p0 <- MatCorPlot(env.dat = env1,tabOTU = tabOTU1)
p0

numpoint2改变link部分点形,并标记显著性，去除不显著的元素，标记显著的R值。

p1 <- MatCorPlot(
  env.dat = env1,
  tabOTU = tabOTU1,
  numpoint2 = 22,
  sig = F,
  lacx = "right", # 改变位置x
  lacy = "bottom",# 改变位置y
  onlysig = F # 是否只显示显著的link线
  
)

p1

curvature 修改link的弯曲程度。越小线越直。去除数字，标记为星号，其次换point形状为新形势

p2 <- MatCorPlot(
  env.dat = env1,
  tabOTU = tabOTU1,
  siglabel = T, # 使用星号标记显著的关系
  shownum = F, # 不显示数字
  numpoint = NULL,# 不使用ggplot吗，默认点
  numpoint2 = 13,#link使用的point为pch = 13
  numsymbol = 27,# 使用扩展点。
  curvature = 0,
  lacx = "right",
  lacy = "bottom",
  p.thur = 0.3,
  sig = F,
  onlysig = F
)
p2

我们更换图形角度和去除不显著的连线，这里不显著设置为0.3，保证去除一部分线，模拟数据测试用，不作为参考值。

p4 <- MatCorPlot(
  env.dat = env1,
  tabOTU = tabOTU1,
  distance = "bray",
  method = "metal",
  x = F,
  y = F,
  diag = T,
  sig = TRUE,
  siglabel = FALSE,
  shownum = TRUE,
  numpoint = 22,
  numsymbol = NULL,
  lacx = "left",
  lacy = "bottom",
  range = 1,
  p.thur = 0.3,
  onlysig = T
  
)
p4

即使是两个群落

tabOTU1 = list(otu1 = otu,otu2 = otu)
data(env1)
p5 <- MatCorPlot(
  env.dat = env1,
  tabOTU = tabOTU1,
  distance = "bray",
  method = "metal",
  x = F,
  y = F,
  diag = T,
  sig = TRUE,
  siglabel = FALSE,
  shownum = TRUE,
  numpoint = NULL,
  numsymbol = 27,
  lacx = "right",
  lacy = "bottom",
  range = 0.5,
  p.thur = 0.3,
  onlysig = F
  
)
p5

最多三个群落（细菌，真菌，古菌）

我就设置了三个群落，如果后面有需要，可以设置成为无限个，但是也没有意义。后面如果有时间，我会升级成6个群落，这就差不多了，不然也太乱了。

tabOTU1 = list(otu1 = otu,otu2 = otu,otu3 = otu)
p6 <- MatCorPlot(
  env.dat = env1,
  tabOTU = tabOTU1,
  distance = "bray",
  method = "metal",
  x = F,
  y = F,
  diag = T,
  sig = TRUE,
  siglabel = FALSE,
  shownum = TRUE,
  numpoint = NULL,
  numsymbol = 27,
  lacx = "left",
  lacy = "bottom",
  range = 0.5,
  p.thur = 0.3,
  onlysig = F
  
)

p6

library(patchwork)

(p0 |  p1 | p2) / (p4 | p5 | p6)

对于颜色和填充的改变

library(RColorBrewer)#调色板调用包
m=brewer.pal(9,"YlGn")
library(scales)
m=c("#E41A1C" ,"#377EB8")
show_col(m)

p + scale_colour_manual(values=c("red","blue"))p + scale_colour_manual(values = m) + scale_fill_continuous()

# p + scale_colour_manual(values = m) + scale_fill_distiller(palette="RdYlBu")

p + scale_colour_manual(values = m) + scale_fill_distiller(palette="PRGn")

p + scale_colour_manual(values = m) + scale_fill_distiller(palette="BrBG")(p0 | p7 | p8) / (p0 | p9 | p10)

猜你喜欢

系列教程：微生物组入门 Biostar 微生物组  宏基因组

专业技能：学术图表 高分文章 生信宝典 不可或缺的人

一文读懂：宏基因组 寄生虫益处 进化树

必备技能：提问 搜索  Endnote

文献阅读 热心肠 SemanticScholar Geenmedical

扩增子分析：图表解读 分析流程 统计绘图

在线工具：16S预测培养基 生信绘图

科研经验：云笔记  云协作 公众号

编程模板: Shell  R Perl

生物科普:  肠道细菌 人体上的生命 生命大跃进  细胞暗战 人体奥秘  

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外5000+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。PI请明示身份，另有海内外微生物相关PI群供大佬合作交流。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍未解决群内讨论，问题不私聊，帮助同行。

学习16S扩增子、宏基因组科研思路和分析实战，关注“宏基因组”

点击阅读原文，跳转最新文章目录阅读