175_微生物组Beta多样性可视化¶
一句话概述¶
Beta多样性衡量不同样本之间微生物群落组成的差异，通过距离矩阵和降维可视化（PCoA/NMDS）将多维群落差异映射到二维空间，是微生物组研究中比较不同处理/环境群落结构差异的核心方法。
核心知识点表格¶
| 知识点 | 说明 |
|---|---|
| Beta多样性定义 | 样本间群落组成差异的度量 |
| 常用距离 | Bray-Curtis（丰度）、UniFrac（系统发育）、Jaccard（存在/缺失） |
| Weighted UniFrac | 考虑丰度的系统发育距离，对优势物种敏感 |
| Unweighted UniFrac | 仅考虑存在与否的系统发育距离，对稀有物种敏感 |
| PCoA | 主坐标分析，保持样本间距离关系的线性降维 |
| NMDS | 非度量多维缩放，基于排序保持距离的非线性降维 |
| PERMANOVA | 置换多元方差分析（adonis2），检验分组是否显著 |
| betadisper | 检验组间方差齐性，PERMANOVA的前提假设检验 |
步骤详解¶
第一步：理解Beta多样性距离¶
白话解释：把每个样本想象成空间中的一个点，样本间微生物组成越相似，点之间的距离就越近。Beta多样性分析就是计算这些"距离"，然后用图形展示出来。
技术细节：不同距离度量关注的生态学信息不同。Bray-Curtis是最常用的生态距离，考虑物种的相对丰度差异；Jaccard只看物种有没有，不看多少；UniFrac还加入了物种间的进化关系。
library(vegan)
library(phyloseq)
# 模拟一个OTU表（行=样本，列=OTU）
set.seed(42)
otu_table <- matrix(rpois(500, lambda = 10), nrow = 20, ncol = 25)
rownames(otu_table) <- paste0("Sample", 1:20)
colnames(otu_table) <- paste0("OTU", 1:25)
# 计算不同距离矩阵
dist_bray <- vegdist(otu_table, method = "bray") # Bray-Curtis
dist_jaccard <- vegdist(otu_table, method = "jaccard") # Jaccard
dist_euclidean <- vegdist(otu_table, method = "euclidean") # Euclidean（不推荐直接用）
# 如果有phyloseq对象，可以计算UniFrac
# dist_unifrac <- UniFrac(physeq, weighted = TRUE)
# dist_uunifrac <- UniFrac(physeq, weighted = FALSE)
第二步：PCoA降维可视化¶
白话解释：距离矩阵是一个n×n的大表，没法直接看。PCoA把这个表"压缩"到二维或三维的图上，让我们能直观地看到哪些样本组成相似（距离近）或不同（距离远）。
技术细节：PCoA（Principal Coordinates Analysis）是对距离矩阵进行特征分解，找到能最大程度保留原始距离信息的低维坐标。每个轴的解释方差百分比反映了该维度捕获的信息量。
library(ape)
library(ggplot2)
# 元数据
metadata <- data.frame(
Sample = paste0("Sample", 1:20),
Group = rep(c("Control", "Treatment"), each = 10),
Timepoint = rep(c("T1", "T2"), 10)
)
# PCoA分析
pcoa_result <- pcoa(dist_bray)
# 提取坐标和解释方差
pcoa_df <- data.frame(
PC1 = pcoa_result$vectors[, 1],
PC2 = pcoa_result$vectors[, 2],
Group = metadata$Group,
Timepoint = metadata$Timepoint
)
var_explained <- round(pcoa_result$values$Relative_eig[1:2] * 100, 1)
# 绘制PCoA图
ggplot(pcoa_df, aes(x = PC1, y = PC2, color = Group, shape = Timepoint)) +
geom_point(size = 3, alpha = 0.8) +
stat_ellipse(aes(group = Group), level = 0.95, linetype = "dashed") +
labs(
x = paste0("PCoA1 (", var_explained[1], "%)"),
y = paste0("PCoA2 (", var_explained[2], "%)"),
title = "Beta Diversity (Bray-Curtis)"
) +
theme_minimal() +
theme(legend.position = "right")
第三步：NMDS降维可视化¶
白话解释：NMDS是另一种降维方法，它不像PCoA那样追求精确保留距离值，而是保持样本间距离的"排名顺序"。有时候NMDS能比PCoA更好地展示复杂的群落模式。
技术细节：NMDS通过迭代优化将样本映射到低维空间，使得低维空间中的距离排序尽可能接近原始距离排序。stress值衡量拟合质量：<0.05表示极好，<0.1表示好，<0.2表示可接受，>0.2表示较差。
# NMDS分析
set.seed(123)
nmds_result <- metaMDS(otu_table, distance = "bray", k = 2, trymax = 100)
# 查看stress值
cat("Stress:", nmds_result$stress, "\n")
# 提取坐标
nmds_df <- data.frame(
NMDS1 = nmds_result$points[, 1],
NMDS2 = nmds_result$points[, 2],
Group = metadata$Group
)
# 绘制NMDS图
ggplot(nmds_df, aes(x = NMDS1, y = NMDS2, color = Group)) +
geom_point(size = 3) +
stat_ellipse(level = 0.95) +
annotate("text", x = Inf, y = Inf,
label = paste("Stress =", round(nmds_result$stress, 3)),
hjust = 1.1, vjust = 1.5, size = 4) +
theme_bw() +
labs(title = "NMDS (Bray-Curtis)")
第四步：统计检验（PERMANOVA）¶
白话解释：光看图还不够，我们需要统计检验来判断分组之间的差异是否"真实"而非随机。PERMANOVA就是专门用于此的置换检验方法。
技术细节：PERMANOVA（adonis2）将距离矩阵的总变异分解为组间和组内部分，通过置换检验计算p值。R²值表示分组因素解释的变异比例。注意PERMANOVA假设组间方差齐性，需用betadisper先检验。
# PERMANOVA检验
set.seed(999)
adonis_result <- adonis2(dist_bray ~ Group * Timepoint, data = metadata, permutations = 999)
print(adonis_result)
# 方差齐性检验
bd <- betadisper(dist_bray, metadata$Group)
permutest(bd, pairwise = TRUE)
# p > 0.05 表示方差齐性，PERMANOVA结果可靠
# 多组成对比较（如果有>2组）
library(pairwiseAdonis)
pairwise.adonis(dist_bray, metadata$Group, p.adjust.m = "BH")
# ANOSIM检验（替代方法）
anosim_result <- anosim(dist_bray, metadata$Group, permutations = 999)
summary(anosim_result)
第五步：使用phyloseq进行整合分析¶
library(phyloseq)
# 创建phyloseq对象
OTU <- otu_table(otu_table, taxa_are_rows = FALSE)
META <- sample_data(metadata)
rownames(META) <- metadata$Sample
physeq <- phyloseq(OTU, META)
# phyloseq内置的排序方法
ord_pcoa <- ordinate(physeq, method = "PCoA", distance = "bray")
ord_nmds <- ordinate(physeq, method = "NMDS", distance = "bray")
# 使用plot_ordination
p <- plot_ordination(physeq, ord_pcoa, color = "Group", shape = "Timepoint") +
geom_point(size = 3) +
stat_ellipse(aes(group = Group)) +
theme_minimal() +
ggtitle("PCoA - Bray-Curtis")
print(p)
第六步：Python实现（scikit-bio + matplotlib）¶
import numpy as np
import pandas as pd
from skbio import DistanceMatrix
from skbio.stats.ordination import pcoa
from skbio.stats.distance import permanova
from scipy.spatial.distance import braycurtis, pdist, squareform
import matplotlib.pyplot as plt
# 计算Bray-Curtis距离矩阵
otu = np.random.poisson(10, (20, 25))
sample_ids = [f"Sample{i}" for i in range(1, 21)]
bc_distances = squareform(pdist(otu, metric='braycurtis'))
dm = DistanceMatrix(bc_distances, ids=sample_ids)
# PCoA
pcoa_result = pcoa(dm)
pcoa_df = pcoa_result.samples[['PC1', 'PC2']]
pcoa_df['Group'] = ['Control']*10 + ['Treatment']*10
# 可视化
fig, ax = plt.subplots(figsize=(8, 6))
for group, color in [('Control', 'blue'), ('Treatment', 'red')]:
mask = pcoa_df['Group'] == group
ax.scatter(pcoa_df.loc[mask, 'PC1'], pcoa_df.loc[mask, 'PC2'],
c=color, label=group, s=60, alpha=0.7)
ax.set_xlabel(f"PC1 ({pcoa_result.proportion_explained[0]*100:.1f}%)")
ax.set_ylabel(f"PC2 ({pcoa_result.proportion_explained[1]*100:.1f}%)")
ax.legend()
ax.set_title("PCoA (Bray-Curtis)")
plt.tight_layout()
plt.savefig("beta_diversity_pcoa.png", dpi=300)
# PERMANOVA
grouping = pd.Series(['Control']*10 + ['Treatment']*10, index=sample_ids)
result = permanova(dm, grouping, permutations=999)
print(result)
实战命令速查¶
# QIIME2中的Beta多样性
# qiime diversity beta --i-table table.qza --p-metric braycurtis --o-distance-matrix bray.qza
# qiime diversity pcoa --i-distance-matrix bray.qza --o-pcoa pcoa.qza
# qiime emperor plot --i-pcoa pcoa.qza --m-metadata-file meta.tsv --o-visualization emperor.qzv
# R一行完成PCoA+可视化
library(microeco); library(magrittr)
# dataset$cal_betadiv(unifrac = FALSE) %>% dataset$plot_betadiv()
# 快速计算多种距离
lapply(c("bray","jaccard","horn","morisita"), function(m) vegdist(otu_table, method = m))
面试常问点¶
Q1: Bray-Curtis和UniFrac距离的区别是什么？ A: Bray-Curtis基于物种丰度计算差异，不考虑物种间的进化关系。UniFrac结合系统发育树信息，亲缘关系近的物种差异被降低权重。如果两个样本有不同物种但这些物种亲缘关系很近，UniFrac距离会小，Bray-Curtis距离可能较大。
Q2: PCoA和NMDS如何选择？ A: PCoA是线性方法，计算快速确定性，适合大多数情况。NMDS是非线性方法，对复杂的非线性关系表现更好，但需要迭代且结果可能不唯一。实际中建议两者都做，结果一致则更可靠。
Q3: PERMANOVA的R²值如何解读？ A: R²表示分组因素解释的总变异比例。例如R²=0.15表示分组解释了15%的群落变异。微生物组研究中R²通常在0.05-0.30之间。R²小但p值显著说明存在真实但较小的差异。
Q4: PERMANOVA显著但PCoA图上看不出分离怎么办？ A: 这通常是因为组间差异虽然统计显著但效应量小，或者差异主要在PCoA未展示的高阶轴上。可以查看PCoA的第3、4轴，或使用dbRDA约束排序方法直接展示分组相关的变异。
Q5: 为什么不直接用PCA分析微生物组数据？ A: PCA假设变量间呈线性关系，使用欧氏距离，不适合物种丰度数据（含大量零值、高度偏态）。PCoA可以使用生态距离，NMDS可以处理非线性关系，都比PCA更适合微生物组数据。
Q6: betadisper检验不通过怎么办？ A: betadisper不通过意味着组间方差不齐，PERMANOVA的p值可能被膨胀。可以报告betadisper结果作为补充，或使用其他方法如ANOSIM（对方差齐性假设更宽松）。
易错点¶
- 忘记数据标准化：直接用原始counts计算Bray-Curtis会被测序深度主导，应先做相对丰度或rarefaction
- PCoA负特征值：某些距离（如UniFrac）可能产生负特征值，需用correction参数校正
- NMDS的stress值过高：stress>0.2时应增加维度（k=3），或检查数据质量
- PERMANOVA嵌套设计遗漏：纵向数据需在PERMANOVA中指定strata参数控制重复测量
- 椭圆图的误解：PCoA图上的椭圆是95%置信椭圆，不代表统计显著性
补充知识¶
距离度量选择指南¶
| 距离 | 适用场景 | 特点 |
|---|---|---|
| Bray-Curtis | 通用，最常用 | 考虑丰度，0-1范围 |
| Jaccard | 关注物种组成 | 仅考虑有/无 |
| Weighted UniFrac | 关注优势物种的系统发育差异 | 需要进化树 |
| Unweighted UniFrac | 关注稀有物种的系统发育差异 | 需要进化树 |
| Aitchison | 组成型数据 | 基于CLR变换+欧氏距离 |
约束排序方法¶
当想直接展示环境因子对群落的影响时，可以使用约束排序方法（如db-RDA）：