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)