From 41fe33034bc863f3a3a542861e2e63a913145b30 Mon Sep 17 00:00:00 2001 From: lalo-caballero Date: Wed, 20 Mar 2024 09:31:24 +0000 Subject: [PATCH] optional align in dt --- R/align-GCIMSDataset.R | 16 ++++++++++------ R/align-GCIMSSample.R | 37 +++++++++++++++++++++++++++---------- 2 files changed, 37 insertions(+), 16 deletions(-) diff --git a/R/align-GCIMSDataset.R b/R/align-GCIMSDataset.R index 9050bf0..20c2e54 100644 --- a/R/align-GCIMSDataset.R +++ b/R/align-GCIMSDataset.R @@ -4,7 +4,9 @@ #' 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 method_rt Method for alignment, should be "ptw" or "pow" +#' if pow is selected the package "pow must be installed, to do so visit: +#' https://github.com/sipss/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] @@ -12,7 +14,7 @@ setMethod( "align", "GCIMSDataset", - function(object, method = "pow", align_ip = TRUE, ...) { + function(object, method_rt = "ptw", align_dt = TRUE, align_ip = TRUE, ...) { tis_matrix <- getTIS(object) ric_matrix <- getRIC(object) dt <- dtime(object) @@ -23,7 +25,8 @@ setMethod( rt = rt, tis_matrix = tis_matrix, ric_matrix = ric_matrix, - method = method, + method_rt = method_rt, + align_dt = align_dt, align_ip = align_ip, ...) @@ -76,9 +79,9 @@ setMethod( } -alignParams <- function(dt, rt, tis_matrix, ric_matrix, method, align_ip, ...) { +alignParams <- function(dt, rt, tis_matrix, ric_matrix, method_rt, align_dt, align_ip, ...) { # Optimize ret time alignment parameters: - if (method == "ptw"){ + if (method_rt == "ptw"){ ref_ric_sample_idx <- ptw::bestref(ric_matrix)$best.ref } else { ref_ric_sample_idx <- pow::select_reference(ric_matrix) @@ -101,7 +104,8 @@ alignParams <- function(dt, rt, tis_matrix, ric_matrix, method, align_ip, ...) { ric_ref = ric_ref, ric_ref_rt = rt, ip_ref_s = ip_ref_s, - method = method, + method_rt = method_rt, + align_dt = align_dt, align_ip = align_ip, ...) } diff --git a/R/align-GCIMSSample.R b/R/align-GCIMSSample.R index e16692a..08794b3 100644 --- a/R/align-GCIMSSample.R +++ b/R/align-GCIMSSample.R @@ -4,19 +4,24 @@ #' @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 method_rt 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, ip_ref_s, method, align_ip, ...){ + function(object, rip_ref_ms, ric_ref, ric_ref_rt, ip_ref_s, method_rt, align_dt, 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 (align_dt){ + object <- alignDt(object, rip_ref_ms = rip_ref_ms) + } else{ + object@proc_params$align$dt_kcorr <- 1 + } + 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") } @@ -25,15 +30,20 @@ methods::setMethod( if (align_ip){ object <- alignRt_ip(object, ip_ref_s = ip_ref_s) } - if (method == "ptw") { + if (method_rt == "ptw") { object <- alignRt_ptw(object, ric_ref = ric_ref, ric_ref_rt = ric_ref_rt) + } else if(method_rt == "pow"){ + object <- alignRt_pow(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") + if (!"align" %in% names(object@proc_params)) { + object@proc_params$align <- list() } - object <- alignRt_pow(object, ric_ref = ric_ref, ric_ref_rt = ric_ref_rt, ...) + object@proc_params$align$rt_min_s <- ric_ref_rt[1] + object@proc_params$align$rt_max_s <- ric_ref_rt[length(ric_ref_rt)] + object@proc_params$align$rt_poly_coefs <- c(0, 1) } + 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") } @@ -138,19 +148,20 @@ methods::setMethod( ric_ref_rt, lambdas = pracma::logspace(-2, 4, 31), p = 10, - max_it = 1000) { + max_it = 5000, + lambda1 = 10^6) { 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") + result_val <- pow:::val(ric, ric_ref, W, iv, lambdas, fom ="rms", lambda1 = lambda1) 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) + w <- pow::pow(ric, best_lambda, ric_ref, max_it = max_it, lambda1= lambda1) int <- intensity(x) inter <- t(apply(int, 1, pow::interpolation, w = w, return = FALSE)) @@ -159,6 +170,12 @@ methods::setMethod( x@data <- inter[,sel] + if (!"align" %in% names(x@proc_params)) { + x@proc_params$align <- list() + } + x@proc_params$align$warp <- w + #x@proc_params$align$rt_min_s <- x@retention_time[1] + #x@proc_params$align$rt_max_s <- x@retention_time[length(x@retention_time)] x@proc_params$align$rt_poly_coefs <- c(0, 1) return(x)