在R中实现表达矩阵标准化及两组t检验获取p值的方法

一、为什么edgeR::cpm(Counts, log10 = T)会出现异常大数值？ #

• 最常见原因：原始Counts矩阵维度搞反——若样本是行、ENST是列，按行计算cpm时，每个样本的总counts极高，分摊到ENST的数值会异常膨胀。先确认矩阵结构：ENST应为行，样本为列。
• 其次，直接对cpm取log10时，原始counts为0会生成-Inf这类极端值（绝对值大）。edgeR的cpm函数默认log转换不会自动加伪计数，需手动指定参数规避。

二、适合RNA-seq Counts的标准化方法 #

针对两组各3个重复的实验设计，推荐以下两种方案：

1. edgeR TMM标准化（差异分析首选） #

TMM（Trimmed Mean of M-values）是RNA-seq领域校正测序深度偏差的经典方法：

library(edgeR)
# 构建DGEList对象，确保ENST为行、样本为列
dge <- DGEList(counts = Counts)
# 过滤低表达基因（可选，建议保留至少3个样本中cpm>1的基因）
keep <- filterByExpr(dge, group = group)
dge <- dge[keep, , keep.lib.sizes = FALSE]
# 执行TMM标准化
dge <- calcNormFactors(dge, method = "TMM")
# 生成带伪计数的log2(cpm)，避免log(0)
log2_cpm <- cpm(dge, log = TRUE, prior.count = 3)

2. 带伪计数的log10(cpm)（若坚持用log10转换） #

必须添加伪计数避免极端值，两种实现方式：

# 方式1：手动加伪计数后计算cpm再log10
log10_cpm <- log10(cpm(Counts + 1))
# 方式2：用edgeR参数自动处理
log10_cpm <- cpm(Counts, log = TRUE, log.base = 10, prior.count = 1)

三、对每个ENST执行两组t检验并导出p值 #

假设分组信息：group <- factor(c("con1", "con1", "con1", "con2", "con2", "con2"))（对应样本a_25至a_30）

步骤1：准备标准化后的表达矩阵 #

以TMM标准化得到的log2_cpm为例，也可使用上述log10_cpm。

步骤2：批量执行t检验 #

用apply遍历每个ENST（矩阵行）：

# 定义分组索引
con1_idx <- which(group == "con1")
con2_idx <- which(group == "con2")

# 逐行执行t检验并提取p值
p_values <- apply(log2_cpm, 1, function(x) {
  t_test <- t.test(x[con1_idx], x[con2_idx])
  return(t_test$p.value)
})

# 整理为包含ENST编号和p值的数据框
result_df <- data.frame(ENST = rownames(log2_cpm), p_value = p_values, stringsAsFactors = FALSE)

步骤3：导出结果到CSV文件 #

write.csv(result_df, "ENST_t_test_pvalues.csv", row.names = FALSE)

额外建议 #

• 若标准化后数据不符合正态分布，可改用非参数检验，将t.test替换为wilcox.test即可。
• 务必做多重检验校正，计算FDR校正后的p值：

result_df$fdr <- p.adjust(result_df$p_value, method = "fdr")
write.csv(result_df, "ENST_t_test_pvalues_with_fdr.csv", row.names = FALSE)

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