获取t检验t值用于GSEA分析的R实现方案咨询

一、批量获取每个基因的t值 #

1. For循环实现 #

假设你的数据框expr_data包含分组列group（取值为"normal"/"tumor"）和所有基因的表达列，用循环逐个基因做t检验：

# 提取所有基因名
genes <- colnames(expr_data)[colnames(expr_data) != "group"]
# 创建存储结果的数据框
t_results <- data.frame(gene = genes, t_value = NA, p_value = NA)

# 循环处理每个基因
for(i in seq_along(genes)){
  current_gene <- genes[i]
  # 做两组t检验
  test_res <- t.test(expr_data[[current_gene]] ~ expr_data$group)
  # 提取t值和p值
  t_results$t_value[i] <- test_res$statistic
  t_results$p_value[i] <- test_res$p.value
}

2. tidyverse + broom 简洁批量处理 #

用长格式数据配合分组统计，代码更简洁易读：

library(tidyverse)
library(broom)

# 把宽格式数据转成长格式（每行对应一个样本的一个基因表达值）
long_expr <- expr_data %>%
  pivot_longer(cols = -group, names_to = "gene", values_to = "expression")

# 按基因分组做t检验，提取结果
t_results <- long_expr %>%
  group_by(gene) %>%
  do(tidy(t.test(expression ~ group, data = .))) %>%
  select(gene, statistic, p.value) %>%
  rename(t_value = statistic)

处理完后可以直接整理成fgsea需要的命名向量：

gene_t_vector <- setNames(t_results$t_value, t_results$gene)

3. limma包批量分析（推荐） #

limma是微阵列差异分析的标准工具，结果更稳健，同时能直接输出t值：

library(limma)

# 构建设计矩阵
design <- model.matrix(~ 0 + group, data = expr_data)
colnames(design) <- levels(expr_data$group)

# 拟合线性模型
fit <- lmFit(expr_data[, genes], design)

# 设置对比：肿瘤组 vs 正常组
contrast_mat <- makeContrasts(tumor - normal, levels = design)
fit_contrast <- contrasts.fit(fit, contrast_mat)
fit_contrast <- eBayes(fit_contrast)

# 提取所有基因的结果，包含t值
limma_res <- topTable(fit_contrast, coef = 1, number = Inf)
# 整理为fgsea可用的命名向量
gene_t_limma <- setNames(limma_res$t, rownames(limma_res))

二、关于limma分析后做GSEA是否更优的问题 #

• 若能获取原始未归一化数据：优先用limma流程——先做背景校正、标准化（如分位数归一化），再结合RUV或limma自带的removeBatchEffect处理批次效应，最后做差异分析得到t值用于GSEA。limma的经验贝叶斯校正能提升统计效力，结果更可靠。
• 若只有已RUV归一化的数据：直接用limma对归一化后的数据做差异分析（如上述代码），得到的t值用于GSEA，比普通t检验的结果更稳健。
• 批次效应注意点：如果原始数据可获取，建议在原始数据层面完成批次校正；如果只有归一化后的数据，确认RUV已有效消除批次效应后，再进行后续分析即可。

内容的提问来源于stack exchange，提问作者Rui123