diff --git a/dev/app.R b/dev/app.R index 4b97b26..3d99621 100644 --- a/dev/app.R +++ b/dev/app.R @@ -15,10 +15,13 @@ library(purrr) library(gslnls) library(tidyverse) library(ggplot2) +library(ggExtra) +library(gtable) library(reshape2) library(openxlsx) library(DT) library(ggpubr) +library(grid) library(gridExtra) library(drc) library(twopartm) @@ -583,17 +586,16 @@ server <- function(input, output, session) { box( title = "Upload multiple worksheets", status = "warning", solidHeader = TRUE, width = 12, "Please upload your EXCEL file here", fileInput("MiFile", "", accept = ".xlsx") - ) + ), + sliderInput("dilslider", "Adjust the dilutions(+-change in %)", min = -100,max=100, value=0, step=1, round=0), + checkboxInput("fixupper","Fix highest concentration (if unticked, the center is fixed)",FALSE) ) ), mainPanel( tabsetPanel( id = "tabs", tabPanel("4pl", - box( - title = "ANOVA table", status = "primary", solidHeader = TRUE, width = 12, - tableOutput("Anovatab") - ), + column(12, # h3("Confidence intervals"), @@ -604,42 +606,44 @@ server <- function(input, output, session) { plotOutput("sigPlotREF"), plotOutput("sigPlotTEST"), #selectInput(inputId = "scenario", label = "Select an 'optimal' scenario:", choices = c("scenario 2", "scenario 3", "scenario 6", "steep slope")) - - plotOutput("plotfordilutions"), - # h4("in grey: most extreme bend point lines of theoretical samples with 50% and 200% potency"), - # sliderInput("dilslider", "Adjust the dilutions(+-change in %)", min = -100, max = 100, value = 0, step = 1, round = 0), - # checkboxInput("fixupper", "Fix highest concentration (if unticked, the center is fixed)", FALSE), - # h5("Dilution factors"), - # tableOutput("adjlogdil"), - # "Short guidance: wider dilution ranges increase the CIs of rel. potency, and decrease the CIs of upper and lower asymptote ratios, as well as Hill's slope ratios", br(), - # "Narrower dilution ranges decrease the CIs of rel. potency, and increase the CIs of upper and lower asymptote ratios, ands Hill's slope ratios", - ), - column( - 3, - h3("Bend points"), - tableOutput("bps"), - tableOutput("extremebps"), - h4("Explanation of the plots") - ) + )), + tabPanel("Dilution slider", + h4("Adjust the dilutions"), + plotOutput("plotfordilutions"), + h5("Dilution factors"), + box( + title = "Adjusted dilutions", status = "primary", solidHeader = TRUE, width = 12, collapsible = TRUE, + tableOutput("adjlogdil") + ), + "Short guidance: wider dilution ranges increase the CIs of rel. potency, and decrease the CIs of upper and lower asymptote ratios, as well as Hill's slope ratios", br(), + "Narrower dilution ranges decrease the CIs of rel. potency, and increase the CIs of upper and lower asymptote ratios, ands Hill's slope ratios", + ), tabPanel("Histograms", h4("Histograms of parameters"), plotOutput("histCIs"), - plotOutput("histEC50REF"), - plotOutput("histLasREF"), - plotOutput("histUasREF"), + column(6, + plotOutput("histEC50REF"), + plotOutput("histLasREF"), + plotOutput("histUasREF") + ), + column(6, + plotOutput("histEC50TEST"), + plotOutput("histLasTEST"), + plotOutput("histUasTEST") + ) ), tabPanel( "Report", - h4("Settings for report") - ) + h4("Settings for report")) + ) - ) - ) - ) - ) + ) # main panel + ) # sidebar Layout + ) # tabpanel + ) # navbarPage }) @@ -696,7 +700,7 @@ server <- function(input, output, session) { reset(id = "") # from shinyjs package }) - +#### input optim XL file ---- observe({ if (!is.null(input$MiFile)) { MinFile <- input$MiFile @@ -2114,7 +2118,7 @@ server <- function(input, output, session) { #### Dilutions Simulator ---- - output$plotfordilutions <- renderPlot({ + observe({ if (!is.null(Dat$Mws)) { @@ -2125,55 +2129,55 @@ server <- function(input, output, session) { for (N_WS in 1:length(AllXL)) { - datWS <- as.data.frame(AllXL[[N_WS]]) - - cn <- colnames(datWS) - logI <- grep("log|ln", cn) - logDoseI <- grep("log_dose", cn) - if (length(logI) > 0 & length(logDoseI) == 0) { - datWS$log_dose <- datWS[, logI] - datWS2 <- datWS[, -logI] - CORro <- cor(datWS$log_dose, datWS[, 3]) - } else if (length(logI) == 0 & length(logDoseI) == 0) { - Ind <- grep(".ilution|.ose|.onc", cn) - datWS$log_dose <- log(datWS[, Ind]) - CORro <- cor(datWS[, Ind], datWS[, 3]) - datWS2 <- datWS[, -Ind] - } else if (length(logI) > 0 & length(logDoseI) > 0) { - datWS2 <- datWS - CORro <- cor(datWS[, logI], datWS[, 3]) - } - Dat$datWS2 <- datWS2 - - FITs <- Fitting_FUNC(datWS2, TransFlag = F) - - pot_est <- FITs[[3]] - potEstL[[N_WS]] <- pot_est - potU_est <- FITs[[4]] - # unrestricted - SU_mu <- FITs[[2]] - URMcoefs1 <- SU_mu$coefficients - URMcoefs <- t(matrix(unlist(URMcoefs1[,1]))) - URMcoefs_ <- cbind(AllSheets[[N_WS]], URMcoefs) - URMcoefsL[[N_WS]] <- URMcoefs_ - - SU_mr <- FITs[[1]] - RMcoefs1 <- SU_mr$coefficients - RMcoefs <- t(matrix(unlist(RMcoefs1[,1]))) - RMcoefs_ <- cbind(AllSheets[[N_WS]], RMcoefs) - RMcoefsL[[N_WS]] <- RMcoefs_ - - - - X <- seq(min(datWS2$log_dose), max(datWS2$log_dose), 0.1) - sigRef <- URMcoefs[1,1] + (URMcoefs[1,3]-URMcoefs[1,1])/(1+exp(URMcoefs[1,2]*(URMcoefs[1,4]-X))) - sigTest1 <- URMcoefs[1,5] + (URMcoefs[1,7]-URMcoefs[1,5])/(1+exp(URMcoefs[1,6]*(URMcoefs[1,4] - URMcoefs[1,8]-X))) - #browser() - dfPlotsigRef <- data.frame(X=X, sigRef = sigRef, Sheet = AllSheets[[N_WS]]) - dfPlotsigTest <- data.frame(X=X, sigTest = sigTest1, Sheet = AllSheets[[N_WS]]) - - if (!exists("SIGrefDF")) SIGrefDF <- dfPlotsigRef else SIGrefDF <- rbind(SIGrefDF, dfPlotsigRef) - if (!exists("SIGtestDF")) SIGtestDF <- dfPlotsigTest else SIGtestDF <- rbind(SIGtestDF,dfPlotsigTest) + datWS <- as.data.frame(AllXL[[N_WS]]) + + cn <- colnames(datWS) + logI <- grep("log|ln", cn) + logDoseI <- grep("log_dose", cn) + if (length(logI) > 0 & length(logDoseI) == 0) { + datWS$log_dose <- datWS[, logI] + datWS2 <- datWS[, -logI] + CORro <- cor(datWS$log_dose, datWS[, 3]) + } else if (length(logI) == 0 & length(logDoseI) == 0) { + Ind <- grep(".ilution|.ose|.onc", cn) + datWS$log_dose <- log(datWS[, Ind]) + CORro <- cor(datWS[, Ind], datWS[, 3]) + datWS2 <- datWS[, -Ind] + } else if (length(logI) > 0 & length(logDoseI) > 0) { + datWS2 <- datWS + CORro <- cor(datWS[, logI], datWS[, 3]) + } + Dat$datWS2 <- datWS2 + + FITs <- Fitting_FUNC(datWS2, TransFlag = F) + + pot_est <- FITs[[3]] + potEstL[[N_WS]] <- pot_est + potU_est <- FITs[[4]] + # unrestricted + SU_mu <- FITs[[2]] + URMcoefs1 <- SU_mu$coefficients + URMcoefs <- t(matrix(unlist(URMcoefs1[,1]))) + URMcoefs_ <- cbind(AllSheets[[N_WS]], URMcoefs) + URMcoefsL[[N_WS]] <- URMcoefs_ + + SU_mr <- FITs[[1]] + RMcoefs1 <- SU_mr$coefficients + RMcoefs <- t(matrix(unlist(RMcoefs1[,1]))) + RMcoefs_ <- cbind(AllSheets[[N_WS]], RMcoefs) + RMcoefsL[[N_WS]] <- RMcoefs_ + + + + X <- seq(min(datWS2$log_dose), max(datWS2$log_dose), 0.1) + sigRef <- URMcoefs[1,1] + (URMcoefs[1,3]-URMcoefs[1,1])/(1+exp(URMcoefs[1,2]*(URMcoefs[1,4]-X))) + sigTest1 <- URMcoefs[1,5] + (URMcoefs[1,7]-URMcoefs[1,5])/(1+exp(URMcoefs[1,6]*(URMcoefs[1,4] - URMcoefs[1,8]-X))) + #browser() + dfPlotsigRef <- data.frame(X=X, sigRef = sigRef, Sheet = AllSheets[[N_WS]]) + dfPlotsigTest <- data.frame(X=X, sigTest = sigTest1, Sheet = AllSheets[[N_WS]]) + + if (!exists("SIGrefDF")) SIGrefDF <- dfPlotsigRef else SIGrefDF <- rbind(SIGrefDF, dfPlotsigRef) + if (!exists("SIGtestDF")) SIGtestDF <- dfPlotsigTest else SIGtestDF <- rbind(SIGtestDF,dfPlotsigTest) } #for N_WS @@ -2205,18 +2209,83 @@ server <- function(input, output, session) { x_UA <- max(X); x_LA <- min(X) } else { x_UA <- min(X); x_LA <- max(X) } + #browser() + BoxDF <- data.frame(EC50REF = EC50REF, EC50TEST = EC50TEST, LasREF = LasREF, UasREF = UasREF) + p1 <- ggplot(SIGrefDF, aes(x=X, y=sigRef, col=as.factor(Sheet))) + geom_line() + annotate("text", label="x", x=x_UA, y=UasREF, alpha=0.2) + annotate("text", label="o", x=x_LA, y=LasREF, alpha=0.2) + geom_vline(xintercept = EC50REF, alpha = 0.2) + + scale_x_continuous(expand = c(0, 0)) + + scale_y_continuous(expand = c(0, 0)) + + expand_limits(y = c(min(SIGrefDF$sigRef) - 0.1 * diff(range(SIGrefDF$sigRef)), + max(SIGrefDF$sigRef) + 0.1 * diff(range(SIGrefDF$sigRef)))) + + expand_limits(x = c(min(SIGrefDF$X) - 0.1 * diff(range(SIGrefDF$X)), + max(SIGrefDF$X) + 0.1 * diff(range(SIGrefDF$X)))) + xlab("dilutions") + - ggtitle("Plot of all calculated reference fits (unrestricted model, in gray vertical lines: EC50)") + + #ggtitle("Plot of all calculated reference fits (unrestricted model, in gray vertical lines: EC50)") + theme_bw() + theme(axis.text = element_text(face = "bold", size = 15), - plot.title = element_text(size = 15, face = "bold")) - - output$sigPlotREF <- renderPlot({ p1 }) + plot.title = element_text(size = 15, face = "bold"), + plot.margin = unit(c(0.2, 0.2, 0.5, 0.5), "lines")) + # Horizontal marginal boxplot - to appear at the top of the chart + pBox_hor <- ggplot( BoxDF, aes(x = factor(1), y = EC50REF)) + + geom_boxplot(outlier.colour = NA) + + geom_jitter(position = position_jitter(width = 0.05)) + + scale_y_continuous(expand = c(0, 0)) + + expand_limits(y = c(min(SIGrefDF$X) - 0.1 * diff(range(SIGrefDF$X)), + max(SIGrefDF$X) + 0.1 * diff(range(SIGrefDF$X)))) + + coord_flip() + + theme_bw() + + theme(axis.text = element_blank(), + axis.title = element_blank(), + axis.ticks = element_blank(), + plot.margin = unit(c(1, 0.2, -0.5, 0.5), "lines")) + + # Vertical marginal boxplot - to appear at the right of the chart + pBox_ver <- ggplot(BoxDF, aes(x = factor(1), y = UasREF)) + + geom_boxplot(outlier.colour = NA) + + geom_jitter(position = position_jitter(width = 0.05)) + + scale_y_continuous(expand = c(0, 0)) + + expand_limits(y = c(min(SIGrefDF$sigRef) - 0.1 * diff(range(SIGrefDF$sigRef)), + max(SIGrefDF$sigRef) + 0.1 * diff(range(SIGrefDF$sigRef)))) + + theme_bw() + + theme(axis.text = element_blank(), + axis.title = element_blank(), + axis.ticks = element_blank(), + plot.margin = unit(c(0.2, 1, 0.5, -0.5), "lines")) + + #browser() + gt1 <- ggplot_gtable(ggplot_build(p1)) + gt2 <- ggplot_gtable(ggplot_build(pBox_hor)) + gt3 <- ggplot_gtable(ggplot_build(pBox_ver)) + + # Get maximum widths and heights + maxWidth <- unit.pmax(gt1$widths[2:3], gt2$widths[2:3]) + maxHeight <- unit.pmax(gt1$heights[4:5], gt3$heights[4:5]) + + # Set the maximums in the gtables for gt1, gt2 and gt3 + gt1$widths[2:3] <- as.list(maxWidth) + gt2$widths[2:3] <- as.list(maxWidth) + + gt1$heights[4:5] <- as.list(maxHeight) + gt3$heights[4:5] <- as.list(maxHeight) + # Create a new gtable + gt <- gtable(widths = unit(c(7, 1), "null"), height = unit(c(1, 7), "null")) + + # Instert gt1, gt2 and gt3 into the new gtable + gt <- gtable_add_grob(gt, gt1, 2, 1) + gt <- gtable_add_grob(gt, gt2, 1, 1) + gt <- gtable_add_grob(gt, gt3, 2, 2) + + # grid.rect(x = 0.5, y = 0.5, height = 0.995, width = 0.995, default.units = "npc", + # gp = gpar(col = "black", fill = NA, lwd = 1)) + # And render the plot + grid.newpage() +#browser() + + output$sigPlotREF <- renderPlot({ grid.draw(gt) }) Dat$sigPlotREF <- p1 # @@ -2267,35 +2336,44 @@ server <- function(input, output, session) { output$histUasREF <- renderPlot({ hist(UasREF, col="darkturquoise", border="white",main = 'Histogram of upper asymptotes REF') }) + + output$histEC50TEST <- renderPlot({ + hist(EC50TEST, col="steelblue", border="white", main = 'Histogram of EC50TEST') + }) + output$histLasTEST <- renderPlot({ + hist(LasTEST, col="violet", border="white",main = 'Histogram of lower asymptotes TEST') + }) + output$histUasTEST <- renderPlot({ + hist(UasTEST, col="darkturquoise", border="white",main = 'Histogram of upper asymptotes TEST') + }) Dat$histEC50REF <- hist(EC50REF, col="steelblue", border="white", main = 'Histogram of EC50REF') Dat$histLasREF <- hist(LasREF, col="violet", border="white", main = 'Histogram of EC50REF') Dat$histUasREF <- hist(UasREF, col="darkturquoise", border="white", main = 'Histogram of EC50REF') - # dils <- tab$log_dose - # min_y <- min(tab[, 1:3]) - # max_y <- max(tab[, 1:3]) - # - # if (input$fixupper) { - # dils_av <- dils - max(dils) - # dils_av_ <- dils_av * (input$dilslider / 100 + 1) - # dils2 <- round(dils_av_ + max(dils), 4) - # dilfactors <- 1 / exp(dils2 - lag(dils2)) - # } else { - # if (!is.null(Dat$cfordils)) { - # av <- Dat$cfordils - # } else { - # av <- (min(dils) + max(dils)) / 2 - # } - # dils_av <- dils - av - # dils_avsc <- dils_av * (input$dilslider / 100 + 1) - # dils2 <- dils_avsc + av - # dilfactors <- 1 / exp(dils2 - lag(dils2)) - # } + tab <- AllXL[[1]] + dils <- tab$log_dose + min_y <- min(tab[, 1:2]) + max_y <- max(tab[, 1:2]) + + if (input$fixupper) { + dils_av <- dils - max(dils) + dils_av_ <- dils_av * (input$dilslider / 100 + 1) + dils2 <- round(dils_av_ + max(dils), 4) + dilfactors <- 1 / exp(dils2 - lag(dils2)) + } else { + if (!is.null(EC50TEST)) { + av <- mean(EC50TEST, na.rm = TRUE) + } else { + av <- (min(dils) + max(dils)) / 2 + } + dils_av <- dils - av + dils_avsc <- dils_av * (input$dilslider / 100 + 1) + dils2 <- dils_avsc + av + dilfactors <- 1 / exp(dils2 - lag(dils2)) + } - - - #Dat$newDils <- dils2 + Dat$newDils <- dils2 #sigmoid <- sigmoid() @@ -2333,111 +2411,65 @@ server <- function(input, output, session) { # # pl_df <- cbind(dils2, SAMPLE50, SAMPLE200) # - # - # output$adjlogdil <- renderTable({ - # adjlogdilfactors <- round(dilfactors, 3) - # adjlogdils <- round(dils2, 3) - # adjdils <- round(exp(dils2), 3) - # DilsTable <- data.frame( - # "adjusted ln(dilutions)" = adjlogdils, - # "adjusted ln_dilution_factors" = adjlogdilfactors, - # "adjusted dilutions" = adjdils - # ) - # DilsTable - # }) - # if (!is.null(Dat$p2)) { - # p2 <- Dat$p2 - # p_dil <- p2 + - # annotate("pointrange", x = dils2, y = rep(min_y, length(dils2)), xmin = min(dils2), xmax = max(dils2)) + - # annotate("text", x = dils2, y = rep(min_y + (max_y - min_y) * 0.05, length(dils2)), label = as.character(round(dils2, 3))) + - # annotate("text", - # x = dils2[-1] + (max(dils2) - min(dils2)) * 0.05, - # y = rep(min_y + (max_y - min_y) * 0.1, length(dils2[-1])), - # label = as.character(round(dilfactors[-1], 3)) - # ) + - # geom_line( - # data = as.data.frame(pl_df), aes(x = dils2, y = SAMPLE50), color = "grey15", linetype = 2, - # inherit.aes = F - # ) + - # geom_line( - # data = as.data.frame(pl_df), aes(x = dils2, y = SAMPLE200), color = "grey15", linetype = 2, - # inherit.aes = F - # ) + - # geom_vline(xintercept = c(Xbend50, Xbend200), col = "grey15", linetype = 2) + - # { - # if (input$scenario == "scenario 6") { - # annotate("pointrange", - # x = optdils2, y = rep(min_y + (max_y - min_y) * 0.2, length(optdils2)), - # xmin = min(optdils2), xmax = max(optdils2), color = "seagreen" - # ) - # } - # } + - # { - # if (input$scenario == "scenario 6") { - # annotate("text", - # x = optdils2, y = rep(min_y + (max_y - min_y) * 0.25, length(optdils2)), - # label = as.character(round(optdils2, 3)), color = "seagreen" - # ) - # } - # } + - # { - # if (input$scenario == "scenario 2") { - # annotate("pointrange", - # x = optdils, y = rep(min_y + (max_y - min_y) * 0.2, length(optdils)), - # xmin = min(optdils), xmax = max(optdils), color = "seagreen" - # ) - # } - # } + - # { - # if (input$scenario == "scenario 2") { - # annotate("text", - # x = optdils, y = rep(min_y + (max_y - min_y) * 0.25, length(optdils)), - # label = as.character(round(optdils, 3)), color = "seagreen" - # ) - # } - # } + - # { - # if (input$scenario == "scenario 3") { - # annotate("pointrange", - # x = optdils_3, y = rep(min_y + (max_y - min_y) * 0.2, length(optdils_3)), - # xmin = min(optdils_3), xmax = max(optdils_3), color = "seagreen" - # ) - # } - # } + - # { - # if (input$scenario == "scenario 3") { - # annotate("text", - # x = optdils_3, y = rep(min_y + (max_y - min_y) * 0.25, length(optdils_3)), - # label = as.character(round(optdils_3, 3)), color = "seagreen" - # ) - # } - # } + - # { - # if (input$scenario == "steep slope") { - # annotate("pointrange", - # x = optdils3, y = rep(min_y + (max_y - min_y) * 0.2, length(optdils3)), - # xmin = min(optdils3), xmax = max(optdils3), color = "seagreen" - # ) - # } - # } + - # { - # if (input$scenario == "steep slope") { - # annotate("text", - # x = optdils3, y = rep(min_y + (max_y - min_y) * 0.25, length(optdils3)), - # label = as.character(round(optdils3, 3)), color = "seagreen" - # ) - # } - # } + - # annotate("text", - # x = optdils[1], y = (max_y + min_y) * 0.5, - # label = paste("in green: optimal \n dilutions acc. to Whitepaper\n", input$scenario), color = "seagreen", - # size = 14 / .pt, fontface = "bold" - # ) - # } - #print(p_dil) - # } - } + output$adjlogdil <- renderTable({ + adjlogdilfactors <- round(dilfactors, 3) + adjlogdils <- round(dils2, 3) + adjdils <- round(exp(dils2), 3) + DilsTable <- data.frame( + "adjusted ln(dilutions)" = adjlogdils, + "adjusted ln_dilution_factors" = adjlogdilfactors, + "adjusted dilutions" = adjdils + ) + DilsTable + }) + + if (!is.null(p2)) { + #p2 <- Dat$p2 + p_dil <- p2 + + annotate("pointrange", x = dils2, y = rep(min_y, length(dils2)), xmin = min(dils2), xmax = max(dils2)) + + annotate("text", x = dils2, y = rep(min_y + (max_y - min_y) * 0.05, length(dils2)), label = as.character(round(dils2, 3))) + + annotate("text", + x = dils2[-1] + (max(dils2) - min(dils2)) * 0.05, + y = rep(min_y + (max_y - min_y) * 0.1, length(dils2[-1])), + label = as.character(round(dilfactors[-1], 3))) + # geom_line( + # data = as.data.frame(pl_df), aes(x = dils2, y = SAMPLE50), color = "grey15", linetype = 2, + # inherit.aes = F + # ) + + # geom_line( + # data = as.data.frame(pl_df), aes(x = dils2, y = SAMPLE200), color = "grey15", linetype = 2, + # inherit.aes = F + # ) + + # geom_vline(xintercept = c(Xbend50, Xbend200), col = "grey15", linetype = 2) + + # { if (input$scenario == "scenario 6") { + # annotate("pointrange", + # x = optdils2, y = rep(min_y + (max_y - min_y) * 0.2, length(optdils2)), + # xmin = min(optdils2), xmax = max(optdils2), color = "seagreen" + # ) + # } + # } + + # { + # if (input$scenario == "scenario 6") { + # annotate("text", + # x = optdils2, y = rep(min_y + (max_y - min_y) * 0.25, length(optdils2)), + # label = as.character(round(optdils2, 3)), color = "seagreen" + # ) + # } + # } + + + + # annotate("text", + # x = optdils[1], y = (max_y + min_y) * 0.5, + # label = paste("in green: optimal \n dilutions acc. to Whitepaper\n", input$scenario), color = "seagreen", + # size = 14 / .pt, fontface = "bold" + # ) + output$plotfordilutions <- renderPlot({ + print(p_dil) + }) + + } # if (!is.null(p2)) + } # if !is.null Dat$Mws })