Limma ebayes returns a list with Affy pd.hugene.2.0.st arrays
2
0
Entering edit mode
Paulo Nuin ▴ 200
@paulo-nuin-3012
Last seen 6.9 years ago

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()
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      
ebayes expression affymetrix pd.hugene.2.0.st limma • 1.1k views
3
Entering edit mode
@james-w-macdonald-5106
Last seen 9 days ago
United States

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) ADD COMMENT 1 Entering edit mode You might also want to exclude all the control probesets on this array. Currently, the best way to do this is load(paste0(package.path("pd.hugene.2.0.st"), "/extdata/netaffxTranscript.rda")) annot <- pd(netaffxTranscript) ## check this all.equal(as.character(annot[,1]), featureNames(eset)) ## if not all equal annot <- annot[match(featureNames(eset), as.character(annot[,1]),] ind <- annot$category == "main"

fit2 <- fit2[ind,]

topTable(fit2,1)

0
Entering edit mode

Thanks. I had another procedure to annotate, didn't add to the question. Always good to have more ways to do. Thanks for pointing out the B in eBayes and the other tips.

2
Entering edit mode
@gordon-smyth
Last seen 2 minutes ago
WEHI, Melbourne, Australia

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.

0
Entering edit mode

In this case, the files were provided to me this way and numbered in order, so I ended up with a "lazy loading" procedure. My first time with this type of arrays, and first time with oligo. Thanks for the answer.