用limma包进行多个分组的差异分析_多组别limma差异分析

数据集GSE75380
根据描述，是一个DNA芯片的数据
16年的文章，算是挺久了
一共4个组，分别是control, si-sox7, si-sox17, double-knockdown

主要是想梳理一下limma的步骤

1. 载入包

library(idmap1)
library(AnnoProbe)
library(GEOmirror)
library(GEOquery)
library(Biobase)

1. 从GEO数据库获取数据

gset <- geoChina("GSE75380")
gset

1. 从下载的数据里获取基因表达矩阵

eSet <- gset[[1]]
exprSet <- exprs(eSet)

4.基因id转换和注释

eSet@annotation
checkGPL(eSet@annotation)
 
ids <- idmap(eSet@annotation)
dat <- filterEM(exprSet,ids)

1. 获取和添加分组信息

dat <- dat[order(rownames(dat)),]
pd <- pData(eSet)
library(stringr)
group_list=str_split(pd$title,' ',simplify = T)[,1]
table(group_list)

1. 保存以便后续使用

save(dat,group_list,file = 'step1-output.Rdata')

1. 检查矩阵，归一化处理

boxplot(dat,las=2)
dat <- log2(dat+1)

1. 差异分析

library(limma)

# 设定分组
condition <- factor(group_list, levels = c("Control","Sox7","Sox17","Double"),ordered = F)# 这里注意，默认是按字母顺序排列，所以要强行设定一个我自己想要的顺序
condition
table(condition)

# 设定差异比较矩阵 **这里注意了，经常绕不清楚的地方来了**
# 这是需要声明差异比较矩阵的方法
design <- model.matrix(~0+condition)
colnames(design) = levels(factor(condition))
rownames(design) = colnames(dat)
design
 
fit=lmFit(dat,design)
cont.matrix=makeContrasts('Sox7-Control',levels = design)
fit2=contrasts.fit(fit,cont.matrix)
fit2=eBayes(fit2)
options(digits = 4)
a <- topTable(fit2,adjust='BH')

image.png

# 现在是不需要声明差异比较矩阵的方法
design1=model.matrix(~factor(condition))
fit1=lmFit(dat,design1)
fit1=eBayes(fit1)
options(digits = 4)
b <- topTable(fit1,coef=2,adjust='BH')

image.png

实际上是完全一样的！！！！

不信的话可以多试试
法一的Sox7-Control，Sox17-Control，Double-Control分别对应
法二的coef=2，coef=3，coef=4！

作者：小狼小狼_e211
链接：https://www.jianshu.com/p/527134c81a83
来源：简书
著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。