675 微生物组多变量统计分析¶

一句话概述：PERMANOVA检验组间差异是否显著，dbRDA/CCA/RDA揭示环境因素如何驱动群落变化——它们是微生物组研究中β多样性分析的核心统计方法。

核心知识点速查表¶

知识点	关键内容
PERMANOVA	基于置换的多变量方差分析，检验组间差异
dbRDA	基于距离的冗余分析，约束排序
CCA	典型对应分析，适合单峰响应数据
RDA	冗余分析，适合线性响应数据
ANOSIM	组间vs组内距离比较（被PERMANOVA取代）
envfit	将环境变量拟合到排序空间

一、为什么需要多变量分析？（白话解释）¶

打个比方：微生物组有几百上千个物种，就像一张有几百列的数据表。你不可能一列一列去比较（多重检验会崩溃）。多变量分析把这几百列"压缩"成少数几个维度，然后在这个压缩空间里检验组间是否有差异。

方法分类： - 非约束排序：PCoA/NMDS——"让数据自己说话"，不考虑分组/环境信息 - 约束排序：CCA/RDA/dbRDA——"用环境变量解释群落变化" - 假设检验：PERMANOVA/ANOSIM——"组间差异是否统计上显著"

二、PERMANOVA（核心检验方法）¶

# PERMANOVA分析
library(vegan)        # 生态学统计核心包
library(phyloseq)     # 微生物组数据

# ============ 准备数据 ============
otu <- as(otu_table(ps), "matrix")       # 提取OTU表
if (!taxa_are_rows(ps)) otu <- t(otu)    # 确保物种在行
otu_t <- t(otu)                          # 转置（样本在行）
meta <- as(sample_data(ps), "data.frame")  # 元数据

# 计算Bray-Curtis距离矩阵
bc_dist <- vegdist(otu_t, method = "bray")  # Bray-Curtis距离

# ============ 1. 基本PERMANOVA ============
# 检验疾病组和对照组的微生物组是否有显著差异
perm1 <- adonis2(
  bc_dist ~ Group,                       # 距离 ~ 分组
  data = meta,                           # 元数据
  permutations = 999                     # 999次置换
)
print(perm1)
# 关键输出：
# - Df: 自由度
# - SumOfSqs: 平方和
# - R2: 效应量（解释的方差比例）
# - F: F统计量
# - Pr(>F): p值（基于置换）

# ============ 2. 多因素PERMANOVA ============
# 控制混杂因素
perm2 <- adonis2(
  bc_dist ~ Group + Age + BMI + Sex,     # 多个因素
  data = meta,
  permutations = 999,
  by = "margin"                          # 边际检验（推荐）
  # by="margin": 每个因素控制其他因素后的独立效应
  # by="terms": 顺序检验（结果受变量顺序影响）
)
print(perm2)

# ============ 3. 嵌套设计PERMANOVA ============
# 例如：多个医院的患者，患者嵌套在医院内
perm3 <- adonis2(
  bc_dist ~ Hospital / Group,            # 医院/分组（嵌套）
  data = meta,
  permutations = 999,
  strata = meta$Hospital                 # 限制置换在医院内进行
)
print(perm3)

# ============ 4. 配对PERMANOVA ============
# 纵向研究：同一个人治疗前后
perm4 <- adonis2(
  bc_dist ~ Timepoint,                   # 时间点
  data = meta,
  permutations = 999,
  strata = meta$Subject_ID              # 限制置换在个体内
)
print(perm4)

# ============ 5. 方差齐性检验（betadisper） ============
# PERMANOVA假设组间方差齐性——必须检验！
bd <- betadisper(bc_dist, meta$Group)    # 方差齐性检验
bd_test <- permutest(bd, permutations = 999)  # 置换检验
print(bd_test)
# 如果显著 → 组间方差不齐 → PERMANOVA结果可能有偏
# 需要报告betadisper结果作为补充

# 可视化方差齐性
plot(bd, main = "方差齐性检验 (betadisper)")  # 绘图
boxplot(bd, main = "到组质心的距离")           # 箱线图

三、约束排序（dbRDA/CCA/RDA）¶

# 约束排序分析
library(vegan)        # 核心包

# ============ 1. dbRDA（距离基冗余分析） ============
# 适用于任何距离矩阵，最推荐的约束排序方法
dbrda_result <- dbrda(
  bc_dist ~ Group + Age + BMI,           # 距离 ~ 环境因素
  data = meta                            # 元数据
)
summary(dbrda_result)                    # 查看结果

# 提取解释的方差
cat("约束轴解释的总方差:",
    sum(dbrda_result$CCA$eig) /
    dbrda_result$tot.chi * 100, "%\n")   # 约束轴的方差比例

# 显著性检验
anova_dbrda <- anova(dbrda_result, by = "axis", permutations = 999)
print(anova_dbrda)                       # 每个轴的显著性

anova_terms <- anova(dbrda_result, by = "terms", permutations = 999)
print(anova_terms)                       # 每个因素的显著性

# dbRDA可视化
plot(dbrda_result, display = c("sites", "bp"),  # 样本点+环境箭头
     type = "n")                         # 先画空图
points(dbrda_result, display = "sites",
       col = as.numeric(factor(meta$Group)),  # 按组着色
       pch = 19, cex = 1.2)             # 点样式
text(dbrda_result, display = "bp",
     col = "blue", cex = 0.8)           # 环境变量箭头
legend("topright", legend = levels(factor(meta$Group)),
       col = 1:2, pch = 19)             # 图例

# ============ 2. CCA（典型对应分析） ============
# 适用于物种丰度数据（单峰响应）
# 注意：CCA需要原始计数或相对丰度，不用距离矩阵
cca_result <- cca(
  otu_t ~ Group + Age + BMI,            # OTU表 ~ 环境因素
  data = meta                            # 元数据
)
summary(cca_result)

# CCA显著性检验
anova(cca_result, by = "terms", permutations = 999)

# ============ 3. RDA（冗余分析） ============
# 适用于线性响应的数据（如Hellinger变换后）
otu_hell <- decostand(otu_t, method = "hellinger")  # Hellinger变换
rda_result <- rda(
  otu_hell ~ Group + Age + BMI,          # 变换后数据 ~ 环境
  data = meta
)
summary(rda_result)

# 前向选择最重要的环境变量
# 从空模型开始，逐步添加变量
rda0 <- rda(otu_hell ~ 1, data = meta)  # 空模型
rda_full <- rda(otu_hell ~ ., data = meta)  # 全模型
fwd_sel <- ordiR2step(
  rda0, scope = formula(rda_full),       # 搜索范围
  direction = "forward",                 # 前向选择
  permutations = 999                     # 置换检验
)
print(fwd_sel$anova)                     # 选择的变量和顺序

# ============ 4. envfit（环境向量拟合） ============
# 将环境变量拟合到非约束排序空间
pcoa <- cmdscale(bc_dist, k = 2, eig = TRUE)  # PCoA
env_fit <- envfit(
  pcoa, meta[, c("Age", "BMI", "Shannon", "pH")],  # 环境变量
  permutations = 999                     # 置换检验
)
print(env_fit)
# 输出R²和p值：哪些环境变量与群落结构显著相关

# envfit可视化
plot(pcoa$points, pch = 19,
     col = as.numeric(factor(meta$Group)))  # PCoA散点
plot(env_fit, p.max = 0.05)              # 只画显著的箭头

四、高级分析与可视化¶

# Python实现PERMANOVA和排序可视化
import numpy as np   # 数值计算
import pandas as pd  # 数据处理
from skbio.stats.distance import permanova  # PERMANOVA
from skbio import DistanceMatrix  # 距离矩阵
import matplotlib.pyplot as plt  # 绑图

# 1. 使用scikit-bio做PERMANOVA
def run_permanova(distance_matrix, metadata, group_col, permutations=999):
    """运行PERMANOVA检验"""
    # 创建DistanceMatrix对象
    dm = DistanceMatrix(distance_matrix,
                        ids=metadata.index.tolist())  # 距离矩阵
    # 运行PERMANOVA
    result = permanova(
        dm,
        metadata[group_col],              # 分组
        permutations=permutations         # 置换次数
    )
    print(f"PERMANOVA结果:")
    print(f"  检验统计量: {result['test statistic']:.4f}")
    print(f"  p值: {result['p-value']:.4f}")
    print(f"  样本量: {result['sample size']}")
    return result

# 2. PCoA + PERMANOVA结合可视化
def plot_pcoa_with_permanova(distance_matrix, metadata, group_col,
                              output="pcoa_permanova.png"):
    """PCoA可视化 + PERMANOVA检验结果"""
    from sklearn.manifold import MDS  # 多维缩放

    # MDS降维（等同于PCoA）
    mds = MDS(n_components=2, dissimilarity="precomputed",
              random_state=42)  # 2维
    coords = mds.fit_transform(distance_matrix)  # 降维

    # 运行PERMANOVA
    dm = DistanceMatrix(distance_matrix,
                        ids=metadata.index.tolist())
    perm_result = permanova(dm, metadata[group_col],
                            permutations=999)

    # 可视化
    fig, ax = plt.subplots(figsize=(10, 8))
    groups = metadata[group_col].unique()  # 所有组
    colors = plt.cm.Set1(np.linspace(0, 1, len(groups)))  # 颜色

    for i, group in enumerate(groups):
        mask = metadata[group_col] == group  # 该组样本
        ax.scatter(coords[mask, 0], coords[mask, 1],
                  c=[colors[i]], label=group,
                  s=80, alpha=0.7, edgecolors="black")  # 散点

        # 画置信椭圆
        from matplotlib.patches import Ellipse
        cov = np.cov(coords[mask, 0], coords[mask, 1])  # 协方差
        vals, vecs = np.linalg.eigh(cov)  # 特征值分解
        angle = np.degrees(np.arctan2(vecs[1, 1], vecs[0, 1]))
        width, height = 2 * 2 * np.sqrt(vals)  # 95%椭圆
        ellip = Ellipse(xy=(coords[mask, 0].mean(), coords[mask, 1].mean()),
                        width=width, height=height, angle=angle,
                        facecolor=colors[i], alpha=0.15)
        ax.add_patch(ellip)

    # 添加PERMANOVA结果
    r2 = perm_result['test statistic']   # 实际上是pseudo-F
    p_val = perm_result['p-value']       # p值
    ax.text(0.02, 0.98,
            f"PERMANOVA: p={p_val:.3f}",
            transform=ax.transAxes, fontsize=11,
            verticalalignment='top',
            bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))

    ax.set_xlabel("PCoA1")               # x轴
    ax.set_ylabel("PCoA2")               # y轴
    ax.set_title("PCoA + PERMANOVA")     # 标题
    ax.legend(loc="lower right")          # 图例
    plt.tight_layout()
    plt.savefig(output, dpi=150)
    print(f"图已保存: {output}")

# 3. 变异分解（Variation Partitioning）
def variation_partitioning_plot():
    """变异分解——将总变异分解为各因素的贡献"""
    # 在R中用vegan::varpart()实现：
    # varpart_result <- varpart(otu_t, ~ Group, ~ Age + BMI, ~ Diet, data=meta)
    # plot(varpart_result)
    #
    # Python中可用Venn图表示各因素解释的方差比例
    print("变异分解示例（需在R中运行）:")
    print("  varpart(Y, ~Group, ~Age+BMI, ~Diet)")
    print("  → 各因素独立解释 + 共同解释的方差比例")

五、方法选择决策树¶

你的数据是什么？
│
├─ 距离矩阵（Bray-Curtis/UniFrac等）
│   ├─ 检验组间差异 → PERMANOVA (adonis2)
│   ├─ 约束排序 → dbRDA
│   └─ 检验方差齐性 → betadisper
│
├─ 物种丰度表（原始计数/相对丰度）
│   ├─ 梯度长度>4 SD → CCA（单峰响应）
│   ├─ 梯度长度<3 SD → RDA（线性响应）
│   └─ 不确定 → 用DCA先判断梯度长度
│
└─ 要做什么？
    ├─ 检验假设（组间有差异？）→ PERMANOVA
    ├─ 找驱动因素 → dbRDA/CCA + envfit
    ├─ 变量选择 → ordiR2step（前向选择）
    └─ 分解变异来源 → varpart（变异分解）

常见报错与解决¶

报错	原因	解决方案
adonis2报"negative eigenvalues"	距离矩阵不满足欧几里得性质	用`sqrt`或`cailliez`校正，或改用dbRDA
PERMANOVA显著但PCoA图上看不出差异	效应量小但样本量大	报告R²（效应量），关注实际差异大小
betadisper显著	组间方差不齐	PERMANOVA结果需谨慎解读，补充报告
CCA轴特征值为负	数据不适合CCA	改用dbRDA或对数据做Hellinger变换后用RDA
envfit所有变量都不显著	环境变量与群落结构无关或样本量太小	检查变量类型是否正确，增加样本量

速查表¶

# 多变量分析流程
1. 计算距离矩阵（Bray-Curtis/UniFrac）
2. PCoA/NMDS非约束排序（探索性分析）
3. PERMANOVA检验组间差异
4. betadisper检验方差齐性
5. dbRDA/CCA约束排序（找驱动因素）
6. envfit拟合环境变量
7. varpart变异分解

# PERMANOVA关键参数
permutations = 999  # 至少999次
by = "margin"       # 推荐边际检验
strata =            # 配对/嵌套设计必须设

# CCA vs RDA vs dbRDA选择
DCA梯度>4 SD → CCA
DCA梯度<3 SD → RDA (+ Hellinger变换)
有距离矩阵 → dbRDA（最通用）

# 结果报告标准
PERMANOVA: F值, R², p值
约束排序: 约束轴解释方差%, 各因素显著性
envfit: R², p值

面试高频问题¶

Q1：PERMANOVA和ANOSIM有什么区别？ A：PERMANOVA基于距离矩阵的方差分解（类似ANOVA），直接输出效应量R²，支持多因素和协变量；ANOSIM基于距离排名比较组间vs组内距离，只支持单因素，没有效应量。PERMANOVA更强大、更灵活，现在基本替代了ANOSIM。但两者都假设方差齐性，需用betadisper检验。

Q2：PERMANOVA的by="margin"和by="terms"有什么区别？ A：by="terms"按公式中变量的顺序依次检验，前面的变量会影响后面的结果（类似Type I SS）；by="margin"控制其他所有变量后检验每个变量的独立效应（类似Type III SS）。实际研究中推荐用by="margin"，因为结果不受变量顺序影响。

Q3：dbRDA、CCA、RDA怎么选？ A：关键看数据类型和响应关系：(1) 有距离矩阵→用dbRDA（最通用）；(2) 原始丰度数据+单峰响应→CCA；(3) 原始丰度数据+线性响应→RDA（需Hellinger变换）。可以先用DCA看梯度长度：>4 SD用CCA，<3 SD用RDA。不确定时用dbRDA最安全。

Q4：PERMANOVA显著但PCoA图上两组重叠很大，怎么解释？ A：这是正常的——(1) PERMANOVA检测的是质心（中心位置）差异，即使分布重叠，质心不同就能显著；(2) 看R²，如果R²=0.02说明只解释了2%的变异，大部分变异来自个体差异；(3) 也可能是betadisper不齐导致的假阳性——一定要检验方差齐性。