[BioC] limma contrasts with multiple groups

[Tarca, Adi]

Hi Jim,
Thanks for the suggestion. Now I see that actually the call
Treat2=(a+b)-(c+d), in in makeContrasts does:
 mean(a)+mean(b)-[mean(c)+mean(d)]
I thought it would pool together all samples from group a and b and then
compare them with the pooled group c and d.

Thanks again,
Adi


-----Original Message-----
From: James MacDonald [mailto:jmacdon at med.umich.edu] 
Sent: Sunday, September 07, 2008 11:21 AM
To: Tarca, Adi; bioconductor at stat.math.ethz.ch
Subject: Re: [BioC] limma contrasts with multiple groups

Hi Adi,

I think you will see the problem if you look at your contrasts matrix.
In the second column you should have a 0.5 or a -0.5 for the samples you
want to compare, but instead you will have a 1 or -1. You need to change
the makeContrasts call to read like

Treat2=(a+b)/2-(c+d)/2

In order to get what you want - see the limma User's Guide for more
info.

Best,

Jim




James W. MacDonald, M.S.
Biostatistician
Affymetrix and cDNA Microarray Core
University of Michigan Cancer Center
1500 E. Medical Center Drive
7410 CCGC
Ann Arbor MI 48109
734-647-5623
>>> "Tarca, Adi" <atarca at med.wayne.edu> 09/05/08 6:25 PM >>>

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

_______________________________________________
Bioconductor mailing list
Bioconductor at stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/bioconductor
Search the archives:
http://news.gmane.org/gmane.science.biology.informatics.conductor

**********************************************************
Electronic Mail is not secure, may not be read every day, and should not
be used for urgent or sensitive issues