[BioC] limma separate channel analysis: extremely low consensus intraspotcorrelation

[Thierry Janssens]

Dear Bioconductor,
Dear Dr. Smyth,

I have hybridizid  two color labelled cRNA from five treatments to 10 
custom made 8x15k Agilent oligo arrays, in a saterated/unconnected 
design. I am looking for differentially expressed genes between the five 
treatments. QC reports from the Agilent Feature Extraction Software are 
good (e.g. nice observed/expected R^2 around 0.99). For separate channel 
analysis you need to calculate the consensus correlation of intraspot 
correlation between red and green intensities. This is in my case very 
low (0.06922669). On one of the threads from the past I read that this 
should be 0.8-0.9. What can be the reason for this low correlation 
coefficents?  I am convinced that the function intraSpotCorrelation 
calculates the R and G values from the M value, so I do not have to do 
RG.MA(), right? I can not get rid of the intuition that this function 
calculates the correlation between the M and A values, then my observed 
consensus correlation would make sense.
When I manually set the consensus.correlation to, let's say, 0.9 I do 
get sets of differentially expressed genes (but of course this is not 
the way because I want to know what is going on ...!).

For the interested reader, my script is below here. I would be very 
grateful to get some feedback on this topic. Chapter 9 of the limma 
usersguide, concerning this topic, is rather short.

kind regards,

Thierry Janssens

library(limma)
library(statmod)
limmaUsersGuide(view=TRUE)
setwd("C:/Documents and Settings/thierry/My Documents/Data/Swantje 
Staaden/BioC")
filenames <- readTargets("filelist.txt")
filenames <- as.vector(filenames[, 2])
RGlist <- read.maimages(files = filenames, source = "agilent")
# ...

#remove spike-in  probes
j <- RGlist$genes$ControlType == 0

#normalize within arrays
RGbc <- backgroundCorrect(RGlist, method = "edwards", offset = 30)
MA <- normalizeWithinArrays(RGbc[j, ], method ="loess")

# MA plotjes

pdf("MAplotsnorm.pdf")
for(i in 1:10) {plotMA(MA, array = i)
abline(h=0, col = "red")}
dev.off()

#normalize between arrays
MAbet <- normalizeBetweenArrays(MA, method = "Aquantile")
pdf("densityplots.pdf")
plotDensities(MAbet)
dev.off()

#sort on duplo
o <- order(MA$genes$ProbeUID)
MAsorted <- MA[o,]
o <- order(MAbet$genes$ProbeUID)
MAbetsorted <- MAbet[o,]
r <- 0
for(q in seq(1, nrow(MAbetsorted), 3)) {
   r <- as.numeric((identical(MAbetsorted$genes$probeUID[q], 
MAbetsorted$genes$probeUID[q+1]))
   && (identical(MAbetsorted$genes$probeUID[q], 
MAbetsorted$genes$probeUID[q+2])) ) + r
}
r
# r should be 5069

# Separate channel analysis in limma
MAbetsortedav <- avedups(MAbetsorted, ndups = 3, spacing =1)
targets <- readTargets("filelist.txt")
targetstest <- targetsA2C(targets)
u <- unique(targetstest$Target)
f <- factor(targetstest$Target, levels=u)
design1 <- model.matrix(~0+f)
colnames(design1) <- u
corfit <- intraspotCorrelation(MAbetsortedav, design1)
corfit$consensus.correlation

 > 0.06922669

# ... etc and the script continues ...
 > targetstest
     channel.col Archive      Filename Target
1.1            1    File  SSArray1.txt      B
1.2            2    File  SSArray1.txt      A
2.1            1    File  SSArray2.txt      C
2.2            2    File  SSArray2.txt      B
3.1            1    File  SSArray3.txt     AC
3.2            2    File  SSArray3.txt      C
4.1            1    File  SSArray4.txt     AB
4.2            2    File  SSArray4.txt     AC
5.1            1    File  SSArray5.txt      A
5.2            2    File  SSArray5.txt     AB
6.1            1    File  SSArray6.txt      C
6.2            2    File  SSArray6.txt      A
7.1            1    File  SSArray7.txt     AB
7.2            2    File  SSArray7.txt      C
8.1            1    File  SSArray8.txt      B
8.2            2    File  SSArray8.txt     AB
9.1            1    File  SSArray9.txt     AC
9.2            2    File  SSArray9.txt      B
10.1           1    File SSArray10.txt      A
10.2           2    File SSArray10.txt     AC
 > design1
   B A C AC AB
1  1 0 0  0  0
2  0 1 0  0  0
3  0 0 1  0  0
4  1 0 0  0  0
5  0 0 0  1  0
6  0 0 1  0  0
7  0 0 0  0  1
8  0 0 0  1  0
9  0 1 0  0  0
10 0 0 0  0  1
11 0 0 1  0  0
12 0 1 0  0  0
13 0 0 0  0  1
14 0 0 1  0  0
15 1 0 0  0  0
16 0 0 0  0  1
17 0 0 0  1  0
18 1 0 0  0  0
19 0 1 0  0  0
20 0 0 0  1  0
attr(,"assign")
[1] 1 1 1 1 1
attr(,"contrasts")
attr(,"contrasts")$f
[1] "contr.treatment"

-- 
Thierry K.S. Janssens
Vrije Universiteit Amsterdam
Faculty of Earth and Life Sciences
Institute of Ecological Science
Department of Animal Ecology,
De Boelelaan 1085
1081 HV AMSTERDAM, The Netherlands
Phone: +31 (0)20-5989147
Fax: +31 (0)20-5987123
thierry.janssens at ecology.falw.vu.nl