Hi there,

I'm having problems trying to get the limma package to analyse a two-colour dye-swap microarray experiment. The experiment itself is a polysome fractionation on RNA samples both before and after treatment (labelled as C50 below). I want to be able to estimate the change in the polysome ratio between the conditions and assess significance using limma. Each sample is fractionated and in dye-swap pairs.

Essentially the targets table looks like:

Array Cy3      Cy5
1     polyCtrl subCtrl
2     polyC50  subC50
3     subC50   polyC50
4     subCtrl  polyCtrl
... and 4 others (a shuffled biological replicate of the above)

Arrays 1 and 4 are a dye-swap pair, as are 2 and 3 (and the corresponding pairs in the other arrays)

To analyse this, I set up my design like so (following section 11.3 of the limma user guide) as I'm interested in the change in the poly/sub ratio between Ctrl and C50 conditions, taking a possible dye effect into consideration:

> design <- cbind(Dye                  =1,
                polyC50.over.subC50  =c(0,1,-1,0,1,0,-1,0),
                polyCtrl.over.subCtrl=c(1,0,0,-1,0,-1,0,1))

> corfit <-   duplicateCorrelation(MA.filt,design,block=biolrep) # MA.filt is MAList object with extracted, background-corrected, normalised, filtered log2(data)

> fit <- lmFit(MA.filt,design,cor=corfit$consensus.correlation)

and then extract the contrast I'm interested in 

> cont.matrix <- makeContrasts(polyC50.over.subC50-polyCtrl.over.subCtrl,levels=design)

> fit2 <- contrasts.fit (fit, cont.matrix)
> fit2 <- eBayes(fit2)

However, the estimate of the contrast (the log2 ratio of the polysome ratios) is not equal to what I calculate manually using the log2-intensities myself (using rowMeans in the table). Am I analysing incorrectly or are the results from limma correct in some way I'm not getting?

Edit; I probably should add that it's the dye-swap technicality I don't think I'm doing correctly - for a one colour array, it's straightforward to build the right design matrix and extract the contrast in a way that "matches up" with the manual calculation (using rowMeans on the log2-intensity table)

Are you sure you've computed the log2-fold changes correctly? Using rowMeans directly isn't appropriate, because the M-values for the dye-swaps need to be flipped around:

aveM.C50 <- (MA.filt[,2] + MA.filt[,5] - MA.filt[,3] - MA.filt[,7])/4 
aveM.Ctrl <- (MA.filt[,1] + MA.filt[,8] - MA.filt[,4] - MA.filt[,6])/4
logFC <- aveM.C50 - aveM.Ctrl

This should give pretty close values to what topTable reports.