[BioC] limma contrasts with multiple groups
James MacDonald
jmacdon at med.umich.edu
Sun Sep 7 17:21:25 CEST 2008
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
More information about the Bioconductor
mailing list