R表达矩阵处理例子

BiocManager::install("CLL")
suppressPackageStartupMessages(library(CLL))
data(sCLLex)
exprSet=exprs(sCLLex)   ##sCLLex是依赖于CLL这个package的一个对象
samples=sampleNames(sCLLex)
pdata=pData(sCLLex)
group_list=as.character(pdata[,2]) ##这个也是存储在sCLLex对象里面的信息，描述了样本的状态，需要根据它来进行样本分类
dim(exprSet)
#可以看到我们得到了表达矩阵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可以看到各个芯片直接数据还算整齐，可以进行差异比较
#然后我们制作分组矩阵
library(limma)
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(exprSet)
design
#这个design就是我们制作好的分组矩阵，需要根据我们下载的芯片数据的实验设计方案来，此处例子是CLL疾病的探究，22个样本分成了两组，你们自己的数据只需要按照同样的方法制作即可！ 
#最后一个差异比较矩阵，此处不能省略，即使这里的数据只分成了两组，所以肯定是它们两组直接进行比较啦！如果分成多于两个组，那么就需要制作更复杂的差异比较矩阵。
contrast.matrix<-makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design)
contrast.matrix ##这个矩阵声明，我们要把progres.组跟stable进行差异分析比较


#使用limma这个包来进行差异分析
##step1
fit <- lmFit(exprSet,design)
##step2
fit2 <- contrasts.fit(fit, contrast.matrix) ##这一步很重要，大家可以自行看看效果
fit2 <- eBayes(fit2)  ## default no trend !!!
##eBayes() with trend=TRUE
##step3
tempOutput = topTable(fit2, coef=1, n=Inf)
nrDEG = na.omit(tempOutput) 
#write.csv(nrDEG2,"limma_notrend.results.csv",quote = F)
head(nrDEG)
最后输出的nrDEG就是我们的差异分析结果！ 至于结果该如何解释，大家可以仔细阅读说明书啦！