610_微生物组样本量估算¶
一句话概述： 微生物组研究的样本量估算是实验设计的关键一步，需要根据效应量、检验功效和统计方法选择合适的样本数，避免"样本不够白做"或"样本太多浪费钱"。
核心知识点速查表¶
| 概念 | 白话解释 |
|---|---|
| 统计功效（Power） | 当真实差异存在时，你能检测到它的概率，通常要求≥80% |
| 效应量（Effect size） | 两组之间差异的大小，越大越容易检测到 |
| I类错误（α） | 假阳性率，通常设为0.05（5%的概率误报） |
| II类错误（β） | 假阴性率，Power = 1 - β |
| PERMANOVA | 基于距离矩阵的非参数多元方差分析 |
| ω² | PERMANOVA中的效应量指标（调整后的R²） |
| Dirichlet-Multinomial | 描述微生物组数据分布的统计模型 |
| Pilot study | 预实验，用来估计参数为正式实验做准备 |
一、为什么微生物组的样本量估算特别难？¶
白话解释：普通临床试验算样本量，比如"比较两组血压差异"，就是两个数字比大小。但微生物组数据有几个"坑"：
- 高维度：几百上千个物种同时存在，不是一个变量而是一个"军团"
- 零膨胀：很多物种在很多样本里检测不到（大量0值）
- 组成性数据：相对丰度加起来等于100%，一个物种增加别的就减少
- 个体差异大：每个人的肠道菌群都不一样，噪音很大
- 多重比较：检验几百个物种，p值要校正，需要更大样本量
二、常用样本量估算方法¶
2.1 基于PERMANOVA的功效分析（micropower包）¶
# === 安装micropower包 ===
# install.packages("devtools") # 如果没安装devtools
devtools::install_github("brendankelly/micropower") # 从GitHub安装
library(micropower) # 加载包
# === 方法1：已知效应量，估算需要的样本量 ===
# omega2：效应量（类似R²，0.01=小,0.06=中,0.14=大）
# alpha：显著性水平
# power：目标功效
result <- power_anova(
omega2 = 0.05, # 效应量（中等偏小）
alpha = 0.05, # 显著性水平
groups = 2, # 组数（比如健康vs疾病）
nsim = 1000 # 模拟次数（越多越准但越慢）
)
print(result) # 打印需要的样本量
# === 方法2：从已有距离矩阵估算 ===
# 如果你有预实验数据的距离矩阵
library(vegan) # 加载vegan包
# 假设已有OTU表和分组信息
otu_table <- read.csv("otu_table.csv", row.names = 1) # 读取OTU表
metadata <- read.csv("metadata.csv") # 读取样本信息
# 计算Bray-Curtis距离
dist_matrix <- vegdist(otu_table, method = "bray") # 计算距离矩阵
# 运行PERMANOVA看效应量
perm_result <- adonis2(
dist_matrix ~ Group, # 按分组检验
data = metadata, # 元数据
permutations = 999 # 置换次数
)
print(perm_result) # 查看R²和p值
# R²就是效应量的粗略估计
# 一般 R² = 0.01-0.05 为小效应
# R² = 0.05-0.15 为中等效应
# R² > 0.15 为大效应
2.2 基于Dirichlet-Multinomial模型的样本量估算¶
# === 安装HMP包（用于DM模型的功效分析）===
install.packages("HMP") # 从CRAN安装
library(HMP) # 加载包
# === 基于DM模型的样本量计算 ===
# 需要两组的丰度向量（从预实验获得）
# pi1：健康组的平均相对丰度
# pi2：疾病组的平均相对丰度
pi1 <- c(0.3, 0.2, 0.15, 0.1, 0.05, 0.2) # 健康组丰度
pi2 <- c(0.1, 0.3, 0.2, 0.15, 0.1, 0.15) # 疾病组丰度
# 计算所需样本量
# 注意：这是简化示例，实际需要从预实验获取参数
result <- MC.Xdc.statistics(
group.data = list(pi1, pi2), # 两组的丰度
n = c(30, 30), # 每组起始样本数
alpha = 0.05 # 显著性水平
)
2.3 使用Evident估算（基于大数据库，Python/QIIME2）¶
# === 安装Evident（QIIME2插件）===
pip install evident # pip安装
# 或者QIIME2插件方式
qiime dev refresh-cache # 刷新QIIME2缓存
# === Python中使用Evident ===
import evident # 导入evident
# Evident可以从大型公开数据库（如American Gut Project）
# 挖掘效应量，用于新研究的功效分析
# 基本用法示例
from evident import UnivariateEffect # 导入单变量效应
# 加载特征表和元数据
effect = UnivariateEffect(
feature_table, # 特征表（OTU/ASV表）
metadata, # 元数据（包含分组信息）
"group_column" # 分组列名
)
# 估算不同样本量下的功效
power_results = effect.power_analysis(
alpha=0.05, # 显著性水平
total_observations=[20, 40, 60, 80, 100] # 测试不同样本量
)
三、不同研究设计的样本量经验值¶
3.1 经验参考表¶
| 研究类型 | 效应量 | 建议最低样本量（每组） | 备注 |
|---|---|---|---|
| 16S扩增子-α多样性 | 大（d>0.8） | 15-20 | Cohen's d |
| 16S扩增子-β多样性 | 中（ω²~0.05） | 30-50 | PERMANOVA |
| 16S扩增子-差异物种 | 小-中 | 50-100 | 多重比较校正后 |
| 宏基因组-功能差异 | 中 | 40-80 | 更多特征需更多样本 |
| 纵向研究 | 中 | 20-30人×多时间点 | 考虑丢失率 |
| 多中心研究 | 小 | 100+ | 批次效应需更多样本 |
3.2 常见疾病研究的效应量参考¶
# === 从文献中获取效应量的经验值 ===
# 以下数据来自大规模微生物组研究的荟萃分析
effect_sizes <- data.frame(
disease = c("IBD", "CRC", "T2D", "Obesity", "Depression"), # 疾病
alpha_div_effect = c(0.8, 0.4, 0.3, 0.2, 0.15), # α多样性效应量
beta_div_R2 = c(0.10, 0.03, 0.02, 0.01, 0.01), # β多样性R²
min_samples = c(30, 50, 80, 100, 150) # 建议最低样本量/组
)
print(effect_sizes) # 打印效应量参考表
# 注意：发表的效应量往往被高估（发表偏倚）
# 建议用发表效应量的50%来估算样本量
四、实战：完整的样本量估算流程¶
#!/usr/bin/env Rscript
# 微生物组研究样本量估算完整脚本
# === 第1步：加载必要的包 ===
library(vegan) # 群落生态学分析
library(pwr) # 基础功效分析
# === 第2步：基于预实验数据估算效应量 ===
# 读取预实验OTU表
otu <- read.csv("pilot_otu_table.csv", row.names = 1) # OTU表
meta <- read.csv("pilot_metadata.csv") # 元数据
# 计算α多样性
shannon <- diversity(otu, index = "shannon") # Shannon指数
# 估算α多样性的效应量（Cohen's d）
group1 <- shannon[meta$Group == "Control"] # 对照组Shannon
group2 <- shannon[meta$Group == "Disease"] # 疾病组Shannon
cohens_d <- abs(mean(group1) - mean(group2)) /
sqrt((sd(group1)^2 + sd(group2)^2) / 2) # 计算Cohen's d
cat("Cohen's d for Shannon diversity:", cohens_d, "\n")
# === 第3步：计算α多样性所需样本量 ===
pwr_result <- pwr.t.test(
d = cohens_d, # 效应量
sig.level = 0.05, # 显著性水平
power = 0.8, # 目标功效80%
type = "two.sample" # 双样本t检验
)
cat("每组需要样本数（α多样性）:", ceiling(pwr_result$n), "\n")
# === 第4步：估算β多样性所需样本量 ===
dist_bc <- vegdist(otu, method = "bray") # Bray-Curtis距离
perm <- adonis2(dist_bc ~ Group, data = meta, permutations = 999)
r2 <- perm$R2[1] # 提取R²
cat("PERMANOVA R²:", r2, "\n")
# 经验公式：β多样性的样本量通常是α多样性的1.5-2倍
n_beta <- ceiling(pwr_result$n * 1.5) # 粗略估算
cat("每组建议样本数（β多样性）:", n_beta, "\n")
# === 第5步：考虑多重检验校正 ===
n_taxa <- ncol(otu) # 物种数量
# Bonferroni校正后需要更大样本量
# 经验：物种数每增加10倍，样本量增加约30%
correction_factor <- 1 + 0.3 * log10(n_taxa) # 校正系数
n_corrected <- ceiling(pwr_result$n * correction_factor) # 校正后样本量
cat("校正多重比较后每组样本数:", n_corrected, "\n")
# === 第6步：汇总建议 ===
cat("\n=== 样本量建议汇总 ===\n")
cat("α多样性比较：每组", ceiling(pwr_result$n), "个样本\n")
cat("β多样性比较：每组", n_beta, "个样本\n")
cat("差异物种分析：每组", n_corrected, "个样本\n")
cat("最终建议（取最大值）：每组",
max(ceiling(pwr_result$n), n_beta, n_corrected), "个样本\n")
cat("建议再加20%冗余（考虑样本丢失）\n")
五、功效分析可视化¶
# === 画功效曲线图 ===
library(ggplot2) # 加载ggplot2
# 不同样本量下的功效
sample_sizes <- seq(10, 200, by = 10) # 样本量范围
powers <- sapply(sample_sizes, function(n) {
pwr.t.test(n = n, d = 0.5, sig.level = 0.05)$power # 计算功效
})
# 创建数据框
power_df <- data.frame(
n = sample_sizes, # 样本量
power = powers # 功效
)
# 画图
ggplot(power_df, aes(x = n, y = power)) + # 映射x和y
geom_line(color = "blue", linewidth = 1.2) + # 画线
geom_hline(yintercept = 0.8, linetype = "dashed", color = "red") + # 80%功效线
annotate("text", x = 150, y = 0.82, label = "Power = 80%", color = "red") +
labs(
title = "功效分析曲线", # 标题
x = "每组样本量", # x轴标签
y = "统计功效" # y轴标签
) +
theme_minimal() # 简洁主题
ggsave("power_curve.pdf", width = 8, height = 6) # 保存图片
六、常见报错与解决¶
| 报错/问题 | 原因 | 解决方案 |
|---|---|---|
| 效应量太小，算出需要>1000样本 | 组间差异确实很小 | 考虑是否值得做这个研究 |
| 预实验样本量不够估计效应量 | pilot study太小 | 至少每组10个样本做预实验 |
| PERMANOVA R²波动很大 | 样本量不足导致估计不稳定 | 多次置换取平均，增加预实验样本 |
| 发表的效应量用来算样本量后太少 | 发表偏倚导致效应量被高估 | 用发表效应量的50%来计算 |
| 不同指标算出的样本量差异大 | 不同分析目标要求不同 | 取最大值作为最终样本量 |
七、面试高频题¶
Q1：微生物组研究为什么要做样本量估算？¶
答： 样本量不够会导致统计功效不足，真实的差异检测不到（假阴性），研究白做。样本量过多则浪费资源。微生物组数据因为高维度、零膨胀和个体差异大，比普通临床研究需要更多样本。一般建议统计功效达80%以上。
Q2：什么是PERMANOVA，为什么用它来做功效分析？¶
答： PERMANOVA是基于距离矩阵的非参数方差分析，特别适合微生物组数据。它不要求数据正态分布，能处理高维数据。在微生物组中常用Bray-Curtis或UniFrac距离来衡量群落差异，PERMANOVA检验这些距离在不同组间是否显著不同。用PERMANOVA做功效分析就是模拟不同样本量下能否检测到组间距离差异。
Q3：效应量小的微生物组研究怎么设计？¶
答： 如果效应量很小（如肥胖与菌群的关系，R²约0.01），可以：(1) 增加样本量，至少每组100个以上；(2) 控制混杂变量（饮食、药物等）减少噪音；(3) 使用配对设计或纵向设计提高检验效率；(4) 聚焦功能通路而非单个物种，可能效应量更大；(5) 多中心合作增加样本量。
Q4：你在毕业设计中是如何考虑样本量的？¶
答： 在2型糖尿病肠道菌群项目中，我使用了公开数据集。在分析前，通过预实验数据的PERMANOVA计算了效应量（R²），确认样本量是否足够。同时参考了已发表T2D微生物组研究的样本量（通常每组50-100），确保分析有足够的统计功效。对于机器学习模型，还需要额外考虑训练集/测试集划分的样本量需求。
参考资料：Kelly et al., Bioinformatics 2015（micropower）| La Rosa et al., PLOS ONE 2012（HMP）| Mattiello et al., Bioinformatics 2016 | Evident: Rahman et al., Genes 2023 | Nature Mucosal Immunology 2022 综述