mantel test可视化，别再只知道ggcor！

医学和生信笔记

发布于 2022-11-15 10:43:49

7860

发布于 2022-11-15 10:43:49

文章被收录于专栏：医学和生信笔记

在几年前出现了一个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, union
library(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))