【R画图学习5】mantel_test和相关性组合图

今天学习我在看微生物相关的文章的时候经常见到的一个图，这个图可以看出有两部分，一部分是各个环境因子的相关性图，另外一个是mantel_test的结果。所以我们就来学习一下如何得到下面的这样子的图。

先来看下怎么理解这个图：

一开始看到这个图，可能觉得挺酷炫，但是并不知道这个图有啥意义！所以作图之前一定要搞清楚你想表达什么，highlight什么，展示什么样子的结果，做这个图啥含义，这一点至关重要！

Mantel test 是对两个矩阵相关关系的检验，之所以抛开相关系数这样一种方法，是因为相关系数只能处理两列数据之间的相关性，而在面对两个矩阵之间的相关性时就束手无策。另外，普通的相关分析只能检验单个环境因子与单个微生物类群(例如单个属)之间的相关性，而基于距离矩阵的Mantel test可以检验单个或一组环境因子与整个微生物群落的相关性。Mantel test的相关性越大，p值越小，则说明环境因子对微生物群落的影响越大。这也是我自己目前看到过的paper中最多的应用场景。

当然，这个也有很多的应用场景：

- 微生物群落与生态环境变量之间的相关性；

- 人体微生物区与某疾病程度的相关性；

- 不同药物组合处理疾病后，微生物组成结构与病情改善之间的相关性；

Mantel test，顾名思义，是一种检验。既然是检验就得有原假设，它的原假设是两个矩阵之间没有相关关系。检验过程如下：

- 两个矩阵都对应展开，变成两列，计算相关系数（理论上什么相关系数都可以计算，但常用pearson相关系数），

- 然后其中一列或两列同时随机打乱，再计算一个值，反复打乱成千上万次，看打乱前的r值，在打乱后所得r值分布中的位置，如果跟随机置换得到的结果较近，则不大相关，反之则为显著。

今天的学习需要一下几个包：

library(tidyverse)

library(ggcor)

library(vegan)

先查看一下包里自带的示例数据：

varespec 描述了24块样地中44种植物的丰度信息

varechem 描述了这24块样地土壤的14个理化参数

当然，同样的，对于微生物的丰度我们也可以采用类似的数据形式呈现，并绘图。

data("varechem", package = "vegan")

varechem[1:5,1:10]

data("varespec", package = "vegan")

varespec[1:5,1:10]

dim(varechem)

dim(varespec)

我们先来做mantel test，这个的输入应该是2个矩阵，但是样品名应该是一样的。

mantel <- mantel_test(varespec, varechem,

                      spec.select = list(Spec01 = 1:7,

                                        Spec02 = 8:18,

                                        Spec03 = 19:37,

                                        Spec04 = 38:44))

大家在做检验时一定要根据实际情况，把spec.select设置好，比如这里的话，就是把输入的44个植物的峰度信息分成了4组，对应于varespec的44列。paper中常见的也有算细菌，真菌，代谢组和个环境理化指标的关联强度的。

从mantel的结果文件来看，主要包含4个group分别和每个环境理化指标的相关性指数r以及对应的p-value。

因为最终的图中对于r和pvalue是分了区间的，所以我们把r和pvalue变成离散变量。

mantel <- mantel %>% mutate(rd = cut(r, breaks = c(-Inf, 0.2, 0.4, Inf),

          labels = c("< 0.2", "0.2 - 0.4", ">= 0.4")),

          pd = cut(p.value, breaks = c(-Inf, 0.01, 0.05, Inf),

          labels = c("< 0.01", "0.01 - 0.05", ">= 0.05")))

其中：%>%是R中经常用的管道符。mutate是增加一个新的变量（我们新增加了rd和pd）。cut就是把一个连续的量离散化。

下面我们再学习一下绘制左上角相关性的图，主要是基于ggcor实现的。

quickcor(varechem, type = "lower")+geom_square()  

通过这个命令，我们就计算并绘制了varechem中各环境因子之间的相关系数，type控制绘制的区域和位置（默认的话就是全矩阵，指定的话，可以是lower，也可以是upper），geom_square()  控制绘制的单元格里面用方块展示，经常用的还有geom_circle2()，geom_star()等。

quickcor(varechem, type = "upper")+geom_circle2()

换一种展示方式。

quickcor(varechem)+geom_circle2()    //这个展示的就是一个对称的全矩阵。

也可以一半用颜色展示，一半用具体的数字展示。

quickcor(varechem, cor.test=TRUE)+

geom_square(data=get_data(type="lower",show.diag=FALSE)) +

geom_mark(data=get_data(type="upper",show.diag=FALSE),size=2.5)+

geom_abline(slope=-1,intercept=15)

回到我们原先要画的图上，学会了画相关性的图，下面就是要把mantel_test的结果放上去了。在这包中，是通过anno_link添加上去的，和ggplot挺像的，都是不断的望上面添加。我们设置是让颜色跟着pd也就是pvalue变化，线条的大小跟着r来变化。

quickcor(varechem, type = "lower") +

 geom_square() +

 anno_link(aes(colour = pd, size = rd), data = mantel)

线的颜色看起来好像有点太粗了，我们调整一下线的粗细。

quickcor(varechem, type = "lower") +

geom_square() +

 anno_link(aes(colour = pd, size = rd), data = mantel) +

 scale_size_manual(values = c(0.5, 1, 2))

另外哦我们想要把legend的pd和rd放在一起，所以需要调整legend的顺序。图例的调整一般是通过guides来调整的。

quickcor(varechem, type = "lower") +

geom_square() +

anno_link(aes(colour = pd, size = rd), data = mantel) +

scale_size_manual(values = c(0.5, 1, 2))+

guides(size = guide_legend(title = "Mantel's r",order = 2), 

        colour = guide_legend(title = "Mantel's p",order = 1),

        fill = guide_colorbar(title = "Pearson's r", order = 3))

这样我们就调整了图例的顺序。

别人paper中，好像是曲线，我们也可以换成曲线。主要是在anno_link中的curvature参数调整。

quickcor(varechem, type = "lower") +

  geom_square() +

  anno_link(aes(colour = pd, size = rd), data = mantel,curvature = -0.2) +

  scale_size_manual(values = c(0.5, 1, 2))+

  guides(size = guide_legend(title = "Mantel's r",order = 2), 

        colour = guide_legend(title = "Mantel's p",order = 1),

        fill = guide_colorbar(title = "Pearson's r", order = 3))

比如我的审美点的话，我还不太喜欢这个默认色系，所以可能还需要自己修改色系。

quickcor(varechem, type = "lower") +

geom_square() +

 anno_link(aes(colour = pd, size = rd), data = mantel,curvature = -0.2) +

 scale_size_manual(values = c(0.5, 1, 2))+

 scale_colour_manual(values = c("#d85c01", "#29d300", "#A2A2A288")) +

  guides(size = guide_legend(title = "Mantel's r",order = 2), 

        colour = guide_legend(title = "Mantel's p",order = 1),

        fill = guide_colorbar(title = "Pearson's r", order = 3))

我通过scale_colour_manual修改了线条的颜色。

我们还可能想换相关系数的色系，这个量是个连续变量。一般情况下，修改颜色属性：gradientn -- 渐变色； manul -- 分类变量颜色。对于scale系列的函数，我也掌握的比较迷糊，查了一下资料，大概分为这几类：

scale_alpha_*() 【设置透明度】

scale_color_*() 或 scale_colour_*() 【设置边框/散点颜色】

scale_fill_*（) 【设置填充颜色和图例相关内容】

scale_linetype_*() 【设置线条样式】

scale_shape_*() 【设置散点样式】

scale_size_*() 【设置散点/文本大小】

scale_radius_*() 【设置散点半径大小】

scale_x_*() 【设置横坐标相关参数】

scale_y_*() 【设置纵坐标相关参数】

所以从资料可以看出，对于颜色类的通知一般是通过scale_fill_*系列来实现的。

ls("package:ggplot2", pattern = "^scale_fill.+") 

可以看出，即便fill系列也是很多类的。具体每一类的用法我也不是很熟悉，只能边用边查了。资料显示适合连续型变量的有continuous和gradient等系列，好像scale_fill_distiller() 也可以。

对于数据为非因子型的填充色映射，ggplot2自动使用“continuous”类型颜色标尺表示连续颜色空间。

quickcor(varechem, type = "lower") +

  geom_square() +

  anno_link(aes(colour = pd, size = rd), data = mantel,curvature = -0.2) +

  scale_size_manual(values = c(0.5, 1, 2))+

  scale_colour_manual(values = c("#d85c01", "#29d300", "#A2A2A288")) +

  guides(size = guide_legend(title = "Mantel's r",order = 2), 

        colour = guide_legend(title = "Mantel's p", order = 1),

        fill = guide_colorbar(title = "Pearson's r", order = 3))+

  scale_fill_continuous(low = "blue", high = "red", space = "rgb")

我们换了蓝红色系。

连续填充色设置函数还有scale_fill_gradient，scale_fill_gradient2和 scale_fill_gradientn，其中scale_fill_gradient的用法和作用和scale_fill_continuous完全相同（其实ggplot2早期版本连续颜色标尺默认使用scale_fill_gradient，没有scale_fill_continuous函数；后者可能是H.W头脑清楚以后加进去的，相当于前者的别名）。scale_fill_gradient2增加了中间点和中间颜色的设置

quickcor(varechem, type = "lower") +

  geom_square() +

  anno_link(aes(colour = pd, size = rd), data = mantel,curvature = -0.2) +

  scale_size_manual(values = c(0.5, 1, 2))+

  scale_colour_manual(values = c("#d85c01", "#29d300", "#A2A2A288")) +

  scale_fill_gradient2(low = "blue", high = "red", mid = "white", midpoint = 0)+

  guides(size = guide_legend(title = "Mantel's r",order = 2), 

        colour = guide_legend(title = "Mantel's p", order = 1),

        fill = guide_colorbar(title = "Pearson's r", order = 3))

用起来貌似比continuous方便些，因为可以设置多个颜色区间。

我们也可以修改性状为自己想要的：

quickcor(varechem, type = "lower", show.diag = FALSE) +

  geom_star(n=6) + 

  anno_link(aes(colour = pd, size = rd), data = mantel,curvature = -0.2) + 

  scale_size_manual(values = c(0.5, 1, 2))+

  scale_fill_gradient2(low = "blue", high = "red", mid = "white", midpoint = 0)+

  scale_colour_manual(values = c("#d85c01", "#29d300", "#A2A2A288")) +

  guides(size = guide_legend(title = "Mantel's r",order = 2), 

        colour = guide_legend(title = "Mantel's p", order = 1), 

        fill = guide_colorbar(title = "Pearson's r", order = 3))