QuantileFirstPrinciples.R

1. ===============================================================================
2. Generating Example Data for quantile normalization
3. ===============================================================================

y <- cbind(ProbeSet = rep(1:2, each=3), Probe = rep(1:3, 2),

          Chip1 = c(7,3,2,4,10,12), Chip2 = c(9,5,6,8,11,10),
          Chip3 = c(19,14,11,8,16,15))

library(limma) ynorm <- cbind(y[,1:2],normalizeQuantiles(y[,3:5], ties=T))

1. ===============================================================================
2. Generating Example Data for median polish
3. ===============================================================================

y <- cbind(probe1=c(18,13,15,19), probe2=c(11,7,6,15), probe3=c(8,5,7,12),

          probe4=c(21,16,16,18), probe5=c(4,7,6,5))

E <- y

1. ===============================================================================
2. Median polish by first principles (3x), convergence is fast
3. ===============================================================================

FUN <- median

1. Sweep out medians on rows

meds <- apply(E, 1, FUN) E <- sweep(E, 1, STATS=meds)

1. Sweep out medians on columns

meds <- apply(E, 2, FUN) E <- sweep(E, 2, STATS=meds)

1. Sweep out medians on rows

meds <- apply(E, 1, FUN) E <- sweep(E, 1, STATS=meds)

1. Sweep out medians on columns

meds <- apply(E, 2, FUN) E <- sweep(E, 2, STATS=meds)

1. Sweep out medians on rows

meds <- apply(E, 1, FUN) E <- sweep(E, 1, STATS=meds)

1. Sweep out medians on columns

meds <- apply(E, 2, FUN) E <- sweep(E, 2, STATS=meds)

1. see ?medpolish for R function

medpolish(y, trace.iter=F)$res

1. ===============================================================================
2. Fitted values
3. ===============================================================================

mu <- y - E

1. ===============================================================================
2. Y_ijk = mu_i.k + alpha_.jk + Eijk -Gene chip specific probe affinities mu_i.k
3. i = 1..I chips, j = 1..J probes, k = 1..G genes
4. ===============================================================================

1. mu_i.k calculation

mu.i <- apply(mu,1, mean)

1. alpha_.jk calculation

alpha <- y - mu.i - E alpha.j <- apply(alpha,2, mean)

1. ===============================================================================
2. Mean instead of median decomposition using anova, aov
3. ===============================================================================

dset <- as.data.frame(cbind(y = c(y), probe = rep(1:5, each=4), array = rep(1:4, 5))) dset[,"probe"] <- factor(dset[,"probe"]) dset[,"array"] <- factor(dset[,"array"])

fit <- aov(y ~ probe + array, data=dset) summary(fit)

matrix(fit$resid, nc=5)