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

###################################################
### code chunk number 1: install (eval = FALSE)
###################################################
## source("http://bioconductor.org/biocLite.R")
## biocLite("qusage")


###################################################
### code chunk number 2: setup
###################################################
options(width=70)
options(SweaveHooks=list(fig=function()
             par(mar=c(5.1, 4.1, 1.1, 2.1))))

library(qusage)
data("fluExample")
data("GeneSets")


###################################################
### code chunk number 3: exprSimple
###################################################
dim(eset)
eset[1:5,1:5]


###################################################
### code chunk number 4: labelsSimple
###################################################
labels = c(rep("t0",17),rep("t1",17))
labels


###################################################
### code chunk number 5: contrastSimple
###################################################
contrast = "t1-t0"


###################################################
### code chunk number 6: geneSetSimple
###################################################
ISG.geneSet[1:10]


###################################################
### code chunk number 7: qsSimple
###################################################
qs.results = qusage(eset, labels, contrast, ISG.geneSet)


###################################################
### code chunk number 8: pValSimple
###################################################
pdf.pVal(qs.results)


###################################################
### code chunk number 9: plotCodeSimple (eval = FALSE)
###################################################
## plot(qs.results, xlab="Gene Set Activation")


###################################################
### code chunk number 10: plotSimple
###################################################
plot(qs.results, xlab="Gene Set Activation")


###################################################
### code chunk number 11: read.MSIG.geneSets (eval = FALSE)
###################################################
##   MSIG.geneSets = read.gmt("c2.cp.kegg.v4.0.symbols.gmt")


###################################################
### code chunk number 12: MSIG.geneSets
###################################################
  summary(MSIG.geneSets[1:5])
  MSIG.geneSets[2]

  qs.results.msig = qusage(eset, labels, contrast, MSIG.geneSets)
  numPathways(qs.results.msig)


###################################################
### code chunk number 13: multiHypCorr
###################################################
  p.vals = pdf.pVal(qs.results.msig)
  head(p.vals)
  q.vals = p.adjust(p.vals, method="fdr")
  head(q.vals)


###################################################
### code chunk number 14: qusage.Rnw:143-144
###################################################
  options(width=100)


###################################################
### code chunk number 15: qsTable
###################################################
  qsTable(qs.results.msig, number=10)


###################################################
### code chunk number 16: qusage.Rnw:149-150
###################################################
  options(width=70)


###################################################
### code chunk number 17: plotCIs
###################################################
par(mar=c(10,4,1,2))
plot(qs.results.msig)


###################################################
### code chunk number 18: pairVector
###################################################
  pairs = c(1:17,1:17)
  pairs


###################################################
### code chunk number 19: twoWayLabels
###################################################
  resp
  twoWay.labels = paste(resp, labels, sep=".")
  twoWay.labels


###################################################
### code chunk number 20: twoWayRun
###################################################
  sx.results = qusage(eset, twoWay.labels, "sx.t1-sx.t0", MSIG.geneSets, pairVector=pairs)
  asx.results = qusage(eset, twoWay.labels, "asx.t1-asx.t0", MSIG.geneSets, pairVector=pairs)


###################################################
### code chunk number 21: twoCurvePval
###################################################
  p.vals = twoCurve.pVal(sx.results,asx.results)
  head(p.vals)


###################################################
### code chunk number 22: twoWayDensityCurvesCode
###################################################
  plotDensityCurves(asx.results,path.index=125,col="black",xlim=c(-0.2,0.5),main="CYTOSOLIC DNA SENSING")
  plotDensityCurves(sx.results,path.index=125,col="red",add=T)
  abline(v=0, lty=2)
  legend("topleft", legend=c("asx","sx"),lty=1,col=1:2)


###################################################
### code chunk number 23: twoWayDensityCurves
###################################################
  plotDensityCurves(asx.results,path.index=125,col="black",xlim=c(-0.2,0.5),main="CYTOSOLIC DNA SENSING")
  plotDensityCurves(sx.results,path.index=125,col="red",add=T)
  abline(v=0, lty=2)
  legend("topleft", legend=c("asx","sx"),lty=1,col=1:2)


###################################################
### code chunk number 24: plotGSD
###################################################
plotGeneSetDistributions(sx.results,asx.results,path.index=125)


###################################################
### code chunk number 25: plotCIsGenes
###################################################
par(mar=c(8,4,1,2))
plotCIsGenes(sx.results,path.index=125)


###################################################
### code chunk number 26: prettyPDF
###################################################
##select the best 5 gene sets
p.vals = pdf.pVal(qs.results.msig)
topSets = order(p.vals)[1:5]

##plot 
plotDensityCurves(qs.results.msig, path.index=topSets, col=1:5, 
                  lwd=2, xlab=expression(log[2](Activity)))


###################################################
### code chunk number 27: plotPDF
###################################################
##select the best 5 gene sets
p.vals = pdf.pVal(qs.results.msig)
topSets = order(p.vals)[1:5]

##plot 
plotDensityCurves(qs.results.msig, path.index=topSets, col=1:5, 
                  lwd=2, xlab=expression(log[2](Activity)))


###################################################
### code chunk number 28: createCDF
###################################################
##get coordinates
coord.list = plotDensityCurves(qs.results, plot=FALSE)
coords = coord.list[[1]]
x=coords$x
##Calculate the CDF as the integral of the PDF at each point.
y = cumsum(coords$y)/sum(coords$y)

plot(x,y, type="l", xlim=calcBayesCI(qs.results,low=0.001)[,1])


###################################################
### code chunk number 29: plotCDF
###################################################
##get coordinates
coord.list = plotDensityCurves(qs.results, plot=FALSE)
coords = coord.list[[1]]
x=coords$x
##Calculate the CDF as the integral of the PDF at each point.
y = cumsum(coords$y)/sum(coords$y)

plot(x,y, type="l", xlim=calcBayesCI(qs.results,low=0.001)[,1])


###################################################
### code chunk number 30: plotISG_GSD
###################################################
plotGeneSetDistributions(qs.results)


###################################################
### code chunk number 31: plotGSDHighlight
###################################################
path = 125 ##the CYTOSOLIC_DNA_SENSING pathway
path.gene.i = sx.results$pathways[[path]] ##the genes are stored as indexes
path.genes = rownames(eset)[path.gene.i]

cols = rep("#33333399",length(path.genes))
##make the IFNA genes red, and the IFNB genes blue.
cols[path.genes=="CXCL10"] = "#FF0000"

plotGeneSetDistributions(sx.results,asx.results,path.index=125, 
                         groupLabel=c("Symptomatic","Asymptomatic"),
                         colorScheme=cols, barcode.col=cols)