diff --git a/DESCRIPTION b/DESCRIPTION index eee56db..f199317 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -107,9 +107,11 @@ Suggests: knitr, labeling, magick, + Matrix, plotly, png, pow, + pracma, rmarkdown, scales, testthat (>= 3.0.0), diff --git a/NAMESPACE b/NAMESPACE index fe727e2..d67a944 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -54,6 +54,7 @@ export(plotRIC) export(plotRaw) export(plotTIS) export(plot_interactive) +export(prealign) export(read_mea) export(realize) export(rtime) @@ -71,10 +72,6 @@ exportMethods("description<-") exportMethods("intensity<-") exportMethods("peaks<-") exportMethods(align) -exportMethods(alignDt) -exportMethods(alignRt_ip) -exportMethods(alignRt_pow) -exportMethods(alignRt_ptw) exportMethods(baseline) exportMethods(decimate) exportMethods(description) @@ -91,6 +88,7 @@ exportMethods(peaks) exportMethods(plot) exportMethods(plotRIC) exportMethods(plotTIS) +exportMethods(prealign) exportMethods(rtime) exportMethods(smooth) exportMethods(updateObject) diff --git a/R/aaa-AllGenerics.R b/R/aaa-AllGenerics.R index 128ef9a..ef47568 100644 --- a/R/aaa-AllGenerics.R +++ b/R/aaa-AllGenerics.R @@ -110,35 +110,7 @@ setGeneric("align", function(object, ...) standardGeneric("align")) #' @param ... Further arguments, possibly used by downstream methods. #' @return The object, modified #' @export -setGeneric("alignDt", function(object, ...) standardGeneric("alignDt")) - -#' @describeIn GCIMS-generics Align an object in retention time -#' -#' @param object An object to align -#' @param ... Further arguments, possibly used by downstream methods. -#' @return The object, modified -#' @export -setGeneric("alignRt_ptw", function(object, ...) standardGeneric("alignRt_ptw")) - -#' @describeIn GCIMS-generics Align an object in retention time -#' -#' @param object An object to align -#' @param ... Further arguments, possibly used by downstream methods. -#' @return The object, modified -#' @export -setGeneric("alignRt_pow", function(object, ...) standardGeneric("alignRt_pow")) - -#' @describeIn GCIMS-generics Align an object in retention time -#' -#' @param object An object to align -#' @param ... Further arguments, possibly used by downstream methods. -#' @return The object, modified -#' @export -setGeneric("alignRt_ip", function(object, ...) standardGeneric("alignRt_ip")) - -# #' @importMethodsFrom ProtGenerics alignRt -# #' @export -# setGeneric("alignRt", getGeneric("alignRt", package = "ProtGenerics")) +setGeneric("prealign", function(object, ...) standardGeneric("prealign")) #' @importMethodsFrom Biobase pData #' @export diff --git a/R/align-GCIMSDataset.R b/R/align-GCIMSDataset.R index 2a6678e..0832ad7 100644 --- a/R/align-GCIMSDataset.R +++ b/R/align-GCIMSDataset.R @@ -4,47 +4,83 @@ #' Parametric Time Warping correction in retention time #' @param object A [GCIMSDataset] object, modified in-place #' @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: +#' if pow is selected the package "pow" must be installed, to do so visit: #' 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 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, reference_sample_idx = NULL, ploynomial_order_ptw = 5, ...) { - tis_matrix <- getTIS(object) - ric_matrix <- getRIC(object) - dt <- dtime(object) - rt <- rtime(object) - - align_params <- alignParams( - dt = dt, - rt = rt, - tis_matrix = tis_matrix, - ric_matrix = ric_matrix, - 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( - name = "align", - fun = align, - params = align_params, - fun_extract = .align_fun_extract, - fun_aggregate = .align_fun_aggregate - ) - object$appendDelayedOp(delayed_op) - # We recompute these, but maybe we could just reset them to zero... - object$extract_dtime_rtime() - object$extract_RIC_and_TIS() + function(object, + method_rt = "ptw", + align_dt = TRUE, + align_ip = TRUE, + reference_sample_idx = NULL, + ...) { + if (align_dt | align_ip){ + tis_matrix <- getTIS(object) + ric_matrix <- getRIC(object) + dt <- dtime(object) + rt <- rtime(object) + + # Optimize drift time alignment parameters: + rip_position <- apply(tis_matrix, 1L, which.max) + rip_ref_idx <- round(stats::median(rip_position, na.rm = TRUE)) + rip_ref_ms <- dt[rip_ref_idx] + + # IP alignment parameters + mins <- apply(ric_matrix, 1L ,which.min) + rt_ref <- rt[1 : (length(rt) - (max(mins) - min(mins)))] + min_start <- min(mins) - 1 + + delayed_op <- DelayedOperation( + name = "prealign", + fun = prealign, + params = list(align_dt = align_dt, + align_ip = align_ip, + rip_ref_ms = rip_ref_ms, + min_start = min_start, + rt_ref = rt_ref), + fun_extract = .align_fun_extract, + fun_aggregate = .align_fun_aggregate + ) + object$appendDelayedOp(delayed_op) + object$extract_dtime_rtime() + object$extract_RIC_and_TIS() + } + + if (method_rt == "ptw" | method_rt == "pow"){ + rt <- rtime(object) + ric_matrix <- getRIC(object) + # Optimize ret time alignment parameters: + if (is.null(reference_sample_idx)){ + if (method_rt == "pow"){ + reference_sample_idx <- pow::select_reference(ric_matrix) + } else { + reference_sample_idx <- ptw::bestref(ric_matrix)$best.ref + } + } + # Select reference RIC + ric_ref <- as.numeric(ric_matrix[reference_sample_idx, ]) + + delayed_op <- DelayedOperation( + name = "align", + fun = align, + params = list(method_rt = method_rt, + ric_ref = ric_ref, + ric_ref_rt = rt, + ...), + fun_extract = .align_fun_extract, + fun_aggregate = .align_fun_aggregate + ) + object$appendDelayedOp(delayed_op) + object$extract_dtime_rtime() + object$extract_RIC_and_TIS() + } invisible(object) } ) @@ -52,7 +88,7 @@ setMethod( .align_fun_extract <- function(x) { list( dt_kcorr = x@proc_params$align$dt_kcorr, - rt_poly_coefs = x@proc_params$align$rt_poly_coefs + w = x@proc_params$align$w ) } @@ -61,63 +97,27 @@ setMethod( ds$align <- list() } ds$align[["dt_kcorr"]] <- purrr::map_dbl(extracted_obj, "dt_kcorr") - poly_order <- purrr::map_int( - extracted_obj, - function(obj) { - length(obj$rt_poly_coefs) - 1L - } - ) - ds$align[["rt_poly_order"]] <- poly_order - poly_coefs <- matrix(0.0, nrow = length(extracted_obj), ncol = max(poly_order) + 1L) - poly_coefs[,2] <- 1 - for (i in seq_along(extracted_obj)) { - coefs <- extracted_obj[[i]]$rt_poly_coefs - poly_coefs[i, seq_along(coefs)] <- coefs - } - dimnames(poly_coefs) <- list( - SampleID = sampleNames(ds), - PolyOrder = paste0("Order_", seq_len(ncol(poly_coefs)) - 1L) - ) - ds$align[["rt_poly_coefs"]] <- poly_coefs + # poly_order <- purrr::map_int( + # extracted_obj, + # function(obj) { + # length(obj$rt_poly_coefs) - 1L + # } + # ) + # ds$align[["rt_poly_order"]] <- poly_order + # poly_coefs <- matrix(0.0, nrow = length(extracted_obj), ncol = max(poly_order) + 1L) + # poly_coefs[,2] <- 1 + # for (i in seq_along(extracted_obj)) { + # coefs <- extracted_obj[[i]]$rt_poly_coefs + # poly_coefs[i, seq_along(coefs)] <- coefs + # } + # dimnames(poly_coefs) <- list( + # SampleID = sampleNames(ds), + # PolyOrder = paste0("Order_", seq_len(ncol(poly_coefs)) - 1L) + # ) + # ds$align[["rt_poly_coefs"]] <- poly_coefs ds } - -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 (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, ]) - - # Optimize drift time alignment parameters: - rip_position <- apply(tis_matrix, 1L, which.max) - rip_ref_idx <- round(stats::median(rip_position, na.rm = TRUE)) - rip_ref_ms <- dt[rip_ref_idx] - - # IP alignment parameters - - 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, - 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, - ...) -} - #' Plots to interpret alignment results #' #' @@ -152,26 +152,26 @@ alignPlots <- function(object) { tibble::as_tibble(reshape2::melt(dt_diff, value.name = "correction_ms")) } - get_corrected_rt <- function(rt_before, rt_poly_coefs) { - rt_after <- eval_poly(rt_before, rt_poly_coefs) - rt_before_mat <- matrix(rt_before, nrow = nrow(rt_poly_coefs), ncol = length(rt_before), byrow = TRUE) - rt_diff <- rt_after - rt_before_mat - dimnames(rt_diff) <- list( - SampleID = rownames(rt_poly_coefs), - ret_time_s = rt_before - ) - tibble::as_tibble(reshape2::melt(rt_diff, value.name = "correction_s")) - } + # get_corrected_rt <- function(rt_before, rt_poly_coefs) { + # rt_after <- eval_poly(rt_before, rt_poly_coefs) + # rt_before_mat <- matrix(rt_before, nrow = nrow(rt_poly_coefs), ncol = length(rt_before), byrow = TRUE) + # rt_diff <- rt_after - rt_before_mat + # dimnames(rt_diff) <- list( + # SampleID = rownames(rt_poly_coefs), + # ret_time_s = rt_before + # ) + # tibble::as_tibble(reshape2::melt(rt_diff, value.name = "correction_s")) + # } object$realize() rt <- rtime(object) dt <- dtime(object) - rt_diff <- get_corrected_rt(rt, object$align$rt_poly_coefs) - - rt_diff_plot <- ggplot2::ggplot(rt_diff) + - ggplot2::geom_line(ggplot2::aes(x = .data$ret_time_s, y = .data$correction_s, group = .data$SampleID, color = .data$SampleID)) + - ggplot2::labs(x = "Retention time (s)", y = "Ret. time correction (s)", color = "SampleID") + # rt_diff <- get_corrected_rt(rt, object$align$rt_poly_coefs) + # + # rt_diff_plot <- ggplot2::ggplot(rt_diff) + + # ggplot2::geom_line(ggplot2::aes(x = .data$ret_time_s, y = .data$correction_s, group = .data$SampleID, color = .data$SampleID)) + + # ggplot2::labs(x = "Retention time (s)", y = "Ret. time correction (s)", color = "SampleID") dt_diff <- get_corrected_dt(dt, object$align$dt_kcorr) @@ -189,7 +189,7 @@ alignPlots <- function(object) { ggplot2::labs(x = "SampleID", y = "Multiplicative factor (drift time correction, unitless)") list( - rt_diff_plot = rt_diff_plot, + #rt_diff_plot = rt_diff_plot, dt_diff_plot = dt_diff_plot, dt_kcorr_plot = dt_kcorr_plot ) diff --git a/R/align-GCIMSSample.R b/R/align-GCIMSSample.R index 9b0c75a..477bfc1 100644 --- a/R/align-GCIMSSample.R +++ b/R/align-GCIMSSample.R @@ -1,55 +1,84 @@ -#' Align a GCIMSSample object, in drift and retention time +#' Align a GCIMSSample object, in retention time #' @param object A [GCIMSSample] object -#' @param rip_ref_ms The reference position of the Reactant Ion Peak in the dataset (in ms) +#' @param method_rt Method for alignment, should be "ptw" or "pow" #' @param ric_ref The reference Reverse Ion Chromatogram #' @param ric_ref_rt The retention times corresponding to `ric_ref` -#' @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, min_start, rt_ref, method_rt, align_dt, align_ip, ploynomial_order_ptw, ...){ + function(object, method_rt, ric_ref, ric_ref_rt, ...){ + 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 (!"align" %in% names(object@proc_params)) { + object@proc_params$align <- list() + } + 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$w <- ric_ref_rt + } + if (all(is.na(object@data))) { - cli_abort("All the data matrix of {description(object)} are missing values. Align is impossible") + cli_abort("After aligning retention time, all the data matrix of {description(object)} are missing values. This should not happen") + } + + rt_range <- c( + object@proc_params$align$rt_min_s, + object@proc_params$align$rt_max_s + ) + object <- subset(object, rt_range = rt_range) + + if (all(is.na(object@data))) { + cli_abort("After aligning and subsetting all the data matrix of {description(object)} are missing values. This should not happen") } + object + }) +#' Align a GCIMSSample object, in drift time and to the injection point in retention time +#' @param object A [GCIMSSample] object +#' @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 rip_ref_ms The reference position of the Reactant Ion Peak in the dataset (in ms) +#' @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 +#' @return The modified [GCIMSSample] +#' @export +methods::setMethod( + "prealign", "GCIMSSample", + function(object, align_dt, align_ip, rip_ref_ms, min_start, rt_ref){ + 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) { object <- alignDt(object, rip_ref_ms = rip_ref_ms) } else{ + if (!"align" %in% names(object@proc_params)) { + object@proc_params$align <- list() + } object@proc_params$align$dt_kcorr <- 1 + object@proc_params$align$dt_min_ms <- min(object@drift_time) + object@proc_params$align$dt_max_ms <- max(object@drift_time) } 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") } - # Align in retention time + # Align to injection point in retention time if (align_ip){ 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, ploynomial_order_ptw) - } else if(method_rt == "pow"){ - object <- alignRt_pow(object, ric_ref = ric_ref, ric_ref_rt = ric_ref_rt, ...) } else { if (!"align" %in% names(object@proc_params)) { object@proc_params$align <- list() } - 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") + object@proc_params$align$rt_min_s <- min(object@retention_time) + object@proc_params$align$rt_max_s <- max(object@retention_time) + object@proc_params$align$w <- rep(1, length(object@retention_time)) } dt_range <- c( @@ -72,50 +101,46 @@ methods::setMethod( #' @param rip_ref_ms The position of the RIP in ms #' @return The modified [GCIMSSample] #' @export -methods::setMethod( - "alignDt", "GCIMSSample", - function(object, rip_ref_ms) { - tis <- getTIS(object) - dt <- dtime(object) - rip_pos_ms <- dt[which.max(tis)] - Kcorr <- rip_ref_ms/rip_pos_ms - - int_mat <- intensity(object) - dt_corr <- Kcorr * dt - for (j in seq_len(ncol(int_mat))) { - int_mat[, j] <- signal::interp1(dt_corr, int_mat[, j], dt, extrap = TRUE) - } - intensity(object) <- int_mat - - # Save kcorr and non-extrapolated dt range: - dt_beg <- dt[1] - dt_end <- dt[length(dt)] - dt_corr_beg <- dt_corr[1] - dt_corr_end <- dt_corr[length(dt_corr)] - - if (dt_corr_beg > dt_beg) { - dt_min_ms <- dt[dt >= dt_corr_beg][1] - } else { - dt_min_ms <- dt_beg - } - - if (dt_corr_end < dt_end) { - dt_max_ms <- utils::tail(dt[dt <= dt_corr_end], 1L) - } else { - dt_max_ms <- dt_end - } - - - if (!"align" %in% names(object@proc_params)) { - object@proc_params$align <- list() - } - object@proc_params$align$dt_kcorr <- Kcorr - object@proc_params$align$dt_min_ms <- dt_min_ms - object@proc_params$align$dt_max_ms <- dt_max_ms - # Return the sample - object - }) - +alignDt <- function(object, rip_ref_ms) { + tis <- getTIS(object) + dt <- dtime(object) + rip_pos_ms <- dt[which.max(tis)] + Kcorr <- rip_ref_ms/rip_pos_ms + + int_mat <- intensity(object) + dt_corr <- Kcorr * dt + for (j in seq_len(ncol(int_mat))) { + int_mat[, j] <- signal::interp1(dt_corr, int_mat[, j], dt, extrap = TRUE) + } + intensity(object) <- int_mat + + # Save kcorr and non-extrapolated dt range: + dt_beg <- dt[1] + dt_end <- dt[length(dt)] + dt_corr_beg <- dt_corr[1] + dt_corr_end <- dt_corr[length(dt_corr)] + + if (dt_corr_beg > dt_beg) { + dt_min_ms <- dt[dt >= dt_corr_beg][1] + } else { + dt_min_ms <- dt_beg + } + + if (dt_corr_end < dt_end) { + dt_max_ms <- utils::tail(dt[dt <= dt_corr_end], 1L) + } else { + dt_max_ms <- dt_end + } + + if (!"align" %in% names(object@proc_params)) { + object@proc_params$align <- list() + } + object@proc_params$align$dt_kcorr <- Kcorr + object@proc_params$align$dt_min_ms <- dt_min_ms + object@proc_params$align$dt_max_ms <- dt_max_ms + # Return the sample + object +} #' Align a GCIMSSample in retention time with a multiplicative correction #' @param object A [GCIMSSample] object @@ -123,15 +148,21 @@ methods::setMethod( #' @param rt_ref retention time reference #' @return The modified [GCIMSSample] #' @export -methods::setMethod( - "alignRt_ip","GCIMSSample", - function(object, min_start, rt_ref) { - ric <- getRIC(object) - injection_point <- which.min(ric) - object@retention_time <- rt_ref - object@data <- object@data[, (injection_point - min_start):((injection_point - min_start)+length(rt_ref)-1)] - object - }) +alignRt_ip <- function(object, min_start, rt_ref) { + ric <- getRIC(object) + injection_point <- which.min(ric) + object@retention_time <- rt_ref + object@data <- object@data[, (injection_point - min_start):((injection_point - min_start)+length(rt_ref)-1)] + + if (!"align" %in% names(object@proc_params)) { + object@proc_params$align <- list() + } + + object@proc_params$align$rt_min_s <- rt_ref[1] + object@proc_params$align$rt_max_s <- rt_ref[length(rt_ref)] + object@proc_params$align$w <- rep(1, length(rt_ref)) + object +} #' Align a GCIMSSample in retention time with parametric optimized warping #' @param object A [GCIMSSample] object @@ -143,161 +174,146 @@ methods::setMethod( #' @param lambda1 Regularization parameter for second derivative of warp #' @return The modified [GCIMSSample] #' @export -methods::setMethod( - "alignRt_pow", "GCIMSSample", - function(object, - ric_ref, - ric_ref_rt, - lambdas = pracma::logspace(-2, 4, 31), - p = 10, - max_it = 5000, - lambda1 = 10^6) { - ric <- getRIC(object) - 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", 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, lambda1= lambda1) - - int <- intensity(object) - inter <- t(apply(int, 1, pow::interpolation, w = w, return = FALSE)) - sel <- !is.na(inter[1,]) - object@retention_time <- ric_ref_rt[sel] - object@data <- inter[,sel] - - - if (!"align" %in% names(object@proc_params)) { - object@proc_params$align <- list() - } - object@proc_params$align$warp <- w - #object@proc_params$align$rt_min_s <- object@retention_time[1] - #object@proc_params$align$rt_max_s <- object@retention_time[length(object@retention_time)] - object@proc_params$align$rt_poly_coefs <- c(0, 1) - - return(object) - }) +alignRt_pow <- function(object, + ric_ref, + ric_ref_rt, + lambdas = pracma::logspace(-2, 4, 31), + p = 10, + max_it = 5000, + lambda1 = 10^6) { + ric <- getRIC(object) + 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", 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, lambda1= lambda1) + + int <- intensity(object) + inter <- t(apply(int, 1, pow::interpolation, w = w, return = FALSE)) + sel <- !is.na(inter[1,]) + object@retention_time <- ric_ref_rt[sel] + object@data <- inter[,sel] + + if (!"align" %in% names(object@proc_params)) { + object@proc_params$align <- list() + } + object@proc_params$align$w <- w + object@proc_params$align$rt_min_s <- object@retention_time[1] + object@proc_params$align$rt_max_s <- object@retention_time[length(object@retention_time)] + return(object) +} #' Align a GCIMSSample in retention time using parametric time warping #' @param object 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 +#' @param ploynomial_order maximum order of the polynomial to be used by default 5 #' @return The modified [GCIMSSample] #' @export -methods::setMethod( - "alignRt_ptw", "GCIMSSample", - function(object, ric_ref, ric_ref_rt, ploynomial_order_ptw) { - optimize_polynomial_order <- function(ric_sample, ric_ref) { - 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") - for (j in seq_along(poly_orders)) { - poly <- rep(0, poly_orders[j]) - poly[2] <- 1 - corr[j + 1L] <- stats::cor( - ric_ref, - as.numeric(ptw::ptw(ref = ric_ref, samp = ric_sample, init.coef = poly)$warped.sample[1,]), - use = "complete.obs" - ) - } - - ### RIGHT NOW IF THIS IS UNCOMMENTED 1ST ORDER POLYNOMIAL IS GOING TO BE SELECTED - # # Initialize index: - # idx_max <- idx_sel <- idx_zero <- idx_sign <- length(poly_orders) + 1L - # # Check when correlation decreases for the first time or does not change when increasing degree of the polynomial. - # diff_corr_i <- diff(corr) - # if (all(sign(diff_corr_i) == 1)) { - # # Correlation is always increasing (don't do anything) - # } else if (any(diff_corr_i == 0)) { - # # Correlation is equal for at least for different degrees of the polynomial. - # idx_zero <- which(diff_corr_i == 0)[1] - # } else if (any(sign(diff_corr_i) == -1)) { - # # Correlation decreases at least for two consecutive degrees of the polynomial. - # idx_sign <- which(sign(diff_corr_i) == -1)[1] - # } - # # Combine indexes and look for the minimum - # idx_combine <- c(idx_sign, idx_zero) - # if (any(idx_combine < idx_max)) { - # # Select the minimum index in which the correlation is still increasing - # idx_sel <- min(idx_combine) - # } - correction_type_options[which.max(corr)] +alignRt_ptw <- function(object, ric_ref, ric_ref_rt, ploynomial_order = 5) { + optimize_polynomial_order <- function(ric_sample, ric_ref) { + correction_type_options <- seq.int(0, ploynomial_order) + poly_orders <- seq.int(1, ploynomial_order) + xi <- seq_len(length(ric_sample)) + corr <- numeric(length(poly_orders) + 1L) + corr[1] <- stats::cor(ric_ref, ric_sample, use = "complete.obs") + for (j in seq_along(poly_orders)) { + poly <- rep(0, poly_orders[j]) + poly[2] <- 1 + corr[j + 1L] <- stats::cor( + ric_ref, + as.numeric(ptw::ptw(ref = ric_ref, samp = ric_sample, init.coef = poly)$warped.sample[1,]), + use = "complete.obs" + ) } - # if (missing(ric_ref) && !missing(y)) { - # cli_warn("Please provide ric_ref as a named argument") - # ric_ref <- y + ### RIGHT NOW IF THIS IS UNCOMMENTED 1ST ORDER POLYNOMIAL IS GOING TO BE SELECTED + # # Initialize index: + # idx_max <- idx_sel <- idx_zero <- idx_sign <- length(poly_orders) + 1L + # # Check when correlation decreases for the first time or does not change when increasing degree of the polynomial. + # diff_corr_i <- diff(corr) + # if (all(sign(diff_corr_i) == 1)) { + # # Correlation is always increasing (don't do anything) + # } else if (any(diff_corr_i == 0)) { + # # Correlation is equal for at least for different degrees of the polynomial. + # idx_zero <- which(diff_corr_i == 0)[1] + # } else if (any(sign(diff_corr_i) == -1)) { + # # Correlation decreases at least for two consecutive degrees of the polynomial. + # idx_sign <- which(sign(diff_corr_i) == -1)[1] # } + # # Combine indexes and look for the minimum + # idx_combine <- c(idx_sign, idx_zero) + # if (any(idx_combine < idx_max)) { + # # Select the minimum index in which the correlation is still increasing + # idx_sel <- min(idx_combine) + # } + correction_type_options[which.max(corr)] + } - ric <- getRIC(object) - rt <- rtime(object) - - if (identical(rt, ric_ref_rt)) { - ric_ref_interp <- ric_ref - } else { - ric_ref_interp <- signal::interp1(x = ric_ref_rt, y = ric_ref, xi = rt, extrap = FALSE) - } - - poly_order <- optimize_polynomial_order(ric, ric_ref_interp) - - if (!"align" %in% names(object@proc_params)) { - object@proc_params$align <- list() - } - - - if (poly_order == 0L) { - object@proc_params$align$rt_poly_order <- 2 - object@proc_params$align$rt_poly_coefs <- c(0, 1) - object@proc_params$align$rt_min_s <- rt[1] - object@proc_params$align$rt_max_s <- rt[length(rt)] - return(object) - } - # Correct retention time axis using the reference RIC and sample RIC. - poly <- rep(0, poly_order + 1L) - poly[2] <- 1 - align_result <- ptw::ptw(ref = ric_ref_interp, samp = ric, init.coef = poly) - rt_corr_idx <- align_result$warp.fun[1,] - # Interpolate data using the correction - int_mat <- intensity(object) - rt_idx <- seq_len(ncol(int_mat)) - for (i in seq_len(nrow(int_mat))) { - int_mat[i,] <- signal::interp1(x = rt_corr_idx, y = int_mat[i,], xi = rt_idx, extrap = FALSE) - } - - intensity(object) <- int_mat - - rt_beg <- rt_idx[1] - rt_end <- rt_idx[length(rt_idx)] - rt_corr_beg <- rt_corr_idx[1] - rt_corr_end <- rt_corr_idx[length(rt_corr_idx)] - - if (rt_corr_beg > rt_beg) { - rt_min_idx <- rt_idx[rt_idx >= rt_corr_beg][1] - } else { - rt_min_idx <- rt_beg - } - rt_min_s <- rt[rt_min_idx] + ric <- getRIC(object) + rt <- rtime(object) + if (identical(rt, ric_ref_rt)) { + ric_ref_interp <- ric_ref + } else { + ric_ref_interp <- signal::interp1(x = ric_ref_rt, y = ric_ref, xi = rt, extrap = FALSE) + } - if (rt_corr_end < rt_end) { - rt_max_idx <- utils::tail(rt_idx[rt_idx <= rt_corr_end], 1L) - } else { - rt_max_idx <- rt_end - } - rt_max_s <- rt[rt_max_idx] + poly_order <- optimize_polynomial_order(ric, ric_ref_interp) - object@proc_params$align$rt_poly_order <- poly_order - object@proc_params$align$rt_poly_coefs <- as.numeric(align_result$warp.coef[1L,]) - object@proc_params$align$rt_min_s <- rt_min_s - object@proc_params$align$rt_max_s <- rt_max_s + if (!"align" %in% names(object@proc_params)) { + object@proc_params$align <- list() + } - object - }) + if (poly_order == 0L) { + object@proc_params$align$w <- rep(1,length(rt)) + object@proc_params$align$rt_min_s <- rt[1] + object@proc_params$align$rt_max_s <- rt[length(rt)] + return(object) + } + # Correct retention time axis using the reference RIC and sample RIC. + poly <- rep(0, poly_order + 1L) + poly[2] <- 1 + align_result <- ptw::ptw(ref = ric_ref_interp, samp = ric, init.coef = poly) + rt_corr_idx <- align_result$warp.fun[1,] + # Interpolate data using the correction + int_mat <- intensity(object) + rt_idx <- seq_len(ncol(int_mat)) + for (i in seq_len(nrow(int_mat))) { + int_mat[i,] <- signal::interp1(x = rt_corr_idx, y = int_mat[i,], xi = rt_idx, extrap = FALSE) + } + + intensity(object) <- int_mat + + rt_beg <- rt_idx[1] + rt_end <- rt_idx[length(rt_idx)] + rt_corr_beg <- rt_corr_idx[1] + rt_corr_end <- rt_corr_idx[length(rt_corr_idx)] + + if (rt_corr_beg > rt_beg) { + rt_min_idx <- rt_idx[rt_idx >= rt_corr_beg][1] + } else { + rt_min_idx <- rt_beg + } + rt_min_s <- rt[rt_min_idx] + + + if (rt_corr_end < rt_end) { + rt_max_idx <- utils::tail(rt_idx[rt_idx <= rt_corr_end], 1L) + } else { + rt_max_idx <- rt_end + } + rt_max_s <- rt[rt_max_idx] + + object@proc_params$align$w <- rt_corr_idx + object@proc_params$align$rt_min_s <- rt_min_s + object@proc_params$align$rt_max_s <- rt_max_s + + object +} diff --git a/man/GCIMS-generics.Rd b/man/GCIMS-generics.Rd index 145aca0..cd7d481 100644 --- a/man/GCIMS-generics.Rd +++ b/man/GCIMS-generics.Rd @@ -10,10 +10,7 @@ \alias{filterDt} \alias{decimate} \alias{align} -\alias{alignDt} -\alias{alignRt_ptw} -\alias{alignRt_pow} -\alias{alignRt_ip} +\alias{prealign} \alias{estimateBaseline} \alias{baseline} \alias{baseline<-} @@ -36,13 +33,7 @@ decimate(object, ...) align(object, ...) -alignDt(object, ...) - -alignRt_ptw(object, ...) - -alignRt_pow(object, ...) - -alignRt_ip(object, ...) +prealign(object, ...) estimateBaseline(object, ...) @@ -80,12 +71,6 @@ The object, modified The object, modified -The object, modified - -The object, modified - -The object, modified - The object, with a baseline estimated The baseline of the object @@ -116,13 +101,7 @@ to an existing generics-only package if you need so. \item \code{align()}: Align an object -\item \code{alignDt()}: Align an object in drift time - -\item \code{alignRt_ptw()}: Align an object in retention time - -\item \code{alignRt_pow()}: Align an object in retention time - -\item \code{alignRt_ip()}: Align an object in retention time +\item \code{prealign()}: Align an object in drift time \item \code{estimateBaseline()}: Estimate the baseline in an object diff --git a/man/align-GCIMSDataset-method.Rd b/man/align-GCIMSDataset-method.Rd index c2252ec..c2fb99d 100644 --- a/man/align-GCIMSDataset-method.Rd +++ b/man/align-GCIMSDataset-method.Rd @@ -10,7 +10,6 @@ align_dt = TRUE, align_ip = TRUE, reference_sample_idx = NULL, - ploynomial_order_ptw = 5, ... ) } @@ -18,17 +17,15 @@ \item{object}{A \link{GCIMSDataset} object, modified in-place} \item{method_rt}{Method for alignment, should be "ptw" or "pow" -if pow is selected the package "pow must be installed, to do so visit: +if pow is selected the package "pow" must be installed, to do so visit: https://github.com/sipss/pow} \item{align_dt}{if \code{TRUE} the drift time axis will be aligned using a multiplicative correction} -\item{align_ip}{if TRUE a multiplicative correction will be done in retention time before applying the other algorithm} +\item{align_ip}{if \code{TRUE} a multiplicative correction will be done in retention time before applying the other algorithm} \item{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} -\item{ploynomial_order_ptw}{One number, by default 2, the maximum order of the polynomial for the parametric time warping alignment.} - \item{...}{additional parameters for POW alignment} } \value{ diff --git a/man/align-GCIMSSample-method.Rd b/man/align-GCIMSSample-method.Rd index 5068c4a..0552052 100644 --- a/man/align-GCIMSSample-method.Rd +++ b/man/align-GCIMSSample-method.Rd @@ -2,48 +2,24 @@ % Please edit documentation in R/align-GCIMSSample.R \name{align,GCIMSSample-method} \alias{align,GCIMSSample-method} -\title{Align a GCIMSSample object, in drift and retention time} +\title{Align a GCIMSSample object, in retention time} \usage{ -\S4method{align}{GCIMSSample}( - object, - rip_ref_ms, - ric_ref, - ric_ref_rt, - min_start, - rt_ref, - method_rt, - align_dt, - align_ip, - ploynomial_order_ptw, - ... -) +\S4method{align}{GCIMSSample}(object, method_rt, ric_ref, ric_ref_rt, ...) } \arguments{ \item{object}{A \link{GCIMSSample} object} -\item{rip_ref_ms}{The reference position of the Reactant Ion Peak in the dataset (in ms)} +\item{method_rt}{Method for alignment, should be "ptw" or "pow"} \item{ric_ref}{The reference Reverse Ion Chromatogram} \item{ric_ref_rt}{The retention times corresponding to \code{ric_ref}} -\item{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} - -\item{rt_ref}{retention time reference for alignment to injection point} - -\item{method_rt}{Method for alignment, should be "ptw" or "pow"} - -\item{align_dt}{if \code{TRUE}, align the drift time axis using a multiplicative correction} - -\item{align_ip}{if TRUE a multiplicative correction will be done in retention time before applying the other algorithm} - -\item{ploynomial_order_ptw}{One number, by default 2, the maximum order of the polynomial for the parametric time warping alignment.} - \item{...}{Additional arguments passed on to the alignment method.} } \value{ The modified \link{GCIMSSample} } \description{ -Align a GCIMSSample object, in drift and retention time +Align a GCIMSSample object, in retention time } diff --git a/man/alignDt-GCIMSSample-method.Rd b/man/alignDt.Rd similarity index 77% rename from man/alignDt-GCIMSSample-method.Rd rename to man/alignDt.Rd index 457f1ab..b00b9bb 100644 --- a/man/alignDt-GCIMSSample-method.Rd +++ b/man/alignDt.Rd @@ -1,10 +1,10 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/align-GCIMSSample.R -\name{alignDt,GCIMSSample-method} -\alias{alignDt,GCIMSSample-method} +\name{alignDt} +\alias{alignDt} \title{Align a GCIMSSample in drift time with a multiplicative correction} \usage{ -\S4method{alignDt}{GCIMSSample}(object, rip_ref_ms) +alignDt(object, rip_ref_ms) } \arguments{ \item{object}{A \link{GCIMSSample} object} diff --git a/man/alignRt_ip-GCIMSSample-method.Rd b/man/alignRt_ip.Rd similarity index 79% rename from man/alignRt_ip-GCIMSSample-method.Rd rename to man/alignRt_ip.Rd index 219d5f6..b2d0a99 100644 --- a/man/alignRt_ip-GCIMSSample-method.Rd +++ b/man/alignRt_ip.Rd @@ -1,10 +1,10 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/align-GCIMSSample.R -\name{alignRt_ip,GCIMSSample-method} -\alias{alignRt_ip,GCIMSSample-method} +\name{alignRt_ip} +\alias{alignRt_ip} \title{Align a GCIMSSample in retention time with a multiplicative correction} \usage{ -\S4method{alignRt_ip}{GCIMSSample}(object, min_start, rt_ref) +alignRt_ip(object, min_start, rt_ref) } \arguments{ \item{object}{A \link{GCIMSSample} object} diff --git a/man/alignRt_pow-GCIMSSample-method.Rd b/man/alignRt_pow.Rd similarity index 88% rename from man/alignRt_pow-GCIMSSample-method.Rd rename to man/alignRt_pow.Rd index 4ff9619..0b84fb3 100644 --- a/man/alignRt_pow-GCIMSSample-method.Rd +++ b/man/alignRt_pow.Rd @@ -1,10 +1,10 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/align-GCIMSSample.R -\name{alignRt_pow,GCIMSSample-method} -\alias{alignRt_pow,GCIMSSample-method} +\name{alignRt_pow} +\alias{alignRt_pow} \title{Align a GCIMSSample in retention time with parametric optimized warping} \usage{ -\S4method{alignRt_pow}{GCIMSSample}( +alignRt_pow( object, ric_ref, ric_ref_rt, diff --git a/man/alignRt_ptw-GCIMSSample-method.Rd b/man/alignRt_ptw.Rd similarity index 67% rename from man/alignRt_ptw-GCIMSSample-method.Rd rename to man/alignRt_ptw.Rd index 9c22aa1..64a21bc 100644 --- a/man/alignRt_ptw-GCIMSSample-method.Rd +++ b/man/alignRt_ptw.Rd @@ -1,10 +1,10 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/align-GCIMSSample.R -\name{alignRt_ptw,GCIMSSample-method} -\alias{alignRt_ptw,GCIMSSample-method} +\name{alignRt_ptw} +\alias{alignRt_ptw} \title{Align a GCIMSSample in retention time using parametric time warping} \usage{ -\S4method{alignRt_ptw}{GCIMSSample}(object, ric_ref, ric_ref_rt, ploynomial_order_ptw) +alignRt_ptw(object, ric_ref, ric_ref_rt, ploynomial_order = 5) } \arguments{ \item{object}{A \link{GCIMSSample} object} @@ -13,7 +13,7 @@ \item{ric_ref_rt}{The retention times corresponding to \code{ric_ref}} -\item{ploynomial_order_ptw}{maximum order of the polynomial to be used} +\item{ploynomial_order}{maximum order of the polynomial to be used by default 5} } \value{ The modified \link{GCIMSSample} diff --git a/man/prealign-GCIMSSample-method.Rd b/man/prealign-GCIMSSample-method.Rd new file mode 100644 index 0000000..684caa5 --- /dev/null +++ b/man/prealign-GCIMSSample-method.Rd @@ -0,0 +1,27 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/align-GCIMSSample.R +\name{prealign,GCIMSSample-method} +\alias{prealign,GCIMSSample-method} +\title{Align a GCIMSSample object, in drift time and to the injection point in retention time} +\usage{ +\S4method{prealign}{GCIMSSample}(object, align_dt, align_ip, rip_ref_ms, min_start, rt_ref) +} +\arguments{ +\item{object}{A \link{GCIMSSample} object} + +\item{align_dt}{if \code{TRUE}, align the drift time axis using a multiplicative correction} + +\item{align_ip}{if \code{TRUE} a multiplicative correction will be done in retention time before applying the other algorithm} + +\item{rip_ref_ms}{The reference position of the Reactant Ion Peak in the dataset (in ms)} + +\item{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} + +\item{rt_ref}{retention time reference for alignment to injection point} +} +\value{ +The modified \link{GCIMSSample} +} +\description{ +Align a GCIMSSample object, in drift time and to the injection point in retention time +}