Nothing
R/gpr_draw.R
In nsgp: Non-Stationary Gaussian Process Regression
Defines functions
.ggplotcolors
.drawShading
.drawSamples
.plotls
plot.gpsimple
plot.gp
plot.gppack
Documented in
plot.gp
plot.gppack
plot.gpsimple
#
# tool functions
#
# returns a set of colors used by ggplot2, they are pretty good
.ggplotcolors <- function(n=6, h=c(0, 360) +15){
if ((diff(h)%%360) < 1) h[2] <- h[2] - 360/n
hcl(h = (seq(h[1], h[2], length = n)), c = 100, l = 65)
}
# draw the variance shading around mean
.drawShading <- function(x, mean, std, color, noisestd=NULL,
minalpha=0.05, maxalpha=0.55, sigma=2, plotgradient=TRUE) {
r <- col2rgb(color)[1]/255
g <- col2rgb(color)[2]/255
bl <- col2rgb(color)[3]/255
pupper95 = mean + sigma*std
plower95 = mean - sigma*std
# draw posterior 95% CI
lines(x=c(x,rev(x)), y=c(plower95,rev(pupper95)), col=rgb(r,g,bl,maxalpha))
# also draw noise curves
if (!is.null(noisestd)) {
npupper95 = mean + sigma*(sqrt(std^2 + noisestd^2))
nplower95 = mean - sigma*(sqrt(std^2 + noisestd^2))
lines(x=c(x,rev(x)), y=c(nplower95,rev(npupper95)), col=rgb(r,g,bl,1), lty=2, lwd=1.3)
}
# draw shading color by variance
if (plotgradient) {
for (i in 1:length(x)) {
if (i == 1) {
xprev <- x[1]
ylprev <- plower95[1]
yuprev <- pupper95[1]
xnext <- mean(c(x[i],x[i+1]))
ylnext <- mean(c(plower95[i], plower95[i+1]))
yunext <- mean(c(pupper95[i], pupper95[i+1]))
} else if (i == length(x)) {
xprev <- mean(c(x[i],x[i-1]))
ylprev <- mean(c(plower95[i], plower95[i-1]))
yuprev <- mean(c(pupper95[i], pupper95[i-1]))
xnext <- x[i]
yunext <- pupper95[i]
ylnext <- plower95[i]
} else {
# mean time point between current and previous
xprev <- mean(c(x[i],x[i-1]))
xnext <- mean(c(x[i],x[i+1]))
ylprev <- mean(c(plower95[i], plower95[i-1]))
ylnext <- mean(c(plower95[i], plower95[i+1]))
yuprev <- mean(c(pupper95[i], pupper95[i-1]))
yunext <- mean(c(pupper95[i], pupper95[i+1]))
}
polygon(x= c(xprev, xnext, xnext, xprev),
y= c(ylprev, ylnext, yunext, yuprev),
density=NA, col= rgb(r,g,bl, max(minalpha,maxalpha-1.1*std[i])), border=NA)
}
} else { # draw simple polygon
polygon(x= c(x,rev(x)),
y= c(mean-2*std,rev(mean)+2*rev(std)),
density=NA, col= rgb(r,g,bl, (minalpha+maxalpha)/2), border=NA)
}
}
#' @import MASS
.drawSamples <- function(timepoints, mean, cov, N=3) {
sample = mvrnorm(N, mean, cov)
for (i in 1:N) {
lines(x= timepoints, y= sample[i,], col = "grey", lw=1)
points(timepoints, sample[i,], col='black', pch=20, cex=0.6)
}
}
# add lengthscale on right y-axis
.plotls = function(obj) {
par(new=T)
plot(obj$targets$x, .bexp(obj$targets$x, obj$params$l, obj$params$lmin, obj$params$c), axes=F, ylab="", ylim=c(0,6),
t='l', lty=1, lwd=1.2, xaxs='i', col='blue', xlab='')
axis(4, at=c(1,3,5), col.axis='blue', mgp=c(3,0.5,0))
mtext('ell(t)', side=4, line=1.5, col='blue', cex=1.2)
}
#' @keywords internal
#' @title Plot a simple gaussian process
#' @description Plots a gaussian process. Several boolean parameters for modifying the
#' plot. By default plots the data, posterior mean and 95\% interval.
#' @param x the gp-object
#' @param y placeholder variable
#' @param plotdata plot the data (default)
#' @param plotmean plot the GP mean (default)
#' @param plotcov plot the GP covariances (default)
#' @param plotnoise plot the observational noise (default)
#' @param samples plot N samples from the GP
#' @param sigma variance level to plot
#' @param title plot title
#' @param legend plot legend
#' @param plotgradient use gradient graphics
#' @param ... ...
plot.gpsimple = function(x, y=NULL, plotdata=TRUE, plotmean=TRUE, plotcov=TRUE, plotnoise=FALSE,
samples=0, sigma=2, title=NULL, legend=FALSE, plotgradient=TRUE,...) {
obj=x
xmin <- min(obj$x)
xmax <- max(obj$x)
if (!hasArg("ylim") & plotnoise == TRUE) {
ymax = max(sigmaupper(obj,T,sigma))
ymin = min(sigmalower(obj,T,sigma))
ylim = c(ymin-0.2,ymax+0.2)
} else if (!hasArg("ylim") & plotnoise == FALSE) {
ymax = max(sigmaupper(obj,F,sigma))
ymin = min(sigmalower(obj,F,sigma))
ylim = c(ymin-0.2,ymax+0.2)
} else {
ylim = list(...)$ylim
}
if (!hasArg("xlab")) { xlab = 'Time' } else { xlab = list(...)$xlab }
if (!hasArg("ylab")) { ylab = 'Value' } else { ylab = list(...)$ylab }
colors <- .ggplotcolors(n=3)[1]
plot(NA, xlim=c(xmin,xmax), ylim=ylim, xlab='', ylab='', xaxs='i', main=title, mgp=c(3,0.7,0))
mtext(xlab, side=1, line=2, cex=1.2)
mtext(ylab, side=2, line=1.8, cex=1.2)
if (plotcov && plotnoise) {
.drawShading(obj$x, obj$mean, sqrt(diag(obj$cov)), colors, obj$noisestd, plotgradient=plotgradient, sigma=sigma)
} else if (plotcov) {
.drawShading(obj$x, obj$mean, sqrt(diag(obj$cov)), colors, plotgradient=plotgradient, sigma=sigma)
}
if (plotmean) { lines(obj$x, obj$mean, col=colors, lwd=3) }
if (plotdata) { points(obj$x.obs, obj$y.obs, bg=colors, pch=21) }
if (samples > 0) { .drawSamples(obj$x, obj$mean, obj$cov, samples) }
if (legend) {
legend('topright', c('Mean', paste('sigma=',sigma,' CI',sep=''),paste('sigma=',sigma,' noisy CI',sep='')), lty=c(1,1,2), lwd=c(3,1,2), col=c('red','red','grey'))
}
}
#' @export
#' @title Plot a gaussian process
#' @description Plots a gaussian process. Several boolean parameters for modifying the
#' plot. By default plots the data, posterior mean and 95\% interval.
#' @param x the gp-object
#' @param y placeholder variable
#' @param plotdata plot the data (default)
#' @param plotmean plot the GP mean (default)
#' @param plotcov plot the GP covariances (default)
#' @param plotnoise plot the observational noise (default)
#' @param samples plot N samples from the GP
#' @param sigma variance level to plot
#' @param title plot title
#' @param legend plot legend
#' @param plotgradient use gradient graphics
#' @param plotls plot lengthscale function
#' @param ... ...
#' @examples
#' # read toy data
#' data(toydata)
#'
#' \dontrun{can take several minutes
#' # perform one-sample regression
#' res = gpr2sample(toydata$ctrl$x, toydata$ctrl$y, seq(0,22,0.1))
#'
#' # pre-computed model for toydata
#' data(toygps)
#' res = toygps$ctrlmodel
#'
#' # basic plot with data, estimated mean and 95\%
#' plot(res)
#'
#' # don't plot the data, plot some samples drawn from the learned gp
#' plot(res, plotdata=FALSE, samples=3)}
plot.gp <- function(x, y=NULL, plotdata=TRUE, plotmean=TRUE, plotcov=TRUE, plotnoise=FALSE,
samples=0, sigma=2, title=NULL, legend=FALSE, plotgradient=TRUE, plotls=FALSE, ...) {
obj=x
xmin <- min(obj$x)
xmax <- max(obj$x)
if (!hasArg("ylim") & plotnoise == TRUE) {
ymax = max(sigmaupper(obj,T,sigma))
ymin = min(sigmalower(obj,T,sigma))
ylim = c(ymin-0.2,ymax+0.2)
} else if (!hasArg("ylim") & plotnoise == FALSE) {
ymax = max(sigmaupper(obj,F,sigma))
ymin = min(sigmalower(obj,F,sigma))
ylim = c(ymin-0.2,ymax+0.2)
} else {
ylim = list(...)$ylim
}
if (!hasArg("xlab")) { xlab = 'Time' } else { xlab = list(...)$xlab }
if (!hasArg("ylab")) { ylab = 'Value' } else { ylab = list(...)$ylab }
colors <- .ggplotcolors(n=3)[1]
plot(NA, xlim=c(xmin,xmax), ylim=ylim, xlab='', ylab='', xaxs='i', main=title, mgp=c(3,0.7,0))
mtext(xlab, side=1, line=2, cex=1.2)
mtext(ylab, side=2, line=1.8, cex=1.2)
if (plotcov && plotnoise) {
.drawShading(obj$targets$x, obj$targets$pmean, obj$targets$pstd, colors, obj$targets$noisestd, plotgradient=plotgradient, sigma=sigma)
} else if (plotcov) {
.drawShading(obj$targets$x, obj$targets$pmean, obj$targets$pstd, colors, plotgradient=plotgradient, sigma=sigma)
}
if (plotmean) { lines(obj$targets$x, obj$targets$pmean, col=colors, lwd=3) }
if (plotdata) { points(obj$x, obj$y, bg=colors, pch=21) }
if (samples > 0) { .drawSamples(obj$targets$x, obj$targets$pmean, obj$cov, samples) }
if (legend) {
legend('topright', c('Mean', paste('sigma=',sigma,' CI',sep=''),paste('sigma=',sigma,' noisy CI',sep='')), lty=c(1,1,2), lwd=c(3,1,2), col=c('red','red','grey'))
}
# draw also the lengthscale function
if (plotls) { .plotls(obj) }
}
#' @export
#' @title Plots several GP's simultaneously
#' @description Plots the GP's corresponding to the control and case data,
#' as well as the null model. Visualizes the log likelihood ratios between
#' the null and individual models. Several boolean parameters for modifying
#' the plot. By default plots the data, posterior mean and 95\% interval for
#' CASE and CONTROL.
#' @details The threshold \code{thr} is the logarithmic likelihood ratio between
#' null and control+case models. The default value 1 hence corresponds to
#' a likelihood ratio of 2.72.
#' @param x the \code{gppack}-object
#' @param y placeholder variable
#' @param plotdata plot the data (default)
#' @param plotmeans plot the GP mean (default)
#' @param plotcovs plot the GP covariances (default)
#' @param plotnoises plot the observational noise (default)
#' @param plotnull plots also the null model
#' @param plotratios plots the ratios, choices are \code{emll}, \code{mll}, \code{npc}, \code{pc}
#' @param thr ratio threshold
#' @param samples plot N samples from the GP
#' @param sigma variance level to plot
#' @param title plot title
#' @param legend plot legend
#' @param plotgradient use gradient graphics
#' @param ... ...
#' @examples
#' # read toy data
#' data(toydata)
#'
#' \dontrun{can take several minutes
#' # perform two-sample regression
#' res = gpr2sample(toydata$ctrl$x, toydata$ctrl$y, toydata$case$x, toydata$case$y, seq(0,22,0.1))
#'
#' # pre-computed model for toydata
#' data(toygps)
#' res = toygps
#'
#' # basic plot
#' plot(res)
#'
#' # plot also the null model, don't plot data, means or noise
#' plot(res, plotnull=TRUE, plotdata=FALSE, plotmeans=FALSE, plotnoise=FALSE)}
plot.gppack <- function(x, y=NULL, plotdata=TRUE, plotmeans=TRUE, plotcovs=TRUE, plotnoises=FALSE,
plotnull=FALSE, plotratios='emll', thr=1,
samples=0, sigma=2, title=NULL, legend=FALSE, plotgradient=TRUE, ...) {
obj=x
x.targets = obj$ctrlmodel$targets$x
K <- length(x.targets) # timepoint count
xmin <- min(x.targets)
xmax <- max(x.targets)
if (!hasArg("xlab")) { xlab = 'Time' } else { xlab = list(...)$xlab }
if (!hasArg("ylab")) { ylab = 'Value' } else { ylab = list(...)$ylab }
colors <- .ggplotcolors(n=3)
colors = colors[c(3,2,1)]
if (plotcovs & plotnoises) {
maxcase = max(sigmaupper(obj$casemodel,T,sigma))
maxctrl = max(sigmaupper(obj$ctrlmodel,T,sigma))
mincase = min(sigmalower(obj$casemodel,T,sigma))
minctrl = min(sigmalower(obj$ctrlmodel,T,sigma))
ymax = max(c(maxctrl,maxcase))
ymin = min(c(minctrl,mincase))
} else if (plotcovs) {
ymax <- max(max(c(sigmaupper(obj$ctrlmodel,F,sigma), sigmaupper(obj$casemodel,F,sigma))))
ymin <- min(min(c(sigmalower(obj$casemodel,F,sigma), sigmalower(obj$casemodel,F,sigma))))
} else {
ymax = max(max(c(obj$ctrlmodel$y, obj$casemodel$y)))
ymin = min(min(c(obj$ctrlmodel$y, obj$casemodel$y)))
}
################### upper
if (plotratios == "emll" || plotratios == 'npc' || plotratios == 'pc' || plotratios == 'mll') {
# Use base-graphics
layout(matrix(c(1,2), 2, 1, byrow=TRUE), heights=c(0.15,0.85), widths=1)
par(mar=c(0,4,1,1), mgp=c(3,0.8,0))
plot(-1, xlim=c(xmin, xmax), ylim=c(0, max(2,2*thr)), ylab="", xaxt="n",
yaxt="n",xaxs="i", yaxs="i", main=title, cex.main=1.2)
mtext(text="log\nratio", side=2, line=1.4, cex=1.4)
axis(2, at=c(0, thr, max(2, 2*thr)), las=2, cex.axis=.9, tck=0.06)
axis(4, tck=0.06, at=c(1, thr, 3), labels=NA)
axis(1, tck=0.06, labels=NA)
axis(3, tck=0.06, labels=NA)
lines(x.targets, obj$ratios[[plotratios]], lwd=1.5)
polygon(c(x.targets,rev(x.targets)),
c(rep(0,K), rev(obj$ratios[[plotratios]])),
density=NA, border=NA, col=colors[3])
if (!is.null(thr)) {
lines(x=x.targets, y=rep(thr,K), col="grey", lwd=1)
}
}
if (!hasArg("ylim")) { ylim = c(ymin-0.2,ymax+0.2) } else { ylim = list(...)$ylim }
################### lower
par(mar=c(5,4,0,1))
plot(NA, xlim=c(xmin,xmax), ylim=ylim, xlab='', ylab='', xaxt="n", yaxt="n", xaxs="i", yaxs="i")
mtext(text=xlab, side=1, line=2.5, cex=1.4)
mtext(text=ylab, side=2, line=2.5, cex=1.4)
axis(2, las=2, cex.axis=.9, tck=0.015)
axis(1, tck=0.015)
axis(4, labels=NA, tck=0.015)
axis(3, tck=0.015, labels=NA)
if (legend==TRUE) {
if (plotnoises && plotnull) {
legend('topright', c('Control posterior', 'Case posterior', 'Shared posterior', ' noisy posterior'), pch=c(20,20,NA,NA), col=c(colors[1],colors[3],colors[2],'grey45'), lty=c(1,1,1,2))
} else if (!plotnoises && plotnull) {
legend('topright', c('Control posterior', 'Case posterior', 'Shared posterior'), pch=c(20,20,NA), col=c(colors[1],colors[3],colors[2]), lty=c(1,1,1))
} else if (plotnoises && !plotnull) {
legend('topright', c('Control posterior', 'Case posterior', ' noisy posterior'), pch=c(20,20,NA), col=c(colors[1],colors[3],'gray45'), lty=c(1,1,2))
} else {
legend('topright', c('Control posterior', 'Case posterior'), pch=c(20,20), col=c(colors[1],colors[3]), lty=c(1,1))
}
}
if (!is.null(thr) && (plotratios == "emll" || plotratios == 'npc' || plotratios == 'pc' || plotratios == 'mll')) {
i <- 1
while(i < K) {
if (obj$ratios[[plotratios]][i]>thr) {
j <- 1
while(obj$ratios[[plotratios]][i+j]>thr && i+j < K-1){
j <- j+1
}
xprev <- mean(c(x.targets[i],x.targets[i-1]))
xnext <- mean(c(x.targets[i+j],x.targets[i+j+1]))
polygon(x=c(xprev, xnext, xnext, xprev),
y=c(-100,-100,100,100),
density=NA, border=NA, col=rgb(0.0,0.0,0.0,0.10))
i <- i+j+1
}
else { i <- i+1 }
}
}
if (plotcovs) {
if (plotnoises) {
.drawShading(x.targets, obj$ctrlmodel$targets$pmean, obj$ctrlmodel$targets$pstd, colors[1], obj$ctrlmodel$targets$noisestd, plotgradient=plotgradient)
.drawShading(x.targets, obj$casemodel$targets$pmean, obj$casemodel$targets$pstd, colors[3], obj$casemodel$targets$noisestd, plotgradient=plotgradient)
if (plotnull) {
.drawShading(x.targets, obj$nullmodel$targets$pmean, obj$nullmodel$targets$pstd, colors[2], obj$nullmodel$targets$noisestd, plotgradient=plotgradient)
}
} else {
.drawShading(x.targets, obj$ctrlmodel$targets$pmean, obj$ctrlmodel$targets$pstd, colors[1], plotgradient=plotgradient)
.drawShading(x.targets, obj$casemodel$targets$pmean, obj$casemodel$targets$pstd, colors[3], plotgradient=plotgradient)
if (plotnull) {
.drawShading(x.targets, obj$nullmodel$targets$pmean, obj$nullmodel$targets$pstd, colors[2], plotgradient=plotgradient)
}
}
}
if (plotdata) {
points(obj$ctrlmodel$x, obj$ctrlmodel$y, bg=colors[1], pch=21, cex=1.4)
points(obj$casemodel$x, obj$casemodel$y, bg=colors[3], pch=21, cex=1.4)
}
if (plotmeans) {
lines(x.targets, obj$ctrlmodel$targets$pmean, col=colors[1], lwd=3.5)
lines(x.targets, obj$casemodel$targets$pmean, col=colors[3], lwd=3.5)
if (plotnull) {
lines(x.targets, obj$nullmodel$targets$pmean, col=colors[2], lwd=3.5)
}
}
}
#
#
#
#
#
#
#
#
#
#
#
# # draw error bars
# .errorBars = function(x, y, noise, sigma=2, bar=0.1) {
# segments(x,y-sigma*noise,x,y+sigma*noise)
# segments(x-bar,y-sigma*noise,x+bar,y-sigma*noise)
# segments(x-bar,y+sigma*noise,x+bar,y+sigma*noise)
# }
#
#
# # Draw a colored line(or area), Red/Green/Black(or transparent) based on (mean +/- 2*std)'s value
# .drawRGB <- function(timepoints, mean, std, ylimits, lwd=3, shift=0){
# # Shift value is used to draw multiple curves in a single plot by shifting it's y value.
# color <- ifelse( (mean+2*std)<(1+shift),'forestgreen',ifelse( (mean-2*std)>(1+shift),'red','black'))
# # Colored lines
# # s <- seq(length(timepoints)-1)
# # segments(timepoints[s], mean[s],timepoints[s+1], mean[s+1], col= color, lwd=lwd)
# # segments(timepoints[s], 1+shift,timepoints[s+1], 1+shift , col= color, lwd=lwd/2)
# # Colored areas
# i <- 1
# while(i < (length(timepoints)-1)){
# actualcol <- color[i]
# j <- 1
# while(i+j <= (length(timepoints)-1) && color[i+j] == actualcol){
# j <- j+1
# }
# if(actualcol!="black"){
# polygon(x=c(timepoints[i], timepoints[i+j], timepoints[i+j], timepoints[i]),
# y=c(ylimits[1]+shift,ylimits[1]+shift,ylimits[2]+shift,ylimits[2]+shift),
# density=NA, border=NA,
# col=rgb(col2rgb(actualcol)[1]/255,col2rgb(actualcol)[2]/255,col2rgb(actualcol)[3]/255,0.25))
# }
# i <- i+j+1
# }
# }
#
# # Adds \n in a text
# .wrap.it <- function(x, len)
# {
# sapply(x, function(y) paste(strwrap(y, len),
# collapse = "\n"),
# USE.NAMES = FALSE)
# }
#
# ############################################################
# #################### PLOTTING FUNCTIONS ####################
# ############################################################
#
#
#
#
# # Create a set of gp plots. Either the whole data in a PDF (several pages) or just the first page in R.
# # The "names" argument require the full csv file since it shows the prot name, id and gene name.
# plot.gpset <- function(x, names=NULL, nrow=6, ncol=6, file=TRUE, filename="gpsetPlot"
# , plotdata=FALSE, plotnoise=FALSE, colored=FALSE, uniformvar=FALSE, lwd=3, A4=TRUE){
# # nrow=max(1,min(50,nrow))
# ncol=max(1,min(50,ncol))
# ext = ".pdf"
# filename = paste0(filename,ext)
# obj = x
# objSize = length(obj)
# if(!colored){
# color <- .ggplotcolors(n=3)[1]
# colorLine <- color
# } else{
# color <- "gray"
# colorLine <- "black"
# }
# txt <- 30 # White area where text is displayed
# ysize <- 50 # Height of the plot
# minSpace <- 0.4 # Minimum space between 2 rows of curves
# if(A4){
# lwd <- 1.5
# ysize <- 30
# # Settings without protNames
# ncol <- 5
# # Settings with protNames
# # ncol <- 4
# # txt <- 70
# # minSpace <- 0.7
# }
# count <- 1
# # Order the data with the biggest plots first
# ord<-c()
# for(i in 1:length(obj)){
# ord<-c(ord,max((1-min(obj[[i]]$targets$plower95)),(max(obj[[i]]$targets$pupper95)-1)))
# }
# obj <- obj[order(ord,decreasing=T)]
# if(!is.null(names)){
# if(nrow(names)==length(ord)){
# geneName <- names[order(ord,decreasing=T),"symbol"]
# protAcc <- names[order(ord,decreasing=T),"uniprot"]
# protName <- names[order(ord,decreasing=T),"prot_name"]
# }
# }
# ord <- sort(ord,decreasing=T)
#
# # Start of the plotting.
# if(file){
# if(A4){
# pdf(filename,8.27,11.69) # A4 size in inches
# } else {
# pdf(filename,ncol*2,ysize)
# }
# }
# k <- 1
# while(k<=objSize){
# par(mar=c(0,4,0,0))
# plot=plot(NA,xlim=c(0,24*ncol+txt*(ncol-1)),ylim=c(0.5,ysize),type="n",axes=FALSE,ann=FALSE)
# height <- ysize
# while(height>=1){
# if(k<=objSize){
# height <- height-max(minSpace,round(ord[k]+0.1,digits=1))
# }
# if(height>=1){
# # cat("ligne ",count,"\n")
# count <- count+1
# for(j in 0:ncol-1){
# if(k<=objSize){
# # cat(k,"\n")
# # Draw axis segments
# # if(!colored){
# segments(0+j*(24+txt),height+1,24+j*(24+txt),height+1,lwd=0.5)
# # }
# segments(0+j*(24+txt),height+1+0.05,0+j*(24+txt),height+1-0.05,lwd=0.5)
# segments(7+j*(24+txt),height+1+0.05,7+j*(24+txt),height+1-0.05,lwd=0.5)
# segments(14+j*(24+txt),height+1+0.05,14+j*(24+txt),height+1-0.05,lwd=0.5)
# segments(21+j*(24+txt),height+1+0.05,21+j*(24+txt),height+1-0.05,lwd=0.5)
#
# # Draw axis names
# # Version without protNames
# # lab <- paste0(protAcc[k],"\n",geneName[k])
# # Version with protNames
# if(!is.null(names)){
# pName <- protName[k]
# if(nchar(toString(pName))>30){
# pName <- paste0(substring(pName,1,15),"\n",substring(pName,16,min(27,nchar(toString(pName)))),"...")
# } else if(nchar(toString(pName))>15) {
# pName <- paste0(substring(pName,1,15),"\n",substring(pName,16,nchar(toString(pName))))
# }
# lab <- paste0(protAcc[k],"-",geneName[k],"\n",pName)
#
# axis(side=2, las=1, at=height+1, pos=j*(24+txt), labels=lab,cex.axis=0.4)
# }
#
# # Draw plotnoise
# if (plotnoise) {
# .drawShading(obj[[k]]$targets$x + j*(24+txt), obj[[k]]$targets$pmean + height
# , obj[[k]]$targets$pstd, color, obj[[k]]$targets$noisestd, uniformvar=uniformvar)
# } else {
# .drawShading(obj[[k]]$targets$x + j*(24+txt), obj[[k]]$targets$pmean + height
# , obj[[k]]$targets$pstd, color, uniformvar=uniformvar)
# }
# # Draw plotmean
# if(colored){
# ymax = max(1.3,max(obj[[k]]$targets$npupper95))
# ymin = min(1/1.3,min(obj[[k]]$targets$nplower95))
# ylim = c(ymin,ymax)
# .drawRGB(obj[[k]]$targets$x + j*(24+txt), obj[[k]]$targets$pmean + height
# , obj[[k]]$targets$pstd,lwd=lwd,ylimits=ylim,shift=height)
# }
# lines(obj[[k]]$targets$x + j*(24+txt), obj[[k]]$targets$pmean + height,
# col=colorLine, lwd=lwd)
# k <- k+1
# } else {
# height <- 0
# k <- objSize+1
# }
#
# }
# if(k<=objSize){
# height <- height-max(minSpace,round(ord[k-ncol]+0.1,digits=1))
# }
# }
# }
# if(!file) { k <- objSize+1 }
# }
# if(file){
# dev.off()
# }
# }
#
#
# # Create a set of gppack plots. Either the whole data in a PDF (several pages) or just the first page in R.
# plot.gppackset <- function(obj, names=NULL, nrow=5, ncol=7, file=TRUE, filename=NULL,
# plotdata=FALSE, plotnoise=FALSE) {
# N = length(obj)
# nx = round(0.8*sqrt(N))
# ny = ceiling(N/nx)
#
# if (!is.null(filename)) {
# pdf(filename, nx*1.5, ny*1.5)
# }
# get( getOption( "device" ) )()
# split.screen(figs = c(nx,ny))
# # par(mfrow=c(nx,ny))
# for (i in 1:N) {
# screen(i)
# plot(obj[[i]])
# }
# close.screen(all=TRUE)
# if (!is.null(filename)) {
# dev.off()
# }
#
# k <- 1
# while(k<=N){
# split.screen(c(nx,ny))
# for(i in 0:(nx-1)){
# for(j in 1:ny){
# if(k<=N){
# screen(i*ny+j)
# plot(obj[[k]])
# if(length(names)==N){
# mtext(names[k],side=1)
# }
# k <- k+1
# }
# }
# }
# close.screen(all=T)
# }
# }
#
#
#
# ###############################################
# #################### QUERY ####################
# ###############################################
#
# # Query the plots and informations of a gene
# query <- function(target, gdata, pdata, ggps, pgps, gGO=NULL, pGO=NULL,colored=F){
# grow = match(target, gdata[,"symbol"])
# prow = match(target, pdata[,"symbol"])
# if(is.na(prow)&&is.na(grow)){
# prow = match(target, pdata[,"uniprot"])
# }
#
# hgncID <- geneName <- refseq <- protAccess <- protName <- ""
# if(!is.na(grow)){
# hgncID <- gdata[grow,"hgnc"]
# geneName <- gdata[grow,"symbol"]
# refseq <- gdata[grow,"refseq"]
# }
# if(!is.na(prow)){
# hgncID <- pdata[prow,"hgnc"]
# geneName <- pdata[prow,"symbol"]
# protAccess <- pdata[prow,"uniprot"]
# protName <- pdata[prow,"prot_name"]
# }
#
# if(geneName==""){
# cat("NO RESULT")
# }else{
# mainSplit <- rbind(c(0,1,0.7,1),c(0,1,0,0.7))
# split.screen(mainSplit)
#
# ########## INFOS
# screen(1)
# par(mar=c(0,0,0,0))
# plot(NA,xlim=c(1,10),ylim=c(4,0.5),axes=F,ann=F)
# txt <- paste0("[GENE] ",geneName)
# text(1,0.8,txt,pos=4,cex=1.2)
# if(!is.na(grow)){
# txt <- paste0(hgncID," - ",refseq)
# } else {
# txt <- hgncID
# }
# text(1,1.3,txt,pos=4)
# if(nchar(toString(protName))>40){
# #sub <- paste("-",substr(toString(protName),41,nchar(toString(protName))))
# #protName <- paste(substr(toString(protName),1,40),"-")
# protName <- paste0(substr(toString(protName),1,40),"-\n-"
# ,substr(toString(protName),41,nchar(toString(protName))))
# } #else { sub <- "" }
# if(!is.na(prow)){
# txt <- paste0("[PROTEIN] ",protName)
# } else {
# txt <- ""
# }
# text(6,0.8,txt,pos=4,cex=1.2)
# text(6,1.3,protAccess,pos=4)
#
# if(!is.na(grow) && !is.null(gGO)){ # GENE GO
# txt <- "Process :"
# text(1,1.7,txt,pos=4,cex=0.85)
# txt <- .wrap.it(toString(gGO[[grow]][["BPterm"]]),170)
# if(nchar(txt)>480){ txt <- paste0(substring(txt,1,480),"...")}
# text(1,2.2,txt,pos=4,cex=0.85)
# txt <- "Function :"
# text(1,2.7,txt,pos=4,cex=0.85)
# txt <- .wrap.it(toString(gGO[[grow]][["MFterm"]]),170)
# if(nchar(txt)>300){ txt <- paste0(substring(txt,1,300),"...")}
# text(1,3.1,txt,pos=4,cex=0.85)
# txt <- "Component :"
# text(1,3.6,txt,pos=4,cex=0.85)
# txt <- .wrap.it(toString(gGO[[grow]][["CCterm"]]),170)
# if(nchar(txt)>300){ txt <- paste0(substring(txt,1,300),"...")}
# text(1,3.9,txt,pos=4,cex=0.85)
# } else if(!is.na(prow) && !is.null(pGO)){ # PROT GO
# txt <- "Process :"
# text(1,1.7,txt,pos=4,cex=0.85)
# txt <- .wrap.it(toString(pGO[[prow]][["BPterm"]]),170)
# if(nchar(txt)>480){ txt <- paste0(substring(txt,1,480),"...")}
# text(1,2.2,txt,pos=4,cex=0.85)
# txt <- "Function :"
# text(1,2.7,txt,pos=4,cex=0.85)
# txt <- .wrap.it(toString(pGO[[prow]][["MFterm"]]),170)
# if(nchar(txt)>300){ txt <- paste0(substring(txt,1,300),"...")}
# text(1,3.05,txt,pos=4,cex=0.85)
# txt <- "Component :"
# text(1,3.6,txt,pos=4,cex=0.85)
# txt <- .wrap.it(toString(pGO[[prow]][["CCterm"]]),170)
# if(nchar(txt)>300){ txt <- paste0(substring(txt,1,300),"...")}
# text(1,3.95,txt,pos=4,cex=0.85)
# }
#
# ########## PLOTS
# screen(2)
# if(!is.na(grow) && !is.na(prow)){
# plotscr <- split.screen(c(1,2))
# screen(plotscr[1])
# plot(ggps[[grow]])
# screen(plotscr[2])
# plot(pgps[[prow]],colored=colored)
# } else if(!is.na(grow)) {
# plot(ggps[[grow]])
# } else if(!is.na(prow)) {
# plot(pgps[[prow]],colored=colored)
# }
# close.screen(all=T)
# }
# }
#
#
#
#
# plotemll = function(res) {
# plot(timepoints, rep(1,175), ylim=c(-3,1.5), t='l')
# lines(timepoints, res$nullmodel$targets$emll, col='green', lty=2)
# lines(timepoints, res$casemodel$targets$emll, col='blue', lty=2)
# lines(timepoints, res$ctrlmodel$targets$emll, col='red', lty=2)
# lines(timepoints, 0.5*res$ctrlmodel$targets$emll+0.5*res$casemodel$targets$emll-res$nullmodel$targets$emll, col='orange', lty=2)
#
# ctrl.f.mean = res$ctrlmodel$targets$pmean - mean(res$ctrlmodel$y)
# ctrl.k.m = res$ctrlmodel$ekernel
# ctrl.f.cov = res$ctrlmodel$cov
# ctrl.k.y = res$ctrlmodel$kernel
# ctrl.noise = res$ctrlmodel$target$noisestd
# ctrl.emll <- -0.5*ctrl.f.mean^2/diag(ctrl.k.m) -0.5*log(diag(ctrl.k.m)) -0.5*log(2*pi)
#
# case.f.mean = res$casemodel$targets$pmean - mean(res$casemodel$y)
# case.k.m = res$casemodel$ekernel
# case.f.cov = res$casemodel$cov
# case.emll <- -0.5*case.f.mean^2/diag(case.k.m) -0.5*log(diag(case.k.m)) -0.5*log(2*pi)
# null.f.mean = res$nullmodel$targets$pmean - mean(res$nullmodel$y)
# null.k.m = res$nullmodel$ekernel
# null.f.cov = res$nullmodel$cov
# null.emll <- -0.5*null.f.mean^2/diag(null.k.m) -0.5*log(diag(null.k.m)) -0.5*log(2*pi)
#
#
# lines(timepoints, ctrl.f.mean, lty=2)
# lines(timepoints, diag(ctrl.k.m))
#
# plot(timepoints, rep(1,175), ylim=c(-2,1), t='l')
# lines(timepoints, null.emll, col='green')
# lines(timepoints, case.emll, col='blue')
# lines(timepoints, ctrl.emll, col='red')
# lines(timepoints, 0.5*case.emll+0.5*ctrl.emll-null.emll, col='orange')
# }
#
Try the nsgp package in your browser
Any scripts or data that you put into this service are public.
nsgp documentation built on May 2, 2019, 9:19 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions .ggplotcolors .drawShading .drawSamples .plotls plot.gpsimple plot.gp plot.gppack
Documented in plot.gp plot.gppack plot.gpsimple
#
# tool functions
#
# returns a set of colors used by ggplot2, they are pretty good
.ggplotcolors <- function(n=6, h=c(0, 360) +15){
if ((diff(h)%%360) < 1) h[2] <- h[2] - 360/n
hcl(h = (seq(h[1], h[2], length = n)), c = 100, l = 65)
}
# draw the variance shading around mean
.drawShading <- function(x, mean, std, color, noisestd=NULL,
minalpha=0.05, maxalpha=0.55, sigma=2, plotgradient=TRUE) {
r <- col2rgb(color)[1]/255
g <- col2rgb(color)[2]/255
bl <- col2rgb(color)[3]/255
pupper95 = mean + sigma*std
plower95 = mean - sigma*std
# draw posterior 95% CI
lines(x=c(x,rev(x)), y=c(plower95,rev(pupper95)), col=rgb(r,g,bl,maxalpha))
# also draw noise curves
if (!is.null(noisestd)) {
npupper95 = mean + sigma*(sqrt(std^2 + noisestd^2))
nplower95 = mean - sigma*(sqrt(std^2 + noisestd^2))
lines(x=c(x,rev(x)), y=c(nplower95,rev(npupper95)), col=rgb(r,g,bl,1), lty=2, lwd=1.3)
}
# draw shading color by variance
if (plotgradient) {
for (i in 1:length(x)) {
if (i == 1) {
xprev <- x[1]
ylprev <- plower95[1]
yuprev <- pupper95[1]
xnext <- mean(c(x[i],x[i+1]))
ylnext <- mean(c(plower95[i], plower95[i+1]))
yunext <- mean(c(pupper95[i], pupper95[i+1]))
} else if (i == length(x)) {
xprev <- mean(c(x[i],x[i-1]))
ylprev <- mean(c(plower95[i], plower95[i-1]))
yuprev <- mean(c(pupper95[i], pupper95[i-1]))
xnext <- x[i]
yunext <- pupper95[i]
ylnext <- plower95[i]
} else {
# mean time point between current and previous
xprev <- mean(c(x[i],x[i-1]))
xnext <- mean(c(x[i],x[i+1]))
ylprev <- mean(c(plower95[i], plower95[i-1]))
ylnext <- mean(c(plower95[i], plower95[i+1]))
yuprev <- mean(c(pupper95[i], pupper95[i-1]))
yunext <- mean(c(pupper95[i], pupper95[i+1]))
}
polygon(x= c(xprev, xnext, xnext, xprev),
y= c(ylprev, ylnext, yunext, yuprev),
density=NA, col= rgb(r,g,bl, max(minalpha,maxalpha-1.1*std[i])), border=NA)
}
} else { # draw simple polygon
polygon(x= c(x,rev(x)),
y= c(mean-2*std,rev(mean)+2*rev(std)),
density=NA, col= rgb(r,g,bl, (minalpha+maxalpha)/2), border=NA)
}
}
#' @import MASS
.drawSamples <- function(timepoints, mean, cov, N=3) {
sample = mvrnorm(N, mean, cov)
for (i in 1:N) {
lines(x= timepoints, y= sample[i,], col = "grey", lw=1)
points(timepoints, sample[i,], col='black', pch=20, cex=0.6)
}
}
# add lengthscale on right y-axis
.plotls = function(obj) {
par(new=T)
plot(obj$targets$x, .bexp(obj$targets$x, obj$params$l, obj$params$lmin, obj$params$c), axes=F, ylab="", ylim=c(0,6),
t='l', lty=1, lwd=1.2, xaxs='i', col='blue', xlab='')
axis(4, at=c(1,3,5), col.axis='blue', mgp=c(3,0.5,0))
mtext('ell(t)', side=4, line=1.5, col='blue', cex=1.2)
}
#' @keywords internal
#' @title Plot a simple gaussian process
#' @description Plots a gaussian process. Several boolean parameters for modifying the
#' plot. By default plots the data, posterior mean and 95\% interval.
#' @param x the gp-object
#' @param y placeholder variable
#' @param plotdata plot the data (default)
#' @param plotmean plot the GP mean (default)
#' @param plotcov plot the GP covariances (default)
#' @param plotnoise plot the observational noise (default)
#' @param samples plot N samples from the GP
#' @param sigma variance level to plot
#' @param title plot title
#' @param legend plot legend
#' @param plotgradient use gradient graphics
#' @param ... ...
plot.gpsimple = function(x, y=NULL, plotdata=TRUE, plotmean=TRUE, plotcov=TRUE, plotnoise=FALSE,
samples=0, sigma=2, title=NULL, legend=FALSE, plotgradient=TRUE,...) {
obj=x
xmin <- min(obj$x)
xmax <- max(obj$x)
if (!hasArg("ylim") & plotnoise == TRUE) {
ymax = max(sigmaupper(obj,T,sigma))
ymin = min(sigmalower(obj,T,sigma))
ylim = c(ymin-0.2,ymax+0.2)
} else if (!hasArg("ylim") & plotnoise == FALSE) {
ymax = max(sigmaupper(obj,F,sigma))
ymin = min(sigmalower(obj,F,sigma))
ylim = c(ymin-0.2,ymax+0.2)
} else {
ylim = list(...)$ylim
}
if (!hasArg("xlab")) { xlab = 'Time' } else { xlab = list(...)$xlab }
if (!hasArg("ylab")) { ylab = 'Value' } else { ylab = list(...)$ylab }
colors <- .ggplotcolors(n=3)[1]
plot(NA, xlim=c(xmin,xmax), ylim=ylim, xlab='', ylab='', xaxs='i', main=title, mgp=c(3,0.7,0))
mtext(xlab, side=1, line=2, cex=1.2)
mtext(ylab, side=2, line=1.8, cex=1.2)
if (plotcov && plotnoise) {
.drawShading(obj$x, obj$mean, sqrt(diag(obj$cov)), colors, obj$noisestd, plotgradient=plotgradient, sigma=sigma)
} else if (plotcov) {
.drawShading(obj$x, obj$mean, sqrt(diag(obj$cov)), colors, plotgradient=plotgradient, sigma=sigma)
}
if (plotmean) { lines(obj$x, obj$mean, col=colors, lwd=3) }
if (plotdata) { points(obj$x.obs, obj$y.obs, bg=colors, pch=21) }
if (samples > 0) { .drawSamples(obj$x, obj$mean, obj$cov, samples) }
if (legend) {
legend('topright', c('Mean', paste('sigma=',sigma,' CI',sep=''),paste('sigma=',sigma,' noisy CI',sep='')), lty=c(1,1,2), lwd=c(3,1,2), col=c('red','red','grey'))
}
}
#' @export
#' @title Plot a gaussian process
#' @description Plots a gaussian process. Several boolean parameters for modifying the
#' plot. By default plots the data, posterior mean and 95\% interval.
#' @param x the gp-object
#' @param y placeholder variable
#' @param plotdata plot the data (default)
#' @param plotmean plot the GP mean (default)
#' @param plotcov plot the GP covariances (default)
#' @param plotnoise plot the observational noise (default)
#' @param samples plot N samples from the GP
#' @param sigma variance level to plot
#' @param title plot title
#' @param legend plot legend
#' @param plotgradient use gradient graphics
#' @param plotls plot lengthscale function
#' @param ... ...
#' @examples
#' # read toy data
#' data(toydata)
#'
#' \dontrun{can take several minutes
#' # perform one-sample regression
#' res = gpr2sample(toydata$ctrl$x, toydata$ctrl$y, seq(0,22,0.1))
#'
#' # pre-computed model for toydata
#' data(toygps)
#' res = toygps$ctrlmodel
#'
#' # basic plot with data, estimated mean and 95\%
#' plot(res)
#'
#' # don't plot the data, plot some samples drawn from the learned gp
#' plot(res, plotdata=FALSE, samples=3)}
plot.gp <- function(x, y=NULL, plotdata=TRUE, plotmean=TRUE, plotcov=TRUE, plotnoise=FALSE,
samples=0, sigma=2, title=NULL, legend=FALSE, plotgradient=TRUE, plotls=FALSE, ...) {
obj=x
xmin <- min(obj$x)
xmax <- max(obj$x)
if (!hasArg("ylim") & plotnoise == TRUE) {
ymax = max(sigmaupper(obj,T,sigma))
ymin = min(sigmalower(obj,T,sigma))
ylim = c(ymin-0.2,ymax+0.2)
} else if (!hasArg("ylim") & plotnoise == FALSE) {
ymax = max(sigmaupper(obj,F,sigma))
ymin = min(sigmalower(obj,F,sigma))
ylim = c(ymin-0.2,ymax+0.2)
} else {
ylim = list(...)$ylim
}
if (!hasArg("xlab")) { xlab = 'Time' } else { xlab = list(...)$xlab }
if (!hasArg("ylab")) { ylab = 'Value' } else { ylab = list(...)$ylab }
colors <- .ggplotcolors(n=3)[1]
plot(NA, xlim=c(xmin,xmax), ylim=ylim, xlab='', ylab='', xaxs='i', main=title, mgp=c(3,0.7,0))
mtext(xlab, side=1, line=2, cex=1.2)
mtext(ylab, side=2, line=1.8, cex=1.2)
if (plotcov && plotnoise) {
.drawShading(obj$targets$x, obj$targets$pmean, obj$targets$pstd, colors, obj$targets$noisestd, plotgradient=plotgradient, sigma=sigma)
} else if (plotcov) {
.drawShading(obj$targets$x, obj$targets$pmean, obj$targets$pstd, colors, plotgradient=plotgradient, sigma=sigma)
}
if (plotmean) { lines(obj$targets$x, obj$targets$pmean, col=colors, lwd=3) }
if (plotdata) { points(obj$x, obj$y, bg=colors, pch=21) }
if (samples > 0) { .drawSamples(obj$targets$x, obj$targets$pmean, obj$cov, samples) }
if (legend) {
legend('topright', c('Mean', paste('sigma=',sigma,' CI',sep=''),paste('sigma=',sigma,' noisy CI',sep='')), lty=c(1,1,2), lwd=c(3,1,2), col=c('red','red','grey'))
}
# draw also the lengthscale function
if (plotls) { .plotls(obj) }
}
#' @export
#' @title Plots several GP's simultaneously
#' @description Plots the GP's corresponding to the control and case data,
#' as well as the null model. Visualizes the log likelihood ratios between
#' the null and individual models. Several boolean parameters for modifying
#' the plot. By default plots the data, posterior mean and 95\% interval for
#' CASE and CONTROL.
#' @details The threshold \code{thr} is the logarithmic likelihood ratio between
#' null and control+case models. The default value 1 hence corresponds to
#' a likelihood ratio of 2.72.
#' @param x the \code{gppack}-object
#' @param y placeholder variable
#' @param plotdata plot the data (default)
#' @param plotmeans plot the GP mean (default)
#' @param plotcovs plot the GP covariances (default)
#' @param plotnoises plot the observational noise (default)
#' @param plotnull plots also the null model
#' @param plotratios plots the ratios, choices are \code{emll}, \code{mll}, \code{npc}, \code{pc}
#' @param thr ratio threshold
#' @param samples plot N samples from the GP
#' @param sigma variance level to plot
#' @param title plot title
#' @param legend plot legend
#' @param plotgradient use gradient graphics
#' @param ... ...
#' @examples
#' # read toy data
#' data(toydata)
#'
#' \dontrun{can take several minutes
#' # perform two-sample regression
#' res = gpr2sample(toydata$ctrl$x, toydata$ctrl$y, toydata$case$x, toydata$case$y, seq(0,22,0.1))
#'
#' # pre-computed model for toydata
#' data(toygps)
#' res = toygps
#'
#' # basic plot
#' plot(res)
#'
#' # plot also the null model, don't plot data, means or noise
#' plot(res, plotnull=TRUE, plotdata=FALSE, plotmeans=FALSE, plotnoise=FALSE)}
plot.gppack <- function(x, y=NULL, plotdata=TRUE, plotmeans=TRUE, plotcovs=TRUE, plotnoises=FALSE,
plotnull=FALSE, plotratios='emll', thr=1,
samples=0, sigma=2, title=NULL, legend=FALSE, plotgradient=TRUE, ...) {
obj=x
x.targets = obj$ctrlmodel$targets$x
K <- length(x.targets) # timepoint count
xmin <- min(x.targets)
xmax <- max(x.targets)
if (!hasArg("xlab")) { xlab = 'Time' } else { xlab = list(...)$xlab }
if (!hasArg("ylab")) { ylab = 'Value' } else { ylab = list(...)$ylab }
colors <- .ggplotcolors(n=3)
colors = colors[c(3,2,1)]
if (plotcovs & plotnoises) {
maxcase = max(sigmaupper(obj$casemodel,T,sigma))
maxctrl = max(sigmaupper(obj$ctrlmodel,T,sigma))
mincase = min(sigmalower(obj$casemodel,T,sigma))
minctrl = min(sigmalower(obj$ctrlmodel,T,sigma))
ymax = max(c(maxctrl,maxcase))
ymin = min(c(minctrl,mincase))
} else if (plotcovs) {
ymax <- max(max(c(sigmaupper(obj$ctrlmodel,F,sigma), sigmaupper(obj$casemodel,F,sigma))))
ymin <- min(min(c(sigmalower(obj$casemodel,F,sigma), sigmalower(obj$casemodel,F,sigma))))
} else {
ymax = max(max(c(obj$ctrlmodel$y, obj$casemodel$y)))
ymin = min(min(c(obj$ctrlmodel$y, obj$casemodel$y)))
}
################### upper
if (plotratios == "emll" || plotratios == 'npc' || plotratios == 'pc' || plotratios == 'mll') {
# Use base-graphics
layout(matrix(c(1,2), 2, 1, byrow=TRUE), heights=c(0.15,0.85), widths=1)
par(mar=c(0,4,1,1), mgp=c(3,0.8,0))
plot(-1, xlim=c(xmin, xmax), ylim=c(0, max(2,2*thr)), ylab="", xaxt="n",
yaxt="n",xaxs="i", yaxs="i", main=title, cex.main=1.2)
mtext(text="log\nratio", side=2, line=1.4, cex=1.4)
axis(2, at=c(0, thr, max(2, 2*thr)), las=2, cex.axis=.9, tck=0.06)
axis(4, tck=0.06, at=c(1, thr, 3), labels=NA)
axis(1, tck=0.06, labels=NA)
axis(3, tck=0.06, labels=NA)
lines(x.targets, obj$ratios[[plotratios]], lwd=1.5)
polygon(c(x.targets,rev(x.targets)),
c(rep(0,K), rev(obj$ratios[[plotratios]])),
density=NA, border=NA, col=colors[3])
if (!is.null(thr)) {
lines(x=x.targets, y=rep(thr,K), col="grey", lwd=1)
}
}
if (!hasArg("ylim")) { ylim = c(ymin-0.2,ymax+0.2) } else { ylim = list(...)$ylim }
################### lower
par(mar=c(5,4,0,1))
plot(NA, xlim=c(xmin,xmax), ylim=ylim, xlab='', ylab='', xaxt="n", yaxt="n", xaxs="i", yaxs="i")
mtext(text=xlab, side=1, line=2.5, cex=1.4)
mtext(text=ylab, side=2, line=2.5, cex=1.4)
axis(2, las=2, cex.axis=.9, tck=0.015)
axis(1, tck=0.015)
axis(4, labels=NA, tck=0.015)
axis(3, tck=0.015, labels=NA)
if (legend==TRUE) {
if (plotnoises && plotnull) {
legend('topright', c('Control posterior', 'Case posterior', 'Shared posterior', ' noisy posterior'), pch=c(20,20,NA,NA), col=c(colors[1],colors[3],colors[2],'grey45'), lty=c(1,1,1,2))
} else if (!plotnoises && plotnull) {
legend('topright', c('Control posterior', 'Case posterior', 'Shared posterior'), pch=c(20,20,NA), col=c(colors[1],colors[3],colors[2]), lty=c(1,1,1))
} else if (plotnoises && !plotnull) {
legend('topright', c('Control posterior', 'Case posterior', ' noisy posterior'), pch=c(20,20,NA), col=c(colors[1],colors[3],'gray45'), lty=c(1,1,2))
} else {
legend('topright', c('Control posterior', 'Case posterior'), pch=c(20,20), col=c(colors[1],colors[3]), lty=c(1,1))
}
}
if (!is.null(thr) && (plotratios == "emll" || plotratios == 'npc' || plotratios == 'pc' || plotratios == 'mll')) {
i <- 1
while(i < K) {
if (obj$ratios[[plotratios]][i]>thr) {
j <- 1
while(obj$ratios[[plotratios]][i+j]>thr && i+j < K-1){
j <- j+1
}
xprev <- mean(c(x.targets[i],x.targets[i-1]))
xnext <- mean(c(x.targets[i+j],x.targets[i+j+1]))
polygon(x=c(xprev, xnext, xnext, xprev),
y=c(-100,-100,100,100),
density=NA, border=NA, col=rgb(0.0,0.0,0.0,0.10))
i <- i+j+1
}
else { i <- i+1 }
}
}
if (plotcovs) {
if (plotnoises) {
.drawShading(x.targets, obj$ctrlmodel$targets$pmean, obj$ctrlmodel$targets$pstd, colors[1], obj$ctrlmodel$targets$noisestd, plotgradient=plotgradient)
.drawShading(x.targets, obj$casemodel$targets$pmean, obj$casemodel$targets$pstd, colors[3], obj$casemodel$targets$noisestd, plotgradient=plotgradient)
if (plotnull) {
.drawShading(x.targets, obj$nullmodel$targets$pmean, obj$nullmodel$targets$pstd, colors[2], obj$nullmodel$targets$noisestd, plotgradient=plotgradient)
}
} else {
.drawShading(x.targets, obj$ctrlmodel$targets$pmean, obj$ctrlmodel$targets$pstd, colors[1], plotgradient=plotgradient)
.drawShading(x.targets, obj$casemodel$targets$pmean, obj$casemodel$targets$pstd, colors[3], plotgradient=plotgradient)
if (plotnull) {
.drawShading(x.targets, obj$nullmodel$targets$pmean, obj$nullmodel$targets$pstd, colors[2], plotgradient=plotgradient)
}
}
}
if (plotdata) {
points(obj$ctrlmodel$x, obj$ctrlmodel$y, bg=colors[1], pch=21, cex=1.4)
points(obj$casemodel$x, obj$casemodel$y, bg=colors[3], pch=21, cex=1.4)
}
if (plotmeans) {
lines(x.targets, obj$ctrlmodel$targets$pmean, col=colors[1], lwd=3.5)
lines(x.targets, obj$casemodel$targets$pmean, col=colors[3], lwd=3.5)
if (plotnull) {
lines(x.targets, obj$nullmodel$targets$pmean, col=colors[2], lwd=3.5)
}
}
}
#
#
#
#
#
#
#
#
#
#
#
# # draw error bars
# .errorBars = function(x, y, noise, sigma=2, bar=0.1) {
# segments(x,y-sigma*noise,x,y+sigma*noise)
# segments(x-bar,y-sigma*noise,x+bar,y-sigma*noise)
# segments(x-bar,y+sigma*noise,x+bar,y+sigma*noise)
# }
#
#
# # Draw a colored line(or area), Red/Green/Black(or transparent) based on (mean +/- 2*std)'s value
# .drawRGB <- function(timepoints, mean, std, ylimits, lwd=3, shift=0){
# # Shift value is used to draw multiple curves in a single plot by shifting it's y value.
# color <- ifelse( (mean+2*std)<(1+shift),'forestgreen',ifelse( (mean-2*std)>(1+shift),'red','black'))
# # Colored lines
# # s <- seq(length(timepoints)-1)
# # segments(timepoints[s], mean[s],timepoints[s+1], mean[s+1], col= color, lwd=lwd)
# # segments(timepoints[s], 1+shift,timepoints[s+1], 1+shift , col= color, lwd=lwd/2)
# # Colored areas
# i <- 1
# while(i < (length(timepoints)-1)){
# actualcol <- color[i]
# j <- 1
# while(i+j <= (length(timepoints)-1) && color[i+j] == actualcol){
# j <- j+1
# }
# if(actualcol!="black"){
# polygon(x=c(timepoints[i], timepoints[i+j], timepoints[i+j], timepoints[i]),
# y=c(ylimits[1]+shift,ylimits[1]+shift,ylimits[2]+shift,ylimits[2]+shift),
# density=NA, border=NA,
# col=rgb(col2rgb(actualcol)[1]/255,col2rgb(actualcol)[2]/255,col2rgb(actualcol)[3]/255,0.25))
# }
# i <- i+j+1
# }
# }
#
# # Adds \n in a text
# .wrap.it <- function(x, len)
# {
# sapply(x, function(y) paste(strwrap(y, len),
# collapse = "\n"),
# USE.NAMES = FALSE)
# }
#
# ############################################################
# #################### PLOTTING FUNCTIONS ####################
# ############################################################
#
#
#
#
# # Create a set of gp plots. Either the whole data in a PDF (several pages) or just the first page in R.
# # The "names" argument require the full csv file since it shows the prot name, id and gene name.
# plot.gpset <- function(x, names=NULL, nrow=6, ncol=6, file=TRUE, filename="gpsetPlot"
# , plotdata=FALSE, plotnoise=FALSE, colored=FALSE, uniformvar=FALSE, lwd=3, A4=TRUE){
# # nrow=max(1,min(50,nrow))
# ncol=max(1,min(50,ncol))
# ext = ".pdf"
# filename = paste0(filename,ext)
# obj = x
# objSize = length(obj)
# if(!colored){
# color <- .ggplotcolors(n=3)[1]
# colorLine <- color
# } else{
# color <- "gray"
# colorLine <- "black"
# }
# txt <- 30 # White area where text is displayed
# ysize <- 50 # Height of the plot
# minSpace <- 0.4 # Minimum space between 2 rows of curves
# if(A4){
# lwd <- 1.5
# ysize <- 30
# # Settings without protNames
# ncol <- 5
# # Settings with protNames
# # ncol <- 4
# # txt <- 70
# # minSpace <- 0.7
# }
# count <- 1
# # Order the data with the biggest plots first
# ord<-c()
# for(i in 1:length(obj)){
# ord<-c(ord,max((1-min(obj[[i]]$targets$plower95)),(max(obj[[i]]$targets$pupper95)-1)))
# }
# obj <- obj[order(ord,decreasing=T)]
# if(!is.null(names)){
# if(nrow(names)==length(ord)){
# geneName <- names[order(ord,decreasing=T),"symbol"]
# protAcc <- names[order(ord,decreasing=T),"uniprot"]
# protName <- names[order(ord,decreasing=T),"prot_name"]
# }
# }
# ord <- sort(ord,decreasing=T)
#
# # Start of the plotting.
# if(file){
# if(A4){
# pdf(filename,8.27,11.69) # A4 size in inches
# } else {
# pdf(filename,ncol*2,ysize)
# }
# }
# k <- 1
# while(k<=objSize){
# par(mar=c(0,4,0,0))
# plot=plot(NA,xlim=c(0,24*ncol+txt*(ncol-1)),ylim=c(0.5,ysize),type="n",axes=FALSE,ann=FALSE)
# height <- ysize
# while(height>=1){
# if(k<=objSize){
# height <- height-max(minSpace,round(ord[k]+0.1,digits=1))
# }
# if(height>=1){
# # cat("ligne ",count,"\n")
# count <- count+1
# for(j in 0:ncol-1){
# if(k<=objSize){
# # cat(k,"\n")
# # Draw axis segments
# # if(!colored){
# segments(0+j*(24+txt),height+1,24+j*(24+txt),height+1,lwd=0.5)
# # }
# segments(0+j*(24+txt),height+1+0.05,0+j*(24+txt),height+1-0.05,lwd=0.5)
# segments(7+j*(24+txt),height+1+0.05,7+j*(24+txt),height+1-0.05,lwd=0.5)
# segments(14+j*(24+txt),height+1+0.05,14+j*(24+txt),height+1-0.05,lwd=0.5)
# segments(21+j*(24+txt),height+1+0.05,21+j*(24+txt),height+1-0.05,lwd=0.5)
#
# # Draw axis names
# # Version without protNames
# # lab <- paste0(protAcc[k],"\n",geneName[k])
# # Version with protNames
# if(!is.null(names)){
# pName <- protName[k]
# if(nchar(toString(pName))>30){
# pName <- paste0(substring(pName,1,15),"\n",substring(pName,16,min(27,nchar(toString(pName)))),"...")
# } else if(nchar(toString(pName))>15) {
# pName <- paste0(substring(pName,1,15),"\n",substring(pName,16,nchar(toString(pName))))
# }
# lab <- paste0(protAcc[k],"-",geneName[k],"\n",pName)
#
# axis(side=2, las=1, at=height+1, pos=j*(24+txt), labels=lab,cex.axis=0.4)
# }
#
# # Draw plotnoise
# if (plotnoise) {
# .drawShading(obj[[k]]$targets$x + j*(24+txt), obj[[k]]$targets$pmean + height
# , obj[[k]]$targets$pstd, color, obj[[k]]$targets$noisestd, uniformvar=uniformvar)
# } else {
# .drawShading(obj[[k]]$targets$x + j*(24+txt), obj[[k]]$targets$pmean + height
# , obj[[k]]$targets$pstd, color, uniformvar=uniformvar)
# }
# # Draw plotmean
# if(colored){
# ymax = max(1.3,max(obj[[k]]$targets$npupper95))
# ymin = min(1/1.3,min(obj[[k]]$targets$nplower95))
# ylim = c(ymin,ymax)
# .drawRGB(obj[[k]]$targets$x + j*(24+txt), obj[[k]]$targets$pmean + height
# , obj[[k]]$targets$pstd,lwd=lwd,ylimits=ylim,shift=height)
# }
# lines(obj[[k]]$targets$x + j*(24+txt), obj[[k]]$targets$pmean + height,
# col=colorLine, lwd=lwd)
# k <- k+1
# } else {
# height <- 0
# k <- objSize+1
# }
#
# }
# if(k<=objSize){
# height <- height-max(minSpace,round(ord[k-ncol]+0.1,digits=1))
# }
# }
# }
# if(!file) { k <- objSize+1 }
# }
# if(file){
# dev.off()
# }
# }
#
#
# # Create a set of gppack plots. Either the whole data in a PDF (several pages) or just the first page in R.
# plot.gppackset <- function(obj, names=NULL, nrow=5, ncol=7, file=TRUE, filename=NULL,
# plotdata=FALSE, plotnoise=FALSE) {
# N = length(obj)
# nx = round(0.8*sqrt(N))
# ny = ceiling(N/nx)
#
# if (!is.null(filename)) {
# pdf(filename, nx*1.5, ny*1.5)
# }
# get( getOption( "device" ) )()
# split.screen(figs = c(nx,ny))
# # par(mfrow=c(nx,ny))
# for (i in 1:N) {
# screen(i)
# plot(obj[[i]])
# }
# close.screen(all=TRUE)
# if (!is.null(filename)) {
# dev.off()
# }
#
# k <- 1
# while(k<=N){
# split.screen(c(nx,ny))
# for(i in 0:(nx-1)){
# for(j in 1:ny){
# if(k<=N){
# screen(i*ny+j)
# plot(obj[[k]])
# if(length(names)==N){
# mtext(names[k],side=1)
# }
# k <- k+1
# }
# }
# }
# close.screen(all=T)
# }
# }
#
#
#
# ###############################################
# #################### QUERY ####################
# ###############################################
#
# # Query the plots and informations of a gene
# query <- function(target, gdata, pdata, ggps, pgps, gGO=NULL, pGO=NULL,colored=F){
# grow = match(target, gdata[,"symbol"])
# prow = match(target, pdata[,"symbol"])
# if(is.na(prow)&&is.na(grow)){
# prow = match(target, pdata[,"uniprot"])
# }
#
# hgncID <- geneName <- refseq <- protAccess <- protName <- ""
# if(!is.na(grow)){
# hgncID <- gdata[grow,"hgnc"]
# geneName <- gdata[grow,"symbol"]
# refseq <- gdata[grow,"refseq"]
# }
# if(!is.na(prow)){
# hgncID <- pdata[prow,"hgnc"]
# geneName <- pdata[prow,"symbol"]
# protAccess <- pdata[prow,"uniprot"]
# protName <- pdata[prow,"prot_name"]
# }
#
# if(geneName==""){
# cat("NO RESULT")
# }else{
# mainSplit <- rbind(c(0,1,0.7,1),c(0,1,0,0.7))
# split.screen(mainSplit)
#
# ########## INFOS
# screen(1)
# par(mar=c(0,0,0,0))
# plot(NA,xlim=c(1,10),ylim=c(4,0.5),axes=F,ann=F)
# txt <- paste0("[GENE] ",geneName)
# text(1,0.8,txt,pos=4,cex=1.2)
# if(!is.na(grow)){
# txt <- paste0(hgncID," - ",refseq)
# } else {
# txt <- hgncID
# }
# text(1,1.3,txt,pos=4)
# if(nchar(toString(protName))>40){
# #sub <- paste("-",substr(toString(protName),41,nchar(toString(protName))))
# #protName <- paste(substr(toString(protName),1,40),"-")
# protName <- paste0(substr(toString(protName),1,40),"-\n-"
# ,substr(toString(protName),41,nchar(toString(protName))))
# } #else { sub <- "" }
# if(!is.na(prow)){
# txt <- paste0("[PROTEIN] ",protName)
# } else {
# txt <- ""
# }
# text(6,0.8,txt,pos=4,cex=1.2)
# text(6,1.3,protAccess,pos=4)
#
# if(!is.na(grow) && !is.null(gGO)){ # GENE GO
# txt <- "Process :"
# text(1,1.7,txt,pos=4,cex=0.85)
# txt <- .wrap.it(toString(gGO[[grow]][["BPterm"]]),170)
# if(nchar(txt)>480){ txt <- paste0(substring(txt,1,480),"...")}
# text(1,2.2,txt,pos=4,cex=0.85)
# txt <- "Function :"
# text(1,2.7,txt,pos=4,cex=0.85)
# txt <- .wrap.it(toString(gGO[[grow]][["MFterm"]]),170)
# if(nchar(txt)>300){ txt <- paste0(substring(txt,1,300),"...")}
# text(1,3.1,txt,pos=4,cex=0.85)
# txt <- "Component :"
# text(1,3.6,txt,pos=4,cex=0.85)
# txt <- .wrap.it(toString(gGO[[grow]][["CCterm"]]),170)
# if(nchar(txt)>300){ txt <- paste0(substring(txt,1,300),"...")}
# text(1,3.9,txt,pos=4,cex=0.85)
# } else if(!is.na(prow) && !is.null(pGO)){ # PROT GO
# txt <- "Process :"
# text(1,1.7,txt,pos=4,cex=0.85)
# txt <- .wrap.it(toString(pGO[[prow]][["BPterm"]]),170)
# if(nchar(txt)>480){ txt <- paste0(substring(txt,1,480),"...")}
# text(1,2.2,txt,pos=4,cex=0.85)
# txt <- "Function :"
# text(1,2.7,txt,pos=4,cex=0.85)
# txt <- .wrap.it(toString(pGO[[prow]][["MFterm"]]),170)
# if(nchar(txt)>300){ txt <- paste0(substring(txt,1,300),"...")}
# text(1,3.05,txt,pos=4,cex=0.85)
# txt <- "Component :"
# text(1,3.6,txt,pos=4,cex=0.85)
# txt <- .wrap.it(toString(pGO[[prow]][["CCterm"]]),170)
# if(nchar(txt)>300){ txt <- paste0(substring(txt,1,300),"...")}
# text(1,3.95,txt,pos=4,cex=0.85)
# }
#
# ########## PLOTS
# screen(2)
# if(!is.na(grow) && !is.na(prow)){
# plotscr <- split.screen(c(1,2))
# screen(plotscr[1])
# plot(ggps[[grow]])
# screen(plotscr[2])
# plot(pgps[[prow]],colored=colored)
# } else if(!is.na(grow)) {
# plot(ggps[[grow]])
# } else if(!is.na(prow)) {
# plot(pgps[[prow]],colored=colored)
# }
# close.screen(all=T)
# }
# }
#
#
#
#
# plotemll = function(res) {
# plot(timepoints, rep(1,175), ylim=c(-3,1.5), t='l')
# lines(timepoints, res$nullmodel$targets$emll, col='green', lty=2)
# lines(timepoints, res$casemodel$targets$emll, col='blue', lty=2)
# lines(timepoints, res$ctrlmodel$targets$emll, col='red', lty=2)
# lines(timepoints, 0.5*res$ctrlmodel$targets$emll+0.5*res$casemodel$targets$emll-res$nullmodel$targets$emll, col='orange', lty=2)
#
# ctrl.f.mean = res$ctrlmodel$targets$pmean - mean(res$ctrlmodel$y)
# ctrl.k.m = res$ctrlmodel$ekernel
# ctrl.f.cov = res$ctrlmodel$cov
# ctrl.k.y = res$ctrlmodel$kernel
# ctrl.noise = res$ctrlmodel$target$noisestd
# ctrl.emll <- -0.5*ctrl.f.mean^2/diag(ctrl.k.m) -0.5*log(diag(ctrl.k.m)) -0.5*log(2*pi)
#
# case.f.mean = res$casemodel$targets$pmean - mean(res$casemodel$y)
# case.k.m = res$casemodel$ekernel
# case.f.cov = res$casemodel$cov
# case.emll <- -0.5*case.f.mean^2/diag(case.k.m) -0.5*log(diag(case.k.m)) -0.5*log(2*pi)
# null.f.mean = res$nullmodel$targets$pmean - mean(res$nullmodel$y)
# null.k.m = res$nullmodel$ekernel
# null.f.cov = res$nullmodel$cov
# null.emll <- -0.5*null.f.mean^2/diag(null.k.m) -0.5*log(diag(null.k.m)) -0.5*log(2*pi)
#
#
# lines(timepoints, ctrl.f.mean, lty=2)
# lines(timepoints, diag(ctrl.k.m))
#
# plot(timepoints, rep(1,175), ylim=c(-2,1), t='l')
# lines(timepoints, null.emll, col='green')
# lines(timepoints, case.emll, col='blue')
# lines(timepoints, ctrl.emll, col='red')
# lines(timepoints, 0.5*case.emll+0.5*ctrl.emll-null.emll, col='orange')
# }
#
Try the nsgp package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.