################################################### ### chunk number 1: loadpacks ################################################### library(Biobase) library(ctest) library(multtest) library(bioclabs)

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

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

################################################### ### chunk number 4: scores ################################################### Index1 <- which(eset3$population==0) Index2 <- which(eset3$population==1) scores <- esApply(eset3,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: MvA ################################################### plot(scores[,1],scores[,2],xlab="A",ylab="M",pch=".") #points(scores[genenames,1],scores[genenames,2],pch=16,col=1:16) text(scores[genenames,1],scores[genenames,2],genenames,col=1:16) abline(h=c(-1,1)) ##fc = 2

################################################### ### chunk number 6: volcano ################################################### plot(scores[,2],abs(scores[,3]),xlab="M",ylab="t.stat",pch=".") points(scores[genenames,2],abs(scores[genenames,3]),pch=16,col=1:16) abline(v=c(-1,1)) a <- qt(.975,4) abline(h=a)

################################################### ### chunk number 7: volcanocloseup ################################################### plot(scores[,2],abs(scores[,3]),xlab="M", ylab="t.stat",pch=".",xlim=c(-1.5,1.5),ylim=c(0,15)) text(scores[genenames,2],abs(scores[genenames,3]),genenames,col=1:16) abline(v=c(-1,1)); a <- qt(.975,4); abline(h=a)

################################################### ### chunk number 8: pvalues ################################################### sum(scores[,"p.value"]<=0.05) sum(scores[,"p.value"]<=0.01)

################################################### ### chunk number 9: adjustedpvalues ################################################### tmp <- mt.rawp2adjp(scores[,4]) adj.p.values <- tmp$adjp[order(tmp$index),-1] scores <- cbind(scores,adj.p.values)

################################################### ### chunk number 10: maxT ################################################### tmp <- mt.maxT(exprs(eset3),eset3$population) scores <- cbind(scores,tmp$adjp[order(tmp$index)]) colnames(scores)[12] <- "maxT"

################################################### ### chunk number 11: results ################################################### round(scores[genenames,],2)

################################################### ### chunk number 12: ranks ################################################### m.ranks <- rank(-abs(scores[,2])) ##by fold change names(m.ranks) <- rownames(scores) t.ranks <- rank(-abs(scores[,3])) ##by fold change names(t.ranks) <- rownames(scores) cbind(m.ranks,t.ranks)[genenames,]