用“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)

详细链接：https://www.jianshu.com/p/11534a814405