今天学习我在看微生物相关的文章的时候经常见到的一个图，这个图可以看出有两部分，一部分是各个环境因子的相关性图，另外一个是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))
网友评论