一、PERMDISP 结果逐段解读（中文）
1. 样本信息
[INFO] Samples used for PERMDISP: 225
1:78 2:21 3:69 4:10 5:33 negative control:14与 PERMANOVA 完全相同的 225 个样本（5 个患者组 + NTC），组间样本量不均衡——这正是需要做 PERMDISP 的原因。
2. 全局 PERMDISP：显著 → 各组离散度不齐
- F(5, 219) = 15.198，置换检验 p = 1×10⁻⁴（9,999 次置换）。
- 含义：至少有一组的”组内离散度”（样本到本组中心点的平均距离）与其他组不同，即”多变量离散度同质性”假设被违反。
- 这直接回应了 Nicole 的顾虑：由于设计不均衡，PERMANOVA 的显著性不能自动全部归因于组成（location）差异，必须结合下面的成对结果逐对判断。
3. 各组平均离散度（到中心点距离）
| 组 | 均值 | SD | 解读 |
|---|---|---|---|
| 5 | 0.5696 | 0.047 | 离散度最高，且组内非常一致（SD 小）→ 该组样本普遍”彼此不像” |
| 3 | 0.5315 | 0.094 | 第二高，组内异质性强 |
| 1 | 0.4720 | 0.113 | 中等 |
| 4 | 0.4477 | 0.107 | 中等偏低 |
| 2 | 0.4378 | 0.084 | 较低 |
| NTC | 0.3269 | 0.178 | 均值最低（阴性对照背景一致），但 SD 最大 → 内部混合了”很紧”的正常 NTC 和”很散”的异常 NTC（NTC_3/6/13） |
4. 成对 PERMDISP（BH 校正后）
不显著（离散度齐性，4 对）：1–2 (p_BH=0.229)、1–4 (0.555)、2–4 (0.790)、4–NTC (0.093)。 → 这 4 对的 PERMANOVA 若显著，可以干净地解释为群落组成（位置）差异。
显著（离散度不齐，11 对）：1–3、1–5、1–NTC、2–3、2–5、2–NTC、3–4、3–5、3–NTC、4–5、5–NTC（p_BH = 0.0003–0.045）。 → 这些对的 PERMANOVA 显著性可能部分由离散度差异驱动，需在论文中注明谨慎解释，并用 taxon 水平方法（DESeq2 / ANCOM-BC2）交叉验证。 → 规律：几乎所有显著对都涉及 Group 3、Group 5 或 NTC（即”高异质组 vs 低异质组”的对比）。
5. 一句话结论
全局与成对 PERMDISP 表明组间离散度确实不齐（尤其 Group 3/5 vs Group 2/4/NTC），因此 PERMANOVA 中涉及这些组的比较应表述为”组成差异可能伴随离散度差异”；而 1–2、1–4、2–4、4–NTC 四对则可作纯组成差异解释。
二、为什么有两个”全局 PERMDISP”？保留哪一个？
anova(disp_glob) |
permutest(disp_glob, 9999) |
|
|---|---|---|
| 原理 | 经典参数 ANOVA，用 F 分布查表求 p | 置换检验，随机打乱组标签 9,999 次求 p |
| 假设 | 要求残差独立、正态 | 无分布假设 |
| 问题 | “到中心点的距离”由同一距离矩阵算出，彼此不独立、也常非正态 → p 值偏乐观（此处 8.15×10⁻¹³ 过小） | 无此问题，是 vegan / Anderson (2006) 推荐做法 |
| F 值 | 15.198（相同） | 15.198（相同） |
结论：只保留置换检验（permutest）。两个输出的 F 值完全一样，区别仅在 p 值的来源；置换 p 值才有效、可报告（写作 p = 1×10⁻⁴, 9,999 permutations；置换检验的 p 下限为 1/10000）。
修改后的代码（删掉 ANOVA 表）
# --- 3) GLOBAL PERMDISP (betadisper + permutation test ONLY) ---
disp_glob <- betadisper(bc_dist, group = grp) # type = "centroid"
perm_glob <- permutest(disp_glob, permutations = 9999) # 唯一的全局检验
cat("```text\n")
cat("[Global PERMDISP - permutation test, 9999 perms]\n")
print(perm_glob)
cat("```\n")
cat("\n[Mean distance to group centroid (multivariate dispersion) per group]\n")
disp_mean <- tapply(disp_glob$distances, grp, mean)
disp_sd <- tapply(disp_glob$distances, grp, sd)
disp_tab <- data.frame(Group = names(disp_mean),
Mean_dist_to_centroid = round(as.numeric(disp_mean), 4),
SD_dist = round(as.numeric(disp_sd), 4))
rownames(disp_tab) <- NULL
print(disp_tab)即：删除原来的 print(anova(disp_glob)) 及其标题行，其余不变。报告中统一写：Global PERMDISP: F(5,219) = 15.20, p = 1×10⁻⁴ (9,999 permutations)。
一、PERMDISP（多变量离散度同质性检验）详解
注：正确拼写为 PERMDISP（PERMutational analysis of multivariate DISPersions），在 vegan 包中由
betadisper()实现，也称 “homogeneity of multivariate dispersions test”。
1. 为什么 PERMANOVA 之后还必须做 PERMDISP？
PERMANOVA（adonis2）检验的是：各组的”组中心点”（centroid，即多元空间中的位置/location）是否不同。但它的伪 F 统计量对另一种差异同样敏感——组内离散度（dispersion，即组内样本彼此之间的变异程度）。
也就是说，一个显著的 PERMANOVA 结果可能有两种来源：
- 位置效应（location effect）：组的平均群落组成确实不同（这是我们想要的生物学结论）；
- 离散效应（dispersion effect）：组的平均组成可能相同，但某一组内部样本间差异极大（例如某组混入了高度异质的样本），导致组间距离被拉大。
当各组样本量不均衡时（您的数据：Group1=78、Group2=21、Group3=69、Group4=10、Group5=33、NTC=14），PERMANOVA 对离散度差异尤其敏感，这正是 Nicole 提出此顾虑的原因。PERMDISP 的作用就是把这两种效应拆开。
2. PERMDISP 的原理（逐步）
- 计算距离矩阵：与 PERMANOVA 完全相同的矩阵（Hellinger 转化后的 Bray–Curtis 距离）。
- 求每组的中心点：在多元主坐标空间（PCoA 空间）中计算每组的 centroid（默认）或 spatial median（
type = "median"，对离群值更稳健）。 - 计算每个样本到其组中心点的距离
d_i：这个距离就是该样本的”多变量离散度”。组内平均距离大 = 该组内部异质性强。 - 比较各组的平均距离：形式上相当于对
d_i做单因素 ANOVA，但由于距离之间不独立、分布非正态，不使用 F 分布查表，而是用置换检验（permutation test，通常 9,999 次） 得到 p 值。
3. 全局检验 + 成对检验
- 全局 PERMDISP：H₀ = 所有组的离散度相同。
betadisper(bc_dist, group)→permutest(disp, permutations = 9999)或anova(disp)。
- 成对 PERMDISP：对每一组两两组合（6 个组 → 15 对），只取这两组的样本重新做 betadisper + 置换检验，再对 15 个 p 值做 BH 校正。
- 同时报告每组的平均 centroid 距离，便于描述哪一组内部更异质。
4. 结果解读矩阵
| PERMANOVA | PERMDISP | 结论 |
|---|---|---|
| 显著 | 不显著 | 组间差异来自群落组成（位置），结论可靠 ✅ |
| 显著 | 显著 | 差异可能部分/全部来自离散度，需谨慎解释，并在论文中注明受影响的比较 ⚠️ |
| 不显著 | 显著 | 组成均值相似但变异性不同（本身也是有意义的生物学发现） |
二、ANCOM-BC2 结合流行率过滤（Prevalence Filtering）——最详细解释
1. 根本问题：16S 数据是”组成型数据”（compositional data）
测序得到的不是绝对菌量，而是相对丰度（比例）。数学上：
观测reads数 O_ij = 真实绝对丰度 A_ij × 采样分数 N_jN_j（样本 j 的总测序深度/采样分数）未知且与生物学无关；- 每个样本的所有比例之和恒等于 1（sum constraint）→ 一个 taxon 比例升高会”挤压”其他所有 taxon 的比例，产生虚假负相关；
- 后果：某 taxon 看起来”上调”，可能仅仅是因为另一个 taxon 大幅减少，而非它自己真的增多。
2. 为什么 DESeq2 不够
- DESeq2 源自 RNA-seq：其 size factor（median-of-ratios）假设”大多数基因不变”，且 counts 反映绝对表达量；微生物组中该假设常被违反（群落结构整体位移很常见）。
- 微生物组数据高度稀疏（多数 taxon 在多数样本中为 0）。仅出现在 1–2 个样本中的 taxon 会产生极大的 log2FoldChange（由单个样本驱动）——这正是 Nicole 说 “some of the DESeq2 fold changes are very large and may be driven by sparse taxa” 的原因。
- DESeq2 不建模组成型约束。因此需要一个组成型方法（compositional method） 作为互补/验证。
3. ANCOM 家族的演进
| 方法 | 核心思想 | 局限 |
|---|---|---|
| ANCOM (Mandal 2015) | 对每对 taxon 做 log-ratio 检验（W 统计量） | 保守、计算慢、无 FDR 控制 |
| ANCOM-BC (Lin & Peddada 2020) | 直接估计并校正采样分数（bias correction），在”近似绝对丰度”尺度上做检验 | 对零值、小样本、协变量支持有限 |
| ANCOM-BC2 (Lin & Peddada) | 统一线性回归框架 + 改进的零值处理 + 正规 FDR 控制 + 可含协变量/连续变量 + 结构零检验 | 仍需配合流行率过滤使用效果最佳 |
4. ANCOM-BC2 的偏倚校正核心思想
模型（log 尺度）：
log(O_ij) = log(N_j) + log(A_ij) + ε_ij- ANCOM-BC2 通过约束优化/迭代算法估计每个样本的未知偏倚项 log(N_j)（假设：大多数 taxon 在组间无差异——与 DESeq2 的假设类似，但实现于 log-线性模型中）；
- 去除偏倚后，对每个 taxon 拟合线性回归：
y = μ + β·group + ε，检验β = 0； - 输出的 lfc（log fold change）位于偏倚校正后的尺度（近似绝对丰度尺度），因此比相对丰度的 fold change 更不易受组成效应扭曲；
- p 值经多重检验校正得到 q 值（FDR），q < 0.05 判为显著。
5. 零值（zeros）处理
- 抽样零（sampling zero）：菌真实存在但因深度不够未检出；
- 结构零（structural zero）：该菌在某组中真实不存在；
- ANCOM-BC2 通过
zero_cut参数（默认 0.90）识别结构零：若某 taxon 在某一组中零的比例超过阈值，则视为结构零，单独做”存在/缺失”差异检验；其余零通过伪计数（pseudo-count）策略处理，保证 log 运算可行。
6. 流行率过滤（Prevalence Filtering）——本次要求的核心
定义：taxon i 的流行率 = （检出该 taxon（count > 0）的样本数）÷（总样本数）。
在您的流程中的具体做法：
- 仅在患者样本（Group1–5 的并集，不含 NTC/PC）上计算流行率——因为对照的污染谱与患者不同，混入会扭曲过滤；
- 主分析：保留流行率 ≥ 10% 的 taxon（
prev_cutoff = 0.10，或手动rowMeans(otu > 0) >= 0.10后prune_taxa）； - 敏感性分析：保留流行率 ≥ 20% 的 taxon。
为什么要过滤（四个理由）：
- 稳定性：只出现在 1–2 个样本中的稀疏 taxon 方差估计极不稳定，fold change 由个别样本驱动（直接回应 Nicole 的顾虑）；
- 统计功效：减少被检验的 taxon 数 → 减轻多重检验负担 → 提高剩余 taxon 的检出功效；
- 去伪信号：低流行率 taxon 更可能是污染、index hopping 或测序错误的产物；
- 模型收敛与 FDR 控制的可靠性：稀疏 taxon 会破坏回归残差假设。
为什么还要做 20% 的敏感性分析：
- 10% 是”宽松”阈值，20% 是”严格”阈值；
- 在两个阈值下都显著的 taxon = 稳健（robust）taxon；
- 若两阈值结果差异巨大，说明结论对过滤阈值敏感，需在论文中如实说明；
- 过滤过严可能丢失”真实但罕见”的重要 taxon，因此双阈值互为印证。
7. 输出与解读
ancombc2() 返回的 res 表每个 taxon 一行，关键列：
lfc_GroupX：组 X 相对参考组（ref_level）的 log fold change（校正后尺度）；p_GroupX/q_GroupX：p 值与 FDR 校正 q 值；q < 0.05 判显著；- 方向：lfc 符号 = 上调/下调（与您 DESeq2 中
relevel的参考组保持一致，便于直接比较方向）。
与 DESeq2 的一致性（concordance）评估：
- 对每个两两比较：n(ANCOM-BC2 显著)、n(DESeq2 显著)、交集大小、方向一致的比例；
- 稳健清单 = ANCOM-BC2(10%) ∩ ANCOM-BC2(20%) ∩ DESeq2 三者交集；
- 若某 taxon 仅 DESeq2 显著而 ANCOM-BC2 不显著 → 提示其可能为稀疏/组成效应驱动的假阳性；反之则提示 DESeq2 功效不足。
8. 对应到您 Rmd 代码的每一步
# 1) 在患者样本上计算流行率
prev_patient <- rowMeans(otu_mat[, patient_ids, drop = FALSE] > 0)
ps_prev10 <- prune_taxa(prev_patient >= 0.10, ps_filt) # 主分析
ps_prev20 <- prune_taxa(prev_patient >= 0.20, ps_filt) # 敏感性分析
# 2) 每个两两比较：ref_level 与 DESeq2 的 relevel 保持一致
out <- ancombc2(data = sub, formula = "Group", ref_level = g2, alpha = 0.05)
# 3) 10 个比较 × 2 个阈值 = 20 次运行；导出 xlsx/csv
# 4) 一致性表 + 稳健 taxon 表（三交集）三、一句话总结
- PERMDISP：检验”组内变异是否齐性”，用以证明 PERMANOVA 的显著性反映的是群落组成差异而非离散度差异；
- ANCOM-BC2 + 流行率过滤：在校正测序偏倚与组成效应的前提下做差异丰度检验，先用 10%/20% 流行率剔除由稀疏样本驱动的假阳性 taxon，再与 DESeq2 结果交叉验证，得到稳健的差异 taxon 清单。