################################################### ### chunk number 1: loadpacks ################################################### library(Biobase) library(golubEsets) library(ctest) library(multtest) library(bioclabs) library(modreg) library(wavethresh) library(Milan)

################################################### ### chunk number 2: loaddata ################################################### data("eset12") genenames <- colnames(pData(eset12))[-1]

################################################### ### chunk number 3: phenodata ################################################### pData(eset12)

################################################### ### chunk number 4: scores ################################################### Index1 <- which(eset12$population==0) Index2 <- which(eset12$population==1) scores <- esApply(eset12,1,function(x){ tmp <- t.test(x[Index2],x[Index1],var.equal=TRUE) c(mean(tmp$estimate),-diff(tmp$estimate),tmp$statistic,tmp$p.value) }) scores <- t(scores) ##put genes back in rows colnames(scores) <- c("A","M","t.stat","p.value")

################################################### ### chunk number 5: dwt ################################################### yy<-scores[1:8192,2] plot(acf(yy))

################################################### ### chunk number 6: dwt2 ################################################### wds1 <- wd(yy) par(mfcol = c(4, 2),mgp = c(5, 0.4, 0)) plot(acf(accessD(wds1,level=11))) plot(acf(accessD(wds1,level=10))) plot(acf(accessD(wds1,level=9))) plot(acf(accessD(wds1,level=8))) plot(acf(accessD(wds1,level=7))) plot(acf(accessD(wds1,level=6))) plot(acf(accessD(wds1,level=5))) plot(acf(accessD(wds1,level=4)))

################################################### ### chunk number 7: bayes ################################################### y<-scores[,3] fit<- ebayesthresh(y, prior = "cauchy", a = NA, bayesfac = T,sdev = NA, verbose = F, threshrule = "median") aux<- abs(fit) auxsv <- aux > 0 cat("Number of signicant genes: \n") sum(as.integer(auxsv)) rownames(scores[aux>0,])

################################################### ### chunk number 8: results ################################################### genenames

################################################### ### chunk number 9: loaddatag ################################################### data(golub) # the golub data set classlabel<-golub.cl dat=golub o <- rep(T,nrow(dat)) for(i in 1:nrow(dat)) { if(max(dat[i,])>20) {o[i]<-F} } dat <- dat[o,] n <- ncol(dat) m <- nrow(dat) y=golub.cl y=y+1 #makes two-sample t-statistics ttest <- function(x,y,f=0) { m1 <- apply(x[,y==1],1,mean) m2 <- apply(x[,y==2],1,mean) v1 <- apply(x[,y==1],1,var) v2 <- apply(x[,y==2],1,var) s1 <- sqrt(v1/sum(y==1)) s2 <- sqrt(v2/sum(y==2)) tt <- (m2-m1)/(sqrt(s1^2 + s2^2)+f) return(tt=tt, ss=(sqrt(s1^2 + s2^2)+f)) } tt <- ttest(dat,y)$tt hist(tt,breaks=100,col="red")

################################################### ### chunk number 10: pvalues ################################################### tt0 <- 0 set.seed(123) B <- 100 for(i in 1:B) { v <- sample(y) tt0 <- c(tt0,ttest(dat,v)$tt) } tt0 <- tt0[-1]

#form p-values att <- abs(tt) att0 <- abs(tt0) v <- c(rep(T,m),rep(F,mB)) v <- v[rev(order(c(att,att0)))] u <- 1:length(v) w <- 1:m p <- (u[v==TRUE]-w)/(Bm) p <- p[rank(-att)] hist(p,breaks=100,col="blue")

################################################### ### chunk number 11: adjustedqvalues ################################################### qobj <- qvalue(p) #qobj$qval contains q-values #qobj$pval contains p-values #qobj$pi0 contains pi0 estimate

#calculate fold-change (on log base 2 scale) m1 <- apply(dat[,y==1],1,mean) m2 <- apply(dat[,y==2],1,mean) fc <- (m2-m1)

qval <- rep(NA,length(o)) pval <- rep(NA,length(o)) fchng <- rep(NA,length(o)) qval[o] <- qobj$q pval[o] <- qobj$pval fchng[o] <- fc for(i in 1:length(o)) { cat(c(round(qval[i],6),pval[i],round(fchng[i],3),\n),file="qvalues.txt",append=T) }

################################################### ### chunk number 12: Figure1 ################################################### hist(p,nclass=20,main=" ",xlab="Density of observed p-values",ylab=" ",freq=F,ylim=c(0,4)) abline(h=1) abline(h=0.67,lty=2)

################################################### ### chunk number 13: Figure2a ################################################### plot(tt[order(tt)],qobj$q[order(tt)],type="l",xlim=c(-8,6),xlab="t-statistics",ylab="q-values",main="a")

################################################### ### chunk number 14: Figure2b ################################################### p2 <- qobj$pval[order(qobj$pval)] q2 <- qobj$qval[order(qobj$pval)] plot(p2,q2,type="l",xlab="p-values",ylab="q-values",main="b")

################################################### ### chunk number 15: Figure2c ################################################### plot(q2,1:m,type="l",xlim=c(0,0.5),ylim=c(0,3051),xlab="q-value",ylab="number of significant genes",main="c")

################################################### ### chunk number 16: Figure2d ################################################### plot(1:500,q2[1:500]*(1:500),type="l",xlab="number of significant genes",ylab="number of expected false positives",main="d")

################################################### ### chunk number 17: Figure3 ################################################### lam <- seq(0,0.95,0.01) pi0 <- rep(0,length(lam)) for(i in 1:length(lam)) { pi0[i] <- mean(p>lam[i])/(1-lam[i]) } spi0 <- smooth.spline(lam,pi0,df=3,w=(1-lam)) plot(lam,pi0,xlab=expression(lambda),ylab=expression(hat(pi)[0](lambda))) lines(spi0)