-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1 I have an experiment where samples have 1 of 3 phenotype: H, P and S. All pairwise comparisons need to be made, but the P vs H is the most important. As a result we have the following Agilent 2-colour hybridisations in a loop design (there are 4 biological reps of each): P --------> H H --------> S S --------> P In addition to these, we also have conducted 4 x P ---> H hybridisations as technical rep dye-swaps: P <-------- H Thus we have 16 arrays. I'm trying to figure out how to analyse this in limma, which untill recently I was only familiar with analysing Affy chips. Before I go into code, I'd like to ask: 1) How best to utilise the technical replicate dye-swaps to estimate/remove dye-bias. Can I just use loess normalisation or should I include the dye-effect in the linear model? Anyway, onto some existing code I have. RG <- read.maimages(files=targets$FileName, source="agilent", names=targets$Name) biolrep <- c(1,1,2,3,4,4,5,6,7,7,8,9,10,10,11,12) H <- c(1,-1,-1,0,1,-1,-1,0,1,-1,-1,0,1,-1,-1,0) S <- c(0,0,1,-1,0,0,1,-1,0,0,1,-1,0,0,1,-1) design <- matrix(c(H, S), length(biolrep)) colnames(design) <- c("H", "S") cont.matrix <- cbind( "H vs P" = c(1,0), "S vs P" = c(0,1), "H vs S" = c(1,-1) ) MA.l <- normalizeWithinArrays(RG, method="loess") MA.l.Aq - normalizeBetweenArrays(MA.l, method="Aquantile") corfit <- duplicateCorrelationMA.l.Aq, design, ndups=1, block=biolrep) fit <- lmFitMA.l.Aq, design, block=biolrep, cor=corfit$consensus) fit2 <- contrasts.fit(fit, cont.matrix) fit2 <- eBayes(fit2) result <- decideTests(fit2, adjust.method="BH", p.value=0.05) vennDiagram(result, cex=0.8, main=paste("Genes with BH corrected p-values <= ", 0.05, sep=""), mar=c(0.1,0.2,2,0.2)) If I plot M vs M for a pair of technical replicate dye-swaps all the points, if there are no dye-effects, should centre around the origin 0,0. For my data I can create 4 such plots: # For the raw data ################## MA <- normalizeWithinArrays(RG, method="none") png(file = "Raw_data_dye-bias.png", type = "cairo1", width=10, height=10, units="in", res=300) par(mfrow=c(2,2)) plot(MA$M[,1], MA$M[,2], col=rgb(0,0,0,0.1), pch=20) plot(MA$M[,5], MA$M[,6], col=rgb(0,0,0,0.1), pch=20) plot(MA$M[,9], MA$M[,10], col=rgb(0,0,0,0.1), pch=20) plot(MA$M[,13], MA$M[,14], col=rgb(0,0,0,0.1), pch=20) dev.off() # For the raw data ################## MA <- normalizeWithinArrays(RG, method="loess") png(file = "Loess_normalised_dye-bias.png", type = "cairo1", width=10, height=10, units="in", res=300) par(mfrow=c(2,2)) plot(MA$M[,1], MA$M[,2], col=rgb(0,0,0,0.1), pch=20) plot(MA$M[,5], MA$M[,6], col=rgb(0,0,0,0.1), pch=20) plot(MA$M[,9], MA$M[,10], col=rgb(0,0,0,0.1), pch=20) plot(MA$M[,13], MA$M[,14], col=rgb(0,0,0,0.1), pch=20) dev.off() The images in both these files (see: http://www.bioinf.watsonhaigh.net/pub/tmp/) show a trend (-ve slope) in the data. With all this in mind, I have several questions: 1) Am I correct in thinking the loess normalisation hasn't worked well given that the M vs M plots for tech. rep. dye-swaps show a trend? If there is a large number of affected genes, would this explain the large number of DE genes (~4000 for H vs P) picked up by the limma analysis above? 2) Can I include the dye-effect in the model given that I don't have tech. rep. dye-swaps for all direct comparisons, only H vs P samples? If so, should this work: design <- cbind(Dye=1, design) cont.matrix <- cbind( "H vs P" = c(0,1,0), "S vs P" = c(0,0,1), "H vs S" = c(0,1,-1) ) corfit <- duplicateCorrelationMA.l.Aq, design, ndups=1, block=biolrep) fit <- lmFitMA.l.Aq, design, block=biolrep, cor=corfit$consensus) fit2 <- contrasts.fit(fit, cont.matrix) fit2 <- eBayes(fit2) Cheers, Nathan - -- - -------------------------------------------------------- Dr. Nathan S. Watson-Haigh OCE Post Doctoral Fellow CSIRO Livestock Industries Queensland Bioscience Precinct St Lucia, QLD 4067 Australia Tel: +61 (0)7 3214 2922 Fax: +61 (0)7 3214 2900 Web: http://www.csiro.au/people/Nathan.Watson-Haigh.html - -------------------------------------------------------- -----BEGIN PGP SIGNATURE----- Version: GnuPG v1.4.9 (MingW32) Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org iEYEARECAAYFAkoJ+PAACgkQ9gTv6QYzVL64zgCgnYVJ8GTinzp1koBAk8hB71IL nKwAnAxngPgyHVw1RSTq4LdaqVj696xY =Zul7 -----END PGP SIGNATURE-----