置换多元方差分析PERMANOVA+主坐标分析PCoA 算法与代码复现

置换多元方差分析 Adonis

（PERMANOVA）别称Adonis

一般而言，多元方差分析的前提假设(多元正态分布等)通常难以满足。

所以需要鲁棒性更强、对数据分布依赖更少的检验方法。Adonis是其中代表性的方法。

通过一个样本间的距离矩阵或相似性矩阵构建ANOVA分析类似的统计量，然后对每组的观测结果进行随机置换来计算显著性P-value。

Fig 1 | 置换多元方差分析计算方式

黄色部分是组间差异，蓝色部分是组内差异，每个值是计算出的距离（一般为bary方法）。

Fig 2 | p值计算方法

SSa是组间方差平方和，SSw是组内方差平方和，通过F统计量从而得到P值。

PCoA 主坐标分析

原理类似主成分分析PCA，即通过特征向量和特征值的分解进行数据降维以便进行可视化，一般与Adonis一起使用。相对于主成分分析，两者的区别为PCA是基于样本的相似系数矩阵(如欧式距离)来寻找主成分，而PCoA是基于距离矩阵(欧式距离以外的其他距离)来寻找主坐标。

Fig 3 主坐标分析

计算距离后，维度会减一，从而2个样本可在一维进行展现，3个样本可在二维进行展现……依次类推。如果是n个样本，则需要进行降维，将高维展现到二维或三维。

在这里，使用PCoA进行降维，因为与置换多元方差分析使用同样的距离计算方式。

代码复现

library(vegan)
library(tidyverse)
#置换多元方差分析
library(pairwiseAdonis)
library(ggpubr)
data(dune) #vegan包中附带的数据，与农业生产有关
data(dune.env)
otu<-t(dune) #测序得到的结果

# Alpha diversity计算方法
alpha_vegan<-function(otu) {
  shannon <- vegan::diversity(otu, index = "shannon")
  simpson <- vegan::diversity(otu, index = "simpson")
  invsimpson <- vegan::diversity(otu, index = "invsimpson")
  richness <-  rowSums(otu > 0)
  alpha <- data.frame(shannon,simpson,invsimpson,richness)
  return (alpha)
}

species<-colnames(dune)
metadata<-dune.env

#A1: 土壤厚度信息 a numeric vector of thickness of soil A1 horizon.
#Moisture: 湿度等级信息，分4个等级，1 < 2 < 4 < 5.
#Management: 分组信息，不同的管理方式 a factor with levels: BF (Biological farming), HF (Hobby farming), NM (Nature Conservation Management), and SF (Standard Farming).
#Use: 一个分组信息 an ordered factor of land-use with levels: Hayfield < Haypastu < Pasture.
#Manure: 一个分组信息，0 < 1 < 2 < 3 < 4.


alpha_species <- alpha_vegan(otu)
alpha_species$sample_id <- rownames(alpha_species)


dune_dist <- vegdist(dune, method="bray", binary=F) #计算距离，使用bray方法



#PcoA计算
dune_pcoa <- cmdscale(dune_dist, k=3, eig=T)
#cmdscale函数用于执行经典多维尺度分析（Classical Multidimensional Scaling，简称MDS），
#也称为主坐标分析。
#dune_dist是输入的距离矩阵，
#k=3表示我们希望将数据降维到三维空间，
#eig=T表示我们希望返回特征值。
dune_pcoa_points <- as.data.frame(dune_pcoa$points)
sum_eig <- sum(dune_pcoa$eig)
eig_percent <- round(dune_pcoa$eig/sum_eig*100,1)

colnames(dune_pcoa_points) <- paste0("PCoA", 1:3)

dune_pcoa_result <- cbind(dune_pcoa_points, dune.env)

# Adonis置换多元方差分析
dune.div <- adonis2(dune ~ Management, data = dune.env, permutations = 999, method="bray")

dune_adonis <- paste0("adonis R2: ",round(dune.div$R2,2), "; P-value: ", dune.div$`Pr(>F)`)




library(ggalt) #用ggalt和ggplot2进行可视化
#推荐ggalt包对ggplot2有一些其他的补充
ggplot(dune_pcoa_result, aes(x=PCoA1, y=PCoA2, color=Management, group = Management)) +
  labs(x=paste("PCoA 1 (", eig_percent[1], "%)", sep=""),
       y=paste("PCoA 2 (", eig_percent[2], "%)", sep=""),
       title=dune_adonis) +
  geom_point(size=5) + 
  geom_encircle(aes(fill=Management), alpha = 0.1, show.legend = F) +
  theme_classic() + coord_fixed(1)