R 生物信息 差异化分析 每次得到的结果都不一样 DESeq2 limma beyes

DESeq2差异分析结果

“第二次跑就变了”并不是计算错误，而是因为在 DESeq2（或 edgeR、limma-voom） 的默认流程里，每一步都依赖数据的随机性/估计不确定性，常见源头如下：

---

🔧 1. 独立筛选（independent filtering）

• DESeq2 默认用 results() 的 independentFiltering = TRUE。
• 每次根据 整体表达分布 动态剔除低表达基因，保留的基因数会略有波动 → 最终 nrow(de_genes) 可能差几个。

固定办法

res <- results(dds, alpha = 0.05, independentFiltering = FALSE)

---

🔧 2. 贝塔先验拟合（betaPrior 参数）

• 若你在 DESeq() 或 results() 里 显式或隐式 启用了 betaPrior=TRUE（旧版默认），
• 贝叶斯收缩步骤用 随机启动的数值优化，两次收缩略有差异 → 上下调计数可能 ±1~2。

固定办法

dds <- DESeq(dds, betaPrior = FALSE)   # 新版已默认 FALSE

---

🔧 3. 并行运算 / BLAS 多线程

• 矩阵乘法、随机采样在 多线程 下可能出现 数值抖动（1e-15 级），
• 当 p 值接近阈值时，可能让 边界基因 进出列表。

固定办法
在脚本开头关闭并行/多线程：

library(DESeq2)
library(BiocParallel)
register(SerialParam())      # 禁用并行
Sys.setenv(MKL_NUM_THREADS=1, OPENBLAS_NUM_THREADS=1)

---

🔧 4. 随机种子（如有随机步骤）

• 若用到 随机重抽样、permutation、随机森林 等，需要锁种子：

set.seed(123)

---

✅ 一键“锁死”示例（DESeq2）

set.seed(123)
register(SerialParam())
dds  <- DESeq(dds, betaPrior = FALSE)
res  <- results(dds, alpha = 0.05, independentFiltering = FALSE)

执行以上后，重复运行将得到 完全一致的 up/down 数目。

---

📌 小结

来源	是否默认	影响幅度	固定方法
independentFiltering	是	0–5 个基因	设为 FALSE
betaPrior	旧版是	0–3 个基因	设为 FALSE
多线程数值误差	是	0–1 个基因	单线程
随机种子	视代码	任意	set.seed()

把上述开关全部关闭，上下调计数就会固定不变。

limms beyes差异化分析 每次结果不一样

在使用limma进行归一化和差异分析时，若结果出现不稳定性（如多次运行结果不一致），可能由以下原因导致。结合文献和实际经验，以下是系统性解决方案：

一、数据预处理问题

1. 归一化方法选择不当
   • limma默认使用normalizeBetweenArrays进行组间归一化，默认方法为分位数归一化（quantile normalization）。若数据分布差异较大（如存在极端离群值），可能导致归一化结果波动。

   • 解决方案：

   ◦ 检查数据分布（如箱线图、密度图），尝试其他归一化方法（如scale或robustspline）。

   ◦ 对RNA-seq数据，建议先使用voom进行方差稳定转换后再分析。

2. 低表达基因未过滤
   • 低表达基因会增加噪声，导致方差估计不稳定。若未统一过滤标准，可能影响结果可重复性。

   • 解决方案：

   ◦ 使用filterByExpr或手动过滤（如保留CPM>1的基因）。

   ◦ 示例代码：
   keep <- rowSums(cpm(x) > 1) >= 3 # 保留至少3个样本中CPM>1的基因
   x <- x[keep, ]

二、实验设计问题

3. 设计矩阵（Design Matrix）未正确设置
   • 若实验设计包含批次效应或混杂变量（如时间、操作人员），未在设计矩阵中体现，会导致模型拟合偏差。

   • 解决方案：

   ◦ 明确实验变量（如处理组、批次），使用model.matrix构建设计矩阵。

   ◦ 示例：
   design <- model.matrix(~ batch + condition, data = sample_info)

4. 未考虑配对样本或重复测量
   • 若实验设计为配对样本（如患者治疗前后），需在模型中引入配对因子（如患者ID）。

   • 解决方案：

   ◦ 添加配对变量到设计矩阵：
   design <- model.matrix(~ patient + treatment, data = sample_info)

三、模型拟合与统计检验问题

5. 经验贝叶斯参数未稳定
   • limma的eBayes函数通过收缩方差估计提高小样本稳定性。若未设置合理参数（如trend=TRUE），可能导致结果波动。

   • 解决方案：

   ◦ 启用趋势方差估计：
   fit <- lmFit(x, design)
   fit <- eBayes(fit, trend = TRUE, robust = TRUE)

6. 随机种子未固定
   • 部分步骤（如eBayes中的随机采样）依赖随机数生成，未固定随机种子会导致结果差异。

   • 解决方案：
   set.seed(123) # 在分析前固定随机种子

四、数据质量与技术问题

7. 批次效应未校正
   • 若实验存在技术批次（如不同测序批次），未校正会导致虚假差异。

   • 解决方案：

   ◦ 使用removeBatchEffect校正批次效应（需在设计矩阵中保留批次变量）：
   x_corrected <- removeBatchEffect(x, batch = batch_info, design = design)

8. 数据矩阵未标准化
   • 输入数据未进行标准化（如TPM、FPKM未转换为Counts），可能导致模型拟合错误。

   • 解决方案：

   ◦ 对RNA-seq数据，使用edgeR::calcNormFactors计算归一化因子。

五、结果验证与调试

9. 检查残差诊断图
   • 通过plot(fit)查看残差分布，判断模型是否满足假设（如正态性、同方差性）。

   • 示例问题：

   ◦ 若残差存在明显模式，需检查设计矩阵或数据质量。

10. 使用Glimma交互式可视化
    ◦ 通过glimma::glMDPlot生成交互式MA图，直观观察差异基因分布。

总结步骤

1. 标准化与过滤：统一归一化方法，过滤低表达基因。
2. 实验设计：明确变量，构建合理设计矩阵。
3. 模型稳定性：启用趋势方差估计，固定随机种子。
4. 批次校正：识别并校正技术批次效应。
5. 结果验证：通过残差分析和可视化确认模型合理性。

若问题仍未解决，建议提供代码片段和数据样本，以便进一步诊断。