第七章 微生物组数据的探索性分析.2
Redundancy Analysis (RDA)：与三种基本排序方法相对应，vegan包有三种版本的约束排序：冗余分析(RDA，与主成分分析相关)、主坐标约束分析(CAP，与度量尺度相关)和约束对应分析(CCA，与对应分析相关)。函数rda()中的RDA基于欧几里得距离，并结合了多元回归和PCA。在RDA中，每个规范轴都是所有解释变量的线性组合。正则轴的数量对应于解释变量的数量，或者更准确地说，对应于自由度的数量。您可以使用两种不同的语法通过函数rda()从vegan语法：Matrix和公式语法计算RDA。矩阵语法最简单：只需列出用逗号分隔的对象名称，如下所示。RDA = rda(Y, X, W),其中Y是响应矩阵(类群组成)，X是解释矩阵(环境因子)，W是协变量的可选矩阵。公式语法如下所示：
RDA= rda（Y~var 1+factor A+var 2 * var 3+condition（var 4）;data= XW）。公式的左侧Y是响应矩阵(类群组成)；约束或解释变量在公式的右侧，包括定量变量var 1、分类变量因子A、var 2和var 3之间的相互作用项以及可以分离出来的协变量var 4。解释变量和协变量在Object XW数据集中，必须是数据框。假设检验在PCA中不是这样的。然而，对于两个数据集Y和X，在RDA中，我们可以检验它们之间没有线性关系的零假设。我们使用Smoker数据集来说明RDA的用法。首先，从GUniFrac包加载数据，并使用函数select()from dplyr package创建解释变量的子集。
> library(GUniFrac)
> data(throat.otu.tab)
> data(throat.meta)
> library(dplyr)
> throat_meta <- select(throat.meta, SmokingStatus, Age, Sex, PackYears)
出于RDA的目的，我们使用Hellinger变换来变换分类群数据：
> abund_hell <- decostand (throat.otu.tab, 'hell')
如果提供了环境变量矩阵，则由函数rda()计算RDA()(如果未提供，则将像我们在第7.4.1节中所做的那样计算PCA)。
> rda_hell<- rda(abund_hell * ., throat_meta)
> summary (rda_hell)
我们在这里部分输出
Call:
rda(formula = abund_hell ~ SmokingStatus + Age + Sex + PackYears, data =
throat_meta)
Partitioning of variance:
Inertia Proportion
Total 0.46271 1.0000
Constrained 0.04823 0.1042
Unconstrained 0.41448 0.8958
Eigenvalues, and their contribution to the variance
Importance of components:
RDA1 RDA2 RDA3 RDA4 PC1 PC2 PC3
Eigenvalue 0.02393 0.01330 0.005969 0.005039 0.0683 0.05946 0.03554
Proportion Explained 0.05172 0.02874 0.012900 0.010890 0.1476 0.12849 0.07681
Cumulative Proportion0.05172 0.08045 0.093350 0.104240 0.2518 0.38034 0.45715
Accumulated constrained eigenvalues
Importance of components:
RDA1 RDA2 RDA3 RDA4
Eigenvalue 0.02393 0.0133 0.005969 0.005039
Proportion Explained 0.49612 0.2757 0.123740 0.104460
Cumulative Proportion 0.49612 0.7718 0.895540 1.000000
Scaling 2 for species and site scores
* Species are scaled proportional to eigenvalues
* Sites are unscaled: weighted dispersion equal on all dimensions
* General scaling constant of scores: 2.285815
Species scores
RDA1 RDA2 RDA3 RDA4 PC1 PC2
4695 1.921e-03 -1.129e-02 -7.088e-03 -2.362e-03 -7.858e-04 -5.495e-04
2983 2.385e-03 -1.850e-03 -4.702e-03 4.402e-03 8.961e-04 -3.710e-03
(…)
Site scores (weighted sums of species scores)
RDA1 RDA2 RDA3 RDA4 PC1 PC2
ESC_1.1_OPL -0.444360 -0.479468 -0.17464 -0.351952 0.081583 -0.135377
ESC_1.3_OPL 0.601039 -0.226340 -0.10553 -0.275939 0.183698 -0.077295
(…)
Site constraints (linear combinations of constraining variables)
RDA1 RDA2 RDA3 RDA4 PC1 PC2
ESC_1.1_OPL -0.30617 -0.220346 0.030881 -0.328152 0.081583 -0.135377
ESC_1.3_OPL 0.32421 0.055762 -0.132139 0.107388 0.183698 -0.077295
(…)
Biplot scores for constraining variables
RDA1 RDA2 RDA3 RDA4 PC1 PC2
SmokingStatusSmoker 0.9064 -0.3399 -0.1354 0.21095 0 0
Age 0.2275 -0.0442 0.7482 0.62172 0 0
SexMale 0.4961 0.8360 0.1026 0.21097 0 0
PackYears 0.6576 -0.1515 0.7376 0.02219 0 0
Centroids for factor constraints
RDA1 RDA2 RDA3 RDA4 PC1 PC2
SmokingStatusNonSmoker -0.2502 0.09382 0.03738 -0.05823 0 0
SmokingStatusSmoker 0.2860 -0.10723 -0.04273 0.06655 0 0
SexFemale -0.1995 -0.33619 -0.04127 -0.08484 0 0
SexMale 0.1074 0.18103 0.02222 0.04568 0 0
四个受约束的排序轴(命名为RDA1到RDA4)与4个环境变量(SmokingStatus、Age、Sex、PackYears)相关。不受约束的轴被命名为PC轴。将总方差分为约束方差和非约束方差。约束方差由约束轴(即环境变量)解释；而非约束方差由非约束轴解释(即不受环境因素解释的方差)。约束方差是Y矩阵由解释变量解释的方差的量。它相当于多元回归中有偏的、未调整的R2。方差分割表显示，这4个环境因子对总方差的解释比例为4.8%(0.04823)。函数coef()检索每个规范轴上每个解释变量的规范系数(相当于回归系数)。
> coef(rda_hell)
RDA1 RDA2 RDA3 RDA4
SmokingStatusSmoker 0.177040 -0.120402 -0.187134 0.15184
Age -0.004452 -0.002161 0.002997 0.01757
SexMale 0.092005 0.267873 -0.027579 0.01856
PackYears 0.005257 -0.000520 0.011633 -0.01596
我们可以使用函数RsquareAdj()从排序结果中提取R2和校正后的R2的值。
> RsquareAdj(rda_hell)
$r.squared
[1] 0.1042
$adj.r.squared
[1] 0.03909
该函数返回两个R2：一个是普通R2(r.squared)，另一个是调整后的R2 adj(adj.r.squared)。请注意，R2调整总是低于R2，并且差异随着解释变量数量的增加而增大。R2调整可以是负的，这意味着解释变量比相同数量的随机生成变量解释的变化要小。对残余轴应用Kaiser-Guttman准则。
> rda_hell$CA$eig[rda_hell$CA$eig
> mean(rda_hell$CA$eig)]
PC1 PC2 PC3 PC4 PC5 PC6 PC7
0.068299 0.059456 0.035542 0.024040 0.017406 0.015601 0.014898
PC8 PC9 PC10 PC11 PC12 PC13
0.012182 0.010314 0.009654 0.009245 0.008760 0.007745
接下来，我们绘制RDA的结果。根据您感兴趣的排序，您可以选择绘制类型1比例或类型2比例。如果您主要对样本排序感兴趣，则比例1是您最合适的选择。如果您主要对分类单元的排序感兴趣，那么，比例2是您最合适的选择。类型1缩放强调样本之间的关系。基本特征(Borcard等人。Legendre和Legendre 2012)是样本，充当响应变量(列)的质心，样本点之间的距离表示它们的v2距离。解释基于：(1)彼此接近的样本点在其相对频率方面可能相对相似；(2)表示分类或定性变量状态的质心附近的样本更有可能具有该物种/分类单元的状态；(3)样本点在表示定量解释变量的向量上的直角投影近似于该样本实现的变量值。类型2测量强调响应变量之间的关系。其基本特征是响应变量(柱)作为样本的质心，响应变量点之间的距离表示它们的v2距离。该解释基于：(1)彼此接近的物种/分类点沿样本可能具有相似的相对频率；(2)响应变量(物种/分类单元)越接近表示分类解释变量的状态的质心，该响应变量在该状态下越有可能具有较高的值；以及(3)表示响应变量(物种/分类单元)的点在表示解释变量的箭头上的直角投影表示最大值(最佳)的位置；以及(3)表示响应变量(物种/分类单元)的点在表示解释变量的箭头上的直角投影表示最大值(最佳)的位置。类型2测量强调响应变量之间的关系。其基本特征是响应变量(柱)作为样本的质心，响应变量点之间的距离表示它们的v2距离。该解释基于：(1)彼此接近的物种/分类点沿样本可能具有相似的相对频率；(2)响应变量(物种/分类单元)越接近表示分类解释变量的状态的质心，该响应变量在该状态下越有可能具有较高的值；以及(3)表示响应变量(物种/分类单元)的点在表示解释变量的箭头上的直角投影表示最大值(最佳)的位置；以及(3)表示响应变量(物种/分类单元)的点在表示解释变量的箭头上的直角投影表示最大值(最佳)的位置。在受约束的排序中，还有一个额外的源数据可用于绘图：解释性(环境)变量。因此，有三种不同的数据可供绘制：样本、响应变量和解释变量。如果您显示来自三个来源的所有数据，您就有了“TriPlot”。如果您显示两个源数据，则有‘bilot’。以下R代码生成类型2缩放三线图：
> plot(rda_hell, display=c("sp", "lc", "cn"),main="Triplot RDA - scaling 2")
在上述代码中，“sp”代表物种(display=“sp”)，“cn”代表约束或解释变量(display=“cn”)。vegan包中有两个样本得分：物种/分类群的加权总和(DISPLAY=“WA”)和拟合样本得分或“LC得分\”(DISPLAY=“Lc”)，即解释变量的线性组合，您可以选择其中之一显示在一个曲线图中。下列R代码用于添加箭头以显示物种/分类群。SCORES()函数的参数“CHOICES=c(1，2)”是要绘制的轴。默认值绘制的是轴1和轴2。因此，我们可以省略代码“CHOICES=c(1，2)”(图7.27)。
> taxa_scores <- scores(rda_hell, choices=c(1,2), display="sp")
> arrows(0, 0, taxa_scores[,1], taxa_scores[,2], length=0, lty=1, col="red")
在上面的三线图RDA-Scaling 2中，绿色空心三角形表示样本，蓝色十字表示分类解释变量(例如，男性、女性或吸烟者、不吸烟者)的状态，蓝色箭头表示定量解释变量(这里是PackYears和Age)，箭头指示其增加的方向，物种/分类群显示为红色加号。为了检验由解释变量解释的吸烟者社区数据变化的显著性，我们可以通过vegan包中的函数anova.cca()进行蒙特卡洛排列测试。此函数可以测试全局模型(默认)、所有轴(BY=“AXIS”)、单个解释变量(BY=“TERMS”)、第一个约束轴(FIRST=TRUE)或在删除模型中所有其他变量的变量(BY=“margin”)后由单个解释变量解释的变量的显著性。我们通过分别测试全局模型、每个轴和每个解释变量来说明它的性能。为确保每次运行时获得相同的置换测试结果，请设置相同的seed.。
> set.seed (123)
下面的R代码测试全局模型的重要性。参数“step”指定排列的最小数量。（置换检验）
> anova(rda_hell, step=1000)
Permutation test for rda under reduced model
Permutation: free
Number of permutations: 999
Model: rda(formula = abund_hell ~ SmokingStatus + Age + Sex + PackYears, data = throat_meta)
Df Variance F Pr(>F)
Model 4 0.048 1.6 0.006 **
Residual 55 0.414
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
这项全局模型在统计上具有重要意义。现在让我们测试一下每个轴。
> anova(rda_hell, by="axis", step=1000)
Permutation test for rda under reduced model
Marginal tests for axes
Permutation: free
Number of permutations: 999
Model: rda(formula = abund_hell ~ SmokingStatus + Age + Sex + PackYears, data
= throat_meta)
Df Variance F Pr(>F)
RDA1 1 0.024 3.18 0.001 ***
RDA2 1 0.013 1.76 0.032 *
RDA3 1 0.006 0.79 0.735
RDA4 1 0.005 0.67 0.886
Residual 55 0.414
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
我们可以看到第一个和第二个轴是重要的。现在让我们测试每个解释变量的重要性。
> anova(rda_hell, by="terms", step=1000)
Permutation test for rda under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999
Model: rda(formula = abund_hell ~ SmokingStatus + Age + Sex + PackYears, data
= throat_meta)
Df Variance F Pr(>F)
SmokingStatus 1 0.022 2.86 0.001 ***
Age 1 0.006 0.75 0.787
Sex 1 0.014 1.92 0.017 *
PackYears 1 0.007 0.88 0.601
Residual 55 0.414
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
我们可以看到，SmokingStatus和Sex都很重要。最后，我们使用正向选择来减少进入分析的解释变量的数量，同时优化它们所解释的变量。目前，RDA中提供了三个用于正向选择的函数：ordistep()、ordiR2step()和forward.sel()，前两个函数来自素食包。第三个可以从“adespace”包中获得。这些函数对变量选择使用不同的标准。函数ordistep()使用AIC准则和蒙特卡罗置换检验中的p值进行变量比较，函数ordiR2step()使用R2 adj。函数forward.sel()使用预先选择的α的显著性级别作为变量选择标准。尽管函数forward.sel()设置参数的逻辑不同，但与其他两个函数相比，返回的结果基本上与ordiR2step()相同。因此，我们只说明前两个函数。过程ordistep()适用于函数rda()、cca()或cmdscale()。OrdiR2step()只能应用于rda()和capscale()，但不能应用于cca()，因为cca()不返回r2adj。以下R代码使用函数ordistep()运行正向选择。还提供了分步选择和向后选择的选项。此功能允许使用因子。
> step_forward <- ordistep(rda(abund_hell ~ 1, data=throat_meta), scope=formula(rda_hell), direction="forward", pstep=1000)
Start: abund_hell ~ 1
Df AIC F Pr(>F)
+ SmokingStatus 1 -46.1 2.83 0.005 **
+ Sex 1 -45.3 2.01 0.015 *
+ PackYears 1 -45.1 1.80 0.045 *
+ Age 1 -44.1 0.83 0.720
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Step: abund_hell ~ SmokingStatus
Df AIC F Pr(>F)
+ Sex 1 -46.0 1.87 0.025 *
+ PackYears 1 -45.0 0.87 0.650
+ Age 1 -44.9 0.74 0.820
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Step: abund_hell ~ SmokingStatus + Sex
Df AIC F Pr(>F)
+ PackYears 1 -45.0 0.85 0.66
+ Age 1 -44.9 0.81 0.69
此过程选择变量SmokingStatus和Sex。以下向前选择使用函数ordiR2step()。此过程包含新变量的默认设置基于R2调整及其与全局模型(包含所有变量)。如果新变量不重要或包括该新变量的模型的R2调整将超过全局模型的R2调整，则停止选择。
> step_forward <- ordiR2step(rda(abund_hell ~ 1, data=throat_meta), scope=formula(rda_hell), direction="forward",pstep=1000)
Step: R2.adj= 0
Call: abund_hell ~ 1
R2.adjusted
<All variables> 0.039094
+ SmokingStatus 0.030093
+ Sex 0.016765
+ PackYears 0.013326
<none> 0.000000
+ Age -0.002834
Df AIC F Pr(>F)
+ SmokingStatus 1 -46.1 2.83 0.002 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Step: R2.adj= 0.03009
Call: abund_hell ~ SmokingStatus
R2.adjusted
+ Sex 0.04449
<All variables> 0.03909
<none> 0.03009
+ PackYears 0.02786
+ Age 0.02574
使用函数ordiR2step()通过向前选择选择变量吸烟状态。运行正向选择程序后，我们拟合最终的简约RDA模型，并用全局模型、每个轴和每个变量进行测试：
> rda_final<-rda(abund_hell~SmokingStatus+Sex, data=throat_meta)
> anova(rda_final, step=1000)
Permutation test for rda under reduced model
Permutation: free
Number of permutations: 999
Model: rda(formula = abund_hell ~ SmokingStatus + Sex, data = throat_meta)
Df Variance F Pr(>F)
Model 2 0.036 2.37 0.001 ***
Residual 57 0.427
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> anova(rda_final, by="axis", step=1000)
Permutation test for rda under reduced model
Marginal tests for axes
Permutation: free
Number of permutations: 999
Model: rda(formula = abund_hell ~ SmokingStatus + Sex, data = throat_meta)
Df Variance F Pr(>F)
RDA1 1 0.023 3.01 0.001 ***
RDA2 1 0.013 1.74 0.031 *
Residual 57 0.427
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> anova(rda_final, by="terms", step=1000)
Permutation test for rda under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999
Model: rda(formula = abund_hell ~ SmokingStatus + Sex, data = throat_meta)
Df Variance F Pr(>F)
SmokingStatus 1 0.022 2.87 0.003 **
Sex 1 0.014 1.87 0.022 *
Residual 57 0.427
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Constrained Correspondence Analysis (CCA)：CCA也称为典范对应分析。自1986年引入以来，它一直是群落生态学中最流行的排序方法之一，并被微生物组研究人员所接受。就像RDA与PCA相关一样，CAP与PCoA相关，CCA与CA相关。CCA共享CA的基本属性，并将它们组合成受约束的排序。Cca由函数cca()执行。它的算法是基于Legendre和Legendre的：它保留了样本之间的v2距离，分类群被表示为Triplot中的点(Borcard等人。(2011年)。计算出的v2距离对约束变量进行加权线性回归，拟合值通过奇异值分解(SvD)进行对应分析。因此，它是RDA的加权形式。与RDA一样，CCA也有两种语法。一种是简单(默认)矩阵语法：cca（X，Y，Z）。其中，X=社区数据矩阵或数据框，必须给出；Y=约束矩阵或数据框，通常为环境变量(可省略)；如果提供矩阵Y，则用于进行CCA，如果不提供，则计算CA。Z=条件矩阵或数据框(也可以省略)。如果提供矩阵Z，则效果将从社区矩阵中分割出来。另一种是公式语法：cca （formula，data，na:action = na:fail，subset =NULL ）。其中，公式是典型的模型公式：公式左侧必须是社区数据矩阵(X)；右侧定义约束模型：约束变量可以包含有序或无序因子、变量间的相互作用和变量的函数；条件变量可以在特定的函数条件下给出，用于分析前分割出的条件变量(协变量)。在本节中，我们将使用吸烟者数据通过vegan包中的函数cca()进行CCA。我们在前面提到过，要进行CCA，必须提供环境变量矩阵，否则cca()函数将计算CA。回想一下，吸烟者数据有两个数据集：throat.otu.tab(社区丰度数据框)和throat.meta(元数据，包括两个二进制变量SmokingStatus和Sex，以及两个连续变量Age和PackYears)。要运行CCA，请现在加载vegan，并且不像我们在RDA中所做的那样，使用Hellinger方法转换社区丰富的数据。否则，无法计算X2距离，也无法解释结果。在下面的CCA中，吸烟者丰度数据受喉咙中的四个环境变量的约束。元数据包括SmokingStatus、Age、Sex和PackYears。
> smoker_cca <- cca(throat.otu.tab ~ ., throat_meta)
> smoker_cca
Call: cca(formula = throat.otu.tab ~ SmokingStatus + Age+ Sex + PackYears, data = throat_meta)
Inertia Proportion Rank
Total 4.9330 1.0000
Constrained 0.3782 0.0767 4
Unconstrained 4.5548 0.9233 55
Inertia is mean squared contingency coefficient
Eigenvalues for constrained axes:
CCA1 CCA2 CCA3 CCA4
0.1517 0.0912 0.0781 0.0572
Eigenvalues for unconstrained axes:
CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8
0.457 0.358 0.325 0.299 0.237 0.187 0.165 0.160
(Showed only 8 of all 55 unconstrained eigenvalues)
第一节之后的“Call”：给出“mean squared contingency coefficients（均方偶合系数）”的分析。它们在CCA中的作用与总方差在RDA中的作用相同。由于社区矩阵被转换为v2距离，因此条目是权变系数。矩阵进行加权回归前的总变异为4.9330，这是可以解释的变异。加权回归后解释的群落矩阵中的变异为0.3782；这是将由CCA中的轴来解释的变异。回归残差的方差是4.5548；这是CCA中的轴不能解释的方差，它可以受到CA的影响。在这里，我们注意到CCA不是很成功，因为CCA中只有0.3782/4.9330或0.0767(见列比例和行约束)的数据总变异被所有四个变量捕获。可以通过Summary()函数找到由特定轴解释的方差：
> summary(smoker_cca)
Call:
cca(formula = throat.otu.tab ~ SmokingStatus + Age + Sex + PackYears, data =
throat_meta)
Partitioning of mean squared contingency coefficient:
Inertia Proportion
Total 4.933 1.0000
Constrained 0.378 0.0767
Unconstrained 4.555 0.9233
Eigenvalues,and their contribution to the mean squared contingency
coefficient
Importance of components:
CCA1 CCA2 CCA3 CCA4 CA1 CA2
Eigenvalue 0.1517 0.0912 0.0781 0.0572 0.4567 0.3581
Proportion Explained 0.0307 0.0185 0.0158 0.0116 0.0926 0.0726
Cumulative Proportion 0.0307 0.0492 0.0651 0.0767 0.1692 0.2418
(…)
Accumulated constrained eigenvalues
Importance of components:
CCA1 CCA2 CCA3 CCA4
Eigenvalue 0.152 0.0912 0.0781 0.0572
Proportion Explained 0.401 0.2411 0.2066 0.1512
Cumulative Proportion 0.401 0.6422 0.8488 1.0000
(…)
Biplot scores for constraining variables
CCA1 CCA2 CCA3 CCA4 CA1 CA2
SmokingStatusSmoker 0.831 0.452 0.0597 0.3181 0 0
Age 0.202 -0.179 -0.7738 0.5729 0 0
SexMale 0.683 -0.730 0.0444 0.0112 0 0
PackYears 0.615 0.163 -0.7704 0.0303 0 0
Centroids for factor constraints
CCA1 CCA2 CCA3 CCA4 CA1 CA2
SmokingStatusNonSmoker -0.746 -0.406 -0.0531 -0.28569 0 0
SmokingStatusSmoker 0.926 0.503 0.0660 0.35459 0 0
SexFemale -0.922 0.985 -0.0599 -0.01516 0 0
SexMale 0.505 -0.540 0.0328 0.00831 0 0
前四个约束轴(CCA1、CCA2、CCA3、CCA4)解释的方差分别为3.07、1.85、1.58和1.16%。第一无约束轴(CA1)解释的方差为9.26%。“Accumulated constrained eigenvalues（累积约束特征值）”部分提供了与投影相关的特征值。因为我们的环境数据框架中有四个变量，所以有四个受约束的特征值。如下图所示，第一个轴约占约束变化的40%，第二个轴占24%，第三个轴占21%，第四个轴占15%。接下来，我们使用plot()函数绘制类型1缩放的CCA结果。在该图中，我们希望显示线性约束或“LC scores”(display = “lc”)和因子变量级别的质心(display = “cn”)。
plot(smoker_cca,scaling=1,display = c("lc","cn"),main="Biplot CCA-scaling 1")
类型1缩放显示四组示例，其中吸烟者组链接到PackYears。在这个分析中，第一个轴与PackYears的增加有关，而第二个轴与年龄的减少有关(图7.28)。
与RDA类似，让我们对全局模型、每个轴和每个解释变量进行排列测试。为确保获得相同的置换检验结果，每次运行时，请设置相同的seed。
> set.seed (123)
下面的R代码测试全局模型的重要性。
> anova(smoker_cca, step=1000)
Permutation test for cca under reduced model
Permutation: free
Number of permutations: 999
Model: cca(formula = throat.otu.tab ~ SmokingStatus + Age + Sex + PackYears,
data = throat_meta)
Df ChiSquare F Pr(>F)
Model 4 0.38 1.14 0.2
Residual 55 4.55
这一结果在统计上并不显著。然后测试每个轴的显著性。
> anova (smoker_cca, by = 'axis',step=1000)
Permutation test for cca under reduced model
Marginal tests for axes
Permutation: free
Number of permutations: 999
Model: cca(formula = throat.otu.tab ~ SmokingStatus + Age + Sex + PackYears,
data = throat_meta)
Df ChiSquare F Pr(>F)
CCA1 1 0.15 1.83 0.004 **
CCA2 1 0.09 1.10 0.350
CCA3 1 0.08 0.94 0.611
CCA4 1 0.06 0.69 0.959
Residual 55 4.55
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
我们可以看到第一个轴是显著的。现在让我们测试每个解释变量的显著性。
> anova (smoker_cca, by = 'terms')
Permutation test for cca under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999
Model: cca(formula = throat.otu.tab ~ SmokingStatus + Age + Sex + PackYears,
data = throat_meta)
Df ChiSquare F Pr(>F)
SmokingStatus 1 0.13 1.56 0.009 **
Age 1 0.07 0.89 0.751
Sex 1 0.10 1.27 0.104
PackYears 1 0.07 0.85 0.768
Residual 55 4.55
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
我们可以看到，SmokingStatus非常显著，p值为0.009。最后，让我们使用vegan包中的函数ordistep()进行前向选择。正如我们在RDA中所述，函数ordistep()适用于函数rda()、cca()或cmdscale()。变量的比较基于AIC标准和蒙特卡罗排列检验的p值。
> ordistep(cca(throat.otu.tab ~ 1,data=throat_meta),scope=formula(smoker_cca), direction="forward",pstep=1000)
Start: throat.otu.tab ~ 1
Df AIC F Pr(>F)
+ SmokingStatus 1 539.05 1.5638 0.015 *
+ Sex 1 539.17 1.4378 0.025 *
+ PackYears 1 539.34 1.2775 0.225
+ Age 1 539.73 0.8914 0.795
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Step: throat.otu.tab ~ SmokingStatus
Df AIC F Pr(>F)
+ Sex 1 539.71 1.2838 0.075 .
+ PackYears 1 540.04 0.9688 0.555
+ Age 1 540.12 0.8877 0.760
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Call: cca(formula = throat.otu.tab ~ SmokingStatus, data = throat_meta)
Inertia Proportion Rank
Total 4.93301 1.00000
Constrained 0.12951 0.02625 1
Unconstrained 4.80350 0.97375 58
Inertia is mean squared contingency coefficient
Eigenvalues for constrained axes:
CCA1
0.12951
Eigenvalues for unconstrained axes:
CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8
0.4840 0.3834 0.3322 0.3025 0.2486 0.1884 0.1774 0.1614
(Showed only 8 of all 58 unconstrained eigenvalues)
> ordistep(cca(throat.otu.tab ~ 1, data=throat_meta),scope=formula(smoker_cca), direction="forward", pstep=1000)
Start: throat.otu.tab ~ 1
Df AIC F Pr(>F)
+ SmokingStatus 1 539 1.56 0.005 **
+ Sex 1 539 1.44 0.045 *
+ PackYears 1 539 1.28 0.190
+ Age 1 540 0.89 0.710
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Step: throat.otu.tab ~ SmokingStatus
Df AIC F Pr(>F)
+ Sex 1 540 1.28 0.08 .
+ PackYears 1 540 0.97 0.54
+ Age 1 540 0.89 0.78
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Call: cca(formula = throat.otu.tab ~ SmokingStatus, data = throat_meta)
Inertia Proportion Rank
Total 4.9330 1.0000
Constrained 0.1295 0.0263 1
Unconstrained 4.8035 0.9737 58
Inertia is mean squared contingency coefficient
Eigenvalues for constrained axes:
CCA1
0.1295
Eigenvalues for unconstrained axes:
CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8
0.484 0.383 0.332 0.302 0.249 0.188 0.177 0.161
(Showed only 8 of all 58 unconstrained eigenvalues)
所选变量为SmokingStatus。最后的模型安装在下面。
> (smoker_cca_final <- cca(throat.otu.tab ~ SmokingStatus, data=throat_meta))
> anova.cca(smoker_cca_final, step=1000)
Call: cca(formula = throat.otu.tab ~ SmokingStatus, data = throat_meta)
Inertia Proportion Rank
Total 4.9330 1.0000
Constrained 0.1295 0.0263 1
Unconstrained 4.8035 0.9737 58
Inertia is mean squared contingency coefficient
Eigenvalues for constrained axes:
CCA1
0.1295
Eigenvalues for unconstrained axes:
CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8
0.484 0.383 0.332 0.302 0.249 0.188 0.177 0.161
(Showed only 8 of all 58 unconstrained eigenvalues)
> anova.cca(smoker_cca_final, step=1000)
Permutation test for cca under reduced model
Permutation: free
Number of permutations: 999
Model: cca(formula = throat.otu.tab ~ SmokingStatus, data = throat_meta)
Df ChiSquare F Pr(>F)
Model 1 0.13 1.56 0.007 **
Residual 58 4.80
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
吸烟状况达到同样的显著性水平，p值为0.007。然而，节俭模式得到了回报，留下了更大的剩余。
Constrained Analysis of Principal Coordinates (CAP)：CAP(也称为vegan包装中的约束邻近分析)，是一种类似于RDA的排序方法。它只是对主坐标分析(或公制多维尺度)结果的冗余分析。CAP允许非欧几里得的不同指数，例如曼哈顿距离或Bray-Curtis距离。如果将欧几里得距离指定为排序方法，则结果将与RDA相同。vegan包中的函数capscale()用于实现CAP。它需要一个相异矩阵作为输入数据集，可以使用函数vegdist()、dist()或产生相似矩阵的任何其他方法来计算。它需要两个步骤：首先使用cmdscale()函数对相异矩阵进行排序，然后使用RDA对结果进行分析。与既可以使用矩阵语法又可以使用公式语法的RDA不同，函数capscale()只能使用公式语法调用。下面列出了函数capscale()的一种用法。capscale （formula， data， distance =“bray”，dfun = vegdist ）。其中，formula是在RDA()和CCA()中定义的典型模型公式。公式的左侧必须是社区数据矩阵(框架)或相异矩阵，可以从函数vegdist()或dist()估计。如果公式的左侧是数据矩阵(框)而不是相异度矩阵，则必须提供相异度(或距离)指标作为距离的输入。公式的右侧定义了约束。约束变量可以是连续变量、因素、交互作用项或用作定义要划分的变量的特殊项条件。数据是包含模型公式右侧变量的数据框。 dfun是使用的距离或相异函数。基本的CAP可以使用以下代码来完成。约束变量包括两个二进制变量(SmokingStatus和Sex)和两个连续变量(PackYears和Age)。特殊的条件(年龄)是用来分割年龄效应的。
> throat_cap <- capscale(throat.otu.tab ~ SmokingStatus + Sex + PackYears +Condition(Age), throat_meta, dist="bray")
> throat_cap
Call: capscale(formula = throat.otu.tab ~ SmokingStatus + Sex + PackYears +
Condition(Age), data = throat_meta, distance = "bray")
Inertia Proportion Eigenvals Rank
Total 14.0932 1.0000 14.5356
Conditional 0.2020 0.0143 0.2084 1
Constrained 1.2438 0.0883 1.2623 3
Unconstrained 12.6473 0.8974 13.0648 46
Imaginary -0.4425 13
Inertia is squared Bray distance
Eigenvalues for constrained axes:
CAP1 CAP2 CAP3
0.731 0.376 0.155
Eigenvalues for unconstrained axes:
MDS1 MDS2 MDS3 MDS4 MDS5 MDS6 MDS7 MDS8
2.278 1.994 1.144 0.925 0.749 0.615 0.534 0.473
(Showed only 8 of all 46 unconstrained eigenvalues)
前三个轴称为“CAP1”、“CAP2”和“CAP3”，然后是原始MDS。我们可以看到，三个约束变量(SmokingStatus、Sex和PackYears)可以解释整个数据集总变异的8.83%。
> anova(throat_cap)
Permutation test for capscale under reduced model
Permutation: free
Number of permutations: 999
Model: capscale(formula = throat.otu.tab ~ SmokingStatus + Sex + PackYears +
Condition(Age), data = throat_meta, distance = "bray")
Df SumOfSqs F Pr(>F)
Model 3 1.24 1.8 0.004 **
Residual 55 12.65
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
该模型具有统计学意义，p值为0.004。默认的Plot()函数生成以下图形，它不是信息性的(图7.29)。
> plot(throat_cap)
由于我们对吸烟者和不吸烟者之间的不同之处很感兴趣，所以让我们在这里提取组信息。
> groups <- throat.meta$SmokingStatus
> groups
[1] NonSmoker Smoker Smoker Smoker Smoker Smoker
[7] NonSmoker NonSmoker NonSmoker NonSmoker Smoker NonSmoker
[13] NonSmoker Smoker NonSmoker Smoker Smoker NonSmoker
[19] NonSmoker NonSmoker NonSmoker NonSmoker NonSmoker NonSmoker
[25] NonSmoker NonSmoker NonSmoker NonSmoker NonSmoker NonSmoker
[31] NonSmoker NonSmoker NonSmoker Smoker NonSmoker NonSmoker
[37] NonSmoker NonSmoker NonSmoker Smoker NonSmoker Smoker
[43] NonSmoker Smoker Smoker Smoker Smoker Smoker
[49] Smoker Smoker Smoker Smoker Smoker Smoker
[55] Smoker Smoker Smoker NonSmoker Smoker Smoker
Levels: NonSmoker Smoker
在下面的一串R代码中，函数 plot()生成一个空的CAP排序图；函数points() 将点添加到由该图创建的排序图(低级绘图函数)。函数ordispider()通过将组的各个成员与组质心连接来创建蜘蛛图。函数ordiellse()像信封一样用椭圆包围组内的点云(图7.30)。
> plot(throat_cap, type="n")
> points(throat_cap, col=as.numeric(as.factor(groups)), pch=as.numeric(as.factor(groups)))
> ordispider(throat_cap, groups, lty=2, col="grey", label=T)
> ordiellipse(throat_cap, groups, lty=2, col="grey", label=F)
与RDA和CCA类似，让我们对全局模型、每个轴和每个解释变量进行置换检验。为确保每次运行时获得相同的置换检验结果，请设置相同的seed。
> set.seed (123)
以下代码测试全局模型的重要性。
> anova(throat_cap, step=1000)
Permutation test for capscale under reduced model
Permutation: free
Number of permutations: 999
Model: capscale(formula = throat.otu.tab ~ SmokingStatus + Sex + PackYears +
Condition(Age), data = throat_meta, distance = "bray")
Df SumOfSqs F Pr(>F)
Model 3 1.24 1.8 0.003 **
Residual 55 12.65
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
其结果具有统计学意义，p值为0.003。然后测试每个轴的显著性。
> anova(throat_cap, by="axis", step=1000)
Permutation test for capscale under reduced model
Marginal tests for axes
Permutation: free
Number of permutations: 999
Model: capscale(formula = throat.otu.tab ~ SmokingStatus + Sex + PackYears +
Condition(Age), data = throat_meta, distance = "bray")
Df SumOfSqs F Pr(>F)
CAP1 1 0.73 3.08 0.001 ***
CAP2 1 0.38 1.58 0.073 .
CAP3 1 0.16 0.65 0.869
Residual 55 13.06
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
我们可以看到，第一个轴是显着的，p值为0.001，第二个轴是边际显着的，p值为0.073。现在让我们测试每个解释变量的显著性。
> anova(throat_cap, by="terms", step=1000)
Permutation test for capscale under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999
Model: capscale(formula = throat.otu.tab ~ SmokingStatus + Sex + PackYears +
Condition(Age), data = throat_meta, distance = "bray")
Df SumOfSqs F Pr(>F)
SmokingStatus 1 0.63 2.72 0.003 **
Sex 1 0.44 1.89 0.035 *
PackYears 1 0.18 0.79 0.682
Residual 55 12.65
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
我们可以看到，吸烟状态和性别都是显著的，p值分别为0.003和0.035。现在，让我们使用vegan包中的函数ordistep()进行前向选择。正如我们在RDA和CCA中所述，该函数基于AIC标准和蒙特卡罗置换测试的p值比较变量。
> step_forward <- ordistep(capscale(throat.otu.tab ~ 1, data=throat_meta), scope=formula(throat_cap), direction="forward", pstep=1000)
Start: throat.otu.tab ~ 1
Df AIC F Pr(>F)
+ SmokingStatus 1 712 2.05 0.035 *
+ Sex 1 712 1.68 0.090 .
+ PackYears 1 712 1.27 0.215
+ Condition(Age) 1 713 0.00
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Step: throat.otu.tab ~ SmokingStatus
Df AIC F Pr(>F)
+ Sex 1 712 1.53 0.14
+ PackYears 1 713 0.58 0.78
+ Condition(Age) 1 713 0.00
所选变量为吸烟状态。最后的模型安装在下面。
> cap_final<- capscale(throat.otu.tab ~ SmokingStatus, throat_meta, dist="bray")
> anova(cap_final, step=1000)
Permutation test for capscale under reduced model
Permutation: free
Number of permutations: 999
Model: capscale(formula = throat.otu.tab ~ SmokingStatus, data = throat_meta,
distance = "bray")
Df SumOfSqs F Pr(>F)
Model 1 0.65 2.8 0.004 **
Residual 58 13.44
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
吸烟状况达到显著水平，p值为0.004，但简约模型获得了较大的残差。
7.5 Summary and Discussion
在本章中，我们使用鼠和人类数据集来说明对微生物组数据的探索性分析。我们首先使用Phyloseq软件包对richness、 abundance bar, heatmap, network and phylogenetic tree这五个图进行了说明。然后，介绍了生态学和微生物群研究中常用的几类聚类方法，重点阐述了四种聚类方法(single linkage agglomerative, complete linkage agglomerative, average linkage agglomerative, and Ward’s minimum variance)。我们还简要描述了聚类、排序和距离测度之间的关系。最后，我们说明了最常见的无约束和受约束的排序：PCA、PCoA、NMDS、CA、RDA、CCA和CAP。无约束坐标的特点被认为是“探索性的”。然而，约束排序的能力不仅仅是探索性的数据分析，还可以成为假设检验。介绍了无约束坐标和约束坐标的不同特点。我们通过使用专门开发的微生物群普查数据分析软件包“phyloseq”和生态学中常用的软件包(如vegan软件包)来说明对微生物组数据的探索性分析和假设检验(在约束坐标方面)。扩展函数Cleanplot.pca()和evlot()是很有吸引力的，评价排序性能的两个标准(Kaiser-Guttman criterion and broken-stick model）有助于排序分析的评价。
