用“limma”进行差异表达分析
1.概述
| 矩阵 | 函数 |
|---|---|
| 表达矩阵 | lmFit |
| 分组矩阵 | eBayes |
| 差异比较矩阵 | 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",las=2)
如果这些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",quote = F)
