#The goal of this function is to calculate the LOB/LOD of the data provided in the data frame.
#The function returns a new data frame containing the value of the LOB/LOD
#' @export
plot_quantlim <- function(spikeindata, quantlim_out, alpha, dir_output, xlim_plot){
if(is.null(quantlim_out)){
print("Assay fit was incorrectly calculated by linear_quantlim or nonlinear_quantlim and cannot be plotted")
return(NULL)
}
#percentile of the prediction interval considered
if(missing(alpha)){
alpha = 5/100;
}
if( alpha >= 1 | alpha <= 0){
print("incorrect specified value for alpha, 0 < alpha < 1")
return(NULL)
}
expdata = spikeindata
datain = quantlim_out
#Define some colors here for the plots:
black1 <- '#000000'; orange1 <- "#E69F00"; blue1 <- "#56B4E9"; green1 <- "#009E73"; yellow1 <- "#F0E442"; blue2 <- "#0072B2"; red1 <- "#D55E00"; pink1 <- "#CC79A7";
cbbPalette <- c(black1, orange1, blue1, green1, yellow1, blue2 ,red1 , pink1)
#Rename variables for the function:
names(datain)[names(datain) == 'CONCENTRATION'] <- 'C'
names(expdata)[names(expdata) == 'CONCENTRATION'] <- 'C'
names(expdata)[names(expdata) == 'INTENSITY'] <- 'I'
names(datain)[names(datain) == 'MEAN'] <- 'mean'
names(datain)[names(datain) == 'LOW'] <- 'low'
names(datain)[names(datain) == 'UP'] <- 'up'
#Remove NA and infinite numbers from spike in data:
expdata <- expdata[!is.na(expdata$I) & !is.na(expdata$C),]
expdata <- expdata[!is.infinite(expdata$I) & !is.infinite(expdata$C),]
#Extract actual data points for plotting:
Cdata = expdata$C
Idata = expdata$I
tmp_blank <- expdata[expdata$C == 0,]
n_blank = length(unique(tmp_blank$I))
noise = mean(tmp_blank$I)
var_noise = var(tmp_blank$I)
## BLANK VARIANCE ZERO PATCH pt1
unique_c = sort(unique(expdata$C)); # list unique concentration points
iii = 1
for (j in unique_c){
data_f <- subset(expdata, C == j)
## patch to account for when the blank has values of 0
if (min(data_f$I) > 0){
nonzeroC <- iii # save the first nonzero concentration point
n_blank <- length(unique(data_f$I)) # set the number of blanks to this
noise <- mean(data_f$I) # set the noise to the mean of this nonzero point
var_noise <- var(data_f$I)
var_blank <- var_noise
tmp_blank <- data_f
break
}
iii = iii +1
}
# PATCH pt2 to account for situations where the blank is all quantitative values of 0
# IF VARIANCE ZERO, ITERATE THROUGH INCREASING CONCENTRATION UNTIL VARIANCE NOT ZERO
#if(var_noise <= 0){
# print("Variance of noise is <=0! Attempting to find nonzero variance in upper concentration levels.")
# var_blank = var_v[nonzeroC]
# var_noise = var_blank
# print(paste("Found! Blank variance set to", var_blank))
#}
## NEED TO FIX THIS TO REFLECT MODIFIED nonlinear_quantlim
fac = qt(1-alpha,n_blank - 1)*sqrt(1+1/n_blank)
#upper bound of noise prediction interval
up_noise = noise + fac* sqrt(var_noise)
rel_size = 2.5; rel_size_2 = 1.8; lw = 1; pw = 2.5;
xaxis_orig_2 = datain$C
tmp_all = datain
LOQ_pred = datain$LOD[1]
LOD_pred = datain$LOB[1]
lower_Q_pred = datain$low
upper_Q_pred = datain$up
mean_bilinear <- datain$mean
if(LOD_pred >= 0){
y_LOD_pred = up_noise
}
if(LOQ_pred >= 0){
y_LOQ_pred = up_noise
}
filename = paste(dir_output,'/',datain$NAME[1],'_',datain$METHOD[1],'_overall.svg',sep='')
svg(filename)
if(LOQ_pred > xaxis_orig_2[3]){
C_max = xaxis_orig_2[min(which(abs(LOQ_pred - xaxis_orig_2) == min(abs(LOQ_pred - xaxis_orig_2))) +1, length(xaxis_orig_2))] }else{
C_max = xaxis_orig_2[which(abs(mean(xaxis_orig_2) - xaxis_orig_2) == min(abs(mean(xaxis_orig_2) - xaxis_orig_2)))]
}
low_p <- paste(alpha*100,"%",sep = "")
high_p <- paste(100-alpha*100,"%", sep = "")
upp_noise <- paste(high_p," upper bound of noise")
low_pred <- paste(low_p , 'percentile of predictions')
p1 <- ggplot() + .theme_complete_bw()
p1 <- p1 + geom_point(data=data.frame(Cdata,Idata) , aes(x=Cdata,y=Idata),size =pw*1.5)
p1 <- p1 + geom_line(data=data.frame(x=xaxis_orig_2,y=mean_bilinear,col='mean prediction', lt = 'mean'),aes(x=x,y=y, color=col), size = lw) #, color = 'black' #as.factor(c('mean bootstrap'))
p1 <- p1 + geom_ribbon(data=data.frame(x=xaxis_orig_2,ymin =lower_Q_pred , ymax =upper_Q_pred), aes(x=x, ymin=ymin, ymax=ymax), fill = red1, alpha = 0.3)
p1 <- p1 + geom_line(data = data.frame(x=xaxis_orig_2,ymin =lower_Q_pred, col = low_pred, lt = 'Int') , aes(x=x,y=ymin, color = col), size = lw)
p1 <- p1 + geom_line(data=data.frame(x=xaxis_orig_2, ymax = rep(up_noise,length(xaxis_orig_2)),
col = "95% upper bound of noise"), aes(x=x, y=ymax, color = col), size = lw)
p1 <- p1 + scale_alpha_continuous(guide = 'none')
p1 <- p1 + xlab('Protein copies on column') + ylab('Integrated peak area') # LKP - CHANGED AXIS LABELS
p1 <- p1 + theme(axis.text.x = element_text(size = rel(rel_size))) + theme(axis.text.y = element_text(size = rel(rel_size)))
p1 <- p1 + theme(axis.title.x = element_text(size = rel(rel_size))) + theme(axis.title.y = element_text(size = rel(rel_size)))
p1 <- p1 + theme(axis.title.y=element_text(vjust=0.7))
p1 <- p1 + scale_color_manual(values = c(orange1, blue1, red1),
labels = c(low_pred,
upp_noise,
"mean prediction"))
#p1 <- p1 + scale_colour_manual(values=c('mean prediction'= red1,'upper 95% prediction' = orange1, '5% percentile of predictions' = orange1, '95% upper bound of noise' = blue1)) + theme(legend.title = element_blank())
#colour_scales <- setNames(c('mean prediction','upper 95% prediction','5% percentile of predictions'),c('dasssta','messssan',"rrrr"))
#p1 <- p1 + scale_colour_manual(values = colour_scales)
p1 <- p1 + theme(legend.title = element_blank()) + theme(legend.position = c(0.05, 0.5), legend.justification = c(0, 0), legend.text=element_text(size=rel(rel_size_2)))
LOD_y = mean_bilinear[which(abs(xaxis_orig_2 - LOD_pred) == min(abs(xaxis_orig_2 - LOD_pred)) )]
p1 <- p1 + geom_point(data=data.frame(x= LOD_pred, y= y_LOD_pred, shape='LOD'),aes(x=x,y=y,shape = shape, guide=FALSE),colour="purple",size=5)
LOQ_y = lower_Q_pred[which(abs(up_noise - lower_Q_pred) == min(abs(up_noise - lower_Q_pred)) )]
p1 <- p1 + geom_point(data=data.frame(x= LOQ_pred, y= y_LOQ_pred, shape='LOQ'),aes(x=x,y=y,shape = shape, guide=FALSE),colour=orange1,size=5)
LOD_string <- paste('LOB=',round(LOD_pred, digits=1),sep=''); LOQ_string <- paste('LOD=',round(LOQ_pred, digits=1),sep='');
p1 <- p1 + guides(colour = guide_legend(order = 1), linetype = guide_legend(order = 1), shape = guide_legend(order = 2)) + guides(shape=FALSE)
p1 <- p1 + ggtitle(paste(datain$NAME,'\n', LOD_string,', ',LOQ_string,sep="")) + theme(plot.title = element_text(size = 20))
print(p1);
dev.off()
#produce a second plot showing a zoomed view:
if(missing(xlim_plot)){ #missing argument for the x limit in the function, pick a x limit that is close to the LOD/LOQ:
if(LOQ_pred > 0){
xlim = LOQ_pred*3.0
} else{
xlim = unique(Cdata)[4]
}
} else{
xlim = xlim_plot
}
filename = paste(dir_output,'/',datain$NAME[1],'_',datain$METHOD[1],'_zoom.svg',sep='')
svg(filename)
Idata = subset(Idata, Cdata < xlim)
Cdata = subset(Cdata, Cdata < xlim)
lower_Q_pred = subset(lower_Q_pred, xaxis_orig_2 < xlim)
upper_Q_pred = subset(upper_Q_pred, xaxis_orig_2 < xlim)
mean_bilinear = subset(mean_bilinear, xaxis_orig_2 < xlim)
xaxis_orig_2 = subset(xaxis_orig_2, xaxis_orig_2 < xlim)
p1 <- ggplot() + .theme_complete_bw()
p1 <- p1 + geom_point(data=data.frame(Cdata,Idata) , aes(x=Cdata,y=Idata),size =pw*1.5)
p1 <- p1 + geom_line(data=data.frame(x=xaxis_orig_2,y=mean_bilinear,col='mean prediction', lt = 'mean'),aes(x=x,y=y, color=col), size = lw) #, color = 'black' #as.factor(c('mean bootstrap'))
p1 <- p1 + geom_ribbon(data=data.frame(x=xaxis_orig_2,ymin =lower_Q_pred , ymax =upper_Q_pred), aes(x=x, ymin=ymin, ymax=ymax), fill = red1, alpha = 0.3)
p1 <- p1 + geom_line(data = data.frame(x=xaxis_orig_2,ymin =lower_Q_pred, col = 'lower 95% prediction', lt = 'Int') , aes(x=x,y=ymin, color = col), size = lw)
p1 <- p1 + geom_line(data=data.frame(x=xaxis_orig_2, ymax = rep(up_noise,length(xaxis_orig_2)),
col = "95% upper bound of noise"), aes(x=x, y=ymax, color = col), size = lw)
p1 <- p1 + scale_alpha_continuous(guide = 'none')
p1 <- p1 + xlab('Protein copies on column') + ylab('Integrated peak area') # LKP - CHANGED AXIS LABELS
p1 <- p1 + theme(axis.text.x = element_text(size = rel(rel_size))) + theme(axis.text.y = element_text(size = rel(rel_size)))
p1 <- p1 + theme(axis.title.x = element_text(size = rel(rel_size))) + theme(axis.title.y = element_text(size = rel(rel_size)))
p1 <- p1 + theme(axis.title.y=element_text(vjust=0.7))
p1 <- p1 + scale_color_manual(values = c(blue1,orange1, red1),
labels = c(upp_noise,low_pred,
"mean prediction"))
#p1 <- p1 + scale_colour_manual(values=c('mean prediction'= red1,'upper 95% prediction' = orange1, 'lower 95% prediction' = orange1, '95% upper bound of noise' = blue1)) + theme(legend.title = element_blank())
#p1 <- p1 + scale_linetype_manual(values=c('dashed','dashed','solid','solid'))
p1 <- p1 + theme(legend.title = element_blank()) + theme(legend.position = c(0.05, 0.5), legend.justification = c(0, 0), legend.text=element_text(size=rel(rel_size_2)))
LOD_y = mean_bilinear[which(abs(xaxis_orig_2 - LOD_pred) == min(abs(xaxis_orig_2 - LOD_pred)) )]
p1 <- p1 + geom_point(data=data.frame(x= LOD_pred, y= y_LOD_pred, shape='LOD'),aes(x=x,y=y,shape = shape, guide=FALSE),colour="purple",size=5)
LOQ_y = lower_Q_pred[which(abs(up_noise - lower_Q_pred) == min(abs(up_noise - lower_Q_pred)) )]
p1 <- p1 + geom_point(data=data.frame(x= LOQ_pred, y= y_LOQ_pred, shape='LOQ'),aes(x=x,y=y,shape = shape, guide=FALSE),colour=orange1,size=5)
LOD_string <- paste('LOB=',round(LOD_pred, digits=1),sep=''); LOQ_string <- paste('LOD=',round(LOQ_pred, digits=1),sep='');
p1 <- p1 + guides(colour = guide_legend(order = 1), linetype = guide_legend(order = 1), shape = guide_legend(order = 2)) + guides(shape=FALSE)
p1 <- p1 + ggtitle(paste(datain$NAME,'\n', LOD_string,', ',LOQ_string,sep="")) + theme(plot.title = element_text(size = 20))
print(p1);
dev.off()
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.