gprege_quick.R
In gprege: Gaussian Process Ranking and Estimation of Gene Expression time-series

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

###################################################
### code chunk number 1: gprege_quick.Rnw:31-32
###################################################
  options(width = 60)


###################################################
### code chunk number 2: gprege_quick.Rnw:47-50 (eval = FALSE)
###################################################
##   if (!requireNamespace("BiocManager", quietly=TRUE))
##       install.packages("BiocManager")
##   BiocManager::install("gprege")


###################################################
### code chunk number 3: gprege_quick.Rnw:57-58
###################################################
  library(gprege)


###################################################
### code chunk number 4: gprege_quick.Rnw:63-64 (eval = FALSE)
###################################################
##   data(FragmentDellaGattaData)


###################################################
### code chunk number 5: gprege_quick.Rnw:69-73 (eval = FALSE)
###################################################
##   # Download full Della Gatta dataset.
##   load(url('https://github.com/alkalait/gprege-datasets/raw/master/R/DellaGattaData.RData'))
##   # OR
##   data(FullDellaGattaData)


###################################################
### code chunk number 6: gprege_quick.Rnw:81-92
###################################################
  # Download full Della Gatta dataset.
  ## con <- url('https://github.com/alkalait/gprege-datasets/raw/master/R/DellaGattaData.RData')
  ## attempts = 0
  ## while(!exists('DGdata') && attempts < 3) {
  ##  try(load(con),TRUE) #  close.connection(con)
  ##  attempts = attempts + 1
  ## }
  data(FullDellaGattaData)
  # Timepoints / GP inputs.
  tTrue = matrix(seq(0,240,by=20), ncol=1)
  gpregeOptions <- list()


###################################################
### code chunk number 7: gprege_quick.Rnw:95-98 (eval = FALSE)
###################################################
##   data(FullDellaGattaData)
##   # Set index range so that only a top few targets suggested by TSNI be explored.
##   gpregeOptions$indexRange <- which(DGatta_labels_byTSNItop100)[1:2]


###################################################
### code chunk number 8: gprege_quick.Rnw:101-106
###################################################
  # Load installed fragment-dataset.
  data(FragmentDellaGattaData)
  # Fragment dataset is comprised of top-ranked targets suggested by TSNI.
  # Explore only the first few.
  gpregeOptions$indexRange <- 1:2


###################################################
### code chunk number 9: gprege_quick.Rnw:109-130
###################################################
  # Explore individual profiles in interactive mode.
  gpregeOptions$explore <- TRUE
  # Exhaustive plot resolution of the LML function.
  gpregeOptions$exhaustPlotRes <- 30
  # Exhaustive plot contour levels.
  gpregeOptions$exhaustPlotLevels <- 10
  # Exhaustive plot maximum lengthscale.
  gpregeOptions$exhaustPlotMaxWidth <- 100
  # Noisy ground truth labels: which genes are in the top 786 ranks of the TSNI ranking.
  gpregeOptions$labels <- DGatta_labels_byTSNItop100 
  # SCG optimisation: maximum number of iterations.
  gpregeOptions$iters <- 100
  # SCG optimisation: no messages.
  gpregeOptions$display <- FALSE
  # Matrix of different hyperparameter configurations as rows:
  # [inverse-lengthscale   percent-signal-variance   percent-noise-variance].
  gpregeOptions$inithypers <-
    matrix( c(  1/1000,	1e-3,	0.999,
                1/30,	0.999,	1e-3,
                1/80,	2/3, 1/3
  ), ncol=3, byrow=TRUE)


###################################################
### code chunk number 10: gprege_quick.Rnw:133-134
###################################################
  gpregeOutput<-gprege(data=exprs_tp63_RMA,inputs=tTrue,gpregeOptions=gpregeOptions)


###################################################
### code chunk number 11: gprege_quick.Rnw:136-156
###################################################
  for (i in 1:length(gpregeOptions$indexRange)) {
    cat("\\begin{figure}[H]")
    cat("\\centering")
#     cat("\\subfigure[]{")
    cat("\\includegraphics[width=0.7\\linewidth]{",  paste("gpPlot",i,".pdf",sep=""), "}\n\n",  sep="")
#     cat("\\label{fig:seqfig", i, "}", sep = "")
    cat("\\caption{GP fit with different initialisations on profile \\#", gpregeOptions$indexRange[i], ".}", sep="")
    cat("\\label{fig:gpPlot", i, "}", sep = "")
    cat("\\end{figure}")
    
    cat("\\begin{figure}[H]")
    cat("\\centering")
    cat("\\subfigure[]{")
    cat("\\includegraphics[page=1,width=.7\\linewidth]{",paste("exhaustivePlot",i,".pdf",sep=""),"}}\n\n", sep="")
    cat("\\subfigure[]{")
    cat("\\includegraphics[page=2,width=.7\\linewidth]{",paste("exhaustivePlot",i,".pdf",sep=""),"}}\n\n", sep="")
    cat("\\caption{Profile \\#", gpregeOptions$indexRange[i], " : (a) Log-marginal likelihood (LML) contour. (b) GP fit with maximum LML hyperparameters from the exhaustive search.}", sep = "")
    cat("\\label{fig:exPlot", i, "}", sep = "")
    cat("\\end{figure}")
  }


###################################################
### code chunk number 12: gprege_quick.Rnw:163-185
###################################################
  # Load fragment dataset.
  data(FragmentDellaGattaData)
  data(DGdat_p63)
  # Download BATS rankings (Angelini, 2007)
  # Case 1: Delta error prior, case 2: Inverse Gamma error prior,
  #  case 3: Double Exponential error prior.
  BATSranking = matrix(0, length(DGatta_labels_byTSNItop100), 3)
  ##for (i in 1:3){
  ##  # Read gene numbers.
  ##  tmp = NULL
  ##  con <- url(paste('https://github.com/alkalait/gprege-datasets/raw/master/R/DGdat_p63_case',i,'_GL.txt',sep=''))
  ##  while(is.null(tmp)) try(tmp <- read.table(con, skip=1), TRUE)
  ##  genenumbers <- as.numeric(lapply( as.character(tmp[,2]), function(x) x=substr(x,2,nchar(x))))
  ##  # Sort rankqings by gene numbers.
  ##  BATSranking[,i] <- tmp[sort(genenumbers, index.return=TRUE)$ix, 4]
  ##}
  genenumbers <- as.numeric(lapply(as.character(DGdat_p63_case1_GL[,2]), function(x) x=substr(x,2,nchar(x))))
  BATSranking[,1] <- DGdat_p63_case1_GL[sort(genenumbers, index.return=TRUE)$ix, 4]
  genenumbers <- as.numeric(lapply(as.character(DGdat_p63_case2_GL[,2]), function(x) x=substr(x,2,nchar(x))))
  BATSranking[,2] <- DGdat_p63_case2_GL[sort(genenumbers, index.return=TRUE)$ix, 4]
  genenumbers <- as.numeric(lapply(as.character(DGdat_p63_case3_GL[,2]), function(x) x=substr(x,2,nchar(x))))
  BATSranking[,3] <- DGdat_p63_case3_GL[sort(genenumbers, index.return=TRUE)$ix, 4]


###################################################
### code chunk number 13: gprege_quick.Rnw:189-196
###################################################
  # The smaller a BATS-rank metric is, the better the rank of the gene
  # reporter. Invert those rank metrics to compare on a common ground
  # with gprege.
  BATSranking = 1/BATSranking
  compareROC(output=gpregeOutput$rankingScores,
          groundTruthLabels=DGatta_labels_byTSNItop100,
          compareToRanking=BATSranking)


###################################################
### code chunk number 14: gprege_quick.Rnw:208-237 (eval = FALSE)
###################################################
##   # Download full Della Gatta dataset.
##   #con <- url('https://github.com/alkalait/gprege-datasets/raw/master/R/DellaGattaData.RData')
##   #while(!exists('DGdata')) try(load(con),TRUE); close.connection(con)
##   data(FullDellaGattaData)
##   # Timepoints / GP inputs.
##   tTrue = matrix(seq(0,240,by=20), ncol=1)
##   gpregeOptions <- list()
##   # Explore individual profiles in interactive mode.
##   gpregeOptions$explore <- FALSE
##   # Exhaustive plot resolution of the LML function.
##   gpregeOptions$exhaustPlotRes <- 30
##   # Exhaustive plot contour levels.
##   gpregeOptions$exhaustPlotLevels <- 10
##   # Exhaustive plot maximum lengthscale.
##   gpregeOptions$exhaustPlotMaxWidth <- 100
##   # Noisy ground truth labels: which genes are in the top 786 ranks of the TSNI ranking.
##   gpregeOptions$labels <- DGatta_labels_byTSNItop100 
##   # SCG optimisation: maximum number of iterations.
##   gpregeOptions$iters <- 100
##   # SCG optimisation: no messages.
##   gpregeOptions$display <- FALSE
##   # Matrix of different hyperparameter configurations as rows:
##   # [inverse-lengthscale   percent-signal-variance   percent-noise-variance].
##   gpregeOptions$inithypers <-
##     matrix( c(  1/1000,  1e-3,	0.999,
##                 1/8,	0.999,	1e-3,
##                 1/80,	2/3, 1/3
##   ), ncol=3, byrow=TRUE)
##   gpregeOutput<-gprege(data=exprs_tp63_RMA,inputs=tTrue,gpregeOptions=gpregeOptions)