[BioC] Limma: mean for M values
Ido M. Tamir
tamir at imp.univie.ac.at
Wed Dec 7 17:54:23 CET 2005
Hi,
please apologize my ignorance, but:
I have the impression that limma is calculating the means
of M values by averaging them after adjusting the signs of the
values to the design matrix.
But shouldn't it be something like:
log2( sum(2^M)/ncol(M))
since these values are on a log 2 scale?
A second question is how do I specify contrasts when I would like
to have the genes that are regulated in the same direction in A vs wt
and B vs wt (genes that are down or up in both cases).
Thank you very much in advance for explaining this to me
Sincerely,
Ido Tamir
library("limma")
ap <- 5
am <- -9
bp <- 5
bm <- -7
nu <- 3
rglist <- list()
M <- matrix(c(rep(ap,nu),rep(ap,nu),rep(am,nu),rep(bp,nu),rep(bm,nu)),ncol=5)
A <- matrix(c(rep(12,nu*5)),ncol=5)
names <- c( "A1","A2","A3","B1","B2")
rglist$M <- M
rglist$A <- A
colnames(rglist$M) <- names
colnames(rglist$A) <- names
rglist$genes <- c(paste("g",1:nrow(M),sep=""))
targets <- data.frame( names=names, Cy3=c("wt", "wt", "A", "wt","B"),
Cy5=c("A","A","wt","B","wt"))
design <- modelMatrix(targets, ref="wt")
fit <- lmFit(rglist, design)
cont.matrix <- makeContrasts( A,B, levels=design)
fit2 <- contrasts.fit(fit,cont.matrix)
fit2 <- eBayes(fit2)
coef <- colnames( fit2$coefficients )
print(coef)
for(ci in 1:length(coef)){
tab <- topTable(fit2,coef=coef[ci], number=nrow(rglist$M),
adjust="none",sort.by="M")
print( paste("limma mean", tab$M[1]))
print(paste("A", mean(c( ap, ap, -am ))))
print(paste("B", mean(c( bp, -bm ))))
}
More information about the Bioconductor
mailing list