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