非平衡嵌套配对设计RNA-seq分析：Limma/DESeq2结果差异咨询

研究背景 #

我有一套RNA-seq数据集，包含3组巨噬细胞样本：成人组（n=6）、足月婴儿组（n=5）、早产婴儿组（n=3），每组样本均包含免疫刺激处理与未处理的配对样本。研究目标是分析组内处理效应，以及组间处理/未处理样本的表达差异，但分别用Limma（结合voom与duplicateCorrelation）和DESeq2（手动构建非满秩模型矩阵）流程分析后，差异表达结果差异极大。

Limma分析代码及结果 #

分析代码 #

#SummarizedExp to DGEList
x <- SE2DGEList(se)
#Introducing the "new" variable
x[["samples"]] <- x[["samples"]] %>% mutate(new_col = paste0(Group, Treatment, SEP = ""))
#Filter dataset 
keep.exprs <- edgeR::filterByExpr(x, group=x[["samples"]][["new_col"]])
x <- x[keep.exprs,, keep.lib.sizes=FALSE]
x <- calcNormFactors(x, method = "TMM")
#Design model matrix 
Treat <- factor(x[["samples"]][["new_col"]])
design <- model.matrix(~0+Treat)
colnames(design) <- levels(Treat)
#Run corfit and voom twice
v <- voom(x, design)
dupcor <- duplicateCorrelation(v,design,block=x[["samples"]][["Index"]])
v <- voom(x, design, block=x[["samples"]][["Index"]], correlation=dupcor$consensus)
corfit <- duplicateCorrelation(v, design, block=x[["samples"]][["Index"]])
#fit/contrasts/eBayes
fit <- lmFit(v,design,block=x[["samples"]][["Index"]],correlation=corfit$consensus)
cm <- makeContrasts(
  Adult_Effect = adultstimulated-adultcontrol,
  Preterm_Effect = pretermstimulated-pretermcontrol,
  Term_Effect = termstimulated-termcontrol,
  cont_AT = termcontrol-adultcontrol,
  cont_AP = pretermcontrol-adultcontrol,
  cont_PT = pretermcontrol-termcontrol,
  trtm_AT = termstimulated-adultstimulated,
  trtm_AP = pretermstimulated-adultstimulated,
  trtm_PT = pretermstimulated-termstimulated,
  levels=design)
fit <- contrasts.fit(fit, cm)
fit <- eBayes(fit)
ab <- decideTests(fit)
summary(ab)

分析结果 #

Adult_Effect Preterm_Effect Term_Effect cont_AT cont_AP cont_PT stim_AT stim_AP stim_PT
Down             34              0           0      38      41       0      43     180       0
NotSig        16990          17061       17063   16998   16970   17064   16968   16638   17064
Up               40              3           1      28      53       0      53     246       0

DESeq2分析代码及结果 #

分析代码 #

#Following the section "matrix not full rank" from the vignette = creating model matrix "on my own"
metadata$Treatment = relevel(metadata$Treatment, "mock")
m1 <- model.matrix(~0 + Group + Group:Index_2 + Group:Treatment, metadata)
all.zero <- apply(m1, 2, function(x) all(x==0))
all.zero
idx <- which(all.zero)
m1 <- m1[,-idx]
#Initialize DESeq2dataset from summarized experiment dataset (se)
dds <- DESeqDataSet(se, design = m1)
#filter
smallestGroupSize <- 6
keep <- rowSums(counts(dds) >= 10) >= smallestGroupSize
dds <- dds[keep,]
#Run DESeq2
dds <- DESeq(dds)
#extract e.g. effects of stimulation on adult cells
x <- results(dds, contrast = list("Groupadult.Treatmentstimulated"))
summary(x)

分析结果 #

out of 16075 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 540, 3.4%
LFC < 0 (down)     : 560, 3.5%
outliers [1]       : 0, 0%
low counts [2]     : 1247, 7.8%
(mean count < 10)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

核心问题 #

1. 针对研究目标，整体分析思路是否存在概念错误？
2. 使用Limma手册9.7节的voom+duplicateCorrelation两次迭代流程是否正确？
3. DESeq2中手动构建非满秩模型矩阵的方式是否合理，代码参数有无误用？

问题解答 #

1. 整体分析思路的概念检查 #

研究目标清晰，针对配对样本设计的分析方向正确：聚焦组内处理配对差异、组间同处理状态样本差异。但两个工具的模型设定逻辑存在明显偏差，这是结果差异悬殊的核心诱因，具体问题看以下工具实现细节。

2. Limma流程的正确性验证 #

你的Limma流程存在两处关键错误：

• 模型矩阵设计冗余矛盾：将Group与Treatment拼接成Treat因子后用~0+Treat建模，同时又通过block=Index控制配对效应，属于重复控制配对信息。正确做法是用包含分组、处理及交互项的模型，结合block控制配对，比如：

  design <- model.matrix(~0 + Group*Treatment, data=x$samples)
  

  再搭配block=x$samples$Index和duplicateCorrelation，才能准确区分分组效应、处理效应、交互效应，同时合理控制配对重复测量的相关性。当前错误设计稀释了处理效应的统计效力，导致组内处理效应差异基因极少。
• voom迭代流程冗余：Limma手册9.7节推荐的流程是：先无block运行voom得到初始权重→用duplicateCorrelation估计相关性→带block和correlation重新运行voom，不需要第二次duplicateCorrelation。你多跑的一次corfit属于冗余操作，但不是结果差异的核心原因。

3. DESeq2非满秩模型的合理性检查 #

你的DESeq2模型构建存在明显逻辑错误：

• 模型公式完全不符合配对设计：~0 + Group + Group:Index_2 + Group:Treatment的公式逻辑混乱，配对样本的正确模型应该是~ Index + Group*Treatment，其中Index为配对个体标识，用于控制个体间差异，同时估计分组、处理及交互效应。当前错误模型未正确控制配对效应，导致假阳性升高，出现大量差异基因。
• 手动构建非满秩矩阵完全没必要：DESeq2会自动处理模型秩的问题，你参考的“matrix not full rank”章节仅针对特殊复杂设计，你的配对设计属于标准情况，无需手动构建模型矩阵并删除零列。
• 过滤标准不一致：Limma用filterByExpr，DESeq2用自定义的计数过滤规则，两者过滤后的基因集不同，也会加剧结果差异，但不是核心原因。

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