Merge pull request #19 from sipss/mix-alignments

Mix alignments

DESCRIPTION

@@ -106,15 +106,17 @@ Suggests:
     httr,
     knitr,
     labeling,
+    magick,
     plotly (>= 4.10.2),
+    png,
     pow,
     rmarkdown,
     scales,
     testthat (>= 3.0.0),
     viridisLite
 biocViews: Software, Preprocessing, Visualization, Classification,
  Cheminformatics, Metabolomics, DataImport
-RoxygenNote: 7.2.3
+RoxygenNote: 7.3.1
 VignetteBuilder: knitr
 Config/testthat/edition: 3
 Roxygen: list(markdown = TRUE)

NAMESPACE

@@ -22,9 +22,9 @@ export(add_peaklist_rect)
 export(align)
 export(alignDt)
 export(alignPlots)
-export(alignRt)
-export(align_ip)
-export(align_pow)
+export(alignRt_ip)
+export(alignRt_pow)
+export(alignRt_ptw)
 export(baseline)
 export(clusterPeaks)
 export(create_annotations_table)
@@ -58,6 +58,7 @@ export(read_mea)
 export(realize)
 export(rtime)
 export(sampleNames)
+export(show_progress_bar)
 export(smooth)
 export(updateObject)
 exportClasses(DelayedOperation)
@@ -71,9 +72,9 @@ exportMethods("intensity<-")
 exportMethods("peaks<-")
 exportMethods(align)
 exportMethods(alignDt)
-exportMethods(alignRt)
-exportMethods(align_ip)
-exportMethods(align_pow)
+exportMethods(alignRt_ip)
+exportMethods(alignRt_pow)
+exportMethods(alignRt_ptw)
 exportMethods(baseline)
 exportMethods(decimate)
 exportMethods(description)
@@ -113,7 +114,6 @@ importMethodsFrom(Biobase,sampleNames)
 importMethodsFrom(BiocGenerics,updateObject)
 importMethodsFrom(ProtGenerics,"intensity<-")
 importMethodsFrom(ProtGenerics,"peaks<-")
-importMethodsFrom(ProtGenerics,alignRt)
 importMethodsFrom(ProtGenerics,filterRt)
 importMethodsFrom(ProtGenerics,intensity)
 importMethodsFrom(ProtGenerics,peaks)

R/aaa-AllGenerics.R

@@ -88,7 +88,6 @@ setGeneric("plotRIC", function(object, ...) standardGeneric("plotRIC"))
 #' @export
 setGeneric("filterDt", function(object, ...) standardGeneric("filterDt"))
 
-
 #' @describeIn GCIMS-generics Decimate an object
 #'
 #' @param object An object to decimate
@@ -97,7 +96,6 @@ setGeneric("filterDt", function(object, ...) standardGeneric("filterDt"))
 #' @export
 setGeneric("decimate", function(object, ...) standardGeneric("decimate"))
 
-
 #' @describeIn GCIMS-generics Align an object
 #'
 #' @param object An object to align
@@ -106,30 +104,41 @@ setGeneric("decimate", function(object, ...) standardGeneric("decimate"))
 #' @export
 setGeneric("align", function(object, ...) standardGeneric("align"))
 
-#' @describeIn GCIMS-generics Align to injection point in retention time
+#' @describeIn GCIMS-generics Align an object in drift time
 #'
-#' @param object An object to align
+#' @param x An object to align
 #' @param ... Further arguments, possibly used by downstream methods.
 #' @return The object, modified
 #' @export
-setGeneric("align_ip", function(object, ...) standardGeneric("align_ip"))
+setGeneric("alignDt", function(x, ...) standardGeneric("alignDt"))
 
-#' @describeIn GCIMS-generics Align in retention time with POW algorithm
+#' @describeIn GCIMS-generics Align an object in retention time
 #'
-#' @param object An object to align
+#' @param x An object to align
 #' @param ... Further arguments, possibly used by downstream methods.
 #' @return The object, modified
 #' @export
-setGeneric("align_pow", function(object, ...) standardGeneric("align_pow"))
+setGeneric("alignRt_ptw", function(x, ...) standardGeneric("alignRt_ptw"))
 
-#' @describeIn GCIMS-generics Align an object in drift time
+#' @describeIn GCIMS-generics Align an object in retention time
+#'
+#' @param x An object to align
+#' @param ... Further arguments, possibly used by downstream methods.
+#' @return The object, modified
+#' @export
+setGeneric("alignRt_pow", function(x, ...) standardGeneric("alignRt_pow"))
+
+#' @describeIn GCIMS-generics Align an object in retention time
 #'
 #' @param x An object to align
-#' @param y Another object to use as reference
 #' @param ... Further arguments, possibly used by downstream methods.
 #' @return The object, modified
 #' @export
-setGeneric("alignDt", function(x, y, ...) standardGeneric("alignDt"))
+setGeneric("alignRt_ip", function(x, ...) standardGeneric("alignRt_ip"))
+
+# #' @importMethodsFrom ProtGenerics alignRt
+# #' @export
+# setGeneric("alignRt", getGeneric("alignRt", package = "ProtGenerics"))
 
 #' @importMethodsFrom Biobase pData
 #' @export
@@ -153,11 +162,6 @@ setGeneric("sampleNames<-", getGeneric("sampleNames<-", package = "Biobase"))
 #' @export
 setGeneric("updateObject", getGeneric("updateObject", package = "BiocGenerics"))
 
-
-#' @importMethodsFrom ProtGenerics alignRt
-#' @export
-setGeneric("alignRt", getGeneric("alignRt", package = "ProtGenerics"))
-
 #' @importMethodsFrom ProtGenerics filterRt
 #' @export
 setGeneric("filterRt", getGeneric("filterRt", package = "ProtGenerics"))

R/align-GCIMSDataset.R

@@ -4,12 +4,15 @@
 #' Parametric Time Warping correction in retention time
 #'
 #' @param object A [GCIMSDataset] object, modified in-place
+#' @param method Method for alignment, should be "ptw" or "pow"
+#' @param align_ip if TRUE a multiplicative correction will be done in retention time before applying the other algorithm
+#' @param ... additional parameters for POW alignment
 #' @return The modified [GCIMSDataset]
 #' @export
 setMethod(
   "align",
   "GCIMSDataset",
-  function(object) {
+  function(object, method = "pow", align_ip = TRUE, ...) {
     tis_matrix <- getTIS(object)
     ric_matrix <- getRIC(object)
     dt <- dtime(object)
@@ -19,8 +22,10 @@ setMethod(
       dt = dt,
       rt = rt,
       tis_matrix = tis_matrix,
-      ric_matrix = ric_matrix
-    )
+      ric_matrix = ric_matrix,
+      method = method,
+      align_ip = align_ip,
+      ...)
 
     delayed_op <- DelayedOperation(
       name = "align",
@@ -70,30 +75,14 @@ setMethod(
   ds
 }
 
-#' Find the  best retention time reference chromatogram.
-#' @noRd
-#'
-#' @description This function provides the index corresponding to the reference Reactant Ion Chromatogram (RIC) to correct
-#' misalignments in retention time.
-#' @param rics A matrix. Each row correspond to a different RIC. There are as many RICs as samples.
-#' @return  An Integer number that indicates the reference sample.
-#' @examples
-#' rics <- rbind(
-#'   dnorm(1:100, mean=50, sd =1),
-#'   dnorm(1:100, mean=51, sd =1),
-#'   dnorm(1:100, mean=52, sd =1)
-#' )
-#' find_reference_ric(rics) == 2L
-find_reference_ric <- function(rics){
-  ref_ric_sample_idx <- ptw::bestref(rics)$best.ref
-  return(ref_ric_sample_idx)
-}
-
 
-
-alignParams <- function(dt, rt, tis_matrix, ric_matrix) {
+alignParams <- function(dt, rt, tis_matrix, ric_matrix, method, align_ip, ...) {
   # Optimize ret time alignment parameters:
-  ref_ric_sample_idx <- find_reference_ric(ric_matrix)
+  if (method == "ptw"){
+    ref_ric_sample_idx <- ptw::bestref(ric_matrix)$best.ref
+  } else {
+    ref_ric_sample_idx <- pow::select_reference(ric_matrix)
+  }
   # Select reference RIC
   ric_ref <- as.numeric(ric_matrix[ref_ric_sample_idx, ])
 
@@ -102,7 +91,19 @@ alignParams <- function(dt, rt, tis_matrix, ric_matrix) {
   rip_ref_idx <- round(stats::median(rip_position, na.rm = TRUE))
   rip_ref_ms <- dt[rip_ref_idx]
 
-  list(rip_ref_ms = rip_ref_ms, ric_ref = ric_ref, ric_ref_rt = rt)
+  # 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]
+
+  list(rip_ref_ms = rip_ref_ms,
+       ric_ref = ric_ref,
+       ric_ref_rt = rt,
+       ip_ref_s = ip_ref_s,
+       method = method,
+       align_ip = align_ip,
+       ...)
 }
 
 #' Plots to interpret alignment results

R/align-GCIMSSample.R

@@ -3,19 +3,37 @@
 #' @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 method Method for alignment, should be "ptw" or "pow"
+#' @param align_ip if TRUE a multiplicative correction will be done in retention time before applying the other algorithm
 #' @return The modified [GCIMSSample]
 #' @export
 methods::setMethod(
   "align", "GCIMSSample",
-  function(object, rip_ref_ms, ric_ref, ric_ref_rt){
+  function(object, rip_ref_ms, ric_ref, ric_ref_rt, ip_ref_s, method, align_ip, ...){
     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
     object <- alignDt(object, rip_ref_ms = rip_ref_ms)
     if (all(is.na(object@data))) {
       cli_abort("After aligning drift times, all the data matrix of {description(object)} are missing values. This should not happen")
     }
-    object <- alignRt(object, ric_ref = ric_ref, ric_ref_rt = ric_ref_rt)
+
+    # Align in retention time
+    if (align_ip){
+      object <- alignRt_ip(object, ip_ref_s = ip_ref_s)
+    }
+    if (method == "ptw") {
+      object <- alignRt_ptw(object, ric_ref = ric_ref, ric_ref_rt = ric_ref_rt)
+    } else {
+      if(method != "pow"){
+        cli_warn("{method} is not a valid method, using POW instead")
+      }
+      object <- alignRt_pow(object, ric_ref = ric_ref, ric_ref_rt = ric_ref_rt, ...)
+    }
+
     if (all(is.na(object@data))) {
       cli_abort("After aligning drift and retention times, all the data matrix of {description(object)} are missing values. This should not happen")
     }
@@ -37,23 +55,17 @@ methods::setMethod(
 
 #' Align a GCIMSSample in drift time with a multiplicative correction
 #' @param x A [GCIMSSample] object
-#' @param y Not used
 #' @param rip_ref_ms The position of the RIP in ms
 #' @return The modified [GCIMSSample]
 #' @export
 methods::setMethod(
   "alignDt", "GCIMSSample",
-  function(x, y, rip_ref_ms) {
-    if (missing(rip_ref_ms) && !missing(y)) {
-      cli_warn("Please provide rip_ref_ms as a named argument")
-      rip_ref_ms <- y
-    }
+  function(x, rip_ref_ms) {
     tis <- getTIS(x)
     dt <- dtime(x)
     rip_pos_ms <- dt[which.max(tis)]
     Kcorr <- rip_ref_ms/rip_pos_ms
 
-    dt <- dtime(x)
     int_mat <- intensity(x)
     dt_corr <- Kcorr * dt
     for (j in seq_len(ncol(int_mat))) {
@@ -91,18 +103,76 @@ 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
+#' @return The modified [GCIMSSample]
+#' @export
+methods::setMethod(
+  "alignRt_ip","GCIMSSample",
+  function(x, ip_ref_s) {
+    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
+    x
+  })
+
+#' Align a GCIMSSample in retention time with parametric optimized warping
+#' @param x A [GCIMSSample] object
+#' @param ric_ref The reference Reverse Ion Chromatogram
+#' @param ric_ref_rt The retention times corresponding to `ric_ref`
+#' @return The modified [GCIMSSample]
+#' @export
+methods::setMethod(
+  "alignRt_pow", "GCIMSSample",
+  function(x,
+           ric_ref,
+           ric_ref_rt,
+           lambdas = pracma::logspace(-2, 4, 31),
+           p = 10,
+           max_it = 1000) {
+    ric <- getRIC(x)
+    v <- rep(1, length(ric_ref))
+    iv <- seq(2, length(ric_ref) - 1, by = p)
+    v[iv] <- 0L
+    W <- Matrix::Diagonal(x = v)
+    result_val <- pow:::val(ric, ric_ref, W, iv, lambdas, fom ="rms")
+    e_ix <- result_val$e
+    ti_ix <- result_val$ti
+    e_ix[ti_ix == 1] <- NA
+    lambdas[ti_ix == 1] <- NA
+    best_lambda <- lambdas[which.min(e_ix)]
+    w <- pow::pow(ric, best_lambda, ric_ref, max_it = max_it)
+
+    int <- intensity(x)
+    inter <- t(apply(int, 1, pow::interpolation, w = w, return = FALSE))
+    sel <- !is.na(inter[1,])
+    x@retention_time <- ric_ref_rt[sel]
+    x@data <- inter[,sel]
+
+
+    x@proc_params$align$rt_poly_coefs <- c(0, 1)
+
+    return(x)
+  })
+
 #' Align a GCIMSSample in retention time using parametric time warping
 #' @param x A [GCIMSSample] object
-#' @param y Not used
 #' @param ric_ref The reference Reverse Ion Chromatogram
 #' @param ric_ref_rt The retention times corresponding to `ric_ref`
 #' @return The modified [GCIMSSample]
-#' @importMethodsFrom ProtGenerics alignRt
 #' @export
 methods::setMethod(
-  "alignRt",
-  signature = c(x = "GCIMSSample", y = "ANY"),
-  function(x, y, ric_ref, ric_ref_rt) {
+  "alignRt_ptw", "GCIMSSample",
+  function(x, ric_ref, ric_ref_rt) {
     optimize_polynomial_order <- function(ric_sample, ric_ref) {
       correction_type_options <- seq.int(0, 5)
       poly_orders <- seq.int(1, 5)
@@ -142,10 +212,10 @@ methods::setMethod(
       correction_type_options[which.max(corr)]
     }
 
-    if (missing(ric_ref) && !missing(y)) {
-      cli_warn("Please provide ric_ref as a named argument")
-      ric_ref <- y
-    }
+    # if (missing(ric_ref) && !missing(y)) {
+    #   cli_warn("Please provide ric_ref as a named argument")
+    #   ric_ref <- y
+    # }
 
     ric <- getRIC(x)
     rt <- rtime(x)