Add findpeaks for GCIMSSpectrum and update alignment options

R/aaa-class-GCIMSSpectrum.R

@@ -11,6 +11,8 @@
 #' @slot retention_time_idx The index or indices used to get the intensity
 #' @slot retention_time_s The retention times corresponding to the retention time indices.
 #' @slot description A string with a description (used as plot title, useful e.g. to know the sample it came from)
+#' @slot peaks  A data frame with peaks, use [findPeaks()] to fill it, or [peaks()] to set/get it
+#' @slot peaks_debug_info A list with arbitrary debug information from [findPeaks()]
 #'
 #' @export
 methods::setClass(
@@ -21,7 +23,9 @@ methods::setClass(
     baseline = "numericOrNULL",
     retention_time_idx = "numeric",
     retention_time_s = "numeric",
-    description = "character"
+    description = "character",
+    peaks = "DataFrameOrNULL",
+    peaks_debug_info = "listOrNULL"
   )
 )
 

R/align-GCIMSDataset.R

@@ -8,13 +8,15 @@
 #' https://github.com/sipss/pow
 #' @param align_dt if `TRUE` the drift time axis will be aligned using a multiplicative correction
 #' @param align_ip if TRUE a multiplicative correction will be done in retention time before applying the other algorithm
+#' @param reference_sample_idx One number, the index of the sample to use as reference for the alignment in retention time, if NULL the reference will be calculated automatically depending on the method
+#' @param ploynomial_order_ptw One number, by default 2, the maximum order of the polynomial for the parametric time warping alignment.
 #' @param ... additional parameters for POW alignment
 #' @return The modified [GCIMSDataset]
 #' @export
 setMethod(
   "align",
   "GCIMSDataset",
-  function(object, method_rt = "ptw", align_dt = TRUE, align_ip = TRUE, ...) {
+  function(object, method_rt = "ptw", align_dt = TRUE, align_ip = TRUE, reference_sample_idx = NULL, ploynomial_order_ptw = 5, ...) {
     tis_matrix <- getTIS(object)
     ric_matrix <- getRIC(object)
     dt <- dtime(object)
@@ -28,6 +30,8 @@ setMethod(
       method_rt = method_rt,
       align_dt = align_dt,
       align_ip = align_ip,
+      ref_ric_sample_idx = reference_sample_idx,
+      ploynomial_order_ptw = ploynomial_order_ptw,
       ...)
 
     delayed_op <- DelayedOperation(
@@ -79,12 +83,14 @@ setMethod(
 }
 
 
-alignParams <- function(dt, rt, tis_matrix, ric_matrix, method_rt, align_dt, align_ip, ...) {
+alignParams <- function(dt, rt, tis_matrix, ric_matrix, method_rt, align_dt, align_ip, ref_ric_sample_idx, ploynomial_order_ptw, ...) {
   # Optimize ret time alignment parameters:
-  if (method_rt == "ptw"){
-    ref_ric_sample_idx <- ptw::bestref(ric_matrix)$best.ref
-  } else {
-    ref_ric_sample_idx <- pow::select_reference(ric_matrix)
+  if (is.null(ref_ric_sample_idx)){
+    if (method_rt == "pow"){
+      ref_ric_sample_idx <- pow::select_reference(ric_matrix)
+    } else {
+      ref_ric_sample_idx <- ptw::bestref(ric_matrix)$best.ref
+    }
   }
   # Select reference RIC
   ric_ref <- as.numeric(ric_matrix[ref_ric_sample_idx, ])
@@ -96,17 +102,19 @@ alignParams <- function(dt, rt, tis_matrix, ric_matrix, method_rt, align_dt, ali
 
   # IP alignment parameters
 
-  ip_position <- apply(ric_matrix, 1L ,which.min)
-  ip_ref_idx <- round(stats::median(ip_position, na.rm = TRUE))
-  ip_ref_s <- rt[ip_ref_idx]
+  mins <- apply(ric_matrix, 1L ,which.min)
+  rt_ref <- rt[1 : (length(rt) - (max(mins) - min(mins)))]
+  min_start <- min(mins) - 1
 
   list(rip_ref_ms = rip_ref_ms,
        ric_ref = ric_ref,
        ric_ref_rt = rt,
-       ip_ref_s = ip_ref_s,
+       min_start = min_start,
+       rt_ref = rt_ref,
        method_rt = method_rt,
        align_dt = align_dt,
        align_ip = align_ip,
+       ploynomial_order_ptw = ploynomial_order_ptw,
        ...)
 }
 
@@ -177,10 +185,8 @@ alignPlots <- function(object) {
       x = names(object$align$dt_kcorr),
       y = object$align$dt_kcorr
     )
-  ) + ggplot2::geom_col(ggplot2::aes(x = .data$x, y = .data$y)) +
-    ggplot2::labs(x = "SampleID", y = "Multiplicative factor (drift time correction, unitless)") +
-    ggplot2::coord_flip()
-
+  ) + ggplot2::geom_point(ggplot2::aes(x = .data$x, y = .data$y)) +
+    ggplot2::labs(x = "SampleID", y = "Multiplicative factor (drift time correction, unitless)")
 
   list(
     rt_diff_plot = rt_diff_plot,

R/align-GCIMSSample.R

@@ -3,22 +3,24 @@
 #' @param rip_ref_ms The reference position of the Reactant Ion Peak in the dataset (in ms)
 #' @param ric_ref The reference Reverse Ion Chromatogram
 #' @param ric_ref_rt The retention times corresponding to `ric_ref`
-#' @param ip_ref_s Injection point reference for retention time
+#' @param min_start minimun injection point, to calculate where to begin the spectrums and cut as few points as posible, to be used in injection point alignment
+#' @param rt_ref retention time reference for alignment to injection point
 #' @param method_rt Method for alignment, should be "ptw" or "pow"
 #' @param align_dt if `TRUE`, align the drift time axis using a multiplicative correction
 #' @param align_ip if TRUE a multiplicative correction will be done in retention time before applying the other algorithm
+#' @param ploynomial_order_ptw One number, by default 2, the maximum order of the polynomial for the parametric time warping alignment.
 #' @param ... Additional arguments passed on to the alignment method.
 #' @return The modified [GCIMSSample]
 #' @export
 methods::setMethod(
   "align", "GCIMSSample",
-  function(object, rip_ref_ms, ric_ref, ric_ref_rt, ip_ref_s, method_rt, align_dt, align_ip, ...){
+  function(object, rip_ref_ms, ric_ref, ric_ref_rt, min_start, rt_ref, method_rt, align_dt, align_ip, ploynomial_order_ptw, ...){
     if (all(is.na(object@data))) {
       cli_abort("All the data matrix of {description(object)} are missing values. Align is impossible")
     }
 
     # Align in drift time
-    if (align_dt){
+    if (align_dt) {
       object <- alignDt(object, rip_ref_ms = rip_ref_ms)
     } else{
       object@proc_params$align$dt_kcorr <- 1
@@ -30,10 +32,10 @@ methods::setMethod(
 
     # Align in retention time
     if (align_ip){
-      object <- alignRt_ip(object, ip_ref_s = ip_ref_s)
+      object <- alignRt_ip(object, min_start = min_start, rt_ref = rt_ref)
     }
     if (method_rt == "ptw") {
-      object <- alignRt_ptw(object, ric_ref = ric_ref, ric_ref_rt = ric_ref_rt)
+      object <- alignRt_ptw(object, ric_ref = ric_ref, ric_ref_rt = ric_ref_rt, ploynomial_order_ptw)
     } else if(method_rt == "pow"){
       object <- alignRt_pow(object, ric_ref = ric_ref, ric_ref_rt = ric_ref_rt, ...)
     } else {
@@ -117,23 +119,17 @@ methods::setMethod(
 
 #' Align a GCIMSSample in retention time with a multiplicative correction
 #' @param x A [GCIMSSample] object
-#' @param ip_ref_s The position of the injection point in s
+#' @param min_start minimun injection point, to calculate where to begin the spectrums and cut as few points as posible
+#' @param rt_ref retention time reference
 #' @return The modified [GCIMSSample]
 #' @export
 methods::setMethod(
   "alignRt_ip","GCIMSSample",
-  function(x, ip_ref_s) {
+  function(x, min_start, rt_ref) {
     ric <- getRIC(x)
-    rt <- rtime(x)
-    ip_pos_s <- rt[which.min(ric)]
-    Kcorr <- ip_ref_s/ip_pos_s
-
-    int_mat <- intensity(x)
-    rt_corr <- Kcorr * rt
-    for (j in seq_len(nrow(int_mat))) {
-      int_mat[j,] <- signal::interp1(rt_corr, int_mat[j, ], rt, extrap = TRUE)
-    }
-    intensity(x) <- int_mat
+    injection_point <- which.min(ric)
+    x@retention_time <- rt_ref
+    x@data <- x@data[, (injection_point - min_start):((injection_point - min_start)+length(rt_ref)-1)]
     x
   })
 
@@ -191,14 +187,15 @@ methods::setMethod(
 #' @param x A [GCIMSSample] object
 #' @param ric_ref The reference Reverse Ion Chromatogram
 #' @param ric_ref_rt The retention times corresponding to `ric_ref`
+#' @param ploynomial_order_ptw maximum order of the polynomial to be used
 #' @return The modified [GCIMSSample]
 #' @export
 methods::setMethod(
   "alignRt_ptw", "GCIMSSample",
-  function(x, ric_ref, ric_ref_rt) {
+  function(x, ric_ref, ric_ref_rt, ploynomial_order_ptw) {
     optimize_polynomial_order <- function(ric_sample, ric_ref) {
-      correction_type_options <- seq.int(0, 5)
-      poly_orders <- seq.int(1, 5)
+      correction_type_options <- seq.int(0, ploynomial_order_ptw)
+      poly_orders <- seq.int(1, ploynomial_order_ptw)
       xi <- seq_len(length(ric_sample))
       corr <- numeric(length(poly_orders) + 1L)
       corr[1] <- stats::cor(ric_ref, ric_sample, use = "complete.obs")

R/findPeaks-GCIMSChromatogram.R

@@ -20,26 +20,3 @@ setMethod(
     object
   }
 )
-
-#' Peak detection for a GCIMSSpectrum
-#' @param object A [GCIMSSpectrum] object
-#' @inheritDotParams findPeaksImpl1D -x -y
-#' @return The modified [GCIMSSpectrum], with a peak list
-#' @family GCIMSSpectrum
-#' @export
-setMethod(
-  "findPeaks",
-  "GCIMSSpectrum",
-  function(object, ...) {
-    dt <- dtime(object)
-    intens <- intensity(object)
-    peak_list_and_debug_info <- findPeaksImpl1D(
-      x = dt,
-      y = intens,
-      ...
-    )
-    object@peaks_debug_info <- peak_list_and_debug_info$debug_info
-    peaks(object) <- peak_list_and_debug_info$peak_list
-    object
-  }
-)

R/findPeaks-GCIMSSpectrum.R

@@ -0,0 +1,22 @@
+#' Peak detection for a GCIMSSpectrum
+#' @param object A [GCIMSSpectrum] object
+#' @inheritDotParams findPeaksImpl1D -x -y
+#' @return The modified [GCIMSSpectrum], with a peak list
+#' @family GCIMSSpectrum
+#' @export
+setMethod(
+  "findPeaks",
+  "GCIMSSpectrum",
+  function(object, ...) {
+    dt <- dtime(object)
+    intens <- intensity(object)
+    peak_list_and_debug_info <- findPeaksImpl1D(
+      x = dt,
+      y = intens,
+      ...
+    )
+    object@peaks_debug_info <- peak_list_and_debug_info$debug_info
+    peaks(object) <- peak_list_and_debug_info$peak_list
+    object
+  }
+)

R/peaks-GCIMSSpectrum.R

@@ -0,0 +1,36 @@
+#' @describeIn GCIMSSpectrum-class Get the peak list
+#' @importMethodsFrom ProtGenerics peaks
+#' @return A data frame with the peaks in the spectrum
+#' @export
+setMethod(
+  "peaks",
+  "GCIMSSpectrum",
+  function(object) {
+    if (is.null(object@peaks)) {
+      stop("Please run findPeaks() on the spectrum first")
+    }
+    if(nrow(object@peaks)==0){
+      p <- NULL
+    }else{
+      p <- tibble::as_tibble(cbind(SampleID = object@description, object@peaks))
+    }
+    return(p)
+  }
+)
+
+#' @describeIn GCIMSSpectrum-class Set the peak list
+#' @param value A data frame with the peak list
+#' @return The GCIMSSpectrum object
+#' @importMethodsFrom ProtGenerics "peaks<-"
+#' @export
+setMethod(
+  "peaks<-",
+  "GCIMSSpectrum",
+  function(object, value) {
+    value <- methods::as(value, "DataFrame")
+    # We don't want this column:
+    value$SampleID <- NULL
+    object@peaks <- value
+    object
+  }
+)