221.Beta多样性PCoA和NMDS排序

时间  2020-07-03 标签 221.beta beta 多样性 pcoa nmds 排序

Beta多样性与PCoA和NMDS排序

本节做者：文涛，刘永鑫php

版本1.0.4，更新日期：2020年6月27日html

本项目永久地址：https://github.com/YongxinLiu/MicrobiomeStatPlot ，本节目录 221BetaPCoA，包含R markdown(*.Rmd)、Word(*.docx)文档、测试数据和结果图表，欢迎广大同行帮忙审核校对、并提修改意见。提交反馈的三种方式：1. 公众号文章下方留言；2. 下载Word文档使用审阅模式修改和批注后，发送至微信(meta-genomics)或邮件(metagenome@126.com)；3. 在Github中的Rmd文档直接修改并提交Issue。审稿人请在创做者登记表 https://www.kdocs.cn/l/c7CGfv9Xc 中记录我的信息、时间和贡献，以避免专著发表时遗漏。react

做为高通量测序的表明之一，扩增子目早已成为表征微生物群落的最主流手段，在后续的数据处理，生物信息学分析中最基础也是最重要的分析就是群落多样性分析（alpha多样性和beta多样性）。今天咱们来学习的就是群落bate多样性分析中最重要的-非限制性排序分析。ios

基本概念

β多样性(Beta diversity)

β多样性又称生境间的多样性(between-habitat diversity)，是指沿环境梯度不一样生境群落之间物种组成的相异性或物种沿环境梯度的更替速率，用于研究群落之间的种多度关系，例如：物种更替或物种组成的差别。git

群落的Beta多样性分析包括非限制性排序（如PCoA，NMDS等）、层次聚类、限制性排序等，且均以群落类似或距离为基础计算。github

非限制性排序和层次聚类并是不独立的，下面这张图表示的就是非限制性排序和层次聚类的关系：express

图1. 排序与聚类的关系。a. 距离矩阵，b. 聚类结果，c. 排序结果，d. 聚类结果叠加排序结果。图片来源 http://mb3is.megx.net/gustame/home编程

将聚类分析的结果与诸如非度量多维尺度分析(Non-metric multidimensional scaling, NMDS)产生的排序之类的结果结合，在图中，聚类和NMDS结果叠加在左下方面板中。在此示例中，结果展示出现一致性：汇集在一块儿的对象也彼此接近。c#

类似性和距离

生态类似性（Ecological resemblance）以计算不一样样本群落组成类似程度或距离（相异程度）为基础，是处理多元生态数据的基本方法之一。在群落数据的分析中，经常使用其反映Beta多样性。微信

如在物种数据的分析中，对于两个群落，若它们共享相同的物种，而且全部物种的丰度也一致，那么这两个群落就具备最高的类似程度（或最低距离0）。生态学数据分析中的不少统计方法都以样方之间的类似性或距离为基础，例如上述提到的Beta多样性分析中的聚类、排序等，即便对于PCA实质上在计算时基于欧几里得（euclidean）距离考虑的。

若两个对象在各属性上越近似，那么它们的类似性就越高。对于群落数据，这些属性通常就是物种组成，或者环境属性等。一般使用物种组成数据，依据类似性指数（similarity indices）判断群落类似性，范围由0（两个群落不共享任何物种）到1（两个群落的物种类型和丰度彻底一致）。全部类似性指数都可以转换为距离指数，转化公式为“距离指数 = 1 – 类似性指数”的关系。

（1）能够转化为类似性指数的距离指数，例如定量数据的相异百分率（也称为Bray-Curtis距离）等。两者相互转换的公式一般表示为D=1-S或S=1-D，其中S是类似性指数，D为距离指数。

（2）没法转化为类似性指数的距离指数，例如欧几里得距离、卡方距离。

距离指数（distance indices）或称距离测度（distance measures），与类似性指数相反，距离数值越大代表群落间差别越大。存在多种距离类型，例如欧几里得（Euclidean）距离、Bray-Curtis距离、UniFrac距离等。对于物种组成数据，距离指数的最小值为0（两个群落的物种类型和丰度彻底一致），最大取值取决于距离类型和数据自己。

类似性或距离的衡量标准有不少种，Legendre于1998编写的“Numerical Ecology”一书中的“Ecological resemblance”章节内容列出大约30种方法，

常见的类似性/距离指数

Jaccard:

Jaccard类似性指数（Jaccard similarity index）将两个样方共享的物种数量（a）除以两个样方中出现的全部物种的总和（a + b + c，其中b和c是仅在第一个和第二个样方中出现的物种数量）。计算公式以下：

其中，y1j和y2j分别是对象1和2中元素j的数值。如果群落物种数据，那么y1j和y2j即分别是样方1和2中物种j的丰度。p是物种数（样方-物种矩阵中的物种数）

Bray-curtis距离(Bray-curtis distance):

Bray-curtis距离（Bray-curtis dissimilarity）其计算公式以下：

欧几里得距离(Euclidean distance):

欧几里得距离是多变量分析中常用的一种距离，如在线性排序方法PCoA、CCA。计算公式以下：

其中p是物种数（样方-物种矩阵中的物种数），y1j和y2j表示两个样方中对应的物种多度。

可是在物种数据的分析中，欧几里得距离却表现不佳。由于它将“双零”现象视做相同存在的方式处理，会缩小两个共享不多物种的群落之间的距离。双零”是指在计算群落类似性（或距离）时，所比较的两个样方中缺失某些物种的状况。具体在群落中一个物种在两个样方内同时缺失，并不能成为这两个样方具备组成类似的依据，由于引发缺失的缘由可能彻底不一样,其次在物种矩阵内，不可解释的双零的数量取决于物种的数量，所以也会随着检测到的稀有种数量的增长而显著增长。

若在群落距离计算过程当中使用欧几里得距离，能够先对原始物种数据进行数据转化（常见的如弦转化、Hellinger转化等），而后再使用转化后的数据计算欧几里得距离。尽管弦距离、Hellinger距离等然是对称指数的范畴，可是相较于使用原始物种丰度数据所得的欧几里得距离，弦距离、Hellinger距离的优点体如今存在距离的“上限”，下降了欧几里得距离对“物种丰度”的敏感性，有效减小了“双零”问题致使的偏差。可是咱们一般选择使用非对称的Bray-curtis距离等。

Bray-curtis距离的取值范围范围由0（两个群落的物种类型和丰度彻底一致）到1（两个群落不共享任何物种），所以它也能够直接经过“1 – 距离指数=类似性指数”转化为类似性指数（上文提到的“类似百分率”）。Bray-curtis距离适用于群落物种数据分析的缘由在于它是一个非对称指数，可有效忽略双零。

Unifrac距离(Unifrac distance):

Unifrac距离，它经常使用于微生物群落的数据中（例如，16S扩增子测序）。Bray-curtis距离仅考虑了物种的存在与否及其丰度，没有考虑物种之间的进化关系，距离0表示两个群落的物种组成结构彻底一致。在Unifrac距离中，除了关注考虑了物种的存在与否及其丰度外，还将物种之间的进化关系考虑在内，距离0更侧重于表示两个群落的进化分类彻底一致。

在16S扩增子测序中，根据16S序列组成构建OTUs进化树，OTUs之间存在进化上的联系，所以不一样OTUs之间的（系统发育）距离实际上有“亲远”之分。将系统发育树和OTUs丰度数据共同考虑到距离计算就是Unifrac距离。而其它非进化距离，忽略了OTUs之间的进化关系，认为OTUs间的关系平等。

Unifrac距离分为非加权Unifrac距离（Unweighted unifrac distance）和加权Unifrac距离（Weighted unifrac distance）。两种的主要区别是否考虑了物种的丰度。非加权Unifrac距离只考虑了物种有无的变化，不关注物种丰度，若两个微生物群落间存在的物种种类彻底一致，则距离为0；加权Unifrac距离同时考虑物种有无和物种丰度的变化，若两个微生物群落间存在的物种种类及丰度彻底一致，则距离为0。

关于Unifrac距离的计算方法，详见 http://scikit-bio.org/docs/latest/generated/skbio.diversity.beta.html。

排序

排序过程是将样品或物种排列在必定的空间，在一个低维空间中，使类似的样品或物种距离相近，相异的样品或物种距离较远。也就是说排序能够揭示微生物-环境间的生态关系，下降维数，减小坐标轴的数目，使排序轴可以反映必定的生态梯度。常见的方法有：PCA、PCoA、CA、DCA、NMDS、RDA、CCA等等。

PCoA排序

主坐标分析（PCoA；也称为度量多维标度）展现在低维欧氏空间中的对象间（非）类似性。PCoA不使用原始数据，而是使用（相异）类似度矩阵做为输入。
从概念上讲，它与主成分分析（PCA）和对应分析（CA）类似，后者分别保留对象之间的欧几里得距离和χ2（卡方）距离。可是，PCoA能够保留任何（距离）度量产生的距离，从而能够更灵活地处理复杂的生态学数据。

若是对解析理解有困难，能够结合下图理解。假如你是一本养花工具宣传册的摄影师，你正在拍摄一个水壶。水壶是三维的，可是照片是二维的，为了更全面的把水壶展现给客户，你须要从不一样角度拍几张图片。下图是你从水壶背面，正面，正上方，斜上方的照片。

图2. 主坐标分析的通俗解析。图片来源 https://blog.csdn.net/HLBoy_happy

咱们看到斜上方的照片能最好的展现咱们观察的特征。咱们的PCoA分析的第1/2主轴的结果就相似于此图。

PCoA和PCA的不一样之处：PCA是基于OTU 表两两样品间欧式距离计算，而PCoA是基于两两样品之间的任何一种距离距离计算，即有更多的选择，若是PCoA 也使用欧式距离，则PCA和PCoA的分析结果是同样的。

另外，PCoA是基于距离矩阵，它的排序的目的是将N个样品排列在必定的空间，使得样品间的空间差别与原始距离矩阵保持一致，这类排序方法也称做多维标定排序（Multi—dimensional scaling）。若是排序依赖于相异系数的数值，就叫有度量多维标定法（metric multi—dimensional scaling）因此说PCoA分析也叫有度量多维标定法；若是排序仅仅决定于相异系数的大小顺序（秩次排序），则称为无度量多维标定法（Non—Metric Multi—Dimensional Scaling；NMDS）。

NMDS 排序

非度量多维尺度法是一种将多维空间的研究对象（样本或变量）简化到低维空间进行定位、分析和归类，同时又保留对象间原始关系的数据分析方法。适用于没法得到研究对象间精确的类似性或相异性数据，仅能获得他们之间等级关系数据的情形。其基本特征是将对象间的类似性或相异性数据当作点间距离的单调函数，在保持原始数据次序(秩)关系的基础上，用新的相同次序的数据列替换原始数据进行度量型多维尺度分析。换句话说，当资料不适合直接进行变量型多维尺度分析时，对其进行变量变换，再采用变量型多维尺度分析，对原始资料而言，就称之为非度量型多维尺度分析。其特色是根据样品中包含的物种信息，以点的形式反映在多维空间上，而对不一样样品间的差别程度，则是经过点与点间的距离体现的，最终得到样品的空间定位点图。

NMDS过程是迭代的，而且分几个步骤进行：

  • 在多维空间中定义群落的原始位置；

  • 指定下降维度的数量（一般为2个维度）；

  • 二维构造样本的初始配置；

  • 该初始配置下的距离相对于观察到的（测量的）距离进行回归；

  • 根据回归肯定应力(stress)或二维构造与预测值之间的差别；

若是应力较高，则按减少应力的方向从新定位2维中的点，而后重复进行直到应力低于某个阈值。经验法则：应力<0.05可很好地表示尺寸减少，<0.1很是好，<0.2 还不错，而应力<0.3 有待提升。

附加说明：最终结果可能会因初始配置（一般是随机的）和迭代次数而有所不一样，所以建议屡次运行NMDS并尽量减下降应力值。

首先，NMDS须要距离矩阵或相异矩阵。原始欧几里得距离并非达到此目的的理想方法：它们对总丰度敏感，所以即便物种的标识不一样，也可能将具备类似数量物种的站点(site)视为类似物种。它们对物种的缺失也很敏感，所以能够将缺乏相同物种数的站点视为类似物种。

所以，生态学家使用Bray-Curtis相异性计算，该计算具备许多理想属性：

  • 它不随单位的变化而变化

  • 它不受添加/删除两个群落中不存在的物种的影响

  • 它不受新增群落的影响

  • 它能够识别总丰度的差别

实例分析

例1.两地点重复的两组PCoA

本文于2019年6月5日发表在Nature Biotechnology 杂志(37卷第6期)，并选为当期封面文章。点击查看中文解读

两图并列展现两组间明显的微生物组差别且在不一样地点可重复。不一样组采用着色配置信椭圆突出组间差别。

c.基于Bray-Cutis距离的主坐标轴分析(PCoA)代表籼粳稻的根系微生物组在第一主轴分开(P < 0.001，PERMANOVA采用Adonis函数置换检验)。椭圆包括亚种68%的数据。
d. 基于Bray-Cutis距离的PCoA在地块2中结果代表籼粳稻根系微生物组也在第一主轴分开。

c, Unconstrained PCoA (for principal coordinates PCo1 and PCo2) with Bray–Curtis distance showing that the root microbiota of indica separate from those of japonica in field I in the first axis (P < 0.001, permutational multivariate analysis of variance (PERMANOVA) by Adonis). Ellipses cover 68% of the data for each rice subspecies. d, Unconstrained PCoA with Bray–Curtis distance showing that the root microbiota of indica separate from those of japonica in field II in the first axis (P < 0.001, PERMANOVA by Adonis).

结果：

咱们发现不一样水稻亚种根系微生物组成存在差别。基于Bray-Curtis距离的非限制性主坐标轴分析(PCoA)代表籼粳稻在地块1的微生物组成明显造成两大类，且在第一主轴分开(图1c；附图2)，代表水稻亚种分化是微生物组变异的最主要影响因素。同时也观察到了因为地块2的土壤不一样，在地块2存在微生物组的变化(附图3；附表1)。但在两块地块中，籼粳稻显著分开保持一致(图1d；附图3)

We found that the composition of the root bacterial microbiotadiffered in rice subspecies. Unconstrained principal coordinateanalysis (PCoA) of Bray–Curtis distance revealed that theroot microbiota of indica and japonica in field I formed two distinctclusters, which separated along the first coordinate axis(Fig. 1c and Supplementary Fig. 2), indicating that the largest sourceof variation in the rice root microbiota was proximity to the subspeciationbetween indica and japonica. As expected, the root microbiotain field II differed from that in field I due to soil differences(Supplementary Fig. 3 and Supplementary Table 1), but the separationof root microbiota between indica and japonica varieties wasconsistent in the two locations (Fig. 1d and Supplementary Figs. 2and 3).

例2. PCoA时间序列

本文是刘永鑫博士负责分析，于2018年发表在中国科学的一篇文章封面文章，详细描述了水稻田间全生育期根系微生物组的变化规律，发表2年被引31次。详见：手把手带你重现菌群封面文章图表。本文对图1中的B/C子图为例进行说明和点评。

田间水稻微生物组随生育时间变化。以水稻日本晴和IR24为材料,并分别种植于昌平和上庄两地，CP表明北京昌平农场，SZ表明北京海淀上庄。B-C. 主坐标轴分析(PCoA)展现水稻微生物组随时间变化，其中微生物群落结构主要在第1/2轴上随时间变化(B)，而不一样土壤类型主要在第3轴上明显分开(C)

Figure 1 The rice root microbiota in fields shift over rice residence time in the field. B−C, Principle coordinate analysis showing that the rice microbiota shifts with rice residence time in the field and developmental stages in the first axis (B) and separated by geographical locations in the third axis (C).

结果：在全部样品的Bray-Curtis距离的主坐标分析（PCoA）中，土体土壤样品汇集在一块儿，而且水稻根样品在田间和发育阶段的第一坐标轴上沿着水稻的生长时间从土壤开始移动（图1B），代表稻田在田间的停留时间和发育阶段是影响根系微生物组成的主要因素。另外，尽管根微生物组在第三轴上被地理位置清楚地分开了，但水稻的生长时间和根系微生物组的动态变化在两个不一样的地点显示出一致的趋势（图1C）。

In Principle Coordinate Analysis (PCoA) of Bray-Curtis distance from all samples, bulk soil samples clustered together, and rice root samples shifted farfrom the soil across rice residence time in the field and developmental stages in the first coordinate axis (Figure 1B),indicating that rice residence time in the field and developmental stage are main factors influencing the root microbiota composition. Additionally, although the root microbiotawere clearly separated by geographic location in the third axis, the rice residence time and development dependent shift of the root microbiota showed consistent trends in the two separate fields (Figure 1C).

  • 总结

  1. 图1B/C是基于Bray-Curtis距离进行的PCoA分析，采用散点图展现，并按时间顺序填充彩虹色(比单色过渡明显，但对色盲人群不友好，有些杂志不接受)，按不一样生态位和地点设置形状，信息较丰富；通常人类颜色区分明显，把颜色赋予想要表达的第一变量，如本文的时间变量，形态分配给次要因素；

  2. 图1B展现PCo1/2轴，组间最大差别为不一样生态位与时间梯度上的变化，但 不一样地点间是没法很好区分时，咱们还须要继续探索其余主坐标轴。本文在图1C展现PCo1/3轴，可进一步看到1轴的差别与时间变化一致，而3轴能够很好分开不一样地点。

例3.NMDS分析不一样食物昆虫组肠道菌群

本文由荷兰皇家科学院生态研究所的S. Emilia Hannula和中科院遗传发育所朱峰研究员于2019年8月发表于Nature Communications (https://doi.org/10.1038/s41467-019-09284-w)。揭示了食叶昆虫微生物群落来源于土壤而不是取食植物。中文解读详见：[Nature子刊：植食昆虫微生物组来自土壤](https://mp.weixin.qq.com/s/uiXqcGZEt3QX-V49r88J2w)

图a-d 表明了植物群落对土壤、毛虫肠道、根系、植物叶片细菌群落的影响。图e-h表明了植物群落对土壤、毛虫肠道、根系、植物叶片真菌群落的影响。NMDS分析基于Bray-Curtis类似性，二维应力值介于0.11-0.18之间。草地植被相关的群落使用亮绿色表示。非禾本草本植物/阔叶草(forb)植被群落使用青绿色点表示。草地和阔叶草植被混合群落使用深绿色表示。每幅图中小点表明样品，大点表明每组样品的中心点。图中的标识为置换检验结果。a，e表明土壤微生物群落。b，f表明食用离体叶片和植株的毛虫肠道微生物。c，g表明植物根系微生物群落，d，h表明叶微生物群落。

Plant community identity effects on bacterial a–d and fungal (e–h) communities in caterpillars, leaves, roots, and soil. NMDS plots are presented based on Bray–Curtis similarity. The 2D stress value for each panel ranges between 0.11–0.18. Soils originating from grass communities are presented with light green symbols, soils from forb communities with turquoise symbols and soils from mixed grass and forb communities with dark green symbols. In each panel, smaller symbols depict individual samples, centroids are depicted with larger markers. Significance of the plant community treatment effect based on a PERMANOVA is also presented in each panel. a, e represent the composition of microbiomes in soils, b, f microbiomes in caterpillars both on intact plants and on detached leaves. c, g microbiomes in roots and d, h microbiomes in leaves.

结果：

咱们经过两个独立的平行试验，研究了田间植物群落对土壤中微生物群落组成、蒲公英和在这些植物上放养的毛虫的影响。植被群落改变了土壤细菌和真菌群落，可是使人惊讶的是并无改变蒲公英根系和叶片微生物组成 (图3c, d, g, h)。可是咱们却检测到了不一样植物群落对毛虫微生物群落的影响，但这只有在以完整植株为食的毛虫中检测到。

We investigated the legacy effects created by field-grown plant communities, on the composition of microbial communities in soils, dandelions grown in those soils, and caterpillars reared on these plants, in two parallel assays. The composition of the plant community (fast- and slow-growing grasses or forbs) that conditioned the soils that were used, influenced the fungal and bacterial community structure in these soils (Fig. 3a, e). Surprisingly, this did not alter the root- or leaf -associated microbiomes in the dandelion plants that were growing in these soils (Fig. 3c, d, g, h). However, we did detect these soil-derived plant community effects in caterpillar microbiomes, but only when the caterpillars were fed on intact plants (Fig. 3b, f), suggesting that, even though they are plant feeders, the caterpillars had been in direct contact with the soil.

例4.NMDS分析组间的功能基因类群

瑞士EAWAG研究所，西湖大学鞠峰教授于2019年发表于The ISME Journal的成果，发现污水厂抗性组受细菌组成和基因交换驱动且出水中抗性表达活跃(https://doi.org/10.1038/s41396-018-0277-8)。全文解读详见：[ISME：污水中抗性基因在细菌群落和基因交换双重驱动下在出水中活跃表达](https://mp.weixin.qq.com/s/dN0_iQkDSFS42TcFOjKXIA)。

污水厂不一样处理部位抗性基因组的组成与细菌群落组成相关。a-c NMDS分析描述了不一样部位之间基于ARG（a）、BRG（b）、MRG（c）组成的Bary-curtis距离。

Resistome composition correlates with bacterial community composition and phylogeny across wastewater treatment compart-ments. a–c Non-metric multidimensional scaling plots depict Bray-Curtis distances between treatment compartments based on relative abundance of antibiotic (a), biocide (b), and metal (c) resistance genes in the metagenomes.

结果：

细菌抗性组系统发育结构。为了测试在咱们的数据集中是否存在这种状况，咱们使用排序方法来跟踪抗性组（图5）的结构变化。不管分析是基于抗性组、杀菌剂和金属抗性基因的丰度指标（图5a–c），样品始终分为三个主要类别。

Bacterial phylogeny structures soil resistomes. To test if this was the case in our dataset, we used ordination to follow structural variations in the resistomes (Fig. 5) both between and within treatment compartments. The samples consistently clustered into three main groups by treatment compartment with bioreactor samples closely clustered together, whether the analysis was based on abundance metrics of antibiotic, biocide, and metal resistance genes (Fig. 5a–c).

PCoA/NMDS实战

关于更多本项目中示例文件的下载，R包安装的内容，请参考以前的章节：

安装和载入R包
if (!requireNamespace("devtools", quietly=TRUE))
    install.packages("devtools")
library(devtools)
if (!requireNamespace("amplicon", quietly=TRUE))
    install_github("microbiota/amplicon")
suppressWarnings(suppressMessages(library(amplicon)))
主坐标轴分析 PCoA

主坐标轴分析(principal coordinate analysis, PCoA)

在amplicon包中有beta_pcoa函数能够快速绘制PCoA散点图，并按组着色和添加68%的置信椭圆

本次绘制使用函数内置数据演示，查看函数帮助，打问题(?)+函数名，如?beta_pcoa

# 使用内置数据，输入距离矩阵、元数据和分组绘制PCoA
(p=beta_pcoa(beta_bray_curtis, metadata, "Group"))
# 保存位图和矢量图，分别用于预览和出版
ggsave(paste0("p1.PCoA.bray.jpg"), p, width=89, height=56, units="mm")
ggsave(paste0("p1.PCoA.bray.pdf"), p, width=89, height=56, units="mm")

图1. 散点图展现基于Bray-Curtis距离的Beta多样性PCoA。点表明样本，颜色表明分组，并按每组添加68%置信度的椭圆方便组间比较，图中展现主坐标分析的前两轴，解析率见坐标轴括号中。

本次测试数据来自刘永鑫博士负责分析并于2019年发表于Science的文章(即上图展现的内置数据)，讨论了基因型对菌群的影响。详见宏基因组公众号详细解读-Science：拟南芥三萜化合物特异调控根系微生物组

咱们再演示从文件读取距离矩阵和元数据，数据位于Data/Science2019目录，本次须要元数据(metadata.txt)和Beta多样性距离矩阵(alpha/unifrac.txt)两个输入文件(注：距离矩阵这里是由USEARCH -beta_div生成，将在扩增子流程部分详细介绍，也可由vegan包计算生成)。

# USEARCH可选距离矩阵bray_curtis、unifrac、unifrac_binary、jaccard、manhatten、euclidean
# 设置距离矩阵相似，本次使用unifrac
distance_type="unifrac"
# 读取距离矩阵并预测前3行3列，再读取元数据
distance_mat=read.table(paste0("../Data/Science2019/beta/",distance_type,".txt"), header=T, row.names=1, sep="\t", comment.char="")
distance_mat[1:3, 1:3]
metadata=read.table("../Data/Science2019/metadata.txt", header=T, row.names=1, sep="\t", comment.char="", stringsAsFactors=F)
# PCoA散点图，按metadata的Group列着色，添加标签，PCo1/3
(p=beta_pcoa(distance_mat, metadata, groupID="Group", ellipse=T, label=T, PCo=13))
# 保存8：5的半版图
ggsave(paste0("p2.PCoA.unifrac.jpg"), p, width=89, height=56, units="mm")
ggsave(paste0("p2.PCoA.unifrac.pdf"), p, width=89, height=56, units="mm")

图2. 基于Unifrac距离的PCoA。看到PCo 1/2的解析率比前面Bray-Curtis距离结果有提升，代表在Unifrac距离前两主轴通常能够解析更高比例的差别。因为Unifrac考虑进化距离，通常样本/组间差别会进一步缩小。

有时咱们更想知识组间是否存在显著差别，使用?beta_pcoa_stat查看函数帮助，使用距离矩阵指定分组，对所有组别两两差别使用adonis函数进行检测。

# 使用adonis检测组件差别，注意是两两检测，而且将检测结果保存到当前路径下。
beta_pcoa_stat(distance_mat, metadata, "Group", "beta_pcoa_stat.txt")
# 结果文件默认见beta_pcoa_stat.txt
beta_pcoa_stat(dis_mat=distance_mat, metadata=metadata, groupID="Group", pairwise=F, pairwise_list="../Data/Science2019/compare.txt")

输入文件compare.txt，即两组比较列表，制表符分隔。

KO    WT
OE    WT

结果文件beta_pcoa_stat.txt，计算时间和两组比较P值

Sat Jun 27 22:35:03 2020
KO    WT    0.0174982501749825
OE    WT    0.0096990300969903
非度量多维尺度NMDS

咱们将会用到BetaDiv函数，这个函数依赖phyloseq能够计算目前主流的降维排序方法，包括DCA, CCA, RDA, NMDS, MDS, PCoA, PCA, LDA，t-sne，而且结合了群落差别分析，为你们带来相对全面的beta多样性分析。咱们下面以NMDS为例演示函数的用法。?BetaDiv 显示帮助

# 安装Bioconductor的R包phyloseq
if (!requireNamespace("BiocManager", quietly=TRUE))
    install.packages("BiocManager")
suppressWarnings(suppressMessages(library(BiocManager)))
if (!requireNamespace("phyloseq", quietly=TRUE))
    BiocManager::install("phyloseq")
library(phyloseq)

# 输入抽平标准化的特征表、元数据、分组列名、距离类型、降维和统计方法
result=BetaDiv(otu=otutab_rare, map=metadata, group="Group",
                 dist="bray", method="NMDS", Micromet="adonis")
# 返回结果列表：标准图，数据，标签图，成对比较结果，总体结果

#提取排序散点图(结果列表中的1)
(p=result[[1]])
ggsave(paste0("p3.NMDS.bray.jpg"), p, width=89, height=56, units="mm")
ggsave(paste0("p3.NMDS.bray.pdf"), p, width=89, height=56, units="mm")

图3. NMDS分析样本微生物群落结构，按组着色，stress值显示于左上角。

# 提取出图坐标
plotdata=result[[2]]
plotdata[1:3,1:3]

# 提取带标签排序散点图
(p=result[[3]])
ggsave(paste0("p4.NMDS.bray.label.jpg"), p, width=89, height=56, units="mm")
ggsave(paste0("p4.NMDS.bray.label.pdf"), p, width=89, height=56, units="mm")

图4. NMDS分析样本微生物群落结构，添加样本标签。

# 提取两两比较差别检测结果
(pair=result[[4]])

# 提取所有组总体差别检测结果
(Mtest=result[[5]])
了解PhyloSeq对象

输入数据除了支持特征表、元数据+分组；还支持phyloseq对象。

咱们将特征表和元数据转换为PhyloSeq对象(简称ps)

# 指定目标分组列为Group，做为默认分组
metadata$Group=metadata[["Group"]]
# 输入特征表和元数据为PhyloSeq对象
ps=phyloseq(otu_table(as.matrix(otutab),taxa_are_rows=TRUE),
            sample_data(metadata),phy_tree(tree))

固然，除了经常使用的adonis置换检验，可选anosim/MRPP差别显著性检验方法。

result=BetaDiv(ps=ps, dist="bray", method ="NMDS", Micromet ="anosim")
result[[5]]

参考文献

刘尧，Beta多样性和生态类似性，科学网， http://wap.sciencenet.cn/home.php?mod=space&uid=3406804&do=blog&id=1195182

HLBayes，通俗理解PCA降维做用，CSDN，https://blog.csdn.net/HLBoy_happy/article/details/77146012

Xiao-Tao Jiang, Xin Peng, Guan-Hua Deng, Hua-Fang Sheng, Yu Wang, Hong-Wei Zhou & Nora Fung-Yee Tam. (2013). Illumina Sequencing of 16S rRNA Tag Revealed Spatial Variations of Bacterial Communities in a Mangrove Wetland. Microbial Ecology 66, 96-104, doi: https://doi.org/10.1007/s00248-013-0238-8

Jingying Zhang, Yong-Xin Liu, Na Zhang, Bin Hu, Tao Jin, Haoran Xu, Yuan Qin, Pengxu Yan, Xiaoning Zhang, Xiaoxuan Guo, Jing Hui, Shouyun Cao, Xin Wang, Chao Wang, Hui Wang, Baoyuan Qu, Guangyi Fan, Lixing Yuan, Ruben Garrido-Oter, Chengcai Chu & Yang Bai. (2019). NRT1.1B is associated with root microbiota composition and nitrogen use in field-grown rice. Nature Biotechnology 37, 676-684, doi: https://doi.org/10.1038/s41587-019-0104-4

Jingying Zhang, Na Zhang, Yong-Xin Liu, Xiaoning Zhang, Bin Hu, Yuan Qin, Haoran Xu, Hui Wang, Xiaoxuan Guo, Jingmei Qian, Wei Wang, Pengfan Zhang, Tao Jin, Chengcai Chu & Yang Bai. (2018). Root microbiota shift in rice correlates with resident time in the field and developmental stage. Science China Life Sciences 61, 613-621, doi: https://doi.org/10.1007/s11427-018-9284-4

S. Emilia Hannula, Feng Zhu, Robin Heinen & T. Martijn Bezemer. (2019). Foliar-feeding insects acquire microbiomes from the soil rather than the host plant. Nature Communications 10, 1254, doi: https://doi.org/10.1038/s41467-019-09284-w

Feng Ju, Karin Beck, Xiaole Yin, Andreas Maccagnan, Christa S. McArdell, Heinz P. Singer, David R. Johnson, Tong Zhang & Helmut Bürgmann. (2019). Wastewater treatment plant resistomes are shaped by bacterial composition, genetic exchange, and upregulated expression in the effluent microbiomes. The ISME Journal 13, 346-360, doi: https://doi.org/10.1038/s41396-018-0277-8

责编：刘永鑫 中科院遗传发育所

版本更新历史

1.0.0，文涛，初稿

1.0.1，刘永鑫，主题限定为非限制排序，添加实例和整理代码

1.0.2，席娇，文字修改

1.0.3，文涛，删减背景中非必需理论，校对全文

1.0.4，刘永鑫，添加PCoA背景示意图，整理参考文献

猜你喜欢

10000+：菌群分析 宝宝与猫狗 梅毒狂想曲 提DNA发Nature Cell专刊 肠道指挥大脑

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

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

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

必备技能：提问 搜索  Endnote

文献阅读 热心肠 SemanticScholar Geenmedical

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

16S功能预测   PICRUSt  FAPROTAX  Bugbase Tax4Fun

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

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

编程模板: Shell  R Perl

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

写在后面

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

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

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

联系我们 沪ICP备13005482号-10