191_微生物组常见统计错误¶

# 错误示例：直接对相对丰度做Pearson相关
# 这会产生虚假的负相关！
rel_abundance <- matrix(c(0.3, 0.3, 0.4,
                           0.5, 0.2, 0.3,
                           0.1, 0.4, 0.5), nrow = 3, byrow = TRUE)
cor(rel_abundance)  # 会显示负相关，但这是数学必然而非生物学真实

# 正确方法：使用CLR变换后再计算相关
library(compositions)
clr_data <- clr(rel_abundance)
cor(clr_data)  # 更可靠的相关性估计

# 或使用SparCC（专门处理组成性数据的相关分析）
# sparcc(otu_table)

# 错误做法：对每个物种做t检验，不校正p值
p_values <- apply(otu_table, 2, function(x) {
  t.test(x ~ group)$p.value
})
sum(p_values < 0.05)  # 可能很多"显著"，但多数是假阳性

# 正确做法：BH校正
p_adjusted <- p.adjust(p_values, method = "BH")
sum(p_adjusted < 0.05)  # 校正后的显著结果

# 更好的做法：使用专门的差异丰度工具
library(ANCOMBC)
result <- ancombc2(
  data = phyloseq_obj,
  fix_formula = "group",
  p_adj_method = "BH",
  alpha = 0.05
)
# ANCOM-BC内部处理了组成性偏差和多重检验

# 错误1：直接用原始counts比较
# t.test(raw_counts[group == "A", "OTU1"], raw_counts[group == "B", "OTU1"])
# 结果主要反映测序深度差异

# 错误2：简单相对丰度可能引入组成性偏差
relative <- sweep(otu_table, 1, rowSums(otu_table), "/")

# 推荐方法1：Rarefaction（有争议但仍广泛使用）
library(vegan)
set.seed(42)
otu_rarefied <- rrarefy(otu_table, min(rowSums(otu_table)))

# 推荐方法2：CLR变换
library(compositions)
# 加伪计数处理零值
otu_pseudo <- otu_table + 0.5
clr_transformed <- clr(otu_pseudo)

# 推荐方法3：CSS标准化（metagenomeSeq）
library(metagenomeSeq)
MRobj <- newMRexperiment(t(otu_table))
MRobj <- cumNorm(MRobj, p = cumNormStatFast(MRobj))
css_normalized <- MRcounts(MRobj, norm = TRUE, log = TRUE)

# 推荐方法4：DESeq2的VST/rlog
library(DESeq2)
dds <- DESeqDataSetFromMatrix(t(otu_table), colData = metadata, design = ~ group)
dds <- estimateSizeFactors(dds, type = "poscounts")
vsd <- varianceStabilizingTransformation(dds)

# 零值问题演示
otu_example <- c(0, 0, 0, 0, 0, 5, 8, 12, 0, 0)
# 这个物种只在3个样本中检出

# 错误：直接做log变换（log(0)=-Inf）
# log2(otu_example)  # 产生-Inf

# 方法1：加伪计数（简单但有偏）
log2(otu_example + 1)  # 最常用，但1的选择是任意的

# 方法2：零值替换（如Bayesian multiplicative replacement）
library(zCompositions)
otu_imputed <- cmultRepl(otu_table, method = "CZM", output = "p-counts")

# 方法3：使用零膨胀模型
library(pscl)
# 零膨胀负二项模型
zinb_model <- zeroinfl(abundance ~ group | 1, dist = "negbin", data = df)

# 方法4：使用原生支持零值的差异分析工具
# ANCOM-BC, DESeq2, ALDEx2 都内部处理了零值

# 错误1：对非正态数据使用t检验
# 微生物丰度数据几乎从不是正态分布
shapiro.test(otu_table[, 1])  # 检查正态性
# 应使用非参数检验
wilcox.test(otu_table[, 1] ~ group)  # Wilcoxon检验

# 错误2：对相关丰度数据使用Pearson相关
# 应使用SparCC或基于CLR的相关
# sparcc_result <- sparcc(otu_table, iter = 20, xiter = 10)

# 错误3：对多组比较使用多次t检验
# 应使用Kruskal-Wallis
kruskal.test(abundance ~ group, data = df)
# 事后多重比较
library(dunn.test)
dunn.test(df$abundance, df$group, method = "bh")

# 错误4：PERMANOVA不检查方差齐性
# 必须先做betadisper
library(vegan)
dist_matrix <- vegdist(otu_table, method = "bray")
bd <- betadisper(dist_matrix, metadata$group)
permutest(bd)  # p>0.05 才能信任PERMANOVA结果

# 2024年推荐的差异丰度分析方法

# 1. ANCOM-BC2（推荐首选）
library(ANCOMBC)
result <- ancombc2(
  data = phyloseq_obj,
  fix_formula = "group + age + sex",  # 支持协变量
  p_adj_method = "BH",
  pseudo_sens = TRUE,
  alpha = 0.05
)
significant <- result$res[result$res$q_group < 0.05, ]

# 2. ALDEx2（组成性数据理论基础最严格）
library(ALDEx2)
aldex_result <- aldex(t(otu_table), conditions = group,
                       mc.samples = 128, test = "t", effect = TRUE)
sig_aldex <- aldex_result[aldex_result$we.eBH < 0.05, ]

# 3. MaAsLin2（灵活，支持多种模型）
library(Maaslin2)
maaslin_result <- Maaslin2(
  input_data = otu_table,
  input_metadata = metadata,
  output = "maaslin2_output",
  fixed_effects = c("group"),
  normalization = "TSS",
  transform = "LOG"
)

# 推荐策略：运行2-3种方法，取交集

# FDR校正
p.adjust(p_values, method = "BH")

# CLR变换
compositions::clr(otu_table + 0.5)

# ANCOM-BC快速运行
ancombc2(data = physeq, fix_formula = "group", p_adj_method = "BH")

# ALDEx2快速运行
aldex(counts, conditions, mc.samples = 128, test = "t")

# 正态性检验
shapiro.test(x)  # n<5000