### R code from vignette source 'PharmacoGx.Rnw'
### Encoding: ISO8859-1

###################################################
### code chunk number 1: setup (eval = FALSE)
###################################################
## options(keep.source=TRUE)


###################################################
### code chunk number 2: install-pkg (eval = FALSE)
###################################################
## source('http://bioconductor.org/biocLite.R')
## biocLite('PharmacoGx')


###################################################
### code chunk number 3: loadlib (eval = FALSE)
###################################################
## library(PharmacoGx)


###################################################
### code chunk number 4: download-example (eval = FALSE)
###################################################
## ## Example
## availablePSets()
## GDSC <- downloadPSet("GDSC") 


###################################################
### code chunk number 5: download-sig (eval = FALSE)
###################################################
## ## Example
## CMAP.sigs <- downloadPertSig("CMAP")


###################################################
### code chunk number 6: inconsistencies (eval = FALSE)
###################################################
##   library(Biobase)
##   library(PharmacoGx)
##   data("GDSCsmall")
##   data("CCLEsmall")
##   commonGenes <- intersect(fNames(GDSCsmall, "rna"),
##                            fNames(CCLEsmall,"rna"))
##   common <- intersectPSet(list('CCLE'=CCLEsmall,
##                                'GDSC'=GDSCsmall),
##                           intersectOn=c("cell.lines", "drugs"), strictIntersect=TRUE)
## 
##   
##   GDSC.auc <- summarizeSensitivityProfiles(
##                 pSet=common$GDSC,
##                 sensitivity.measure='auc_published', 
##                 summary.stat="median",
##                 verbose=FALSE)
##   CCLE.auc <- summarizeSensitivityProfiles(
##                 pSet=common$CCLE,
##                 sensitivity.measure='auc_published', 
##                 summary.stat="median",
##                 verbose=FALSE)
##   
##   GDSC.ic50 <- summarizeSensitivityProfiles(
##                 pSet=common$GDSC, 
##                 sensitivity.measure='ic50_published', 
##                 summary.stat="median",
##                 verbose=FALSE)
##   CCLE.ic50 <- summarizeSensitivityProfiles(
##                 pSet=common$CCLE, 
##                 sensitivity.measure='ic50_published', 
##                 summary.stat="median",
##                 verbose=FALSE)
##   
##   GDSCexpression <- summarizeMolecularProfiles(common$GDSC, 
##                                         cellNames(common$GDSC),
##                                         mDataType="rna",
##                                         features=commonGenes,
##                                         verbose=FALSE)
##   CCLEexpression <- summarizeMolecularProfiles(common$CCLE, 
##                                          cellNames(common$CCLE),
##                                          mDataType="rna",
##                                          features=commonGenes,
##                                          verbose=FALSE)
##   gg <- fNames(common[[1]], 'rna')
##   cc <- cellNames(common[[1]])
##   
##   ge.cor <- sapply(cc, function (x, d1, d2) {
##     return (stats::cor(d1[ , x], d2[ , x], method="spearman",
##                 use="pairwise.complete.obs"))
##   }, d1=exprs(GDSCexpression), d2=exprs(CCLEexpression))
##   ic50.cor <- sapply(cc, function (x, d1, d2) {
##     return (stats::cor(d1[, x], d2[ , x], method="spearman",
##                 use="pairwise.complete.obs"))
##   }, d1=GDSC.ic50, d2=CCLE.ic50)
##   auc.cor <- sapply(cc, function (x, d1, d2) {
##     return (stats::cor(d1[ , x], d2[ , x], method="spearman",
##                 use="pairwise.complete.obs"))
##   }, d1=GDSC.auc, d2=CCLE.auc)
##   
##   w1 <- stats::wilcox.test(x=ge.cor, y=auc.cor,
##                            conf.int=TRUE, exact=FALSE)
##   w2 <- stats::wilcox.test(x=ge.cor, y=ic50.cor,
##                            conf.int=TRUE, exact=FALSE)
##   yylim <- c(-1, 1)
##   ss <- sprintf("GE vs. AUC = %.1E\nGE vs. IC50 = %.1E",
##                 w1$p.value, w2$p.value)
##   boxplot(list("GE"=ge.cor,
##                "AUC"=auc.cor,
##                "IC50"=ic50.cor),
##           main="Concordance between cell lines",
##           ylab=expression(R[s]),
##           sub=ss,
##           ylim=yylim,
##           col="lightgrey",
##           pch=20,
##           border="black")
## 


###################################################
### code chunk number 7: fig2
###################################################

  library(Biobase)
  library(PharmacoGx)
  data("GDSCsmall")
  data("CCLEsmall")
  commonGenes <- intersect(fNames(GDSCsmall, "rna"),
                           fNames(CCLEsmall,"rna"))
  common <- intersectPSet(list('CCLE'=CCLEsmall,
                               'GDSC'=GDSCsmall),
                          intersectOn=c("cell.lines", "drugs"), strictIntersect=TRUE)

  
  GDSC.auc <- summarizeSensitivityProfiles(
                pSet=common$GDSC,
                sensitivity.measure='auc_published', 
                summary.stat="median",
                verbose=FALSE)
  CCLE.auc <- summarizeSensitivityProfiles(
                pSet=common$CCLE,
                sensitivity.measure='auc_published', 
                summary.stat="median",
                verbose=FALSE)
  
  GDSC.ic50 <- summarizeSensitivityProfiles(
                pSet=common$GDSC, 
                sensitivity.measure='ic50_published', 
                summary.stat="median",
                verbose=FALSE)
  CCLE.ic50 <- summarizeSensitivityProfiles(
                pSet=common$CCLE, 
                sensitivity.measure='ic50_published', 
                summary.stat="median",
                verbose=FALSE)
  
  GDSCexpression <- summarizeMolecularProfiles(common$GDSC, 
                                        cellNames(common$GDSC),
                                        mDataType="rna",
                                        features=commonGenes,
                                        verbose=FALSE)
  CCLEexpression <- summarizeMolecularProfiles(common$CCLE, 
                                         cellNames(common$CCLE),
                                         mDataType="rna",
                                         features=commonGenes,
                                         verbose=FALSE)
  gg <- fNames(common[[1]], 'rna')
  cc <- cellNames(common[[1]])
  
  ge.cor <- sapply(cc, function (x, d1, d2) {
    return (stats::cor(d1[ , x], d2[ , x], method="spearman",
                use="pairwise.complete.obs"))
  }, d1=exprs(GDSCexpression), d2=exprs(CCLEexpression))
  ic50.cor <- sapply(cc, function (x, d1, d2) {
    return (stats::cor(d1[, x], d2[ , x], method="spearman",
                use="pairwise.complete.obs"))
  }, d1=GDSC.ic50, d2=CCLE.ic50)
  auc.cor <- sapply(cc, function (x, d1, d2) {
    return (stats::cor(d1[ , x], d2[ , x], method="spearman",
                use="pairwise.complete.obs"))
  }, d1=GDSC.auc, d2=CCLE.auc)
  
  w1 <- stats::wilcox.test(x=ge.cor, y=auc.cor,
                           conf.int=TRUE, exact=FALSE)
  w2 <- stats::wilcox.test(x=ge.cor, y=ic50.cor,
                           conf.int=TRUE, exact=FALSE)
  yylim <- c(-1, 1)
  ss <- sprintf("GE vs. AUC = %.1E\nGE vs. IC50 = %.1E",
                w1$p.value, w2$p.value)
  boxplot(list("GE"=ge.cor,
               "AUC"=auc.cor,
               "IC50"=ic50.cor),
          main="Concordance between cell lines",
          ylab=expression(R[s]),
          sub=ss,
          ylim=yylim,
          col="lightgrey",
          pch=20,
          border="black")

    



###################################################
### code chunk number 8: load-cmap
###################################################
  library(PharmacoGx)
  require(xtable)
  data(CMAPsmall)
  drug.perturbation <- drugPerturbationSig(CMAPsmall, 
                                           mDataType="rna",
                                           verbose=FALSE)
  data(HDAC_genes)
  
  res <- apply(drug.perturbation[,,c("tstat", "fdr")],
               2, function(x, HDAC){ 
	    return(connectivityScore(x=x, 
	                             y=HDAC[,2,drop=FALSE], 
	                             method="gsea", nperm=100))
	}, HDAC=HDAC_genes)
  
  rownames(res) <- c("Connectivity", "P Value")
  res <- t(res)
  res <- res[order(res[,1], decreasing=TRUE),]
  xtable(res, 
    caption='Connectivity Score results for HDAC inhibitor gene signature.')


###################################################
### code chunk number 9: biomarkers
###################################################

  data(CCLEsmall)
  features <- fNames(CCLEsmall, "rna")[
                          which(featureInfo(CCLEsmall,
                                            "rna")$Symbol == "NQO1")]
  sig.rna <- drugSensitivitySig(pSet=CCLEsmall, 
                            mDataType="rna", 
                            drugs=c("17-AAG"), 
                            features=features, 
                            sensitivity.measure="auc_published", 
                            molecular.summary.stat="median", 
                            sensitivity.summary.stat="median",
                            verbose=FALSE)
  sig.mut <- drugSensitivitySig(pSet=CCLEsmall, 
                            mDataType="mutation", 
                            drugs=c("PD-0325901"), 
                            features="BRAF", 
                            sensitivity.measure="auc_published", 
                            molecular.summary.stat="and", 
                            sensitivity.summary.stat="median",
                            verbose=FALSE)
  sig <- rbind(sig.rna, sig.mut)
  rownames(sig) <- c("17-AAG + NQO1","PD-0325901 + BRAF")
  colnames(sig) <- dimnames(sig.mut)[[3]]
  xtable(sig, caption='P Value of Gene-Drug Association')


###################################################
### code chunk number 10: sessionInfo
###################################################
toLatex(sessionInfo())


