[BioC] limma and technical rep dye-swaps

[Nathan S. Watson-Haigh]

-----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.

<code>
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 <- duplicateCorrelation(MA.l.Aq, design, ndups=1, block=biolrep)
fit <- lmFit(MA.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))
</code>

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:

<code>
# 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()
</code>

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:
<code>
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 <- duplicateCorrelation(MA.l.Aq, design, ndups=1, block=biolrep)
fit <- lmFit(MA.l.Aq, design, block=biolrep, cor=corfit$consensus)
fit2 <- contrasts.fit(fit, cont.matrix)
fit2 <- eBayes(fit2)
</code>

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-----