细菌上清液RNA-seq分析中edgeR的estimateCommonDisp报错invalid 'tol' value

实验背景 #

• 细菌上清液RNA-seq分析，包含4组（组1、2、3、4），组2为3次重复，其余组为4次重复，共15个样本。
• 使用edgeR的分析流程如下：

annotation <- data.frame("GeneID" = cds$locus_tag, "chr" = cds$V1 , "start" = cds$V4 , "end" = cds$V5, "strand" = cds$V7)

bams <- dir(path = "./data/input/reads/bams", pattern = "*.BAM", full.names = T, recursive = F)
ann <- annotation
fc_SE <- featureCounts(files = bams, annot.ext = ann, isPairedEnd = T, nthreads = 20)
write.csv(fc_SE$counts, "./data/input/counts.csv", quote = F, row.names = T)
write.csv(fc_SE$stat, "./data/input/stats.csv", quote = F, row.names = T)

counts <- read.csv("./data/input/counts.csv", header = T)
rownames(counts) <- counts[ , 1]

x <- counts[ , 2:16]

group <- factor(c("3", "3", "3", "3", "1", "1", "1", "1", "4", "4", "4", "4", "2", "2", "2"))

y <- DGEList(counts = x, group = group)
design <- model.matrix(~0 + group, data = y$samples)
keep <- filterByExpr(y)
y <- y[keep, keep.lib.sizes = F]
y <- calcNormFactors(y)
y <- estimateCommonDisp(y, design)

报错信息 #

运行至y <- estimateCommonDisp(y, design)时出现错误：

Error in optimize(commonCondLogLikDerDelta, interval = c(1e-04, 100/(100 + :
invalid 'tol' value

临时解决：移除design参数调用estimateCommonDisp(y)可正常运行，但担心影响分析准确性。

解决方案 #

1. 检查过滤后基因数量 #

报错常因过滤后剩余基因过少，导致dispersion估计缺乏足够统计量。运行以下代码查看保留的基因数：

length(keep)

如果基因数少于500，调整filterByExpr的阈值，保留更多有表达的基因：

keep <- filterByExpr(y, min.count=3, min.total.count=10)
y <- y[keep, keep.lib.sizes = F]

2. 替换为GLM-based的dispersion估计函数 #

estimateCommonDisp是旧版的dispersion估计方法，当使用设计矩阵时，edgeR官方更推荐使用estimateGLMCommonDisp，它对复杂设计的兼容性更好，可直接指定tol参数避免报错：

# 替换原estimateCommonDisp行
y <- estimateGLMCommonDisp(y, design, tol=1e-4)

后续补充trended和tagwise dispersion估计，这是edgeR标准流程：

y <- estimateGLMTrendedDisp(y, design)
y <- estimateGLMTagwiseDisp(y, design)

3. 调整设计矩阵形式 #

当前使用的是无截距模型~0 + group，换成带截距的模型尝试：

design <- model.matrix(~group, data=y$samples)
y <- estimateCommonDisp(y, design)

无截距模型在组样本数不均时，可能引发内部优化函数的参数异常。

4. 排查样本质量 #

• 查看各样本总reads数，排除低质量样本：

colSums(y$counts)

如果某样本总reads远低于其他组，考虑移除该样本后重新分析。

• 绘制表达分布箱线图，检查组间整体表达是否存在极端差异：

boxplot(log2(y$counts + 1), las=2)

极端差异可能导致dispersion估计的数值计算问题。

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