R/offsets.R

In MoseleyBioinformaticsLab/ScanCentricPeakCharacterization: Functionality for Characterizing Peaks in Mass Spectrometry in a Scan-Centric Manner

default_offset_fit_function = function(x, y){
  loess_frame = data.frame(x = x, y = y)
  loess_frame = loess_frame[!(is.na(x) | is.na(y)), ]
  loess_fit = stats::loess(y ~ x, data = loess_frame, span = 1.5, control = stats::loess.control(surface = "direct"))
  loess_fit
}


#' offset predict function
#'
#' The offset predictor using loess
#'
#' @param model the model to use
#' @param x the new values
#'
#' @export
#'
#' @return numeric
default_offset_predict_function = function(model, x){
  stats:::predict.loess(model, newdata = x)
}

correctly_round_numbers = function(number_of_things, fraction){
  assertthat::assert_that(number_of_things > 0)
  assertthat::assert_that(fraction > 0)
  if (fraction < 1) {
    floor_value = floor(number_of_things * fraction)
  } else {
    floor_value = fraction
  }
  if (floor_value < 1) {
    floor_value = 1
  }
  floor_value
}

#' correct offsets
#'
#' Given a MasterPeakList object and the MultiScansPeakList that generated it,
#' correct the m/z values using offset predictions
#'
#' @param master_peak_list the MasterPeakList object of correspondent peaks
#' @param multi_scan_peaklist the MultiScansPeakList to be corrected
#' @param min_scan what is the minimum number of scans a peak should be in to be used for correction.
#'
#' @return list
#' @export
default_correct_offset_function = function(master_peak_list, multi_scan_peaklist, min_scan = 0.1){

  n_col = ncol(master_peak_list$scan_mz)
  n_min_scan = correctly_round_numbers(n_col, min_scan)
  correspond_peaks = master_peak_list$count_notna() >= n_min_scan

  scan_indices = master_peak_list$scan_indices

  multi_scan_peaklist_corrected = multi_scan_peaklist$clone(deep = TRUE)

  offset_models = vector(mode = "list", length = length(scan_indices))

  for (iscan in seq(1, n_col)) {
    scan_offset = master_peak_list$scan_mz[correspond_peaks, iscan] -
      master_peak_list$master[correspond_peaks]
    offset_model = master_peak_list$offset_fit_function(master_peak_list$master[correspond_peaks], scan_offset)

    use_scan = scan_indices[iscan]
    tmp_peaks = multi_scan_peaklist_corrected$peak_list_by_scans[[use_scan]]$peak_list
    use_mz = tmp_peaks$ObservedMZ
    offset_predict = master_peak_list$offset_predict_function(offset_model, use_mz)
    corrected_mz = use_mz - offset_predict
    tmp_peaks$ObservedMZ = corrected_mz
    multi_scan_peaklist_corrected$peak_list_by_scans[[use_scan]]$peak_list = tmp_peaks
    offset_model$scan_index = scan_indices[iscan]
    offset_models[[iscan]] = offset_model
  }

  list(multi_scan_peaklist = multi_scan_peaklist_corrected,
       models = offset_models)
}

#' create loess data.frame
#'
#' Given a loess model, creates a data.frame suitable for plotting via ggplot2
#'
#' @param loess_model the model object generated by loess
#'
#' @return data.frame
#' @export
loess_to_df = function(loess_model){
  out_frame = data.frame(x = rep(loess_model$x[,1], 2),
                          y = c(loess_model$y, loess_model$fitted),
                          which = rep(c("original", "fitted"), each = loess_model$n),
                          stringsAsFactors = FALSE)
  if (!is.null(loess_model$scan_index)) {
    out_frame$scan_index = loess_model$scan_index
  }
  out_frame
}