数据集GSE75380
根据描述，是一个DNA芯片的数据
16年的文章，算是挺久了
一共4个组，分别是control, si-sox7, si-sox17, double-knockdown
主要是想梳理一下limma的步骤
- 载入包
library(idmap1)
library(AnnoProbe)
library(GEOmirror)
library(GEOquery)
library(Biobase)
- 从GEO数据库获取数据
gset <- geoChina("GSE75380")
gset
- 从下载的数据里获取基因表达矩阵
eSet <- gset[[1]]
exprSet <- exprs(eSet)
4.基因id转换和注释
eSet@annotation
checkGPL(eSet@annotation)
ids <- idmap(eSet@annotation)
dat <- filterEM(exprSet,ids)
- 获取和添加分组信息
dat <- dat[order(rownames(dat)),]
pd <- pData(eSet)
library(stringr)
group_list=str_split(pd$title,' ',simplify = T)[,1]
table(group_list)
- 保存以便后续使用
save(dat,group_list,file = 'step1-output.Rdata')
- 检查矩阵，归一化处理
boxplot(dat,las=2)
dat <- log2(dat+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！
