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@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