Improve default peak detection default parameters

MassSpecWavelet::peakDetectionCWT() has a lot of parameters.

Some of those parameters depend on the chosen scales, and/or the expected peak widths or signal lengths.

This commit defines defaultPeakDetectionCWTParams() which is a function from the GCIMS package that sets reasonable default values for those parameters, suitable
for GCIMS data.

Author: zeehio

DESCRIPTION

@@ -1,9 +1,9 @@
 Package: GCIMS
 Type: Package
 Title: Pre-processing for GC-IMS data
-Version: 0.1.1
+Version: 0.2.0
 Encoding: UTF-8
-Date: 2024-04-29
+Date: 2025-05-03
 Authors@R: c(
   person(
     given = "Sergio",
@@ -118,7 +118,7 @@ Suggests:
     viridisLite
 biocViews: Software, Preprocessing, Visualization, Classification,
  Cheminformatics, Metabolomics, DataImport
-RoxygenNote: 7.3.1
+RoxygenNote: 7.3.2
 VignetteBuilder: knitr
 Config/testthat/edition: 3
 Roxygen: list(markdown = TRUE)

NAMESPACE

@@ -30,6 +30,7 @@ export(clusterPeaks)
 export(create_annotations_table)
 export(cubic_root_trans)
 export(decimate)
+export(defaultPeakDetectionCWTParams)
 export(description)
 export(download_three_ketones_dataset)
 export(dtime)

NEWS.md

@@ -1,3 +1,9 @@
+# GCIMS 0.2.0 (2025-05-03)
+
+- Deprecate `peakDetectionCWTParams = list(exclude0scaleAmpThresh = TRUE)`. Please
+  remove that argument, since the new default is better. For further details on
+  the new default please check the vignette or `?GCIMS::defaultPeakDetectionCWTParams`.
+
 # GCIMS 0.1.1 (2024-04-29)
 
 - Fix undefined variable on RIP saturation detection (#21)

R/findPeaks.R

@@ -67,11 +67,10 @@ interleave_peaks_with_zeros <- function(peak_idx, zero_idx, signal_idx) {
   out
 }
 
-detect_peaks_and_zeros_one_signal <- function(x, x_zeros, idx, save_debugging, prep_wav, peakDetectionCWTParams) {
+detect_peaks_and_zeros_one_signal <- function(x, x_zeros, idx, save_debugging, peakDetectionCWTParams) {
   peakdet <- rlang::exec(
     MassSpecWavelet::peakDetectionCWT,
     ms = x,
-    scales = prep_wav,
     !!!peakDetectionCWTParams
   )
   if (is.null(x_zeros)) {
@@ -97,8 +96,11 @@ detect_peaks_and_zeros_one_signal <- function(x, x_zeros, idx, save_debugging, p
   )
 }
 
-detect_peaks_and_zeros <- function(xmat, rowwise, scales, peakDetectionCWTParams = list(exclude0scaleAmpThresh = TRUE),
-                                   signals_to_save_extra = NULL, xmat_zeros = NULL) {
+detect_peaks_and_zeros <- function(
+    xmat, rowwise,
+    peakDetectionCWTParams,
+    signals_to_save_extra = NULL, xmat_zeros = NULL
+) {
   if (rowwise) {
     direction <- "rowwise"
     signal_length <- ncol(xmat)
@@ -129,7 +131,10 @@ detect_peaks_and_zeros <- function(xmat, rowwise, scales, peakDetectionCWTParams
     }
 
   }
-  prep_wav <- MassSpecWavelet::prepareWavelets(signal_length, scales = scales)
+  peakDetectionCWTParams$scales <- MassSpecWavelet::prepareWavelets(
+    signal_length,
+    scales = peakDetectionCWTParams$scales
+  )
   save_debugging <- rep(FALSE, num_signals)
   if (!is.null(signals_to_save_extra)) {
     save_debugging[signals_to_save_extra] <- TRUE
@@ -145,22 +150,9 @@ detect_peaks_and_zeros <- function(xmat, rowwise, scales, peakDetectionCWTParams
       save_debugging = save_debugging
     ),
     detect_peaks_and_zeros_one_signal,
-    prep_wav = prep_wav,
     peakDetectionCWTParams = peakDetectionCWTParams,
     .progress = ifelse(show_progress_bar(), glue("Detecting peaks and zeros ({direction})"), FALSE)
   )
-  # peaks_and_zeros_and_extra <- mapply(
-  #   FUN = detect_peaks_and_zeros_one_signal,
-  #   x = signals,
-  #   x_zeros = signals_zeros,
-  #   idx = seq_len(num_signals),
-  #   save_debugging = save_debugging,
-  #   MoreArgs = list(
-  #     prep_wav = prep_wav,
-  #     peakDetectionCWTParams = peakDetectionCWTParams
-  #   ),
-  #   SIMPLIFY = FALSE
-  # )
 
   peaks_and_zeros <- purrr::map_dfr(peaks_and_zeros_and_extra, "peaks_and_zeros")
   stopifnot(ncol(peaks_and_zeros) == 4L)
@@ -391,12 +383,10 @@ prep_wav_from_peakwidths <- function(axis, peakwidth_min, peakwidth_max) {
 #' 5. We merge similar ROIs using a threshold on the intersection over union
 #' 6. Get some ROI metrics and return.
 #'
-#' For the MassSpecWavelet-based peak detection, the `scales` are computed based on
-#' the requested peak widths. Besides, the scales, further tuning beyond the
-#' MassSpecWavelet defaults is possible through the `*_peakDetectionCWTParams` argument.
-#' By default, the only change we introduce is the `exclude0scaleAmpThresh = TRUE`
-#' which is a reasonable peak detection setting not enabled in MassSpecWavelet
-#' for backwards compatibility reasons.
+#' For the MassSpecWavelet-based peak detection, custom parameters are calculated
+#' with our [GCIMS::defaultPeakDetectionCWTParams()], feel free to call that
+#' function and customize the `SNR.Th` or any other parameter as shown in the
+#' vignette.
 #'
 #' @keywords internal
 #' @param drift_time The drift time vector
@@ -418,15 +408,50 @@ findPeaksImpl <- function(
     rt_length_s = 8,
     dt_peakwidth_range_ms = c(0.15, 0.4),
     rt_peakwidth_range_s = c(10, 25),
-    dt_peakDetectionCWTParams = list(exclude0scaleAmpThresh = TRUE),
-    rt_peakDetectionCWTParams = list(exclude0scaleAmpThresh = TRUE),
+    dt_peakDetectionCWTParams = defaultPeakDetectionCWTParams(x = drift_time, peakwidth_range_xunits = dt_peakwidth_range_ms),
+    rt_peakDetectionCWTParams = defaultPeakDetectionCWTParams(x = retention_time, peakwidth_range_xunits = rt_peakwidth_range_s),
     dt_extension_factor = 0,
     rt_extension_factor = 0,
     exclude_rip = FALSE,
     iou_overlap_threshold = 0.2,
     debug_idx = list(dt = NULL, rt = NULL)
 ) {
 
+  #### Handle older dt_peakDetectionCWTParams. Remove in 2026-01-01.
+  old_peakDetectionCWTParams <- list(exclude0scaleAmpThresh = TRUE)
+  if (identical(dt_peakDetectionCWTParams, old_peakDetectionCWTParams)) {
+    cli::cli_warn(c(
+      "Deprecated dt_peakDetectionCWTParams value",
+      "!" = "Passing the argument {.code dt_peakDetectionCWTParams = list(exclude0scaleAmpThresh = TRUE)} is deprecated since GCIMS v0.2.1.",
+      "i" = "Please remove the {.arg dt_peakDetectionCWTParams} argument. The new default value is better."
+    ))
+    dt_peakDetectionCWTParams = defaultPeakDetectionCWTParams(
+      x=drift_time,
+      peakwidth_range_xunits = dt_peakwidth_range_ms
+    )
+  }
+  if (identical(rt_peakDetectionCWTParams, old_peakDetectionCWTParams)) {
+    cli::cli_warn(c(
+      "Deprecated rt_peakDetectionCWTParams value",
+      "!" = "Passing the argument {.code rt_peakDetectionCWTParams = list(exclude0scaleAmpThresh = TRUE)} is deprecated since GCIMS v0.2.1.",
+      "i" = "Please remove the {.arg rt_peakDetectionCWTParams} argument. The new default value is better."
+    ))
+    rt_peakDetectionCWTParams = defaultPeakDetectionCWTParams(
+      x=retention_time,
+      peakwidth_range_xunits = rt_peakwidth_range_s
+    )
+  }
+  #### Code handling older dt_peakDetectionCWTParams ends here.
+
+    if (!is.null(dt_peakwidth_range_ms) && !identical(attr(dt_peakDetectionCWTParams, "peakwidth_range_xunits"), dt_peakwidth_range_ms)) {
+      stop("Please provide a new dt_peakDetectionCWTParams created with the given dt_peakwidth_range_ms")
+    }
+
+    if (!is.null(rt_peakwidth_range_s) && !identical(attr(rt_peakDetectionCWTParams, "peakwidth_range_xunits"), rt_peakwidth_range_s)) {
+      stop("Please provide a new rt_peakDetectionCWTParams created with the given rt_peakwidth_range_s")
+    }
+
+
     debug_info <- list()
     if (!is.null(debug_idx$rt) && length(debug_idx$rt) > 0) {
       debug_info$rt <- list()
@@ -459,31 +484,9 @@ findPeaksImpl <- function(
     rm(deriv2)
 
     # 5. Peaks and Zero-crossings
-    scales_rt <- prep_wav_from_peakwidths(
-      retention_time,
-      peakwidth_min = rt_peakwidth_range_s[1L],
-      peakwidth_max = rt_peakwidth_range_s[2L]
-    )
-    scales_dt <- prep_wav_from_peakwidths(
-      drift_time,
-      peakwidth_min = dt_peakwidth_range_ms[1L],
-      peakwidth_max = dt_peakwidth_range_ms[2L]
-    )
-
-    if (verbose) {
-      cli_inform(
-        message = c(
-          "Using the following scales",
-          "i" = "Retention time scales: {glue::glue_collapse(scales_rt$scales, sep = ', ')}",
-          "i" = "Drift time scales: {glue::glue_collapse(scales_dt$scales, sep = ', ')}"
-        )
-      )
-    }
-
     peaks_and_zeros_drt_extra <- detect_peaks_and_zeros(
       xmat = int_mat,
       rowwise = TRUE,
-      scales = scales_rt,
       peakDetectionCWTParams = rt_peakDetectionCWTParams,
       signals_to_save_extra = debug_idx$rt,
       xmat_zeros = drt
@@ -497,7 +500,6 @@ findPeaksImpl <- function(
     peaks_and_zeros_ddt_extra <- detect_peaks_and_zeros(
       xmat = int_mat,
       rowwise = FALSE,
-      scales = scales_dt,
       peakDetectionCWTParams = dt_peakDetectionCWTParams,
       signals_to_save_extra = debug_idx$dt,
       xmat_zeros = ddt
@@ -567,6 +569,77 @@ findPeaksImpl <- function(
 }
 
 
+#' Get the default peakDetectionCWT parameters for GCIMS peak detection
+#'
+#' @description
+#'
+#' Our selection of defaults for [peakDetectionCWT()] consist of the following:
+#'
+#' \itemize{
+#' \item{`exclude0scaleAmpThresh=TRUE`: Does not include the zeroth scale in the default amplitude threshold.}
+#' \item{`nearbyPeak=TRUE`: This is already a default, but since we customize nearbyWinSize we prefer to be explicit.}
+#' \item{`nearbyWinSize`: This is set to the maximum desired peak width in points (instead of the arbitrary 150 from [MassSpecWavelet::identifyMajorPeaks()])}
+#' \item{`ridgeLength`: We set this to the second calculated scale, that would correspond to a peak with shorter width in our scale preparation.}
+#' \item{`excludeBoundariesSize`: Instead of an arbitrary number, we exclude half of the peak minimum width}
+#' \item{`SNR.Th`: We keep the default (3) but we explicitly specify it because we consider this a sensitive parameter}
+#' \item{`minNoiseLevel`: [MassSpecWavelet::identifyMajorPeaks()] uses a default of 0.001 and we keep that default}
+#' \item{`scales`: We provide here the [MassSpecWavelet::prepareWavelets()]
+#'       output, created with scales representative of the peak widths. You can
+#'       provide just a numeric vector with scales, but if you provide the prepared
+#'       wavelets you get some extra performance.}
+#'
+#' }
+#'
+#' @keywords internal
+#' @param x The x axis, to determine units
+#' @param peakwidth_range_xunits A vector of length 2 with the minimum and maximum peak width.
+#'
+#' @returns
+#' A list with the default parameters we will use with [MassSpecWavelet::peakDetectionCWT()].
+#' @export
+#'
+#' @examples
+#' chromat1 <- GCIMSChromatogram(
+#'   retention_time = seq(from = 0, to = 10, length.out = 200),
+#'   intensity = 1:200
+#' )
+#'
+#' peakDetectionCWTParams <- defaultPeakDetectionCWTParams(
+#'   x = rtime(chromat1),
+#'   peakwidth_range_xunits = c(5, 25)
+#' )
+#' peakDetectionCWTParams$SNR.Th <- 2.0
+#' peakDetectionCWTParams$scales <- MassSpecWavelet::prepareWavelets(
+#'   x = length(rtime(chromat1)),
+#'   scales = c(1, seq(2, 30, 2), seq(32, 64, 4))
+#' )
+#' peakDetectionCWTParams
+defaultPeakDetectionCWTParams <- function(x, peakwidth_range_xunits) {
+  step <- x[2] - x[1]
+  peakwidth_min_pts <- units_to_points(peakwidth_range_xunits[1L], step)
+  peakwidth_max_pts <- units_to_points(peakwidth_range_xunits[2L], step)
+
+  prep_wavelets <- prep_wav_from_peakwidths(
+    x,
+    peakwidth_min = peakwidth_range_xunits[1L],
+    peakwidth_max = peakwidth_range_xunits[2L]
+  )
+
+  peakDetectionCWTParams <- list(
+    exclude0scaleAmpThresh = TRUE,
+    nearbyPeak=TRUE, # This is already the default
+    nearbyWinSize=peakwidth_max_pts,
+    ridgeLength=prep_wavelets$scales[2],
+    excludeBoundariesSize=max(1, floor(peakwidth_min_pts/2)),
+    SNR.Th=3,
+    minNoiseLevel = 0.001,
+    scales = prep_wavelets
+  )
+  attr(peakDetectionCWTParams, "peakwidth_range_xunits") <- peakwidth_range_xunits
+  peakDetectionCWTParams
+}
+
+
 #---------------#
 #   FUNCTIONS   #
 #---------------#

R/findPeaksImpl1D.R

@@ -92,20 +92,20 @@ rois_to_peaklist_1D <- function(ROIs, x, y, deriv2) {
 #' 4. We merge similar ROIs using a threshold on the 1D-intersection over union
 #' 5. Get some ROI metrics and return.
 #'
-#' For the MassSpecWavelet-based peak detection, the `scales` are computed based on
-#' the requested peak widths. Besides, the scales, further tuning beyond the
-#' MassSpecWavelet defaults is possible through the `peakDetectionCWTParams` argument.
-#' By default, the only change we introduce is the `exclude0scaleAmpThresh = TRUE`
-#' which is a reasonable peak detection setting not enabled in MassSpecWavelet
-#' for backwards compatibility reasons.
+#' For the MassSpecWavelet-based peak detection, the `scales` and some other parameters
+#' are customized based on the requested peak widths. Still, if you want to further
+#' customize other peak detection sensitivity parameters, you can start calling
+#' [defaultPeakDetectionCWTParams()] and modify its output. Then pass that output
+#' to this function using `peakDetectionCWTParams=`.
+#'
 #'
 #' @keywords internal
 #' @param x The x axis, to determine units
 #' @param y The vector where to find peaks, of dimensions `length(x)`
 #' @param verbose If `TRUE` information will be printed on screen
 #' @param length_in_xunits Length of the filter used to compute the second derivative. See details.
 #' @param peakwidth_range_xunits A vector of length 2 with the minimum and maximum peak width. See details.
-#' @param peakDetectionCWTParams Additional parameters to [MassSpecWavelet::peakDetectionCWT()]. See details.
+#' @param peakDetectionCWTParams Parameters to [MassSpecWavelet::peakDetectionCWT()]. See details.
 #' @param extension_factor A number to extend the ROIs beyond their default size
 #' @param iou_overlap_threshold A number, between 0 and 1. Pairs of ROIs with an intersection over union larger than this threshold are merged.
 #' @param debug If TRUE, return as well the debug information
@@ -115,11 +115,32 @@ findPeaksImpl1D <- function(
     verbose = FALSE,
     length_in_xunits = 0.07,
     peakwidth_range_xunits = c(0.15, 0.4),
-    peakDetectionCWTParams = list(exclude0scaleAmpThresh = TRUE),
+    peakDetectionCWTParams = defaultPeakDetectionCWTParams(x=x, peakwidth_range_xunits = peakwidth_range_xunits),
     extension_factor = 0,
     iou_overlap_threshold = 0.2,
     debug = FALSE
 ) {
+
+  #### Handle older peakDetectionCWTParams. Remove in 2026-01-01.
+  old_peakDetectionCWTParams <- list(exclude0scaleAmpThresh = TRUE)
+  if (identical(peakDetectionCWTParams, old_peakDetectionCWTParams)) {
+    cli::cli_warn(c(
+      "Deprecated peakDetectionCWTParams value",
+      "!" = "Passing the argument {.code peakDetectionCWTParams = list(exclude0scaleAmpThresh = TRUE)} is deprecated since GCIMS v0.2.1.",
+      "i" = "Please remove the {.arg peakDetectionCWTParams} argument. The new default value is better."
+    ))
+    peakDetectionCWTParams = defaultPeakDetectionCWTParams(
+      x=x,
+      peakwidth_range_xunits = peakwidth_range_xunits
+    )
+  }
+  #### Code handling older peakDetectionCWTParams ends here.
+
+
+  if (!is.null(peakwidth_range_xunits) && !identical(attr(peakDetectionCWTParams, "peakwidth_range_xunits"), peakwidth_range_xunits)) {
+    stop("Please provide a new peakDetectionCWTParams created with the given peakwidth_range_xunits")
+  }
+
   x_step <- x[2L] - x[1L]
   x_length_pts <- units_to_points(length_in_xunits, x_step, must_odd = TRUE)
   spec_length <- length(y)
@@ -129,24 +150,10 @@ findPeaksImpl1D <- function(
   sg_filt <- signal::sgolay(p = 2L, n = x_length_pts, m = 2)
   deriv2 <- sgolayfilt(y, sg_filt)
   # 5. Peaks and Zero-crossings
-  scales <- prep_wav_from_peakwidths(
-    x,
-    peakwidth_min = peakwidth_range_xunits[1L],
-    peakwidth_max = peakwidth_range_xunits[2L]
-  )
-
-  if (verbose) {
-    cli_inform(
-      message = c(
-        "i" = "Using the following scales: {glue::glue_collapse(scales$scales, sep = ', ')}"
-      )
-    )
-  }
 
   peaks_and_zeros_extra <- detect_peaks_and_zeros(
     xmat = matrix(y, ncol = 1),
     rowwise = FALSE,
-    scales = scales,
     peakDetectionCWTParams = peakDetectionCWTParams,
     signals_to_save_extra = if (debug) 1L else NULL,
     xmat_zeros = matrix(deriv2, ncol = 1)