PAYATHON 2026

edgeR 中饮食与提取批次完全混杂，设计矩阵不满秩怎么办？

数据阿明

结论

报错与提取日期的安排有关：所有 ND 样本都在 date_1 提取，所有 WD 样本都在 date_2 提取。diet 与 block_extract 完全混杂，现有数据无法分别估计饮食效应和提取批次效应。

因此，当前模型不能同时将这两个因素作为独立效应，也无法通过 removeBatchEffect 可靠地得到“去除批次后的饮食差异”。可以使用满秩设计继续估计离散度并进行适当的组间比较，但跨饮食的差异仍包含潜在的批次影响。不能因为当前 estimateDisp 报错就直接跳过离散度估计。

为什么设计矩阵不满秩？

给出的 designBlock 中，dietWD 列与提取批次指示列完全相同：

ND 样本：dietWD = 0，批次指示值 = 0
WD 样本：dietWD = 1，批次指示值 = 1

模型中的相关部分可以写成：

dietWD * beta_diet + batch * beta_batch
= dietWD * (beta_diet + beta_batch)

数据只能确定两个系数之和，无法确定它们各自的值。无论同一批次内有多少个重复样本，都不能解决这个不可辨识的问题。

可以检查设计矩阵的秩：

c(
  rank = qr(designBlock)$rank,
  columns = ncol(designBlock)
)

题目展示的矩阵有 5 列，秩为 4。错误中的 not of full rank 表示列之间存在线性依赖。

报错列出某个无法估计的批次系数，不代表只有那个日期有问题。改变参考水平或列顺序，可能会改变报错中列出的系数，但不会消除混杂。

为什么 removeBatchEffect 也会报错？

当前调用通过 designBlock 包含了批次列，又通过 batch = block_extract 要求函数加入批次项，造成了重复建模。

即使从 design 中删除批次列，只保留包含 diet 的生物学设计，函数仍无法同时保留饮食差异并删除与之完全重合的批次差异。

Partial NA coefficients for 16443 probe(s) 表示这些表达特征的拟合中存在无法估计的系数。不能据此认定这 16443 个特征都有批次效应，也不应直接将输出视为有效校正后的结果。

edgeR 差异表达分析应使用计数数据和可辨识的设计矩阵。removeBatchEffect 通常用于 log-expression 的可视化或探索分析；它的输出不应替代 edgeR 模型中的计数输入。

如何继续进行 edgeR 分析？

如果当前重点是分别在每种饮食内部比较 KO 与 CT，可以使用四组设计。它与原先包含 diet、treatment 和交互作用的模型表达相同的四组均值，对比也更直观。

下面假设 diet 和 treatment 的样本顺序与 dgef 计数矩阵的列顺序一致，且 dgef 已完成适用的过滤和归一化步骤：

library(edgeR)

stopifnot(
  length(diet) == ncol(dgef),
  length(treatment) == ncol(dgef)
)

group <- factor(
  paste(diet, treatment, sep = "_"),
  levels = c("ND_CT", "ND_KO", "WD_CT", "WD_KO")
)
stopifnot(!anyNA(group))

design <- model.matrix(~ 0 + group)
colnames(design) <- levels(group)
stopifnot(qr(design)$rank == ncol(design))

y <- estimateDisp(dgef, design, robust = TRUE)
fit <- glmQLFit(y, design, robust = TRUE)

# ND 内：KO - CT
res_ND <- glmQLFTest(
  fit,
  contrast = c(-1, 1, 0, 0)
)

# WD 内：KO - CT
res_WD <- glmQLFTest(
  fit,
  contrast = c(0, 0, -1, 1)
)

# 两种饮食下处理效应的差异：
# (WD_KO - WD_CT) - (ND_KO - ND_CT)
res_interaction <- glmQLFTest(
  fit,
  contrast = c(1, -1, -1, 1)
)

topTags(res_ND)
topTags(res_WD)
topTags(res_interaction)

这些 contrast 的位置对应显式指定的列顺序：ND_CT、ND_KO、WD_CT、WD_KO。如果修改列顺序，必须同步修改 contrast。

矩阵满秩只能确认模型参数在代数上可估计，不能证明批次影响已经消除，也不能证明实验设计足以支持某种因果解释。

哪些结果可以怎样解释？

比较可以支持的解释与限制
ND 内的 KO 与 CT两组都在 date_1 提取。若批次对两组产生相同的加性影响，这部分批次效应不会进入处理对比。
WD 内的 KO 与 CT两组都在 date_2 提取，同样需要批次对两种 treatment 的影响一致这一前提。
WD 与 ND饮食和提取日期同时改变，无法将观察到的差异单独归因于饮食。
两种饮食下 KO−CT 的差异该对比在模型中可估计，但解释为 diet
交互作用，需要假设批次不会改变 treatment 效应。

最后一项需要谨慎解释。如果 date_1 与 date_2 对 KO 和 CT 的影响方式不同，批次与 treatment 的交互作用也可能表现为饮食与 treatment 的交互作用，当前数据无法区分二者。

删除 block_extract 后模型能够拟合，只表示模型不再尝试拆分不可辨识的参数，不表示数据中的批次影响已经被去除。

estimateDisp 失败后能否跳过？

当前错误来自设计矩阵不满秩，应先修正设计矩阵，再重新运行 estimateDisp。

edgeR 的计数模型需要离散度来描述超出 Poisson 变异的生物学变异。某些分析可以使用另有充分依据的离散度，但题目没有提供这种依据，因此不能随意指定一个值，也不能把报错调用中的 dispersion = 0.05 当作后续分析应采用的离散度。

题目中四组都有重复样本，可以使用上面的满秩设计重新估计离散度。估计质量仍取决于实际数据和前处理。

如果目标是独立估计饮食效应，应怎么办？

现有数据缺少区分 diet 与提取日期所需的信息。换参考水平、增加相同结构的重复，或改用其他批次校正方法，都无法补上这部分信息。

后续实验应让同一提取批次包含不同 diet 的样本，并合理安排 treatment，使饮食和提取批次不再完全对应。对于当前数据，应明确报告两者完全混杂，并据此限制跨饮食比较的解释。

备注：内容仅供参考。