Question: Limma ebayes returns a list with Affy pd.hugene.2.0.st arrays
0
gravatar for Paulo Nuin
4.8 years ago by
Paulo Nuin200
Canada
Paulo Nuin200 wrote:

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      
ADD COMMENTlink modified 4.8 years ago by Gordon Smyth39k • written 4.8 years ago by Paulo Nuin200
Answer: Limma ebayes returns a list with Affy pd.hugene.2.0.st arrays
3
gravatar for James W. MacDonald
4.8 years ago by
United States
James W. MacDonald51k wrote:

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 COMMENTlink written 4.8 years ago by James W. MacDonald51k
1

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)

 

ADD REPLYlink written 4.8 years ago by James W. MacDonald51k

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.

ADD REPLYlink written 4.8 years ago by Paulo Nuin200
Answer: Limma ebayes returns a list with Affy pd.hugene.2.0.st arrays
2
gravatar for Gordon Smyth
4.8 years ago by
Gordon Smyth39k
Walter and Eliza Hall Institute of Medical Research, Melbourne, Australia
Gordon Smyth39k wrote:

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.

ADD COMMENTlink modified 4.8 years ago • written 4.8 years ago by Gordon Smyth39k

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.

ADD REPLYlink written 4.8 years ago by Paulo Nuin200
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 415 users visited in the last hour