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 )))) }