1/3
2/3
3/3
【R 语言可视化】Mantel_test 的实现与绘图
东瓶西镜未有时
2022年09月20日 01:55

在查阅了相关的资料后，对它进行了深度的分解发现这种图主要由左右两部分组成：

左边部分进行的是 Mantel test 分析右边是相同矩阵不用因子或参数进行的相关性分析，并绘制了热力图。通过编程语言将 Mantel test 分析与热力图进行结合，既展示了同一矩阵内不同因素的相关性，又展示了不同矩阵的相关性。

一、什么是 Mantel test 分析？它能用来干什么？

传统的相关性检验只能就两列数据，或者一组数据进行相关关系的检验。随着对科学问题更加深入的探讨，如何对两组数据进行相关性分析呢？Mantel test 的提出就可以解决这个问题。其次 Mantel tests 是确定两组距离测度矩阵之间的相关性测试方法，目的是判断矩阵 A 中的样本距离与矩阵 B 中的样本距离是否相关。

Mantel test 使用的前提是两个矩阵，矩阵 A 和 B 要有相同的行数，例如，矩阵 A 有 n*120（列 * 行）数据，矩阵 B 有 m*120（列 * 行）数据，此时可以进行 Mantel test。Mantel test 假设两个矩阵之间不存在相关性，如果数据结果中发现 P 值显著，说明存在相关性。并且矩阵 A 和矩阵 B 中样本样本之间的距离具有一定关系（同增同减，即矩阵 A 的样本距离增加，则矩阵 B 对应的样本距离也会增加）。

更加具体和详细的内容可以参考 https://doi.org/10.1111/2041-210x.12018（Dismantling the Mantel tests）

二、具体 Mantel test 分析过程是怎样的？

在上文中简单的介绍了 Mantel test 拟解决的问题与一些需要注意的地方。在这里简单说一下 Mantel test 过程：

  1. 矩阵 A 和矩阵 B 对应展开；

  2. 计算相关系数，在没有具体设置时，默认是 pearson 系数；

  3. 置换其中的一列或者两列数据在计算一个值，在进行 N 次循环后

  4. 观察实际的 R 值和 R 值分布的位置（实际的 R 值跟置换后的 R 值位置接近，则不相关）。

三、绘图之前还需要做什么准备？

1. 理论基本了解，在绘图之前需要对上文中提到的图进行深度的剖析，即绘图流程。

在经过分析后发现，此图会分为这么几部分：

（PS. 将进行 Mantel test 的两组数据分为矩阵 A 和矩阵 B 其中矩阵 B 需要绘制热力图）

对矩阵 B 做相关性分析、绘制矩阵 B 的热力图并保留 1/2。对矩阵 A 和矩阵 B 做 Mantel test 分析然后进行组合。

2. 在了解绘图流程后，还需要安装绘图软件相关的东西。

这里用的绘图软件组合是 R 语言和 R studio 软件以及比较关键的一个工具 R tools。关于 R 语言以及 R 语言以及相关软件的安装可以参考：

可能是最好的R语言安装指南​（Ps. 并不是手懒，是因为视频讲解真的比手写快很多）

在安装好 R 语言和 Rstudio 后，如果有一定其他语言的编程基础最好，加上注释可以很好的编写 R 语言的代码。如果没有 R 语言的基础想更加系统的学习 R 语言可以参考 B 站上的相关教学视频。（Ps. 如果只是想画这个图可以参考我下面的案例）

在这里需要安装必须的几个 R 包：dplyr、linkKEt、ggplot2 包

需要在 Console 输入代码进行安装：

输入：

代码块
R
自动换行
复制代码
install.packages("ggplot2") install.packages("dplyr") devtools::install_github("Hy4m/linKEt"，force = TRUE) 
PackageVersion（"linKEt"）
复制成功

其中 linKEt 包可以再 github 上查询到源代码以及具体使用的一些案例 

https://github.com/Hy4m/linkET（膜拜大佬）

🔶*理论、编程语言、编程软件以及必备的包已经安装完毕下面就可以开始进行绘图。🚩

➡首先还是先把数据进行分类好并储存成两个.csv 文件。切记保存的哪个 csv 文件是用来绘制热力图！！⬅

1. 调用 R 包并读取数据：

代码块
R
自动换行
复制代码
library(dplyr)
library(linkET)
library(ggplot2)
speciese <- read.csv(file.choose(),header = T,row.names = NULL)
env <- read.csv(file.choose(),header = T,row.names = NULL)
复制成功

读取数据主要利用 red.csv 指令。这个指令的利用方式有很多种，为了避免大家在修改路径时因为复制粘贴的原因导致路径会一直报错，所以选择 fil.choose（）的方式，可以利用交互的方式选择事前准备好的 CSV 文件。（特别说明，准备好的 CSV 文件要和正在编程的 R 语言程序在同一个目录下，即同一个文件夹内，如下图）

在执行 speciese <- read.csv (file.choose (),header = T,row.names = NULL) env <- read.csv (file.choose (),header = T,row.names = NULL) 意思是读取选择的 csv 文件并命名为 speciese 和 env

读取后会在 R stuido 右侧显示

2. 进行 Mantel test 分析

代码块
R
自动换行
复制代码
mantel01 <- mantel_test(speciese, env, 
                 spec_select = list(Y1 = 1:7,
                               Y2 = 8:18,     
                               Y3 = 19:37, 
                               Y4 = 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),                 
              labels = c("< 0.01", "0.01 - 0.05", ">= 0.05")))
复制成功

将读取好的 csv 文件，进行 mantel test。讲运行结果进行了进一步的划分。运行结果如下：

相比较上图会出现一个新的文件，点开后会看到具体分析的数据。

在运行此程序时会在 R studio 下方出现两个红色的警告：

     `mantel_test()` using 'bray&#​39; dist method for 'spec&#​39;.      `mantel_test()` using 'euclidean&#​39; dist method for 'env&#​39;.

可以不用管这个，不影响后面的程序运行。

代码块
R
自动换行
复制代码
spec_select = list(Y1 = 1:7,                    
              Y2 = 8:18,                    
              Y3 = 19:37,                    
              Y4 = 38:44))
复制成功

这语句代表的是，将矩阵 A 按照列进行分类，其中 1-7 分为第一类，8-18 分为第二类依次类推。后面可以根据自己的数据进行灵活应用修改，比如说只有四列数据就想分成四类则可以写成：

代码块
R
自动换行
复制代码
spec_select = list(Y1 = 1,                    
              Y2 = 2,                    
              Y3 = 3,                    
              Y4 = 4 ))
复制成功

**Pay Attention **💡在 R 语言中 %>% 是一个管道函数（来自于 dplyr 包）意思是将这一步的结果直接传参给下一步。

代码块
R
自动换行
复制代码
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),                  
              labels = c("< 0.01", "0.01 - 0.05", ">= 0.05")))
复制成功

将 Mantel test 运算结果进行分类，rd = cut (r, breaks = c (-Inf, 0.2, 0.4, Inf),labels = c ("< 0.2&#​34;, "0.2 - 0.4&#​34;, ">= 0.4&#​34;)),// 将 R 值分为 < 0.2,0.2-0.4,>=0.4 三类。Label 是标签的意思。下面 pd 的语句同理。

3. 绘图

代码块
R
自动换行
复制代码
qcorrplot(correlate(env),          
          type = "lower",          
          diag = FALSE,) +  
   geom_square() + 
   geom_couple(aes(colour = pd, size = rd),data = mantel01, curvature = 0.1, 
            node.colour = c("blue", "blue"),             
            node.fill = c("grey", "grey"),             
            node.size = c(3.5, 2.5), ) +  
    scale_fill_gradientn(colours = RColorBrewer::brewer.pal(11, "RdBu"),
                   limits = c(-1, 1),                      
                   breaks = seq(-1,1,0.5)) +  
    scale_size_manual(values = c(0.5, 1, 1.5, 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 = 1.5), 
                                 order = 1),        
           fill = guide_colorbar(title = "Pearson's r", order = 3))
复制成功

绘图主要指令是 qcorrlot 指令，在此段程序中它分为这么几部分：

（1）热力图的绘制

代码块
R
自动换行
复制代码
qcorrplot(correlate(env), 
          type = "lower",
          diag = FALSE,
          ) +
  geom_square()
复制成功

correlate（绘制热力图的矩阵），type=“lower” 保留下三角 ，diag=FALSE 去除对角线的填充。

修改一下：type=“upper” 保留上三角，diag=TRUE。保留对角线填充

geom_square () // 热力图里以方块进行填充。

（2）绘制 Mantel test

代码块
R
自动换行
复制代码
qcorrplot(correlate(env), 
          type = "lower",
          diag = FALSE,
          ) +
  geom_square() +
  geom_couple(aes(colour = pd, size = rd),data = mantel01, curvature = 0.1,
              node.colour = c("blue", "blue"),
              node.fill = c("grey", "grey"),
              node.size = c(3.5, 2.5),
              )
复制成功

解析一下：geom_couple (aes (colour = pd, size = rd),data = mantel01, curvature = 0.1）//colour 是指将上面代码中的 pd 也就是 p 值用连线的颜色表示，size 表示上面代码中的 rd 也就是 R 值用连线的粗细表示，data=mantel01 数据来源为 manltel01 文件。    geom_couple (node.colour = c ("blue&#​34;, "blue&#​34;),node.fill = c ("grey&#​34;, "grey&#​34;), node.size = c (3.5, 2.5),) // 用来设置 Y1-Y4 N 到 pH 连接点的设置 它并不是成对出现的。

上图中橘色的的设置为 Y1-Y4 的连接点设置，blue 外圈颜色，grey 内部颜色，3.5 是大小。蓝色部分同理设置的是 N-pH 连接点。

至此其实已经完成了 Mantel test 和热力图的组合图。

但是还需要在进行亿点点的优化（头秃.jpg）。

因此在此基础上，再加上如下代码：

代码块
R
自动换行
复制代码
scale_fill_gradientn(colours = RColorBrewer::brewer.pal(11, "RdBu"),
                       limits = c(-1, 1),
                       breaks = seq(-1,1,0.5)) +
  scale_size_manual(values = c(0.5, 1, 1.5, 2)) +
  scale_colour_manual(values = color_pal(3))
复制成功

scale_fill_gradientn (colours = RColorBrewer::brewer.pal (11, "RdBu&#​34;),limits = c (-1, 1), breaks = seq (-1,1,0.5))// 用于设置热力图的标签，颜色设置，图例的设置。

在这个基础上在加上以下代码：

代码块
R
自动换行
复制代码
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 = 1.5), 
                               order = 1),
         fill = guide_colorbar(title = "Pearson's r", order = 3))
复制成功

这段代码可以实现对图例标签的改变：

（3）保存图片

代码块
R
自动换行
复制代码
ggsave("Mantel test.png",width = 8,height = 6)
复制成功

可以自定义需要保存图片的名字。格式以及图片的长和宽。

DONE！

至此已经完成了 Mantel test 和绘图。

（PS.  R studio 可以选择执行部分代码，需要按住鼠标左键选择一部分代码然后点击 RUN 进行执行，如下）大佬可自动省略该提醒

全部代码如下：

代码块
R
自动换行
复制代码
library(dplyr)
library(linkET)
library(ggplot2)
speciese <- read.csv(file.choose(),header = T,row.names = NULL)
env <- read.csv(file.choose(),header = T,row.names = NULL)

mantel01 <- mantel_test(speciese, env,
                        spec_select = list(Y1 = 1:7,
                                           Y2 = 8:18,
                                           Y3 = 19:37,
                                           Y4 = 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),
                  labels = c("< 0.01", "0.01 - 0.05", ">= 0.05")))
qcorrplot(correlate(env), 
          type = "lower",
          diag = FALSE,
          ) +
  geom_square() +
  geom_couple(aes(colour = pd, size = rd),data = mantel01, curvature = 0.1,
              node.colour = c("blue", "blue"),
              node.fill = c("grey", "grey"),
              node.size = c(3.5, 2.5),
              ) +
  scale_fill_gradientn(colours = RColorBrewer::brewer.pal(11, "RdBu"),
                       limits = c(-1, 1),
                       breaks = seq(-1,1,0.5)) +
  scale_size_manual(values = c(0.5, 1, 1.5, 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 = 1.5), 
                               order = 1),
         fill = guide_colorbar(title = "Pearson's r", order = 3))
ggsave("Mantel test.png",width = 8,height = 6)
复制成功

最终效果图及文件夹变化图：

写在最后：

目前正在进行进一步的调整与开发，关于基于 R 语言可视化的 mantel test 和绘图已经很成熟。我也只是进行了初步的调用和调试。后期还会在此基础上修改或增加其他代码，尽量以更加简单的方式进行绘图作业。感谢每个点开这篇文章的大佬们。我们下篇博客见！

（ps.文章中的代码和测试数据集已经上传到github上，包括后面对代码的一些维护也在github中进行。https://github.com/DongPingXiJin/Mantel-test）

可视化 R语言 mantel test
cv18682729