生物统计与实验设计（Experimental Design for Bioinformatics）¶

样本量(n) ←→ 效应量(d) ←→ 显著性水平(α) ←→ 功效(1-β)

# 宏基因组示例项目中用的 StratifiedKFold
# 保证每一折中T2D和健康样本的比例与整体一致
from sklearn.model_selection import StratifiedKFold
skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)

        混杂变量（年龄）
         /         \
        ↓           ↓
  暴露因素(T2D) → 结果(菌群变化)
      ↑
   这条关系可能是假的

预算固定 = 总reads数固定

方案A：3个样本 × 100M reads/样本 = 300M reads
方案B：10个样本 × 30M reads/样本 = 300M reads

→ 方案B几乎总是更优（检测到更多差异基因/菌）

错误设计：
  批次1 → 全部健康样本
  批次2 → 全部T2D样本

正确设计：
  批次1 → 5个健康 + 5个T2D
  批次2 → 5个健康 + 5个T2D

#!/usr/bin/env Rscript
# ============================================================
# 样本量计算实操 —— pwr 包（Power Analysis）
# pwr包版本: 1.3-0，基于Cohen (1988) 的方法
# 核心思路：给定效应量、显著性水平、功效，计算需要多少样本
# ============================================================

# 安装（如果没有的话）
# install.packages("pwr")
library(pwr)  # 加载pwr包

# ============================================================
# 场景1：两组t检验的样本量计算
# 问题：比较T2D组和健康组的某个菌属丰度，需要多少样本？
# ============================================================

# 参数设定：
#   d = 0.8   → 大效应量（两组均值差 = 0.8倍标准差，对应约2倍丰度差异）
#   sig.level = 0.05  → 显著性水平5%（犯I类错误的概率）
#   power = 0.8       → 统计功效80%（检测到真实差异的概率）
#   type = "two.sample" → 两组独立样本
#   alternative = "two.sided" → 双侧检验

result_t <- pwr.t.test(
  d = 0.8,             # 效应量（Cohen's d）
  sig.level = 0.05,    # 显著性水平
  power = 0.8,         # 目标功效
  type = "two.sample", # 两组独立样本t检验
  alternative = "two.sided"  # 双侧检验
)

cat("===== 场景1：两组t检验样本量 =====\n")
print(result_t)
cat(sprintf("结论：每组至少需要 %d 个样本\n", ceiling(result_t$n)))
# ceiling()向上取整，因为样本数必须是整数

# ============================================================
# 场景2：不同效应量下的样本量比较
# 实际意义：效应越小，需要越多样本才能检测到
# ============================================================

cat("\n===== 场景2：效应量与样本量的关系 =====\n")
cat("（固定 α=0.05, power=0.8, two-sample, two-sided）\n\n")

effect_sizes <- c(0.2, 0.3, 0.5, 0.8, 1.0, 1.2)  # 从小到大的效应量
for (d in effect_sizes) {
  res <- pwr.t.test(d = d, sig.level = 0.05, power = 0.8,
                    type = "two.sample", alternative = "two.sided")
  cat(sprintf("  d = %.1f → 每组需要 %3d 个样本\n", d, ceiling(res$n)))
}

# ============================================================
# 场景3：已知样本量，计算能达到的功效
# 问题：该 T2D项目只有每组10个样本，功效够不够？
# ============================================================

cat("\n===== 场景3：已有样本量的功效评估 =====\n")
result_power <- pwr.t.test(
  n = 10,              # 每组10个样本
  d = 0.8,             # 假设大效应
  sig.level = 0.05,    # α = 0.05
  type = "two.sample",
  alternative = "two.sided"
)
cat(sprintf("每组10个样本、d=0.8时的功效：%.1f%%\n", result_power$power * 100))
# 如果power < 80%，说明样本量不够

# ============================================================
# 场景4：ANOVA的样本量计算（多组比较）
# 问题：比较3组（健康/轻度T2D/重度T2D）需要多少样本？
# ============================================================

cat("\n===== 场景4：ANOVA(3组)样本量 =====\n")
result_anova <- pwr.anova.test(
  k = 3,               # 3个组
  f = 0.4,             # 中-大效应量（Cohen's f）
  sig.level = 0.05,    # α = 0.05
  power = 0.8          # 目标功效80%
)
cat(sprintf("3组ANOVA(f=0.4)：每组需要 %d 个样本\n", ceiling(result_anova$n)))

# ============================================================
# 场景5：相关性分析的样本量
# 问题：检测菌属丰度与HbA1c的Spearman相关(r=0.3)需要多少样本？
# ============================================================

cat("\n===== 场景5：相关性检验样本量 =====\n")
result_cor <- pwr.r.test(
  r = 0.3,             # 预期相关系数（中等相关）
  sig.level = 0.05,
  power = 0.8,
  alternative = "two.sided"
)
cat(sprintf("检测r=0.3的相关性：需要 %d 个样本\n", ceiling(result_cor$n)))

#!/usr/bin/env Rscript
# ============================================================
# 批次效应校正实操
# 方法1：sva::ComBat() —— 经验贝叶斯方法，最常用
# 方法2：limma::removeBatchEffect() —— 线性模型方法
# ============================================================

# 安装（Bioconductor包）
# if (!require("BiocManager")) install.packages("BiocManager")
# BiocManager::install(c("sva", "limma"))

library(sva)    # ComBat函数所在的包
library(limma)  # removeBatchEffect函数所在的包

# ============================================================
# 1. 模拟有批次效应的表达数据
# ============================================================
set.seed(42)

n_genes <- 1000   # 模拟1000个基因
n_samples <- 20   # 20个样本

# 模拟真实表达矩阵（行=基因，列=样本）
true_expr <- matrix(rnorm(n_genes * n_samples, mean = 8, sd = 2),
                    nrow = n_genes, ncol = n_samples)
rownames(true_expr) <- paste0("Gene", 1:n_genes)
colnames(true_expr) <- paste0("Sample", 1:n_samples)

# 定义分组信息
condition <- factor(c(rep("Healthy", 10), rep("T2D", 10)))

# 定义批次信息（正确设计：每个批次中各组各5个）
batch <- factor(c(rep("Batch1", 5), rep("Batch2", 5),
                  rep("Batch1", 5), rep("Batch2", 5)))

# 添加批次效应（Batch2的所有基因系统性偏高1.5）
batch_effect <- matrix(0, nrow = n_genes, ncol = n_samples)
batch_effect[, batch == "Batch2"] <- rnorm(n_genes * sum(batch == "Batch2"),
                                            mean = 1.5, sd = 0.5)

# 添加真实生物学差异（T2D组前50个基因上调）
bio_effect <- matrix(0, nrow = n_genes, ncol = n_samples)
bio_effect[1:50, condition == "T2D"] <- 2  # 前50个基因在T2D中上调

# 最终观测数据 = 真实表达 + 批次效应 + 生物学差异
observed_expr <- true_expr + batch_effect + bio_effect

cat("===== 数据概况 =====\n")
cat(sprintf("基因数: %d, 样本数: %d\n", nrow(observed_expr), ncol(observed_expr)))
cat(sprintf("批次: %s\n", paste(levels(batch), collapse = ", ")))
cat(sprintf("分组: %s\n", paste(levels(condition), collapse = ", ")))

# ============================================================
# 2. ComBat 校正
# ============================================================
cat("\n===== ComBat 批次效应校正 =====\n")

# 构建模型矩阵（保留我们感兴趣的生物学差异）
mod <- model.matrix(~ condition)  # 保护condition的效应不被去除

# 运行ComBat
# par.prior=TRUE: 使用参数化经验贝叶斯（推荐，尤其是小样本）
combat_expr <- ComBat(
  dat = observed_expr,   # 输入表达矩阵（行=基因，列=样本）
  batch = batch,         # 批次标签
  mod = mod,             # 保护的生物学模型（不去除condition效应）
  par.prior = TRUE       # 参数化先验（适合小样本，TRUE=默认推荐）
)

cat("ComBat校正完成\n")
cat(sprintf("校正前Batch1均值: %.2f, Batch2均值: %.2f\n",
            mean(observed_expr[, batch == "Batch1"]),
            mean(observed_expr[, batch == "Batch2"])))
cat(sprintf("校正后Batch1均值: %.2f, Batch2均值: %.2f\n",
            mean(combat_expr[, batch == "Batch1"]),
            mean(combat_expr[, batch == "Batch2"])))

# ============================================================
# 3. limma::removeBatchEffect 校正
# ============================================================
cat("\n===== limma::removeBatchEffect 校正 =====\n")

# design参数：告诉函数哪些效应要保留
design <- model.matrix(~ condition)

# removeBatchEffect：去除batch效应，保留condition效应
limma_expr <- removeBatchEffect(
  x = observed_expr,     # 输入表达矩阵
  batch = batch,         # 要去除的批次变量
  design = design        # 要保留的设计矩阵
)

cat("limma校正完成\n")
cat(sprintf("校正后Batch1均值: %.2f, Batch2均值: %.2f\n",
            mean(limma_expr[, batch == "Batch1"]),
            mean(limma_expr[, batch == "Batch2"])))

# ============================================================
# 4. PCA可视化验证校正效果
# ============================================================
cat("\n===== PCA验证（看样本是否不再按批次聚类） =====\n")

# 校正前PCA
pca_before <- prcomp(t(observed_expr), scale. = TRUE)
cat("校正前PC1解释方差:", 
    sprintf("%.1f%%", summary(pca_before)$importance[2,1] * 100), "\n")

# 校正后PCA（ComBat）
pca_after <- prcomp(t(combat_expr), scale. = TRUE)
cat("ComBat校正后PC1解释方差:", 
    sprintf("%.1f%%", summary(pca_after)$importance[2,1] * 100), "\n")

# ============================================================
# 5. 差异分析中的正确做法：模型中纳入批次
# ============================================================
cat("\n===== 差异分析中的批次处理（推荐做法） =====\n")
cat("在DESeq2/edgeR/limma中，不要先去批次再分析，\n")
cat("而是在模型公式中直接加入批次：\n")
cat("  DESeq2:  design = ~ batch + condition\n")
cat("  edgeR:   design <- model.matrix(~ batch + condition)\n")
cat("  limma:   design <- model.matrix(~ batch + condition)\n")
cat("\n这样统计模型会同时估计批次效应和处理效应，更稳健。\n")
cat("ComBat/removeBatchEffect的输出主要用于PCA可视化和热图展示。\n")