微生物群落组成的研究方法可以分为两大类：分类多样性分析和微生物群落组成的多元分析。多变量分析包括各种多变量技术，如聚类和(无约束和有约束的)排序，以及各组之间的假设检验差异。无约束排序虽然涉及到事后假设，但本质上属于探索性分析。约束排序是一种假设检验。在这一章中，我们将使用各种图形技术来探索分类学的多样性，并使用聚类和排序技术来探索微生物组的组成。

7.1 小鼠和人类的数据集

VDR−/−小鼠数据集：我们将继续使用VDR−/−小鼠数据集，小鼠肠道菌群数据，从粪便和盲肠粪便样本中采集。这里使用的是粪便样本。

吸烟人数数据集：第二个数据集来自Charlson等人。其中包括吸烟对上呼吸道微生物群的影响的研究。原始数据集包含来自咽喉和鼻子微生物群的样本，以及来自身体两侧。本章使用的数据集来自左侧喉部微生物群。包括60名受试者(32名非吸烟者和28名吸烟者)。数据集包括三个数据： abundance count, tree, and meta-data。适用于图解树木标绘和约束排序分析。

7.2  探索性分析与图形总结

微生物组数据可以通过各种图表进行探索。在这里，我们使用上述两个数据集说明了五个常用的图：richness, abundance bar, heatmap, network and phylogenetic tree。这些是由Phyloseq软件包生成的。Phyloseq软件包是一个用于导入、存储、分析和以图形方式显示复杂的系统发育测序数据的工具。此软件包中使用的输入数据可以是OTUS或丰富的计数数据。该软件包使用先进/灵活的图形系统(Ggplot2)轻松生成出版质量的数据图形。

Plot Richness：可以使用phyloseq包中的函数plot_richness()通过图形汇总估计的alpha多样性。虽然它的名字暗示着绘制“丰富度”，通常指的是绘制样本或环境中物种/类群/OTUS的总数，但实际上该功能不仅绘制丰富度，还生成观察到的和其他估计的多样性的数字。首先，我们需要加载phyloseq和ggplot2软件包，以及vdr−/−小鼠数据集。

> library(phyloseq)

> library(ggplot2)

> abund_table=read.csv("VdrFecalGenusCounts.csv",row.names=1,check.names=FALSE)as

0wertyuiop

> abund_table<-t(abund_table)

使用phyloseq包时的关键步骤是构建phyloseq类对象。下面的R代码块使用构造函数phyloseq()构建一个名为physeq的phyloseq类对象。Phyloseq类对象是从其组件数据构建的：OTU表、示例数据、分类表和Phylo树。作为实验级对象，必须提供两个或多个组件数据对象。论据的顺序并不重要。

> meta_table <-data.frame(row.names=rownames(abund_table),t(as.data.frame(strsplit(rownames(abund_table),"_"))))

> meta_table$Group <- with(meta_table,ifelse(as.factor(X2)%in% c(11,12,13,14,15),c("Vdr-/-"), c("WT")))

> #Convert the data to phyloseq format

> OTU = otu_table(as.matrix(abund_table), taxa_are_rows = FALSE)

> SAM = sample_data(meta_table)

> physeq<-merge_phyloseq(phyloseq(OTU),SAM)

> physeq

phyloseq-class experiment-level object

otu_table() OTU Table:  [ 248 taxa and 8 samples ]

sample_data() Sample Data:  [ 8 samples by 4 sample variables ]

在构建Physeq对象之后，我们可以通过ggplot2包使用函数plot_richness()来绘制观察到的和估计的alpha差异度。

> plot_richness(physeq, x = "Group", color = "Group")

输入数据“physeq”是必需的，它是phyloseq类，或者可替换地，是OTU表类。可选参数“x”是映射到水平轴的变量；x可以是字符串，也可以是向量。默认值为“Samples”，它会将每个样本的名称映射到绘图中的单独水平位置。在本例中，x=“Group”会将组成员身份映射到x轴。参数颜色也是可选的。它将指定映射到不同颜色的示例变量。与参数x类似，它可以是单个字符串或向量。在α分集估计中，噪声可以被修剪；然而，由于许多丰富度估计，甚至是“观察到的”丰富度，都高度依赖于单例的数量。因此，如果您想要有意义的结果，您必须使用未修剪的数据集。图7.1由上述R代码生成。您还可以选择要绘制的Alpha多样性度量。例如，下面的R代码仅绘制Chao1和Shannon多样性。

> plot_richness(physeq, measures = c("Chao1", "Shannon"),x = "Group", color= "Group")

Plot Abundance Bar：Phyloseq函数lot_bar()功能强大且灵活。它的主要目的是快速、轻松地创建提供信息的汇总图形，显示实验样本之间分类群丰度的差异。下面给出一种用法：plot_bar(physeq, x=“Sample”, y =“Abundance”, fill=NULL, title=NULL,facet_grid=NULL)

“x”和“y”参数都是可选的字符串。数据中的变量将分别映射到x轴和y轴。通常，“y”参数将是“丰度”，以定量显示每个OTU/组的丰度值。另外两个复杂的自定义可选参数是“Fill”和“facet_grid”。“Fill”选项是一个字符串，用于指定哪个示例变量将用于映射到条的填充颜色。“facet_grid”选项是一个公式对象，用于指定您想要在“ggplot2”图形中显示的刻面。让我们先调整一下默认主题。

> theme_set(theme_bw())

没有给定任何参数的默认条形图将绘制每个样本分别映射到x轴，而丰度值映射到y轴。每个OTU/样品的丰度值按从大到小的顺序堆叠，由一条细水平线分隔。

> plot_bar(physeq)

以下R代码用于填充每个样品所属的VDR、−/−和WT组的不同颜色。

> plot_bar(physeq, fill="Group")

在图7.2中，所有丰度值都被绘制出来，如果你只想知道前五个最丰富的细菌，那么你可以使用下面的R代码块(图7.3)。

> TopNGenus <- names(sort(taxa_sums(physeq), TRUE)[1:5])

> Top5Genus <- prune_taxa(TopNGenus,physeq)

> plot_bar(Top5Genus, fill="Group", facet_grid=*Group)

到目前为止，我们使用变量group演示了“fill”和“facet_grid”选项。实际上，要可视化VDRPhyloseq和WT组之间的样本丰度差异，更有意义的方法是使用“属”填充此处的颜色，这需要在−/−对象中包含分类数据。然而，我们没有这项研究的分类学数据。感兴趣的读者可以使用下面的示例代码进行自己的学习

> plot_bar(Top5Genus, fill="Genus", facet_grid=*Group）

Plot Heatmap：Rajaram和Oono已经展示了如何使用排序方法而不是层次聚类分析(Rajaram)创建热图来组织行和列。基于排序的排序方法比层次聚类方法更适合于微生物组数据的表示。Phyloseq在线教程中有许多有用的Phyloseq热图图形示例。感兴趣的读者可以参考这些示例热图来完成自己的工作。在这里，我们重点说明基于NMDS和PCA排序方法的plot_heatmap()函数。plot_heatmap()函数的一种用法如下。

plot_heatmap(physeq, method = “NMDS”, distance = “bray”, sample.label =NULL, taxa.label = NULL, low = “#000033”, high = “#66CCFF”, na.value =“black”)

输入数据参数“physeq”是必需的。它是一个phyloseq类对象(Otu_Table)。方法和距离增加都是可选的。排序方法用于组织热图。生态距离法是用于排序的字符串。“sample.label”和“taxa.label”都是字符串，可分别用于标记样本(水平)轴和重新标记分类/物种/OTU(垂直)轴。低增长和高增长都是字符串，并且都是可选的。它们用于选择R中支持的颜色选项。R理解600多种颜色。您可以在R中键入颜色()以检查名称。

> colors()

可以在R Cookbook(http：//www.cookbook-r.com/graph/Colors_(ggplot2)/)上找到这些颜色的表格摘要。

在热图中，零值被视为NA，并设置为“黑色”以表示背景颜色。low选项用于表示最低非零值，low参数的默认值为深蓝色“#000033”；high选项用于表示最高值，默认值为“#66CCFF”。在热图图中，太多的分类群(在本例中是属)会使图形过于拥挤而看不清。我们需要限制要划分的分类群的数量。下面的R代码选择最丰富的前五个属来制作热图。

> TopNGenus <- names(sort(taxa_sums(physeq), TRUE)[1:5])

> Top5Genus <- prune_taxa(TopNGenus, physeq)

> plot_heatmap(Top5Genus)

下面的R码使用NMDS排序法和Bray-Curtis距离法来绘制前五个属

> (p <- plot_heatmap(Top5Genus, "NMDS", "bray"))

可以尝试不同颜色的热图以满足日志要求。以下是一些可以使用的示例。

> plot_heatmap(Top5Genus, "NMDS", "bray", low="#000033", high="#CCFF66")

> plot_heatmap(Top5Genus, "NMDS", "bray", low="#000033", high="#FF3300")

> plot_heatmap(Top5Genus, "NMDS", "bray", low="#000033", high="#66CCFF")

> plot_heatmap(Top5Genus, "NMDS", "bray", low="#66CCFF", high="#000033",na.value="white")

可以尝试使用生态距离和坐标的不同组合。例如，下面的热图在默认的Bray-Curtis距离上使用PCoA排序(图7.4)。

> plot_heatmap(Top5Genus, "PCoA", "bray")

Plot Network：Phyloseq包中有两个函数可以使用“ggplot2”绘制微生物组网络：plot_network()和plot_net()。如果使用函数 plot_network()，则其第一个参数需要是一个字形对象，并且网络本身应该是使用字形包表示。可以通过phyloseq包使用make_network()函数创建网络对象(参数g)。函数plot_net() 是对plot_network()的性能和接口修订，它的第一个/主要参数是phyloseq-class实例。plot_network()和plot_net()的用法示例如下：plot_network(g, physeq=NULL, type=“samples”, color=“Group”, shape=“Group”)。plot_net(physeq, distance = “bray”, type = “samples”, maxdist = 0.7, color =NULL, shape = NULL)。在plot_network()中，g是必需的igraph类对象，可以由函数make_network()创建，也可以直接由igraph-package创建。可选参数physeq是g所基于的phyloseq类对象。类型选项指示主参数g中表示的网络是samples还是taxa/OTUS。默认值为“Samples”。颜色和形状这两个选项都是可选的。它们用于点的颜色映射和形状映射。plot_net()函数中所需的参数physeq是您想要表示为网络的phyloseq类对象。距离选项是距离方法或已计算的距离类。默认值为“Bray”。选项maxdist表示两个顶点之间与图形中的边连接的最大距离值。默认值为0.7。以下R代码基于默认距离方法“Jaccard”和连接的节点之间的最大距离创建基于图形的网络0.8。“Group”用于颜色和形状映射，以可视化VDR−/−和WT小鼠样本的结构(图7.5)。

> ig <- make_network(physeq, max.dist=0.8)

> plot_network(ig, physeq, color="Group", shape="Group")

上图显示了一些有趣的结构，两个子图分别由来自Vdr−/−和WT小鼠的样本组成。此外，样本之间似乎存在相关性。与 plot_network()函数相比，较新的plot_net() 函数不需要单独的make_network()函数调用，也不需要单独的igraph对象。以下代码基于连接的节点之间的最大距离0.5创建网络。由于我们对VDR−/−和WT鼠标样本的结构很感兴趣，所以我们使用“组”对点进行颜色映射和形状映射(图7.6)。

> plot_net(physeq, maxdist = 0.5, color = "Group", shape="Group")

Plot Phylogenetic Tree：Phyloseq包中的函数plot_tree()旨在方便以图形方式研究系统发育树以及样本数据。对于丰富度较高的样本的系统发育测序，由此函数生成的树将很难阅读和解释。由Phyloseq软件包的作者提出的一个粗略的“经验法则”是使用数据子集，每个小区的OTU不超过200。下面给出此函数的一种用法：plot_tree(physeq, method = “sampledodge”, color = NULL, shape = NULL,ladderize = FALSE).需要输入数据“physeq”。数据应该是系统级的，并且至少包含一个系统发育树成分。您想要绘制和注释系统发育树的正是这些数据。“method”选项是要使用的注释方法的名称。如果观察到来自该分类单元的个体，则默认的“sampledodge”将在树叶旁边绘制点，并为每个样本绘制一个单独的点。“color”和“shape”参数都是可选的。它们提供Physeq中的变量名称，以分别映射到点颜色和点形状。选项“ladderize”用于指定是否对树进行阶梯化，即在绘制之前根据节点的封闭子树的深度对节点进行重新排序。默认值为False，不应用阶梯化。当为true或“right”时，使用“right”阶梯化。当设置为“左”时，将应用“左”阶梯。Smokers数据集用于说明绘制系统发育树。数据来自GUniFrac软件包。让我们首先加载这个包，并使数据可供使用。

> library(GUniFrac)

> data(throat.otu.tab)

> data(throat.tree)

> data(throat.meta)

然后，我们需要通过phyloseq包构建phyloseq类对象，如下所示。

> library(phyloseq)

> #Convert the data to phyloseq format

> OTU = otu_table(as.matrix(throat.otu.tab), taxa_are_rows = FALSE)

> SAM = sample_data(throat.meta)

> TRE <-throat.tree

> physeq<-merge_phyloseq(phyloseq(OTU),SAM,TRE) )

检查数据集中有多少个分类群/OTU。

> ntaxa(physeq)

[1] 856

有856个OTU。太多了，不能用来在系统发育树上作注释。下面的prune_taxa()函数只修剪数据集中的前50个OTU。我们使用这前50个OTU来绘制和注释系统发育树。

> physeq = prune_taxa(taxa_names(physeq)[1:50], physeq)

以下R代码将颜色映射到吸烟状态变量(图7.7)。

> plot_tree(physeq, ladderize = "left", color = "SmokingStatus")

下面的R代码将颜色和形状都映射到吸烟状态变量(图7.8)。

> plot_tree(physeq, ladderize = "left", color = "SmokingStatus",shape "SmokingStatus")

树的上方是垂直方向的树。它是数据到图形的笛卡尔映射。在文献中，经常使用辐射树。我们还可以使用与垂直定向树中相同的映射来创建具有ggplot2的放射状树。不同的是，这次的映射使用的是极坐标(图7.9)。

> plot_tree(physeq,color = "SmokingStatus",shape = "SmokingStatus", ladderize = "left") + coord_polar(theta = "y")

7.3 聚类

①聚类、距离和排序简介：聚类(或分类)和排序是微生物组研究人员和群落生态学家经常使用的两种主要的多变量方法。在某种程度上，这两种方法是相辅相成的。聚类的目的是将样本归入(可能是分层的)类，以降低数据的复杂性(维度)；它可以减少一个维度(x轴)上的所有样本。然而，由于大多数社区数据是连续的，所以二维或三维的数据比一维的数据更容易解释。有了排序，社区数据可以简化为二维或三维。因此，微生物组研究人员和群落生态学家通常希望得到排序。在微生物群和生态学研究中，已经发展了许多基于排序的多变量方法。在进行聚类和排序时，需要通过距离方法提供距离度量。

②聚类：文献中提供了几类聚类方法：顺序或同时算法、聚集或分裂算法、一元方法与综合方法、分层方法与非分层方法、概率方法与非概率方法。在这些分类中，Ward的层次聚类法和k-均值分割法是微生物群研究中常用的方法。分层聚类结果通常用树状图表示。我们将介绍几种不同的聚类方法，包括“单链接聚集聚类”、“完全链接聚集聚类”、“平均链接聚集聚类”和“Ward最小方差聚类”。我们将使用您之前遇到的带有粪便样本的Vdr−/−小鼠数据集来说明聚类分析。这里，用Bray-Curtis距离来说明样本的分类。还会应用其他距离。如果vegan包未加载，请立即加载。让我们首先使用函数decstand()标准化丰度表，并使用vegan包中的函数vegdist()计算所有样本对之间的Bray-Curtis差异。

> abund_table_norm <- decostand(abund_table, "normalize")

> bc_dist<- vegdist(abund_table_norm , method = "bray")

考虑到已经计算出Bray-Curtis不同点，我们现在将层次聚类函数hclust()与四种不同的聚类算法一起应用-““average”, “complete”, “single”链接方法和Ward的聚类方法。

Single Linkage Agglomerative Clustering：“Agglomerate”是指聚集成一个簇；单链接聚集聚类也称为最近邻排序。在每个步骤中，聚类将包含最短成对距离(或最大相似性)的两个样本组合在一起。聚类的结果可以可视化为树状图。导致单一连锁聚类的树状图的缺点通常是样品的链状。由于单个元素彼此接近，所以形成的簇被强迫在一起，而每个簇中的许多元素可能彼此非常远。单链接群集的另一个相关缺点是很难从数据分区的角度进行解释。这是单链接集群的一个真正的缺点，因为我们对数据分区感兴趣。

> cluster_single <- hclust (bc_dist, method = 'single')

> plot(cluster_single)

Complete Linkage Agglomerative Clustering：完全链接群集也称为最远邻居群集。它允许一个样本(或一组)仅在彼此最远的距离处与另一样本(或组)凝聚。因此，两个组的所有成员都链接在一起。完全连锁聚类避免了单一连锁方法链接样本的缺点。完全连锁往往会找到许多单独的小群体。

> cluster_complete <- hclust (bc_dist, method = 'complete')

> plot(cluster_complete)

Average Linkage Agglomerative Clustering：平均链接聚类允许样本在该样本与该簇的所有成员之间的距离的平均值处被分组到该簇。两个群集以一个群集的所有成员与另一个群集的所有成员之间的距离的平均值连接。产生树状图的聚类看起来有点介于单个和完全连锁聚类之间。

> cluster_average <- hclust (bc_dist, method = 'average')

> plot(cluster_average)

Ward’s Minimum Variance Clustering：Ward的最小方差聚类(又名Ward‘s聚类)最初由Ward(1963)提出，它是基于最小二乘线性模型准则的。该方法最小化样本之间距离的平方和(即方差分析的平方误差)的簇内和。基本上，它将聚类分析视为方差问题的分析，而不是使用距离度量或关联度量。这种方法最适合于数量变量，而不是二元变量。Ward的聚类是通过找出两个子集的对来实现的在导致合并后总的簇内方差最小增加的每个步骤上的簇。该增加是群集中心之间的加权平方距离。虽然初始聚类距离被定义为Ward聚类中点之间的平方欧几里得距离，但其他距离方法也可以产生有意义的结果。

> cluster_ward <- hclust (bc_dist, method = 'ward.D2')

> plot(cluster_ward)

使用par (mfrow = c(2,2))将结果一起绘制到一个图表中，以创建一个包含两行和两个面板的图表。

> par (mfrow = c(2,2))

> plot(cluster_single)

> plot(cluster_complete)

> plot(cluster_average)

> plot(cluster_ward)

恢复默认窗口面板。

> par (mfrow = c(1,1))

这四个数据图之间的比较表明，完全链接和Ward‘s聚类生成了相同的分区簇，并且在按照我们感兴趣的方向划分数据方面具有更好的性能(样本11、12、13、14和15来自vdr−/−小鼠，而22、23和24来自野生型小鼠)。

7.4  排序

排序的主要目标被认为是“探索性的”，随着规范对应分析(CCA)的引入，排序已经超越了单纯的“探索性”分析，也变成了假设检验。两个变量的数据结构通常通过样本的散点图来揭示。微生物组数据是多维的，一般有两个以上的变量。微生物组数据集可以被视为位于每个变量或物种(或OTUS/分类群)定义一个维度的空间中的样本(对象)的集合。因此，维度与变量或物种(或OTU/分类群)一样多。例如，在我们的VDR−/−小鼠粪便数据集中，有8个样本，248属。因此，数据维度是248。对于具有n个变量的数据集，需要绘制的散点图数量为n(n−1)/2。在我们的例子中，248个属将需要(248247)/2=30,628个散点图。如此大量的散点图对于了解数据结构不是很有意义，而且工作起来也很繁琐。排序主要致力于在低维空间中尽可能忠实地表示样本和物种(或OTUS/分类群)关系。这一目标是可取的，因为群落数据是混合了噪声的多维数据，低维可能理想且通常代表对物种(或OTUS/分类群)-环境关系的重要且直观的解释。为了一个nXP包含n个主题和p个变量的p个数据集，则可以将数据结构视为p维空间中的n个主题(表示为点的群集)。排序的主要目的是在减少数量的正交(即，独立)轴上表示多个样本(对象)，其中轴的总数小于或等于样本(对象)的数量。排序轴的重要性按顺序递减。排序的第一个轴解释了数据集中的最大变化，然后是第二个轴，然后是第三个轴，依此类推。排序图对于可视化样本(对象)之间的相似性特别有用。例如，在贝塔多样性的背景下，排序空间中距离较近的样本的物种组合比排序空间中相距较远的样本具有更相似的物种组合。排序方法包括约束排序和非约束排序，具体取决于排序轴是否受环境因素(变量)的约束。顾名思义，在约束排序中，排序轴受环境因素的约束：样本在排序中的位置受环境变量的约束。在无约束排序中，排序轴不受环境因素的约束。换言之，约束排序在构建排序时直接使用环境变量；而在非约束排序中，环境变量被输入到事后分析中。从假设检验的角度来看，无约束排序分析本身主要是描述性的方法，并不真正涉及假设检验。多变量数据。虽然它涉及使用环境变量来解释轴的假设，例如在vegan包中使用函数envfit()，但无约束排序很简单。它分析一个数据矩阵；它的目标是从一组减少的正交轴中揭示图中的主要数据结构。相反，约束排序是一种“假设驱动”排序：一种假设检验方法，直接检验关于环境因素对物种(或OTUS/分类群)组成的影响的假设。约束排序与多元线性模型的关系是这样的：左侧的“因变量”(或社区变量)作为响应，右侧的“独立”变量(或约束)作为解释因素。因此，约束排序是不对称的。

Principal Component Analysis (PCA)：就vegan而言，在我们的案例中，环境变量是遗传条件(VDR、−/−和WT)。我们想知道VDR基因缺失是否解释了属组成中的β多样性。我们进行主成分分析，探讨群落属组成(β多样性)的变化是否由遗传因素引起。通过主成分分析，基于轴1上的属A、轴2上的属B、轴3上的属C的丰度来绘制样本，直到在非常高维空间中绘制了n个样本。N个样本创建(n−1)个PC：PC1是穿过所有这些样本创建的空间的第一条直线，PC2是垂直于PC1的第二条直线，因此是第三条PC3，直到PC(n−1)。个人电脑的重要性按顺序递减。PC1是最重要的PC，也是所有样本中变化最大的。PC2解释了第二多的变异，PC3解释了第三多的变异，依此类推，(n-1)PC解释了样本的最小变异。可以使用几个R函数(包括预安装的stats包中的prcomp()、vegan包中的rda()、labdsv包中的PCA())来执行PCA。在这里，我们使用函数rda()。两个扩展函数：evlot()和BioDiversityR软件包中的PCAsignance()可以改善地块。evlot()通过使用Keiser-Guttman准则和折杆模型，提供了确定排序轴重要性的直观方法。函数的作用是：计算PCA轴的折线模型。当使用函数rda()执行PCA时，通过不指定环境数据矩阵(即，group variable)，该函数执行无约束排序PCA。在微生物组数据分析中，绝对丰度计数不合适，因为过大的数值会对分析产生太大的影响。因此，我们需要对大量的读取数据进行标准化处理，然后才能进行分析。在这里，我们使用带有total方法的函数decstand()来标准化读取。

> stand_abund_table <- decostand(abund_table, method = "total")

> PCA <-rda(stand_abund_table)

> PCA

Call: rda(X = stand_abund_table)

Inertia Rank

Total          0.0408   

Unconstrained  0.0408    7

Inertia is variance

Eigenvalues for unconstrained axes:

PC1    PC2    PC3    PC4    PC5    PC6    PC7

0.02029 0.01279 0.00371 0.00173 0.00143 0.00077 0.00010

在vegan的语言中，“Inertia"是数据中“variation”的总称。在这种情况下，整个数据集的总变异为0.0408，第一个轴解释了总变异的60.4%(0.02029/0.0408=0.4973)。总变异是分析矩阵中每个属的变异之和。以下R代码检查总方差：

> sum (apply (stand_abund_table, 2, var))

[1] 0.04082

然后，让我们使用函数bilot()绘制图表。显示选项“species”是OTUS/TATA的vegan包标签。默认值为“sites”(样品标签)(图7.14)。

> biplot(PCA, display = 'species')

上面由bilot()绘制的图表只是绘制属的箭头，这并不能提供信息。更具信息量的情节是使用函数ordilot()将属和样本分数都绘制为质心，如下所示(图7.15)：

> ordiplot(PCA, display = "sites", type = "text")

在上述增强中，type=“text”或“t”将文本标签添加到地物(默认设置仅添加点)。作为替代方案，我们可以使用François Gillet和DanielBorcard编写的函数leanplot.pca()将绘制PCA结果绘制到两个图中，并且根据比例的不同而有所不同。此函数用于从PCA中“rda”类的对象或vegan的rda()函数的RDA结果中绘制两个双线图(比例1和比例2)。这在Borcard等人的书“Numerical Ecology with R”中提供了函数。我们首先运行函数leanplot.pca()，如下所示。

"cleanplot.pca" <- function(res.pca, ax1=1, ax2=2, point=FALSE,

ahead=0.07, cex=0.7)

{

# A function to draw two biplots (scaling 1 and scaling 2) from an object

# of class "rda" (PCA or RDA result from vegan's rda() function)

#

# License: GPL-2

# Authors: Francois Gillet & Daniel Borcard, 24 August 2012

require("vegan")

par(mfrow=c(1,2))

p <- length(res.pca$CA$eig)

# Scaling 1: "species" scores scaled to relative eigenvalues

sit.sc1 <- scores(res.pca, display="wa", scaling=1, choices=c(1:p))

spe.sc1 <- scores(res.pca, display="sp", scaling=1, choices=c(1:p))

plot(res.pca, choices=c(ax1, ax2), display=c("wa", "sp"), type="n",

main="PCA - scaling 1", scaling=1)

if (point)

{

points(sit.sc1[,ax1], sit.sc1[,ax2], pch=20)

text(res.pca, display="wa", choices=c(ax1, ax2), cex=cex,

pos=3, scaling=1)

}

else

{

text(res.pca, display="wa", choices=c(ax1, ax2), cex=cex,

scaling=1)

}

text(res.pca, display="sp", choices=c(ax1, ax2), cex=cex, pos=4,

col="red", scaling=1)

arrows(0, 0, spe.sc1[,ax1], spe.sc1[,ax2], length=ahead, angle=20,

col="red")

pcacircle(res.pca)

# Scaling 2: site scores scaled to relative eigenvalues

sit.sc2 <- scores(res.pca, display="wa", choices=c(1:p))

spe.sc2 <- scores(res.pca, display="sp", choices=c(1:p))

plot(res.pca, choices=c(ax1,ax2), display=c("wa","sp"), type="n",

main="PCA - scaling 2")

if (point) {

points(sit.sc2[,ax1], sit.sc2[,ax2], pch=20)

text(res.pca, display="wa", choices=c(ax1 ,ax2), cex=cex, pos=3)

}

else

{

text(res.pca, display="wa", choices=c(ax1, ax2), cex=cex)

}

text(res.pca, display="sp", choices=c(ax1, ax2), cex=cex, pos=4,

col="red")

arrows(0, 0, spe.sc2[,ax1], spe.sc2[,ax2], length=ahead,

angle=20, col="red")

}

"pcacircle" <- function (pca)

{

# Draws a circle of equilibrium contribution on a PCA plot

# generated from a vegan analysis.

# vegan uses special constants for its outputs, hence

# the 'const' valuebelow.

eigenv <- pca$CA$eig

p <- length(eigenv)

n <- nrow(pca$CA$u)

tot <- sum(eigenv)

const <- ((n - 1) * tot)^0.25

radius <- (2/p)^0.5

radius <- radius * const

symbols(0, 0, circles=radius, inches=FALSE, add=TRUE, fg=2)

}

然后调用PCA，如下所示(图7.16)。

> cleanplot.pca (PCA)

“Scaling”是将排序结果投影到减少的图形显示空间中的方式。函数leanplot.pca()生成两个PCA缩放。左图中的PCA-Scaling 1是距离双线图，它关注的是样本之间的距离。如果您主要想了解样本之间的关系，请选择Scaling 1，Scaling 1的基本功能如下：特征向量按单位长度缩放，双曲线图中样本(对象)之间的距离是它们的欧几里得距离的近似值多维空间。然而，变量(在这种情况下是细菌)向量之间的夹角是没有意义的。这个圆叫做平衡贡献圆，代表变量的平衡贡献。对于给定的轴组合，向量大于圆周半径的变量可以有把握地解释为最重要的细菌，而向量小于平衡贡献圆半径的变量对给定的约化空间贡献很小。右图中的PCA-Scaling 2是相关双线图。它绘制了变量(在本例中为细菌)之间的相关性。如果您的主要兴趣集中在变量之间的关系上，请选择Scaling 2。Scaling 2的基本功能是：每个特征向量都被缩放到其特征值的平方根；向量的长度近似变量(细菌)的标准差；变量(在这种情况下是细菌)之间的角度反映了它们的相关性：角度的余弦近似变量(细菌)之间的相关性；然而，双线图中样本(对象)之间的距离不是它们在多维空间中的欧几里得距离的近似。

Principal Coordinate Analysis (PCoA):PCoA也称为度量多维缩放。PCoA是一种灵活的排序技术，其允许用户选择几乎任何距离度量(例如，Jaccard、Bray-Curtis、欧几里德等)。与PCA一样，PCoA使用特征值来衡量一组返回的正交轴的重要性。通过确定每个特征向量和特征值，对矩阵进行降维。通过缩放每个特征向量来获得主坐标。当根据样本之间的欧几里得距离计算PCoA时，产生的结果与基于相同数据集的协方差矩阵计算的PCA相同(如果使用1尺度)。R函数，包括vegan包中的cmdscale()和APE包中的pcoA()，可以执行pcoA。对于vegan，输入数据可以通过函数vegdist()计算(default为Bray-Curtis Disimility)，排序图可以使用函数ordilot()绘制。也可以使用APE包中的函数biplot.pcoa()绘制排序图。在本例中，我们使用cmdscale()函数和相同的vdr−/−小鼠排泄物数据来执行pCoA。该函数需要一个相似矩阵作为输入数据。首先，让我们使用函数vegdist()计算Bray-Curtis相异度，并将其命名为：

> bc_dist <-vegdist(abund_table, "bray")

然后，我们将显式设置k=2(要返回的维数的默认值)和eig = TRUE (保存特征值)。

> PCoA <- cmdscale (bc_dist, eig = TRUE,k = 2)

> PCoA

$points

[,1]      [,2]

5_15_drySt-28F  0.35060  0.007035

1_11_drySt-28F -0.17119  0.137244

2_12_drySt-28F  0.02928  0.115931

3_13_drySt-28F -0.17013  0.126964

4_14_drySt-28F  0.27190  0.088394

7_22_drySt-28F  0.13943 -0.258936

8_23_drySt-28F -0.25091 -0.157891

9_24_drySt-28F -0.19898 -0.058740

$eig

[1]  3.779e-01  1.517e-01 1.170e-01  8.855e-02  2.771e-02

[6]  2.167e-02 -2.515e-17 -4.493e-03

$x

NULL

$ac

[1] 0

$GOF

[1] 0.6713 0.6751

函数cmdscale()的作用是：生成输出列表。第一个输出点包含每个降维中每个样本的坐标。第二个输出EIG包含特征值。最后三个输出与分析的其他选项有关，我们在此不会涉及这些选项。下面的R码块用于检查由PCoA的前两个轴解释的数据集中的百分比变化。

> explainedvar1 <- round(PCoA$eig[1] / sum(PCoA$eig), 2) * 100

> explainedvar1

[1] 48

> explainedvar2 <- round(PCoA$eig[2] / sum(PCoA$eig), 2) * 100

> explainedvar2

[1]19

> sum_eig <- sum(explainedvar1, explainedvar2)

> sum_eig

[1] 67

第一个轴解释了48%的数据变化，第二个轴解释了19%的数据变化。因此，这两个轴解释了数据的大量变化(总计67%)。有两个标准来评估前几个PCoA轴是否捕获了不成比例的大量总解释变异：(1) Kaiser-Guttman criterion and (2) broken-stick model。Kaiser-Guttman准则规定，“与前几个轴相关的特征值应该大于所有特征值的平均值。”根据折杆模型的判据，将与前几个轴相关的特征值与broken-stick model模型的期望值进行了比较。broken-stick model假设特征值的总和随着PCoA轴的有序顺序递减。我们使用这两个标准通过以下曲线图评估PCoA的性能。

> # D efine Plot Parameters

> par(mar = c(5, 5, 1, 2) + 0.1)

> # Plot Eigenvalues

> plot(PCoA$eig, xlab = "PCoA", ylab = "Eigenvalue",  las = 1, cex.lab = 1.5, pch = 16)

> # Add Expectation based on Kaiser-Guttman criterion and Broken Stick Model

> abline(h = mean(PCoA$eig), lty = 2, lwd = 2, col = "blue")

> b_stick <- bstick(8, sum(PCoA$eig))

> lines(1:8, b_stick, type = "l", lty = 4, lwd = 2, col = "red")

> # Add Legend

> legend(“topright", legend = c("Avg Eigenvalue", "Broken-Stick"), lty = c(2, 4), bty = "n", col = c("blue", "red"))

图7.17显示，与前三个轴相关的特征值大于所有特征值的平均值，也大于折杆模型的预期。在评估PCoA输出之后，接下来我们将为两个PCoA轴创建排序图(图7.18)。

> # D efine Plot Parameters

> par(mar = c(5, 5, 1, 2) + 0.1)

> # Initiate Plot

> plot(PCoA$points[ ,1], PCoA$points[ ,2], ylim = c(-0.5, 0.5),  xlab = paste("PCoA 1 (", explainedvar1, "%)", sep = ""), ylab = paste("PCoA 2 (", explainedvar2, "%)", sep = ""), pch = 5, cex = 1.0, type = "n", cex.lab = 1.0, cex.axis = 1.2,axes = FALSE)

> # Add Axes

> axis(side = 1, labels = T, lwd.ticks = 2, cex.axis = 1.2, las = 1)

> axis(side = 2, labels = T, lwd.ticks = 2, cex.axis = 1.2, las = 1)

> abline(h = 0, v = 0, lty = 3)

> box(lwd = 2)

> # Add Points & Labels

> points(PCoA$points[ ,1], PCoA$points[ ,2],  pch = 19, cex = 3, bg = "blue", col = "blue")

> text(PCoA$points[ ,1], PCoA$points[ ,2],  labels = row.names(PCoA$points))

基本排序图使我们可以看到样本是如何彼此分开的。在我们的例子中，由于不同老鼠属的丰度不同，样品沿着PCoA轴分离。一个合乎逻辑的后续问题是询问是什么类型的数据集驱动了观察到的点之间的差异。我们能识别并可视化PCoA中这些有影响力的属吗？我们可以使用BioDiversityR包中的add.spec.core()函数获得此信息。首先，相对丰度计算如下：

> fecalREL <- abund_table

> for(i in 1:nrow(abund_table)){ fecalREL[i, ] = abund_table[i, ] / sum(abund_table[i, ]) }

然后计算属级分数并加到图中(图7.19)

> require("BiodiversityR")

> PCoA <- add.spec.scores(PCoA,fecalREL,method = "pcoa.scores")

> text(PCoA$cproj[ ,1], PCoA $cproj[ ,2],   labels = row.names(PCoA$cproj), col = "black")

函数add.spec.scores()的作用是：确定沿着PCoA轴的每个毒素之间的相关性。这是识别有影响力的类群的一种更有效的定量方法。

> Genus_corr <- add.spec.scores(PCoA, cecalREL, method = "cor.scores")$cproj

为了识别和提取重要的分类单元，我们需要定义一个相关系数截止值，例如0.70。然后从任一维打印出相关系数大于0.70的分类单元。

> corrcut <- 0.7

> genus_corr <- add.spec.scores(PCoA, fecalREL, method = "cor.scores")$cproj

> import_genus <- genus_corr[abs(genus_corr[, 1]) >= corrcut |abs(genus_corr[, 2]) >= corrcut, ]

沿着PCoA轴具有大于或等于0.7的相关性的12个重要属输出如下：

> import_genus[complete.cases(import_genus),]

Dim1    Dim2

Tannerella            0.82046  0.16847

Lactobacillus        -0.30640 -0.86651

Helicobacter          0.74876  0.11758

Paraprevotella        0.74631 -0.08582

Bacillus              0.19803 -0.77265

Pedobacter          -0.24607  0.87028

Limibacter          -0.08244  0.70361

Mycoplasma            0.72268  0.22945

Slackia              0.09773 -0.75775

Fluviicola            0.24247 -0.71075

Caldicellulosiruptor  0.24247 -0.71075

Anaerophaga          0.24247 -0.71075

最后，我们使用vegan包中的envfit()函数对这些相关性的轴上的一般丰度进行排列测试。

> envfit(PCoA, fecalREL, perm = 999)

部分输出如下：

Tannerella 0.951 0.308 0.70 0.044 *

Lactobacillus          -0.219 -0.976    0.84  0.013 *

Parasutterella          -0.754 -0.657    0.63  0.031 *

Porphyromonas            0.359 -0.933    0.56  0.072 .

Bacillus                0.160 -0.987    0.64  0.056 .

Pedobacter              -0.176  0.984    0.82  0.015 *

Slackia                  0.081 -0.997    0.58  0.085 .

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Permutation: free

Number of permutations: 999

q值小于0.05的属如下图所示(图7.20)：

> fit <- envfit(PCoA, fecalREL, perm = 999)

> plot(fit, p.max = 0.05, col = "red")

Non-metric Multidimensional Scaling (NMDS)：NMDS是PCoA分析的非度量替代方案。NMDS已经被认为是一种很好的排序方法，因为它使用具有生态学意义的方法来衡量群落差异和样本之间的任何距离(不同)作为输入。因此，这是社区排序中推荐的方法。NMDS分析的重点是将样本点的相对位置投影到低维排序空间(两个或三个轴)。vegan包中的函数metaMDS()执行NMDS分析。为简化起见，NMDS分析算法总结如下：首先，它使用vegdist()函数获取足够的不同度量；然后，它使用随机启动配置运行NMDS多次，并通过函数Procrstes()比较结果，然后在找到两次类似的最小应力解决方案后停止。最后，它缩放和旋转解决方案，并使用函数wacores()将物种(或OTUS/分类群)得分作为加权平均值添加到配置中。算法完成后，使用PCA旋转最终解以简化其解释。在这种情况下，我们使用vegan和相同的VDR−/−小鼠的粪便数据来说明NMDS分析。首先，我们使用默认设置的metaMDS()函数使用vegdist()函数获得Bray-Curtis相异度度量(默认方法)。它自动转换数据并检查解决方案的健壮性。在metaMDS()函数调用中结合使用了威斯康星双重标准化和SQRT变换。函数调用产生7.57%的应力值。

> bc_nmds <- metaMDS(abund_table, dist = "bray")

Square root transformation

Wisconsin double standardization

Run 0 stress 0.07574

Run 1 stress 0.1562

Run 2 stress 0.1284

Run 3 stress 0.1801

Run 4 stress 0.1284

Run 5 stress 0.1562

Run 6 stress 0.1678

Run 7 stress 0.1393

Run 8 stress 0.1562

Run 9 stress 0.1284

Run 10 stress 0.1393

Run 11 stress 0.1192

Run 12 stress 0.1764

Run 13 stress 0.1192

Run 14 stress 0.1551

Run 15 stress 0.1562

Run 16 stress 0.1562

Run 17 stress 0.07574

… Procrustes: rmse 1.314e-06 max resid 2.644e-06

… Similar to previous best

Run 18 stress 0.1284

Run 19 stress 0.1899

Run 20 stress 0.07574

… Procrustes: rmse 1.503e-06 max resid 2.229e-06

… Similar to previous best

*** Solution reached

> bc_nmds

Call:

metaMDS(comm = abund_table, distance = "bray")

global Multidimensional Scaling using monoMDS

Data:  wisconsin(sqrt(abund_table))

Distance: bray

Dimensions: 2

Stress:  0.07574

Stress type 1, weak ties

Two convergent solutions found after 20 tries

Scaling: centring, PC rotation, halfchange scaling

Species: expanded scores based on 'wisconsin(sqrt(abund_table))'

其次，我们使用ordilot()函数来绘制NMDS结果。默认设置仅将点添加到图中，我们在此图中使用type=‘t’或type=‘text’添加文本标签(图7.21)。

> ordiplot (bc_nmds, type = 't')

要将现场(样本)分数绘制为文本(图7.22)：

> ordiplot(bc_nmds, display = "sites", type = "text")

最后，使用函数StressPlot()绘制Shepards应力曲线图。函数压力图()生成两个图：一个图显示排序距离与观察到的差异(所选群落差异)的关系，并用单调的步进线表示拟合；另一个图绘制拟合优度图，以评估特定样本的NMDS2与NMDS1的排序优度。下列R代码用于将绘图窗口分成两个面板：

> par (mfrow = c(1,2))

函数plot()绘制具有站点(样本)的NMDS排序图：

> stressplot (bc_nmds)

> plot (bc_nmds, display = 'sites', type = 't', main = 'Goodness of fit')

points()函数将反映拟合良好的大小的点相加(较大=拟合较差)(图7.23)。

> points (bc_nmds, display = 'sites', cex = goodness (bc_nmds)*300))

stressplo显示了所得到的m维排序解中样品之间的实际距离与所选的Bray-Curtis相异度所表示的特殊成分相异度之间的关系。拟合优度有两个类相关统计；基于应力的相关性：R2=1−S2(非量度拟合=0.994)以及拟合值与排序距离之间或阶跃线与点之间的相关性：“基于拟合的R2”(线性拟合=0.956)。

Correspondence Analysis (CA)：CA是另一种无约束排序方法。它使用所有排序轴的多维空间中样本之间的卡方距离，并对稀有物种(例如，出现频率较低的多个零点的物种)赋予较高的权重。与PCA和PCoA一样，CA使用特征值来衡量返回的正交轴的重要性。在主成分分析排序图中，分类群是矢量，样本是点，而在CA排序图中，分类群和样本是由点表示的。与PCA类似，CA有两种尺度：尺度1和尺度2。在简化的排序空间中，样本之间的距离(尺度1)近似于它们的卡方距离。例如，任何代表分类单元的点附近的样本都可能包含对该分类单元的高贡献。分类群之间的距离(比例2)也近似于它们的卡方距离。例如，任何靠近代表样本的点的分类单元在该样本中更有可能出现更高的频率。vegan包中的函数cca()可用于执行无约束对应分析。基于Kaiser-Guttman准则或折杆模型，我们使用Cca()函数进行CCA分析，使用evlot()函数选择重要的排序轴。由Borcard等人编写的函数evlot()。这里用来绘制排序对象的特征值和变异百分比。当函数cca()(cca=规范成分分析)用于执行无约束CA时，请不要指定环境矩阵或分组信息。

> fecal_genus_cca=cca(abund_table)

> fecal_genus_cca

Call: cca(X = abund_table)

Inertia Rank

Total          0.502   

Unconstrained  0.502    7

Inertia is mean squared contingency coefficient

122 species (variables) deleted due to missingness

Eigenvalues for unconstrained axes:

CA1    CA2    CA3    CA4    CA5    CA6    CA7

0.2036 0.1256 0.0747 0.0374 0.0322 0.0200 0.0090

数据的总异质性(inertia)为0.502，第一轴捕获了属组成总变异的40.56%(0.2036/0.502=0.4056，其中0.2036是第一轴CA1的特征值，0.502是数据的总异质性)。以下R代码用于绘制排序并在图中显示样品名称(图7.24)：

> plot(fecal_genus_cca, display="sites")

如果您不希望图形过于拥挤，请使用以下R代码仅显示点，而不是图中的示例名称：

> plot(fecal_genus_cca, display="sites", type="p")

下面的排序图显示了排序图中的样本和属的格局(图7.25)：

> ordiplot (fecal_genus_cca)

下面使用函数evlot()进行的CA后分析是为了说明如何确定应该使用哪个CA轴来解释结果。函数evlot()的语法是evlot(Ev)，其中ev是特征值的向量。首先，我们运行函数。

evplot <- function(ev)

{# Broken stick model

n <- length(ev)

bsm <- data.frame(j=seq(1:n), p=0)

bsm$p[1] <- 1/n

for (i in 2:n) bsm$p[i] <- bsm$p[i-1] + (1/(n + 1 - i ) )

bsm$p <- 100*bsm$p/n

# Plot eigenvalues and % of variation for each axis

op <- par(mfrow=c(2,1))

barplot(ev, main="Eigenvalues", col=“bisque", las=2)

abline(h=mean(ev), col="red")

legend("topright", "Average eigenvalue", lwd=1, col=2, bty="n")

barplot(t(cbind(100*ev/sum(ev), bsm$p[n:1])), beside=TRUE,

main="% variation", col=c("bisque",2), las=2)

legend("topright", c("% eigenvalue", "Broken stick model"),

pch=15, col=c("bisque",2), bty="n")

par(op)

}

# Plot eigenvalues and % of variance for each axis

ev <- fecal_genus_cca$CA$eig

windows(title="CA eigenvalues")

evplot(ev)

如PCoA部分所述，为了评估分析的前几个CA轴是否捕获了不成比例的大量总解释变异，使用图形的另外两个后处理方法是 Kaiser-Guttman criterion and  broken-stick model。根据Kaiser-Guttman criterion，与前几个轴相关的特征值应大于所有特征值的平均值；对于折杆模型，将与前几个轴相关的特征值与期望值进行比较。Keiser-Guttman准则和broken-stick model都表明前三个轴是重要的(图7.26)。