#' ms_isothermal_line_innerplot
#'
#' Internal function to generate ggplot objects for isothermal data in a curve plot format
#'
#' @param data isothermal dataset to plot
#' @param legenddata dataset used for condition extraction, at least one protein
#' inside this dataset should contains the full set of experiment conditions
#' @param levelvector a vector of experimental conditions, not complusory for isothermal functions
#' @param nread number of reading channels or sample treatements
#' @param robustfitting whether to apply robust fitting, which is new since v.0.3.7
#' @param minireplicate number of replicates to keep in final data
#' @param withset whether there is set column to perform facet_grid
#' @param loess whether to perform curve fitting using loess model
#' @param dotconnect whether to simply dot connect the readings for each curve
#' @param xtransform, how should the x scale be transform, one value among
#' c("xsqrt","xcubert","xlog10","xlog2","xlinear")
#' @param xinterval a number indicating the numerical interval for linear x-axis
#' @param presetcolor whether to use the pre-defined color scheme
#' @param colorpanel a vector of customizable color scheme provided by the user
#' @param layout a vector indicating the panel layout for multi-panel plots per page
#' @param toplabel textual label at the top of the page
#' @param leftlabel textual label at the left side of the page
#' @param bottomlabel textual label at the bottom of the page
#' @param shadearea a four element vector indicating the xmin, xmax, ymin, ymax for x and y axis for shading
#' @param shadeoutlinecolor color used to shade the area, default value is black
#' @param shadefillcolor color used to shade the area, default value is gray90
#' @param returnplots whether to return the plot objects
#'
#'
#' @import tidyr
#' @import RColorBrewer
#' @import scales
#' @import drc
#' @importFrom plyr . dlply
#' @importFrom grid unit.c
#' @importFrom gridExtra grid.arrange arrangeGrob
#' @import ggplot2
#' @keywords internal
#'
#' @return a list of ggplot2 object
#' @examples \dontrun{
#' }
#'
#'
ms_isothermal_line_innerplot <- function(data, legenddata, levelvector, nread,
robustfitting, minireplicate, withset,
witherrorbar, loess, dotconnect,
xtransform, xinterval,
presetcolor, colorpanel, layout,
toplabel, leftlabel, bottomlabel,
shadearea, shadeoutlinecolor, shadefillcolor,
returnplots, outdir) {
if (xtransform %in% c("xsqrt","xcubert","xlinear")) {
xmin <- as.numeric(names(data)[3])
xmax <- as.numeric(names(data)[(nread+2)])
cat("The x axis scale range is:",c(xmin, xmax),"\n") # scale used for plotting
} else if (xtransform=="xlog10") {
first <- log10(as.numeric(names(data)[3]))
last <- log10(as.numeric(names(data)[(nread+2)]))
xmin <- (floor(first) + ceiling(first))/2
if ( xmin > first ) { xmin <- floor(first) }
if ( xmin == -Inf ) {
# to hardcode zero dose to two log units smaller than the minimal dose
xmin <- floor(log10(as.numeric(names(data)[4])/10))
sub <- as.numeric(names(data)[4])/100
cat("In log scale, the zero was substituted as", sub, "for regression.\n")
names(data) <- gsub("^0$", sub, names(data))
}
xmax <- (floor(last) + ceiling(last))/2
if (xmax < last) { xmax <- ceiling(last) }
cat("The log10 x axis scale range is:",c(xmin, xmax),"\n") # scale used for plotting
} else if (xtransform=="xlog2") {
first <- log2(as.numeric(names(data)[3]))
last <- log2(as.numeric(names(data)[(nread+2)]))
xmin <- (floor(first) + ceiling(first))/2
if ( xmin > first ) { xmin <- floor(first) }
if ( xmin == -Inf ) {
# to hardcode zero dose to two log units smaller than the minimal dose
xmin <- floor(log2(as.numeric(names(data)[4])/10))
sub <- as.numeric(names(data)[4])/100
cat("In log scale, the zero was substituted as", sub, "for regression.\n")
names(data) <- gsub("^0$", sub, names(data))
}
xmax <- (floor(last) + ceiling(last))/2
if (xmax < last) { xmax <- ceiling(last) }
cat("The log2 x axis scale range is:",c(xmin, xmax),"\n") # scale used for plotting
}
nametreatmentvector <- names(data)[3:(nread+2)]
numtreatmentvector <- as.numeric(nametreatmentvector)
cat("The points are as follows:",numtreatmentvector,"\n")
# y <- c(1.0, 0.98, 1.02, 1.1, 1.2, 1.25, 1.3, 1.4, 1.5, 2.1)
# x <- numtreatmentvector
# df <- data.frame(x, y)
# t1 <- drc::drm(y~x, data=df, fct=LL.4())
# print(t1)
if (xtransform %in% c("xsqrt","xcubert","xlinear")) {
treat_min <- xmin
treat_max <- xmax
} else if (xtransform=="xlog10") {
treat_min <- log(min(as.numeric(nametreatmentvector)),10)-0.05 # extra 10% room for buffering
treat_max <- log(max(as.numeric(nametreatmentvector)),10)+0.05
} else if (xtransform=="xlog2") {
treat_min <- log(min(as.numeric(nametreatmentvector)),2)-0.05 # extra 10% room for buffering
treat_max <- log(max(as.numeric(nametreatmentvector)),2)+0.05
}
message("Generating fitted plot file, pls wait...")
if (witherrorbar) {
if (withset) {
setpos <- grep("^set", names(data))
if ( length(setpos)==0 ) {stop("There is no set column to perform facet_grid")}
data_l <- tidyr::gather(data[ ,c(1:(nread+2),setpos)], treatment, reading, -id, -set, -condition)
data_l <- tidyr::separate(data_l, condition, into=c("sample","rep"), sep="\\.")
cdata <- plyr::ddply(data_l, c("id", "set", "sample", "treatment"), summarise,
N = length(reading),
mean = mean(reading),
sd = sd(reading),
se = sd / sqrt(N)
)
} else {
data_l <- tidyr::gather(data[ ,c(1:(nread+2))], treatment, reading, -id, -condition)
data_l <- tidyr::separate(data_l, condition, into=c("sample","rep"), sep="\\.")
cdata <- plyr::ddply(data_l, c("id", "sample", "treatment"), summarise,
N = length(reading),
mean = mean(reading),
sd = sd(reading),
se = sd / sqrt(N)
)
}
# print(head(cdata))
if (length(minireplicate)>0) {
cdata <- subset(cdata, N>=minireplicate)
}
#cdata <- na.omit(cdata) # to remove the measurements without enough replicates
if (withset) {
cdata <- tidyr::complete(cdata, id, set, sample, treatment)
} else {
cdata <- tidyr::complete(cdata, id, sample, treatment)
}
nsample <- length(unique(data_l$sample))
if (presetcolor & length(colorpanel)==0) {
if (nsample<=9) {colorpanel <- brewer.pal(9, "Set1")}
else {stop("The number of conditions in dataset exceeds the preset number of colors, pls provide a vector of colors in colorpanel")}
} else if (length(colorpanel) < nsample){
stop("The number of conditions in dataset exceeds the provided number of colors, pls provide enough colors in colorpanel")
}
if (length(levelvector)) {
cdata$sample <- factor(as.character(cdata$sample), levels=levelvector)
} else {
cdata$sample <- factor(cdata$sample, levels=unique(cdata$sample))
}
cdata$treatment <- as.numeric(as.character(cdata$treatment))
names(colorpanel) <- levels(cdata$sample)
} else { # for individual curves
if (withset) {
setpos <- grep("^set", names(data))
if ( length(setpos)==0 ) {stop("There is no set column to perform facet_grid")}
data_l <- tidyr::gather(data[ ,c(1:(nread+2),setpos)], treatment, reading, -id, -set, -condition)
} else {
data_l <- tidyr::gather(data[ ,c(1:(nread+2))], treatment, reading, -id, -condition)
}
if (length(levelvector)) {
data_l$condition<-factor(as.character(data_l$condition), levels=levelvector)
} else {
data_l$condition<-factor(data_l$condition, levels=sort(unique(data_l$condition)))
}
data_l$treatment <- as.numeric(as.character(data_l$treatment))
# print(head(data_l))
# print(summary(data_l))
# data_l<-data.frame(data_l)
# return(data_l)
ncond <- length(unique(data$condition)) # unique condition names including replicate info
data1 <- tidyr::separate(data, condition, into=c("condition", "replicate"), sep="\\.")
ncondition <- length(unique(data1$condition))
nreplicate <- length(unique(data1$replicate))
if (withset) {
uniquecond <- unique(data1[ ,c("set", "condition", "replicate")])
} else{
uniquecond <- unique(data1[ ,c("condition", "replicate")])
}
row.names(uniquecond) <- NULL
# print(uniquecond)
if (presetcolor & length(colorpanel)==0) {
if (nreplicate==2 & ncond<=8) {colorpanel <- c("blue1", "blue4", "red1", "red4", "green1", "green4", "darkgrey", "black")}
else if (nreplicate==2 & ncond<=12) {colorpanel <- brewer.pal(ncond, "Paired")}
else if (nreplicate==1 & ncond<=8) {colorpanel <- c("black", "red", "blue", "orange", "green", "purple", "yellow", "cyan")}
else if (nreplicate==1 & ncond<=12) {colorpanel <- brewer.pal(ncond, "Set3")}
else if (nreplicate==3 & ncond<=12) {colorpanel <- c("blue1", "blue4", "blueviolet", "red1", "red3", "brown",
"green1", "green3", "forestgreen", "dimgray", "gray30", "black")}
else if (nreplicate==4 & ncond<=8) {colorpanel <- c("#E0E0E0", "#BABABA", "#878787", "#4D4D4D", "#FDDBC7", "#F4A582", "#D6604D", "#B2182B")}
else {stop("The number of conditions in dataset exceeds the preset number of colors, pls provide a vector of colors in colorpanel")}
} else if (length(colorpanel) < ncond){
stop("The number of conditions in dataset exceeds the provided number of colors, pls provide enough colors in colorpanel")
}
#colorpanel <- c("blue1", "blue4", "red1", "red4", "green1", "green4", "yellow", "yellow4", "gray30", "black")
#colorpanel <- brewer_pal("qual")(length(levels(data_l$condition)))
names(colorpanel) <- levels(data_l$condition)
}
# print(colorpanel)
zz <- file(paste0(outdir,"/","curve_fitting_record.txt"), open="wt")
sink(zz, type = "message")
on.exit(sink(type="message"))
plotting1 <- function(d1) { # individual line plot function
q <- ggplot()
if (loess) {
q <- ggplot(d1, aes(x=treatment, y=reading, group=condition, colour=condition)) +
geom_point(aes(colour=condition), shape=20, size=2) +
geom_smooth(method="loess", aes(colour=condition), se=FALSE, size=0.5) +
scale_colour_manual(drop=FALSE, values=colorpanel)
} else if (dotconnect) {
q <- ggplot(d1, aes(x=treatment, y=reading, group=condition, colour=condition)) +
geom_point(aes(colour=condition), shape=20, size=2) +
geom_line(aes(colour=condition), size=0.5) +
scale_colour_manual(drop=FALSE, values=colorpanel)
} else {
# q <- ggplot(d1, aes(x=treatment, y=reading, group=condition, colour=condition)) +
# geom_point(aes(colour=condition), shape=20, size=2)
# if (withset) { q <- q + facet_grid(set~., drop=FALSE) }
# q <- q + geom_smooth(method=drc::"drm", formula=t1,
# method.args=list(fct=drc::LL.4()), aes(colour=condition), se=FALSE, size=0.5) +
# scale_colour_manual(drop=FALSE, values=colorpanel)
d1 <- droplevels(d1)
d1_list <- split(d1, d1$condition)
for (j in names(d1_list)) { # use levels of conditions, but not length of condition for individual proteins (length(d1_list))
# print(colorpanel[j])
q <- q + geom_point(data=d1_list[[j]], aes(x=treatment, y=reading), colour=colorpanel[j], shape=20, size=2)
d1_list[[j]] <- na.omit(d1_list[[j]])
if (grepl("3.3",R.version.string)) {
if (robustfitting) {
model.drm <- try(drm(reading~treatment, data=d1_list[[j]], fct=LL.4(), robust="median", na.action=na.omit,
control = drmc(noMessage=TRUE)), silent=FALSE)
} else {
model.drm <- try(drm(reading~treatment, data=d1_list[[j]], fct=LL.4(), na.action=na.omit,
control = drmc(noMessage=TRUE)), silent=FALSE)
}
} else {
if (robustfitting) {
model.drm <- try(drm(reading~treatment, data=d1_list[[j]], fct=LL.4(), robust="median", na.action=na.omit,
control = drmc(noMessage=TRUE)), silent=FALSE, outFile=zz)
} else {
model.drm <- try(drm(reading~treatment, data=d1_list[[j]], fct=LL.4(), na.action=na.omit,
control = drmc(noMessage=TRUE)), silent=FALSE, outFile=zz)
}
}
if(class(model.drm) != "try-error"){
if (xtransform %in% c("xsqrt","xcubert","xlinear")) {
mda <- data.frame(treatment=seq(treat_min,treat_max,length.out=100))
mda$reading <- predict(model.drm, mda)
} else if (xtransform=="xlog10") {
mda <- data.frame(treatment=10^(seq(treat_min,treat_max,length.out=100)))
mda$reading <- predict(model.drm,mda)
} else if (xtransform=="xlog2") {
mda <- data.frame(treatment=2^(seq(treat_min,treat_max,length.out=100)))
mda$reading <- predict(model.drm,mda)
}
q <- q + geom_line(data=mda, aes(x=treatment, y=reading), colour=colorpanel[j], size=0.5)
} else {
q <- q + geom_smooth(data=d1_list[[j]], method="lm", formula=y~poly(x,1), aes(x=treatment, y=reading), colour=colorpanel[j], se=FALSE, size=0.5)
}
}
}
if (min(abs(d1$reading),na.rm=T) <= 0.5) {
mi = 0
} else {
mi = 0.5
}
ma <- ceiling(2*max(abs(d1$reading), na.rm=T))/2
if (xtransform=="xlinear") {
if (!length(xinterval)) {
q <- q + coord_cartesian(xlim=c(xmin,xmax))
} else {
q <- q + coord_cartesian(xlim=c(xmin,xmax)) + scale_x_continuous(breaks=seq(xmin, xmax, xinterval))
}
} else if (xtransform=="xlog10") {
q <- q + coord_cartesian(xlim=c(10^xmin,10^xmax*1.5),ylim=c(mi-0.05,max(ma+0.1,2))) +
scale_x_log10(limits=c(10^xmin,10^xmax*1.5), breaks = trans_breaks("log10", function(x) 10^x),
labels = trans_format("log10", math_format(10^.x)))
} else if (xtransform=="xlog2") {
q <- q + coord_cartesian(xlim=c(2^xmin,2^xmax*1.5),ylim=c(mi-0.05,max(ma+0.1,2))) +
scale_x_log2(limits=c(2^xmin,2^xmax*1.5), breaks = trans_breaks("log2", function(x) 2^x),
labels = trans_format("log2", math_format(2^.x)))
} else if (xtransform=="xsqrt") {
q <- q + coord_cartesian(xlim=c(0,xmax),ylim=c(mi-0.05,max(ma+0.1,2))) +
scale_x_continuous(trans=sqrt_trans(), breaks=trans_breaks("sqrt", function(x) x^2))
} else if (xtransform=="xcubert") {
cuberoot_trans = function() trans_new("cuberoot", function(x) x^(1/3), function(x) x^3)
q <- q + coord_cartesian(xlim=c(0,xmax),ylim=c(mi-0.05,max(ma+0.1,2))) +
scale_x_continuous(trans=cuberoot_trans(), breaks=seq(0,xmax,xinterval))
}
if (ma<=3.0) {
q <- q + scale_y_continuous(breaks=seq(mi,ma,0.5))
} else if (ma<=5.0) {
q <- q + scale_y_continuous(breaks=seq(mi,ma,1))
} else {
q <- q + scale_y_continuous(breaks=seq(mi,ma,2))
}
if (length(shadearea)) {
q <- q + geom_rect(mapping=aes(xmin=shadearea[1], xmax=shadearea[2],
ymin=shadearea[3], ymax=shadearea[4]),
color=shadeoutlinecolor, fill=shadefillcolor, alpha=0.1)
}
q <- q + ggtitle(as.character(unique(d1$id))) + labs(x=" ", y=" ")
q <- q + theme_classic() +
theme(
text = element_text(size=10),
strip.text.x = element_text(size=5),
plot.title = element_text(hjust=0.5, size=rel(0.7)),
axis.text.x = element_text(angle=45, hjust=1),
legend.position = "none",
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
strip.background = element_blank(),
strip.text = element_blank(),
axis.line.x = element_line(),
axis.line.y = element_line(),
aspect.ratio = 1
)
return(q)
}
plotting2 <- function(d1) { # plot with error bar function
q <- ggplot()
if (loess) {
q <- ggplot(d1, aes(x=treatment, y=mean, group=sample, colour=sample)) +
geom_point(aes(colour=sample), shape=20, size=2) +
geom_smooth(method="loess", aes(colour=sample), se=FALSE, size=0.5) +
scale_colour_manual(drop=FALSE, values=colorpanel)
q <- q + geom_errorbar(d1, aes(ymin=mean-se, ymax=mean+se, colour=sample), width=.2,
position=position_dodge(.9))
} else if (dotconnect) {
q <- ggplot(d1, aes(x=treatment, y=mean, group=sample, colour=sample)) +
geom_point(aes(colour=sample), shape=20, size=2) +
geom_line(aes(colour=sample), size=0.5) +
scale_colour_manual(drop=FALSE, values=colorpanel)
q <- q + geom_errorbar(d1, aes(ymin=mean-se, ymax=mean+se, colour=sample), width=.2,
position=position_dodge(.9))
} else {
d1 <- droplevels(d1)
d1_list <- split(d1, d1$sample)
for (j in names(d1_list)) {
# use levels of conditions, but not length of condition for individual proteins (length(d1_list))
q <- q + geom_point(data=d1_list[[j]], aes(x=treatment, y=mean), colour=colorpanel[j], shape=20, size=2)
d1_list[[j]] <- na.omit(d1_list[[j]])
if (grepl("3.3",R.version.string)) {
if (robustfitting) {
model.drm <- try(drm(mean~treatment, data=d1_list[[j]], fct=LL.4(), robust="median", na.action=na.omit,
control = drmc(noMessage=TRUE)), silent=FALSE)
} else {
model.drm <- try(drm(mean~treatment, data=d1_list[[j]], fct=LL.4(), na.action=na.omit,
control = drmc(noMessage=TRUE)), silent=FALSE)
}
} else {
if (robustfitting) {
model.drm <- try(drm(mean~treatment, data=d1_list[[j]], fct=LL.4(), robust="median", na.action=na.omit,
control = drmc(noMessage=TRUE)), silent=FALSE, outFile=zz)
} else {
model.drm <- try(drm(mean~treatment, data=d1_list[[j]], fct=LL.4(), na.action=na.omit,
control = drmc(noMessage=TRUE)), silent=FALSE, outFile=zz)
}
}
if (class(model.drm) != "try-error") {
if (xtransform %in% c("xsqrt","xcubert","xlinear")) {
mda <- data.frame(treatment=seq(treat_min,treat_max,length.out=100))
mda$reading <- predict(model.drm, mda)
} else if (xtransform=="xlog10") {
mda <- data.frame(treatment=10^(seq(treat_min,treat_max,length.out=100)))
mda$reading <- predict(model.drm,mda)
} else if (xtransform=="xlog2") {
mda <- data.frame(treatment=2^(seq(treat_min,treat_max,length.out=100)))
mda$reading <- predict(model.drm,mda)
}
q <- q + geom_line(data=mda, aes(x=treatment, y=reading), colour=colorpanel[j], size=0.5)
} else {
q <- q + geom_smooth(data=d1_list[[j]], method="lm", formula=y~poly(x,1),
aes(x=treatment, y=mean), colour=colorpanel[j], se=FALSE, size=0.5)
}
q <- q + geom_errorbar(data=d1_list[[j]], aes(x=treatment, ymin=mean-se, ymax=mean+se),
colour=colorpanel[j], width=.2, # Width of the error bars
position=position_dodge(.9))
}
}
if (min(abs(d1$mean),na.rm=T) <= 0.5) {
mi = 0
} else {
mi = 0.5
}
ma <- ceiling(2*max(abs(d1$mean), na.rm=T))/2
if (xtransform=="xlinear") {
if (!length(xinterval)) {
q <- q + coord_cartesian(xlim=c(xmin,xmax))
} else {
q <- q + coord_cartesian(xlim=c(xmin,xmax)) + scale_x_continuous(breaks=seq(xmin, xmax, xinterval))
}
} else if (xtransform=="xlog10") {
q <- q + coord_cartesian(xlim=c(10^xmin,10^xmax*1.5),ylim=c(mi-0.05,max(ma+0.1,2))) +
scale_x_log10(limits=c(10^xmin,10^xmax*1.5), breaks = trans_breaks("log10", function(x) 10^x),
labels = trans_format("log10", math_format(10^.x)))
} else if (xtransform=="xlog2") {
q <- q + coord_cartesian(xlim=c(2^xmin,2^xmax*1.5),ylim=c(mi-0.05,max(ma+0.1,2))) +
scale_x_log2(limits=c(2^xmin,2^xmax*1.5), breaks = trans_breaks("log2", function(x) 2^x),
labels = trans_format("log2", math_format(2^.x)))
} else if (xtransform=="xsqrt") {
q <- q + coord_cartesian(xlim=c(0,xmax),ylim=c(mi-0.05,max(ma+0.1,2))) +
scale_x_continuous(trans=sqrt_trans(), breaks=trans_breaks("sqrt", function(x) x^2))
} else if (xtransform=="xcubert") {
cuberoot_trans = function() trans_new("cuberoot", function(x) x^(1/3), function(x) x^3)
q <- q + coord_cartesian(xlim=c(0,xmax),ylim=c(mi-0.05,max(ma+0.1,2))) +
scale_x_continuous(trans=cuberoot_trans(), breaks=seq(0,xmax,xinterval))
}
if (ma<=3.0) {
q <- q + scale_y_continuous(breaks=seq(mi,ma,0.5))
} else if (ma<=5.0) {
q <- q + scale_y_continuous(breaks=seq(mi,ma,1))
} else {
q <- q + scale_y_continuous(breaks=seq(mi,ma,2))
}
if (length(shadearea)) {
q <- q + geom_rect(mapping=aes(xmin=shadearea[1], xmax=shadearea[2],
ymin=shadearea[3], ymax=shadearea[4]),
color=shadeoutlinecolor, fill=shadefillcolor, alpha=0.1)
}
q <- q + ggtitle(as.character(unique(d1$id))) + labs(x=" ", y=" ")
q <- q + theme_classic() +
theme(
text = element_text(size=10),
strip.text.x = element_text(size=5),
plot.title = element_text(hjust=0.5, size=rel(0.7)),
axis.text.x = element_text(angle=45, hjust=1),
legend.position = "none",
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
strip.background = element_blank(),
strip.text = element_blank(),
axis.line.x = element_line(),
axis.line.y = element_line(),
aspect.ratio = 1
)
return(q)
}
if (witherrorbar) {
plots <- plyr::dlply(cdata, plyr::.(id), .fun = plotting2)
plotlegend <- ms_isothermal_legend2(legenddata, levelvector,
nread, nsample, presetcolor, colorpanel)
} else {
plots <- plyr::dlply(data_l, plyr::.(id), .fun = plotting1)
plotlegend <- ms_isothermal_legend(legenddata, levelvector,
nread, nreplicate, presetcolor, colorpanel)
}
zz <- paste0(outdir,"/","curve_fitting_record.txt")
if (file.exists(zz)) { file.remove(zz) }
if (returnplots) { return(plots) }
legend <- plotlegend$legend
lheight <- plotlegend$lheight
params <- list(nrow=layout[1], ncol=layout[2])
n <- with(params, nrow*ncol)
## add one page if division is not complete
pages <- length(plots) %/% n + as.logical(length(plots) %% n)
groups <- split(seq_along(plots), gl(pages, n, length(plots)))
pl <- lapply(names(groups), function(i){
gridExtra::grid.arrange(
do.call(arrangeGrob,
c(plots[groups[[i]]], params, top=toplabel,
left=leftlabel, bottom=bottomlabel)),
legend,
ncol = 1,
heights = grid::unit.c(unit(1, "npc") - lheight, lheight))
})
class(pl) <- c("arrangelist", "ggplot", class(pl))
return(pl)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.