[BioC] limma contrasts with multiple groups

[Tarca, Adi]

Hi all,

I am trying to reproduce the results I get with contrasts.fit function
in limma using simple math, but it does not work when the contrast is
made of more than two groups. Here is the code I am using. Any
suggestions are appreciated!

#a simple matrix with 2 genes, 4 groups and 2 samples per group
eset<-matrix(rnorm(16),2,8)
rownames(eset)<-c("g1","g2")
labs=rep(c("a","b","c","d"),each=2)
TS <- factor(labs)

design <- model.matrix(~0+TS)
colnames(design) <- levels(TS)
fit <- lmFit(eset, design)
cont.matrix <- makeContrasts(
 Treat1=a-b,
 Treat2=(a+b)-(c+d),
 levels=design)
fit2 <- contrasts.fit(fit, cont.matrix)
fit2 <- eBayes(fit2)


#check first contrast for first gene
aT1<-topTable(fit2,coef=1, number=60000, adjust="fdr")
aT1
#which I can get also as:
mean(eset[aT1$ID[1],labs=="a"])-mean(eset[aT1$ID[1],labs=="b"])


#check second contrast for first gene
aT1<-topTable(fit2,coef=2, number=60000, adjust="fdr")
aT1
#which I can NOT get like this 
mean(eset[aT1$ID[1],labs%in%c("a","b")])-mean(eset[aT1$ID[1],labs%in%c("
c","d")])

Thanks,
Adi