################################################### ### chunk number 1: loadpacks ################################################### library(Biobase,warn.conflicts = FALSE) library(affy,warn.conflicts = FALSE) library(ctest) library(multtest) library(bioclabs) library(annotate) library(hgu95av2)

################################################### ### chunk number 2: loaddata ################################################### data("spikein")

################################################### ### chunk number 3: histograms ################################################### pops <- pData(spikein)[,1]+2 hist(spikein,col=pops,type="l")

################################################### ### chunk number 4: rma ################################################### eset <- rma(spikein)

################################################### ### chunk number 5: howmanypops ################################################### eset$population Index1 <- which(eset$population==0) Index2 <- which(eset$population==1)

################################################### ### chunk number 6: scores ################################################### scores <- esApply(eset,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) })

################################################### ### chunk number 7: scores2 ################################################### scores <- t(scores) ##put genes back in rows colnames(scores) <- c("a","m","t.test","p.value")

################################################### ### chunk number 8: MvA ################################################### plot(scores[,1],scores[,2],xlab="A",ylab="M",pch=".") abline(h=c(-1,1)) ##fc = 2

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

################################################### ### chunk number 10: pvalues ################################################### sum(scores[,4]<=0.05) sum(scores[,4]<=0.01)

################################################### ### chunk number 11: adjustedpvalues ################################################### tmp <- mt.rawp2adjp(scores[,4],proc="BH") adj.p.values <- tmp$adjp[order(tmp$index),] scores <- cbind(scores,adj.p.values[,-1]) colnames(scores)[5]<- "FDR"

################################################### ### chunk number 12: significant ################################################### Names <- geneNames(eset)[scores[,5]<=0.01] Names <- Names[order(scores[Names,5])]

################################################### ### chunk number 13: annotation ################################################### ll <- multiget(Names,env=hgu95av2LOCUSID) sym <- multiget(Names,env=hgu95av2SYMBOL)

################################################### ### chunk number 14: htmlreport ################################################### res <- data.frame(unlist(sym),signif(scores[Names,],2)) ll.htmlpage(ll,filename="report.html", title="HTML report", othernames=res, table.head=c("Locus ID","Gene Symbol",colnames(scores)), table.center=TRUE)

################################################### ### chunk number 15: truepositives ################################################### true <- colnames(pData(eset))[-1] tp <- sum(Names%in%true) cat(tp,"true positives",length(Names)-tp,"false positives.\n")

################################################### ### chunk number 16: MvA2 ################################################### plot(scores[,1],scores[,2],xlab="A",ylab="M",pch=".",ylim=c(-1,1)) text(scores[Names,1],scores[Names,2],sym,pch=16,col=rainbow(length(Names))) abline(h=c(-1,1)) ##fc = 2 fp <- Names[!Names%in%true] ##lets show the false positves points(scores[fp,1],scores[fp,2],pch=4,cex=4,col="red")

################################################### ### chunk number 17: 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)[true,]