[BioC] Re: limma and affy (fwd)

[Phguardiol at aol.com]

Hi Gordon,
here I am now : with one potential problem in line 15 that maybe I solved ?

1 library(affy)
2 library(limma)
3 data <- ReadAffy()
4 data2 <- rma(data)
5 design <- model.matrix(~ -1+factor(c(1,1,2,2,3,3)))
6 colnames(design) <- c("group1", "group2", "group3")
7 fit <- lm.series(exprs(data2), design)
8 contrast.matrix <- cbind(group2vs1=c(-1,1,0), group3vs2=c(0,-1,1), 
group3vs1=c(-1,0,1))
9 fit2 <- contrasts.fit(fit, contrast.matrix)
10 eb <- ebayes(fit2)
11 clas <- classifyTests(eb$t, design=design, contrasts=contrast.matrix)
12 qqt(eb$t, df=3+eb$df, pch=16, cex=0.1)
13 abline(0,1)

14 options(digits=3)
#then if I do:
15 tab <- toptable(coef="group2vs1", number=30, genelist=geneNames(data2), 
fit=fit,eb=eb, adjust="fdr")
#I get: Error in toptable(coef = "group2vs1", number = 30, genelist = 
geneNames(data2),  :  subscript out of #bounds
#I wonder if it should not be : 
16 tab <- toptable(coef="group2vs1", number=30, genelist=geneNames(data2), 
fit=fit2,eb=eb, adjust="fdr")
#the more so that it works with fit=fit2 ! If so I guess then after I have to 
use fit2 instead of fit in lines 19, 24 ?
17 ord <- order(eb$lods, decreasing=TRUE)
18 top30 <- ord[1:30]
19 M <- fit2$coef
20 A <- apply(exprs(data2), 1, mean)
21 plot(A, M, pch=16, cex=0.1)
22 text(A[top30], M[top30], cex=0.8, col="blue")
23 abline(0,0,col="blue")
24 plot(fit2$coef[,2], eb$lods[,2], pch=16, cex=0.1, xlab="Log Fold Change", 
ylab="Log Odds", main="group1 vs group2")

25 
write.table(clas,file="clastable2vs1b.txt",quote=TRUE,row.names=TRUE,sep="\t")

there if I open the file clastable2vs1b.txt it is like in affy, first line 
has to be moved one cell to the right. No cure for this "bug" ? I d like to have 
instead of 1,2,3 in the first column, the probe set id like in the tab. how 
can I do that ? I tried: 
   clas <- classifyTests(eb$t, design=design, genelist=geneNames(data2), 
contrasts=contrast.matrix)
   but obtained: Error in classifyTests(eb$t, design = design, genelist = 
geneNames(data2),  : unused argument(s) (genelist ...)
   I tried other options but same result and I ll not detail !
thanks
Philippe

	[[alternative HTML version deleted]]