普鲁克分析（Procrustes Analysis）评估物种-环境/功能关联度的一个示例

R包vegan的普鲁克分析（Procrustes Analysis）示例

几天前，同学咨询这篇文献中的分析方法，先来看一下（部分截图）。

（Zhao et al, 2019）

大体意思是，作者测量了活性污泥（AS）中的细菌群落（16S）组成和抗生素抗性基因（ARGs）丰度，然后分别执行了OTUs和ARGs 亚型的非度量多维标度（NMDS）分析，基于获得的两组NMDS坐标执行Procrustes分析以表征各样本的细菌物种丰度组成和抗生素抗性基因丰度组成的潜在一致性。

作者通过Procrustes分析的M²个统计量以及基于999次置换检验获得的显著性p值，指出各样本的细菌物种丰度组成和抗生素抗性基因丰度组成存在很强的关联。

大致搞懂了什么是Procrustes分析之后，本篇作个简介。包括其原理，M²统计量的意义，以及怎样计算（在R里可以使用那个函数执行）。

Procrustes分析简介

Procrustes分析（Procrustes Analysis，普鲁克分析）是一种通过分析形状分布，比较两组数据一致性的方法。数学上来讲，就是不断迭代，寻找标准形状（canonical shape），并利用最小二乘法寻找每个对象形状到这个标准形状的仿射变化方式（Gower, 1975）。该过程也称为最小二乘正交映射（least-squares orthogonal mapping）。

Procrustes分析的过程

简单地说，Procrustes分析基于匹配两个数据集中的对应点（坐标），通过平移、旋转和缩放其中一个数据集中点的坐标以匹配另一数据集中对应点的坐标，并最小化点坐标之间的偏差平方和（表示为M²）。对应点坐标之间的偏差称为矢量残差（vector residuals），小的矢量残差代表了两数据集具有更高的一致性。

A，显示了投影至二维空间中的两数据集，两数据集中样本对应（A-a，B-b，C-c）；

B，平移坐标使质心对齐，此时两数据集具有一个公共质心；

C-D，通过不断迭代旋转并缩放调整数据，使M²最小化。

注：坐标平移、旋转和缩放过程中，对所有的样本点执行统一的操作，因此尽管点的位置发生改变，但不会改变数据集的投影“形状”。

Procrustes分析M²统计量的显著性检验（PROTEST）

Procrustes分析对两个数据集中点的坐标进行描述性总结和图形化比较，尽管M²统计量提供了对一致性的度量（M²越小表示两数据集关联度越高），但是并未评估M²是否比预期的要好（M²是否显著，而不是由偶然所致）。

可通过置换检验的原理实现对M²显著性的检验，称为PROTEST或PROcrustean randomization test。通过在其中一个数据集中随机地对观测值进行置换，同时保持数据集的共变结构，之后使用Procrustes分析计算随机置换后数据集与另一数据集的M²值（称为M_p²）。如此随机置换共N次，并记录原始观测值的M² < M_p²的次数，由此获得p值：

如果p达到显著性水平（如p<0.05），则可以认为原始观测值的M²并非由偶然因素所致，M²显著小于随机数据集的M_p²，两个数据集的期望值比在随机情况下表现出更大的一致性。

Procrustes分析在生物学数据分析中的应用

这是一个在面部几何中广泛应用的算法，例如用于人脸识别。同学们有兴趣的话搜一下了解就可以了，本篇扯人脸识别就偏题了（况且白鱼同学自己就是出了名的脸盲，莫要跟我提人脸识别……）。

生物学数据分析中，Procrustes分析也经常出现。

例如在组学分析中，经常需要分析来自相同样本不同组学数据集之间是否存在相似性关系，就微生物组而言，例如开篇时提到的那个示例，物种组成丰度和ARGs丰度是否存在潜在一致性。当然分析潜在关系的方法有很多，Procrustes分析就是一种很好的选择。

再如群落分析中，常见通过Procrustes分析物种组成与环境的关系，以及物种形态、遗传组成、空间结构、行为特征等更多类型的数据集类似的案例。

考虑到白鱼同学的所学专业，所以下文多以“群落分析”举例描述。

Procrustes分析与Mantel Test的比较

提到群落分析，Procrustes分析与Mantel test都是用于分析物种组成和环境属性关系的常见方法，当然二者的具体关注点还是有区别的，方法各有自身的优点。但如果只聚焦在评估两数据集一致性上，似乎Procrustes分析更直观一些。

M²统计量及其显著性检验p值提供了两个数据集之间一致性的总体度量，同时数据集的图形匹配和相关残差提供了比Mantel test更丰富的信息源。在对应点的坐标匹配度较好时，两个数据集表现出良好的一致性。坐标匹配度越差表明这些点与整体趋势不匹配，这类似于回归分析中残差较大的点，这些点不符合样本的总体趋势。此外，PROTEST的统计功效也被证明优于Mantel test的统计功效（Peres-Neto and Jackson, 2001）。因此，如果两组数据之间存在潜在关系，则Procrustes分析更有能力检测到，并且鉴于结果的图形性质，它还提供了出色的解释性准则。

R包vegan的Procrustes分析示例

Procrustes分析在群落分析中得到广泛应用，因此vegan包（这是一个在群落分析中非常知名的R包）中就提供了Procrustes分析的方法，下文就以vegan包的Procrustes分析为例展示。

其它可用于Procrustes分析的R包，请自行了解。

下文示例数据、R代码等的获取链接（提取码，ij18）：

https://pan.baidu.com/s/1VXqsxicZE4a2VdY6bnuJLQ

示例数据

举一个非常常见的例子，在A、B、C三个环境中（每个地点观测了4个样地）测量了物种组成和当地的环境指标，并期望通过Procrustes分析建立环境和物种的关联。

网盘附件示例数据集，一个物种丰度矩阵（spe_table.txt），一个环境变量矩阵（env_table.txt）。

PCA降维及探索性分析

由于两个数据集的属性不同，直接比较并不合适。因此，分别对两个数据集执行主成分分析（PCA）降维，并提取特征轴的坐标（代表了变量集的线性组合）用于比较。

使用vegan的rda()执行PCA。过程中，将环境变量标准化为均值为0和标准差为1，以消除不同环境测度的量纲差异；对物种变量执行Hellinger转化，避免直接应用欧几里得距离时产生的“马蹄形效应”。（为什么这样做，参考前文主成分分析；当然，对于物种数据的降维，基于Bray-curtis距离的主坐标分析（PCoA）也是一个好的选择，本篇未用它举例，大家可以自行尝试下）

library(vegan)
 
##样方-环境属性矩阵
env <- read.delim('env_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
 
#环境变量的 PCA 需要标准化，详情 ?rda
env_pca <- rda(env, scale = TRUE)
 
##样方-物种丰度矩阵
otu <- read.delim('spe_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
 
#物种数据 Hellinger 预转化，详情 ?decostand
otu_hel <- decostand(otu, method = 'hellinger')
 
#对转化后的物种数据执行 PCA，无需标准化，详情 ?rda
otu_pca <- rda(otu_hel, scale = FALSE)
 
##排序图比较，以 PCA 的 I 型标尺为例
par(mfrow = c(1, 2))
biplot(env_pca, choices = c(1, 2), scaling = 1, 
    main = '环境组成的PCA', col = c('red', 'blue'))
biplot(otu_pca, choices = c(1, 2), scaling = 1, 
    main = '物种组成的PCA', col = c('red', 'blue'))

同时也可借助PCA作为探索性分析，看两个数据集中样方的排序坐标是否具有一致性。比如说，在排序图中不同分组都区分的很明显，表明环境与物种存在协同性；如果PCA结果中两个数据集比较后找不到规律，可能表明它们之间的关联度较低。

在示例的两个数据集各自的PCA排序图中，样方的排布还是非常明显的，A、B、C三组的样本都能够分开。尽管具体的坐标有区别，但这无妨，因为Procrustes分析可以通过平移、旋转和缩放坐标，最小化最组坐标之间的偏差平方和（M²），描述二者的关联程度。

Procrustes分析度量两数据集一致性

接下来，将上述两个PCA结果中的样方坐标提取出来，作为Procrustes分析的输入。vegan中，procrustes()执行Procrustes分析。

关于procrustes()中参数X和Y的指定，环境数据的PCA坐标指定于X，物种组成的PCA坐标指定于Y，因为后续Procrustes分析中旋转和缩放操作是对Y而言的，将Y匹配于X。这种顺序可以表示物种依赖于环境的逻辑，代表了将物种匹配于环境，最好不要反过来。

注：procrustes()参数中，当symmetric=FALSE时，将旋转并按比例缩放Y以匹配X，可知这种模式是“非对称”的，X和Y的分配值调换后，Procrustes分析的偏差平方和（M²）也会随之改变。而当symmetric=TRUE时，首先将两组坐标都按比例缩放为单位方差，从而给出更独立于比例的对称统计。这种“对称”模式下，X和Y的分配值调换后，Procrustes分析的偏差平方和（M²）不会发生改变，但注意旋转仍将是非对称的（旋转的仍然是Y）。

#Procrustes 分析
#提取两个 PCA 中的样方排序坐标，均以 I 型标尺为例
site_env <- summary(otu_pca, scaling = 1)$site
site_otu <- summary(otu_pca, scaling = 1)$site
 
#执行 Procrustes 分析，详情 ?procrustes
#以对称分析为例（symmetric = TRUE）
proc <- procrustes(X = env_pca, Y = otu_pca, symmetric = TRUE)
summary(proc)
 
#旋转图
plot(proc, kind = 1, type = 'text')
 
#一些重要的结果提取
names(proc)
 
head(proc$Yrot)  #Procrustes 分析后 Y 的坐标
head(proc$X)  #Procrustes 分析后 X 的坐标
proc$ss  #偏差平方和 M2 统计量
proc$rotation  #通过该值可获得旋转轴的坐标位置

通过平移、旋转和缩放坐标，图中映射在主正交轴中的点是来自环境变量PCA的样方点，映射在斜向正交轴中的点是来自物种组成PCA的样方点，箭头指示了二者中配对的样方。从图形结果中可以看出，环境和物种的潜在关系表现出较好的一致性。

最小化最组坐标之间的平方差之和（M²）为0.2178。

#残差图
plot(proc, kind = 2)
residuals(proc)  #残差值

该图显示了每个配对样方的残差，从底部到顶部的水平线是残差的25％（虚线），50％（实线）和75％（虚线）分位数。如果样方中环境和物种组成更为一致，则旋转图中两个匹配的点距离越近，残差较小，反之环境和物种更无规律，则残差更大。

M²统计量的显著性检验（PROTEST）

上述图形结果显示，环境和物种的潜在关系表现出较好的一致性，同时获得了指示这种一致性的M²统计量。

但是这种一致性是否是显著的，procrustes()没有为我们提供出来。因此，通过另一函数protest()执行置换检验评估Procrustes分析结果的显著性（即PROTEST，见上文概述部分）。

#PROTEST 检验，详情 ?protest
#以 999 次置换为例
#注：protest() 中执行的是对称 Procrustes 分析，X 和 Y 的分配调换不影响 M2 统计量的计算
set.seed(123)
prot <- protest(X = env_pca, Y = otu_pca, permutations = how(nperm = 999))
prot

#重要统计量的提取
names(prot)
prot$signif #p 值
prot$ss #偏差平方和 M2 统计量

999次置换检验后显示p<0.001，结果是非常显著的。

备注：我尝试在procrustes()中添加参数symmetric = TRUE/FALSE，但都返回了对称Procrustes分析的M²统计量，该函数只能执行对称模式的Procrustes分析？不过倒是不影响结果判断，因为对称Procrustes分析显著的话，非对称Procrustes分析肯定也是显著的，因此p值怎样都适用。

ggplot2的一个作图示例

Procrustes分析到这里就大致演示完了。

如果觉得vegan默认的图不方便调整，不妨将上述计算好的坐标导出，通过ggplot2作图，如下展示一个示例。

library(ggplot2)

#提取 Procrustes 分析的坐标
Y <- cbind(data.frame(proc$Yrot), data.frame(proc$X))
X <- data.frame(proc$rotation)

#添加分组信息
group <- read.delim('group.txt', sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
Y$samples <- rownames(Y)
Y <- merge(Y, group, by = 'samples')

#ggplot2 作图
p <- ggplot(Y) +
geom_point(aes(X1, X2, color = groups), size = 1.5, shape = 16) +
geom_point(aes(PC1, PC2, color = groups), size = 1.5, shape = 1) +
scale_color_manual(values = c('red2', 'purple2', 'green3'), limits = c('A', 'B', 'C')) +
geom_segment(aes(x = X1, y = X2, xend = PC1, yend = PC2), arrow = arrow(length = unit(0.1, 'cm')),
color = 'blue', size = 0.3) +
theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'),
legend.key = element_rect(fill = 'transparent')) +
labs(x = 'Dimension 1', y = 'Dimension 2', color = '') +
geom_vline(xintercept = 0, color = 'gray', linetype = 2, size = 0.3) +
geom_hline(yintercept = 0, color = 'gray', linetype = 2, size = 0.3) +
geom_abline(intercept = 0, slope = X[1,2]/X[1,1], size = 0.3) +
geom_abline(intercept = 0, slope = X[2,2]/X[2,1], size = 0.3) +
annotate('text', label = sprintf('M^2 == 0.2178'),
x = -0.21, y = 0.42, size = 3, parse = TRUE) +
annotate('text', label = 'P < 0.001',
x = -0.21, y = 0.38, size = 3, parse = TRUE)

p

#输出图片
ggsave('procrustes.pdf', p, width = 6, height = 5)
ggsave('procrustes.png', p, width = 6, height = 5)

该图和上文中vegan包的默认出图结构是一致的，解读方式参考上文即可。额外添加了样方分组、统计量标识等信息。

参考资料

概念：http://jackson.eeb.utoronto.ca/procrustes-analysis/

R操作：http://john-quensen.com/tutorials/procrustes-analysis/

Gower J C. Generalized procrustes analysis. Psychometrika, 1975, 40, 33-51.

Peres-Neto P R, Jackson D A. How well do multivariate data sets match? The advantages of a Procrustean superimposition approach over the Mantel test. Oecologia, 2001, 129(2):169-178.

Zhao R, Feng J, Liu J, et al. Deciphering of microbial community and antibiotic resistance genes in activated sludge reactors under high selective pressure of different antibiotics. Water Research, 2019, 388-402.

猜你喜欢

系列教程：微生物组入门 Biostar 微生物组  宏基因组

专业技能：学术图表 高分文章 生信宝典 不可或缺的人

一文读懂：宏基因组 寄生虫益处 进化树

必备技能：提问 搜索  Endnote

文献阅读 热心肠 SemanticScholar Geenmedical

扩增子分析：图表解读 分析流程 统计绘图

在线工具：16S预测培养基 生信绘图

科研经验：云笔记  云协作 公众号

编程模板: Shell  R Perl

生物科普:  肠道细菌 人体上的生命 生命大跃进  细胞暗战 人体奥秘  

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外5000+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。PI请明示身份，另有海内外微生物相关PI群供大佬合作交流。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍未解决群内讨论，问题不私聊，帮助同行。

学习16S扩增子、宏基因组科研思路和分析实战，关注“宏基因组”

点击阅读原文，跳转最新文章目录阅读