limma and technical rep dye-swaps
0
0
Entering edit mode
@nathanwatson-haighcsiroau-2863
Last seen 9.6 years ago
-----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-----
affy limma affy limma • 1.1k views

Login before adding your answer.

Traffic: 728 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6