[BioC] Re: limma and affy (fwd)

Gordon Smyth smyth at wehi.edu.au
Sun Jun 22 02:23:13 MEST 2003


At 09:50 PM 21/06/2003, Phguardiol at aol.com wrote:
>Gordon,
>thank you very much for these precious information.
>I have few additional questions...
>once you have done this:
>library(affy)
>library(limma) ----> last version as suggested
>data<-ReadAffy()
>data2<-rma(data)
>design <- model.matrix(~ -1+factor(c(1,1,1,2,2,3,3,3)))
>colnames(design) <- c("group1", "group2", "group3")
>fit <- lm.series(exprs(data2), design)
>contrast.matrix <- cbind(group2vs1=c(-1,1,0), group3vs2=c(0,-1,1), 
>group3vs1=c(-1,0,1))
>fit2 <- contrasts.fit(fit, contrast.matrix)
>eb <- ebayes(fit2)
>
> >>> how do you get the list that you can put for example in Excel...? can 
> I save eb in a way that I can open it later in Excel ? surely a stupid 
> question ! (same for fit2)
> >>> then I tried:
>clas <- classifyTests(eb$t, design=design, contrasts=contrast.matrix)
> >>> and I have the same question again ! in addition is there a way to 
> just get the genes with significant variation: 1 / -1 I guess ?
> >>> then I ran:
>options(digits=3)
>toptable(number=30, genelist=gal, fit=fit, eb=eb, adjust="fdr")
> >>> and of course genelist is not recognized, so I have 2 options: either 
> delete the genelist : toptable(number=30, fit=fit, eb=eb, adjust="fdr") 
> but I d like to have the identification ! how can I include the names 
> from the U133A ?

You can't use 'gal' because you don't have a GenePix Allocation List! You 
need instead

   tab <- toptable(coef="group2vs1",genelist=geneNames(data2), fit=fit, 
eb=eb, adjust="fdr")

Then you can write the results to disk using 'write.table', for example

   write.table(tab,file="toptable2vs1.txt",quote=FALSE,row.names=FALSE,sep="\t")

> >>> then finally I ran
>ord <- order(eb$lods, decreasing=TRUE)
>top30 <- ord[1:30]
>M <- fit$coef
>A <- apply(exprs(data2), 1, mean)
>plot(A, M, pch=16, cex=0.1)
> >>> but here it says: Error in xy.coords(x, y, xlabel, ylabel, log) : x 
> and y lengths differ

Well yes, because fit$coef contains three columns corresponding to your 
three contrasts while A is only one column. You need

   M <- fit$coef[,"group2vs1"]

etc

>I think I m going to stop here for this night !
>thanks for your help
>Philippe

Gordon



More information about the Bioconductor mailing list