R语言做主坐标举例分析

发布时间：2021-11-22 15:54:42 作者：iii
来源：亿速云阅读：310

# R语言做主坐标分析（PCoA）详解与实战案例

主坐标分析（Principal Coordinates Analysis, PCoA）是一种经典的多变量数据降维技术，广泛应用于生态学、微生物组学等领域。本文将详细介绍PCoA原理，并通过R语言实现完整分析流程。

## 一、PCoA基本概念

### 1.1 什么是PCoA？
PCoA（Principal Coordinates Analysis）又称经典多维尺度分析（Classical Multidimensional Scaling），是一种基于距离矩阵的降维方法。与PCA（主成分分析）不同，PCoA可以处理非欧几里得距离，适用于更复杂的数据类型。

### 1.2 核心数学原理
PCoA通过以下步骤实现降维：

1. 计算样本间距离矩阵D
2. 对距离矩阵进行双重中心化：A = -1/2 * D²
3. 计算矩阵A的特征值和特征向量
4. 选择前k个最大特征值对应的特征向量构成新坐标

数学表达式：
\[ B = HAH \]
其中H是中心化矩阵：\( H = I - \frac{1}{n}11^T \)

## 二、R语言实现PCoA

### 2.1 准备环境与数据
首先加载必要的R包：

```R
# 安装必要包（首次需要）
install.packages(c("vegan", "ggplot2", "ape"))

# 加载包
library(vegan)
library(ggplot2)
library(ape)

使用经典微生物组数据集dune作为示例：

data(dune)
head(dune[,1:5])  # 查看前5列

2.2 计算距离矩阵

PCoA分析的关键是选择合适的距离测度。这里使用Bray-Curtis距离：

dist_matrix <- vegdist(dune, method = "bray")

其他常用距离方法： - "euclidean"：欧式距离 - "jaccard"：Jaccard距离 - "manhattan"：曼哈顿距离

2.3 执行PCoA分析

使用cmdscale()函数进行PCoA：

pcoa_result <- cmdscale(dist_matrix, k = 3, eig = TRUE)

参数说明： - k：保留的维度数 - eig：是否返回特征值

2.4 结果解释

查看分析结果：

# 查看前4个样本的坐标
head(pcoa_result$points, 4)

# 特征值（解释方差）
pcoa_result$eig[1:5]

# 计算各轴解释度
explained_var <- pcoa_result$eig / sum(pcoa_result$eig[pcoa_result$eig > 0])
round(explained_var[1:3]*100, 2)

三、结果可视化

3.1 基础散点图

plot(pcoa_result$points[,1], pcoa_result$points[,2],
     xlab = paste0("PCoA1 (", round(explained_var[1]*100,1), "%)"),
     ylab = paste0("PCoA2 (", round(explained_var[2]*100,1), "%)"),
     main = "PCoA Plot")

3.2 使用ggplot2美化

创建更专业的可视化：

# 准备绘图数据
plot_data <- data.frame(
  PCoA1 = pcoa_result$points[,1],
  PCoA2 = pcoa_result$points[,2],
  Sample = rownames(dune)
)

ggplot(plot_data, aes(x = PCoA1, y = PCoA2)) +
  geom_point(size = 3, alpha = 0.7) +
  labs(x = paste0("PCoA1 (", round(explained_var[1]*100,1), "%)"),
       y = paste0("PCoA2 (", round(explained_var[2]*100,1), "%)")) +
  theme_minimal() +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_vline(xintercept = 0, linetype = "dashed")

3.3 添加分组信息

假设我们有环境因子数据：

data(dune.env)
head(dune.env)

# 合并分组信息
plot_data$Management <- dune.env$Management

# 分组可视化
ggplot(plot_data, aes(x = PCoA1, y = PCoA2, color = Management)) +
  geom_point(size = 3) +
  stat_ellipse(level = 0.95) +  # 添加95%置信椭圆
  scale_color_brewer(palette = "Set1") +
  theme_bw()

四、进阶分析技巧

4.1 显著性检验

使用PERMANOVA检验组间差异：

adonis2(dist_matrix ~ Management, data = dune.env)

4.2 环境因子拟合

将环境变量投影到PCoA图中：

env_fit <- envfit(pcoa_result$points[,1:2], dune.env)
plot(env_fit, col = "blue", cex = 0.8)

4.3 三维可视化

当需要展示第三维度时：

library(plotly)
plot_ly(x = pcoa_result$points[,1],
        y = pcoa_result$points[,2],
        z = pcoa_result$points[,3],
        color = dune.env$Management,
        type = "scatter3d",
        mode = "markers")

五、PCoA与PCA的比较

特征	PCoA	PCA
输入数据	距离矩阵	原始数据矩阵
距离限制	支持非欧式距离	仅欧式距离
计算方式	特征值分解	奇异值分解
适用场景	生态学、进化分析	连续变量分析

六、实际应用案例

6.1 微生物β多样性分析

在16S rRNA测序数据分析中，PCoA常用于展示样本间β多样性差异：

# 示例代码（需替换为实际OTU表）
otu_table <- read.csv("otu_table.csv", row.names = 1)
dist_bray <- vegdist(otu_table, method = "bray")
pcoa_microbiome <- cmdscale(dist_bray, eig = TRUE)

6.2 遗传距离可视化

在群体遗传学中展示个体间遗传距离：

genetic_dist <- as.dist(read.table("genetic_distance.txt"))
pcoa_genetic <- cmdscale(genetic_dist, k = 3)

七、常见问题解答

Q1：如何选择PCoA维度数？ A：通常选择肘部法则或累计解释率>70%的维度数

Q2：PCoA结果出现负特征值怎么办？ A：这是非欧式距离的常见现象，可忽略负值或使用修正方法

Q3：样本点重叠严重如何解决？ A：尝试以下方法： 1. 使用不同距离测度 2. 增加jitter扰动 3. 使用三维展示

八、总结

PCoA作为强大的多变量分析工具，在R语言中可通过cmdscale()或vegan包轻松实现。关键步骤包括： 1. 选择合适距离测度 2. 检查特征值解释率 3. 结合统计检验 4. 制作专业可视化

通过本文的完整案例，读者可以快速掌握PCoA分析的R语言实现方法，应用于自己的研究领域。

注意：实际分析时应根据数据特点调整参数，建议配合其他多元统计方法共同验证结果。 “`