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