feat: add adjustRtime and related functions for XcmsExperimentHdf5

R/XcmsExperimentHdf5-functions.R

@@ -331,6 +331,47 @@ NULL
                   drop = TRUE)
 }
 
+#' Replace the retention times of chrom peaks with new values, depending
+#' on the provided rts. This function is used during retention time alignment
+#'
+#' @param id `character(1)` with the ID of the sample
+#'
+#' @param rt_old `numeric` with the original retention times
+#'
+#' @param rt_new `numeric` with the new retention times
+#'
+#' @param ms_level `integer` defining for which MS levels the retention times
+#'     should be adjusted. Ideally for all!
+#'
+#' @param hdf5_file `character(1)` with the name of the HDF5 file.
+#'
+#' @return hdf5_count
+#'
+#' @noRd
+.h5_update_rt_chrom_peaks_sample <- function(id, rt_old, rt_new, ms_level,
+                                             hdf5_file) {
+    ## loop over MS levels
+    cnt <- 0L
+    for (msl in ms_level) {
+        ## read chrom peaks
+        cp <- .h5_read_data(hdf5_file, index = id, name = "chrom_peaks",
+                            ms_level = msl, read_colnames = TRUE,
+                            read_rownames = FALSE)[[1L]]
+        ## adjust chrom peak rt - use .applyRtAdjToChromPeaks for that.
+        cp <- .applyRtAdjToChromPeaks(
+            cbind(cp, sample = rep(1, nrow(cp))), rtraw = list(rt_old),
+            rtadj = list(rt_new))
+        l <- list(cp[, colnames(cp) != "sample", drop = FALSE])
+        names(l) <- id
+        ## replace chrom peaks
+        cnt <- .h5_write_data(hdf5_file, data_list = l, "chrom_peaks",
+                              ms_level = msl, replace = FALSE,
+                              write_colnames = FALSE, write_rownames = FALSE)
+    }
+    cnt
+}
+
+
 ################################################################################
 ##
 ##        FEATURES THINGS

R/XcmsExperimentHdf5.R

@@ -45,9 +45,9 @@
 #' @section Retention time alignment:
 #'
 #' - `adjustRtimePeakGroups()` and `adjustRtime()` with `PeakGroupsParam`:
-#'   parameter `extraPeaks` of `PeakGroupsParam` is ignored. Anchor peaks are
-#'   only defined using the `minFraction` parameter (and eventually, if
-#'   provided, the `subset` parameter).
+#'   parameter `extraPeaks` of `PeakGroupsParam` is **ignored**. Anchor peaks
+#'   are thus only defined using the `minFraction` and the optional `subset`
+#'   parameter.
 #'
 #' @section Correspondence analysis results:
 #'
@@ -72,12 +72,14 @@ setClass("XcmsExperimentHdf5",
                    chrom_peaks_ms_level = "integer", # use to keep track of the
                                         # MS levels for which chrom peaks are
                                         # available
+                   gap_peaks_ms_level = "integer", # gap-filled chrom peaks
                    features_ms_level = "integer"),
          prototype = prototype(
              hdf5_file = character(),
              hdf5_mod_count = 0L,
              sample_id = character(),
              chrom_peaks_ms_level = integer(),
+             gap_peaks_ms_level = integer(),
              features_ms_level = integer()
          ))
 
@@ -410,29 +412,26 @@ setMethod(
         rts[order(rowMedians(rts, na.rm = TRUE)), , drop = FALSE]
     })
 
+#' @rdname hidden_aliases
 setMethod(
     "adjustRtime", signature(object = "XcmsExperimentHdf5",
                              param = "PeakGroupsParam"),
     function(object, param, msLevel = 1L, ...) {
         if (hasAdjustedRtime(object))
             stop("Alignment results already present. Please either remove ",
-                 "them with 'dropAdjustedRtime' in order to perform an ",
-                 "alternative, new, alignment, or use 'applyAdjustedRtime'",
-                 " prior 'adjustRtime' to perform a second round of ",
+                 "them with 'dropAdjustedRtime()' in order to perform an ",
+                 "alternative, new, alignment, or use 'applyAdjustedRtime()'",
+                 " prior 'adjustRtime()' to perform a second round of ",
                  "alignment.")
         if (any(msLevel != 1L))
             stop("Alignment is currently only supported for MS level 1")
         if (!nrow(peakGroupsMatrix(param))) {
-            if (!hasFeatures(object))
+            if (!hasFeatures(object, msLevel = msLevel))
                 stop("No feature definitions present in 'object'. Please ",
                      "perform first a correspondence analysis using ",
-                     "'groupChromPeaks'")
-            peakGroupsMatrix(param) <- adjustRtimePeakGroups(
-                object, param = param)
-            ## Need to implement an `adjustRtimePeakGroups,XcmsExperimentHdf5`.
-
+                     "'groupChromPeaks()'")
+            peakGroupsMatrix(param) <- adjustRtimePeakGroups(object, param)
         }
-        ## LLLLL continue here
         fidx <- as.factor(fromFile(object))
         rt_raw <- split(rtime(object), fidx)
         rt_adj <- .adjustRtime_peakGroupsMatrix(
@@ -446,9 +445,12 @@ setMethod(
         else ph <- list()
         object <- dropFeatureDefinitions(object)
         object@spectra$rtime_adjusted <- unlist(rt_adj, use.names = FALSE)
-        ## LLLLL need to fix this here.
-        object@chromPeaks <- .applyRtAdjToChromPeaks(
-            .chromPeaks(object), rtraw = rt_raw, rtadj = rt_adj)
+        res <- mapply(
+            FUN = .h5_update_rt_chrom_peaks_sample,
+            object@sample_id, rt_raw, rt_adj,
+            MoreArgs = list(ms_level = object@chrom_peaks_ms_level,
+                            hdf5_file = object@hdf5_file))
+        object@hdf5_mod_count <- max(unlist(res, use.names = FALSE))
         xph <- XProcessHistory(
             param = param, type. = .PROCSTEP.RTIME.CORRECTION,
             fileIndex = seq_along(object), msLevel = msLevel)
@@ -457,6 +459,38 @@ setMethod(
         object
 })
 
+#' @rdname hidden_aliases
+setMethod("dropAdjustedRtime", "XcmsExperimentHdf5", function(object) {
+    if (!hasAdjustedRtime(object))
+        return(object)
+    ptype <- vapply(object@processHistory, processType, character(1))
+    nom <- length(ptype) + 1L
+    idx_al <- .match_last(.PROCSTEP.RTIME.CORRECTION, ptype, nomatch = nom)
+    idx_co <- .match_last(.PROCSTEP.PEAK.GROUPING, ptype, nomatch = nom)
+    if (hasChromPeaks(object)) {
+        fidx <- as.factor(fromFile(object))
+        res <- mapply(
+            FUN = .h5_update_rt_chrom_peaks_sample,
+            object@sample_id,
+            split(rtime(object, adjusted = TRUE), fidx),
+            split(rtime(object, adjusted = FALSE), fidx),
+            MoreArgs = list(ms_level = object@chrom_peaks_ms_level,
+                            hdf5_file = object@hdf5_file))
+        object@hdf5_mod_count <- max(unlist(res, use.names = FALSE))
+    }
+    svs <- unique(c(spectraVariables(object@spectra), "mz", "intensity"))
+    object@spectra <- selectSpectraVariables(
+        object@spectra, svs[svs != "rtime_adjusted"])
+    object@processHistory <- dropProcessHistoriesList(
+        object@processHistory, type = .PROCSTEP.RTIME.CORRECTION, num = 1L)
+    if (hasFeatures(object) && idx_co > idx_al) {
+        warning("Had to remove feature definitions along with the adjusted ",
+                "retention times because of the dependency between them.")
+        object <- dropFeatureDefinitions(object)
+    }
+    object
+})
+
 ################################################################################
 ##
 ##        CORRESPONDENCE

R/do_adjustRtime-functions.R

@@ -170,8 +170,8 @@ do_adjustRtime_peakGroups <-
             stop("The retention times in 'peakGroupsMatrix' have to be within",
                  " the retention time range of the experiment!")
         rt <- peakGroupsMatrix
-        message("Performing retention time correction using ", nrow(rt),
-                " peak groups.")
+        message("Performing retention time alignment using ", nrow(rt),
+                " anchor peaks.")
 
         ## Calculate the deviation of each peak group in each sample from its
         ## median
@@ -963,4 +963,4 @@ matchedRtimes  <- function(param){
         stop("The inputs need to be of class 'LamaParama'")
     rtMap  <- param@rtMap
     rtMap
-}
+}

R/functions-XCMSnExp.R

@@ -553,7 +553,12 @@ dropGenericProcessHistory <- function(x, fun) {
 }
 
 .hasFilledPeaks <- function(object) {
-    hasChromPeaks(object) & any(chromPeakData(object)$is_filled, na.rm = TRUE)
+    if (is(object, "XcmsExperimentHdf5")) {
+        length(object@gap_peaks_ms_level) > 0
+    } else {
+        hasChromPeaks(object) &
+            any(chromPeakData(object)$is_filled, na.rm = TRUE)
+    }
 }
 
 #' @description

tests/testthat/test_XcmsExperimentHdf5-functions.R

@@ -610,4 +610,9 @@ test_that(".h5_filter works", {
     expect_equal(.h5_filter(), "NONE")
 })
 
+test_that(".h5_update_rt_chrom_peaks_sample works", {
+    ## Unit test is in text_XcmsExperimentHdf5.R @adjustRtime,XcmsExperimentHdf5
+    expect_true(TRUE)
+})
+
 rm(h5f_full_g)

tests/testthat/test_XcmsExperimentHdf5.R

@@ -1,14 +1,15 @@
 h5f <- tempfile()
 xmse_h5 <- xcms:::.xcms_experiment_to_hdf5(loadXcmsData("faahko_sub2"), h5f)
 h5f_full <- tempfile()
-a <- dropFeatureDefinitions(loadXcmsData("xmse"))
+a <- loadXcmsData("xmse") |>
+    dropFeatureDefinitions()
 xmse_full_h5 <- xcms:::.xcms_experiment_to_hdf5(a, h5f_full)
 
 ## correspondence
 h5f_full_g <- tempfile()
 xmseg_full_h5 <- xcms:::.xcms_experiment_to_hdf5(a, h5f_full_g)
 pdp <- PeakDensityParam(sampleGroups = sampleData(xmseg_full_h5)$sample_group,
-                          minFraction = 0.4, bw = 30)
+                        minFraction = 0.4, bw = 30)
 xmseg_full_h5 <- groupChromPeaks(xmseg_full_h5, pdp, msLevel = 1L)
 
 test_that("XcmsExperimentHdf5 validation works", {
@@ -302,10 +303,82 @@ test_that("adjustRtimePeakGroups works", {
     expect_equal(colnames(apeaks_ref), colnames(apeaks))
 })
 
-test_that("adjustRtime,XcmsExperimentHdf5,PeakGroupsParam works", {
-    object <- xmseg_full_h5
-    msLevel <- 1L
-    param <- PeakGroupsParam(span = 0.4)
+test_that("adjustRtime,XcmsExperimentHdf5 and related function work", {
+    ## Note: using `extraPeaks = 100` because that parameter is not supported
+    ## for XcmsExperimentHdf5
+    p <- PeakGroupsParam(span = 0.4, minFraction = 0.7, subset = c(1, 3, 5, 7),
+                         extraPeaks = 100)
+    ## Define the reference data
+    ref <- loadXcmsData("xmse") |>
+        dropFeatureDefinitions() |>
+        applyAdjustedRtime()
+    res_h5 <- tempfile()
+    res <- xcms:::.xcms_experiment_to_hdf5(ref, res_h5)
+    ## Create a single sample XcmsExperimentHdf5
+    a <- ref[3L]
+    a_h5 <- tempfile()
+    a <- xcms:::.xcms_experiment_to_hdf5(a, a_h5)
+    ## Perform retention time alignment on reference data
+    ref <- ref |>
+        groupChromPeaks(pdp, msLevel = 1L) |>
+        adjustRtime(param = p)
+    rt_raw <- rtime(ref, adjusted = FALSE)
+    rt_raw <- split(rt_raw, fromFile(ref))[[3L]]
+    rt_adj <- rtime(ref, adjusted = TRUE)
+    rt_adj <- split(rt_adj, fromFile(ref))[[3L]]
+    cp_raw <- chromPeaks(a)
+
+    ############################################################################
+    ## .h5_update_rt_chrom_peaks_sample: adjust rt of chrom peaks:
+    cnt <- .h5_update_rt_chrom_peaks_sample(
+        a@sample_id[1L], rt_raw, rt_adj, 1L, a@hdf5_file)
+    expect_equal(cnt, a@hdf5_mod_count + 1L)
+    a@hdf5_mod_count <- cnt
+    cp_adj <- chromPeaks(a)
+    expect_true(all(cp_raw[, "rt"] != cp_adj[, "rt"]))
+    expect_true(all(cp_raw[, "rtmin"] != cp_adj[, "rtmin"]))
+    expect_true(all(cp_raw[, "rtmax"] != cp_adj[, "rtmax"]))
+
+    ############################################################################
+    ## adjustRtime: retention time adjustment
+    ## errors
+    expect_error(adjustRtime(xmseg_full_h5, p), "Alignment results already")
+    expect_error(adjustRtime(res, p, msLevel = 2L), "supported for MS level 1")
+    expect_error(adjustRtime(res, p), "No feature definitions present")
+
+    ## Perform alignment
+    res <- groupChromPeaks(res, pdp, msLevel = 1L)
+    expect_false(hasAdjustedRtime(res))
+    cp_ref_raw <- chromPeaks(res)
+    res <- adjustRtime(res, param = p)
+    expect_true(hasAdjustedRtime(res))
+    expect_equal(rtime(res), rtime(ref))
+    expect_equal(unname(chromPeaks(ref)), unname(chromPeaks(res)))
+    expect_true(all(chromPeaks(res)[, "rt"] != cp_ref_raw[, "rt"]))
+    expect_true(validObject(res))
+
+    ############################################################################
+    ## dropAdjustedRtime: revert retention times
+    cp_ref_adj <- chromPeaks(res)
+    cnt <- res@hdf5_mod_count
+    phl <- length(res@processHistory)
+    res <- dropAdjustedRtime(res)
+    expect_true(res@hdf5_mod_count > cnt)
+    expect_true(length(res@processHistory) < phl)
+    expect_false(hasAdjustedRtime(res))
+    expect_true(all(chromPeaks(res)[, "rt"] != cp_ref_adj[, "rt"]))
+    ref <- dropAdjustedRtime(ref)
+    expect_equal(rtime(ref), rtime(res))
+    expect_equal(unname(chromPeaks(ref)), unname(chromPeaks(res)))
+    res <- dropAdjustedRtime(res)
+    expect_false(hasAdjustedRtime(res))
+
+    unlink(a_h5)
+    unlink(res_h5)
+})
+
+test_that(".hasFilledPeaks works with XcmsExperimentHdf5", {
+    expect_false(.hasFilledPeaks(xmse_h5))
 })
 
 ## test_that(".h5_feature_chrom_peaks_sample works", {