Mantel Test
算法原理讲解、代码实现、绘图可视化
Tips / 提示 这篇博文将是你可以在所有中文搜索引擎里找得到的最详细的关于Mantel Test的中文分析讲解。
起源
在一次课题组师兄汇报的时候，我第一听说了Mantel Test，当时第一眼就被这个漂亮的图形所吸引，所以就想着以后也能用到自己的文章里，便自己花时间了解了下。
全网能找到这张图的资料，最早的时间是2015年，在Science中发表的一篇文章里第一次出现了这张图。
SUNAGAWA S, COELHO L P, CHAFFRON S, et al. Structure and function of the global ocean microbiome[J]. Science, 2015, 348(6237): 1261359.
只能说不愧能发Science，文章里面的每一张图都是经过精心处理的。
介绍
1967年，美国国立卫生研究院（National Institutes of Health）的生物统计学家Nathan Mantel提出了Mantel Test分析方法。在统计学中，传统相关系数只能用于计算分析一个数据矩阵中每两列变量之间的相关性，而在面对两个矩阵之间的相关性时就一筹莫展。Mantel Test的出现克服了传统相关系数的这一缺陷，该方法本质上是分析两个矩阵之间的相关性。自Mantel Test推出以来，其在微生物群落生态学、动植物种群遗传学等领域都得到快速发展与应用，在分析微生物菌群与环境因子之间的相关性方面更是得到了广泛的应用。
在使用Mantel Test分析环境因子与微生物群落结构之间的相关性时，通常对微生物群落OTU数据矩阵使用Bray-Curtis相异度（Bray-Curtis dissimilarity）来计算微生物群落结构之间的差异性；而对环境变量数据矩阵采用欧氏距离（Euclidean Distances）进行差异性计算。
Mantel Test的分析过程主要包括：分别使用各自的距离公式计算两个数据矩阵的距离矩阵，然后将两个距离矩阵进行压缩得到两个压缩距离列，然后计算这两列的相关性（一般都采用皮尔逊pearson相关性指数）；在完成一次计算后，对原数据矩阵中的一列或者两列进行置换，重新计算距离公式以及压缩距离公式，计算新的相关性系数（r值）；经过成千上万次的置换后，观察实际数据的r值在经过多次置换后所得的r值分布中的位置，如果跟随机置换得到的结果站队较近，则说明相关性不显著，相反则说明两个矩阵具有较强的显著相关性。
讲解
下面结合我查阅的一些中外文资料，对Mantel Test做出自己的理解，如果有错误的地方，恳请指正交流（文末有微信）。
图形讲解
右侧上三角
首先来看图形右半部分，这部分大家都很常见，是一个相关性热图，它代表了一个数据矩阵中每两列之间的相关性。而计算相关性的算法一般都选择Pearson相关。
皮尔逊（Pearson）相关（r），它测量两个变量（x和y）之间的线性相关性。它也称为参数相关性检验，因为它取决于数据的分布。仅当x和y来自正态分布时才可以使用它。y = f（x）的图称为线性回归曲线。
举个例子：Oxygen与Temperature的相关性为深蓝色，对比最后边色条可以发现，深蓝色对应的相关性为-0.8~-1.0之间，说明其两者呈现出很强的相关性（相关性强弱看绝对值），且呈现负相关（相关性方向看正负）。
左侧下三角
再来看图形的左侧，左侧反映出来的也是相关性（r值），只不过除了相关性之外还反映出了一个显著性（p值），
在相关性分析中，r值代表相关性，p（sig）值代表显著性。p值是判断r值相关系数是否具有统计学意义，判定标准一般为0.05，如果p值小于0.05说明两者之间具有较强的相关性，即r值有意义。然后再看r值，|r|值越大，相关性越好。正数指正相关，即一个变量随着另一个变量增大而增大；负数指负相关，即一个变量随另一个变量的增大而减小。
可以看出右侧上三角中的r值与左侧下三角中的r值是一样的，具有相同的统计学概念。
Pearson相关系数的一点小补充
Pearson主要是把两个序列看成两个变量，计算两个变量是否线性相关。Pearson公式如：
分子是协方差的计算公式，协方差我们知道，是计算两个变量的相关性，大于0是正相关，小于是负相关，相等不相关。除以分母相当于归一化到[-1,1]之间。所以，Pearson相关系数的计算结果也等于将数据矩阵进行标准化后再求协方差，此时求出的协方差就等于源数据矩阵中各列的相关性。
代码演示：
import pandas as pd # 引包
import seaborn as sns
df = sns.load_dataset('iris') # 加载机器学习经典数据集——鸢尾花数据集
df.info() # 查看数据信息
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 150 entries, 0 to 149
Data columns (total 5 columns):
# Column Non-Null Count Dtype
--- ------ -------------- -----
0 sepal_length 150 non-null float64
1 sepal_width 150 non-null float64
2 petal_length 150 non-null float64
3 petal_width 150 non-null float64
4 species 150 non-null object # 这一列为花的分类，为Object文字对象，需要删除
dtypes: float64(4), object(1)
memory usage: 6.0+ KB
df.pop('species') # 删除分类列species文字列
df # 查看数据集
Out[14]:
sepal_length sepal_width petal_length petal_width
0 5.1 3.5 1.4 0.2
1 4.9 3.0 1.4 0.2
.. ... ... ... ...
[150 rows x 4 columns]
X = (df - df.mean()) / df.std(ddof=0) # 将数据矩阵标准化
SIGMA = (X.T @ X) / X.shape[0] # 计算数据矩阵的协方差矩阵
SIGMA
Out[20]:
sepal_length sepal_width petal_length petal_width
sepal_length 1.000000 -0.117570 0.871754 0.817941
sepal_width -0.117570 1.000000 -0.428440 -0.366126
petal_length 0.871754 -0.428440 1.000000 0.962865
petal_width 0.817941 -0.366126 0.962865 1.000000
df.corr() # 计算数据矩阵中每两列之间的Pearson相关性
Out[21]:
sepal_length sepal_width petal_length petal_width
sepal_length 1.000000 -0.117570 0.871754 0.817941
sepal_width -0.117570 1.000000 -0.428440 -0.366126
petal_length 0.871754 -0.428440 1.000000 0.962865
petal_width 0.817941 -0.366126 0.962865 1.000000由上述代码可以发现，一个数据矩阵的Pearson相关性系数矩阵就等于其协方差矩阵。
算法讲解
Mantel Test算法的讲解使用Python语言实现，编程语言只是工具，与算法概念无关，选择Python只是因为这门语言我比较顺手~
数据准备
正如上面的讲解所说，Mantel Test是对两个数据矩阵之间的相关性进行分析。因此，需要先准备好两个数据矩阵。
举个栗子，我想要对一个微生物数据矩阵、一个环境因子数据矩阵进行分析：
微生物OTU矩阵
环境因子矩阵
注意看，上方就是两个进行检验的矩阵，需要注意的是，两个数据矩阵的行索引应该是相互对应的。
计算距离矩阵
数据准备好之后，要进行的第一步操作便是计算两个数据矩阵各自的距离矩阵。一般情况下，对微生物或者基因数据矩阵采用bc距离公式计算，而对于环境因子数据矩阵采用欧氏距离公式计算。
bc距离计算
计算公式：
详细计算过程：
- min(S_{A,i}, S_{B,i}) = 7+3+4+1+0 = 15
- \sum S_{A,i} = 10+8+4+1+1 = 24
- \sum S_{B,i} = 7+3+4+8+4+0 = 22
- D = 1 - 2*15/(24+22) = 0.3478
D值越小表示二者组成差异小。
也可以直接使用Python科学计算库scipy包中的工具进行计算：
from scipy.spatial.distance import braycurtis as bc
import pandas as pd
OTU = pd.DataFrame([['A', 10, 8, 4, 1, 1], ['B', 7, 3, 8, 4, 0]],
columns=['community', 'OTU1', 'OTU2', 'OTU3', 'OTU4', 'OTU5'])
OTU.set_index('community', inplace=True)
OTU
Out[4]:
OTU1 OTU2 OTU3 OTU4 OTU5
community
A 10 8 4 1 1
B 7 3 8 4 0
bc(OTU.iloc[0,:], OTU.iloc[1,:])
Out[5]: 0.34782608695652173欧式距离
欧式距离就不用介绍了，就是初中就学过的那个欧氏距离，直接自己写下就行了。
非冗余距离公式
经过计算后，得到的距离公式应该是这样的：
这个时候会发现，其实这个距离公式是有瑕疵的，为啥这么说呢？可以看到1→2的距离其实是与2→1的距离一样的，且1→1、2→2的距离都是0。也就是说，这个距离公式中有一半的元素（上三角或者下三角）以及对角线上的元素其实都是无效的数据，因为我们根本不需要它或者用不到。
这个时候我们就需要对距离公式进行压缩，使其成为一个非冗余距离矩阵，应该为：
使其上三角或者下三角全为0，这个时候，这个距离矩阵才是非冗余的。
自己封装了一个可以计算欧式距离或者bc距离公式，并可以得到冗余距离矩阵以及非冗余距离矩阵的函数：
def calc_distance(df, method='euc', redundant=True):
"""
计算欧式距离(Euclidean)和BC(Bray-Curtis)距离
:param df: 输入矩阵
:param method: euc：欧氏距离（默认值）；bc：Bray-Curtis距离；Others：Undefined
:param redundant: 是否返回冗余距离矩阵？冗余矩阵就是上三角和下三角一样的一个矩阵
:return: 返回距离公式
"""
method = method.lower()
distance_df = np.zeros((df.shape[0], df.shape[0])) # 创建一个用于存储距离的矩阵
if redundant: # 计算冗余矩阵
for i in range(len(df)): # 遍历矩阵元素
for j in range(len(df)): # 由于距离矩阵中，对角线代表每个元素到自身的距离（必定为0），且上三角和下三角是重复的，因此精简算法，只对下三角进行填充
if method == 'euc': # 欧氏距离
distance_df[i, j] = np.sqrt(np.sum((df.iloc[i] - df.iloc[j]) ** 2))
elif method == 'bc': # 计算BC距离
# distance_df[i, j] = bc(df.iloc[i], df.iloc[j]) # 引包计算
_ = pd.DataFrame([df.iloc[i], df.iloc[j]])
distance_df[i, j] = 1 - 2 * np.sum(_.min()) / np.sum(np.sum(_))
else: # 报错：未定义的距离计算公式！
raise ValueError(f'Undefined distance formula: {method}!')
else: # 计算非冗余矩阵（算法都一样）
for i in range(1, len(df)): # 遍历矩阵元素
for j in range(i): # 由于距离矩阵中，对角线代表每个元素到自身的距离（必定为0），且上三角和下三角是重复的，因此精简算法，只对下三角进行填充
if method == 'euc': # 欧氏距离
distance_df[i, j] = np.sqrt(np.sum((df.iloc[i] - df.iloc[j]) ** 2))
elif method == 'bc': # 计算BC距离
distance_df[i, j] = bc(df.iloc[i], df.iloc[j])
else: # 报错：未定义的距离计算公式！
raise ValueError(f'Undefined distance formula: {method}!')
return pd.DataFrame(distance_df)距离向量
其实，上一步得到的非冗余距离矩阵还是可以再次压缩，使之成为向量的格式：
这样的话，对两个数据矩阵计算得到的距离矩阵经过一系列的操作后，压缩成为了两个距离向量，这个时候，才完成了距离数据的准备。
在上一步的距离矩阵上，可以得到距离向量的函数。
def condensed_distance_matrices(df):
"""
输入一个距离矩阵，对该矩阵进行压缩
:param df: 一个下三角距离矩阵
:return: 返回压缩后的距离数据（Series）
"""
pointer = [] # 存放距离指向
value = [] # 存放距离
for i in range(1, df.shape[0]): # 遍历矩阵元素
for j in range(i):
pointer.append(str(f'{i}➡{j}'))
value.append(df.iloc[i, j])
return pd.DataFrame({'Pointer': pointer, 'Value': value})相关性的计算
Mantel Test所计算的相关性，其本质上是对两个数据矩阵的距离矩阵相关性进行检验。在得到两个距离向量之后，就可以使用Pearson相关系数来反应两个距离向量之间的相关性。
而且我本人经过检验（R、Python），确实是这么一回事儿。
显著性的计算
按照上述方法，我们已经得到了相关性r值。但是，这个r值是否可靠还需要显著性p值来进行检验。而想要得到p值，就先必须理解统计学里的一个假设-检验（hypothesis testing）的概念。
其基本思路就是：
- 我想要证明数据矩阵A与数据矩阵B具有相关性；
- 那么我先假设数据矩阵A与数据矩阵B之间没有相关性，预先设定的检验水准为0.05，当检验假设为真，但被错误地拒绝的概率，记作α，通常取α=0.05或α=0.01；
- 对样本观察值按相应的公式计算出统计量的大小（例如T值等）；
- 根据统计量的大小及其分布确定检验假设成立的可能性P的大小并判断结果。
而这个假设检验体现在Mantel Test里面就是，按照上述得到r值的方法，对两个源数据矩阵的两行或者两列进行置换，再次计算出一个r值。就这样，进行千千万万次的置换，就可以得到千千万万个r值。这个得到的千千万万个r值，理论上是服从正态分布的。
如果我们第一次计算的r值，能够落在接受原假设的区域上，那就说明原假设为真，即两个数据矩阵之间没有相关性；如果r值落在了拒绝原假设的域上，就说明可以拒绝原假设，即两个数据矩阵之间存在相关性。
举个例子：
选定检验水准α为0.05，如果第一次计算得到的r值为0.1516028059，且该r值在符合正态分布的经过千千万万次置换得到的r值list中可以落到拒绝原假设的域上，那就说明两个数据矩阵之间存在相关性。
至于是单边检验还是双边检验，看了下源代码：
@property
def p(self):
if self._p is None:
if self.tail == "upper": # 右侧检验
self._p = sum(self.correlations >= self.r) / self.perms # self.perms就是千千万万次置换得到的r值得list的长度
elif self.tail == "lower": # 左侧检验
self._p = sum(self.correlations <= self.r) / self.perms
else: # 双边检验
self._p = sum(abs(self.correlations) >= abs(self.r)) / self.perms
return self._pR语言计算
安装包
# 安装包
## install.packages("devtools") # 如果没有安装开发者工具箱的话，需要先安装
devtools::install_github("Hy4m/linkET", force = TRUE) 加载包
# 加载包
library(linkET)
library(ggplot2)加载数据集并进行Mantel检验
# Mantel检验结果
data("varechem", package = "vegan")
data("varespec", package = "vegan")
mantel_test(varespec, varechem)
# A tibble: 14 x 4 # 输出结果
spec env r p
* <chr> <chr> <dbl> <dbl>
1 spec N 0.152 0.029
2 spec P 0.221 0.008
3 spec K 0.213 0.009
4 spec Ca 0.0961 0.138
5 spec Mg 0.0634 0.22
6 spec S 0.137 0.042
7 spec Al 0.278 0.004
8 spec Fe 0.167 0.045
9 spec Mn 0.221 0.014
10 spec Zn 0.0400 0.33
11 spec Mo 0.00330 0.48
12 spec Baresoil 0.145 0.045
13 spec Humdepth 0.188 0.023
14 spec pH -0.0282 0.573举例说明：微生物矩阵与环境因子矩阵中N元素相关性的显著性p值为0.029<0.05，说明其显著性较强，即相关性具有统计学意义；r值为0.152，说明两者呈现出相关性，且为正相关。
绘制Mantel Test相关性图
# Mantel绘图
library(dplyr)
library(ggplot2)
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),
labels = c("< 0.01", "0.01 - 0.05", ">= 0.05")))
qcorrplot(correlate(varechem), type = "lower", diag = FALSE) +
geom_square() +
geom_couple(aes(colour = pd, size = rd), data = mantel, curvature = 0.1) +
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))
----- END -----