##--------------------------------------------------------------- ## densities ##--------------------------------------------------------------- library("vsn") library("geneplotter") data("lymphoma") ll=function(x, thr=2)log2(ifelse(x>thr, x, thr)) multidensity(ll(exprs(lymphoma)), col=brewer.pal(9, "Set1")[ifelse(lymphoma$dye=="Cy3", 3, 1)], lwd=3, main="8 arrays from the lymphoma data (Alizadeh 2000)", xlab=expression(log[2]~intensity)) ##--------------------------------------------------------------- ## quantile normalisation ##--------------------------------------------------------------- library("affy") library("affydata") library("preprocessCore") library("RColorBrewer") library("geneplotter") data("Dilution") nr = apply(exprs(Dilution), 2, rank) nq = normalize.quantiles(exprs(Dilution)) colnames(nq) = sampleNames(Dilution) cols = brewer.pal(9,"Set1") qn1 = function(x, ...) { matplot(nr, x, pch=".", log="y", xlab="rank", col=cols) legend("topleft", colnames(x), fill=cols) } qn2 = function(x, ...) multidensity(log2(x), lwd=2, col=cols, ...) png(file="110628-brixen-microarrays-huber-quantilenorm.png", width=1280, height=960, pointsize=24) par(mfrow=c(2,2)) qn1(exprs(Dilution), main="before") qn1(nq, main="after quantile normalisation") qn2(exprs(Dilution)) qn2(nq) dev.off() ##--------------------------------------------------------------- ## robust regression ##--------------------------------------------------------------- x = runif(10)*10 y = 2 + 2*x + rnorm(length(x)) wm = which.min(x) y[wm] = y[wm]+10 plot(x,y, pch=16) f1 = lm(y~x) f2 = rlm(y~x) cols = c("#E41A1C", "#377EB8") lines(x, predict(f1), lwd=2, col=cols[1]) lines(x, predict(f2), lwd=2, col=cols[2]) legend("bottomright", c("lm", "rlm"), col=cols, lwd=2, lty=1) ##--------------------------------------------------------------- ## two component error model ##--------------------------------------------------------------- pdf(file="110628-brixen-microarrays-huber-twocomponent.pdf", width=8, height=4) par(mfrow=c(1,2)) mean = seq(0, 1000, length=40) a = 1000; b=.01 variance = a + b*mean^2 varmax = variance[length(variance)] plot(mean, variance, type="l", lwd=2, ylim=c(0, varmax)) abline(h=a, col="blue", lty=2) text(0, varmax, substitute(var==a+b~mean^2, list(a=a,b=b)), adj=c(0,1)) plot(mean, sqrt(variance), type="l", lwd=2, ylim=c(0, sqrt(varmax))) abline(h=sqrt(a), col="blue", lty=2) abline(a=0, b=sqrt(b), col="orange", lty=2) dev.off() ##--------------------------------------------------------------- ## Tukey Biweight ##--------------------------------------------------------------- x = seq(1, 100, length=100) m = 50 s = 18 u = (x-m)/s w = ifelse(abs(u)>1, 0, ((1 - u^2)^2)) plot(x, w, main="Tukey Biweight", type="l") ##--------------------------------------------------------------- ## Median Polish ##--------------------------------------------------------------- x1 = c(rnorm(13), 8+rnorm(12)) x2 = x1 x2[1] = x2[1]+8 rg = range(c(x1,x2)) myhist=function(x){ hist(x, xlim=rg, breaks=12, col="orange", main="") med = median(x) trm = mean(x, trim=0.15) abline(v=c(med, trm), col=cols, lwd=2) } par(mfrow=c(1,2)) myhist(x1) myhist(x2)