7. vegan包进行PCoA分析和可视化

简介

PCoA分析，即主坐标分析（principal co-ordinates analysis），是一种非约束性的数据降维分析方法，可用来研究样本的相似性或差异性，与PCA分析类似；但相比于PCA，PCoA以样本距离为整体考虑，更符合生态学数据特征，应用也更为广泛。

PCoA分析，首先对一系列的特征值和特征向量进行排序，然后选择排在前几位的最主要特征值，经过投影后并将其投影在坐标系里，结果相当于是距离矩阵的一个旋转，在低维度空间以最大限度地保留原始样本的距离关系；相似的样本在图形中的距离更为接近，相异的样本距离更远。

示例

PCoA plot showing the difference between sample groups, stratified to only samples originating from participants not receiving any topical treatment. Pvalue corresponds to Adonis PERMANOVA test. Ellipses delineate the 75% prediction areas of samples from each group.

脚本

数据样式：

1. OTU丰度数据就是一般OTU表或注释后的OTU丰度表，每一行为一个OTU，每一列为一个样品。
2. 分组数据为跟样品一一对应的分组数据。

vegan包的分析结果解释：eig记录了PCoA排序结果中，主要排序轴的特征值（再除以特征值总和就是各轴的解释量）；points记录了各样本在各排序轴中的坐标值。

library(readxl)
library(ggplot2)
library(learn)
library(patchwork)
library(tidyverse)
rm(list = ls())
file <-  "C:\\Users\\...total_data\\"

genes_abundance <- read.table(file = paste0(file, "otu_table_g_relative.xls"),
                              header = TRUE, stringsAsFactors = FALSE)
genes_abundance <- genes_abundance[-ncol(genes_abundance)]
str(genes_abundance)
which(duplicated(genes_abundance$Taxonomy) == TRUE)

groups <- read_xls(path = paste0(file, "the_information_of_sample_site.xls"),
                   sheet = 3)

row.names(genes_abundance) <- genes_abundance$Taxonomy
otu <- genes_abundance[-1]
otu <- data.frame(t(otu))
head(otu)
#排序（基于 OTU 丰度表）
library(vegan)
distance <- vegdist(otu, method = 'bray')
pcoa <- cmdscale(distance, k = (nrow(otu) - 1), eig = TRUE)


# 可视化数据提取 ------------------------------------------------


# 提取样本点坐标（points记录了各样本在各排序轴中的坐标值）
# 前两轴
plot_data <- data.frame({pcoa$point})[1:2]

# 提取列名，便于后面操作。
plot_data$Sample_name <- rownames(plot_data)
names(plot_data)[1:2] <- c('PCoA1', 'PCoA2')

# eig记录了PCoA排序结果中，主要排序轴的特征值（再除以特征值总和就是各轴的解释量）
eig = pcoa$eig

#为样本点坐标添加分组信息
plot_data <- merge(plot_data, groups, by = 'Sample_name', all.x = TRUE)

# 绘制主标准轴的第1，2轴
ggplot(data = plot_data, aes(x=PCoA1, y=PCoA2, color=Group3)) +
  geom_point(alpha=.7, size=2) +
stat_chull(fill =NA) +
  labs(x=paste("PCoA 1 (", format(100 * eig[1] / sum(eig), digits=4), "%)", sep=""),
       y=paste("PCoA 2 (", format(100 * eig[2] / sum(eig), digits=4), "%)", sep=""))

Reference

Ring HC, Thorsen J, Saunte DM, Lilje B, Bay L, Riis PT, Larsen N, Andersen LO, Nielsen HV, Miller IM, Bjarnsholt T, Fuursted K, Jemec GB. The Follicular Skin Microbiome in Patients With Hidradenitis Suppurativa and Healthy Controls. JAMA Dermatol. 2017 Sep 1;153(9):897-905. doi: 10.1001/jamadermatol.2017.0904.