微生物组关联分析MWAS（菌群-表型关联分析）¶

一句话概述¶

微生物组全关联分析（Microbiome-Wide Association Study, MWAS）是借鉴GWAS思路，系统性地检测肠道/口腔/皮肤等部位微生物群落组成或功能与宿主表型（疾病、代谢指标、药物反应等）之间关联关系的分析方法，核心工具包括MaAsLin2（多变量回归）、PERMANOVA（群落整体差异）和多种统计建模方法。

核心知识点表格¶

知识点	说明
MWAS	Microbiome-Wide Association Study，微生物组全关联研究
MaAsLin2	Microbiome Multivariable Associations with Linear Models 2
PERMANOVA	Permutational Multivariate Analysis of Variance
Alpha多样性	样本内多样性（Shannon/Simpson/Chao1等）
Beta多样性	样本间多样性（Bray-Curtis/UniFrac/Jaccard等）
差异丰度分析	检测组间差异菌种（LEfSe/DESeq2/ANCOM-BC）
组成数据	Compositional data，相对丰度之和为1的约束数据
CLR变换	Centered Log-Ratio，处理组成数据的标准方法
零膨胀	Zero-inflation，微生物组数据中大量零值
混杂因素	年龄、BMI、饮食、药物等影响菌群的混杂变量
多重检验	同时检测数百个物种需要严格的FDR校正
功能预测	PICRUSt2/HUMAnN3预测菌群功能通路

各步骤详解¶

第一步：微生物组关联分析概念与设计¶

白话解释： MWAS就是对微生物组做"关联研究"——检查成百上千种微生物中哪些与你关心的健康问题（如肥胖、糖尿病、癌症）有统计学关联。这类似于GWAS检查基因组中哪些变异与疾病相关，但对象从基因换成了微生物。挑战在于微生物组数据很"特殊"——相对丰度有组成性约束、大量零值、高度相关等。

技术细节：

微生物组数据特点： 1. 组成性（Compositional）：相对丰度之和为1/100%，一个菌增加必然其他菌减少 2. 零膨胀（Zero-inflated）：大量物种在很多样本中检测不到（0值） 3. 过度离散（Over-dispersed）：方差远大于均值 4. 高维性：物种数远多于样本数 5. 层次结构：从门到种有分类层次 6. 高度相关：共生/竞争关系导致菌种间强相关

研究设计类型： | 设计 | 优点 | 缺点 | 适用 | |------|------|------|------| | 横断面 | 简单快速 | 不能推因果 | 初步探索 | | 纵向 | 可追踪变化 | 费时费力 | 因果推断 | | 干预实验 | 可验证因果 | 最复杂 | 明确机制 | | 病例对照 | 便于回顾 | 回忆偏倚 | 疾病关联 |

第二步：数据预处理¶

白话解释： 16S扩增子测序或宏基因组测序产生的数据需要经过质控、去噪/聚类（得到ASV/OTU表）、物种注释等步骤，最后得到一张"物种×样本"的丰度表。然后需要对这张表做过滤和变换，让它适合统计分析。

# ===== 16S数据预处理（DADA2 + phyloseq） =====
library(dada2)
library(phyloseq)
library(microbiome)

# 假设已通过DADA2得到ASV表、分类表、样本信息
# 创建phyloseq对象
ps <- phyloseq(
  otu_table(asv_table, taxa_are_rows = TRUE),
  tax_table(taxonomy_table),
  sample_data(sample_metadata),
  phy_tree(tree)  # 如有
)

# ===== 质量过滤 =====
# 1. 去除低丰度物种
# 保留至少在10%样本中出现且相对丰度>0.01%的物种
ps_filtered <- filter_taxa(ps, function(x) {
  sum(x > 0) > 0.1 * nsamples(ps)  # 至少10%样本非零
}, prune = TRUE)

# 或使用更严格的过滤
ps_filtered <- prune_taxa(taxa_sums(ps) > 10, ps)  # 总丰度>10

# 2. 去除低测序深度样本
ps_filtered <- prune_samples(sample_sums(ps_filtered) > 1000, ps_filtered)

# 3. 检查数据概览
cat("Samples:", nsamples(ps_filtered), "\n")
cat("Taxa:", ntaxa(ps_filtered), "\n")
cat("Min depth:", min(sample_sums(ps_filtered)), "\n")
cat("Max depth:", max(sample_sums(ps_filtered)), "\n")

# ===== 数据变换 =====
# 方法1：相对丰度（传统，但有组成性问题）
ps_rel <- transform_sample_counts(ps_filtered, function(x) x / sum(x))

# 方法2：CLR变换（推荐，处理组成性）
library(compositions)
# 需要处理零值（加伪计数或使用czm方法）
asv_clr <- microbiome::transform(ps_filtered, "clr")

# 方法3：稀释（rarefaction）到相同深度
ps_rarefied <- rarefy_even_depth(ps_filtered, sample.size = min(sample_sums(ps_filtered)),
                                  rngseed = 42, replace = FALSE)

# 方法4：Variance Stabilizing Transformation (DESeq2)
library(DESeq2)
dds <- phyloseq_to_deseq2(ps_filtered, ~ Group)
dds <- estimateSizeFactors(dds, type = "poscounts")
vst_data <- varianceStabilizingTransformation(dds)

第三步：Alpha和Beta多样性分析¶

白话解释： Alpha多样性衡量"单个样本中微生物的丰富度和均匀度"——是否种类多且分布均匀？Beta多样性衡量"样本之间的群落差异"——两个人的菌群到底有多不同？这是微生物组分析的基础描述性步骤。

# ===== Alpha多样性 =====
library(phyloseq)
library(ggpubr)

# 计算多种Alpha多样性指标
alpha_div <- estimate_richness(ps_rarefied, 
                                measures = c("Observed", "Shannon", "Simpson", "Chao1"))
alpha_div$Group <- sample_data(ps_rarefied)$Group

# 组间比较
# Shannon指数
ggboxplot(alpha_div, x = "Group", y = "Shannon", fill = "Group") +
  stat_compare_means(method = "wilcox.test") +
  labs(title = "Shannon Diversity")

# Chao1丰富度
ggboxplot(alpha_div, x = "Group", y = "Chao1", fill = "Group") +
  stat_compare_means(method = "wilcox.test") +
  labs(title = "Chao1 Richness")

# 多组比较用Kruskal-Wallis
kruskal.test(Shannon ~ Group, data = alpha_div)
# 事后两两比较
pairwise.wilcox.test(alpha_div$Shannon, alpha_div$Group, p.adjust.method = "BH")

# ===== Beta多样性 =====
library(vegan)

# 计算距离矩阵
# Bray-Curtis（基于丰度）
bray_dist <- distance(ps_rel, method = "bray")
# Weighted UniFrac（基于系统发育+丰度）
wunifrac_dist <- distance(ps_rarefied, method = "wunifrac")
# Unweighted UniFrac（基于系统发育+有无）
uunifrac_dist <- distance(ps_rarefied, method = "uunifrac")

# PCoA排序
pcoa_bray <- ordinate(ps_rel, method = "PCoA", distance = bray_dist)
plot_ordination(ps_rel, pcoa_bray, color = "Group") +
  stat_ellipse(level = 0.95) +
  theme_classic() +
  labs(title = "PCoA (Bray-Curtis)")

# NMDS排序
nmds <- ordinate(ps_rel, method = "NMDS", distance = bray_dist)
plot_ordination(ps_rel, nmds, color = "Group") +
  labs(title = paste("NMDS, stress =", round(nmds$stress, 3)))

# ===== PERMANOVA =====
# 检验组间Beta多样性是否显著不同
metadata <- as(sample_data(ps_rel), "data.frame")

# 基本PERMANOVA
permanova_result <- adonis2(
  bray_dist ~ Group,
  data = metadata,
  permutations = 999
)
print(permanova_result)
# R2 = 0.15, P = 0.001 → 分组解释15%的变异，差异显著

# 加入混杂变量的PERMANOVA
permanova_adj <- adonis2(
  bray_dist ~ Age + BMI + Sex + Group,
  data = metadata,
  permutations = 999,
  by = "margin"  # 每个变量独立检验
)
print(permanova_adj)

# 方差齐性检验（PERMANOVA的前提假设）
betadisper_result <- betadisper(bray_dist, metadata$Group)
permutest(betadisper_result)  # P > 0.05说明方差齐性满足

# ===== ANOSIM（另一种群落比较方法） =====
anosim_result <- anosim(bray_dist, metadata$Group, permutations = 999)
print(anosim_result)

第四步：MaAsLin2多变量关联分析¶

白话解释： MaAsLin2是专门为微生物组数据设计的多变量回归方法——它可以在控制混杂因素（年龄、BMI、性别、测序批次等）的情况下，找出哪些微生物与你关心的表型显著关联。相比简单的Wilcoxon检验，它更严谨、更能反映真实关联。

技术细节： MaAsLin2框架： - 支持多种归一化：TSS(相对丰度)、CSS、TMM、CLR - 支持多种变换：LOG、LOGIT、AST(arcsine square root) - 支持多种统计模型：LM(线性模型)、CPLM(复合泊松线性模型)、NEGBIN(负二项)、ZINB(零膨胀负二项) - 内置多重检验校正（BH FDR）

# ===== MaAsLin2 =====
# BiocManager::install("Maaslin2")
library(Maaslin2)

# 准备数据
# 丰度表：样本为行，物种为列
abundance_df <- as.data.frame(t(otu_table(ps_filtered)))

# 元数据：样本为行，变量为列
metadata_df <- as.data.frame(sample_data(ps_filtered))

# 运行MaAsLin2
fit_data <- Maaslin2(
  input_data = abundance_df,
  input_metadata = metadata_df,
  output = "maaslin2_output",
  # 固定效应（要检验的变量 + 混杂变量）
  fixed_effects = c("Group", "Age", "BMI", "Sex"),
  # 随机效应（如纵向数据中的个体）
  # random_effects = c("SubjectID"),
  # 归一化方法
  normalization = "TSS",      # Total Sum Scaling
  # 数据变换
  transform = "LOG",          # log变换
  # 统计模型
  analysis_method = "LM",     # 线性模型
  # 过滤
  min_abundance = 0.0001,     # 最小丰度阈值（默认为0.0，此处设为0.01%）
  min_prevalence = 0.1,       # 最小普遍率（至少10%样本存在，默认0.1）
  # 多重检验
  correction = "BH",
  # 显著性
  max_significance = 0.25,        # q-value阈值（MaAsLin2默认0.25）
  # 参考组
  reference = c("Group,Control"),
  # 核心数
  cores = 4,
  # 图形输出
  plot_heatmap = TRUE,
  plot_scatter = TRUE
)

# 查看结果
results <- read.table("maaslin2_output/all_results.tsv", header = TRUE, sep = "\t")
sig_results <- results[results$qval < 0.25, ]
sig_results <- sig_results[order(sig_results$qval), ]
print(sig_results[, c("feature", "metadata", "coef", "stderr", "pval", "qval")])

# 火山图
library(ggplot2)
group_results <- results[results$metadata == "Group", ]
ggplot(group_results, aes(x = coef, y = -log10(qval))) +
  geom_point(aes(color = qval < 0.25), size = 2) +
  geom_hline(yintercept = -log10(0.25), linetype = "dashed") +
  geom_text_repel(data = subset(group_results, qval < 0.25),
                  aes(label = feature), size = 3) +
  scale_color_manual(values = c("grey", "red")) +
  theme_classic() +
  labs(x = "Coefficient (Effect Size)", y = "-log10(q-value)",
       title = "MaAsLin2: Group Effect on Microbiome")

第五步：差异丰度分析方法对比¶

白话解释： 除了MaAsLin2，还有很多方法可以找"差异菌"。LEfSe是最经典的、ANCOM-BC考虑了组成性偏倚、DESeq2从RNA-seq借用。每种方法有不同的假设和优缺点。

# ===== ANCOM-BC2（推荐，处理组成性偏倚） =====
# BiocManager::install("ANCOMBC")
library(ANCOMBC)

# 运行ANCOM-BC2
ancombc_result <- ancombc2(
  data = ps_filtered,
  fix_formula = "Group + Age + Sex",
  rand_formula = NULL,
  p_adj_method = "holm",
  prv_cut = 0.10,             # 普遍率过滤
  lib_cut = 1000,             # 最小文库大小
  group = "Group",
  struc_zero = TRUE,          # 检测结构性零值
  neg_lb = TRUE,
  alpha = 0.05,
  global = TRUE               # 全局检验
)

# 提取结果（ANCOM-BC2的结果在一个扁平data.frame中）
res <- ancombc_result$res
# diff_列为TRUE表示该taxon在该变量水平上差异显著
sig_taxa <- res[res$diff_GroupDisease == TRUE, ]

# ===== LEfSe（经典方法，适合快速筛选） =====
# 通常在Galaxy或命令行使用
# R中可用microbiomeMarker包
# BiocManager::install("microbiomeMarker")
library(microbiomeMarker)

lefse_result <- run_lefse(
  ps_filtered,
  group = "Group",
  norm = "CPM",
  lda_cutoff = 2,             # LDA score阈值
  wilcoxon_cutoff = 0.05,
  kw_cutoff = 0.05
)

# 绘制LEfSe条形图
plot_ef_bar(lefse_result)

# ===== DESeq2（适合counts数据） =====
library(DESeq2)

# phyloseq转DESeq2
dds <- phyloseq_to_deseq2(ps_filtered, ~ Group)
dds <- DESeq(dds, test = "Wald", fitType = "parametric")
deseq_res <- results(dds, contrast = c("Group", "Disease", "Control"))
deseq_sig <- deseq_res[which(deseq_res$padj < 0.05), ]

# ===== ALDEx2（贝叶斯方法，专为组成数据） =====
library(ALDEx2)

# 运行ALDEx2
aldex_result <- aldex(
  reads = as.data.frame(otu_table(ps_filtered)),
  conditions = sample_data(ps_filtered)$Group,
  mc.samples = 128,
  test = "t",
  effect = TRUE,
  include.sample.summary = FALSE,
  denom = "all",               # CLR变换分母
  verbose = FALSE
)

# 显著结果
aldex_sig <- aldex_result[aldex_result$wi.eBH < 0.05, ]

第六步：功能预测与通路分析¶

白话解释： 知道了哪些菌有变化还不够——我们更想知道这些菌群变化会带来什么功能上的影响。功能预测工具可以根据菌群的物种组成，预测它们携带哪些代谢通路和基因功能，从而理解菌群变化的生物学意义。

# ===== PICRUSt2（16S数据功能预测） =====
# 从ASV序列预测宏基因组功能
conda install -c bioconda picrust2

# 输入：ASV丰度表 + ASV代表序列
picrust2_pipeline.py \
  -s asv_seqs.fasta \
  -i asv_table.biom \
  -o picrust2_output/ \
  -p 16 \
  --stratified \
  --per_sequence_contrib

# 输出包含：
# - EC_metagenome_out/: 酶分类预测
# - KO_metagenome_out/: KEGG直系同源基因预测
# - pathways_out/: MetaCyc通路预测

# ===== HUMAnN3（宏基因组功能注释） =====
# 如果有shotgun宏基因组数据
conda install -c bioconda humann

# 运行HUMAnN3
humann --input metagenome_R1.fastq.gz \
       --output humann3_output/ \
       --threads 16 \
       --nucleotide-database chocophlan_db/ \
       --protein-database uniref_db/

# 合并多个样本
humann_join_tables --input humann3_output/ \
                   --output merged_pathabundance.tsv \
                   --file_name pathabundance

# 归一化
humann_renorm_table --input merged_pathabundance.tsv \
                    --output merged_pathabundance_relab.tsv \
                    --units relab

# ===== R中分析功能预测结果 =====
# 读取PICRUSt2 pathway结果
pathway_abundance <- read.table("picrust2_output/pathways_out/path_abun_unstrat.tsv",
                                 header = TRUE, sep = "\t", row.names = 1)

# 用MaAsLin2分析差异通路
Maaslin2(
  input_data = as.data.frame(t(pathway_abundance)),
  input_metadata = metadata_df,
  output = "pathway_maaslin2/",
  fixed_effects = c("Group", "Age", "BMI"),
  normalization = "TSS",
  transform = "LOG",
  analysis_method = "LM",
  min_abundance = 0.0001,
  min_prevalence = 0.1
)

第七步：多组学整合与因果推断¶

白话解释： 微生物组数据与宿主的代谢组、免疫指标、基因组等多组学数据整合，可以揭示更完整的菌群-宿主互作机制。此外，利用孟德尔随机化等方法，可以从关联推进到因果。

# ===== 微生物组+代谢组联合分析（MOFA2） =====
# BiocManager::install("MOFA2")
library(MOFA2)

# 准备多组学数据
# view1: 微生物组CLR丰度
# view2: 代谢组
# view3: 免疫指标

data_list <- list(
  microbiome = as.matrix(microbiome_clr),
  metabolome = as.matrix(metabolome_data),
  immune = as.matrix(immune_markers)
)

# 创建MOFA对象
mofa_obj <- create_mofa(data_list)
mofa_obj <- prepare_mofa(mofa_obj)
mofa_obj <- run_mofa(mofa_obj)

# 分析因子
plot_variance_explained(mofa_obj)
plot_factor(mofa_obj, factor = 1, color_by = "Group")
plot_weights(mofa_obj, view = "microbiome", factor = 1, nfeatures = 20)

# ===== Spearman相关网络 =====
# 微生物-代谢物相关性
library(Hmisc)

cor_result <- rcorr(
  as.matrix(cbind(microbiome_clr, metabolome_data)),
  type = "spearman"
)

# 过滤显著相关
sig_cors <- which(cor_result$P < 0.05 & abs(cor_result$r) > 0.3, arr.ind = TRUE)

# ===== 孟德尔随机化（MR） =====
# 利用微生物组的SNP（mQTL）作为工具变量
# 推断微生物→疾病的因果方向
# 使用TwoSampleMR包
library(TwoSampleMR)
# 具体实现需要GWAS summary statistics

实战命令（可复制）¶

完整MWAS分析流程¶

# ============================================
# 微生物组关联分析（MWAS）完整流程
# ============================================

# 1. 安装包
BiocManager::install(c("phyloseq", "microbiome", "Maaslin2", "ANCOMBC",
                        "vegan", "ggpubr", "ggrepel"))
library(phyloseq); library(microbiome); library(Maaslin2)
library(ANCOMBC); library(vegan); library(ggpubr)

# 2. 读入数据
ps <- readRDS("phyloseq_object.rds")  # 或从各组件创建

# 3. 数据过滤
ps <- prune_samples(sample_sums(ps) > 1000, ps)
ps <- prune_taxa(taxa_sums(ps) > 10, ps)
ps <- filter_taxa(ps, function(x) sum(x > 0) > 0.1 * nsamples(ps), prune = TRUE)

# 4. Alpha多样性
alpha <- estimate_richness(rarefy_even_depth(ps, rngseed = 42),
                           measures = c("Shannon", "Chao1", "Simpson"))
alpha$Group <- sample_data(ps)$Group
wilcox.test(Shannon ~ Group, data = alpha)

# 5. Beta多样性 + PERMANOVA
ps_rel <- transform_sample_counts(ps, function(x) x/sum(x))
bray <- distance(ps_rel, "bray")
adonis2(bray ~ Group + Age + BMI + Sex, data = as(sample_data(ps), "data.frame"))

# 6. MaAsLin2差异分析
Maaslin2(
  input_data = as.data.frame(t(otu_table(ps))),
  input_metadata = as.data.frame(sample_data(ps)),
  output = "mwas_results/",
  fixed_effects = c("Group", "Age", "BMI", "Sex"),
  normalization = "TSS", transform = "LOG",
  analysis_method = "LM",
  min_abundance = 0.0001, min_prevalence = 0.1,
  reference = c("Group,Control"),
  cores = 4
)

# 7. ANCOM-BC2验证
ancom_res <- ancombc2(ps, fix_formula = "Group + Age + Sex",
                       p_adj_method = "holm", prv_cut = 0.1)

# 8. 汇总一致性结果
maaslin_sig <- read.table("mwas_results/all_results.tsv", header = TRUE, sep = "\t")
maaslin_sig <- maaslin_sig[maaslin_sig$qval < 0.25 & maaslin_sig$metadata == "Group", "feature"]
ancom_sig <- ancom_res$res$taxon[ancom_res$res$diff_GroupDisease == TRUE]
consistent <- intersect(maaslin_sig, ancom_sig)
cat("Consistently significant taxa:", paste(consistent, collapse = ", "))

面试常问点¶

Q1: 为什么微生物组数据不能直接用t检验做差异分析？¶

A: 因为微生物组数据有几个特殊性质：(1) 组成性——相对丰度之和为1，增加一个菌的丰度必然减少其他菌，导致假的负相关（spurious correlations）；(2) 零膨胀——大量零值使正态分布假设不成立；(3) 过度离散——方差远大于均值；(4) 多重检验——同时测试几百个物种需要严格FDR校正；(5) 混杂因素——年龄、饮食、BMI等强烈影响菌群。需要使用专门的工具（MaAsLin2/ANCOM-BC/ALDEx2）处理这些问题。

Q2: MaAsLin2和ANCOM-BC的区别是什么？¶

A: MaAsLin2是通用多变量回归框架，灵活性高（可选多种归一化、变换、模型），但本质上是对变换后的丰度做线性回归，对组成性偏倚处理不够理论化。ANCOM-BC基于严格的组成数据统计理论，显式估计和校正"sampling fraction"偏差（不同样本测序深度的系统差异），统计理论更严格。实践中建议：用两种方法交叉验证，取一致结果更可靠。

Q3: PERMANOVA的R²值很小（如0.05）还有意义吗？¶

A: 有意义。在人类微生物组研究中，单个因素通常只能解释2-10%的变异（因为影响因素太多）。R²=0.05意味着该因素解释了5%的群落变异，如果P值显著（P<0.05），说明这个效应是真实的。但R²也意味着95%的变异由其他因素解释——这提醒我们要全面考虑混杂因素，且不要过度解释单一因素的作用。

Q4: 如何处理微生物组数据中的零值？¶

A: 三种策略：(1) 加伪计数：简单但引入偏倚，通常加0.5或最小非零值的一半；(2) 零膨胀模型：如ZINB（零膨胀负二项），显式建模零的来源（真零vs假零）；(3) 专用变换：如ANCOM系列使用log-ratio不需要处理零值（通过pair-wise ratio避开），ALDEx2用贝叶斯估计替代零。选择取决于零值比例和分析方法的要求。

Q5: 16S和宏基因组数据的MWAS分析有什么不同？¶

A: (1) 分辨率：16S到属/种级别，宏基因组可到种/株级别且直接获得功能信息；(2) 定量准确性：16S有PCR和引物偏倚，宏基因组相对更准确；(3) 功能分析：16S需要PICRUSt2预测（准确性有限），宏基因组直接用HUMAnN3注释；(4) 统计处理：基本相同，但宏基因组的counts数据可能更适合负二项模型。

Q6: 如何控制MWAS中的混杂因素？¶

A: (1) 在MaAsLin2中将混杂变量作为fixed_effects纳入模型；(2) 在PERMANOVA中将混杂变量先于目标变量进入模型（顺序敏感，或用"margin"独立检验）；(3) 匹配设计——病例和对照按年龄、性别、BMI等配对；(4) 排除正在使用抗生素/PPI的样本；(5) 记录并报告饮食信息。最常见的混杂因素：年龄、BMI、性别、饮食、药物（尤其PPI和抗生素）、地理位置、测序批次。

Q7: 如何从关联推进到因果？¶

A: 横断面关联不等于因果。推进方法：(1) 纵向研究：追踪菌群变化是否先于疾病发生；(2) 干预实验：补充/清除特定菌后观察表型变化；(3) 孟德尔随机化：用mQTL作为工具变量推断因果方向；(4) 粪菌移植实验：在动物模型中移植患者菌群看是否重现表型；(5) 中介分析：检验菌群是否是暴露→疾病路径中的中介变量。

易错点¶

1. 不做稀释就比较Alpha多样性¶

问题： 测序深度不同的样本直接比较Shannon/Chao1，高深度样本天然看起来多样性更高。 解决： Alpha多样性分析前必须稀释到相同深度（rarefaction），或使用对深度不敏感的指标（如coverage-based rarefaction）。

2. 把相对丰度当作绝对量解读¶

问题： 一个菌的相对丰度增加不一定说明它的绝对数量增加——可能只是其他菌减少了。 解决： 使用绝对定量方法（如spike-in标准品、qPCR总菌量校正），或用ANCOM-BC等考虑组成性偏倚的方法，或至少在讨论中注明这一局限。

3. PERMANOVA组间方差不齐却不报告¶

问题： PERMANOVA假设组间方差齐性。如果一组样本的组内多样性显著大于另一组（dispersion不同），PERMANOVA的P值不可靠。 解决： 始终做betadisper检验。如果方差不齐，应报告这一情况，可考虑使用PERMANOVA的替代方法或增加样本量。

4. 多种差异分析方法结果不一致时只报告有利的¶

问题： LEfSe、DESeq2、MaAsLin2、ANCOM-BC对同一数据可能给出不同结果。Cherry-picking有利结果是不诚实的。 解决： 报告多种方法的一致性结果（concordance），讨论不一致的原因。以在多种方法中都显著的taxa为高可信结果。

5. 功能预测（PICRUSt2）结果过度解读¶

问题： PICRUSt2基于16S序列预测功能，准确性有限（尤其对新物种或未注释基因），却被当作宏基因组结果解读。 解决： PICRUSt2结果应标注为"predicted functional potential"而非确定的功能。验证应使用真正的宏基因组或代谢组数据。NSTI值高的物种预测不可靠，应过滤。

补充知识¶

常用工具生态¶

分析类型	工具	特点
差异丰度	MaAsLin2	通用多变量回归
差异丰度	ANCOM-BC2	理论严格，处理组成性
差异丰度	ALDEx2	贝叶斯方法
差异丰度	LEfSe	简单快速，适合探索
功能预测	PICRUSt2	16S→功能预测
功能注释	HUMAnN3	宏基因组功能注释
多组学整合	MOFA2	多组学因子分析
网络分析	SpiecEasi	微生物互作网络推断
因果推断	MR	孟德尔随机化

样本量估计¶

检测中等效应的差异菌种通常需要30-50/组
纵向研究可以用更少样本（利用个体内配对信息）
可使用micropower包或模拟法估计样本量

报告规范（STORMS checklist）¶

STORMS (Strengthening The Organizing and Reporting of Microbiome Studies) 是微生物组研究报告规范： - 报告测序平台、引物、数据库版本 - 报告质控过滤标准和保留样本/序列数 - 报告所有检验的多重校正P值 - 提供原始数据到公共数据库（SRA/ENA）