肠型在这里

内容基本上翻译自：https://enterotype.embl.de/index.html

介绍

肠型是一种根据人类肠道微生物群来分层的人的方法。它们来自个体肠道微生物群中不同微生物群的相对丰度（目前处于属级;请参阅维基百科关于"相对物种丰度"的文章，以了解相对丰度的概念）。肠型不是离散类型，如血型，而是在微生物群特征的更高维度空间中稠密的区域。我们假设稠密的区中，个体的很大一部分，它们之间或周围有人烟稀少的区域，但这种分布的原因尚不清楚，也没有公布足够的数据来量化，随着时间的推移，肠道类型的分类和可能开关的鲁棒性。

2011年4月，MetaHIT公布了人类肠道微生物群中肠型的发现。与研究相关的数据被公开，计算过程背后的理论在文章的补充资料中作了解释。然而，确切的命令（在R环境中），使任何人都可以重现文章中的所有数字没有报告在补充。在这里，我们介绍了一个详细的教程，以复制我们在原始文章中的工作，使用R与原始文章中使用的数据集。

可以在此处下载重现本教程结果的 R 代码摘要，但我们强烈建议完成它以及下面提供的详细信息。

数据

本文中有三个数据集。主要结论来自基于Sanger的宏基因组序列，我们在本教程中使用该数据集。处理其他数据集的步骤相同，除非另有提及。注：初始数据集包含 39 个示例。四个日本婴儿样本被排除在外，因为已知婴儿含有非常异质、不稳定和独特的微生物群。两个美国成人样本被排除在外，因为细菌类和疑似技术文物水平异常低。有关详细信息，请参阅文章及其参考。

首先，下载包含 33 个样本的 Sanger 数据集的属丰度表，并将其加载为 R 环境中的表：

data=read.table("MetaHIT_SangerSamples.genus.txt", header=T, row.names=1, dec=".", sep="\t") #读入文件
data=data[-1,] #去除第一行表头

了解属丰度表

属丰度表总结了33个样本的相对属丰度。此表采用 R 表格式 - 这意味着第一行是具有 33 个示例标识符的列标题，以下行为 34 列，每列为属标识符。两个基于基因组的数据集中的第二行，具有属标识符"-1"，对应于无法可靠地分配给已知属的基因组读取分数。在观察到可以可靠地分配的读取部分样本之间存在较大差异（例如，在上表中，19% 到 78% 的读取未分配），我们引入了此未分配分数（而不是忽略所有未可靠地分配给已知属的读取）。对未分配分数进行核算可确保可以更可靠地比较多个样本中给定属的相对丰度。然而，由于此分数不代表真实属，并且是所有属的集团，我们至今没有代表，因此，在估算聚类过程的 Jensen-Shannon 距离 （D） 时，我们不将其视为特征，因此，未分配分数的类似水平不会误认为是真正的遗传相似性。这意味着相对丰度表不会重新规范化 - 当包含"未分配"分数时，它仍然求和为 1，但"未分配"分数不用于计算 Jensen-Shannon 距离。

聚 类

使用JSD距离和Medoids（PAM）聚类算法，根据相对属丰度对样本进行聚类。使用卡林斯基-哈拉巴斯（CH） 指数 评估了星团的最佳数量的结果。我们还通过比较最优聚类的轮廓系数（Rousseeuw，1984年）与从模拟派生的轮廓系数的分布进行比较，从而评估了最优聚类的统计意义（请参阅下面的模拟数据部分）。

距离指标

属丰度剖面（遗传）和OG丰度剖面（功能）被归一化，以产生概率分布（以下简称丰度分布）。我们使用与Jensen-Shannon 区别 （JSD）相关的概率分布距离指标对样本进行分组。Jensen-Shannon 区别是一个积极的确定尺度，满足以下条件：

JSD（a，b）>=0
JSD（a，b）=0如果且仅当a=b

它还是对称的：

JSD（a，b）=JSD（b，a）

但是，它不符合三角不平等条件，被认为是一个真正的指标。

JSD（a，c）=JSD（a，b）=JSD（b，c）

另一方面，JSD的平方根是一个真正的距离指标（Endres & Schindelin，2003年;在Arumugam等人2011年的补充中引证为参考文献59）。这就是我们使用平方根而不是JSD本身的原因。

JSD<- function(x,y) sqrt(0.5 * KLD(x, (x+y)/2) + 0.5 * KLD(y, (x+y)/2))
#我们在丰度分布中添加了 0.000001 的伪计数，以避免公式的分子和/或分母中的零。
KLD <- function(x,y) sum(x * log(x/y))

下面是在 R 中创建 JSD 距离矩阵的完整函数示例：

#Function 1
dist.JSD <- function(inMatrix, pseudocount=0.000001, ...) {
	KLD <- function(x,y) sum(x *log(x/y))
	JSD<- function(x,y) sqrt(0.5 * KLD(x, (x+y)/2) + 0.5 * KLD(y, (x+y)/2))
	matrixColSize <- length(colnames(inMatrix))
	matrixRowSize <- length(rownames(inMatrix))
	colnames <- colnames(inMatrix)
	resultsMatrix <- matrix(0, matrixColSize, matrixColSize)
        
  inMatrix = apply(inMatrix,1:2,function(x) ifelse (x==0,pseudocount,x))

	for(i in 1:matrixColSize) {
		for(j in 1:matrixColSize) { 
			resultsMatrix[i,j]=JSD(as.vector(inMatrix[,i]),
			as.vector(inMatrix[,j]))
		}
	}
	colnames -> colnames(resultsMatrix) -> rownames(resultsMatrix)
	as.dist(resultsMatrix)->resultsMatrix
	attr(resultsMatrix, "method") <- "dist"
	return(resultsMatrix) 
 }

现在，我们可以在数据上应用此功能：

data.dist=dist.JSD(data)

算法

我们使用围绕中心点 （PAM） 聚类算法的来聚类丰度谱文件。PAM 源自基本 k-means 算法，但优点是它支持任何任意距离度量，并且比 k-means 更强劲。它是一个受监督的过程，其中预先确定的群集数作为过程的输入给出，然后将数据分区到这么多群集中。

在 R 中，我们在群集库中使用了pam（）函数

pam(as.dist(x), k, diss=TRUE) # x is a distance matrix and k the number of clusters

pam（）返回不同的对象，但我们只关注$clustering。
下面是一个在 R 中执行 PAM 聚类的示例函数：

 pam.clustering=function(x,k) { # x is a distance matrix and k the number of clusters
                         require(cluster)
                         cluster = as.vector(pam(as.dist(x), k, diss=TRUE)$clustering)
                         return(cluster)
                        }

在确定群集的最佳数量之前，以k=3群集为例，在我们的数据集上运行测试（请参阅下面的部分）：

data.cluster=pam.clustering(data.dist, k=3)

最佳群集数

为了评估数据集最牢固地分区到的群集的最佳数量，我们使用了 Calinski-Harabasz（CH） 索引，该指数在恢复群集数量方面表现出良好的性能。
其中 Bk是聚类之间的平方和 （即所有点i和j之间的平方距离 ，其中 i和j不在同一个聚类） 和Wk是聚类内平方和 （即所有点i和j之间的平方距离， 其中i和j位于同一群集中）。此度量实现了当群集之间的距离比群集内距离大得多时，聚类更加健壮。因此，我们选择了群集k的数量，因此CHk是最大化的。

在这里，我们使用函数索引。G1（）来自库集群Sim：

require(clusterSim)
nclusters = index.G1(t(data), data.cluster, d = data.dist, centrotypes = "medoids")

我们需要计算每个数量的群集k的CH索引，下面是执行该操作的示例代码（并绘制结果）：

nclusters=NULL

for (k in 1:20) { 
	if (k==1) {
		nclusters[k]=NA 
	} else {
		data.cluster_temp=pam.clustering(data.dist, k)
		nclusters[k]=index.G1(t(data),data.cluster_temp,  d = data.dist,
		centrotypes = "medoids")
	}
}

plot(nclusters, type="h", xlab="k clusters", ylab="CH index")

这表明，此特定数据集的最佳群集数为 3（k=3）。

  data.cluster=pam.clustering(data.dist, k=3)

群集证

聚类验证方法可用于评估与基础数据点相关群集的质量。在这里，我们使用轮廓验证技术。
(5)
a（i）是样本i与同一群集中所有其他样本的平均不同（或距离），而b（i）是最接近其他群集中所有对象的平均不同（或距离）。

公式表示-1 =< S（i） =< 1 。与其他群集更靠近其自身群集的样本具有较高的S（i）值，而S（i）接近 0 表示给定的样本位于两个群集之间的某处。较大的负S（i）值表示样本已分配给错误的群集。

为了评估轮廓系数 S（i），我们只需使用群集库中的函数轮廓（），该函数剪影（）用作聚类和距离矩阵对象的输入。要获得群集质量的全局评估，所有数据点上的平均S（i）是一个有用的度量。

obs.silhouette=mean(silhouette(data.cluster, data.dist)[,3])

模拟数据

对于一组给定的N特征载体，每个载体代表基因组样本的遗传（属、植激素等）或功能（基因、正交组、功能模块等）组成，我们模拟了包含相同数特征的N假设基因组样本，从连续体中抽样如下：

模板均匀：每个生成的要素对应于真实数据集中的要素，其跨多个样本的值在真实数据集中该要素的最小和最大丰度之间均匀分布
模板-高斯：每个生成的要素对应于真实数据集中的要素，其跨多个样本的值遵循高斯分布，其平均值和标准差与相应真实数据中观察到的平均值和标准差相同
生成的值然后规范化，使样本总和内的丰度为 1。

然后，对模拟数据进行聚类，并估计每个模拟的轮廓系数。然后，我们将这些与数据集最佳聚类的轮廓系数进行比较，以评估其统计显著性。但是，我们不能排除由于模型选择不当而产生具有统计显著性的结果的可能性，即使最佳聚类结构不能为我们的数据提供良好的模型。

图形解释

肠型可以通过许多不同的方式可视化。其中两种是类间分析（BCA） 和主坐标分析 （PCoA）。BCA 是Arumugam 等人文章中使用的第一个，可用于可视化群集的分类驱动因素（即肠型）。然而，PCoA在生态研究中更为常见，我们目前推荐了这种可视化技术，使用JSD距离矩阵作为输入。

类间分析

进行了类间分析（BCA）以支持聚类并标识肠型的驱动程序。分析使用Ade4包的R完成。在分析之前，在光明网数据集中，如果所有样本的平均丰度低于0.01%，则去除丰度极低的属体，以降低噪声。类间分析是主体元分析与工具变量的一个特定案例：在这里变量是一个定性因子（即肠型聚类）。类间分析使我们能够首先根据每个组的重心找到主要组件，以突出组之间的差异，然后将每个样本与其组联系起来。它是数据的受监督投影，其中预定义类之间的距离（在本例中为肠型）最大化。因此，此投影中的前两个组件不一定与解释数据中最高可变性的前两个组件相同。

下面是消除噪声的函数（即低丰盛的属）：

 noise.removal <- function(dataframe, percent=0.01, top=NULL){
	dataframe->Matrix
	bigones <- rowSums(Matrix)*100/(sum(rowSums(Matrix))) > percent 
	Matrix_1 <- Matrix[bigones,]
	print(percent)
	return(Matrix_1)
 }

我们建议将此功能应用于使用短测序技术生成的数据，如 Illumina 或 Solid。在原始手稿中，我们仅在 Solexa 数据上应用降噪，而不适用于桑格或热测序数据。噪声消除功能可应用如下：

data.denoized=noise.removal(data, percent=0.01)

最后，我们可以执行类间分析，并使用s.class（）函数绘制结果

 obs.pca=dudi.pca(data.frame(t(data)), scannf=F, nf=10)
 obs.bet=bca(obs.pca, fac=as.factor(data.cluster), scannf=F, nf=k-1) 
 s.class(obs.bet$ls, fac=as.factor(data.cluster), grid=F)

此图再现了原始文章中的图 2a（包括下文供比较），而不添加用于发布的修饰。聚类标签（1，2，3）由聚类过程随机分配，它们分别对应于输入型ET3、ET1和ET2。原始文章中已发布的图 2a（下图）是通过编辑原始绘图生成的：（1）添加样本的标签，（2）使用三个不同的符号代替默认点，（3）为不同的肠型对椭圆着色，（4）将每个肠型的驱动程序系投影到 BCA 上。

下面是上面的图，没有代表三个肠型和没有肠型质质连接器的椭圆。

s.class(obs.bet$ls, fac=as.factor(data.cluster), grid=F, cell=0, cstar=0, col=c(4,2,3))

有关ade4包的更多信息和教程，请参阅此处

主体坐标分析
下面是一个简单的代码段，用于使用pco（）执行PCoA，然后使用s.class（）进行绘图：

obs.pcoa=dudi.pco(data.dist, scannf=F, nf=3)
s.class(obs.pcoa$li, fac=as.factor(data.cluster), grid=F)

下面是相同的图，没有代表三个肠型的椭圆，没有肠型质质连接器。

s.class(obs.pcoa$li, fac=as.factor(data.cluster), grid=F, cell=0, cstar=0, col=c(3,2,4))

有关ade4包的详细信息和教程可在此处获取，而本教程中的所有代码都嵌入到BiotypeR中。