[BioC] Re: limma and affy (fwd)

[Rafael A. Irizarry]

if you want to use weights coming from standard errors i would write a
summary function (for expresso) based on rlm as opposed to median polish.
it will take much much longer but the resutls will be more meaningful (in
my opinion)


On Sat, 21 Jun 2003, Gordon Smyth wrote:

> 
> >Date: Sat, 21 Jun 2003 00:37:37 EDT
> >From: Phguardiol at aol.com
> >Subject: limma and affy
> >
> >Sir,
> >I m trying to use limma after having extract background and normalize my data
> >(HU133A chips from Affy). then I have 3 different groups with duplicates ou
> >triplicates and I d like to know what are the genes differentially expressed
> >between these 3 groups.
> >here is what I tried to do so:
> >library(affy)
> >library(limma)
> >data <- ReadAffy()
> >data2<- rma(data)
> >design <- c(1,1,1,2,2,3,3,3)
> >fit<-lm.series(exprs(data2), design, weights=1/se.exprs(data2)^2)
> 
> A couple of problems here. First 'design' has to be the actual design 
> matrix, not just a vector of group identifiers. You can calculate the right 
> matrix using
> 
>     design <- model.matrix(~ -1+factor(c(1,1,1,2,2,3,3,3)))
>     colnames(design) <- c("group1", "group2", "group3")
> 
> Secondly, 'rma' does not always calculate standard errors. In this case you 
> should skip the weights and simply use
> 
>    fit <- lm.series(exprs(data2), design)
> 
> To make all pairwise comparisons between your three groups you can estimate 
> the appropriate contrasts by
> 
>    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)
> 
> If you download LIMMA version 1.1.2 then there are some further 
> possibilities, for example
> 
>    clas <- classifyTests(eb$t, design=design, contrasts=contrast.matrix)
> 
> will use F-tests to classify which of the pair-wise differences between the 
> groups are significant for each gene.
> 
> Gordon
> 
> >eb<-ebayes(fit)
> >and then I got this message: Error in if (out$df.prior == Inf) out$s2.post <-
> >rep(out$s2.prior, length(sigma)) else out$s2.post <- (ifelse(df.residual ==
> >:
> >         missing value where TRUE/FALSE needed
> >I wonder if my design object is correctly define. the way it is constructed
> >is that my first 3 chips are representing group 1, the 2 thereafter are
> >belonging to group 2, and the last 3 to group 3. should it be something else ?
> >thanks for your help
> >yours sincerely
> >Philippe Guardiola
> 
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at stat.math.ethz.ch
> https://www.stat.math.ethz.ch/mailman/listinfo/bioconductor
>