added the positive control normalization. Furthermore, refined the se…

OpenStats/R/Backend_History.R

@@ -158,6 +158,18 @@ eval_entry_V1_2 <- function(entry, DataModelState,
       res$eval(ResultsState, entry[["Result name"]])
       get_result(ResultsState)
     },
+    PrimaryAssay = {
+      res <- primary_assay_V1_2$new(
+        DataModelState$df,
+        DataModelState$formula,
+        entry[["Negative control name"]],
+        entry[["Positive control name"]],
+        backend_communicator_V1_2
+      )
+      res$validate()
+      res$eval(ResultsState)
+      get_result(ResultsState)
+    },
     TTest= {
       res <- t_test_V1_2$new(
         DataModelState$df,

OpenStats/R/Backend_V1_2_Engine.R

@@ -80,7 +80,7 @@ bg_process_V1_2 <- R6::R6Class("bg_process_V1_2",
           if (ResultsState$bgp$warnings != "") {
             ResultsState$bgp$com$print_warn(ResultsState$bgp$warnings)
           }
-          e <- try(env_utils_V1_2$check_rls(ResultsState$all_data, res))
+          e <- try(env_utils_V1_2$check_rls(ResultsState$all_data, res), silent = TRUE)
           if (inherits(e, "try-error")) {
             self$com$print_err(conditionMessage(e))
             self$enable()
@@ -147,7 +147,7 @@ bg_process_V1_2 <- R6::R6Class("bg_process_V1_2",
           )
         }
       } else {
-        self$com$print_err(res)
+        self$com$print_err(attributes(res)$condition$message)
       }
     },
 
@@ -1424,6 +1424,96 @@ dose_response_V1_2 <- R6::R6Class(
 
   )
 )
+primary_assay_V1_2 <- R6::R6Class(
+  "primary_assay_V1_2",
+  public = list(
+    df = NULL,
+    formula = NULL,
+    name_column = NULL,
+    neg_control_name = NULL,
+    pos_control_name = NULL,
+    com = NULL,
+
+    initialize = function(df, formula, neg_control_name, pos_control_name, com = communicator_V1_2) {
+      self$df <- df
+      self$formula <- formula
+      indep <- try({
+        f <- as.character(formula@formula)
+        f[[3L]]
+      }, silent = TRUE)
+      if (inherits(indep, "try-error")) stop("The independent variable is not usable")
+      self$name_column <- indep
+      self$neg_control_name <- neg_control_name
+      self$pos_control_name <- pos_control_name
+      self$com <- com$new()
+    },
+
+    validate = function() {},
+
+    eval = function(ResultsState) {
+      withCallingHandlers(
+        expr = {
+          new_name <- paste0(ResultsState$counter + 1, " Primary Assay")
+          promise_history_entry <- self$create_history(new_name)
+
+          ResultsState$bgp$start(
+            fun = function(df, formula, name_column, neg_control_name, pos_control_name) {
+              indep <- as.character(formula@formula)[3]
+              dep <- as.character(formula@formula)[2]
+              if (pos_control_name != "") {
+                mean_pos <- mean(df[df[, indep] == pos_control_name, dep], na.rm = TRUE)
+                df[, dep] <- df[, dep] - mean_pos
+                df <- df[df[, indep] != pos_control_name, , drop = FALSE] # Remove pos control
+              } else {
+                warning("No positive control is specified the respective normalization will not be conducted")
+              }
+              if (neg_control_name != "") {
+                mean_neg <- mean(df[df[, indep] == neg_control_name, dep], na.rm = TRUE)
+                df[, dep] <- (df[, dep] / mean_neg) * 100
+              } else {
+                stop("You have to define a name for the negative control")
+              }
+              fit <- lm(formula@formula, data = df)
+              emm <- emmeans::emmeans(fit, name_column)
+              res <- emmeans::contrast(emm, method = "trt.vs.ctrl", ref = neg_control_name, adjust = "holm")
+              res <- broom::tidy(res) |> as.data.frame()
+              res <- res[, c(2, 4, 8)]
+              names(res) <- c("name", "Standard Value", "adj. p value")
+              res[["Standard Units"]] <- "%"
+              res <- res[, c(1, 2, 4, 3)]
+              pattern <- paste0("\\s*-\\s*", neg_control_name, "$")
+              res$name <- gsub(pattern, "", res$name)
+              res
+            },
+            args = list(
+              df = self$df, formula = self$formula,
+              name_column = self$name_column,
+              neg_control_name = self$neg_control_name, pos_control_name = self$pos_control_name
+            ),
+            promise_result_name = new_name,
+            promise_history_entry = promise_history_entry,
+            in_background = FALSE, ResultsState = ResultsState
+          )
+
+        },
+        warning = function(warn) {
+          self$com$print_warn(warn$message)
+          invokeRestart("muffleWarning")
+        }
+      )
+    },
+
+    create_history = function(new_name) {
+      list(
+        type = "PrimaryAssay",
+        formula = deparse(self$formula@formula),
+        "Negative control name" = self$neg_control_name,
+        "Positive control name" = self$pos_control_name,
+        "Result name" = new_name
+      )
+    }
+  )
+)
 
 t_test_V1_2 <- R6::R6Class(
   "t_test_V1_2",

OpenStats/R/Getter.R

@@ -159,6 +159,11 @@ get_dose_response <- function(version = VERSION) {
     `1.2` = dose_response_V1_2
   )[[version]]
 }
+get_primary_assay <- function(version = VERSION) {
+  list(
+    `1.2` = primary_assay_V1_2
+  )[[version]]
+}
 get_ttest <- function(version = VERSION) {
   list(
     `1.2` = t_test_V1_2

OpenStats/R/Server_CheckFunctions.R

@@ -63,6 +63,36 @@ check_dose_response <- function(DataModelState) {
 }
 check_correlation <- check_dose_response
 
+check_primary_assay <- function(DataModelState) {
+  if (is.null(DataModelState$df)) return("No data is available")
+  if (!is.data.frame(DataModelState$df)) return("Dataset seems to be malformed")
+
+  formula <- DataModelState$formula
+  if (is.null(formula)) return("You have to define a model in the formula editor")
+  if (!inherits(formula, "LinearFormula")) return("Only linear models are supported.")
+
+  f <- formula@formula
+  response <- all.vars(f[[2]])
+  predictor <- all.vars(f[[3]])
+
+  if (length(response) != 1 || length(predictor) != 1) {
+    return("The formula must have one predictor and one response.")
+  }
+
+  df <- DataModelState$df
+  if (!(response %in% names(df))) {
+    return(sprintf("The response variable '%s' is not a column of the active dataset", response))
+  }
+  if (!(predictor %in% names(df))) {
+    return(sprintf("The predictor variable '%s' is not a column of the active dataset", predictor))
+  }
+  if (!is.numeric(df[[response]])) {
+    return(sprintf("The response variable '%s' must be numeric.", response))
+  }
+
+  return(NULL)
+}
+
 check_assumptions <- function(DataModelState) {
   if (is.null(DataModelState$df)) return("No data is available")
   if (!is.data.frame(DataModelState$df)) return("Dataset seems to be malformed")

OpenStats/R/Server_DoseResponse.R

@@ -71,6 +71,55 @@ DoseResponseServer <- function(id, DataModelState, ResultsState) {
       )
     })
 
+    # Render control group
+    output[["primaryAssayUI"]] <- shiny::renderUI({
+      shiny::req(!is.null(DataModelState$df))
+      shiny::req(is.data.frame(DataModelState$df))
+      shiny::req(inherits(DataModelState$formula, "LinearFormula"))
+
+      message <- try(check_primary_assay(DataModelState), silent = TRUE)
+      if (!is.null(message)) {
+        return(info_div(message))
+      }
+
+      indep <- try({
+        f <- as.character(DataModelState$formula@formula)
+        f[[3L]]
+      }, silent = TRUE)
+      if (inherits(indep, "try-error")) return()
+
+      choices <- unique(DataModelState$df[[indep]])
+      htmltools::div(
+        htmltools::h4("Primary Assay"),
+
+        shiny::selectizeInput(
+          inputId = paste0("DOSERESPONSE-neg_control_name"),
+          label = "Name of the negative control",
+          selected = choices[1],
+          choices = choices,
+          options = list(
+            placeholder = 'Type to search...',
+            maxOptions = 1000
+          )
+        ),
+
+        shiny::selectizeInput(
+          inputId = paste0("DOSERESPONSE-pos_control_name"),
+          label = "Name of the positive control",
+          selected = choices[1],
+          choices = choices,
+          options = list(
+            placeholder = 'Type to search...',
+            maxOptions = 1000
+          )
+        ),
+
+        shiny::actionButton("DOSERESPONSE-primary_assay", "Run primary assay"),
+        message,
+        class = "boxed-output"
+      )
+    })
+
     check_dr <- function() {
       print_req(is.data.frame(DataModelState$df), "The dataset is missing")
       shiny::req(input$substanceNames)
@@ -93,6 +142,16 @@ DoseResponseServer <- function(id, DataModelState, ResultsState) {
       e <- try(run_dr(df, new_name))
     })
 
+    shiny::observeEvent(input$primary_assay, {
+      try({
+        pa <- get_primary_assay()$new(
+          DataModelState$df, DataModelState$formula,
+          input$neg_control_name, input$pos_control_name
+        )
+        pa$eval(ResultsState)
+      })
+    })
+
   })
 
 }

OpenStats/R/UI_DoseResponse.R

@@ -1,12 +1,19 @@
 DoseResponseSidebarUI <- function(id) {
   shiny::tabPanel(
     "Dose Response analysis",
-    shiny::uiOutput(shiny::NS(id, "substanceNamesUI")),
-    shiny::uiOutput(shiny::NS(id, "unitNamesUI")),
-    shiny::sliderInput(shiny::NS(id, "ic_percentage"), "Percentage if IC",
-      min = 1, max = 99, value = 50
+    htmltools::div(
+      htmltools::h4("Dose Response analysis"),
+      shiny::uiOutput(shiny::NS(id, "substanceNamesUI")),
+      shiny::uiOutput(shiny::NS(id, "unitNamesUI")),
+      shiny::sliderInput(shiny::NS(id, "ic_percentage"), "Percentage if IC",
+        min = 1, max = 99, value = 50
+      ),
+      shiny::uiOutput(shiny::NS(id, "DoseResponseUI")),
+      class = "boxed-output"
     ),
-    shiny::uiOutput(shiny::NS(id, "DoseResponseUI"))
+    htmltools::br(),
+    htmltools::br(),
+    shiny::uiOutput(shiny::NS(id, "primaryAssayUI"))
   )
 }
 

development/primary_assay.R

@@ -0,0 +1,9 @@
+df <- read.csv("./test_data/DoseResponseData.csv", sep = ",")
+head(df)
+name_col <- "substance"
+neg_control_name <- "neg"
+dep <- "abs"
+fit <- lm(abs ~ substance, data = df)
+emm <- emmeans::emmeans(fit, "substance")
+res <- emmeans::contrast(emm, method = "trt.vs.ctrl", ref = neg_control_name)
+broom::tidy(res)