Mantel Test分析与绘图

生信学习小达人

已于 2024-08-20 23:55:04 修改

阅读量319

点赞数 14

分类专栏： R 文章标签： linux

于 2024-08-20 21:48:47 首次发布

本文链接：https://blog.csdn.net/hgz2020/article/details/141369033

版权

R 专栏收录该内容

12 篇文章 0 订阅

订阅专栏

1.前言

2.步骤

3.在R语言中，除了mantel_test函数，还有其他几个工具和方法可以用于进行Mantel Test分析：

4.利用ggcor包在进行Mantel Test分析

5.使用ggcor包进行Mantel Test分析

6.两个距离矩阵的行名和列名不完全相同的处理方法

1.前言

Mantel Test是一种统计方法，用于检验两个矩阵之间的相关性。在生态学研究中，它常用来分析群落距离矩阵（比如基于Bray-Curtis距离）和环境变量距离矩阵（如pH值、温度或地理位置）之间的相关性。如果Mantel Test的相关性系数较大，且p值较小，这通常意味着环境因子对微生物群落的影响较大。

2.步骤

在R语言中，进行Mantel Test分析和绘图可以通过以下步骤完成：

加载必要的R包，例如linkET、tidyverse、RColorBrewer等。
准备数据，通常是两个距离矩阵，它们应该有相同的行数和样本标识。
使用mantel_test函数进行Mantel Test分析，该函数可以计算两个矩阵之间的相关性，并给出r值和p值。
对结果进行处理，比如将r值和p值分为不同的区间，以便于后续可视化。
使用qcorrplot函数绘制相关性热图，并使用geom_square和geom_couple添加Mantel Test的结果到热图上。
根据需要调整图表的美观性，包括颜色、大小、标签等。

例如，下面的R代码展示了如何使用linkET包进行Mantel Test分析并绘制热图和网络图：

library(linkET)
library(tidyverse)
library(RColorBrewer)

# 假设varechem和varespec是两个已经准备好的距离矩阵
data("varechem", package = "vegan")
data("varespec", package = "vegan")

# 进行Mantel Test分析
mantel <- mantel_test(varespec, varechem) %>%
  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")))

# 绘制热图和网络图
correlate(varechem) %>%
  qcorrplot(type = "lower", diag = T) +
  geom_square() +
  geom_couple(aes(colour = pd, size = rd), data = mantel, curvature = 0.1) +
  scale_fill_gradientn(colours = RColorBrewer::brewer.pal(9, "RdBu")) +
  scale_size_manual(values = c(0.5, 1, 2)) +
  scale_colour_manual(values = color_pal(3)) +
  labs(fill = "Pearson's correlation",
       size = "Mantel's r value",
       colour = "Mantel's p value")

另一组R代码示例：

rm(list=ls())#好习惯，确保有干净的 R 环境
# setwd("C:/Users/Desktop/take")
library(linkET)
library(ggplot2)
library(ggtext)
library(dplyr)
library(RColorBrewer)
library(cols4all)
library(tidyverse)

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

#计算环境因子相关性系数：
cor2 <- correlate(varechem)
corr2 <- cor2 %>% as_md_tbl() ##在linkET包中，as_md_tbl()可以将相关性矩阵转换成一个"md"表格，这是linkET特有的数据结构，用于存储和操作相关性数据
write.csv(corr2, file = "pearson_correlate(env&env).csv", row.names = TRUE)

head(corr2)
#mantel test:
mantel <- mantel_test(varespec, varechem,
                      mantel_fun = 'mantel', #支持4种："mantel"使用vegan::mantel()；"mantel.randtest"使用ade4::mantel.randtest()；"mantel.rtest"使用ade4::mantel.rtest()；"mantel.partial"使用vegan::mantel.partial()
                      spec_select = list(spec01= 1:1,
                                         spec02=5:5,
                                         spec03 = 7:7
                      )) #这里分组为随机指定，具体实操需按自己的实际数据分组
head(mantel)
write.csv(mantel, file = "mantel_result(bio&env).csv", row.names = TRUE)

#对mantel的r和P值重新赋值（设置绘图标签）：
mantel2 <- mantel %>%
  mutate(r = cut(r, breaks = c(-Inf, 0.25, 0.5, Inf),
                 labels = c("<0.25", "0.25-0.5", ">=0.5")),
         p = cut(p, breaks = c(-Inf, 0.001, 0.01, 0.05, Inf),
                 labels = c("<0.001", "0.001-0.01", "0.01-0.05", ">= 0.05")))
head(mantel2)
#首先，绘制相关性热图(和上文相同):


##############################
p4 <- qcorrplot(cor2,
                grid_col = "#00468BFF",
                "white","#42B540FF",
                grid_size = 0.2,
                type = "upper",
                diag = FALSE) +
  geom_square() +
  scale_fill_gradientn(colours = c("#00468BFF",
                                   "white","#42B540FF"),
                       limits = c(-1, 1))

# 打印出来看看

#添加显著性标签：
p5 <- p4 +
  geom_mark(size = 4,
            only_mark = T,
            sig_level = c(0.05, 0.01, 0.001),
            sig_thres = 0.05,
            colour = 'white')
p5
#在相关性热图上添加mantel连线：
p6 <- p5 +
  geom_couple(data = mantel2,
              aes(colour = p, size = r),
              curvature = nice_curvature())
p6
#继续美化连线：
p7 <- p6 +
  scale_size_manual(values = c(1, 2, 3)) + #连线粗细
  scale_colour_manual(values = c4a('brewer.set2',4)) + #连线配色
  #修改图例：
  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 = 5),
                               order = 1),
         fill = guide_colorbar(title = "Pearson's r", order = 3))+
  theme(
    text = element_text(size = 16, family = "serif"),
    plot.title = element_text(size = 16, colour = "black", hjust = 0.5),
    legend.title = element_text(color = "black", size = 16),
    legend.text = element_text(color = "black", size = 16),
    axis.text.y = element_text(size = 16, color = "black", vjust = 0.5, hjust = 1, angle = 0),
    axis.text.x = element_text(size = 16, color = "black", vjust = 0.5, hjust = 0.5, angle = 0)
  )

p7

2.1代码解读

如果目的是进行Mantel Test，通常需要为每个分组选择多个样本，以便能够分析组内和/或组间的关系。正确的使用方式可能是：

spec_select = list( group1 = c(1:10), # 假设第一组包含前10个样本

group2 = c(11:20), # 第二组包含接下来的10个样本

group3 = c(21:30) # 第三组包含接下来的10个样本 )

#对mantel的r和P值重新赋值（设置绘图标签）： mantel2 <- mantel %>% mutate(r = cut(r, breaks = c(-Inf, 0.25, 0.5, Inf), labels = c("<0.25", "0.25-0.5", ">=0.5")), p = cut(p, breaks = c(-Inf, 0.001, 0.01, 0.05, Inf), labels = c("<0.001", "0.001-0.01", "0.01-0.05", ">= 0.05")))

这段代码是使用dplyr包中的mutate和cut函数来处理Mantel Test的结果，将连续的r值和p值转换为分区间的标签，以便于在绘图时提供更直观的表示。下面是对代码的详细解释：

mantel2 <- mantel %>%：这行代码使用管道操作符%>%将mantel对象传递给mutate函数，并创建一个新的对象mantel2。
mutate(...)：这个函数用于修改数据框，添加新的列或改变现有列。
r = cut(r, breaks = c(-Inf, 0.25, 0.5, Inf), labels = c("<0.25", "0.25-0.5", ">=0.5"))：这里使用cut函数将相关性系数(r值)分为三个区间：小于0.25、0.25到0.5之间、以及大于等于0.5。每个区间都被赋予了一个标签。
p = cut(p, breaks = c(-Inf, 0.001, 0.01, 0.05, Inf), labels = c("<0.001", "0.001-0.01", "0.01-0.05", ">= 0.05"))：同样使用cut函数将p值分为四个区间：小于0.001、0.001到0.01之间、0.01到0.05之间、以及大于等于0.05。每个区间也被赋予了一个标签。

这种处理方式使得原始的连续变量(r值和p值)转换为分类别变量，便于在绘图时通过颜色或形状的不同来区分不同的区间，使得图形输出更加直观易懂。

例如，在使用ggplot2或其他可视化包进行绘图时，您可以使用这些新的标签作为图例的依据，从而在图中直接展示每个区间的类别，而不是原始的数值。这对于展示Mantel Test结果的统计显著性特别有用

p4 <- qcorrplot(cor2,
grid_col = "#00468BFF",
"white","#42B540FF",
grid_size = 0.2,
type = "upper",
diag = FALSE) +
geom_square() +
scale_fill_gradientn(colours = c("#00468BFF",
"white","#42B540FF"),
limits = c(-1, 1))

这段代码是使用qcorrplot函数和ggplot2语法来创建一个相关性热图（heatmap），并通过添加图层来定制其外观。下面是对代码的详细解释：

p4 <- qcorrplot(cor2, ...)：这行代码使用qcorrplot函数创建一个相关性热图，并将其结果赋值给变量p4。cor2是传入的相关性矩阵。
grid_col = "#00468BFF", "white", "#42B540FF"：这个参数为热图中的网格线设置了颜色。这里指定了三种颜色，但正确的语法应该使用向量形式，如grid_col = c("#00468BFF", "white", "#42B540FF")。
grid_size = 0.2：设置网格线的大小。
type = "upper"：指定只显示相关性矩阵的上三角部分。
diag = FALSE：指定不填充对角线上的值。
geom_square()：使用ggplot2的geom_square图层在热图中添加正方形，表示相关性的大小。
scale_fill_gradientn(colours = c("#00468BFF", "white", "#42B540FF"), limits = c(-1, 1))：这个函数调用设置热图中颜色的渐变，从#00468BFF（一种蓝色）经过白色到#42B540FF（一种绿色），并设置颜色范围从-1到1。

3.在R语言中，除了`mantel_test`函数，还有其他几个工具和方法可以用于进行Mantel Test分析：

vegan包：这个包提供了mantel函数，它可以直接用来计算两个距离矩阵之间的Mantel Test统计量。例如，使用mantel函数可以检验基于物种丰度的距离矩阵与基于环境参数的距离矩阵之间的相关性。
ade4包：提供了多种与Mantel Test相关的函数，例如mantel.randtest和mantel.rtest，这些函数可以用于更复杂的Mantel Test变体，如基于随机化的测试。
ggcor包：但通常在R社区中，ggcor包被推荐用于相关性分析，并且可能包含用于Mantel Test的函数或工具。
linkET包：提供了mantel_test函数的另一个实现，可以用来进行Mantel Test分析并绘制热图和网络图展示Mantel Test的相关性。

在进行Mantel Test分析时，您可以根据具体的分析需求和数据类型选择合适的工具和方法。例如，如果您需要检验群落物种组成是否与环境相关，可以使用vegan包中的mantel函数来计算基于Bray-Curtis距离的物种组成矩阵和基于欧几里得距离的环境参数矩阵之间的相关性。如果您需要更高级的分析，如偏Mantel Test，可以考虑使用ade4包中提供的函数。

4.利用`ggcor`包在进行Mantel Test分析

ggcor包在进行Mantel Test分析时具有一些显著的优势，同时也存在一些局限性：

优势：

可视化效果：ggcor包可以创建信息丰富且美观的相关性热图，这些热图可以直观地展示多个变量之间的相关性。
多种相关性图形：提供多种图形选项，例如geom_square(), geom_circle2(), geom_star()等，这为用户提供了多种数据可视化的方式。
简便的Mantel Test结果整合：ggcor包通过anno_link()函数可以方便地将Mantel Test的结果与相关性热图结合起来，提供了一种直观展示Mantel Test结果的方法。

局限性：

安装困难：一些用户报告了安装ggcor包时遇到的困难，可能需要从GitHub上特定用户的仓库进行安装。
文档和社区支持：相比于一些更流行的R包，ggcor可能拥有较少的文档和社区支持，这可能增加了初学者使用该包的难度。
更新和维护：由于ggcor包在GitHub上的最后更新时间可能较长，可能存在一些功能上的滞后或未解决的bug。

总的来说，ggcor包是一个强大的工具，尤其适用于需要复杂相关性热图和Mantel Test结果可视化的高级分析。然而，用户可能需要面对安装和使用上的一些挑战。

5.使用`ggcor`包进行Mantel Test分析

使用ggcor包进行Mantel Test分析的具体步骤如下：

安装和加载ggcor包：如果ggcor包还未安装，可以使用devtools包从GitHub安装：

install.packages("devtools") devtools::install_github("hannet91/ggcor")

然后加载ggcor包：

library(ggcor)
准备数据：确保你有两个距离矩阵，它们应该有相同的行名和列名，代表相同的样本。
进行Mantel Test：使用mantel_test()函数进行Mantel Test分析。这个函数接受两个距离矩阵作为输入，并返回Mantel Test的结果，包括相关性系数(r)和p值。

mantel_result <- mantel_test(matrix1, matrix2)
处理Mantel Test结果：根据需要，你可能想要将Mantel Test的结果(r值和p值)分为不同的区间，以便于后续可视化。可以使用dplyr包中的mutate()和cut()函数来实现：

library(dplyr) mantel_result <- mantel_result %>% 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")))
绘制相关性热图：使用quickcor()函数计算并绘制其中一个矩阵的变量之间的相关性热图，然后使用geom_square()或其他几何对象添加相关性图形：

correlation_plot <- quickcor(matrix1, type = "lower") + geom_square()
将Mantel Test结果添加到热图：使用anno_link()函数将Mantel Test的结果作为链接添加到相关性热图中，可以根据r值和p值调整链接的颜色和大小：

final_plot <- correlation_plot + anno_link(aes(colour = pd, size = rd), data = mantel_result)
调整和展示最终图形：根据需要调整图形的美观性，包括颜色、大小、标签等，并展示最终的图形：

print(final_plot)

请注意，上述步骤是一个通用的流程，具体的函数调用和参数设置可能需要根据你的具体数据和分析需求进行调整。

6.两个距离矩阵的行名和列名不完全相同的处理方法

如果两个距离矩阵的行名和列名不完全相同，你需要确保它们代表相同的样本，并且具有相同的顺序。以下是一些调整步骤：

核对样本标识：检查两个矩阵的行名和列名，列出所有不匹配的样本标识。
重命名样本：使用rownames()和colnames()函数来重命名矩阵的行和列，确保两个矩阵的行名和列名完全一致。例如：

# 假设matrix1和matrix2是两个距离矩阵 # matrix1的行名需要修改为与matrix2一致 rownames(matrix1) <- matrix2_rownames_to_match <- make.names(names_to_match_with_matrix2)
排序样本：如果行名和列名相同，但顺序不一致，可以使用order()函数对行和列进行排序：

# 根据matrix2的行名排序matrix1的行 matrix1 <- matrix1[order(rownames(matrix1)), ] # 如果需要，也可以根据matrix2的列名排序matrix1的列 matrix1 <- matrix1[, order(colnames(matrix1))]
检查和验证：使用all.equal()函数检查两个矩阵的行名和列名是否完全相同：

all.equal(rownames(matrix1), rownames(matrix2)) all.equal(colnames(matrix1), colnames(matrix2))
删除或合并样本：如果两个矩阵中有不匹配的样本，你可能需要删除这些样本或将它们合并。删除样本可以使用逻辑索引或subset()函数：

# 删除matrix1中不存在于matrix2的样本 matrix1 <- matrix1[intersect(rownames(matrix1), rownames(matrix2)), ]
创建新的样本标识：如果需要，可以为两个矩阵创建一个新的共同的样本标识符列表，并用这个列表更新两个矩阵的行名和列名。
重新计算距离矩阵（如果需要）：如果样本被删除或合并，可能需要重新计算距离矩阵以反映这些变化。
保存和备份：在进行任何修改之前，保存原始数据的备份，并在修改后保存更新后的数据。