[BioC] dye swaps of technical replicates and variable numbers of
replicate spots
Ramon Diaz-Uriarte
rdiaz at cnio.es
Thu Aug 21 17:59:18 MEST 2003
> I haven't tried out any technical replicate analyses myself yet and I
> haven't worked out how the details will work. It would appear that one
> can't do it with a fixed effect analysis at the log-ratio level. So you're
> on your own ... good luck!
Thanks anyway. I'll go for it, and let you know.
However, in case it might interest anyone, I think we might be able to use the
fixed effects approach too. This is what I did:
1. First, forget initially about the dye-swaps.
2. The estimate we want is the mean of the biological replicates. The rest of
the design matrix is just to account for the technical replicates part.
3. Set up the design matrix using, for instance, Helmert contrasts (with a
first column of 1s). This way, the term for the intercept is the average of
the biological replicates.
4. If we now use lm.series we get the estimate we are interested in in the
first coeff. and, I think, the right s.e.
5. If we now want to incorporate dye-swap again we just multiply the above
design matrix by the appropriate vector of alternating 1s and -1s.
6. If we use lme we get, of course, slightly different answers (which get
closer the smaller the inter replicate variation; for reasons I don't
understand, however, with lme it seems better to multiply the dye swaps, and
fit a constant term that to fit the alternating (1, -1)).
The following code allows to play with the above (you need to use helmert
contrasts). I used an example with larger number of technical replicates to
better understand the results.
f2 <- function(effect.size = 0,
s.biol = 2, s.tech = 1,
biological.repls = 5, n.tech.reps = 10){
if(getOption("contrasts")[1] != "contr.helmert")
stop("You need to set Helmert contrasts in options")
if(n.tech.reps %% 2) stop("Need even number of n.tech.reps")
## simulate data
tech.rep.effect <- rnorm(biological.repls, mean = 0, sd = s.biol)
tech.rep.effect <- rep(tech.rep.effect,
rep(n.tech.reps, biological.repls))
data.l <- data <-
rnorm(biological.repls * n.tech.reps, mean = 0, sd = s.tech) +
tech.rep.effect + effect.size
dim(data.l) <- c(1, biological.repls * n.tech.reps)
## design matrices, etc.
tech.reps <-
factor(rep(1:biological.repls, rep(n.tech.reps, biological.repls)))
dm1 <- model.matrix(data ~ tech.reps)
dye.swap.codes <- rep(c(-1, 1), (n.tech.reps * biological.repls)/2)
data.sp <- data.l * dye.swap.codes
data.s <- data * dye.swap.codes
## model fitting
tmp1 <- lm.series(data.l, dm1)
tmp2 <- lm.series(data.sp, dm1 * dye.swap.codes)
tmp3 <- summary(lme(data ~ 1, random = ~ 1|tech.reps,
control = list(tolerance = 1e-10, msTol = 1e-10)))
tmp4 <- summary(lme(data.s ~ -1 + dye.swap.codes,
random = ~ 1|tech.reps,
control = list(tolerance = 1e-10, msTol = 1e-10)))
tmp1 <- c(tmp1$coefficients[[1]], tmp1$stdev.unscaled[[1]],
tmp1$sigma, tmp1$sigma *tmp1$stdev.unscaled[[1]])
names(tmp1) <- c("coeff", "unscaled.stdev", "sigma", "scaled se")
tmp2 <- c(tmp2$coefficients[[1]], tmp2$stdev.unscaled[[1]],
tmp2$sigma, tmp2$sigma *tmp2$stdev.unscaled[[1]])
rbind(
lm.series.positivized = tmp1,
lm.series.swap.design = tmp2,
lme.positivized = c(tmp3$coefficients[[1]],
sqrt(tmp3$varFix)/tmp3$sigma, tmp3$sigma, sqrt(tmp3$varFix)),
lme.swap.design = c(tmp4$coefficients[[1]],
sqrt(tmp4$varFix)/tmp4$sigma, tmp4$sigma, sqrt(tmp4$varFix)))
}
Best,
Ramón
--
Ramón Díaz-Uriarte
Bioinformatics Unit
Centro Nacional de Investigaciones Oncológicas (CNIO)
(Spanish National Cancer Center)
Melchor Fernández Almagro, 3
28029 Madrid (Spain)
Fax: +-34-91-224-6972
Phone: +-34-91-224-6900
http://bioinfo.cnio.es/~rdiaz
More information about the Bioconductor
mailing list