Beta多样性分析：一句代码搞定CCA、DCA、DPCoA、MDS、NMDS、PCoA、RDA

 数据分析服务请访问以下链接：

文章发表技术服务，数据分析服务

数据和代码获取：请查看主页个人信息

大家好，今天我将介绍如何使用R语言phyloseq包进行常见的beta多样性分析，并使用ggplot2包进行绘图展示。

相信大家在分析的过程中会遇到一个头疼的点：距离计算方法和beta多样性指数计算方法那么多，我该怎么计算？或者说用哪种指数进行分析后，样本区分度可以达到我的需求？

首先来看下chatGPT对于这七种多样性指数的解释：

1. CCA (Canonical Correspondence Analysis)：CCA是一种用于分析响应变量（通常是物种丰度数据）与解释变量（例如环境变量）之间关系的多变量技术。它旨在识别最大化响应变量与解释变量之间相关性的线性组合。

2. DCA (Detrended Correspondence Analysis)：DCA是一种用于排序和可视化生态数据的方法。通常应用于物种丰度数据，以探索样本或群落之间的差异。DCA试图展现数据中的梯度信息，并且可以显示数据中的非线性关系。

3. DPCoA (Double Principal Coordinate Analysis)：DPCoA是一种多元数据分析方法，类似于PCoA（下面解释）。它用于比较和可视化多个样本或群落之间的差异，并且可以基于任意距离矩阵进行操作。

4. MDS (Multidimensional Scaling)：MDS是一种用于可视化和分析数据的技术。它将样本或群落之间的距离或相似性矩阵转换为低维空间中的坐标，以便在图形上表示它们的相对位置关系。

5. NMDS (Non-metric Multidimensional Scaling)：NMDS是一种在MDS基础上发展而来的技术。与传统的MDS不同，NMDS可以处理非度量（非线性）距离矩阵，使其更适用于生态学中的多样性分析。

6. PCoA (Principal Coordinate Analysis)：PCoA是一种用于分析样本或群落之间差异的技术。它将距离或相似性矩阵转换为低维空间中的坐标，以可视化样本之间的相对位置关系。

7. RDA (Redundancy Analysis)：RDA是一种多元回归技术，常用于分析生态数据中的环境变量对物种组成的解释。它结合了多元线性回归和主成分分析的概念。RDA通过最大化解释物种组成差异的线性组合，将环境变量与生物群落之间的关系可视化。

今天我将介绍使用phyloseq包一次性实现7种Beta多样性分析方法，接下来我们来进行分析和可视化展示：

Step1：数据载入

rm(list=ls())
pacman::p_load(tidyverse,phyloseq,plyr,ggsci,ape)
otu <- read.csv('otu_table.csv', row.names = 1)
tax <- read.csv('tax_table.csv', row.names = 1)
sam <- read.csv('sample_data.csv', row.names = 1)
phy <- read.tree('phy_tree.nwk')

Step2：构建phyloseq对象

GP <- 
  phyloseq(otu_table(as.matrix(otu), taxa_are_rows=TRUE),
           sample_data(sam),
           tax_table(as.matrix(tax)),
           phy_tree(phy))
GP

phyloseq包在微生物组学数据分析的地位我们不用介绍了，相信大家都有一定了解。这种功能强大的R包第一步都需要构建一个数据对象，这里的操作与microeco包的对象构建比较类似：R语言实战：microeco包轻松实现微生物组LEfSe分析。随机森林分析：R包microeco实现微生物组随机森林分析及重要变量的选择 

Step3：基于Bray-Burtis距离的beta多样性

dist = "bray"
ord_meths = c("DCA", "CCA", "RDA", "DPCoA", "NMDS", "MDS", "PCoA")
plist = llply(as.list(ord_meths), function(i, physeq, dist){
  ordi = ordinate(physeq, method=i, distance=dist)
  plot_ordination(physeq, ordi, "samples", color="SampleType")
}, GP, dist)
names(plist) <- ord_meths

pdataframe = ldply(plist, function(x){
  df = x$data[, 1:2]
  colnames(df) = c("Axis_1", "Axis_2")
  return(cbind(df, x$data))
})
names(pdataframe)[1] = "method"

ggplot(pdataframe, aes(Axis_1, Axis_2, color=SampleType, shape=human, fill=SampleType)) +
  geom_point(size = 1) + 
  facet_wrap(~method, scales="free", nrow = 2) + 
  scale_color_npg() +
  scale_fill_npg() +
  theme_bw(base_size = 10)
ggsave('beta on brey.png', width = 8, height = 3.5)

 Step4：基于Unweighted-UniFrac距离的beta多样性

dist = "unifrac"
ord_meths = c("DCA", "CCA", "RDA", "DPCoA", "NMDS", "MDS", "PCoA")
plist = llply(as.list(ord_meths), function(i, physeq, dist){
  ordi = ordinate(physeq, method=i, distance=dist)
  plot_ordination(physeq, ordi, "samples", color="SampleType")
}, GP, dist)
names(plist) <- ord_meths

pdataframe = ldply(plist, function(x){
  df = x$data[, 1:2]
  colnames(df) = c("Axis_1", "Axis_2")
  return(cbind(df, x$data))
})
names(pdataframe)[1] = "method"

ggplot(pdataframe, aes(Axis_1, Axis_2, color=SampleType, shape=human, fill=SampleType)) +
  geom_point(size = 1) + 
  facet_wrap(~method, scales="free", nrow = 2) + 
  scale_color_npg() +
  scale_fill_npg() +
  theme_bw(base_size = 10) + 
  ggtitle("Unweighted-UniFrac distance")
ggsave('beta on unweighted-UniFrac distance.png', width = 8, height = 3.5)

Step5：挑选你心仪的效果图

theme_set(theme_bw(base_size = 7))
cowplot::plot_grid(plist[[4]], plist[[6]], labels = c("a","b"), nrow=1)
ggsave('pic.png', width = 7, height = 2.5)

 关键词“beta多样性”获取数据和代码