## data leuk <- scale(read.table("golubtrain.txt", head = TRUE), center = TRUE, scale = FALSE) ## create z-scores ix.all <- 1:27 ix.aml <- 27 + 1:11 m1 <- length(ix.all) m2 <- length(ix.aml) my.t.test <- function(x) {mm <- t.test(x[ix.all], x[ix.aml], var.equal = TRUE); return(qnorm(pt(mm$stat, mm$par)))} z <- apply(leuk, 1, my.t.test) iz <- order(z) leuk <- leuk[iz,] clone.id <- dimnames(leuk)[[1]] z <- z[iz] n <- nrow(leuk) ### data visualization #rgb.col <- rgb(pmax(seq(-1,1,.2),0), pmax(seq(1,-1,-.2), 0), 0) #image(1:(m1 + m2), 1:n, scale(t(leuk)), col = rgb.col, ann = FALSE) #abline(v = m1 + .5, col = "white", lwd = 0.5) #title(ylab = "genes", xlab = "subjects") #mtext(c("ALL", "AML"), side = 3, at = c(m1/2, m1 + m2/2)) ## discoveries in Lee et al (2003) and Golub et al (1999) lee.list <- c(1882, 760, 2288, 4847, 1144, 1120, 4535, 6218, 6200, 1834, 1630, 5772, 1745, 804, 2354, 3252, 6201, 1685, 6041, 1779, 6855, 173, 2642, 1829, 4107, 697, 229) lee.wt <- c(7.7, 4.8, 2.8, 2.8, 2.7, 2.4, 2.1, 2.0, 1.9, 1.9, 1.9, 1.8, 1.7, 1.7, 1.6, 1.5, 1.3, 1.3, 1.3, 1.3, 1.3, 1.3, 1.2, 1.2, 1.2, 1.2, 1.2) golub.list <- c("U22376", "X59417", "U05259", "M92287", "M31211", "X74262", "D26156", "S50223", "M31523", "L47738", "U32944", "Z15115", "X15949", "X63469", "M91432", "U29175", "Z69881", "U26266", "M31303", "Y08612", "U35451", "M29696", "M13792", "M55150", "X95735", "U50136", "M16038", "U82759", "M23197", "M84526", "Y12670", "M27891", "X17042", "Y00787", "M96326", "U46751", "M80254", "L08246", "M62762", "M28130", "M63138", "M57710", "M69043", "M81695", "X85116", "M19045", "M83652", "X04085") ## PRML multiple-test fit source("prmlmtest.R") z.prml <- prml.mtest.adapt(z, nperm = 10, start = "random", opt.exclude = c("rho")) summary(z.prml) #$par #mu sigma p0 tau gam #0.2315516 1.0408235 0.6326913 3.0516075 0.6700000 # #$logLik #[1] -13202.41 # ## visualization z.fdr <- plot(z.prml, fdr.cut = 0.1) points(z[match(lee.list, iz)], lee.wt / 100, ty = "h", cex = 0.5, col = "darkblue") points(z[charmatch(golub.list, clone.id)], rep(0, length(golub.list)), pch = 19, cex = 0.4, col = "darkred") ## Empirical null estimated by others plot(z.prml, fdr.cut = 0.1) z.jc <- jincai(z) overlay(z.prml, z.jc, col = 2) z.ef <- efron(z) overlay(z.prml, z.ef, col = 3) z.mu <- murali(z) overlay(z.prml, z.mu, col = 4) legend("topleft", c("PRML", "Jin-Cai", "Efron", "Murali"), col = 1:4, bty = "n", lty = 1)