Hello,

I'm writing to point at what might be a bug in limma topTable.

I 've been using limma (v2.2 and now v2.4.13) and 
its duplicateCorrelation function to analyse 4 
Nimblegen expression arrays (WT vs M, 2 different 
days of extraction). 4753 genes are represented 
on the arrays, each by twenty probes duplicated 
in two blocks. For each gene I have used the RMA 
function included in the NimbleScan software 
provided by NGS to calculate an expression call 
for each gene in each block.

I found using the dupCor function greatly 
increases for that experiment the number of genes 
found as differentially expressed between WT and 
mutant strains (437 genes have a FDR<0.05 instead 
of 28 when a single value per gene and slide is 
used for the analysis)

Now the problem: I think the topTable contains an 
erroneous A column when one uses the dupcor 
function. It contains the values from fit$Amean 
ie the mean log intensities per gene and per 
duplicate (which length here is 9506) instead of 
the mean log intensities per gene.

I reproduce below the lines showing how I used this duplicateCorrelation

>  souche<-c(2,2,1,1)
>  jour<-c(1,2,1,2)
>  designSJ2<-model.matrix(~factor(jour)+factor(souche))
>  corfitSJ2<-duplicateCorrelation(myEset2, designSJ2,ndups=2,spacing=1)
# Duplicates are grouped in the expression table
Il y a eu 50 avis ou plus (utilisez warnings() pour voir les 50 premiers)
>  corfitSJ2$cor
[1] -0.2367406

#I was initially confused by this negative 
correlation, since the correlation between the 
intensities of the duplicates is high (R2=0.98). 
Now I suppose this negative correlation comes 
from the fact that the between slides correlation 
is also as high, so that the between slide 
variation must be for many genes smaller than the 
intra-slide variation. Can anybody confirm that 
this is a plausible explanation?

>  fitSJ2<-lmFit(myEset2,designSJ2, ndups=2,spacing=1,correlation=-0.2367406)
>  fitSJ2<-eBayes(fitSJ2)
>  topSJ2.souche<-topTable(fitSJ2,coef=3,n=4753, adjust="BH")
>  topSJ2.souche.byid<-topSJ2[order(as.numeric(topSJ2.souche[,1])),]
>  topSJ2.souche.byid[1:10,3]
  [1] 9.478149 9.398650 7.177368 7.119764 7.903540 
7.731239 8.921696 9.065286 7.391708 7.489187
>  fitSJ2$Amean[1:10]
        1        2        3        4        
5        6        7        8        9       10
9.478149 9.398650 7.177368 7.119764 7.903540 
7.731239 8.921696 9.065286 7.391708 7.489187

Best regards,

Christine
--
-----------------------------------
Dr Christine Oger
PRABI
Rhone-Alpes BioInformatics Center
Université Claude Bernard Lyon1
Bâtiment Gregor Mendel
43 blv du 11 nov 1918
F-69622 Villeurbanne
tel 33 (0)4 72 43 29 52
-----------------------------------
	[[alternative HTML version deleted]]

