### R code from vignette source 'vignettes/PECA/inst/doc/PECA.Rnw'

###################################################
### code chunk number 1: load
###################################################
library(PECA)
library(SpikeIn)
data(SpikeIn133)


###################################################
### code chunk number 2: subset
###################################################
data <- SpikeIn133[,c(1,15,29,2,16,30)]


###################################################
### code chunk number 3: peca
###################################################
peca_results <- PECA_AffyBatch(normalize="true",affy=data)


###################################################
### code chunk number 4: mas
###################################################
library(multtest)
mas5 <- mas5(data)
groups=c(1,1,1,0,0,0)
stats=mt.teststat(exprs(mas5), groups, test="t")
mas5_results=2*(1-pnorm(abs(stats)))


###################################################
### code chunk number 5: rma
###################################################
rma <- rma(data)
groups=c(1,1,1,0,0,0)
stats=mt.teststat(exprs(rma), groups, test="t")
rma_results=2*(1-pnorm(abs(stats)))


###################################################
### code chunk number 6: combine
###################################################
results <- cbind(peca=peca_results[,4],mas5=mas5_results,rma=rma_results,spike=0)
results[,4][rownames(results) %in% colnames(pData(data))] <- 1
head(results)


###################################################
### code chunk number 7: roc
###################################################
library(ROCR)
results_roc <- results
results_roc[,1:3] <- 1-results_roc[,1:3]
pred.peca <- prediction(results_roc[,1],results_roc[,4])
pred.mas5 <- prediction(results_roc[,2],results_roc[,4])
pred.rma <- prediction(results_roc[,3],results_roc[,4])
perf.peca <- performance(pred.peca,"tpr","fpr")
perf.mas5 <- performance(pred.mas5,"tpr","fpr")
perf.rma <- performance(pred.rma,"tpr","fpr")
plot(perf.mas5, lwd=3, col="green")
plot(perf.rma, lwd=3, col="blue", add=TRUE)
plot(perf.peca, lwd=3, col="red", add=TRUE)
legend(0.7,0.2,c('PECA', 'RMA', 'MAS'), col=c('red','blue','green'), lwd=3)