微生物组样本量计算¶
一句话概述¶
微生物组研究的统计功效分析（power analysis）——基于HMP等已有数据估计效应量，利用模拟方法计算达到目标统计功效（通常80%）所需的最小样本量，为微生物组实验设计提供定量依据。
核心知识点总览¶
| 知识点 | 关键内容 | 重要程度 |
|---|---|---|
| 功效分析基本概念 | 效应量/alpha/beta/power关系 | ⭐⭐⭐⭐⭐ |
| 效应量估计 | 从pilot/HMP数据估计微生物效应量 | ⭐⭐⭐⭐⭐ |
| 多样性检验的power | Shannon/Beta多样性的样本量 | ⭐⭐⭐⭐ |
| 差异丰度的power | DESeq2/ANCOM样本量需求 | ⭐⭐⭐⭐ |
| 模拟方法 | 参数/非参数模拟框架 | ⭐⭐⭐⭐ |
| 多重检验校正影响 | FDR对power的稀释效应 | ⭐⭐⭐ |
| 纵向研究设计 | 重复测量的power计算 | ⭐⭐⭐ |
| 实用工具 | HMP数据/micropower/powmic | ⭐⭐⭐ |
各步骤详解¶
第一步：样本量计算的基本原理¶
白话解释： 样本量计算回答一个关键问题：要发现你期望的效应，需要至少多少样本？太少了可能发现不了真实差异（假阴性/type II error），太多则浪费资源。功效分析需要三个输入：(1) 你想检测的最小效应量，(2) 可接受的假阳性率（alpha，通常0.05），(3) 期望的检测力（power，通常0.80）。
技术细节： 微生物组样本量计算的特殊挑战： 1. 高维数据——数百个物种同时检验，多重校正大幅降低power 2. 零膨胀——大量零值使方差难以估计 3. 组成性——相对丰度的约束导致非独立性 4. 高变异性——微生物丰度的个体间变异很大（CV>100%） 5. 非正态分布——通常为过离散的负二项分布
# 基本power analysis框架
# 效应量(d) → 样本量(n) → 统计功效(power)
# 关系：power ↑ when n↑, d↑, alpha↑
# 经典两组比较的样本量公式（正态分布假设）
# n = (Z_alpha/2 + Z_beta)^2 * 2 * sigma^2 / delta^2
# 但微生物组数据违反正态假设，需要模拟方法
第二步：从已有数据估计效应量¶
白话解释： 效应量是你期望在实验中发现的"真实差异的大小"。对新研究，可以从三个来源估计：(1) pilot实验数据；(2) 公共数据（如HMP数据库）的变异参数；(3) 同领域已发表文献报告的效应量。
技术细节：
# === 从HMP/已有数据估计微生物组参数 ===
# 1. 使用HMP数据估计方差
# 下载HMP stool样本数据
library(HMP16SData)
hmp_stool <- HMP16SData::V13() # V1-3区域
# 计算各属的均值和方差（用于后续模拟）
genus_table <- otu_table(hmp_stool)
genus_stats <- data.frame(
mean_abundance = rowMeans(genus_table) / colSums(genus_table),
variance = apply(genus_table, 1, var),
prevalence = rowMeans(genus_table > 0),
cv = apply(genus_table, 1, sd) / apply(genus_table, 1, mean)
)
# 2. 从pilot数据估计效应量
# 假设有pilot数据比较两组
pilot_effect <- function(group1_abundance, group2_abundance) {
# log2 fold change
log2fc <- log2(mean(group2_abundance + 0.5) / mean(group1_abundance + 0.5))
# Cohen's d
pooled_sd <- sqrt((var(group1_abundance) + var(group2_abundance)) / 2)
cohens_d <- (mean(group2_abundance) - mean(group1_abundance)) / pooled_sd
return(list(log2fc = log2fc, cohens_d = cohens_d))
}
# 3. 文献报告的典型效应量
# Beta多样性：PERMANOVA R²通常 0.02-0.10
# Alpha多样性：Cohen's d通常 0.3-0.8
# 单菌丰度变化：log2FC通常 0.5-2.0
第三步：Alpha多样性比较的样本量¶
白话解释： 如果你的主要研究问题是"两组的Shannon多样性是否不同？"——这相当于经典的两组t检验（或Wilcoxon），可以用标准的power计算公式。需要估计多样性指数的组间差异和组内变异。
技术细节：
# === Alpha多样性的Power Analysis ===
library(pwr)
# 从已有数据估计参数
# 假设：Shannon diversity
# Control组均值=3.5, SD=0.5
# 期望Treatment组均值=3.0 (降低0.5)
mean_diff <- 0.5
pooled_sd <- 0.5
cohens_d <- mean_diff / pooled_sd # d = 1.0 (大效应)
# 经典t检验的样本量
power_result <- pwr.t.test(d = cohens_d, sig.level = 0.05,
power = 0.80, type = "two.sample")
cat(sprintf("Required n per group (alpha diversity): %d\n", ceiling(power_result$n)))
# 如果效应量较小
for (d in c(0.3, 0.5, 0.8, 1.0, 1.5)) {
n <- ceiling(pwr.t.test(d = d, sig.level = 0.05, power = 0.80,
type = "two.sample")$n)
cat(sprintf("Cohen's d = %.1f → n = %d per group\n", d, n))
}
# d=0.3 → n=176; d=0.5 → n=64; d=0.8 → n=26; d=1.0 → n=17; d=1.5 → n=9
# Wilcoxon检验的样本量（非参数，约为t检验的1.15倍）
n_wilcoxon <- ceiling(power_result$n * 1.15)
第四步：Beta多样性（PERMANOVA）的样本量¶
白话解释： 如果研究问题是"两组的整体菌群组成是否不同？"（PERMANOVA），需要估计组间距离占总距离的比例（R²/效应量）。R²越小需要越多样本。微生物组研究中PERMANOVA的R²通常较小（0.02-0.10）。
技术细节：
# === PERMANOVA Power Analysis ===
# 方法1：基于模拟的power计算
library(vegan)
library(micropower) # 专用包
# 从pilot数据估计PERMANOVA效应量
pilot_dist <- vegdist(t(pilot_abundance), method = "bray")
pilot_result <- adonis2(pilot_dist ~ group, data = pilot_metadata)
observed_R2 <- pilot_result$R2[1]
cat(sprintf("Pilot PERMANOVA R² = %.4f\n", observed_R2))
# 模拟计算power
# 方法：重采样pilot数据，检验不同n下能否检测到显著差异
simulate_permanova_power <- function(abundance_matrix, group_labels,
sample_sizes, n_sim = 100) {
results <- data.frame()
for (n in sample_sizes) {
sig_count <- 0
for (i in 1:n_sim) {
# 从每组随机抽取n个样本
g1_idx <- sample(which(group_labels == levels(group_labels)[1]), n, replace = TRUE)
g2_idx <- sample(which(group_labels == levels(group_labels)[2]), n, replace = TRUE)
sub_data <- abundance_matrix[, c(g1_idx, g2_idx)]
sub_group <- factor(c(rep("G1", n), rep("G2", n)))
sub_dist <- vegdist(t(sub_data), method = "bray")
perm_result <- adonis2(sub_dist ~ sub_group, permutations = 199)
if (perm_result$`Pr(>F)`[1] < 0.05) sig_count <- sig_count + 1
}
results <- rbind(results, data.frame(n = n, power = sig_count / n_sim))
}
return(results)
}
power_curve <- simulate_permanova_power(
pilot_abundance, pilot_metadata$group,
sample_sizes = c(10, 15, 20, 30, 50, 75, 100),
n_sim = 200
)
# 绘制power曲线
plot(power_curve$n, power_curve$power, type = "b",
xlab = "Sample size per group", ylab = "Power",
main = "PERMANOVA Power Curve")
abline(h = 0.8, lty = 2, col = "red")
第五步：差异丰度分析的样本量¶
白话解释： 检测单个微生物的丰度差异（如DESeq2/ANCOM-BC）面临最大的统计挑战：数百个物种同时检验需要多重校正，每个物种的效应量可能不同，且数据是零膨胀的。需要基于模拟的方法估计样本量。
技术细节：
# === 差异丰度Power Analysis ===
# 方法1：使用powmic包（专为微生物组设计）
# devtools::install_github("lichen-lab/powmic")
library(powmic)
# 基于负二项分布模拟
# 输入：pilot数据的参数（均值、离散度）
power_da <- powmic(
n = c(10, 20, 30, 50, 100), # 每组样本量
pi0 = 0.9, # 真阴性比例（即只有10%的OTU是差异的）
foldchange = 2, # 期望的fold change
zeroinfl = 0.3, # 零膨胀比例
disp = 0.5, # 负二项离散参数
depth = 10000, # 测序深度
fdr_level = 0.1 # FDR控制水平
)
# 方法2：手动模拟
simulate_da_power <- function(n_per_group, n_taxa, n_diff, fold_change,
lib_size = 10000, disp = 0.5, n_sim = 100) {
library(DESeq2)
power_results <- numeric(n_sim)
for (sim in 1:n_sim) {
# 模拟count数据（负二项分布）
base_means <- runif(n_taxa, 10, 500)
diff_idx <- sample(n_taxa, n_diff)
counts_g1 <- matrix(rnbinom(n_taxa * n_per_group,
mu = rep(base_means, n_per_group),
size = 1/disp),
nrow = n_taxa)
counts_g2 <- counts_g1
counts_g2[diff_idx, ] <- matrix(rnbinom(n_diff * n_per_group,
mu = rep(base_means[diff_idx] * fold_change, n_per_group),
size = 1/disp),
nrow = n_diff)
# DESeq2分析
count_data <- cbind(counts_g1, counts_g2)
col_data <- data.frame(group = factor(c(rep("A", n_per_group), rep("B", n_per_group))))
dds <- DESeqDataSetFromMatrix(count_data, col_data, ~group)
dds <- DESeq(dds, quiet = TRUE)
res <- results(dds, alpha = 0.1)
# 计算power = 检测到的真阳性 / 真阳性总数
detected <- which(res$padj < 0.1)
tp <- length(intersect(detected, diff_idx))
power_results[sim] <- tp / n_diff
}
return(mean(power_results))
}
# 计算不同样本量下的power
sample_sizes <- c(5, 10, 15, 20, 30, 50)
powers <- sapply(sample_sizes, function(n) {
simulate_da_power(n, n_taxa = 200, n_diff = 20, fold_change = 2)
})
plot(sample_sizes, powers, type = "b",
xlab = "N per group", ylab = "Power (sensitivity)",
main = "Differential Abundance Power")
abline(h = 0.8, lty = 2, col = "red")
第六步：多重检验对Power的影响¶
白话解释： 微生物组分析通常同时检验数百个物种/OTU，FDR校正后p值阈值实际变得很严格。这意味着：即使对单个物种的检测力80%，经FDR校正后可能只有40-60%的差异物种能被检出。样本量估计必须考虑这种"稀释"效应。
技术细节：
# === 多重检验对power的影响 ===
# 假设检测200个OTU，其中20个是真差异的（10%）
n_tests <- 200
n_true_diff <- 20
pi0 <- (n_tests - n_true_diff) / n_tests # 0.9
# BH-FDR下的有效检测阈值
# 如果对单个检验有alpha=0.05的power=0.8，
# 经FDR校正后真正能检出的比例取决于效应量分布
# 模拟FDR对power的影响
n_sim <- 1000
fdr_powers <- numeric(n_sim)
for (i in 1:n_sim) {
# 生成p值
# 真阴性：uniform(0,1)
# 真阳性：beta分布（偏小的p值）
pvals_null <- runif(n_tests - n_true_diff)
pvals_alt <- rbeta(n_true_diff, 0.5, 5) # 偏向0
pvals_all <- c(pvals_null, pvals_alt)
# FDR校正
padj <- p.adjust(pvals_all, method = "BH")
# 检出的真阳性
detected_alt <- sum(padj[(n_tests - n_true_diff + 1):n_tests] < 0.1)
fdr_powers[i] <- detected_alt / n_true_diff
}
cat(sprintf("Average power after FDR (q<0.1): %.2f\n", mean(fdr_powers)))
# 通常比nominal power低很多
第七步：实用样本量推荐与报告¶
白话解释： 综合以上分析，给出最终的样本量推荐。实际中还需考虑：样本损耗率（预计10-20%样本会因质控不合格被排除）、测序深度的影响、以及预算约束。
技术细节：
# === 综合样本量推荐 ===
# 考虑因素汇总
cat("=== 微生物组样本量推荐 ===\n\n")
cat("研究目标 → 推荐每组最小样本量\n")
cat("-------------------------------------------\n")
cat("1. Alpha多样性差异(d>0.8): 20-30 per group\n")
cat("2. Alpha多样性差异(d=0.5): 50-70 per group\n")
cat("3. Beta多样性差异(R²>0.05): 30-50 per group\n")
cat("4. Beta多样性差异(R²=0.02): 100-200 per group\n")
cat("5. 单菌差异丰度(FC>2, top5%): 30-50 per group\n")
cat("6. 单菌差异丰度(FC=1.5): 100+ per group\n")
cat("7. 机器学习分类器: 50-100 per group\n")
cat("8. 纵向时间序列: 20-30 subjects × 5+ timepoints\n\n")
cat("注意事项：\n")
cat("- 以上为最低推荐，建议增加20%余量应对样本损耗\n")
cat("- 测序深度建议 >10,000 reads/sample (16S) 或 >5M reads (shotgun)\n")
cat("- 实际效应量越小，需要样本量越大\n")
cat("- 多中心研究需考虑批次效应，样本量需进一步增加\n")
# 生成power curve图（用于grant申请）
library(ggplot2)
power_df <- data.frame(
n = rep(c(10, 20, 30, 50, 75, 100, 150), 3),
power = c(0.25, 0.45, 0.60, 0.78, 0.88, 0.93, 0.97, # large effect
0.12, 0.25, 0.38, 0.55, 0.70, 0.80, 0.90, # medium effect
0.07, 0.12, 0.18, 0.30, 0.42, 0.55, 0.70), # small effect
effect = rep(c("Large (R²=0.08)", "Medium (R²=0.04)", "Small (R²=0.02)"), each = 7)
)
ggplot(power_df, aes(x = n, y = power, color = effect)) +
geom_line(size = 1.2) + geom_point(size = 2) +
geom_hline(yintercept = 0.8, linetype = "dashed") +
labs(x = "Sample size per group", y = "Statistical Power",
title = "Power Curve for Microbiome Study Design",
color = "Effect Size") +
theme_minimal() +
scale_y_continuous(limits = c(0, 1))
实战命令速查¶
# 快速power计算
library(pwr)
# Alpha多样性
pwr.t.test(d=0.8, sig.level=0.05, power=0.80)
# PERMANOVA模拟
# 见上文 simulate_permanova_power()
# 差异丰度模拟
# 见上文 simulate_da_power()
面试常问点¶
Q1: 微生物组power analysis为什么比传统研究更难？¶
A: (1) 高维多重检验——FDR校正后有效阈值极低；(2) 零膨胀和高离散——方差难以准确估计；(3) 组成性约束——物种间非独立；(4) 效应量通常小（R²常<5%）；(5) 个体间变异极大（CV>100%）。这些因素叠加使得标准公式不适用，必须用模拟方法。
Q2: 如果没有pilot数据，如何估计效应量？¶
A: (1) 使用公共数据（HMP、American Gut等）估计目标人群的基线变异；(2) 参考同领域文献报告的效应量作为先验；(3) 定义"有意义的最小效应量"——如至少2倍丰度变化才有生物学意义；(4) 做多种效应量假设下的sensitivity analysis。保守起见，应假设较小的效应量。
Q3: 测序深度如何影响样本量需求？¶
A: 更深的测序深度降低了技术噪声，提高了低丰度物种的检测率，从而可以用更少的样本检测到差异。但收益递减——超过10,000 reads/sample(16S)后，增加深度对power的提升很小。样本量的增加通常比测序深度的增加更能提升power。因此宁可"shallow sequencing, more samples"。
Q4: 纵向微生物组研究的样本量如何计算？¶
A: 纵向设计（重复测量）比横截面设计power更高——因为每个个体作为自己的对照，减少了个体间变异。通常20-30个受试者×5+时间点就能检测到中等效应。使用混合效应模型的power计算，需要估计：组内相关系数(ICC)、时间点间自相关、以及时间效应大小。
Q5: Grant申请中如何写微生物组样本量论证（sample size justification）？¶
A: 标准格式：(1) 说明主要研究终点（primary endpoint，如Shannon多样性差异或PERMANOVA R²）；(2) 引用pilot或文献数据给出效应量估计及其来源；(3) 声明显著性水平（α=0.05）和目标功效（80%）；(4) 给出计算方法（模拟/公式）和结果；(5) 加入损耗余量（通常+20%）。附上power curve图更有说服力。
易错点¶
1. 使用经典公式忽略微生物组特殊性¶
直接用t检验公式计算微生物丰度比较的样本量严重低估需求（没考虑零膨胀、过离散、FDR等）。
2. 用单个分析目标的power代表整个研究¶
研究通常有多个分析目标（alpha多样性+beta多样性+差异丰度+功能），应以power最低（样本量需求最大）的分析为主。
3. 忽略样本损耗¶
DNA提取失败、测序深度不足、QC不合格等原因通常导致10-30%样本损耗。计算的n应为最终需要的有效样本量，实际招募应增加余量。
4. 对效应量过于乐观¶
首次研究往往高估效应量（publication bias使文献中效应量偏大）。建议按文献报告的一半作为保守估计。
5. 不区分主要终点和次要终点¶
不能对所有分析目标都要求80% power（会导致样本量爆炸）。应明确primary endpoint，次要终点作为探索性分析。
补充知识¶
样本量计算工具¶
- pwr：经典power计算R包
- micropower：微生物组专用power R包
- powmic：差异丰度power模拟
- HMP package：基于HMP数据的power
- G*Power：通用GUI power计算软件
经验法则参考¶
- 16S研究：每组30-50样本是较安全的起点
- 宏基因组研究：每组20-30样本（信息量更大）
- 临床干预研究：视效应量，通常50-100/组
- 发现性研究：可接受较小样本+后续验证
引用推荐¶
- Kelly et al., Bioinformatics, 2015 (micropower, PERMANOVA power estimation)
- La Rosa et al., PLoS ONE, 2012 (HMP-based power, Dirichlet-Multinomial方法)
- Chen L., Bioinformatics, 2020 (powmic)
- 综述: Knight et al., Nature Reviews Microbiology, 2018
- 综述: Ferdous et al., Mucosal Immunology, 2022 ("The rise to power of the microbiome")