數(shù)據(jù)集GSE75380
根據(jù)描述步悠，是一個DNA芯片的數(shù)據(jù)
16年的文章徙菠，算是挺久了
一共4個組媚值，分別是control, si-sox7, si-sox17, double-knockdown
主要是想梳理一下limma的步驟
- 載入包
library(idmap1)
library(AnnoProbe)
library(GEOmirror)
library(GEOquery)
library(Biobase)
- 從GEO數(shù)據(jù)庫獲取數(shù)據(jù)
gset <- geoChina("GSE75380")
gset
- 從下載的數(shù)據(jù)里獲取基因表達矩陣
eSet <- gset[[1]]
exprSet <- exprs(eSet)
4.基因id轉(zhuǎn)換和注釋
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)
- 保存以便后續(xù)使用
save(dat,group_list,file = 'step1-output.Rdata')
- 檢查矩陣，歸一化處理
boxplot(dat,las=2)
dat <- log2(dat+1)
- 差異分析
library(limma)
# 設(shè)定分組
condition <- factor(group_list, levels = c("Control","Sox7","Sox17","Double"),ordered = F)# 這里注意于置，默認是按字母順序排列穿香，所以要強行設(shè)定一個我自己想要的順序
condition
table(condition)
# 設(shè)定差異比較矩陣 **這里注意了亭引，經(jīng)常繞不清楚的地方來了**
# 這是需要聲明差異比較矩陣的方法
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
# 現(xiàn)在是不需要聲明差異比較矩陣的方法
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宏浩！
