R/calculateCurveFit.R

In CeMOS-Mannheim/MALDIcellassay: Automated MALDI Cell Assays Using Dose-Response Curve Fitting

#' Calculate the fit for a dose-response curve
#'
#' @param intmat  Intensity matrix as generated by
#'                \code{MALDIquant::intensityMatrix()}
#'                with rownames as the respective concentrations of the spectra.
#' @param idx     Numeric vector of the mz indices to perform the fit.
#' @param verbose Logical, print logs to console.
#' @param ...     Additional arguments passed to \code{nplr::nplr()}.
#'
#' @return
#' List of curve fits.
#' @export
#' @importFrom nplr nplr convertToProp
#' @importFrom dplyr mutate arrange %>% .data
#' @importFrom tibble as_tibble
#' @importFrom stats rnorm
#' 
#' @examples
#' data(Blank2022intmat)
#' 
#' # for faster runtime we let it run on 5 peaks only
#' fits <- calculateCurveFit(Blank2022intmat, idx = 1:5)
calculateCurveFit <- function(intmat, idx, verbose = TRUE, ...) {
  if (length(rownames(intmat)) == 0) {
    stop("intmat must have concentrations as rownames!\n")
  }

  current_res <- vector("list", length = length(idx))
  names(current_res) <- colnames(intmat[, idx])
  for (j in 1:length(idx)) {
    if(verbose) {
      cat("fitting", round(as.numeric(colnames(intmat)[j]), 3),
          paste0("(", j, "/", length(idx), ")\n"))  
    }
    
    df <- intmat[, idx[j]] %>%
      as_tibble() %>%
      mutate(conc = rownames(intmat)) %>%
      mutate(conc = as.numeric(.data$conc)) %>%
      arrange(.data$conc)

    # transform concentrations to log10 and handle zero concentrations.
    concLog <- transformConc2Log(df$conc)

    df <- df %>%
      mutate(concLog = concLog)
    resp <- df$value
      #convertToProp(y = df$value)

    model <- tryCatch(expr = {
      suppressWarnings(
        nplr(x = concLog, y = resp, useLog = FALSE, npars = 4, silent = TRUE, ...)
      )
    }, error = function(cond) {
      warning("m/z ", round(as.numeric(colnames(intmat)[j]), 3),
          " failed. Re-trying with npar='all' and additional noise (mean=0, sd=1e-4).\n")
      return(
        tryCatch(expr = {
          resp <- resp + rnorm(length(resp),
                               mean = 0,
                               sd = 1e-4)
        suppressWarnings(
          nplr(x = concLog, y = resp, useLog = FALSE, npars = "all", silent = TRUE, ...)
        )
      }, error = function(cond) {
        warning("m/z ", round(as.numeric(colnames(intmat)[j]), 3),
                " failed again. Forcing flat line.\n")
        resp <- rnorm(n = length(resp), sd = 0.5)
        suppressWarnings( {
          nplr(x = concLog, y = resp, useLog = FALSE, npars = 2, silent = TRUE, ...)
        }
        )
      })
      )
    })

    current_res[[j]] <- list(
      model = model,
      df = df
    )

    res_list <- current_res
  }
  return(res_list)
}