三万字懂得微生物群落差异分析（下） / 开普饭

9.3 用ANOSIM假设检验组间的不同

9.3.1相似性分析介绍(ANOSIM)

ANOSIM仅仅是基于两个距离矩阵之间标准化秩相关的Mantel测试的修正版本。

ANOSIM检验是由Clarke(1993)开发的一种无分布的多元数据分析方法，经常被群落生态学家和微生物组研究人员使用。它是一种非参数程序，用于检验两组或多组样本之间无差异的假设，基于组间和组内相似性的排列检验(Clarke 1993)。它比较了物种(或任何其他分类群)丰度和组成在采样单元(如:就分组因素或实验处理水平而言，与ANOVA分析类似，ANOSIM将组成员资格或处理水平作为因素，并将其建模作为解释变量。相似度的分析基于一个简单的想法:如果测试的组是有意义的，那么组内的样本应该比不同组的样本在组成上更相似。因此，零假设是:在测试组或治疗条件的成员之间没有差异。本实验采用Bray-Curtis相似性测度法。

anosim检验统计量是基于组间和组内的平均秩差。如下:

里面R是检验统计量，是相对群内差异的一个指标。

M=N(N-1)/2 是样品对数。

N是为样本总数(受试者)。

rB是组间相似性排序的平均值

rW是组内相似性排序的平均值。

执行ANOSIM有五个主要步骤:

步骤1:计算距离矩阵。

步骤2:计算秩差值，将最小的不相似值赋1。

步骤三:计算组间和组内秩差的平均值。

步骤4:用上式计算检验统计量R

R被解释为一个相关系数，作为其他类型的相关系数，如Pearson’s相关系数，是“效应大小”的衡量标准。检验统计量是为了检验在原假设下组间是否存在差异。如果零假设是正确的，那么R=0，这意味着平均而言，组间和组内的差异是相同的。当高相似性和低相似性完全混合且与群体没有关系时，就会发生这种情况。如果零假设被拒绝，那么R不等于0，这表明组内所有的样本对比任何来自不同组的样本对更相似。例如，在这种情况下，所有最相似的样本都在样本组内，则R = 1。理论上，R < 0也是有可能的。但实际上这种情况是在生态和微生物组研究中不太可能。极端情况R =-1，表示最相似的样本都不在组内。

步骤五4:检验显著性。

正如在PERMANOVA检验中，检验统计量R的p值是通过排列获得:将样本观察值随机分配给组。然后，将组内和组间的排序相似度与随机生成的相似度进行比较。显著性检验简单地说就是R大于R的观测值的排列R的分数。它是给定R大于或等于这个值的概率。假设检验背后的算法与PERMANOVA相同:如果两组抽样单位在其物种(或其他分类群)组成上确实不同，那么组间的组成差异应该大于组内的组成差异。

9.3.2 用vegan包说明ANOSIM的相似性分析

相似性分析(ANOSIM)提供了一种统计检验两组或多组抽样单位之间是否存在显著差异的方法。相似性分析是通过函数anosim()在vegan包中实现的。该函数假设所有组内的等级差异的中位数和范围都是相等的。输入数据是一个不同矩阵。它可以由函数dist()或vegdist()生成。该功能还具有总结和绘图方法来进行后期建模分析。下面是一个语法示例。

anosim (data, grouping, permutations = 1000, distance = ”bray”)

其中，data为数据矩阵，其中行为样本，列为响应变量。或不同的物体或对称的不同方阵;; grouping = grouping variable (一个因素); permutations = 评估ANOSIM统计量的显著性的排列数; distance = 距离或不相似度量。如果输入数据不具有不相似结构或为对称方阵，则需要指定距离。这里用Vdr小鼠粪便数据来说明ANOSIM测试。首先，我们需要加载数据和尚未加载的vegan包。如前所述，将获得分组信息。

> abund_table=read.csv("VdrFecalGenusCounts.csv",row.names=1, check.names=FALSE) > abund_table<-t(abund_table) > grouping<-data.frame(row.names=rownames(abund_table),t(as.data.frame (strsplit(rownames(abund_table),"_")))) > grouping$Group <- with(grouping,ifelse(as.factor(X2)%in% c (11,12,13,14,15),c("Vdr-/-"), c("WT"))) > grouping<- grouping[,c(4)]

用Bray-Curtis的不同点来匹配ANOSIM

以下R代码使用Bray-Curtis不同矩阵作为输入数据运行ANOSIM。

> library(vegan) > bray<-vegdist(abund_table, "bray") > anosim(bray, grouping,permutations = 1000) Call: anosim(dat = bray, grouping = grouping, permutations = 1000) Dissimilarity: bray ANOSIM statistic R: 0.19 Significance: 0.2 Permutation: free Number of permutations: 1000

从函数anosim()的输出中可以看出，它包括了建模公式以及函数调用中使用的不同方法。由排列值和排列值个数的显著性可知，0.2的p值大于0.05，表明在0.05显著水平下，组内相似度不大于组间相似度。我们可以得出结论，没有证据表明组内样本比随机概率预期的更相似

下面的R代码使用丰度数据帧作为输入数据运行ANOSIM。

> anosim(abund_table, grouping, permutations = 1000, distance = "bray") distance = "bray") Dissimilarity: bray ANOSIM statistic R: 0.19 Significance: 0.19 Permutation: free Number of permutations: 1000

ANOSIM的匹配结果可以由summary().函数概括。

> fit <- anosim(bray, grouping,permutations = 1000) > summary(fit) Call: anosim(dat = bray, grouping = grouping, permutations = 1000) Dissimilarity: bray ANOSIM statistic R: 0.19 Significance: 0.18 Permutation: free Number of permutations: 1000 Upper quantiles of permutations (null model): 90% 95% 97.5% 99% 0.282 0.323 0.538 0.846 Dissimilarity ranks between and within classes: 0% 25% 50% 75% 100% N Between 2 9.5 17 21.50 28 15 Vdr-/- 1 6.5 14 19.75 25 10 WT 5 8.0 11 16.50 22 3

最后我们可以绘出结果（图9.4）:

> plot(fit)

使用Jaccard不同配合ANOSIM. 下面的R代码使用Jaccard方法适合ANOSIM

> anosim(abund_table, grouping, permutations = 1000, distance = "jaccard") Call: anosim(dat = abund_table, grouping = grouping, permutations = 1000, distance= "jaccard") Dissimilarity: jaccard ANOSIM statistic R: 0.19 Significance: 0.2 Permutation: free Number of permutations: 1000

使用Sørensen不同配合ANOSIM. 以下使用Sørensen方法适用于ANOSIM rensen

> fit_S <- anosim(abund_table, grouping, permutations = 1000,distance = "Sørensen") > summary(fit_S) Call: anosim(dat = abund_table, grouping = grouping, permutations = 1000) Dissimilarity: bray ANOSIM statistic R: 0.19 Significance: 0.19 Permutation: free Number of permutations: 1000 Upper quantiles of permutations (null model): 90% 95% 97.5% 99% 0.282 0.344 0.546 0.846 Dissimilarity ranks between and within classes: 0% 25% 50% 75% 100% N Between 2 9.5 17 21.50 28 15 Vdr-/- 1 6.5 14 19.75 25 10 WT 5 8.0 11 16.50 22 3

9.4 Multi-response Permutation Procedures的假设检验

9.4.1 介绍MRPP

多响应置换程序(MRPP)是一种非参数程序，用于检验两组或多组样本之间无差异的假设，基于组间和组内差异的置换检验(Mielke 1984，1991)。测试差异可能会使均值(位置)或组内距离差异(spread) (Warton et al. 2012)。

与PERMANOVA相似，在概念和方法上，MRPP与方差分析(ANOVA)相结合:它比较组内和组间的差异。潜在的原理也是一样的。MRPP检验统计量是基于组间和组内差异的加权均值之差。它如下所示：

式中，delta= MRPP检验统计量，为抽样单位间两两不相似的组内均数的总体加权均值:Ci=ni/N；N=项目总数: n;=组i中项目的数量。

实施MRPP主要有五个步骤: 第一步：一般情况下使用欧几里得距离计算距离矩阵第二步：计算各组di的平均距离。第三步：计算g多组的增量第四步：确定效果大小。统计量A是对效应大小的衡量，它被解释为ANOISIM中的R;它是作为组内同质性与随机期望的比较而得到的。

检验统计量是检验在零假设下组间无差异。如果组内异质性随机等于期望，则A=0:当组内所有项各都相同时，则A= 1。A < 0.1在生态学中很常见;> 0.3在生态学中是相当高的。然而，据我们所知，目前尚无微生物组效应大小的报道标准。

第五步:检验显著性。正如在PERMANOVA和ANOSIM检验中，检验统计量delta的p值是通过Monte Carlo置换得到的。显著性检验是简单地用小于观测到的增量的置换增量的比例，并进行小样本校正。它是在这个或更小的条件下的概率。

9.4.2 用vegan包做MRPP

MRPP提供了一种统计检验两组或两组以上采样单位之间是否存在显著差异的方法。MRPP是通过vegan包中的mrpp()函数实现的。如果输入数据不一致，则可以直接使用。如果它是一个由响应数据结构观察到的矩阵，那么在使用之前需要由函数vegdist()计算不相似性。默认的距离是欧几里得方法，但其他不同的vegdist()也被应用。该函数有总结和绘图方法来进行后期建模分析。下面是一个语法示例:

mrpp(data, grouping, permutations = 1000, distance = “bray”)

其中，data =数据矩阵或数据框架，其中行是观测值，列是结果，或一个不相似对象或对称的不相似方阵。结果可能是单变量的或多变量的;grouping =分组变量(一个因素);permutations =评估MRPP统计量的显著性的置换数:distance = distance or dissimilarity measure如果输入数据没有dissimilarity structure或者是一个对称的方阵，那么需要指定距离。在ANOSIM测试中，我们需要访问数据并加载vegan包。

用Bray-Curtis的不同点来匹配MRPP以下R代码使用Bray-Curtis不同矩阵作为输入数据运行MRPP

> mrpp(bray, grouping,permutations = 1000) Call: mrpp(dat = bray, grouping = grouping, permutations = 1000) Dissimilarity index: bray Weights for groups: n Class means and counts: Vdr-/- WT delta 0.44 0.427 n 5 3 Chance corrected within-group agreement A: 0.0502 Based on observed delta 0.4353 and expected delta 0.4583 Significance of delta: 0.17 Permutation: free Number of permutations: 1000

观测到的和预期的deltas分别为0.4353和0.4583。delta的显著性为0.13，组内一致性A校正的概率为0.0502。由于样本量较小，我们认为在0.05的显著水平下，两个属之间的差异不具有统计学意义。

以下代码使用丰度数据帧作为输入数据运行MRPP:

> mrpp(abund_table, grouping, permutations = 1000, distance = "bray") Call: mrpp(dat = abund_table, grouping = grouping, permutations = 1000, distance= "bray") Dissimilarity index: bray Weights for groups: n Class means and counts: Vdr-/- WT delta 0.44 0.427 n 5 3 Chance corrected within-group agreement A: 0.0502 Based on observed delta 0.4353 and expected delta 0.4583 Significance of delta: 0.15 Permutation: free Number of permutations: 1000

用meandist()函数获得平均距离矩阵函数meandist()计算群内(块)差异(对角线)和簇间(块)差异(非对角线元素)的均值矩阵，以及分组计数的属性n。

> meandist(bray, grouping,permutations = 1000) Vdr-/- WT Vdr-/- 0.4401 0.4766 WT 0.4766 0.4273 attr(,"class") [1] "meandist" "matrix" attr(,"n") grouping Vdr-/- WT 5 3

使用meandist()和 summary()函数获得平均距离和汇总统计我们可以使用meandist()summary()函数来查找组内、组间和总体方法以及他们的差异，以及所有重要的MRPP统计数。选项和分类强度。

> bray_mrpp <- meandist(bray, grouping,permutations = 1000,distance = "bray",weight.type = 1) > summary(bray_mrpp) Mean distances: Average within groups 0.4372 between groups 0.4766 overall 0.4583 Summary statistics: Statistic MRPP A weights n 0.05016 MRPP A weights n-1 0.04900 MRPP A weights n(n-1) 0.04612 Classification strength 0.04131

9.5使用GUniFrac包比较微生物群落

9.5.1 介绍UniFrac，加权UniFrac和广义UniFrac距离矩阵

在多变量微生物组组成分析中，加权和未加权的 UniFrac distance是测量两种微生物组群体的的最常用方法。UniFrac距离，也被称为未加权UniFrac距离，由Lozupone等人在2005年提出(Lozupone and Knight 2005)。它是根据系统发育距离估计微生物组样品之间差异的一种方法。UniFrac距离度量的目标是使来自不同条件的微生物组样本之间的客观比较成为可能。unweighted UniFrac定义如下:

式中，d = unweighted UniFrac distance;A,B=微生物组群落A与B;n=系统发育树的根状分支;bi=i枝长度,piA和piB=群落A和B从分支i下降的类群比例。I()是指示物种是否存在于i支的指示函数。在上式中，群落A和群落B从i支下降的类群比例分别定义为>0和I(piA>0)和I(piB>0)。

按照这种定义，UniFrac度量是通过将两个样本之间不共享的分支长度除以其中一个样本(但不是两个样本都包含的分支长度)来计算的(Lozupone和Knight 2005)。距离为0表示两个样本完全相同，距离为1表示两个样本没有共同类群。unweighted UniFrac distance作为缺席二值检验，只考虑类群的存在和缺席信息;它是检测稀有血统丰度变化的最有效方法。然而，这种距离度量的定义完全忽略了类群丰度信息(Chen et al .2012).2007年，Lozupone等人在原始的未加权方法(Lozupone et al. 2007)的基础上增加了比例加权，因此将这种新的UniFrac方法称为加权UniFrac。加权UniFrac距离利用类群丰度信息，并对丰度差异的分支长度加权;在加权的UniFrad距离度量中，系统发育树的每个分支长度的加权是由两个样本之间的类群丰度比例的差异，而不是只看类群的存在或不存在。加权UniFrac距离(Lozupone et al. 2007)定义为

式中，dW=(标准化)加权UniFrac距离;A、B=微生物组群落A、B;n =系统发育树的根状分支;bi =分支i的长度。

该定义的公式有以下几个特点:(1)以丰度差异为分支长度加权;==(2)即使将丰度数据转化为存在缺失数据，也不能将dW简化为du;(3) dw采用绝对比例差，导致dw值主要由大比例分支决定，对小比例分支丰度变化不太敏感(Chen et al. 2012)。通过向UniFrac距离添加一个比例加权，加权的UniFrac距离减少了低丰度分类群被表示为0或根据采样深度的低计数的问题。在加权的UniFrac中，低丰度分类群的权重更低，因此对度量报告的总距离的影响更小(Wong et al. 2016)。因此，加权的UniFrac可以检测每个谱系中存在的序列数量和分类群的变化(Lozupone et al. 2007);它对检测丰富基因的变化最为敏感(Chen et al. 2012)。然而，未加权或加权的UniFrac距离在检测中等丰富谱系中的变化时可能都不是很有效(Chen et al. 2012)，因为它们对罕见谱系或最丰富的谱系都赋予了过多的权重。因此，Chen等人提出了以下广义UniFrac距离，以统一加权的UniFrac距离和非加权的UniFrac距离。

式中,d^(α)=广义UniFrac距离；α∈[0,1] 是用来控制来自高丰度分支的贡献;A,B=分别是微生物群落A,B;n=生根的系统发育树的分支；bi=分支i的长度；∑ⁿbi(pi^A+pi^B)^(α)=标准因子以确保d^(α)值在[0,1].广义UniFrac 距离包含一个额外的参数α控制家系丰度的权重，所以距离不由高丰度家系决定。α=0.5已经远超于最大值。UniFrac距离的标准化公式有几个特征：（1）利用相对差值|pi^A-pi^B|/（pi^A+pi^B）对大比例分支的权重进行了衰减，权重值在[0,1]区间（2）根据这个公式，如果我们将丰度数据转换为存在/缺失数据，现在dW可以简化为dU，当我们令式（9.11）中α=1，d^(α)减少到dW。

当α=0，我们得到

基于UniFrac距离的伪f统计量定义为

其中tr(.)是矩阵的迹函数，H=X(X^TX)^-1X^T 为设计矩阵X的hat(投影)矩阵，G为高尔居中矩阵，n为样本个数，m为预测因子个数。设dij为群落i与j之间的广义UniFrac距离，且表示为A=(aij)=(-1/2dij²).高尔矩阵定义为G=(I-11/n)A(I-11/n)。未加权、加权 UniFrac距离和广义UniFrac距离的详细的检验统计量及统一公式，读者可以参考作者的原始论文。

9.5.2 母乳数据集

母乳数据集来自两项最近发表的研究(Urbaniak et al.)。2016;Wong等人2016年)。它从哺乳期收集了58个微生物样本加拿大白人女性。母乳是发育中的婴儿的重要细菌来源，并影响新生儿的细菌组成，进而影响到以后生活中的疾病风险。在第一篇论文中，作者使用这组数据比较了早产和出生的、剖腹产(选择性和非选择性)、阴道分娩以及男婴和女婴之间的细菌分布。在第二篇文章中，使用相同的数据集来说明作者自己开发的工具，包括unweighted UniFrac、weighted UniFrac、information UniFrac和ratio UniFrac。在这里，我们使用这个数据集来比较GUniFrac软件包中的微生物群落。

9.5.3

使用GUniFrac包比较微生物群落

GUniFrac包被发展以完成广义的uniFrac距离来比较生物群落，正如PERMANOVA程序， GUniFrac包基于置换估计伪F的统计值的显著性，如果树在一个环境中所占的比例大于偶然的预期值，那么这两个群落就被认为是不同的，在这个包中，加权和不加权UniFrac，以及方差调整加权UniFrac距离也可以实现， UniFrac距离可以被计算出来，并且在GUniFrac包中比较微生物群落，包括d^W,d^(0.5),d⁽⁰⁾,d^(U),以及方差以调整加权UniFrac距离d^VAW，以上部分,我们介绍了d^W,d⁽⁰⁾和d^(U)的公式，d^(0.5)代表一系列距离的中间值，通过以下公式计算出来：

方差调整加权UniFrac距离d^{VAW^{计算公式：}}

其中mi为第i支两个社团的总个体数/reads, m为总个体数/reads。它的开发是为了考虑这样一个事实，即加权UniFrac距离不考虑随机抽样下权重的变化，从而导致较少的功率检测社区之间的差异(Chang et al. 2011)。GUniFrac包依赖于vegan和ape包，作者还建议使用ade4包。一个用法如下所示:

GUniFrac(otu.tab, tree, alpha = c(0, 0.5, 1))

UniFrac测量计算需要两条信息:呕吐表和系统发育树。因此，运行guunifrac包需要两个数据集:一个计数表，如OTU表和系统发育树。接下来，我们将逐步演示如何使用GUniFrac软件包来比较微生物群落。

第一步:加载并读取合适格式的OTU表设置目录的方法与我们在前面章节中所做的相同。“td_OTU_tag_mapped_lineage.txt”数据集包括两条信息:最后一列的OTU计数表和分类。来自原始数据集作者的部分R代码是有用的。我们在这里修改它们，并解释和评论代码的作用。

> setwd("E:/Home/MicrobiomeStatUsingR/Analysis/") > otu_tab <- read.table("td_OTU_tag_mapped_lineage. txt", header=T, sep="\t", row.names=1, comment. char="", check.names=FALSE)

在参数中，row.names = 1 参数指定表中包含行名的列1;comment.char = “ “用于完全关闭注释的解释。它比默认的read.table()快得多。数据集的前6行“td_OTU_tag_mapped_lineage.txt”的内容如下:

> head(otu_tab) S31 S1 S42 S13_T2 S30 S50 S43 S20 S29 S47U S26 S13_T3 S33 S8L 0 38 36 30 14 13 18 27 38 49 7 5251 12 17 35 1 2866 15069 42985 3292 1223 1056 3959 3021 7023 1856 7993 2571 13306 1101 2 437 9831 4628 4231 4473 718 3843 441 4311 2496 5672 4940 4112 4830 3 3356 7407 5355 62 121 14 3616 6184 334 108 31 30 13 56 4 12 478 16 342 23 8 18 97 390 126 4 271 53 7 5 145 238 234 3 56 1 109 156 374 28 0 45 0 1 …… taxonomy 0 Bacteria;Proteobacteria;Gammaproteobacteria;Pasteurellales;Pasteurellaceae;Pasteurella;|93 1 Bacteria;Firmicutes;Bacilli;Bacillales;Staphylococcaceae;Staphylococcus;|77 2 Bacteria;Proteobacteria;Gammaproteobacteria;Pseudomonadales;Pseudomonadaceae;Pseudomonas;|92 3 Bacteria;Proteobacteria;Gammaproteobacteria;Enterobacteriales;Enterobacteriaceae;EscherichiaShigella;|98 4 Bacteria;Proteobacteria;Gammaproteobacteria;Enterobacteriales;Enterobacteriaceae;Klebsiella;|75 5 Bacteria;Actinobacteria;Actinobacteria;Corynebacteriales;Corynebacteriaceae;Corynebacterium;|88

然而，为了正确使用guunifrac包，OTU表应该是一个数字矩阵。下面的R代码用于从原始数据集中删除分类列。

> taxonomy <- otu_tab$taxonomy > otu_tab <- otu_tab[-length(colnames(otu_tab))]

otu表需要通过otu矩阵转置为样本。

> otu_tab <- t(as.matrix(otu_tab)) > head(otu_tab) ...... 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 S31 38 2866 437 3356 12 145 43 0 20 19 33 27 16 21 79 40 30 12 S1 36 15069 9831 7407 478 238 1368 25 306 2821 91 169 323 475 1290 3795 532 186 S42 30 42985 4628 5355 16 234 14491 2 357 153 48 109 136 520 781 617 727 867 S13_T2 14 3292 4231 62 342 3 282 1 188 51 64 168 161 574 557 33 408 7 S30 13 1223 4473 121 23 56 3 0 123 3 13 3 279 207 703 123 617 5 S50 18 1056 718 14 8 1 3 1 20 5 0 1 35 49 119 28 83 5

第二步：使用抽平将样品OTU计数标准化到一个标准测序深度UniFrac的作者推荐使用稀疏数据(Carcer等人，2011年)，并将其用于计算未加权的UniFrac距离矩阵，而非稀疏数据则使用自定义的UniFrac脚本用于加权、信息和比例UniFrac (Wong等人，2016年)。然而，在GUuniFrac软件包手册中，所有的UniFrac距离计算都使用了稀有数据。因为GUuniFrac对不同的测序深度很敏感，为了在相同的基础上比较微生物群，可以使用稀疏法(Chen 2012)。因此，我们遵循包手册;在执行未加权，加权UniFrac和方差调整加权UniFrac距离前细化样本。可以使用vegan package中的rrarefy()或rarefy()函数来实现。

> library(vegan) > otu_tab_rarefy <- rrarefy(otu_tab, min(apply(otu_tab,1,sum)))

步骤3:读取系统发育树GUniFrac包需要一个根树作为输入数据。我们可以使用phangorn包中的函数midpoint()来获取根树。

> library(phangorn) > otu_tree <- midpoint(otu_tree) > otu_tree Phylogenetic tree with 115 tips and 114 internal nodes. Tip labels: 1, 11, 6686, 18, 230, 82, ... Node labels: NA, 71.9, 79.1, 86.9, 87.1, 52.4, ... Rooted; includes branch lengths.

第四步:计算UniFracs现在，可以使用GUuniFrac软件包计算UniFracs

> library(GUniFrac) > #Calculate the UniFracs > unifracs <- GUniFrac(otu_tab_rarefy, otu_tree, alpha=c(0, 0.5, 1)) $unifracs > dw <- unifracs[,, "d_1"] # Weighted UniFrac > du <- unifracs[,, "d_UW"] # Unweighted UniFrac > dv <- unifracs[,, "d_VAW"]# Variance adjusted weighted UniFrac > d0 <- unifracs[,, "d_0"] # GUniFrac with alpha 0 > d5 <- unifracs[,, "d_0.5"]# GUniFrac with alpha 0.5

第五步:进行PERMANOVA比较UniFrac的一项测量为了检验UniFrac测量的假设，需要从元数据中提取群组信息。我们使用read.table()函数读取元数据，如下所示:

> meta_tab<- read.table("metadata.txt", header=T, sep="\t", row.names=1, comment.char="", check.names=FALSE)

匹配OTU表和分组数据，以保留只出现在两个数据集中的样本.

> otu_meta_matched <- match(rownames(meta_tab),rownames(otu_tab)) > otu_meta_matched <- otu_meta_matched[!is.na(otu_meta_matched)] > otu_tab <- otu_tab[otu_meta_matched,] > meta_tab_Ordered <- meta_tab[match(rownames(otu_tab),rownames (meta_tab)),]

下面的R代码使用了vegan package中的adonis()函数来执行 PERMANOVA比较妊娠距离系列的中间距离。感兴趣的读者可以尝试UniFrac的其他距离指标。

> set.seed(123) > adonis(as.dist(d5) ~ meta_tab$Gestation) Call: adonis(formula = as.dist(d5) ~ meta_tab$Gestation) Permutation: free Number of permutations: 999 Terms added sequentially (first to last) Df SumsOfSqs MeanSqs F.Model R2 Pr(>F) meta_tab$Gestation 4 0.4336 0.108391 1.2001 0.08305 0.22 Residuals 53 4.7869 0.090319 0.91695 Total 57 5.2205 1.00000

类似地，性别和奶类型的影响可以被测试如下。

> adonis(as.dist(d5) ~ meta_tab$Gender) Call: adonis(formula = as.dist(d5) ~ meta_tab$Gender) Permutation: free Number of permutations: 999 Terms added sequentially (first to last) Df SumsOfSqs MeanSqs F.Model R2 Pr(>F) meta_tab$Gender 2 0.1864 0.093197 1.0182 0.0357 0.381 Residuals 55 5.0341 0.091529 0.9643 Total 57 5.2205 1.0000 > adonis(as.dist(d5) ~ meta_tab$milktype) Call: adonis(formula = as.dist(d5) ~ meta_tab$milktype) Permutation: free Number of permutations: 999 Terms added sequentially (first to last) Df SumsOfSqs MeanSqs F.Model R2 Pr(>F) meta_tab$milktype 3 0.4019 0.133976 1.5014 0.07699 0.06 . Residuals 54 4.8186 0.089233 0.92301 Total 57 5.2205 1.00000 --- Signif. codes: 0 '***’ 0.001 '**’ 0.01 '*’ 0.05 '.’ 0.1 ' ’ 1

第六步:使用PermanovaG()函数进行PERMANOVA对多种UniFrac测量进行比较

结合多个距离矩阵可以增加假设检验的能力。下面的R代码使用d(0)，d(0.5),d(1)和函数PermanovaG()进行排列多元方差分析:

> PermanovaG(unifracs[, , c("d_0", "d_0.5", "d_1")] ~ meta_tab$Gestation) $aov.tab F.Model p.value meta_tab$Gestation 1.262764 0.263

9.6 概括和讨论

在本章中，我们介绍了多元群落微生物组数据的假设检验及其在R中的逐步实现。我们用来说明的数据来自我们自己的研究(Jin et al. 2015;王等。或公开提供。读者可以使用本章提供的R代码和解释来分析他们自己的微生物组数据。

beta多样性的多元群落分析有两种方法：生态文献中的“原始数据法”与“距离(Mantel)”方法，PERMANOVA，冗余分析(RDA)，基于距离的RDA(db-RDA)，主坐标的规范化分析属于基于原始数据的方法，而Mantel测试是典型的基于距离的方法。ANOSIM是一个基于两个距离矩阵之间标准化秩相关的修正版本的Mantel测试。因此，它也属于基于距离的方法。Mantel试验的应用领域也适用于相似度分析(ANOSIM)。

问题是哪一种方法更适合分析多元群落微生物组数据?在生态文学中，关于这个问题的讨论很多。(McArdle and Anderson 2001; Legendre et al. 2005; Legendre2007; Laliberte 2008; Legendre et al. 2008; Pélissier et al. 2008; Tuomisto and Ruokolainen 2008; Anderson et al. 2011; Legendre and Legendre 2012)

基于距离的方法存在以下问题:(1)它们不能正确划分数据中的变化，不能提供正确的第一类错误率，因此不适合分析beta多样性;(2) 正如 Anderson 等人认为，更有问题的是，使用划分方法对驱动beta多样性模式的底层过程的相对重要性进行直接推断(Anderson等，2011年)。例如，ANOSIM和Mantel测试已被证明不适合测试关于元数据变化的假设检验，因此，建议在分析beta多样性时，应限于分析beta多样性的变化，而不是分析beta多样性。基于原始数据的替代方法被认为为分析beta多样性提供了更合适和更强大的工具。然而，基于距离的方法确实有一些优点，例如，Mantel测试被认为是一种有效的方法来关联两个距离矩阵(Anderson et al.2011)并提供了更大的灵活性，允许使用其他类型的距离函数，如Jaccard或Steinhaus/Bray-Curtis (Legendre et al. 2005)，它可能适合于关于群点(样本)间beta多样性变化的假设 (Legendre et al. 2005)。与简单的Mantel试验相比，生态学研究也指出，部分Mantel方法的解释在一段时间内存在问题.(e.g., Legendre et al. 2005; Anderson et al. 2011).

关于vegan包中的函数的使用，这个包的作者推荐adonis()而不是mrpp()和anosim()。原因在于函数adonis()允许用连续和/或分类预测器解释beta多样性中的方差进行类似方差的测试。然而，mrpp()和anosim()都只处理分类预测因子，它们不如adonis()强大。(Oksanen et al. 2016).

到目前为止，关于多元群落分析的讨论大多局限于生态学等相关研究领域，而在微生物组领域则没有。然而，微生物组学研究已经采用了这些方法和方法，这些方法的优点和局限性也适用于微生物组学研究。

根际互作生物学研究室简介

根际互作生物学研究室是沈其荣教授土壤微生物与有机肥团队下的一个关注于根际互作的研究小组。本小组由袁军副教授带领，主要关注：1.植物和微生物互作在抗病过程中的作用；2 环境微生物大数据整合研究；3 环境代谢组及其与微生物过程研究体系开发和应用。团队在过去三年中在 isme J， Microbiome， PCE，SBB，Horticulture Research等期刊上发表了多篇文章。欢迎关注微生信生物公众号对本研究小组进行了解。

团队工作及其成果（点击查看）

了解交流合作

团队成员邮箱袁军：unyuan@njau.edu.cn；文涛：2018203048@njau.edu.cn
团队公众号：微生信生物添加主编微信，或者后台留言；

加主编微信加入群聊

关于微生信生物你想要的都在这里

微生信生物

References Anderson, M.J. 2001. A new method for non-parametric multivariate analysis of variance. Austral Ecology 26: 32–46. Anderson, M.J., T.O. Crist, et al. 2011. Navigating the multiple meanings of beta diversity: A roadmap for the practicing ecologist. Ecology Letters 14 (1): 19–28. Benjamini, Y., and Y. Hochberg. 1995. Controlling the false discovery rate: A practical and powerful approach to multiple testing. Journal of the Royal Statistical Society Series B 57: 289–300. Benjamini, Y., and D. Yekutieli. 2001. The control of the false discovery rate in multiple testing under dependency. Annals of Statistics 29: 1165–1188. Bonferroni, C.E. 1936. Teoria statistica delle classi e calcolo delle probabilitàby. Pubblicazioni del R Istituto Superiore di Scienze Economiche e Commerciali di Firenze 8: 3–62 Key: citeulike:1778138. Carcer, D.A., S.E. Denman, et al. 2011. Evaluation of subsampling-based normalization strategies for tagged high-throughput sequencing data sets from gut microbiomes. Applied and Environment Microbiology 77 (24): 8795–8798. Chang, Q., Y. Luan, et al. 2011. Variance adjusted weighted UniFrac: A powerful beta diversity measure for comparing communities based on phylogeny. BMC Bioinformatics 12 (1): 118. Chen, J. 2012. GUniFrac: Generalized UniFrac distances. R package version 1.0. https://CRAN.Rproject.org/package=GUniFrac. Chen, J., K. Bittinger, et al. 2012. Associating microbiome composition with environmental covariates using generalized UniFrac distances. Bioinformatics 28 (16): 2106–2113. Clarke, K.R. 1993. Non-parametric multivariate analysis of changes in community structure. Australian Journal of Ecology 18: 117–143. Hochberg, Y. 1988. A sharper Bonferroni procedure for multiple tests of significance. Biometrika 75: 800–803. Holm, S. 1979. A simple sequentially rejective multiple test procedure. Scandinavian Journal of Statistics 6: 65–70. Hommel, G. 1988. A stagewise rejective multiple test procedure based on a modified Bonferroni test. Biometrika 75: 383–386. Jin, D., S. Wu, et al. 2015. Lack of Vitamin D receptor causes dysbiosis and changes the functions of the murine intestinal microbiome. Clinical Therapeutics 37 (5): 996–1009, e1007. Laliberte, E. 2008. Analyzing or explaining beta diversity? comment. Ecology 89 (11): 3232– 3237. Legendre, P. 2007. Studying beta diversity: Ecological variation partitioning by multiple regression and canonical analysis. Journal of Plant Ecology 1 (1): 3–8. Legendre, P., and L. Legendre. 2012. Numerical ecology. Amsterdam: Elsevier Science BV. Legendre, P., and M.J. Fortin. 1989. Spatial pattern and ecological analysis. Vegetatio 80 (2): 107–138. Legendre, P., D. Borcard, et al. 2005. Analyzing beta diversity: Partitioning the spatial variation of community composition data. Ecological Monographs 75 (4): 435–450. Legendre, P., D. Borcard, et al. 2008. Analyzing or explaining beta diversity? comment. Ecology 89 (11): 3238–3244. 9.6 Summary and Discussion 329 Linnenbrink, M., J. Wang, et al. 2013. The role of biogeography in shaping diversity of the intestinal microbiota in house mice. Molecular Ecology 22 (7): 1904–1916. Lozupone, C., and R. Knight. 2005. UniFrac: A new phylogenetic method for comparing microbial communities. Applied and Environmental Microbiology 71 (12): 8228–8235. Lozupone, C.A., M. Hamady, et al. 2007. Quantitative and qualitative beta diversity measures lead to different insights into factors that structure microbial communities. Applied and Environment Microbiology 73 (5): 1576–1585. Mantel, N. 1967. The detection of disease clustering and a generalized regression approach. Cancer Research 27: 209–220. Mantel, N., and R.S. Valand. 1970. A technique of nonparametric multivariate analysis. Biometrics 26 (3): 547–558. McArdle, B.H., and M.J. Anderson. 2001. Fitting multivariate models to community data: A comment on distance-based redundancy analysis. Ecology 82 (1): 290–297. Mielke, P.W. 1984. Meteorological applications of permutation techniques based on distance functions. In Handbook of statistics, vol. 4, ed. P.R. Krishnaiah and P.K. Sen, 813–830. Amsterdam, North-Holland: Elsevier Science Publishers. Mielke Jr., P.W. 1991. The application of multivariate permutation methods based on distance functions in the earth sciences. Earth-Science Reviews 31: 55–71. Oksanen, J., F. Guillaume Blanchet, et al. 2016. Vegan: Community ecology package. R package version 2.4-1. http://CRAN.R-project.org/package=vegan. Pélissier, R., P. Couteron, et al. 2008. Analyzing or explaining beta diversity? comment. Ecology 89 (11): 3227–3232. Smouse, P.E., J.C. Long, et al. 1986. Multiple regression and correlation extensions of the mantel test of matrix correspondence. Systematic Zoology 35 (4): 627–632. Tuomisto, H., and K. Ruokolainen. 2006. Analyzing or explaining beta diversity? Understanding the targets of different methods of analysis. Ecology 87: 2697–2708. Tuomisto, H., and K. Ruokolainen. 2008. Analyzing or explaining beta diversity: Reply. Ecology 89: 3244–3256. Urbaniak, C., M. Angelini, et al. 2016. Human milk microbiota profiles in relation to birthing method, gestation and infant gender. Microbiome 4 (1): 1.

三万字懂得微生物群落差异分析（下）