[BioC] intraSpotCorrelation

[Gordon K Smyth]

Dear Thierry,

I can't comment on whether the low intraspot correlation you observe is a 
feasible value because I don't have any personal experience with the 
Agilent platform you are using.  My earlier comments on the expected size 
of the correlation were written in the context of analysing Agilent Human 
1A oligo arrays, and I don't have good experience with more recent Agilent 
arrays, as they are not used by my collaborators.  There have been some 
suggestions in the literature which might imply that the latest Agilent 
two colour platforms produce data with very small spot correlation 
effects.  If this is true, then it would explain the small correlation you 
see.

My suggestion is that you analyse the data with the estimated correlation, 
even though small.  You definitely should not manually insert an 
artificial value which does not match your data.

Best wishes
Gordon

On Wed, 20 May 2009, Thierry Janssens wrote:

> 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