本文首发于公众号：医学和生信笔记

“

医学和生信笔记，专注R语言在临床医学中的使用，R语言数据分析和可视化。主要分享R语言做医学统计学、meta分析、网络药理学、临床预测模型、机器学习、生物信息学等。

在几年前出现了一个ggcor包，可以用来可视化mantel test的结果，最开始还可以通过cran安装，不过后来也不行了，而且这个包由于一些原因已经停止维护了，最近的更新是2年前了！

但是那张图却一直很风靡。。。其实原作者已经开发了新的包用于可视化mantel test，名字叫linkET，只是由于缺少宣传，大家知道的比较少。

善于搜索一搜就能搜到，我在之前的 可能是最适合初学者的R包安装教程，视频中提到了这个linkET，但是大家不愿意看，一个劲的问封面图是怎么画的，我真是服了。。。

所以，今天专门介绍用于mantel test可视化的linkET包。

安装

首先是安装R包，这个包只能通过github安装，或者下载到本地安装，使用install.packages()100%报错，使用BiocManager::install()也是100%报错！

很多初学者最大的拦路虎绝对是R包安装，有些人宁愿花钱找tb，也不愿意自己学习一下，搞不懂！

# install.packages("devtools")devtools::install_github("Hy4m/linkET", force = TRUE)

使用

一般的相关性分析是用于两列数据之间的，而mantel test 是用于两个矩阵的相关性检验，我在工作中用的很少，做微生物的小伙伴应该用的比较多，做肠道菌群、宏基因组的应该也会用到。

加载R包和数据：

library(linkET)library(dplyr)## ## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':## ##     filter, lag## The following objects are masked from 'package:base':## ##     intersect, setdiff, setequal, unionlibrary(ggplot2)

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

看看这两个数据长什么样，因为数就是图，图就是数！

class(varechem)## [1] "data.frame"class(varespec)## [1] "data.frame"

glimpse(varechem)## Rows: 24## Columns: 14## $ N        <dbl> 19.8, 13.4, 20.2, 20.6, 23.8, 22.8, 26.6, 24.2, 29.8, 28.1, 2…## $ P        <dbl> 42.1, 39.1, 67.7, 60.8, 54.5, 40.9, 36.7, 31.0, 73.5, 40.5, 3…## $ K        <dbl> 139.9, 167.3, 207.1, 233.7, 180.6, 171.4, 171.4, 138.2, 260.0…## $ Ca       <dbl> 519.4, 356.7, 973.3, 834.0, 777.0, 691.8, 738.6, 394.6, 748.6…## $ Mg       <dbl> 90.0, 70.7, 209.1, 127.2, 125.8, 151.4, 94.9, 45.3, 105.3, 11…## $ S        <dbl> 32.3, 35.2, 58.1, 40.7, 39.5, 40.8, 33.8, 27.1, 42.5, 60.2, 3…## $ Al       <dbl> 39.0, 88.1, 138.0, 15.4, 24.2, 104.8, 20.7, 74.2, 17.9, 329.7…## $ Fe       <dbl> 40.9, 39.0, 35.4, 4.4, 3.0, 17.6, 2.5, 9.8, 2.4, 109.9, 4.6, …## $ Mn       <dbl> 58.1, 52.4, 32.1, 132.0, 50.1, 43.6, 77.6, 24.4, 106.6, 61.7,…## $ Zn       <dbl> 4.5, 5.4, 16.8, 10.7, 6.6, 9.1, 7.4, 5.2, 9.3, 9.1, 8.1, 10.2…## $ Mo       <dbl> 0.30, 0.30, 0.80, 0.20, 0.30, 0.40, 0.30, 0.30, 0.30, 0.50, 0…## $ Baresoil <dbl> 43.90, 23.60, 21.20, 18.70, 46.00, 40.50, 23.00, 29.80, 17.60…## $ Humdepth <dbl> 2.2, 2.2, 2.0, 2.9, 3.0, 3.8, 2.8, 2.0, 3.0, 2.2, 2.7, 2.5, 2…## $ pH       <dbl> 2.7, 2.8, 3.0, 2.8, 2.7, 2.7, 2.8, 2.8, 2.8, 2.8, 2.7, 2.9, 2…

glimpse(varespec)## Rows: 24## Columns: 44## $ Callvulg <dbl> 0.55, 0.67, 0.10, 0.00, 0.00, 0.00, 4.73, 4.47, 0.00, 24.13, …## $ Empenigr <dbl> 11.13, 0.17, 1.55, 15.13, 12.68, 8.92, 5.12, 7.33, 1.63, 1.90…## $ Rhodtome <dbl> 0.00, 0.00, 0.00, 2.42, 0.00, 0.00, 1.55, 0.00, 0.35, 0.07, 0…## $ Vaccmyrt <dbl> 0.00, 0.35, 0.00, 5.92, 0.00, 2.42, 6.05, 2.15, 18.27, 0.22, …## $ Vaccviti <dbl> 17.80, 12.13, 13.47, 15.97, 23.73, 10.28, 12.40, 4.33, 7.13, …## $ Pinusylv <dbl> 0.07, 0.12, 0.25, 0.00, 0.03, 0.12, 0.10, 0.10, 0.05, 0.12, 0…## $ Descflex <dbl> 0.00, 0.00, 0.00, 3.70, 0.00, 0.02, 0.78, 0.00, 0.40, 0.00, 0…## $ Betupube <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.02, 0.00, 0.00, 0.00, 0…## $ Vacculig <dbl> 1.60, 0.00, 0.00, 1.12, 0.00, 0.00, 2.00, 0.00, 0.20, 0.00, 0…## $ Diphcomp <dbl> 2.07, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.07, 0…## $ Dicrsp   <dbl> 0.00, 0.33, 23.43, 0.00, 0.00, 0.00, 0.03, 1.02, 0.30, 0.02, …## $ Dicrfusc <dbl> 1.62, 10.92, 0.00, 3.63, 3.42, 0.32, 37.07, 25.80, 0.52, 2.50…## $ Dicrpoly <dbl> 0.00, 0.02, 1.68, 0.00, 0.02, 0.02, 0.00, 0.23, 0.20, 0.00, 0…## $ Hylosple <dbl> 0.00, 0.00, 0.00, 6.70, 0.00, 0.00, 0.00, 0.00, 9.97, 0.00, 0…## $ Pleuschr <dbl> 4.67, 37.75, 32.92, 58.07, 19.42, 21.03, 26.38, 18.98, 70.03,…## $ Polypili <dbl> 0.02, 0.02, 0.00, 0.00, 0.02, 0.02, 0.00, 0.00, 0.00, 0.00, 0…## $ Polyjuni <dbl> 0.13, 0.23, 0.23, 0.00, 2.12, 1.58, 0.00, 0.02, 0.08, 0.02, 0…## $ Polycomm <dbl> 0.00, 0.00, 0.00, 0.13, 0.00, 0.18, 0.00, 0.00, 0.00, 0.00, 0…## $ Pohlnuta <dbl> 0.13, 0.03, 0.32, 0.02, 0.17, 0.07, 0.10, 0.13, 0.07, 0.03, 0…## $ Ptilcili <dbl> 0.12, 0.02, 0.03, 0.08, 1.80, 0.27, 0.03, 0.10, 0.03, 0.25, 0…## $ Barbhatc <dbl> 0.00, 0.00, 0.00, 0.08, 0.02, 0.02, 0.00, 0.00, 0.00, 0.07, 0…## $ Cladarbu <dbl> 21.73, 12.05, 3.58, 1.42, 9.08, 7.23, 6.10, 7.13, 0.17, 23.07…## $ Cladrang <dbl> 21.47, 8.13, 5.52, 7.63, 9.22, 4.95, 3.60, 14.03, 0.87, 23.67…## $ Cladstel <dbl> 3.50, 0.18, 0.07, 2.55, 0.05, 22.08, 0.23, 0.02, 0.00, 11.90,…## $ Cladunci <dbl> 0.30, 2.65, 8.93, 0.15, 0.73, 0.25, 2.38, 0.82, 0.05, 0.95, 2…## $ Cladcocc <dbl> 0.18, 0.13, 0.00, 0.00, 0.08, 0.10, 0.17, 0.15, 0.02, 0.17, 0…## $ Cladcorn <dbl> 0.23, 0.18, 0.20, 0.38, 1.42, 0.25, 0.13, 0.05, 0.03, 0.05, 0…## $ Cladgrac <dbl> 0.25, 0.23, 0.48, 0.12, 0.50, 0.18, 0.18, 0.22, 0.07, 0.23, 0…## $ Cladfimb <dbl> 0.25, 0.25, 0.00, 0.10, 0.17, 0.10, 0.20, 0.22, 0.10, 0.18, 0…## $ Cladcris <dbl> 0.23, 1.23, 0.07, 0.03, 1.78, 0.12, 0.20, 0.17, 0.02, 0.57, 0…## $ Cladchlo <dbl> 0.00, 0.00, 0.10, 0.00, 0.05, 0.05, 0.02, 0.00, 0.00, 0.02, 0…## $ Cladbotr <dbl> 0.00, 0.00, 0.02, 0.02, 0.05, 0.02, 0.00, 0.00, 0.02, 0.07, 0…## $ Cladamau <dbl> 0.08, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0…## $ Cladsp   <dbl> 0.02, 0.00, 0.00, 0.02, 0.00, 0.00, 0.02, 0.02, 0.00, 0.07, 0…## $ Cetreric <dbl> 0.02, 0.15, 0.78, 0.00, 0.00, 0.00, 0.02, 0.18, 0.00, 0.18, 0…## $ Cetrisla <dbl> 0.00, 0.03, 0.12, 0.00, 0.00, 0.00, 0.00, 0.08, 0.02, 0.02, 0…## $ Flavniva <dbl> 0.12, 0.00, 0.00, 0.00, 0.02, 0.02, 0.00, 0.00, 0.00, 0.00, 0…## $ Nepharct <dbl> 0.02, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0…## $ Stersp   <dbl> 0.62, 0.85, 0.03, 0.00, 1.58, 0.28, 0.00, 0.03, 0.02, 0.03, 0…## $ Peltapht <dbl> 0.02, 0.00, 0.00, 0.07, 0.33, 0.00, 0.00, 0.00, 0.00, 0.02, 0…## $ Icmaeric <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.07, 0.00, 0.00, 0…## $ Cladcerv <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0…## $ Claddefo <dbl> 0.25, 1.00, 0.33, 0.15, 1.97, 0.37, 0.15, 0.67, 0.08, 0.47, 1…## $ Cladphyl <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0…

接下来先进行mantel test检验，然后画图即可。

# 计算mantel <- mantel_test(varespec, varechem,                      spec_select = list(Spec01 = 1:7,                                         Spec02 = 8:18,                                         Spec03 = 19:37,                                         Spec04 = 38:44)) %>%   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, breaks = c(-Inf, 0.01, 0.05, Inf), # 对P值进行分割，便于映射颜色                  labels = c("< 0.01", "0.01 - 0.05", ">= 0.05")))## `mantel_test()` using 'bray' dist method for 'spec'.## `mantel_test()` using 'euclidean' dist method for 'env'.

# 画图qcorrplot(correlate(varechem), type = "lower", diag = FALSE) +  geom_square() +  geom_couple(aes(colour = pd, size = rd), # 这行代码是关键              data = mantel,               curvature = nice_curvature()) +

  # 下面就是各种颜色和名称设置  scale_fill_gradientn(colours = RColorBrewer::brewer.pal(11, "RdBu")) +  scale_size_manual(values = c(0.5, 1, 2)) +  scale_colour_manual(values = color_pal(3)) +  guides(size = guide_legend(title = "Mantel's r",                             override.aes = list(colour = "grey35"),                              order = 2),         colour = guide_legend(title = "Mantel's p",                                override.aes = list(size = 3),                                order = 1),         fill = guide_colorbar(title = "Pearson's r", order = 3))

简单方便，快捷画图，出图效果也很好。

当然各种细节都是可以用ggplot2语法修改的，喜欢折腾的可以自己尝试下~

这个图看起来很复杂，但是展示的信息确实很多，不过既然你需要这张图，那我相信你应该理解这张图的意思~

最后，我还是要强调下基础知识的重要性！如果要经常用R，最好还是系统学习下，吃饭的工具怎么能马马虎虎呢？

本文首发于公众号：医学和生信笔记

“

医学和生信笔记，专注R语言在临床医学中的使用，R语言数据分析和可视化。主要分享R语言做医学统计学、meta分析、网络药理学、临床预测模型、机器学习、生物信息学等。