微生物组差异丰度分析：ALDEx2/ANCOM-BC2/MaAsLin3¶

# ====== R语言实现 ======

# 安装ALDEx2
if (!require("BiocManager")) install.packages("BiocManager")  # 安装BiocManager
BiocManager::install("ALDEx2")  # 安装ALDEx2

library(ALDEx2)  # 加载ALDEx2

# 读取数据
# counts: 行=物种/ASV，列=样本，值=原始计数
counts <- read.csv("otu_table.csv", row.names = 1)  # 读取OTU表
groups <- read.csv("metadata.csv")$Group             # 读取分组信息

# 第一步：CLR变换 + 蒙特卡洛采样
# ALDEx2先从贝叶斯Dirichlet先验采样，再做CLR变换
# mc.samples=128表示采样128次取平均
aldex_clr <- aldex.clr(
    reads = counts,        # 输入原始计数矩阵
    conds = groups,        # 分组条件（如"T2D"和"Control"）
    mc.samples = 128,      # 蒙特卡洛采样次数
    denom = "all"          # 使用所有特征作为CLR的分母
)

# 第二步：统计检验
# 两组比较
aldex_test <- aldex.ttest(
    aldex_clr,             # 输入CLR变换后的数据
    paired.test = FALSE    # 非配对检验
)

# 效应量计算
aldex_effect <- aldex.effect(
    aldex_clr              # 计算效应量
)

# 合并所有结果
aldex_results <- data.frame(aldex_test, aldex_effect)  # 合并结果

# 第三步：筛选显著差异的物种
sig_aldex <- aldex_results[
    aldex_results$wi.eBH < 0.05 &    # BH校正后的Wilcoxon p值<0.05
    abs(aldex_results$effect) > 1,     # 效应量绝对值>1
]

cat("ALDEx2发现的显著差异物种数:", nrow(sig_aldex), "\n")  # 打印数量

# 可视化：效应量图
aldex.plot(aldex_results, type = "MA")  # MA图
aldex.plot(aldex_results, type = "MW")  # MW图（效应量 vs 差异）

# 关键列说明：
# we.ep    — Welch t检验的原始p值
# we.eBH   — Welch t检验的BH校正p值
# wi.ep    — Wilcoxon检验的原始p值
# wi.eBH   — Wilcoxon检验的BH校正p值
# effect   — 效应量（衡量差异大小，>1为大效应）
# overlap  — 两组分布的重叠度（0=完全不重叠，1=完全重叠）

# 推荐使用 wi.eBH（Wilcoxon更稳健）+ effect 联合判断

# 安装ANCOM-BC2
BiocManager::install("ANCOMBC")  # 安装

library(ANCOMBC)    # 加载ANCOMBC
library(phyloseq)   # 加载phyloseq（用于数据管理）

# 创建phyloseq对象
otu <- otu_table(as.matrix(counts), taxa_are_rows = TRUE)  # OTU表
meta <- sample_data(metadata)                               # 元数据
tax <- tax_table(as.matrix(taxonomy))                       # 分类信息

ps <- phyloseq(otu, meta, tax)  # 合并为phyloseq对象

# 运行ANCOM-BC2
ancombc2_result <- ancombc2(
    data = ps,                    # 输入phyloseq对象
    fix_formula = "Group",        # 固定效应公式（分组变量）
    p_adj_method = "BH",          # p值校正方法
    prv_cut = 0.10,               # 最低流行度阈值（至少10%样本中出现）
    lib_cut = 1000,               # 最低文库大小
    group = "Group",              # 分组变量名
    struc_zero = TRUE,            # 检测结构零
    neg_lb = TRUE,                # 使用负下界估计
    pseudo_cnt = TRUE,            # 添加伪计数
    pairwise = TRUE,              # 多组两两比较
    global = TRUE                 # 全局检验
)

# 提取结果
res <- ancombc2_result$res  # 差异丰度结果表

# 筛选显著差异的物种
sig_ancombc <- res[res$q_GroupT2D < 0.05, ]  # BH校正q值<0.05

cat("ANCOM-BC2发现的显著差异物种数:", nrow(sig_ancombc), "\n")

# 查看结构零（某组中完全不存在的物种）
struc_zero <- ancombc2_result$zero_ind  # 结构零指示矩阵

# ANCOM-BC2将零值分为三类：
# 1. 异常零(Outlier zeros) → 当作缺失值处理
# 2. 结构零(Structural zeros) → 一组完全不存在，直接标为差异显著
# 3. 采样零(Sampling zeros) → 用伪计数替换
#
# 还会检验结果对伪计数大小的敏感性

# 安装MaAsLin3
# devtools::install_github("biobakery/maaslin3")  # 从GitHub安装
BiocManager::install("maaslin3")  # 或从Bioconductor安装

library(maaslin3)  # 加载MaAsLin3

# 准备输入数据
# input_data: 行=样本，列=物种/通路，值=丰度
# input_metadata: 行=样本，列=元数据变量

# 运行MaAsLin3
maaslin3_results <- maaslin3(
    input_data = "otu_table.tsv",     # 输入丰度表
    input_metadata = "metadata.tsv",  # 输入元数据
    output = "maaslin3_output",       # 输出目录
    fixed_effects = c("Group", "Age", "BMI"),  # 固定效应
    random_effects = c("SubjectID"),  # 随机效应（纵向数据时使用）
    normalization = "TSS",            # 标准化方法：总和缩放
    transform = "LOG",                # 变换方法：对数变换
    analysis_type = "PREVALENCE,ABUNDANCE",  # 同时分析存在性和丰度
    min_prevalence = 0.1,             # 最低流行度10%
    min_abundance = 0.0001,           # 最低丰度阈值
    correction = "BH",               # 多重比较校正
    cores = 4                         # 并行核心数
)

# MaAsLin3输出：
# 1. all_results.tsv — 完整结果表
# 2. significant_results.tsv — 显著结果
# 3. 丰度关联 + 存在性关联分别报告

# MaAsLin3可以通过多种方式处理组成性：

# 方法1：使用实验方法获得的绝对丰度（如qPCR校正）
# 设置 normalization = "NONE"

# 方法2：使用spike-in标准品校正
# 设置 normalization = "NONE"，输入校正后的数据

# 方法3：计算CLR变换（无实验校正时）
# 设置 normalization = "CLR"

# ====== 推荐策略：同时运行多种方法取交集 ======

# 1. 运行ALDEx2
aldex_sig <- run_aldex2(counts, groups)  # 获取显著物种列表

# 2. 运行ANCOM-BC2
ancombc_sig <- run_ancombc2(ps)  # 获取显著物种列表

# 3. 运行MaAsLin3
maaslin_sig <- run_maaslin3(counts, metadata)  # 获取显著物种列表

# 4. 取交集（至少2个方法同时认为显著的物种）
library(VennDiagram)  # 加载Venn图包

# 画Venn图展示重叠
venn.plot <- venn.diagram(
    x = list(
        ALDEx2 = aldex_sig,
        "ANCOM-BC2" = ancombc_sig,
        MaAsLin3 = maaslin_sig
    ),
    filename = "DA_methods_comparison.png",  # 输出文件
    fill = c("red", "blue", "green"),        # 颜色
    alpha = 0.5,                              # 透明度
    main = "差异丰度分析方法比较"             # 标题
)

# 取交集
consensus <- Reduce(intersect, list(aldex_sig, ancombc_sig, maaslin_sig))
cat("三种方法一致的差异物种:", length(consensus), "\n")

# 工具选择
最保守、最可重复 → ALDEx2
多组比较、重复测量 → ANCOM-BC2
协变量多、需要同时分析存在性 → MaAsLin3
建议策略 → 至少跑2种方法取交集

# 预处理要求
ALDEx2: 输入原始计数，不要做任何标准化
ANCOM-BC2: 输入phyloseq对象
MaAsLin3: 输入丰度表（TSV/CSV）+ 元数据表

# 过滤建议
最低流行度：10-20%（至少10%的样本中出现）
最低丰度：0.01-0.1%
移除单例ASV

# 多重比较校正
BH (Benjamini-Hochberg) — 最常用，控制FDR
Bonferroni — 最严格
q < 0.05 为显著阈值

# R包安装命令
BiocManager::install("ALDEx2")
BiocManager::install("ANCOMBC")
BiocManager::install("maaslin3")