R语言tmod包logFC与MSD计算复现困难及方法问询

一、tmod中logFC与MSD的核心计算逻辑 #

logFC计算逻辑 #

tmod的logFC并非直接计算「处理组均值/对照组均值」的log2值，而是基于limma线性模型拟合后的校正系数：

1. 要求输入的表达矩阵必须是log2转换后的数值（如芯片log2信号值、RNA-seq的logCPM）；
2. 内部自动构建设计矩阵，通过lmFit拟合线性模型，再用contrasts.fit指定处理组vs对照组的对比；
3. 最终的logFC是经验贝叶斯校正后的模型系数，会考虑样本间的变异权重、方差收缩，结果和直接计算组均值的log2差异存在偏差（尤其当样本量小、方差异质性大时）。

MSD计算逻辑 #

tmod的MSD（最小显著差异）是组间差异95%置信区间的半宽，计算步骤：

1. 从limma拟合结果中提取经经验贝叶斯校正后的标准误se = 模型未缩放标准差 * 全局收缩后的sigma；
2. 结合残差自由度，用t分布分位数计算：MSD = qt(0.975, df.residual) * se；
3. MSD的意义是：当组间logFC的绝对值大于MSD时，该基因的表达差异在统计学上显著。

二、复现tmod logFC的正确步骤 #

要复现tmod的logFC，必须完全复刻其limma流程，示例代码：

# 前提：expr为log2转换后的表达矩阵（行=基因，列=样本）
# group为分组向量，如c("ctl", "ctl", "trt", "trt")
library(limma)

# 构建设计矩阵（无截距项，明确分组）
design <- model.matrix(~0 + group)
colnames(design) <- c("ctl", "trt")

# limma拟合流程
fit <- lmFit(expr, design)
contrast <- makeContrasts(trt - ctl, levels=design)
fit2 <- contrasts.fit(fit, contrast)
fit2 <- eBayes(fit2)

# 提取tmod同款logFC和MSD
tmod_logFC <- fit2$coefficients[,1]
se <- fit2$stdev.unscaled[,1] * fit2$sigma
tmod_MSD <- qt(0.975, fit2$df.residual) * se

常见复现失败原因：未做log2转换、设计矩阵/对比方向与tmod不一致、直接计算原始均值的log2比值而非模型系数。

三、替代MSD计算方法 #

1. 传统统计方法（无方差收缩） #

适合样本量充足、方差齐性较好的场景：

• 计算每组均值mean_ctl、mean_trt，标准差sd_ctl、sd_trt；
• 计算标准误se = sqrt(sd_ctl²/n_ctl + sd_trt²/n_trt)；
• 用Welch自由度计算MSD：MSD = qt(0.975, df) * se，其中df = (se²)² / ((sd_ctl²/n_ctl)²/(n_ctl-1) + (sd_trt²/n_trt)²/(n_trt-1))。

2. 基于方差收缩的非limma方法 #

• edgeR：用exactTest计算差异，提取table$logFC和table$lfcSE，MSD = qt(0.975, df) * lfcSE；
• DESeq2：拟合模型后提取results(dds)$lfcSE，MSD = qnorm(0.975) * lfcSE（DESeq2采用正态近似）。

3. 非参数MSD计算 #

针对非正态分布数据，用置换检验计算组间差异的95%置信区间，区间半宽即为非参数MSD：

# 示例：对单个基因计算非参数MSD
gene_expr <- expr["GeneX", ]
ctl_expr <- gene_expr[group == "ctl"]
trt_expr <- gene_expr[group == "trt"]

# 置换检验计算差异分布
perm_diff <- replicate(1000, mean(sample(trt_expr)) - mean(sample(ctl_expr)))
# 提取95%置信区间半宽
ci <- quantile(perm_diff, c(0.025, 0.975))
nonparam_MSD <- abs(ci[2] - (mean(trt_expr)-mean(ctl_expr)))

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