197_RNA-seq样本量估计¶

library(RNASeqPower)

# 基本功效分析
# 参数：
# depth: 平均测序深度（每基因的reads数）
# cv: 生物学变异系数（BCV）
# effect: 最小检测的fold change
# alpha: 显著性水平（考虑FDR校正）
# power: 期望的检验力

# 场景1：常规RNA-seq
n1 <- rnapower(
  depth = 20,        # 每基因约20 reads（对应~2000万reads总量）
  cv = 0.4,          # BCV=0.4（人类样本典型值）
  effect = 2,        # 检测2倍变化
  alpha = 0.05/20000, # Bonferroni校正（约20000个基因）
  power = 0.8
)
cat("需要每组样本数:", ceiling(n1), "\n")

# 场景2：不同参数组合
params <- expand.grid(
  depth = c(10, 20, 50, 100),
  cv = c(0.2, 0.4, 0.6),
  effect = c(1.5, 2, 3),
  alpha = 0.05/20000,
  power = 0.8
)

params$n <- apply(params, 1, function(x) {
  ceiling(rnapower(depth = x[1], cv = x[2], effect = x[3],
                    alpha = x[4], power = x[5]))
})

# 查看结果
library(ggplot2)
ggplot(params, aes(x = factor(effect), y = n, fill = factor(cv))) +
  geom_bar(stat = "identity", position = "dodge") +
  facet_wrap(~depth, labeller = label_both) +
  theme_minimal() +
  labs(x = "Fold Change", y = "Samples per Group", fill = "BCV",
       title = "Required Sample Size for RNA-seq")

# install.packages("ssizeRNA")
library(ssizeRNA)

# 基于负二项分布的功效分析
result <- ssizeRNA_single(
  nGenes = 20000,     # 检测的基因数
  pi0 = 0.9,          # 非差异基因比例
  m = 1000,           # 平均reads数
  mu = 10,             # 平均表达量
  disp = 0.1,         # 离散系数
  logfc = log2(2),     # log2 fold change
  fdr = 0.05,          # FDR阈值
  power = 0.8          # 期望检验力
)

cat("推荐每组样本数:", result, "\n")

# 绘制功效曲线
sizes <- 3:20
powers <- sapply(sizes, function(n) {
  # 简化的功效估计
  power_at_n <- 1 - pnorm(qnorm(1 - 0.05/2) - log2(2) / (0.4 / sqrt(n)))
  return(power_at_n)
})

plot(sizes, powers, type = "b", pch = 19,
     xlab = "Samples per Group", ylab = "Power",
     main = "Power Curve for RNA-seq")
abline(h = 0.8, lty = 2, col = "red")
abline(v = sizes[which.min(abs(powers - 0.8))], lty = 2, col = "blue")

library(PROPER)

# 基于DESeq2的模拟
# 使用已有数据估计参数
sim_options <- RNAseq.SimOptions.2grp(
  ngenes = 20000,
  p.DE = 0.1,            # 10%差异基因
  lOD = "300",            # 平均read depth
  lfc = function(x) rnorm(x, mean = 1, sd = 0.5)  # log fold change分布
)

# 运行模拟
sim_result <- runSims(
  Nreps = c(3, 5, 7, 10, 15, 20),  # 测试不同样本量
  sim.opts = sim_options,
  DEmethod = "DESeq2",
  nsims = 50                         # 模拟次数
)

# 提取功效
powers <- comparePower(sim_result, alpha.type = "fdr", alpha.nominal = 0.05,
                        stratify.by = "expr", delta = 0.5)
summaryPower(powers)

# 可视化
plotPower(powers)
plotFDR(powers)

# 快速参考表
sample_size_guide <- data.frame(
  Scenario = c("模式生物，低变异", "人类样本，中等变异", "临床样本，高变异",
                "时间序列实验", "多因素实验", "稀有样本"),
  BCV = c("0.1-0.2", "0.3-0.4", "0.5-0.8", "0.3-0.5", "0.3-0.5", "0.4-0.6"),
  FC_target = c("1.5x", "2x", "2x", "1.5x", "2x", "3x"),
  Min_per_group = c(3, 5, 8, 5, 4, 3),
  Recommended = c(5, 8, 12, 8, 6, 5),
  Note = c("近交系，变异小", "个体差异大", "临床异质性高",
           "需要足够时间点", "每个因素组合", "尽量多收集")
)
print(sample_size_guide)

# RNASeqPower快速计算
library(RNASeqPower)
rnapower(depth=20, cv=0.4, effect=2, alpha=0.05/20000, power=0.8)

# 检验力曲线
sapply(3:20, function(n) rnapower(depth=20, cv=0.4, effect=2, alpha=0.05/20000, n=n))

# PROPER模拟
library(PROPER)
sim <- runSims(Nreps=c(3,5,10), sim.opts=RNAseq.SimOptions.2grp(), DEmethod="DESeq2")