如何进行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