Hi

I am trying to do some differential expression analysis on 18 Affy pd.hugene.2.0.st arrays. The data is divided in two groups and I have a very quick and simple script to the analysis:

library(oligo)
library(human.db0)
library(hugene20sttranscriptcluster.db)
library(limma)
celFiles <- list.celfiles()
affyRaw <- read.celfiles(celFiles)
eset <- rma(affyRaw, target="core")
design <- cbind(SE=c(1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0),RE=c(0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1))
fit <- lmFit(eset, design)
fit2 <- ebayes(fit)

Checking for class(fit) I get the write type of object, but after running ebayes I get a list on fit2. set object is also fine. Any help appreciated.

sessionInfo is as follows:

R version 3.1.2 (2014-10-31)
Platform: x86_64-apple-darwin14.0.0 (64-bit)

locale:
[1] en_CA.UTF-8/en_CA.UTF-8/en_CA.UTF-8/C/en_CA.UTF-8/en_CA.UTF-8

attached base packages:
[1] stats4    parallel  stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] limma_3.22.1                         pd.hugene.2.0.st_3.10.0              hugene20sttranscriptcluster.db_8.2.0 org.Hs.eg.db_3.0.0                  
 [5] RSQLite_1.0.0                        DBI_0.3.1                            human.db0_3.0.0                      AnnotationDbi_1.28.1                
 [9] GenomeInfoDb_1.2.4                   oligo_1.30.0                         Biostrings_2.34.1                    XVector_0.6.0                       
[13] IRanges_2.0.1                        S4Vectors_0.4.0                      Biobase_2.26.0                       oligoClasses_1.28.0                 
[17] BiocGenerics_0.12.1                 

loaded via a namespace (and not attached):
 [1] affxparser_1.38.0     affyio_1.34.0         BiocInstaller_1.16.1  bit_1.1-12            codetools_0.2-9       ff_2.2-13             foreach_1.4.2        
 [8] GenomicRanges_1.18.4  iterators_1.0.7       preprocessCore_1.28.0 splines_3.1.2         tools_3.1.2           zlibbioc_1.12.0      

You are fitting a 'cell means' model, where you are estimating the mean expression for each group. In this situation you will need the additional step of computing the contrast you are interested in.

contrast <- matrix(c(1,-1), ncol = 1, dimnames = list(c("SE","RE"), "SE - RE"))
fit2 <- contrasts.fit(fit, contrast)
fit2 <- eBayes(fit2)
topTable(fit2,1)

If you care to have annotated results, you need to add that to your fit2 object

gns <- select(hugene20sttranscriptcluster.db, featureNames(eset), c("SYMBOL","GENENAME"))
## there are probesets that measure more than one thing, so you have to remove duplicates
gns <- gns[!duplicated(gns[,1]),]
## this is the most naive way to do so. You might want to do something else - it's up to you
fit2$genes <- gns
topTable(fit2,1)

The reason you got a list instead of an MArrayLM fit object is that you used ebayes() instead of eBayes(). Note the difference in capitalization.

To quote from the eBayes documentation page: "eBayes produces an object of class MArrayLM ... ebayes produces an ordinary list".

As James has pointed out, you have tested an incorrect contrast. James has shown you how to use contrast.fit(). Another even simpler way is to define the design matrix more directly. Then there is no need for contrast.fit:

RE <- factor(c(0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1))
design <- model.matrix(~RE)
fit <- lmFit(eset, design)
fit <- eBayes(fit, trend=TRUE)
topTable(fit, coef=2)

BTW, I get worried when I see people read in files using list.celfiles(), which reads files in alphabetical order. Is it really true that the cel files will sort alphabetically in order of treatment group? Maybe this is fine for you, but file names can be arbitrary and can't always be assumed to sort in a sensible order. I like to create a formal targets frame and connect the treatment conditions to the file names directly.