## data hed <- read.table("hedenfalk-breast-data_filt.txt", head = TRUE) ## create z-scores ix1 <- 1 + c(1:6, 18) ix2 <- 1 + c(7:10, 19:22) m1 <- length(ix1) m2 <- length(ix2) hed.12 <- hed[, c(ix2, ix1)] my.t.test <- function(x) {mm <- t.test(x[1:m2], x[m2 + 1:m1]); return(qnorm(pt(mm$stat, mm$par)))} z <- apply(hed.12, 1, my.t.test) iz <- order(z) hed.12 <- hed.12[iz,] clone.id <- hed[iz,1] z <- z[iz] n <- nrow(hed.12) # ## 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(hed.12)), col = rgb.col, ann = FALSE) # abline(v = m2 + .5, col = "white", lwd = 0.5) # title(ylab = "genes", xlab = "subjects") # mtext(c("BRCA2", "BRCA1"), side = 3, at = c(m2/2, m2 + m1/2)) ## discoveries from Lee et al (2003) lee.list <- c(897781, 823940, 26184, 840702, 376516, 47542, 366647, 293104, 28012, 212198, 247818, 26082, 667598, 30093, 73531, 950682, 47681, 46019, 307843, 548957, 788721, 843076, 204897, 812227, 566887, 563598, 324210) lee.wt <- c(8.6, 8.4, 7.8, 7.5, 7.1, 6.9, 6.6, 6.6, 6.2, 6.1, 5.9, 5.5, 5.4, 5.2, 5.1, 5, 5, 4.9, 4.9, 4.8, 4.7, 4.7, 4.7, 4.7, 4.6, 4.6, 4.5) ## 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.007222966 0.962118417 0.462049874 1.739806878 0.670000000 # #$logLik #[1] -5621.06 ## visualization z.fdr <- plot(z.prml, fdr.cut = 0.1) points(z[match(lee.list, clone.id)], lee.wt / 100, ty = "h", cex = 0.5, col = "darkblue") ## 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)