report generation after new setup of repo; linearity XL report improved

R/Global.R

@@ -280,14 +280,14 @@ plot_f <- function(dat, TransFlag = F) { # sigmoid,det_sig,
   s_mr <- MODLS[[1]]
   a <- s_mr$coefficients["a", 1]
   b <- s_mr$coefficients["b", 1]
-  cs <- s_mr$coefficients["cs", 1]
+  c <- s_mr$coefficients["cs", 1]
   d <- s_mr$coefficients["d", 1]
   r <- s_mr$coefficients["r", 1]
 
   log_dose <- unique(all_l$log_dose)
   seq_x <- seq(min(log_dose), max(log_dose), 0.1)
-  SAMPLE <- a + (d - a) / (1 + exp(b * ((cs - r) - seq_x)))
-  REF <- a + (d - a) / (1 + exp(b * ((cs) - seq_x)))
+  SAMPLE <- a + (d - a) / (1 + exp(b * ((c - r) - seq_x)))
+  REF <- a + (d - a) / (1 + exp(b * (c - seq_x)))
 
   s_mu <- MODLS[[2]]
 
@@ -310,16 +310,16 @@ plot_f <- function(dat, TransFlag = F) { # sigmoid,det_sig,
   all_l2$readout[all_l2$readout < 0] <- 0.01
   all_l2$readouttrans <- log(all_l2$readout)
   slopeEC50 <- b * (d - a) / 4
-  Intercept <- a + (d - a) / 2 - b * (d - a) / 4 * cs
+  Intercept <- a + (d - a) / 2 - b * (d - a) / 4 * c
   # browser()
-  Xbendl3 <- cs - (1.5434 / b)
-  Xbendu3 <- cs + (1.5434 / b)
-  XbendlT <- cs - r - (1.5434 / b)
-  XbenduT <- cs - r + (1.5434 / b)
-  XasymplS <- cs - (3 / b)
-  XasympuS <- cs + (3 / b)
-  XasymplT <- cs - r - (3 / b)
-  XasympuT <- cs - r + (3 / b)
+  Xbendl3 <- c - (1.5434 / b)
+  Xbendu3 <- c + (1.5434 / b)
+  XbendlT <- c - r - (1.5434 / b)
+  XbenduT <- c - r + (1.5434 / b)
+  XasymplS <- c - (3 / b)
+  XasympuS <- c + (3 / b)
+  XasymplT <- c - r - (3 / b)
+  XasympuT <- c - r + (3 / b)
   bendpoints <- c(
     bendREF_lower = round(Xbendl3, 3), bendREF_upper = round(Xbendu3, 3),
     bendSAMPLE_lower = round(XbendlT, 3), bendSAMPLE_upper = round(XbenduT, 3),
@@ -362,7 +362,7 @@ plot_f <- function(dat, TransFlag = F) { # sigmoid,det_sig,
     geom_vline(xintercept = c(XbendlT, XbenduT), col = "#C2173F", linetype = 2) +
     geom_vline(xintercept = c(XasymplS, XasympuS), col = "#4545BABB", linetype = 3) +
     geom_vline(xintercept = c(XasymplT, XasympuT), col = "#C2173FBB", linetype = 3) +
-    annotate("text", x = cs, y = a + (d - a) / 2, label = "0", size = 5) +
+    annotate("text", x = c, y = a + (d - a) / 2, label = "0", size = 5) +
     geom_abline(slope = slopeEC50, intercept = Intercept) +
     theme(legend.position = "none")
   Dat$p2 <- p2
@@ -419,16 +419,8 @@ plot_f <- function(dat, TransFlag = F) { # sigmoid,det_sig,
   cte <- Sum_u$coefficients["cs", 1] - Sum_u$coefficients["r", 1]
   dst <- Sum_u$coefficients["ds", 1]
   dte <- Sum_u$coefficients["dt", 1]
-  # } else {
-  #   ast <- det_sig[5]
-  #   ate <- det_sig[6]
-  #   bst <- det_sig[1]
-  #   bte <- det_sig[2]
-  #   cst <- det_sig[7]
-  #   cte <- det_sig[8]
-  #   dst <- det_sig[3]
-  #   dte <- det_sig[4]
-  # }
+ 
+  
   REFu <- ast + (dst - ast) / (1 + exp(bst * (cst - seq_x)))
   SAMPLEu <- ate + (dte - ate) / (1 + exp(bte * (cte - seq_x)))
   pl_df2 <- cbind(seq_x, SAMPLEu, REFu)
@@ -877,7 +869,7 @@ ANOVAlintests <- function(ro_new, circles, Lim, PureErrFlag) {
 #'
 #' PlotLinPLA_FUNC(circle, sigmoid, all_l2, pl_df, indS, indT)
 PlotLinPLA_FUNC <- function(circle, sigmoid, all_l2, pl_df, indS, indT) {
-  # browser()
+   #browser()
   mLin <- gsl_nls(readout ~ (intS + r) * isSample + intS * isRef + k * log_dose,
     data = circle,
     start = list(intS = 0, k = 1, r = 0),

dev/Doc_BioassayLinReport.Rmd

@@ -20,6 +20,8 @@ params:
   REP: NA
   REPlin: NA
   coeffsLin: NA
+  NoP: NA
+  Assay: NA
 author: "Author: `r params$author`"
 title: |
     | ![](logo.png){width=1in}
@@ -51,7 +53,7 @@ coeffsLin <- params$coeffsLin
 
 all_l <- REP$all_l
 circles <- REPlin$circles
-ANOVAXLS <- REP$ANOVAXLS
+#ANOVAXLS <- REP$ANOVAXLS
 SuModAB <- REPlin$SuModAB
 SuModABu <- REPlin$SuModABu
 LinTests <-  REPlin$LinTests
@@ -60,7 +62,8 @@ LinPotTab <- REPlin$LinPotTab
 
 XLdat2 <- REP$XLdat2
 
-
+LinTests1 <- LinTests[,1:3]
+ANOVAlin <- LinTests[,4:ncol(LinTests)]
 
 
 ```
@@ -68,8 +71,8 @@ XLdat2 <- REP$XLdat2
 
 # Introduction
 
-Bioassay potency estimation uses statistical methods to quantify the strength of a biological product or drug by comparing its response to that of a reference standard. Because biological responses are inherently variable, affected by assay conditions, cell systems or organisms, and measurement noise, the 4-parametric logistic regression is used to obtain reliable potency values. The variance for confidence interval calculation is coming from the regression procedure itself and is an excellent predictor for the variability of any future potency determinations.  
-USP<1034> recommends calculation of standard errors of ratios of the parameters using Fieller's theorem [Finney D.J. 1978] or using the "delta" method (for a discussion about the "delta" method see [Ver Hoef 2012]). However, the presented gradient approach using the differences on the log-scale is methematically more stable und thus preferable compared to any ratio approach ([Franz, V.H. 2007]).
+Bioassay potency estimation uses statistical methods to quantify the strength of a biological product or drug by comparing its response to that of a reference standard. Biological responses are inherently variable, affected by assay conditions, cell systems or organisms, and measurement noise. To control this variability, a linear regression approach is used to obtain reliable potency values. Three consecutive dilution steps showing the steepest slope are used for linear fitting.   
+USP<1034> recommends calculation of standard errors of ratios of the parameters using Fieller's theorem [Finney D.J. 1978] or using the "delta" method (for a discussion about the "delta" method see [Ver Hoef 2012]). The present analysis calculated the relative potency with the "delta" method. The formula of the relative potency is in the Appendix.
 
 # Raw data
 
@@ -81,7 +84,7 @@ kable(XLdat2, format = "markdown", caption= "Uploaded data (test and reference)
 
 ```
 
-All data used linerar regression is shown in table 2.
+The linerar regression is calculated on the readout listed in table 2.
 
 ```{r Circles, echo=FALSE, warning=FALSE, results='asis'}
 
@@ -135,21 +138,37 @@ plot_grid(XLplotLin)
 
 ```
 
-The ANOVA of the unconstrained model is listed in table 3.
 
-```{r anovaxls, echo=FALSE, warning=FALSE, results='asis'}
+The relative potency can be read from tbale 3.
 
-kable(ANOVAXLS, format = "markdown", caption= "ANOVA table unrestricted", digits=3)
+```{r LinPotTab, echo=FALSE, warning=FALSE, results='asis'}
+
+kable(LinPotTab, format = "markdown", caption= "Potency table", digits=3)
 
 
 ```
 
-The assay suitability tests are shown in table 4.
+
+
+
+The ANOVA of the unconstrained model is listed in table 4.
+
+```{r anovaxls, echo=FALSE, warning=FALSE, results='asis'}
+
+kable(ANOVAlin, format = "markdown", caption= "ANOVA table unrestricted", digits=3)
+
+RMSE <- sqrt(ANOVAlin[5,4])
+
+```
+
+The standard deviation of the model is `r RMSE`. 
+
+The assay suitability tests are shown in table 5.
 
 
 ```{r SST_ergebn, echo=FALSE, cache=FALSE, warning=FALSE, message=FALSE, tidy=TRUE}
 
-kable(LinTests, row.names = F, format = "markdown", caption="Assay suitability test results", digits=3)
+kable(LinTests1, row.names = F, format = "markdown", caption="Assay suitability test results", digits=3)
 
 
 ```
@@ -159,15 +178,15 @@ The estimate is the p-value of the test.
 F-tests on regression, significance of slopes, and preparation need to have a p-value <0.05 to pass.
 All other tests pass if p-value > 0.05.
 
-  0 ... test passed (for EQ tests: CI within limits);
-  
-  1 ... test failed (for EQ tests CI not within limits);
+  0 ... test passed;
   
+  1 ... test failed);
 
+(NOTE: F-tests are sensitive, when the residual variability of the method is small. On the other hand effects may not be detected if residual variability is high.)
 
 ## Fitting results
 
-The results of the linear fitting procedure for the restricted model is listed in table 5:
+The results of the linear fitting procedure for the restricted model is listed in table 6:
 
 ```{r SumCSSI, echo=FALSE, warning=FALSE, results='asis'}
 
@@ -178,9 +197,10 @@ kable(SuModAB, format = "markdown", caption= "Restricted linear regression (CSSI
 
 CSSI: common slope, separate intercept
 
-The results of the linear fitting procedure for the unrestricted model is listed in table 6.
+The results of the linear fitting procedure for the unrestricted model is listed in table 7.
 
-```{r SuSSSI, echo=FALSE, warning=FALSE, results='asis'}
+
+```{r SumSSSI, echo=FALSE, warning=FALSE, results='asis'}
 
 kable(SuModABu, format = "markdown", caption= "Restricted linear regression (SSSI)", digits=3, row.names = F)
 
@@ -197,7 +217,7 @@ SSSI: separate slope, separate intercept
 ## Potency of linear PLA
 
 $$
- rel Potency = \frac{I_{ref} - I_{test}{k}
+ rel Potency = \frac{I_{ref} - I_{test}}{k}
 $$
 where: I... intercept of reference or test
 k ... common slope

dev/Doc_BioassayReport.Rmd

@@ -93,6 +93,7 @@ kable(XLdat2, format = "markdown", caption= "Uploaded data (test and reference)
 
 ```{r Over_all, echo=FALSE, comment=NA, warning=NA, message=NA}
 
+# browser()
 potFlag <- 0
 if (pottab4plXL["test_result"][[1]][1]==1) potFlag <- 1
 AnalysisFlag <- FALSE

dev/app.R

@@ -548,6 +548,8 @@ server <- function(input, output, session) {
     )
   })
 
+  
+  #### UI wizard ----
   output$wizard <- renderUI({
     navbarPage(
       title = "Dilution setting",
@@ -557,14 +559,12 @@ server <- function(input, output, session) {
           sidebarPanel(
             width = 3,
             fluidRow(
-              column(
-                6,
-                numericInput(
-                  "Limits", p("limit to be >", bsButton("q4", label = "", icon = icon("info"), style = "primary", size = "extra-small")),
-                  bsPopover(id = "q4", title = "", content = "The calculated limits ...")
+              column(6,
+                box(title = "Upload multiple worksheets", status = "warning", solidHeader = T, width = 12, "Please upload your EXCEL file here",
+                  fileInput("MiFile", "", accept = ".xlsx"))
                 )
               )
-            )
+            
           ),
           mainPanel(
             tabsetPanel(
@@ -575,8 +575,7 @@ server <- function(input, output, session) {
                   title = "ANOVA table", status = "primary", solidHeader = T, width = 12,
                   tableOutput("Anovatab")
                 ),
-                column(
-                  4,
+                column(4,
                   h3("Confidence intervals"),
                   tableOutput("CIs"),
                   "The confidence interval table is interaactive for changes in: variability slider (%SD), potency of test-slider,
@@ -584,8 +583,7 @@ server <- function(input, output, session) {
                   tableOutput("optimalDils"),
                   selectInput(inputId = "scenario", label = "Select an 'optimal' scenario:", choices = c("scenario 2", "scenario 3", "scenario 6", "steep slope"))
                 ),
-                column(
-                  5,
+                column(5,
                   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),
@@ -595,8 +593,7 @@ server <- function(input, output, session) {
                   "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,
+                column(3,
                   h3("Bend points"),
                   tableOutput("bps"),
                   tableOutput("extremebps"),
@@ -670,6 +667,22 @@ server <- function(input, output, session) {
   })
 
 
+  observe({
+    if (!is.null(input$MiFile)) {
+      MinFile <- input$MiFile
+      ext <- tools::file_ext(MinFile$name)
+      file.rename(MinFile$datapath, paste(MinFile$datapath, ".xlsx", sep = ""))
+      t.filelocation <- gsub("\\\\", "/", paste(MinFile$datapath, ext, sep = "."))
+      sheets <- openxlsx::getSheetNames(t.filelocation)
+      dat <- lapply(sheets, openxlsx::read.xlsx, xlsxFile = t.filelocation)
+      names(dat) <- sheets
+      Dat$Mws <- dat
+      names(Dat$Mws) <- sheets
+      Dat$Msheets <- sheets
+      Dat$MFileName <- input$MiFile[["name"]]
+    }
+  })
+  
   #### process XLSX file ----
   observe({
     if (!is.null(input$iFile)) {
@@ -1339,7 +1352,7 @@ server <- function(input, output, session) {
     if (is.null(input$PureErr)) {
       return(NULL)
     }
-    if (!is.null(Dat$FITsFlag)) {
+    if (Dat$FITsFlag) {
       return(NULL)
     }
 
@@ -1974,10 +1987,12 @@ server <- function(input, output, session) {
 
   #### 4pl potency table XL ----
   observe({
+    
+   
     if (is.null(Dat$EXCEL)) {
       return(NULL)
     }
-    if (!is.null(Dat$FITsFlag)) {
+    if (Dat$FITsFlag) {
       return(NULL)
     }
     ro_new <- Dat$EXCEL
@@ -2025,7 +2040,7 @@ server <- function(input, output, session) {
     colnames(pottab4_) <- c("model", "potency", "lower95%CI", "upper95%CI", "relative_lower95%CI", "relative_upper95%CI", "test_result")
     row.names(pottab4_) <- NULL
     REP$pottab4plXL <- pottab4_[1:2, ]
-
+#browser()
     output$pottab4plXL <- DT::renderDataTable({
       dat <- datatable(pottab4_[1:2, ],
         rownames = F,
@@ -2055,8 +2070,101 @@ server <- function(input, output, session) {
 
   #### Dilutions Simulator ----
   output$plotfordilutions <- renderPlot({
-    tab <- sim2()
-    tab <- as.data.frame(tab)
+    if (!is.null(Dat$Mws))
+    AllXL <- Dat$Mws
+    AllSheets <- Dat$Msheets
+    
+    
+    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]]
+    potU_est <- FITs[[4]]
+    # unrestricted
+    SU_mu <- FITs[[2]]
+    URMcoeffs <- SU_mu$coefficients
+    
+    X <- seq(min(log(datWS23$log_dose)), max(log(datWS2$log_dose)), 0.1)
+    sigRef <- URMcoefs[1,1] + (URMcoefs1[4,1]-URMcoefs[1,1])/(1+exp(URMcoefs[2,1]*(URMcoefs[3,1]-X)))
+    sigTest1 <- URMcoefs[5,1] + (URMcoefs[8,1]-URMcoefs[5,1])/(1+exp(URMcoefs[6,1]*(URMcoefs[7,1]-X)))
+    
+    dfPlotsigRef <- data.frame(X=X, sigRef = sigRef, Prod = pdfInd)
+    dfPlotsigTest <- data.frame(X=X, sigTest = sigTest1, Prod = AllSheets[[N_WS]])
+    
+    if (!exists("SIGrefDF")) SIGrefDF <- dfPlotsigRef else SIGrefDF <- rbind(SIGrefDF, dfPlotsigRef)
+    if (!exists("SIGtestDF")) SIGtestDF <- dfPlotsigTest else SIGtestDF <- rbind(SIGtestDF,dfPlotsigTest)
+   
+    EC50TEST <- as.numeric(c(URMcoefsDF[,8]))
+    # EC50TEST <- EC50TEST[!EC50TEST %in% boxplot.stats(EC50TEST)$out]
+    EC50REF <- as.numeric(URMcoefsDF[,4])
+    # EC50REF <- EC50REF[!EC50REF %in% boxplot.stats(EC50REF)$out]
+    UasREF <- as.numeric(URMcoefsDF[,5])
+    # UasREF <- UasREF[!UasREF %in% boxplot.stats(UasREF)$out]
+    LasREF <- as.numeric(URMcoefsDF[,2])
+    # LasREF <- LasREF[!LasREF %in% boxplot.stats(LasREF)$out]
+    # 
+    # Dat$URMcoefsDF <- URMcoefsDF
+    # Dat$RestrM <- RestrM
+    # Dat$CalcPot <- CalcPot
+    # 
+    #### sigmoid plots ----
+    Slope <- as.numeric(URMcoefsDF[1,3])
+    # if (Slope > 0) {
+    #   x_UA <- max(X); x_LA <- min(X)
+    # } else { x_UA <- min(X); x_LA <- max(X) }
+    # 
+    # p1 <- ggplot(SIGrefDF, aes(x_X, y=sigRef, col=as.factor(Prod))) +
+    #   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) +
+    #   xlab("dilutions") +
+    #   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 })
+    # 
+    # PLOTS$sigPlotREF <- p1
+    # 
+    # p2 <- ggplot(SIGtestDF, aes(x_X, y=sigTest, col=as.factor(Prod))) +
+    #   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 = EC50TEST, alpha = 0.2) +
+    #   xlab("dilutions") +
+    #   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({ p2 })
+    # 
+    # PLOTS$sigPlotTEST <- p2
+    
+    
     dils <- tab$log_dose
     min_y <- min(tab[, 1:3])
     max_y <- max(tab[, 1:3])
@@ -2077,6 +2185,9 @@ server <- function(input, output, session) {
       dils2 <- dils_avsc + av
       dilfactors <- 1 / exp(dils2 - lag(dils2))
     }
+    } #for N_WS
+    
+    
     
     Dat$newDils <- dils2
 
@@ -2538,6 +2649,8 @@ server <- function(input, output, session) {
         params = list(
           FileName = Dat$FileName,
           author = Dat$Author,
+          NoP = Dat$NoP,
+          Assay = Dat$Assay,
           REP = REP,
           REPlin = REPlin,
           coeffsLin = Dat$coeffs_UN