⽤“limma”进⾏差异表达分析
1.概述
矩阵表达矩阵分组矩阵
函数lmFiteBayes
差异⽐较矩阵topTable
2.读取表达矩阵:
suppressPackageStartupMessages(library(CLL))data(sCLLex)
exprSet=exprs(sCLLex)
samples=sampleNames(sCLLex)pdata=pData(sCLLex)
group_list=as.character(pdata[,2])dim(exprSet)exprSet[1:5,1:5]group_list
得到表达矩阵exprSet,它的列是各个样本名称,⾏是各个探针ID,⼀个纯粹的表达矩阵,必须是数字型的!可以简单地做⼀下该表达矩阵的QC检测:
par(cex = 0.7)
n.sample=ncol(exprSet)
if(n.sample>40) par(cex = 0.5)cols <- rainbow(n.sample*1.2)
boxplot(exprSet, col = cols,main=\"expression value\
如果这些boxplot发现各个芯⽚直接数据还算整齐,则可以进⾏差异⽐较。
3.制作分组矩阵
suppressMessages(library(limma))
design <- model.matrix(~0+factor(group_list))colnames(design)=levels(factor(group_list))rownames(design)=colnames(exprSet)design
design就是分组矩阵,需要根据我们下载的芯⽚数据的实验设计⽅案来调整。
4.制作差异⽐较矩阵
contrast.matrix<-makeContrasts(paste0(unique(group_list), collapse = \"-\"),levels = design)contrast.matrix
以上已经制作好了必要的输⼊数据,下⾯是如何使⽤limma包来进⾏差异分析!
step 1
fit <- lmFit(exprSet,design)step 2
fit2 <- contrasts.fit(fit, contrast.matrix) fit2 <- eBayes(fit2)step 3
tempOutput = topTable(fit2, coef=1, n=Inf)nrDEG = na.omit(tempOutput) head(nrDEG)
write.csv(nrDEG,\"limma_notrend.results.csv\
详细链接: