【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))

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 200,045评论 5 468
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 84,114评论 2 377
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 147,120评论 0 332
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 53,902评论 1 272
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 62,828评论 5 360
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,132评论 1 277
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,590评论 3 390
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,258评论 0 254
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,408评论 1 294
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,335评论 2 317
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,385评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,068评论 3 315
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,660评论 3 303
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,747评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,967评论 1 255
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,406评论 2 346
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 41,970评论 2 341

推荐阅读更多精彩内容