### R code from vignette source 'GSReg.Rnw'

###################################################
### code chunk number 1: start
###################################################
options(width=85)
options(continue=" ")
#rm(list=ls())


###################################################
### code chunk number 2: GSReg.Rnw:166-171
###################################################
library(GSBenchMark)
data(diracpathways)
class(diracpathways)
names(diracpathways)[1:5]
class(diracpathways[[1]])


###################################################
### code chunk number 3: GSReg.Rnw:180-182
###################################################
data(GSBenchMarkDatasets)
print(GSBenchMark.Dataset.names)


###################################################
### code chunk number 4: GSReg.Rnw:188-191
###################################################
DataSetStudy = GSBenchMark.Dataset.names[[9]]
print(DataSetStudy)
data(list=DataSetStudy)


###################################################
### code chunk number 5: GSReg.Rnw:200-202
###################################################
if(sum(apply(is.nan(exprsdata),1,sum)>0))
  exprsdata = exprsdata[-which(apply(is.nan(exprsdata),1,sum)>0),];


###################################################
### code chunk number 6: GSReg.Rnw:206-208
###################################################
genenames = rownames(exprsdata);
genenames[1:10]


###################################################
### code chunk number 7: GSReg.Rnw:238-239
###################################################
library(GSReg)


###################################################
### code chunk number 8: GSReg.Rnw:245-247
###################################################
Nperm = 10
system.time({DIRACAn =GSReg.GeneSets.DIRAC(exprsdata,diracpathways,phenotypes,Nperm=Nperm)})


###################################################
### code chunk number 9: GSReg.Rnw:250-251
###################################################
hist(DIRACAn$pvalues,xlab="pvalue",main="Hist of pvalues applying DIRAC Analysis.")


###################################################
### code chunk number 10: GSReg.Rnw:272-281
###################################################
#Calculating the variance for the pathways
#Calculate how much it takes to calculate the statistics and their p-value for all pathways

system.time({VarAnKendallV = GSReg.GeneSets.EVA(geneexpres=exprsdata,
  pathways=diracpathways, phenotypes=as.factor(phenotypes)) })


names(VarAnKendallV)[[1]]
VarAnKendallV[[1]]


###################################################
### code chunk number 11: GSReg.Rnw:292-309
###################################################
Nperm = 10;
VarAnPerm = vector(mode="list",length=Nperm)
for( i in seq_len(Nperm))
{
  VarAnPerm[[i]] = GSReg.GeneSets.EVA(geneexpres=exprsdata, pathways=diracpathways, 
                                           phenotypes=sample(phenotypes))
}

pvaluesperm = vector(mode="numeric",length=length(VarAnPerm[[1]]))

for( i in seq_along(VarAnPerm[[1]]))
{
  z = sapply(VarAnPerm,function(x) x[[i]]$E1 - x[[i]]$E2)
  pvaluesperm[i] = mean(abs(VarAnKendallV[[i]]$E1-VarAnKendallV[[i]]$E2)<abs(z)) 
}
 zscore = sapply(VarAnKendallV,function(x) x$zscore);
 pvalustat = sapply(VarAnKendallV,function(x) x$pvalue);


###################################################
### code chunk number 12: GSReg.Rnw:314-323
###################################################
 plot(x=abs(zscore),y=pvaluesperm,xlab="|Z-score|",
  ylab="p-value",col="red1",main="p-value comparisons")
 zscorelin = seq(0,6,0.1);
 pvaltheoretic = (1-pnorm(zscorelin))*2
 lines(x=zscorelin,y=pvaltheoretic,type="l",pch=50,lty=5,col="darkblue")

  legend("topright",legend=c("permutation test","U-Stat Estimation"),
         col=c("red","blue"),text.col=c("red","blue"),
  lty=c(NA,1),lwd=c(NA,2.5),pch=c(21,NA))


###################################################
### code chunk number 13: GSReg.Rnw:330-331
###################################################
hist(x=pvalustat,breaks=20,main="P-value Hist of U-Stat",xlim=c(0,1))


###################################################
### code chunk number 14: GSReg.Rnw:344-350
###################################################
plot(x=DIRACAn$pvalues,y=pvalustat,xlab ="DIRAC",
        ylab="EVA",main=sprintf("P-value Comparison corr=%2.2g",cor(x=DIRACAn$pvalues,y=pvalustat)))

lmfit = lm(pvalustat~DIRACAn$pvalues-1)
abline(lmfit)
cor.test(x=DIRACAn$pvalues,y=pvalustat)


###################################################
### code chunk number 15: GSReg.Rnw:354-355
###################################################
cor(x=DIRACAn$pvalues,y=pvalustat)


###################################################
### code chunk number 16: GSReg.Rnw:361-369
###################################################
load("DIRACAn.rda")
significantPathwaysDIRAC = names(DIRACAn$mu1)[which(DIRACAn$pvalues<0.05)];
print(significantPathwaysDIRAC)
mu1 = DIRACAn$mu1[significantPathwaysDIRAC];
mu2 = DIRACAn$mu2[significantPathwaysDIRAC];
significantPathwaysGSV = names(which(pvalustat<0.05));
eta1 = sapply(VarAnKendallV,function(x) x$E1)[significantPathwaysGSV];
eta2 = sapply(VarAnKendallV,function(x) x$E2)[significantPathwaysGSV];


###################################################
### code chunk number 17: GSReg.Rnw:377-381
###################################################
  plot(x=DIRACAn$pvalues,y=pvalustat,xlab ="DIRAC",
        ylab="EVA",main=sprintf("P-value Comparison corr=%2.2g",cor(x=DIRACAn$pvalues,y=pvalustat)))
  lmfit = lm(pvalustat~DIRACAn$pvalues-1)
abline(lmfit)


###################################################
### code chunk number 18: GSReg.Rnw:397-398 (eval = FALSE)
###################################################
## DIRACAn =GSReg.GeneSets.DIRAC(exprsdata,diracpathways,phenotypes,Nperm=1000)


###################################################
### code chunk number 19: GSReg.Rnw:401-411
###################################################
significantPathwaysDIRAC = names(DIRACAn$mu1)[which(DIRACAn$pvalues<0.05)];
mu1 = DIRACAn$mu1[significantPathwaysDIRAC];
mu2 = DIRACAn$mu2[significantPathwaysDIRAC];

#The dysregulated pathways
names(mu1)
plot(x=mu1,y=mu2,
     xlim=c(0,max(mu1,mu2)),ylim=c(0,max(mu1,mu2)),xlab="normal",ylab="disease",
     main="(a)  DIRAC significantly dysregulated pathways")
lines(x=c(0,max(mu1,mu2)),y=c(0,max(mu1,mu2)))


###################################################
### code chunk number 20: GSReg.Rnw:426-435
###################################################
significantPathwaysGSV = names(which(pvalustat<0.05));
eta1 = sapply(VarAnKendallV,function(x) x$E1)[significantPathwaysGSV];
eta2 = sapply(VarAnKendallV,function(x) x$E2)[significantPathwaysGSV];

#The dysregulated pathways
names(eta1)
plot(x=eta1,y=eta2,xlim=c(0,max(eta1,eta2)),ylim=c(0,max(eta1,eta2)),xlab="normal",ylab="disease",
     main="(b) EVA: Dysregulated pathways")
lines(x=c(0,max(eta1,eta2)),y=c(0,max(eta1,eta2)))


###################################################
### code chunk number 21: GSReg.Rnw:441-443
###################################################
print(significantPathwaysGSV)
print(significantPathwaysDIRAC)


###################################################
### code chunk number 22: sessioInfo
###################################################
toLatex(sessionInfo())