引言

在使用limma包进行差异分析的过程中，我们都知道至少需要表达矩阵和分组矩阵两个文件，而在一些例子当中，还出现了一种叫差异比较矩阵的东西，那为什么有些需要有些不需要呢？不需要的会不会得到完全相反的上调下调基因？

其实差异比较矩阵的差距只在于一行代码，是 design <- model.matrix(~Group) 还是 design <- model.matrix(~ 0 + Group) ，那么这个0究竟代表什么含义呢？

过程

根据官方文档 9.2 , 这一段讨论了一个简单的单通道实验，比较了两组老鼠，一组是野生型（Wt），另一组是突变型（Mu）。该实验的目标是识别两组老鼠之间的差异表达基因。为此，提供了两种不同的设计矩阵构建方法。

首先是示例输入矩阵:

FileName Target
File1 WT
File2 WT
File3 Mu
File4 Mu
File5 Mu

可以用R语言进行创建:

1
2
3
4
5
6
7
8
9
10
11
12
## 创建文件名和目标向量
filename <- c("File1", "File2", "File3", "File4", "File5")
target <- c("WT", "WT", "Mu", "Mu", "Mu")
## 使用cbind函数构建矩阵
targets <- data.frame(cbind(filename, target))

##   filename target
## 1    File1     WT
## 2    File2     WT
## 3    File3     Mu
## 4    File4     Mu
## 5    File5     Mu

第一种方法是 实验-对照组比参数化 方法，其中设计矩阵包括突变型和野生型之间差异的系数。设计矩阵是通过为所有样本分配值为1，为突变型组分配值为1，为野生型组分配值为0来创建的。设计矩阵中的第一个系数估计野生型小鼠的平均对数表达，并起到截距的作用，第二个系数估计突变型和野生型之间的差异。

HIere the first coefficient estimates the mean log-expression for wild type mice and plays the role of an intercept. The second coefficient estimates the difference between mutant and wild type.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
Group <- factor(targets$Target, levels=c("WT","Mu"))
design <- model.matrix(~Group)
colnames(design) <- c("WT","MUvsWT")

## r$> design
##   WT MUvsWT
## 1  1      0
## 2  1      0
## 3  1      1
## 4  1      1
## 5  1      1
## attr(,"assign")
## [1] 0 1
## attr(,"contrasts")
## attr(,"contrasts")$Group
## [1] "contr.treatment"

这种方法可以使用 R 中的 lmFit 函数实现。可以使用 eBayes 函数和 topTable 函数来识别不同表达的基因，将系数设置为**“MUvsWT”**。

1
2
3
fit <- lmFit(eset, design)
fit <- eBayes(fit)
topTable(fit, coef = "MUvsWT", adjust = "BH")

第二种方法是组均值参数化方法，其中设计矩阵包括分别为野生型和突变型组分配的系数，并将差异提取为对比。设计矩阵是通过为野生型样本分配值为1，为突变型样本分配值为0，并为突变型样本分配值为1，为野生型样本分配值为0来创建的。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
design <- model.matrix(~0+Group)
colnames(design) <- c("WT", "MU")

## r$> design
##   WT MU
## 1  1  0
## 2  1  0
## 3  0  1
## 4  0  1
## 5  0  1
## attr(,"assign")
## [1] 1 1
## attr(,"contrasts")
## attr(,"contrasts")$Group
## [1] "contr.treatment"

这种方法可以使用 R 中的 makeContrasts 和 contrasts.fit 函数实现。可以使用 eBayes 函数和 topTable 函数来识别不同表达的基因，而不需要指定系数。

1
2
3
4
5
fit <- lmFit(eset, design)
cont.matrix <- makeContrasts(MUvsWT=MU-WT, levels=design)
fit2 <- contrasts.fit(fit, cont.matrix)
fit2 <- eBayes(fit2)
topTable(fit2, adjust="BH")

总之，这一段提供了两种不同的方法，用于创建一个单通道实验的设计矩阵，比较了两组老鼠，野生型和突变型组。这两种方法是处理-对比参数化和组均值参数化方法。这两种方法都可以使用 R 函数实现，可以用于识别两组老鼠之间的不同表达基因。

limma User’s Guide 1
limma User’s Guide 2

结论

因此, 结论是:

仅限两组比较，如已将对照组排在前就可以不要差异比较矩阵，否则将导致结果完全倒转。

原因是 design <- model.matrix(~Group) 会先对需要比较的组进行比较，从第二列开始以对比组填充，而 model.matrix(~ 0 + Group) 只进行分组，不进行比较，如何进行比较由差异比较矩阵和 makeContrasts 函数结果控制。

引用

  1. 关于limma包中model.matrix()的问题-QA-生信技能树
  2. 差异分析是否需要比较矩阵 - basic/makeContrasts.md at master · bioconductor-china/basic · GitHub
  3. limma: Linear Models for Microarray and RNA-Seq Data User’s Guide