[BioC] dye swaps of technical replicates and variable numbers of replicate spots

[Ramon Diaz-Uriarte]

> 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