写在前面【科普】

什么是PCoA？

主坐标分析（Principal Coordinates Analysis，PCoA），也称为经典多维尺度分析（Classical Multidimensional Scaling，CMDS），是一种多变量分析技术，用于将复杂的多维数据映射到较低的维度（通常是二维或三维），以便进行可视化和解释。

工作原理：

• 输入数据：PCoA通常从一个距离矩阵开始，该矩阵表示数据点（如样本、物种、基因等）之间的成对距离或不相似性。
• 特征值分解：PCoA通过特征值分解将距离矩阵转换为坐标矩阵，生成一组主坐标，这些坐标描述了数据点在降低维度后的相对位置。
• 特征值和特征向量：特征值表示每个主坐标解释的总方差（即信息量）的比例。第一个主坐标解释了最大方差，第二个解释了次大方差，依此类推。

可视化：

• PCoA图：PCoA的结果通常会以散点图的形式展示。每个点代表一个数据点（例如一个样本），点之间的距离表示它们在原始多维空间中的相似性或差异性。
• 解释主坐标：通过PCoA图，可以观察到数据点如何聚类（或分离），从而推测不同组之间的差异或相似性。

PCoA的应用场景

PCoA在以下领域中具有广泛应用：

1. 生态学与环境科学：

   • 物种多样性分析：在群落生态学中，PCoA常用于分析不同环境条件下物种群落的相似性。例如，研究不同样地的植物或动物群落结构。
   • 群落构成比较：比较不同样本（如水体或土壤样本）的微生物群落结构，以评估环境影响或污染源的作用。

2. 基因组学与分子生物学：

   • 基因表达数据：在基因表达分析中，PCoA可以帮助研究人员将高维的基因表达数据简化，并可视化不同样本或条件下的基因表达模式。
   • 微生物组分析：分析微生物群落的群落组成数据（如16S rRNA基因测序数据），评估不同群体或环境之间的差异。

3. 人类学与进化生物学：

   • 遗传差异分析：在遗传学研究中，PCoA可以用来分析不同种群或个体的遗传变异，并可视化其分布情况，揭示遗传背景或进化关系。

4. 社会科学与市场研究：

   • 问卷调查数据：PCoA可以用来分析消费者的态度或偏好，帮助研究人员在多维数据中找到主要的差异来源，并进行细分市场分析。

什么时候使用PCoA？

• 数据维度高：当你处理高维数据时（如基因表达数据、微生物组数据），PCoA可以帮助你在降低维度的同时保留尽可能多的信息，并且提供易于解释的图形表示。
• 距离矩阵分析：如果你已经有一个表示样本之间相似性或差异性的距离矩阵，PCoA是将其转化为可视化的理想工具。
• 探索性分析：当你想要发现数据中的潜在模式或结构时，PCoA可以帮助你识别数据中显著的群组或趋势。
• 群体之间的比较：当需要比较不同群体（如实验组与对照组、不同地区的物种群落）之间的异同，并且这些群体在多个变量上存在差异时，PCoA是一种非常有用的工具。

总之，PCoA通过将复杂的多维数据简化为可视化的低维表示，使研究人员能够更直观地理解数据的结构和样本之间的关系。

第一次做

加载相关软件包

library("ggplot2")
library("vegan")

读入实验设计

design = read.table("design.txt", header=T, row.names= 1, sep="\t") 

design

读入bray_curtis的距离矩阵

bray_curtis = read.table("bray_curtis_dm.txt", sep="\t", header=T, row.names = 1)
colnames(bray_curtis) = rownames(bray_curtis)

bray_curtis

将距离矩阵进行主坐标轴分析

pcoa = cmdscale(bray_curtis, k=3, eig=T) # k is dimension, 3 is recommended; eig is eigenvalues
points = as.data.frame(pcoa$points) # get coordinate string, format to dataframme
colnames(points) = c("x", "y", "z") 
eig = pcoa$eig
points = cbind(points, design[match(rownames(points), rownames(design)), ])

plot(points$x,points$y,
     main = "PCOA of ...",
     xlab = paste("PCoA 1 (", format(100 * eig[1] / sum(eig), digits=4), "%)", sep=""),
     ylab = paste("PCoA 2 (", format(100 * eig[2] / sum(eig), digits=4), "%)", sep=""),
     col=c("black"),
     pch=(16)
     )
dev.off()

plot

绘制主标准轴的第1，2轴

p = ggplot(points, aes(x=x, y=y, color=Group)) +
  geom_point(alpha=.7, size=2) + 
  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=""),
       title="bray_curtis PCoA")
p
ggsave("beta_pcoa_bray_curtis.pdf", p, width = 5, height = 3)
ggsave("beta_pcoa_bray_curtis.png", p, width = 5, height = 3)

结果

第二次做

### 清除所有变量
rm(list = ls())

1. 导入所需的库

library(vegan)      # 社区生态学分析：多元分析、物种多样性分析等
library(tidyverse)  # 数据科学工具集：ggplot2、dplyr、tidyr、readr等
library(ggalt)      # 增强ggplot2功能，提供额外的坐标系统和几何对象
library(car)        # 线性回归分析的有用函数
library(ggforce)    # ggplot2扩展，提供新几何对象和坐标系统
library(ggpubr)     # 创建出版质量图形的便捷函数
library(patchwork)  # 组合多个ggplot对象

2. 定义pairwise adonis函数

pairwise.adonis1 <- function(x, factors, p.adjust.m) {
    x <- as.matrix(x)
    co <- as.matrix(combn(unique(factors), 2))
    pairs <- F.Model <- R2 <- p.value <- c()

    # 对每一对factors进行adonis分析
    for (elem in 1:ncol(co)) {
        subset <- factors %in% c(as.character(co[1, elem]), as.character(co[2, elem]))
        ad <- adonis(x[subset, subset] ~ factors[subset], permutations = 999)

        pairs <- c(pairs, paste(co[1, elem], 'vs', co[2, elem]))
        F.Model <- c(F.Model, ad$aov.tab[1, 4])
        R2 <- c(R2, ad$aov.tab[1, 5])
        p.value <- c(p.value, ad$aov.tab[1, 6])
    }

    p.adjusted <- p.adjust(p.value, method = p.adjust.m)
    pairw.res <- data.frame(pairs, F.Model, R2, p.value, p.adjusted)

    return(pairw.res)
}

3. 读取和处理数据

# setwd("~/Desktop")
otu <- read.table("./28_gene_fpkm_filter0.1.xls", row.names = 1, sep = "\t", header = TRUE) %>% as.data.frame()
map <- read.table("./meta.xls", sep = "\t", header = TRUE)
colnames(map)[1] <- "ID"
row.names(map) <- map$ID

# 仅保留map中与otu列名相同的行
idx <- rownames(map) %in% colnames(otu)
map1 <- map[idx,]
otu <- otu[, rownames(map1)]

4. 进行adonis分析并计算统计值

bray_curtis <- vegan::vegdist(t(otu), method = "bray", na.rm = TRUE)
ado <- adonis(bray_curtis ~ map1$Group, permutations = 999, method = "bray")
R2_value <- round(ado$aov.tab[1, "R2"], 3)
p_v_value <- ado$aov.tab[1, "Pr(>F)"]
title <- paste("Adonis: R^2 =", R2_value, "P_value =", p_v_value)

5. 绘制PCoA图

pcoa <- cmdscale(bray_curtis, k = 2, eig = TRUE)
points <- as.data.frame(pcoa$points) %>% rename(x = "V1", y = "V2")
points <- cbind(points, map1[match(rownames(points), map1$ID),])
eig <- pcoa$eig

# 定义颜色和形状
colors <- c("o" = "#00BFFF", "y" = "#FF4500")
shapes <- c("BRO_o_1" = 24, "BRO_o_2" = 22, "BRO_o_3" = 21, "BRO_o_4" = 23,
            "BRO_y_1" = 24, "BRO_y_2" = 22, "BRO_y_3" = 21, "BRO_y_4" = 23)
levels_order <- c("o", "y")
points$Group <- factor(points$Group, levels = levels_order)

# 绘制PCoA图
p1 <- ggplot(points, aes(x = x, y = y, fill = Group, shape = Group)) + 
    geom_point(alpha = .7, size = 3) +  # 增加点的大小以便更清晰地查看
    scale_shape_manual(values = c(24, 22)) +  # 对不同组使用不同的形状
    scale_fill_manual(values = colors) + 
    xlim(-.5, .5) +  # 设置x轴范围
    ylim(-.25, .25) +  # 设置y轴范围
    
    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 = ""),
         title = title) + 
    geom_mark_ellipse(aes(fill = Group, label = Group), 
                      alpha = 0.1, color = "grey", linetype = 3) +  # 使用Group而不是ID进行标记
    theme_bw() + 
    theme(panel.grid.major = element_blank(), 
          panel.grid.minor = element_blank(), 
          axis.text = element_text(color = "black", size = 9))

# 显示图
p1

# 6. 输出pairwise adonis结果并保存为文本文件
pair_bray_adonis <- pairwise.adonis1(bray_curtis, map1$Group, p.adjust.m = "bonferroni")
write.table(as.data.frame(pair_bray_adonis), "table.txt", sep = "\t", quote = FALSE, row.names = FALSE)

# 创建并显示表格
p2 <- ggtexttable(pair_bray_adonis, rows = NULL)
p2

### 保存PCoA图为PNG格式
ggsave(filename = "PCoA_plot.png", plot = p1, width = 12, height = 10, units = "in", dpi = 600)

p1