前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >差异分析分组构建到底谁在前面--关于limma包中model.matrix()的问题

差异分析分组构建到底谁在前面--关于limma包中model.matrix()的问题

原创
作者头像
叶子Tenney
发布2023-03-27 17:50:48
3K2
发布2023-03-27 17:50:48
举报
文章被收录于专栏:叶子的数据科技专栏

引言

在使用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语言进行创建:

代码语言:text
复制
# 创建文件名和目标向量
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来创建的。设计矩阵中的第一个系数估计野生型小鼠的平均对数表达,并起到截距的作用,第二个系数估计突变型和野生型之间的差异。

代码语言:txt
复制
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.
代码语言:text
复制
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”

代码语言:text
复制
fit <- lmFit(eset, design)
fit <- eBayes(fit)
topTable(fit, coef = "MUvsWT", adjust = "BH")

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

代码语言:text
复制
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 中的 makeContrastscontrasts.fit 函数实现。可以使用 eBayes 函数和 topTable 函数来识别不同表达的基因,而不需要指定系数。

代码语言:text
复制
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 1
limma User’s Guide 2
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

原创声明:本文系作者授权腾讯云开发者社区发表,未经许可,不得转载。

如有侵权,请联系 cloudcommunity@tencent.com 删除。

原创声明:本文系作者授权腾讯云开发者社区发表,未经许可,不得转载。

如有侵权,请联系 cloudcommunity@tencent.com 删除。

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 引言
  • 过程
  • 结论
  • 引用
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档