如何进行limma对基因芯片数据基因差异表达分析

>suppressPackageStartupMessages(library(CLL))

> data(sCLLex)

> exprSet=exprs(sCLLex) ##sCLLex是依赖于CLL这个package的一个对象

> samples=sampleNames(sCLLex)

> pdata=pData(sCLLex)

> group_list=as.character(pdata[,2])

> dim(exprSet)

[1] 12625 22

> exprSet[1:5,1:5]

CLL11.CEL CLL12.CEL CLL13.CEL CLL14.CEL CLL15.CEL

1000_at 5.743132 6.219412 5.523328 5.340477 5.229904

1001_at 2.285143 2.291229 2.287986 2.295313 2.662170

1002_f_at 3.309294 3.318466 3.354423 3.327130 3.365113

1003_s_at 1.085264 1.117288 1.084010 1.103217 1.074243

1004_at 7.544884 7.671801 7.474025 7.152482 6.902932

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

> suppressMessages(library(limma))

> design <- model.matrix(~0+factor(group_list))

> colnames(design)=levels(factor(group_list))

> rownames(design)=colnames(exprSet)

> design

progres. stable

CLL11.CEL 1 0

CLL12.CEL 0 1

CLL13.CEL 1 0

CLL14.CEL 1 0

CLL15.CEL 1 0

CLL16.CEL 1 0

CLL17.CEL 0 1

CLL18.CEL 0 1

CLL19.CEL 1 0

CLL20.CEL 0 1

CLL21.CEL 1 0

CLL22.CEL 0 1

CLL23.CEL 1 0

CLL24.CEL 0 1

CLL2.CEL 0 1

CLL3.CEL 1 0

CLL4.CEL 1 0

CLL5.CEL 1 0

CLL6.CEL 1 0

CLL7.CEL 1 0

CLL8.CEL 1 0

CLL9.CEL 0 1

attr(,"assign")

[1] 1 1

attr(,"contrasts")

attr(,"contrasts")$`factor(group_list)`

[1] "contr.treatment"

>contrast.matrix<makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design)

> contrast.matrix

Contrasts

Levels progres.-stable

progres. 1

stable -1

> fit <- lmFit(exprSet,design)

> fit2 <- contrasts.fit(fit, contrast.matrix) ##这一步很重要，大家可以自行看看效果

> fit2 <- eBayes(fit2)

> tempOutput = topTable(fit2, coef=1, n=Inf)

> nrDEG = na.omit(tempOutput)

> head(nrDEG)

logFC AveExpr t P.Value adj.P.Val B

39400_at -1.0284628 5.620700 -5.835799 8.340576e-06 0.03344118 3.233915

36131_at 0.9888221 9.954273 5.771526 9.667514e-06 0.03344118 3.116707

33791_at 1.8301554 6.950685 5.736161 1.048765e-05 0.03344118 3.051940

1303_at -1.3835699 4.463438 -5.731733 1.059523e-05 0.03344118 3.043816

36122_at 0.7801404 7.259612 5.141064 4.205709e-05 0.10619415 1.934581

36939_at 2.5471980 6.915045 5.038301 5.362353e-05 0.11283285 1.736846