From 5147dc71fcce883866df17b5ab5cec9e40d5e390 Mon Sep 17 00:00:00 2001 From: RogerGinBer <45530914+RogerGinBer@users.noreply.github.com> Date: Fri, 2 Dec 2022 10:38:05 +0100 Subject: [PATCH 1/8] Add new Ion-mobility peak picking algorithm Add a new function for ion-mobility peak picking, compatible with `MsExperiment`'s `Spectra` only (`TimsTofBackend`), with a corresponding CentWaveParamIM class: - Create `do_findChromPeaks_IM_centWave`, the API for the peakpicking itself - Create a `IMParam`class that inherits from `Param` and that will serve as superclass for all ion-mobility peak picking Params. - Create a `CentWaveParamIM` that inherits from both `CentWaveParam` and `IMParam` - Add a function dispatch point in .mse_find_chrom_peaks_sample - Create documentation for both API and param class --- NAMESPACE | 1 + R/DataClasses.R | 60 +++++ R/MsExperiment-functions.R | 14 +- R/do_findChromPeaks-functions.R | 211 ++++++++++++++++++ R/functions-Params.R | 27 +++ man/do_findChromPeaks_IM_centWave.Rd | 152 +++++++++++++ man/do_findChromPeaks_centWave.Rd | 1 + ..._findChromPeaks_centWaveWithPredIsoROIs.Rd | 1 + man/do_findChromPeaks_massifquant.Rd | 1 + man/do_findChromPeaks_matchedFilter.Rd | 1 + man/do_findPeaks_MSW.Rd | 1 + man/findChromPeaks-centWave.Rd | 1 + man/findChromPeaks-centWaveIonMobility.Rd | 157 +++++++++++++ man/findChromPeaks-centWaveWithPredIsoROIs.Rd | 1 + man/findChromPeaks-massifquant.Rd | 1 + man/findChromPeaks-matchedFilter.Rd | 1 + man/findChromPeaks.Rd | 1 + man/findPeaks-MSW.Rd | 1 + 18 files changed, 629 insertions(+), 4 deletions(-) create mode 100644 man/do_findChromPeaks_IM_centWave.Rd create mode 100644 man/findChromPeaks-centWaveIonMobility.Rd diff --git a/NAMESPACE b/NAMESPACE index 0061bdc70..2dbb0db9a 100755 --- a/NAMESPACE +++ b/NAMESPACE @@ -447,6 +447,7 @@ export("CentWaveParam", "MassifquantParam", "MSWParam", "CentWavePredIsoParam", + "CentWaveParamIM", "PeakDensityParam", "MzClustParam", "NearestPeaksParam", diff --git a/R/DataClasses.R b/R/DataClasses.R index 810abc1e9..941e3616c 100644 --- a/R/DataClasses.R +++ b/R/DataClasses.R @@ -1298,6 +1298,66 @@ setClass("CentWavePredIsoParam", else TRUE }) +#### Ion mobility peak-picking classes #### + +setClass("IMParam", contains = "Param") + +#' @title Centwave-based ion-mobility peak picking +#' +#' @aliases centWaveIonMobility +#' +#' @description Performs an extension of CentWave peak-picking on LC-IM-MS MS1 +#' data: first it joins all mobility scans into frames and performs .centWave_orig on +#' the summarized LC-MS-like data; then, from each peak, it calculates its mobilogram and +#' performs a second peak-picking on the IM dimension, resolving the peaks. +#' +#' @inheritParams findChromPeaks-centWave +#' +#' @param ppmMerging The maximum mass deviation allowed when grouping individual +#' IM scans into frames. Data points within \code{ppmMerging} ppm will be +#' summed up into a single value and the reported mz will be their weighted +#' average. +#' +#' @param binWidthIM The bin size used when calculating the mobilograms to resolve +#' the peaks into the ion-mobility dimension. Lower values will give better resolution +#' if the data allows it, but can also generate spurious peaks. +#' +#' @details See \code{\link{centWave}} for details on the centWave method. +#' +#' @family peak detection methods +#' +#' @author Roger Gine, Johannes Rainer +#' +#' @seealso The \code{\link{do_findChromPeaks_IM_centWave}} core +#' API function and \code{\link{CentWaveParam}} for the class the +#' \code{CentWaveParamIM} extends. +#' +#' @name findChromPeaks-centWaveIonMobility +NULL + +#' @description The \code{CentWaveParamIM} class allows to specify all +#' settings for +#' Instances should be created with the \code{CentWaveParamIM} +#' constructor. See also the documentation of the +#' \code{\link{CentWaveParam}} for all methods and arguments this class +#' inherits. +#' +#' @slot ppm,peakwidth,snthresh,prefilter,mzCenterFun,integrate,mzdiff,fitgauss,noise,verboseColumns,roiList,firstBaselineCheck,roiScales,ppmMerging,binWidthIM +#' See corresponding parameter above. +#' +#' @rdname findChromPeaks-centWaveIonMobility +setClass("CentWaveParamIM", + contains = c("IMParam", "CentWaveParam"), + slots = c( + ppmMerging = "numeric", + binWidthIM = "numeric" + + ), + prototype = prototype( + ppmMerging = 10, + binWidthIM = 0.02 + )) + setClass("PeakDensityParam", slots = c(sampleGroups = "ANY", bw = "numeric", diff --git a/R/MsExperiment-functions.R b/R/MsExperiment-functions.R index db34ebb33..0471dbba0 100644 --- a/R/MsExperiment-functions.R +++ b/R/MsExperiment-functions.R @@ -3,7 +3,8 @@ MatchedFilterParam = "do_findChromPeaks_matchedFilter", MassifquantParam = "do_findChromPeaks_massifquant", MSWParam = "do_findPeaks_MSW", - CentWavePredIsoParam = "do_findChromPeaks_centWaveWithPredIsoROIs") + CentWavePredIsoParam = "do_findChromPeaks_centWaveWithPredIsoROIs", + CentWaveParamIM = "do_findChromPeaks_IM_centWave") fun <- p2f[class(x)[1L]] if (is.na(fun)) stop("No peak detection function for parameter class ", class(x)[1L]) @@ -62,6 +63,11 @@ #' @noRd .mse_find_chrom_peaks_sample <- function(x, msLevel = 1L, param, ...) { x <- filterMsLevel(x, msLevel) + if(inherits(param, "IMParam")){ + if(!any(c("inv_ion_mobility") %in% Spectra::spectraVariables(x))) # Add any other column name needed + stop("Your Spectra object doesn't contain ion-mobility data") + return(do.call(.param_to_fun(param), args = append(list(x), as(param, "list")))) #Append to avoid concatenating spectra + } pkd <- Spectra::peaksData(x, columns = c("mz", "intensity"), BPPARAM = SerialParam()) vals_per_spect <- vapply(pkd, nrow, integer(1), USE.NAMES = FALSE) @@ -74,6 +80,9 @@ pkd <- do.call(rbind, pkd) if (!length(pkd)) return(NULL) # not returning matrix because of rbind + rts <- rtime(x) + if (is.unsorted(rts)) + stop("Spectra are not ordered by retention time", .call = FALSE) if (inherits(param, "CentWaveParam")) { centroided <- all(centroided(x)) if (is.na(centroided)) { @@ -83,9 +92,6 @@ " works best on data in centroid mode.") } } - rts <- rtime(x) - if (is.unsorted(rts)) - stop("Spectra are not ordered by retention time", .call = FALSE) do.call(.param_to_fun(param), args = c(list(mz = pkd[, 1L], int = pkd[, 2L], scantime = rts, valsPerSpect = vals_per_spect), as(param, "list"))) diff --git a/R/do_findChromPeaks-functions.R b/R/do_findChromPeaks-functions.R index b6d8d5e22..0257460bb 100644 --- a/R/do_findChromPeaks-functions.R +++ b/R/do_findChromPeaks-functions.R @@ -2220,9 +2220,220 @@ do_findPeaks_MSW <- function(mz, int, snthresh = 3, peaklist } +############################################################ +## Ion-mobility peak picking +## +#' @title Core API for Centwave-based ion-mobility peak picking +#' @name do_findChromPeaks_IM_centWave +#' +#' @description Performs an extension of CentWave peak-picking on LC-IM-MS MS1 +#' data. First it joins all scans into frames and performs .centWave_orig on +#' the summarized LC-MS-like data. From each peak, it calculates its mobilogram and +#' performs a second peak-picking on the IM dimension, resolving the peaks. +#' +#' @inheritParams do_findChromPeaks_centWave +#' @inheritParams findChromPeaks-centWaveIonMobility +#' +#' @return A matrix, each row representing an identified peak, with columns: +#' \describe{ +#' \item{mz}{m/z value of the peak at the apex position.} +#' \item{mzmin}{Minimum m/z of the peak.} +#' \item{mzmax}{Maximum m/z of the peak.} +#' \item{rt}{Retention time value of the peak at the apex position.} +#' \item{rtmin}{Minimum retention time of the peak.} +#' \item{rtmax}{Maximum retention time of the peak.} +#' \item{im}{Ion mobility value of the peak at the apex position.} +#' \item{immin}{Minimum ion mobility value of the peak.} +#' \item{immax}{Maximum ion mobility value of the peak.} +#' \item{maxo}{Maximum intensity of the peak.} +#' \item{into}{Integrated (original) intensity of the peak.} +#' \item{intb}{Always \code{NA}.} +#' \item{sn}{Always \code{NA}} +#' } +#' +#' @family core peak detection functions +#' +#' @author Roger Gine, Johannes Rainer +#' +#' @importFrom Spectra peaksData rtime combineSpectra mz +do_findChromPeaks_IM_centWave <- function(spec, + ppm = 25, + peakwidth = c(20, 50), + snthresh = 10, + prefilter = c(3, 100), + mzCenterFun = "wMean", + integrate = 1, + mzdiff = -0.001, + fitgauss = FALSE, + noise = 0, + verboseColumns = FALSE, + roiList = list(), + firstBaselineCheck = TRUE, + roiScales = NULL, + sleep = 0, + extendLengthMSW = FALSE, + ppmMerging = 10, + binWidthIM = 0.01 + ){ + ## Extract frame information + pdata <- peaksData(spec, columns = c("mz", "intensity")) + rt <- rtime(spec) + im <- spec$inv_ion_mobility + + + ## Merging frames into scans and Summarize across IM dimension + message("Collapsing data over IM dimension... ", appendLF = F) + + scans_summarized <- + combineSpectra( + spec, + f = as.factor(spec$frameId), + intensityFun = base::sum, + weighted = TRUE, + ppm = ppmMerging + ) + message("OK") + + ## Peak-picking on summarized data + mzs <- mz(scans_summarized) + valsPerSpect <- lengths(mzs, FALSE) + mz <- unlist(mzs, use.names = FALSE) + int <- unlist(intensity(scans_summarized), use.names = FALSE) + scantime <- sort(unique(rt)) + peaks <- .centWave_orig(mz = mz, int = int, scantime = scantime, + valsPerSpect = valsPerSpect, ppm = ppm, peakwidth = peakwidth, + snthresh = snthresh, prefilter = prefilter, + mzCenterFun = mzCenterFun, integrate = integrate, + mzdiff = mzdiff, fitgauss = fitgauss, noise = noise, + verboseColumns = verboseColumns, roiList = roiList, + firstBaselineCheck = firstBaselineCheck, + roiScales = roiScales, sleep = sleep, + extendLengthMSW = extendLengthMSW) + + ## Resolving peaks across IM dimension + message("Resolving peaks over ion-mobility dimension... ", appendLF = F) + resolved_peaks <- vector("list", nrow(peaks)) + for (i in seq_len(nrow(peaks))) { + current_peak <- peaks[i,] + mobilogram <- .extract_mobilogram(pdata, current_peak, rt, im, binWidthIM) + if (length(mobilogram) == 0) { + warning(i, " mobilogram is empty") + next + } + bounds <- .split_mobilogram(mobilogram) + new_peaks <- data.frame( + mz = current_peak["mz"], + mzmin = current_peak["mzmin"], + mzmax = current_peak["mzmax"], + rt = current_peak["rt"], + rtmin = current_peak["rtmin"], + rtmax = current_peak["rtmax"], + im = vapply(bounds, mean, numeric(1)), + immin = vapply(bounds, min, numeric(1)), + immax = vapply(bounds, max, numeric(1)), + row.names = NULL + ) + resolved_peaks[[i]] <- new_peaks + } + resolved_peaks <- do.call(rbind, resolved_peaks) + message("OK") + + ## Refine and calculate peak parameters + vals <- vector("list", nrow(resolved_peaks)) + for (i in seq(nrow(resolved_peaks))) { + peak <- unlist(resolved_peaks[i, , drop = TRUE]) + + ## Create a EIC for mz, rt and IM ranges + eic <- .extract_EIC_IM(peak, pdata, rt, im) + + if (nrow(eic) == 0 | all(eic[, 2] == 0)) + next + + ## Refine RT bounds + rts <- c(peak["rtmin"], peak["rtmax"]) + apx <- which.max(eic[, 2]) + apx_rt <- eic[apx, 1] + range <- xcms:::descendMin(eic[, 2], apx) + + eic <- eic[range[1]:range[2], , drop = FALSE] + + ## Calculate peak stats + vals[[i]] <- data.frame( + mz = peak["mz"], + mzmin = peak["mzmin"], + mzmax = peak["mzmax"], + rt = apx_rt, + rtmin = min(eic[, 1]), + rtmax = max(eic[, 1]), + im = peak["im"], + immin = peak["immin"], + immax = peak["immax"], + maxo = max(eic[, 2]), + into = sum(eic[, 2]), + intb = NA, + sn = NA + ) + } + resolved_peaks <- do.call(rbind, vals) + resolved_peaks <- + resolved_peaks[resolved_peaks$into > 0, ] #Remove empty peaks + + as.matrix(resolved_peaks) +} + +#' @importFrom MsCoreUtils bin +.extract_mobilogram <- function(pdata, peak, rt, im, binWidthIM = 0.01){ + rtr <- c(peak["rtmin"], peak["rtmax"]) + mzr <- c(peak["mzmin"], peak["mzmax"]) + keep <- dplyr::between(rt, rtr[1], rtr[2]) + if (length(keep) == 0) return() + ims <- im[keep] + ints <- vapply(pdata[keep], xcms:::.aggregate_intensities, + mzr = mzr, INTFUN = sum, na.rm = TRUE, numeric(1)) + if(all(ints == 0)) return() + mob <- MsCoreUtils::bin(x = ints, y = ims, size = binWidthIM, FUN = sum) + mob +} + + +#' @importFrom MassSpecWavelet peakDetectionCWT +.split_mobilogram <- function(mob){ + if(length(mob$x) == 0){return()} + vec <- mob$x + #Add some padding, which will be removed after + padding_size <- 5 + vec <- c(rep(0, padding_size), vec, rep(0, padding_size)) + pks <- MassSpecWavelet::peakDetectionCWT(vec, scales = c(1:7)) + left <- sapply(pks$majorPeakInfo$peakCenterIndex - pks$majorPeakInfo$peakScale, function(x) max(1, x)) + right <- sapply(pks$majorPeakInfo$peakCenterIndex + pks$majorPeakInfo$peakScale, function(x) min(x, length(vec) - padding_size)) + limits <- list() + for (i in seq_along(pks$majorPeakInfo$peakCenterIndex)){ + ranges <- xcms:::descendMinTol(vec, startpos = c(left[i], right[i]), maxDescOutlier = 1) - padding_size + ranges[1] <- min(max(1, ranges[1]), length(mob$mids)) + ranges[2] <- min(ranges[2], length(mob$mids)) + limits[[i]] <- mob$mids[ranges] + } + limits <- limits[vapply(limits, function(x){!any(is.na(x))}, logical(1))] + return(limits) +} +#' @importFrom dplyr between +.extract_EIC_IM <- function(peak, pdata, rt, im){ + rtr <- c(peak["rtmin"], peak["rtmax"]) + mzr <- c(peak["mzmin"], peak["mzmax"]) + imr <- c(peak["immin"], peak["immax"]) + + keep <- dplyr::between(rt, rtr[1], rtr[2]) & dplyr::between(im, imr[1], imr[2]) + rts <- rt[keep] + ints <- vapply(pdata[keep], xcms:::.aggregate_intensities, + mzr = mzr, INTFUN = sum, na.rm = TRUE, numeric(1)) + ints <- vapply(unique(rts), function(x){sum(ints[rts == x])}, numeric(1)) + + cbind(unique(rts), ints) +} + ############################################################ ## MS1 diff --git a/R/functions-Params.R b/R/functions-Params.R index 24e1f7908..e8167a309 100644 --- a/R/functions-Params.R +++ b/R/functions-Params.R @@ -226,6 +226,33 @@ CentWavePredIsoParam <- function(ppm = 25, peakwidth = c(20, 50), snthresh = 10, mzIntervalExtension = mzIntervalExtension, polarity = polarity)) } + +#' @return The \code{CentWaveParamIM} function returns a +#' \code{CentWaveParamIM} class instance with all of the settings +#' specified for the centWave-based peak detection in chromatographic + +#' ion mobility data. +#' +#' @rdname findChromPeaks-centWaveIonMobility +#' +CentWaveParamIM <- function(ppm = 25, peakwidth = c(20, 50), snthresh = 10, + prefilter = c(3, 100), mzCenterFun = "wMean", + integrate = 1L, mzdiff = -0.001, fitgauss = FALSE, + noise = 0, verboseColumns = FALSE, roiList = list(), + firstBaselineCheck = TRUE, roiScales = numeric(), + extendLengthMSW = FALSE, ppmMerging = 10, + binWidthIM = 0.02) { + return(new("CentWaveParamIM", ppm = ppm, peakwidth = peakwidth, + snthresh = snthresh, prefilter = prefilter, + mzCenterFun = mzCenterFun, integrate = as.integer(integrate), + mzdiff = mzdiff, fitgauss = fitgauss, noise = noise, + verboseColumns = verboseColumns, roiList = roiList, + firstBaselineCheck = firstBaselineCheck, roiScales = roiScales, + extendLengthMSW = extendLengthMSW, ppmMerging = ppmMerging, + binWidthIM = binWidthIM)) +} + + + #' @rdname groupChromPeaks PeakDensityParam <- function(sampleGroups = numeric(), bw = 30, minFraction = 0.5, minSamples = 1, diff --git a/man/do_findChromPeaks_IM_centWave.Rd b/man/do_findChromPeaks_IM_centWave.Rd new file mode 100644 index 000000000..d935de6c2 --- /dev/null +++ b/man/do_findChromPeaks_IM_centWave.Rd @@ -0,0 +1,152 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/do_findChromPeaks-functions.R +\name{do_findChromPeaks_IM_centWave} +\alias{do_findChromPeaks_IM_centWave} +\title{Core API for Centwave-based ion-mobility peak picking} +\usage{ +do_findChromPeaks_IM_centWave( + spec, + ppm = 25, + peakwidth = c(20, 50), + snthresh = 10, + prefilter = c(3, 100), + mzCenterFun = "wMean", + integrate = 1, + mzdiff = -0.001, + fitgauss = FALSE, + noise = 0, + verboseColumns = FALSE, + roiList = list(), + firstBaselineCheck = TRUE, + roiScales = NULL, + sleep = 0, + extendLengthMSW = FALSE, + ppmMerging = 10, + binWidthIM = 0.01 +) +} +\arguments{ +\item{ppm}{\code{numeric(1)} defining the maximal tolerated m/z deviation in +consecutive scans in parts per million (ppm) for the initial ROI +definition.} + +\item{peakwidth}{\code{numeric(2)} with the expected approximate +peak width in chromatographic space. Given as a range (min, max) +in seconds.} + +\item{snthresh}{\code{numeric(1)} defining the signal to noise ratio cutoff.} + +\item{prefilter}{\code{numeric(2)}: \code{c(k, I)} specifying the prefilter +step for the first analysis step (ROI detection). Mass traces are only +retained if they contain at least \code{k} peaks with intensity +\code{>= I}.} + +\item{mzCenterFun}{Name of the function to calculate the m/z center of the +chromatographic peak. Allowed are: \code{"wMean"}: intensity weighted +mean of the peak's m/z values, \code{"mean"}: mean of the peak's m/z +values, \code{"apex"}: use the m/z value at the peak apex, +\code{"wMeanApex3"}: intensity weighted mean of the m/z value at the +peak apex and the m/z values left and right of it and \code{"meanApex3"}: +mean of the m/z value of the peak apex and the m/z values left and right +of it.} + +\item{integrate}{Integration method. For \code{integrate = 1} peak limits +are found through descent on the mexican hat filtered data, for +\code{integrate = 2} the descent is done on the real data. The latter +method is more accurate but prone to noise, while the former is more +robust, but less exact.} + +\item{mzdiff}{\code{numeric(1)} representing the minimum difference in m/z +dimension required for peaks with overlapping retention times; can be +negative to allow overlap. During peak post-processing, peaks +defined to be overlapping are reduced to the one peak with the largest +signal.} + +\item{fitgauss}{\code{logical(1)} whether or not a Gaussian should be fitted +to each peak. This affects mostly the retention time position of the +peak.} + +\item{noise}{\code{numeric(1)} allowing to set a minimum intensity required +for centroids to be considered in the first analysis step (centroids with +intensity \code{< noise} are omitted from ROI detection).} + +\item{verboseColumns}{\code{logical(1)} whether additional peak meta data +columns should be returned.} + +\item{roiList}{An optional list of regions-of-interest (ROI) representing +detected mass traces. If ROIs are submitted the first analysis step is +omitted and chromatographic peak detection is performed on the submitted +ROIs. Each ROI is expected to have the following elements specified: +\code{scmin} (start scan index), \code{scmax} (end scan index), +\code{mzmin} (minimum m/z), \code{mzmax} (maximum m/z), \code{length} +(number of scans), \code{intensity} (summed intensity). Each ROI should +be represented by a \code{list} of elements or a single row +\code{data.frame}.} + +\item{firstBaselineCheck}{\code{logical(1)}. If \code{TRUE} continuous +data within regions of interest is checked to be above the first baseline. +In detail, a first rough estimate of the noise is calculated and peak +detection is performed only in regions in which multiple sequential +signals are higher than this first estimated baseline/noise level.} + +\item{roiScales}{Optional numeric vector with length equal to \code{roiList} +defining the scale for each region of interest in \code{roiList} that +should be used for the centWave-wavelets.} + +\item{sleep}{\code{numeric(1)} defining the number of seconds to wait between +iterations. Defaults to \code{sleep = 0}. If \code{> 0} a plot is +generated visualizing the identified chromatographic peak. Note: this +argument is for backward compatibility only and will be removed in +future.} + +\item{extendLengthMSW}{Option to force centWave to use all scales when +running centWave rather than truncating with the EIC length. Uses the "open" +method to extend the EIC to a integer base-2 length prior to being passed to +\code{convolve} rather than the default "reflect" method. See +https://github.com/sneumann/xcms/issues/445 for more information.} + +\item{ppmMerging}{The maximum mass deviation allowed when grouping individual +IM scans into frames. Data points within \code{ppmMerging} ppm will be +summed up into a single value and the reported mz will be their weighted +average.} + +\item{binWidthIM}{The bin size used when calculating the mobilograms to resolve +the peaks into the ion-mobility dimension. Lower values will give better resolution +if the data allows it, but can also generate spurious peaks.} +} +\value{ +A matrix, each row representing an identified peak, with columns: + \describe{ + \item{mz}{m/z value of the peak at the apex position.} + \item{mzmin}{Minimum m/z of the peak.} + \item{mzmax}{Maximum m/z of the peak.} + \item{rt}{Retention time value of the peak at the apex position.} + \item{rtmin}{Minimum retention time of the peak.} + \item{rtmax}{Maximum retention time of the peak.} + \item{im}{Ion mobility value of the peak at the apex position.} + \item{immin}{Minimum ion mobility value of the peak.} + \item{immax}{Maximum ion mobility value of the peak.} + \item{maxo}{Maximum intensity of the peak.} + \item{into}{Integrated (original) intensity of the peak.} + \item{intb}{Always \code{NA}.} + \item{sn}{Always \code{NA}} + } +} +\description{ +Performs an extension of CentWave peak-picking on LC-IM-MS MS1 + data. First it joins all scans into frames and performs .centWave_orig on + the summarized LC-MS-like data. From each peak, it calculates its mobilogram and + performs a second peak-picking on the IM dimension, resolving the peaks. +} +\seealso{ +Other core peak detection functions: +\code{\link{do_findChromPeaks_centWaveWithPredIsoROIs}()}, +\code{\link{do_findChromPeaks_centWave}()}, +\code{\link{do_findChromPeaks_massifquant}()}, +\code{\link{do_findChromPeaks_matchedFilter}()}, +\code{\link{do_findPeaks_MSW}()} +} +\author{ +Roger Gine, Johannes Rainer +} +\concept{core peak detection functions} diff --git a/man/do_findChromPeaks_centWave.Rd b/man/do_findChromPeaks_centWave.Rd index 663b557d1..e2f8a95e9 100644 --- a/man/do_findChromPeaks_centWave.Rd +++ b/man/do_findChromPeaks_centWave.Rd @@ -214,6 +214,7 @@ Ralf Tautenhahn, Christoph B\"{o}ttcher, and Steffen Neumann "Highly \code{\link{centWave}} for the standard user interface method. Other core peak detection functions: +\code{\link{do_findChromPeaks_IM_centWave}()}, \code{\link{do_findChromPeaks_centWaveWithPredIsoROIs}()}, \code{\link{do_findChromPeaks_massifquant}()}, \code{\link{do_findChromPeaks_matchedFilter}()}, diff --git a/man/do_findChromPeaks_centWaveWithPredIsoROIs.Rd b/man/do_findChromPeaks_centWaveWithPredIsoROIs.Rd index 0591e8fad..e5a88b6f6 100644 --- a/man/do_findChromPeaks_centWaveWithPredIsoROIs.Rd +++ b/man/do_findChromPeaks_centWaveWithPredIsoROIs.Rd @@ -222,6 +222,7 @@ For more details on the centWave algorithm see } \seealso{ Other core peak detection functions: +\code{\link{do_findChromPeaks_IM_centWave}()}, \code{\link{do_findChromPeaks_centWave}()}, \code{\link{do_findChromPeaks_massifquant}()}, \code{\link{do_findChromPeaks_matchedFilter}()}, diff --git a/man/do_findChromPeaks_massifquant.Rd b/man/do_findChromPeaks_massifquant.Rd index 0f0b4b1d2..88530dc95 100644 --- a/man/do_findChromPeaks_massifquant.Rd +++ b/man/do_findChromPeaks_massifquant.Rd @@ -204,6 +204,7 @@ detection" \emph{Bioinformatics} 2014, 30(18):2636-43. \code{\link{massifquant}} for the standard user interface method. Other core peak detection functions: +\code{\link{do_findChromPeaks_IM_centWave}()}, \code{\link{do_findChromPeaks_centWaveWithPredIsoROIs}()}, \code{\link{do_findChromPeaks_centWave}()}, \code{\link{do_findChromPeaks_matchedFilter}()}, diff --git a/man/do_findChromPeaks_matchedFilter.Rd b/man/do_findChromPeaks_matchedFilter.Rd index 9932925b6..b362b6f2d 100644 --- a/man/do_findChromPeaks_matchedFilter.Rd +++ b/man/do_findChromPeaks_matchedFilter.Rd @@ -166,6 +166,7 @@ Profiling Using Nonlinear Peak Alignment, Matching, and Identification" \code{\link{matchedFilter}} for the standard user interface method. Other core peak detection functions: +\code{\link{do_findChromPeaks_IM_centWave}()}, \code{\link{do_findChromPeaks_centWaveWithPredIsoROIs}()}, \code{\link{do_findChromPeaks_centWave}()}, \code{\link{do_findChromPeaks_massifquant}()}, diff --git a/man/do_findPeaks_MSW.Rd b/man/do_findPeaks_MSW.Rd index e918b25ba..ac34ab611 100644 --- a/man/do_findPeaks_MSW.Rd +++ b/man/do_findPeaks_MSW.Rd @@ -55,6 +55,7 @@ This is a wrapper around the peak picker in Bioconductor's \code{MassSpecWavelet} package. Other core peak detection functions: +\code{\link{do_findChromPeaks_IM_centWave}()}, \code{\link{do_findChromPeaks_centWaveWithPredIsoROIs}()}, \code{\link{do_findChromPeaks_centWave}()}, \code{\link{do_findChromPeaks_massifquant}()}, diff --git a/man/findChromPeaks-centWave.Rd b/man/findChromPeaks-centWave.Rd index 8157aa4b3..f076cbbef 100644 --- a/man/findChromPeaks-centWave.Rd +++ b/man/findChromPeaks-centWave.Rd @@ -385,6 +385,7 @@ detection in purely chromatographic data. the peak detection. Other peak detection methods: +\code{\link{findChromPeaks-centWaveIonMobility}}, \code{\link{findChromPeaks-centWaveWithPredIsoROIs}}, \code{\link{findChromPeaks-massifquant}}, \code{\link{findChromPeaks-matchedFilter}}, diff --git a/man/findChromPeaks-centWaveIonMobility.Rd b/man/findChromPeaks-centWaveIonMobility.Rd new file mode 100644 index 000000000..5eea3ace1 --- /dev/null +++ b/man/findChromPeaks-centWaveIonMobility.Rd @@ -0,0 +1,157 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/DataClasses.R, R/functions-Params.R +\docType{class} +\name{findChromPeaks-centWaveIonMobility} +\alias{findChromPeaks-centWaveIonMobility} +\alias{centWaveIonMobility} +\alias{CentWaveParamIM-class} +\alias{CentWaveParamIM} +\title{Centwave-based ion-mobility peak picking} +\usage{ +CentWaveParamIM( + ppm = 25, + peakwidth = c(20, 50), + snthresh = 10, + prefilter = c(3, 100), + mzCenterFun = "wMean", + integrate = 1L, + mzdiff = -0.001, + fitgauss = FALSE, + noise = 0, + verboseColumns = FALSE, + roiList = list(), + firstBaselineCheck = TRUE, + roiScales = numeric(), + extendLengthMSW = FALSE, + ppmMerging = 10, + binWidthIM = 0.02 +) +} +\arguments{ +\item{ppm}{\code{numeric(1)} defining the maximal tolerated m/z deviation in +consecutive scans in parts per million (ppm) for the initial ROI +definition.} + +\item{peakwidth}{\code{numeric(2)} with the expected approximate +peak width in chromatographic space. Given as a range (min, max) +in seconds.} + +\item{snthresh}{\code{numeric(1)} defining the signal to noise ratio cutoff.} + +\item{prefilter}{\code{numeric(2)}: \code{c(k, I)} specifying the prefilter +step for the first analysis step (ROI detection). Mass traces are only +retained if they contain at least \code{k} peaks with intensity +\code{>= I}.} + +\item{mzCenterFun}{Name of the function to calculate the m/z center of the +chromatographic peak. Allowed are: \code{"wMean"}: intensity weighted +mean of the peak's m/z values, \code{"mean"}: mean of the peak's m/z +values, \code{"apex"}: use the m/z value at the peak apex, +\code{"wMeanApex3"}: intensity weighted mean of the m/z value at the +peak apex and the m/z values left and right of it and \code{"meanApex3"}: +mean of the m/z value of the peak apex and the m/z values left and right +of it.} + +\item{integrate}{Integration method. For \code{integrate = 1} peak limits +are found through descent on the mexican hat filtered data, for +\code{integrate = 2} the descent is done on the real data. The latter +method is more accurate but prone to noise, while the former is more +robust, but less exact.} + +\item{mzdiff}{\code{numeric(1)} representing the minimum difference in m/z +dimension required for peaks with overlapping retention times; can be +negative to allow overlap. During peak post-processing, peaks +defined to be overlapping are reduced to the one peak with the largest +signal.} + +\item{fitgauss}{\code{logical(1)} whether or not a Gaussian should be fitted +to each peak. This affects mostly the retention time position of the +peak.} + +\item{noise}{\code{numeric(1)} allowing to set a minimum intensity required +for centroids to be considered in the first analysis step (centroids with +intensity \code{< noise} are omitted from ROI detection).} + +\item{verboseColumns}{\code{logical(1)} whether additional peak meta data +columns should be returned.} + +\item{roiList}{An optional list of regions-of-interest (ROI) representing +detected mass traces. If ROIs are submitted the first analysis step is +omitted and chromatographic peak detection is performed on the submitted +ROIs. Each ROI is expected to have the following elements specified: +\code{scmin} (start scan index), \code{scmax} (end scan index), +\code{mzmin} (minimum m/z), \code{mzmax} (maximum m/z), \code{length} +(number of scans), \code{intensity} (summed intensity). Each ROI should +be represented by a \code{list} of elements or a single row +\code{data.frame}.} + +\item{firstBaselineCheck}{\code{logical(1)}. If \code{TRUE} continuous +data within regions of interest is checked to be above the first baseline. +In detail, a first rough estimate of the noise is calculated and peak +detection is performed only in regions in which multiple sequential +signals are higher than this first estimated baseline/noise level.} + +\item{roiScales}{Optional numeric vector with length equal to \code{roiList} +defining the scale for each region of interest in \code{roiList} that +should be used for the centWave-wavelets.} + +\item{extendLengthMSW}{Option to force centWave to use all scales when +running centWave rather than truncating with the EIC length. Uses the "open" +method to extend the EIC to a integer base-2 length prior to being passed to +\code{convolve} rather than the default "reflect" method. See +https://github.com/sneumann/xcms/issues/445 for more information.} + +\item{ppmMerging}{The maximum mass deviation allowed when grouping individual +IM scans into frames. Data points within \code{ppmMerging} ppm will be +summed up into a single value and the reported mz will be their weighted +average.} + +\item{binWidthIM}{The bin size used when calculating the mobilograms to resolve +the peaks into the ion-mobility dimension. Lower values will give better resolution +if the data allows it, but can also generate spurious peaks.} +} +\value{ +The \code{CentWaveParamIM} function returns a + \code{CentWaveParamIM} class instance with all of the settings + specified for the centWave-based peak detection in chromatographic + + ion mobility data. +} +\description{ +Performs an extension of CentWave peak-picking on LC-IM-MS MS1 + data: first it joins all mobility scans into frames and performs .centWave_orig on + the summarized LC-MS-like data; then, from each peak, it calculates its mobilogram and + performs a second peak-picking on the IM dimension, resolving the peaks. + +The \code{CentWaveParamIM} class allows to specify all + settings for + Instances should be created with the \code{CentWaveParamIM} + constructor. See also the documentation of the + \code{\link{CentWaveParam}} for all methods and arguments this class + inherits. +} +\details{ +See \code{\link{centWave}} for details on the centWave method. +} +\section{Slots}{ + +\describe{ +\item{\code{ppm,peakwidth,snthresh,prefilter,mzCenterFun,integrate,mzdiff,fitgauss,noise,verboseColumns,roiList,firstBaselineCheck,roiScales,ppmMerging,binWidthIM}}{See corresponding parameter above.} +}} + +\seealso{ +The \code{\link{do_findChromPeaks_IM_centWave}} core + API function and \code{\link{CentWaveParam}} for the class the + \code{CentWaveParamIM} extends. + +Other peak detection methods: +\code{\link{findChromPeaks-centWaveWithPredIsoROIs}}, +\code{\link{findChromPeaks-centWave}}, +\code{\link{findChromPeaks-massifquant}}, +\code{\link{findChromPeaks-matchedFilter}}, +\code{\link{findChromPeaks}()}, +\code{\link{findPeaks-MSW}} +} +\author{ +Roger Gine, Johannes Rainer +} +\concept{peak detection methods} diff --git a/man/findChromPeaks-centWaveWithPredIsoROIs.Rd b/man/findChromPeaks-centWaveWithPredIsoROIs.Rd index 73cef5e4c..178c16ac2 100644 --- a/man/findChromPeaks-centWaveWithPredIsoROIs.Rd +++ b/man/findChromPeaks-centWaveWithPredIsoROIs.Rd @@ -286,6 +286,7 @@ The \code{\link{do_findChromPeaks_centWaveWithPredIsoROIs}} core the peak detection. Other peak detection methods: +\code{\link{findChromPeaks-centWaveIonMobility}}, \code{\link{findChromPeaks-centWave}}, \code{\link{findChromPeaks-massifquant}}, \code{\link{findChromPeaks-matchedFilter}}, diff --git a/man/findChromPeaks-massifquant.Rd b/man/findChromPeaks-massifquant.Rd index 38a89e8c3..e304b4d22 100644 --- a/man/findChromPeaks-massifquant.Rd +++ b/man/findChromPeaks-massifquant.Rd @@ -411,6 +411,7 @@ The \code{\link{do_findChromPeaks_massifquant}} core API function the peak detection. Other peak detection methods: +\code{\link{findChromPeaks-centWaveIonMobility}}, \code{\link{findChromPeaks-centWaveWithPredIsoROIs}}, \code{\link{findChromPeaks-centWave}}, \code{\link{findChromPeaks-matchedFilter}}, diff --git a/man/findChromPeaks-matchedFilter.Rd b/man/findChromPeaks-matchedFilter.Rd index b622757a5..708e884c8 100644 --- a/man/findChromPeaks-matchedFilter.Rd +++ b/man/findChromPeaks-matchedFilter.Rd @@ -329,6 +329,7 @@ peak detection in purely chromatographic data. the chromatographic peak detection. Other peak detection methods: +\code{\link{findChromPeaks-centWaveIonMobility}}, \code{\link{findChromPeaks-centWaveWithPredIsoROIs}}, \code{\link{findChromPeaks-centWave}}, \code{\link{findChromPeaks-massifquant}}, diff --git a/man/findChromPeaks.Rd b/man/findChromPeaks.Rd index d5efb715c..2dd06480e 100644 --- a/man/findChromPeaks.Rd +++ b/man/findChromPeaks.Rd @@ -91,6 +91,7 @@ chromatographic peaks. \code{\link[=manualChromPeaks]{manualChromPeaks()}} to manually add/define chromatographic peaks. Other peak detection methods: +\code{\link{findChromPeaks-centWaveIonMobility}}, \code{\link{findChromPeaks-centWaveWithPredIsoROIs}}, \code{\link{findChromPeaks-centWave}}, \code{\link{findChromPeaks-massifquant}}, diff --git a/man/findPeaks-MSW.Rd b/man/findPeaks-MSW.Rd index 337afc042..111d7b096 100644 --- a/man/findPeaks-MSW.Rd +++ b/man/findPeaks-MSW.Rd @@ -300,6 +300,7 @@ The \code{\link{do_findPeaks_MSW}} core API function the peak detection. Other peak detection methods: +\code{\link{findChromPeaks-centWaveIonMobility}}, \code{\link{findChromPeaks-centWaveWithPredIsoROIs}}, \code{\link{findChromPeaks-centWave}}, \code{\link{findChromPeaks-massifquant}}, From 92eeec18cd21d7a48c3d288e661f34e3a27ca489 Mon Sep 17 00:00:00 2001 From: RogerGinBer <45530914+RogerGinBer@users.noreply.github.com> Date: Wed, 7 Dec 2022 11:28:16 +0100 Subject: [PATCH 2/8] Change CentWaveParamIM to IMCentWaveParam For better readability and indexing --- NAMESPACE | 2 +- R/DataClasses.R | 8 ++++---- R/MsExperiment-functions.R | 2 +- R/functions-Params.R | 8 ++++---- man/findChromPeaks-centWaveIonMobility.Rd | 16 ++++++++-------- 5 files changed, 18 insertions(+), 18 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 2dbb0db9a..800ab341e 100755 --- a/NAMESPACE +++ b/NAMESPACE @@ -447,7 +447,7 @@ export("CentWaveParam", "MassifquantParam", "MSWParam", "CentWavePredIsoParam", - "CentWaveParamIM", + "IMCentWaveParam", "PeakDensityParam", "MzClustParam", "NearestPeaksParam", diff --git a/R/DataClasses.R b/R/DataClasses.R index 941e3616c..c812ac3e1 100644 --- a/R/DataClasses.R +++ b/R/DataClasses.R @@ -1330,14 +1330,14 @@ setClass("IMParam", contains = "Param") #' #' @seealso The \code{\link{do_findChromPeaks_IM_centWave}} core #' API function and \code{\link{CentWaveParam}} for the class the -#' \code{CentWaveParamIM} extends. +#' \code{IMCentWaveParam} extends. #' #' @name findChromPeaks-centWaveIonMobility NULL -#' @description The \code{CentWaveParamIM} class allows to specify all +#' @description The \code{IMCentWaveParam} class allows to specify all #' settings for -#' Instances should be created with the \code{CentWaveParamIM} +#' Instances should be created with the \code{IMCentWaveParam} #' constructor. See also the documentation of the #' \code{\link{CentWaveParam}} for all methods and arguments this class #' inherits. @@ -1346,7 +1346,7 @@ NULL #' See corresponding parameter above. #' #' @rdname findChromPeaks-centWaveIonMobility -setClass("CentWaveParamIM", +setClass("IMCentWaveParam", contains = c("IMParam", "CentWaveParam"), slots = c( ppmMerging = "numeric", diff --git a/R/MsExperiment-functions.R b/R/MsExperiment-functions.R index 0471dbba0..afa5d1393 100644 --- a/R/MsExperiment-functions.R +++ b/R/MsExperiment-functions.R @@ -4,7 +4,7 @@ MassifquantParam = "do_findChromPeaks_massifquant", MSWParam = "do_findPeaks_MSW", CentWavePredIsoParam = "do_findChromPeaks_centWaveWithPredIsoROIs", - CentWaveParamIM = "do_findChromPeaks_IM_centWave") + IMCentWaveParam = "do_findChromPeaks_IM_centWave") fun <- p2f[class(x)[1L]] if (is.na(fun)) stop("No peak detection function for parameter class ", class(x)[1L]) diff --git a/R/functions-Params.R b/R/functions-Params.R index e8167a309..e429eba02 100644 --- a/R/functions-Params.R +++ b/R/functions-Params.R @@ -227,21 +227,21 @@ CentWavePredIsoParam <- function(ppm = 25, peakwidth = c(20, 50), snthresh = 10, } -#' @return The \code{CentWaveParamIM} function returns a -#' \code{CentWaveParamIM} class instance with all of the settings +#' @return The \code{IMCentWaveParam} function returns a +#' \code{IMCentWaveParam} class instance with all of the settings #' specified for the centWave-based peak detection in chromatographic + #' ion mobility data. #' #' @rdname findChromPeaks-centWaveIonMobility #' -CentWaveParamIM <- function(ppm = 25, peakwidth = c(20, 50), snthresh = 10, +IMCentWaveParam <- function(ppm = 25, peakwidth = c(20, 50), snthresh = 10, prefilter = c(3, 100), mzCenterFun = "wMean", integrate = 1L, mzdiff = -0.001, fitgauss = FALSE, noise = 0, verboseColumns = FALSE, roiList = list(), firstBaselineCheck = TRUE, roiScales = numeric(), extendLengthMSW = FALSE, ppmMerging = 10, binWidthIM = 0.02) { - return(new("CentWaveParamIM", ppm = ppm, peakwidth = peakwidth, + return(new("IMCentWaveParam", ppm = ppm, peakwidth = peakwidth, snthresh = snthresh, prefilter = prefilter, mzCenterFun = mzCenterFun, integrate = as.integer(integrate), mzdiff = mzdiff, fitgauss = fitgauss, noise = noise, diff --git a/man/findChromPeaks-centWaveIonMobility.Rd b/man/findChromPeaks-centWaveIonMobility.Rd index 5eea3ace1..668c7fe9c 100644 --- a/man/findChromPeaks-centWaveIonMobility.Rd +++ b/man/findChromPeaks-centWaveIonMobility.Rd @@ -4,11 +4,11 @@ \name{findChromPeaks-centWaveIonMobility} \alias{findChromPeaks-centWaveIonMobility} \alias{centWaveIonMobility} -\alias{CentWaveParamIM-class} -\alias{CentWaveParamIM} +\alias{IMCentWaveParam-class} +\alias{IMCentWaveParam} \title{Centwave-based ion-mobility peak picking} \usage{ -CentWaveParamIM( +IMCentWaveParam( ppm = 25, peakwidth = c(20, 50), snthresh = 10, @@ -111,8 +111,8 @@ the peaks into the ion-mobility dimension. Lower values will give better resolut if the data allows it, but can also generate spurious peaks.} } \value{ -The \code{CentWaveParamIM} function returns a - \code{CentWaveParamIM} class instance with all of the settings +The \code{IMCentWaveParam} function returns a + \code{IMCentWaveParam} class instance with all of the settings specified for the centWave-based peak detection in chromatographic + ion mobility data. } @@ -122,9 +122,9 @@ Performs an extension of CentWave peak-picking on LC-IM-MS MS1 the summarized LC-MS-like data; then, from each peak, it calculates its mobilogram and performs a second peak-picking on the IM dimension, resolving the peaks. -The \code{CentWaveParamIM} class allows to specify all +The \code{IMCentWaveParam} class allows to specify all settings for - Instances should be created with the \code{CentWaveParamIM} + Instances should be created with the \code{IMCentWaveParam} constructor. See also the documentation of the \code{\link{CentWaveParam}} for all methods and arguments this class inherits. @@ -141,7 +141,7 @@ See \code{\link{centWave}} for details on the centWave method. \seealso{ The \code{\link{do_findChromPeaks_IM_centWave}} core API function and \code{\link{CentWaveParam}} for the class the - \code{CentWaveParamIM} extends. + \code{IMCentWaveParam} extends. Other peak detection methods: \code{\link{findChromPeaks-centWaveWithPredIsoROIs}}, From 857f68180b8baa92a1bbfbb9f3a30344af8bccc4 Mon Sep 17 00:00:00 2001 From: RogerGinBer <45530914+RogerGinBer@users.noreply.github.com> Date: Tue, 13 Dec 2022 20:01:26 +0100 Subject: [PATCH 3/8] Refactor .do_findChromPeaks_IM_centWave - Remove `message` calls from within the function - Split functionality in 3 parts: 1. collapse/merge IM scans, 2. call CentWave on collapsed data, 3. post-process peak matrix to resolve IM peaks. - Add termination `CWT` to `.do_resolve_IM_peaks` and `split_mobilogram` to differenciate them if more methods are added --- R/do_findChromPeaks-functions.R | 60 ++++++++++++++++----------------- 1 file changed, 30 insertions(+), 30 deletions(-) diff --git a/R/do_findChromPeaks-functions.R b/R/do_findChromPeaks-functions.R index 0257460bb..e9a0291d2 100644 --- a/R/do_findChromPeaks-functions.R +++ b/R/do_findChromPeaks-functions.R @@ -2275,52 +2275,53 @@ do_findChromPeaks_IM_centWave <- function(spec, ppmMerging = 10, binWidthIM = 0.01 ){ - ## Extract frame information - pdata <- peaksData(spec, columns = c("mz", "intensity")) - rt <- rtime(spec) - im <- spec$inv_ion_mobility - - - ## Merging frames into scans and Summarize across IM dimension - message("Collapsing data over IM dimension... ", appendLF = F) + ## Merging all scans from the same frame to summarize across IM dimension scans_summarized <- - combineSpectra( + Spectra::combineSpectra( spec, f = as.factor(spec$frameId), intensityFun = base::sum, weighted = TRUE, ppm = ppmMerging ) - message("OK") + Spectra::centroided(scans_summarized) <- TRUE + + ## 1D Peak-picking on summarized data + peaks <- .mse_find_chrom_peaks_sample(scans_summarized, + msLevel = 1L, + param = CentWaveParam(ppm = ppm, peakwidth = peakwidth, + snthresh = snthresh, prefilter = prefilter, + mzCenterFun = mzCenterFun, integrate = integrate, + mzdiff = mzdiff, fitgauss = fitgauss, noise = noise, + verboseColumns = verboseColumns, roiList = roiList, + firstBaselineCheck = firstBaselineCheck, + roiScales = roiScales, + extendLengthMSW = extendLengthMSW)) + + ## 1D Peak-picking, for each individual peak, to resolve across the IM dimension + .do_resolve_IM_peaks_CWT(spec, peaks, binWidthIM) + +} + + +.do_resolve_IM_peaks_CWT <- function(spec, peaks, binWidthIM){ + ## Extract frame information + pdata <- peaksData(spec, columns = c("mz", "intensity")) + rt <- rtime(spec) + im <- spec$inv_ion_mobility - ## Peak-picking on summarized data - mzs <- mz(scans_summarized) - valsPerSpect <- lengths(mzs, FALSE) - mz <- unlist(mzs, use.names = FALSE) - int <- unlist(intensity(scans_summarized), use.names = FALSE) - scantime <- sort(unique(rt)) - peaks <- .centWave_orig(mz = mz, int = int, scantime = scantime, - valsPerSpect = valsPerSpect, ppm = ppm, peakwidth = peakwidth, - snthresh = snthresh, prefilter = prefilter, - mzCenterFun = mzCenterFun, integrate = integrate, - mzdiff = mzdiff, fitgauss = fitgauss, noise = noise, - verboseColumns = verboseColumns, roiList = roiList, - firstBaselineCheck = firstBaselineCheck, - roiScales = roiScales, sleep = sleep, - extendLengthMSW = extendLengthMSW) ## Resolving peaks across IM dimension - message("Resolving peaks over ion-mobility dimension... ", appendLF = F) resolved_peaks <- vector("list", nrow(peaks)) for (i in seq_len(nrow(peaks))) { current_peak <- peaks[i,] mobilogram <- .extract_mobilogram(pdata, current_peak, rt, im, binWidthIM) if (length(mobilogram) == 0) { - warning(i, " mobilogram is empty") + # warning(i, " mobilogram is empty") next } - bounds <- .split_mobilogram(mobilogram) + bounds <- .split_mobilogram_CWT(mobilogram) new_peaks <- data.frame( mz = current_peak["mz"], mzmin = current_peak["mzmin"], @@ -2336,7 +2337,6 @@ do_findChromPeaks_IM_centWave <- function(spec, resolved_peaks[[i]] <- new_peaks } resolved_peaks <- do.call(rbind, resolved_peaks) - message("OK") ## Refine and calculate peak parameters vals <- vector("list", nrow(resolved_peaks)) @@ -2397,7 +2397,7 @@ do_findChromPeaks_IM_centWave <- function(spec, #' @importFrom MassSpecWavelet peakDetectionCWT -.split_mobilogram <- function(mob){ +.split_mobilogram_CWT <- function(mob){ if(length(mob$x) == 0){return()} vec <- mob$x #Add some padding, which will be removed after From 216210f6adb5ef2a68e6bb87e25535146e10919b Mon Sep 17 00:00:00 2001 From: RogerGinBer <45530914+RogerGinBer@users.noreply.github.com> Date: Wed, 21 Dec 2022 22:03:58 +0100 Subject: [PATCH 4/8] feat: Add IM peak-picking dispatch point by chunks Add the dispatch point to `.mse_find_chrom_peaks_chunk` Until now the IM peakpicking was only available through .mse_find_chrom_peaks_sample --- R/MsExperiment-functions.R | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/R/MsExperiment-functions.R b/R/MsExperiment-functions.R index afa5d1393..f5ef5d534 100644 --- a/R/MsExperiment-functions.R +++ b/R/MsExperiment-functions.R @@ -225,6 +225,19 @@ if (lx) f <- factor(x$.SAMPLE_IDX, levels = sidx) else f <- factor(integer(), levels = sidx) + + ## Ion mobility peak-picking dispatch point + if (inherits(param, "IMParam")){ + if (!any(c("inv_ion_mobility") %in% Spectra::spectraVariables(x))) # Add any other column name needed + stop("Spectra object does not seem to have ion-mobility data") + return ( + bplapply(split(x, f), function(spec){ + do.call(.param_to_fun(param), + args = append(list(spec), as(param, "list"))) #Append to avoid concatenating spectra + }, BPPARAM = BPPARAM) + ) + } + ## Check for random number of spectra if they are centroided. NOT all. if (inherits(param, "CentWaveParam")) { cntr <- all(centroided(x[sort(sample(seq_along(x), min(c(100, lx))))])) From 56cc96cb87d5f0630c1e69d95a90b03b4845b486 Mon Sep 17 00:00:00 2001 From: RogerGinBer <45530914+RogerGinBer@users.noreply.github.com> Date: Tue, 3 Jan 2023 23:20:53 +0100 Subject: [PATCH 5/8] feat: extend peak Mz range to compensate for combineSpectra Since we use `combineSpectra` to merge mobility scans into a frame-level scan, the reported mz value *shrinks* to the weigthed average, causing problems downstream when trying to calculate the mobilograms. Therefore, we are expanding by `ppmMerging` the peak mz range after CentWave peak detection Also, add checks for number of rows of output --- R/do_findChromPeaks-functions.R | 18 +++++++++++++----- man/do_findChromPeaks_IM_centWave.Rd | 2 +- 2 files changed, 14 insertions(+), 6 deletions(-) diff --git a/R/do_findChromPeaks-functions.R b/R/do_findChromPeaks-functions.R index e9a0291d2..dd3dc2783 100644 --- a/R/do_findChromPeaks-functions.R +++ b/R/do_findChromPeaks-functions.R @@ -2269,7 +2269,7 @@ do_findChromPeaks_IM_centWave <- function(spec, verboseColumns = FALSE, roiList = list(), firstBaselineCheck = TRUE, - roiScales = NULL, + roiScales = numeric(), sleep = 0, extendLengthMSW = FALSE, ppmMerging = 10, @@ -2298,6 +2298,11 @@ do_findChromPeaks_IM_centWave <- function(spec, firstBaselineCheck = firstBaselineCheck, roiScales = roiScales, extendLengthMSW = extendLengthMSW)) + if (!nrow(peaks)) return() + + #Correcting for the fact that combineSpectra combined close mz values + peaks[,"mzmin"] <- peaks[,"mzmin"] * (1 - ppmMerging * 1e-6) + peaks[,"mzmax"] <- peaks[,"mzmax"] * (1 + ppmMerging * 1e-6) ## 1D Peak-picking, for each individual peak, to resolve across the IM dimension .do_resolve_IM_peaks_CWT(spec, peaks, binWidthIM) @@ -2321,6 +2326,7 @@ do_findChromPeaks_IM_centWave <- function(spec, # warning(i, " mobilogram is empty") next } + bounds <- .split_mobilogram_CWT(mobilogram) new_peaks <- data.frame( mz = current_peak["mz"], @@ -2336,7 +2342,9 @@ do_findChromPeaks_IM_centWave <- function(spec, ) resolved_peaks[[i]] <- new_peaks } + resolved_peaks <- do.call(rbind, resolved_peaks) + if(is.null(resolved_peaks) || !nrow(resolved_peaks)) return() ## Refine and calculate peak parameters vals <- vector("list", nrow(resolved_peaks)) @@ -2381,11 +2389,11 @@ do_findChromPeaks_IM_centWave <- function(spec, as.matrix(resolved_peaks) } -#' @importFrom MsCoreUtils bin +#' @importFrom MsCoreUtils between bin .extract_mobilogram <- function(pdata, peak, rt, im, binWidthIM = 0.01){ - rtr <- c(peak["rtmin"], peak["rtmax"]) - mzr <- c(peak["mzmin"], peak["mzmax"]) - keep <- dplyr::between(rt, rtr[1], rtr[2]) + rtr <- c(peak[["rtmin"]], peak[["rtmax"]]) + mzr <- c(peak[["mzmin"]], peak[["mzmax"]]) + keep <- MsCoreUtils::between(rt, rtr) if (length(keep) == 0) return() ims <- im[keep] ints <- vapply(pdata[keep], xcms:::.aggregate_intensities, diff --git a/man/do_findChromPeaks_IM_centWave.Rd b/man/do_findChromPeaks_IM_centWave.Rd index d935de6c2..3df9bd655 100644 --- a/man/do_findChromPeaks_IM_centWave.Rd +++ b/man/do_findChromPeaks_IM_centWave.Rd @@ -18,7 +18,7 @@ do_findChromPeaks_IM_centWave( verboseColumns = FALSE, roiList = list(), firstBaselineCheck = TRUE, - roiScales = NULL, + roiScales = numeric(), sleep = 0, extendLengthMSW = FALSE, ppmMerging = 10, From ce8daa4069544204509d71c5017a4389e1112cf5 Mon Sep 17 00:00:00 2001 From: RogerGinBer <45530914+RogerGinBer@users.noreply.github.com> Date: Tue, 3 Jan 2023 23:22:49 +0100 Subject: [PATCH 6/8] feat: add getters and setters for IMCentwaveParam class --- R/AllGenerics.R | 4 +++ R/DataClasses.R | 20 ++++++++++++-- R/methods-Params.R | 32 +++++++++++++++++++++++ man/findChromPeaks-centWaveIonMobility.Rd | 31 +++++++++++++++++++++- 4 files changed, 84 insertions(+), 3 deletions(-) diff --git a/R/AllGenerics.R b/R/AllGenerics.R index 0103bc603..18bca35cb 100644 --- a/R/AllGenerics.R +++ b/R/AllGenerics.R @@ -255,6 +255,8 @@ setGeneric("baseValue", function(object, ...) standardGeneric("baseValue")) setGeneric("baseValue<-", function(object, value) standardGeneric("baseValue<-")) setGeneric("binSize", function(object, ...) standardGeneric("binSize")) setGeneric("binSize<-", function(object, value) standardGeneric("binSize<-")) +setGeneric("binWidthIM", function(object, ...) standardGeneric("binWidthIM")) +setGeneric("binWidthIM<-", function(object, value) standardGeneric("binWidthIM<-")) setGeneric("bw", function(object) standardGeneric("bw")) setGeneric("bw<-", function(object, value) standardGeneric("bw<-")) @@ -736,6 +738,8 @@ setGeneric("plotTIC", function(object, ...) standardGeneric("plotTIC")) setGeneric("plotTree", function(object, ...) standardGeneric("plotTree")) setGeneric("ppm", function(object, ...) standardGeneric("ppm")) setGeneric("ppm<-", function(object, value) standardGeneric("ppm<-")) +setGeneric("ppmMerging", function(object, ...) standardGeneric("ppmMerging")) +setGeneric("ppmMerging<-", function(object, value) standardGeneric("ppmMerging<-")) setGeneric("prefilter", function(object, ...) standardGeneric("prefilter")) setGeneric("prefilter<-", function(object, value) standardGeneric("prefilter<-")) setGeneric("present", function(object, class, minfrac) standardGeneric("present")) diff --git a/R/DataClasses.R b/R/DataClasses.R index c812ac3e1..d13d35ab0 100644 --- a/R/DataClasses.R +++ b/R/DataClasses.R @@ -1300,7 +1300,7 @@ setClass("CentWavePredIsoParam", #### Ion mobility peak-picking classes #### -setClass("IMParam", contains = "Param") +setClass("IMParam", contains = "VIRTUAL") #' @title Centwave-based ion-mobility peak picking #' @@ -1356,7 +1356,23 @@ setClass("IMCentWaveParam", prototype = prototype( ppmMerging = 10, binWidthIM = 0.02 - )) + ), + validity = function(object){ + msg <- character() + if (length(object@ppmMerging) != 1 || + object@ppmMerging < 0) { + msg <- c(msg, + "'ppmMerging' should be a positive numeric of length 1") + } + if (length(object@binWidthIM) != 1 || + object@binWidthIM < 0) { + msg <- c(msg, + "'binWidthIM' should be a positive numeric of length 1") + } + if (length(msg)) + msg + else TRUE + }) setClass("PeakDensityParam", slots = c(sampleGroups = "ANY", diff --git a/R/methods-Params.R b/R/methods-Params.R index aa1b1e454..22f7416fd 100644 --- a/R/methods-Params.R +++ b/R/methods-Params.R @@ -915,6 +915,38 @@ setReplaceMethod("polarity", "CentWavePredIsoParam", function(object, value) { return(object) }) +## IMCentWaveParam + +#' @description \code{ppmMerging},\code{ppmMerging<-}: getter and +#' setter for the \code{ppmMerging} slot of the object. +#' +#' @rdname findChromPeaks-centWaveIonMobility +setMethod("ppmMerging", "IMCentWaveParam", function(object){ + return(object@ppmMerging)}) +#' @aliases ppmMerging<- +#' +#' @rdname findChromPeaks-centWaveIonMobility +setReplaceMethod("ppmMerging", "IMCentWaveParam", function(object, value) { + object@ppmMerging <- value + if (validObject(object)) + return(object) +}) + +#' @description \code{binWidthIM},\code{binWidthIM<-}: getter and +#' setter for the \code{binWidthIM} slot of the object. +#' +#' @rdname findChromPeaks-centWaveIonMobility +setMethod("binWidthIM", "IMCentWaveParam", function(object){ + return(object@binWidthIM)}) +#' @aliases binWidthIM<- +#' +#' @rdname findChromPeaks-centWaveIonMobility +setReplaceMethod("binWidthIM", "IMCentWaveParam", function(object, value) { + object@binWidthIM <- value + if (validObject(object)) + return(object) +}) + ############################################################ ## PeakDensityParam diff --git a/man/findChromPeaks-centWaveIonMobility.Rd b/man/findChromPeaks-centWaveIonMobility.Rd index 668c7fe9c..01b28fa2f 100644 --- a/man/findChromPeaks-centWaveIonMobility.Rd +++ b/man/findChromPeaks-centWaveIonMobility.Rd @@ -1,11 +1,18 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/DataClasses.R, R/functions-Params.R +% Please edit documentation in R/DataClasses.R, R/functions-Params.R, +% R/methods-Params.R \docType{class} \name{findChromPeaks-centWaveIonMobility} \alias{findChromPeaks-centWaveIonMobility} \alias{centWaveIonMobility} \alias{IMCentWaveParam-class} \alias{IMCentWaveParam} +\alias{ppmMerging,IMCentWaveParam-method} +\alias{ppmMerging<-,IMCentWaveParam-method} +\alias{ppmMerging<-} +\alias{binWidthIM,IMCentWaveParam-method} +\alias{binWidthIM<-,IMCentWaveParam-method} +\alias{binWidthIM<-} \title{Centwave-based ion-mobility peak picking} \usage{ IMCentWaveParam( @@ -26,6 +33,14 @@ IMCentWaveParam( ppmMerging = 10, binWidthIM = 0.02 ) + +\S4method{ppmMerging}{IMCentWaveParam}(object) + +\S4method{ppmMerging}{IMCentWaveParam}(object) <- value + +\S4method{binWidthIM}{IMCentWaveParam}(object) + +\S4method{binWidthIM}{IMCentWaveParam}(object) <- value } \arguments{ \item{ppm}{\code{numeric(1)} defining the maximal tolerated m/z deviation in @@ -109,6 +124,14 @@ average.} \item{binWidthIM}{The bin size used when calculating the mobilograms to resolve the peaks into the ion-mobility dimension. Lower values will give better resolution if the data allows it, but can also generate spurious peaks.} + +\item{object}{For \code{findChromPeaks}: an + \code{\link{OnDiskMSnExp}} object containing the MS- and all + other experiment-relevant data. + + For all other methods: a parameter object.} + +\item{value}{The value for the slot.} } \value{ The \code{IMCentWaveParam} function returns a @@ -128,6 +151,12 @@ The \code{IMCentWaveParam} class allows to specify all constructor. See also the documentation of the \code{\link{CentWaveParam}} for all methods and arguments this class inherits. + +\code{ppmMerging},\code{ppmMerging<-}: getter and + setter for the \code{ppmMerging} slot of the object. + +\code{binWidthIM},\code{binWidthIM<-}: getter and + setter for the \code{binWidthIM} slot of the object. } \details{ See \code{\link{centWave}} for details on the centWave method. From f8f48262bf3cb1fe7c348feac69ef334dbd07f82 Mon Sep 17 00:00:00 2001 From: RogerGinBer <45530914+RogerGinBer@users.noreply.github.com> Date: Tue, 3 Jan 2023 23:25:34 +0100 Subject: [PATCH 7/8] test: add unit tests and sample IM data Add a compact IM data file for running tests Add tests for `do_findChromPeaks_IM_centWave`, internals `.extract_mobilogram` and `.split_mobilogram_CWT`, as well as for the parameter class `IMCentWaveParam` --- data/IM_test_data.rds | Bin 0 -> 18846 bytes tests/testthat.R | 14 +++ tests/testthat/test_Param_classes.R | 115 ++++++++++++++++++ .../test_do_findChromPeaks-functions.R | 42 +++++++ 4 files changed, 171 insertions(+) create mode 100644 data/IM_test_data.rds diff --git a/data/IM_test_data.rds b/data/IM_test_data.rds new file mode 100644 index 0000000000000000000000000000000000000000..c1542746ebb80d98ea619fd82207c95860540400 GIT binary patch literal 18846 zcmbq)1yClzk|ypxz%aPG4DRp`?(XjVI1KLY?(XjH?(Xg|xVziUd;4}{Z_hU3Ix@Pm zzpl!v&d$iJs;njmhXMKL0e+zkt%W}P#2-c&5<(f4=R!%$CUjoX+MA|UYY_oQOfAt4 zUJ*hOrYCy7QseBeX7JrT6w*DCSf$5LZ}8E(nYn&hP5N_{cXImS&Ncf8-(}&3gRmf{dB2Oc`^u%B%TTF8J+KET#5v(Pdf6QY5-&aM8CHc3m==E z^WDsciMMWd?q@GJT#I;YjhuJ4ZJ-VOwyo~773^~fnw*crAq9VI^V>K&gZq%qRnu)o4ZIu+x=DRN6RZ8 z^4OQXj@d_{>P^-c_r}F75&RJ#*4}O!P8PFch4wZnGkcmgnd0}=2_^{_?l;bZPT2~% z=YX*HuW{(qwQ0C*n%`Y5d7qIh1vLc8(LZp#{@hg&WJ->BWTSnO{G(%kd$)}vzT!xg zkBk<|_`bb@IRaZkEJsFYnllEV2Ppe!`IdnxgO?#GftQfS;K^IBr>ia5 zJDn*Z%-41e?ComFuKYPgPQ=yn_-$j+aVS#QCyYy4y8_ij=-Ge7u2K}t>*M{p*r^lzOUD70hWcn}KHMT9Wn{d{@l zeuS>cbR|l0Jz8i92+D5fM+_9HZNIA2d06m(z%K?Nj3~N<3-6Nu`us%{7Rw*H{HM`{ zPAB}ob$RG+TT%~(dqRdMz(P4@_;1$wH~C+u8pO5e)5$2}q_IM1C+=jSDaAr)|G0UL z5-0A;7W@Y+UT~E@U~eq^Z@9I%z!nKBo9Kh?-2w6U))8U)-?8CZ`~Sw?GT4W5M4|mx z)O&#JPZ#VzM3KgJ8ODX&dH&a$pV2KsxT2x|F6Un&U-}ATVS%6{37c=VtipBO1;e$=5-KRh9}+jCovH!hlrJn@`TOllq26RT#0 zHkM^`hw4U3B_Fpo^5g@Gx_)?ls=hiqp67-qasf& zlwXu>GgN(IwGk+JY4gT<@tAQP&c6#L8m~^j-<@GPDvXcUduaWxb&=O{SkJ#W-pbxe z-UxnGbQ-iKu&!P&jpn{U{=n z5nbldj(F%kPopXpaaboBr+;v|UceWrvkWzL(yK>4{rs#RB%=4mQ|0|*JT6>Ep>2gd zgr7RSCg~JAkmsAde0=uhjS=j+E>b9EL+Z|CR2uEyC!WMRWSlPP?5*x(BbsXY3|P{+ zKi}&Pf==F+!?o<-R`$XvfBb$*u|Qwajq>)KeTNf<+)MVh*`&$G=e*!Fb2F0A*(z$h zN*9FQt#)Gi;c}b5d!P86j>lfqpJ;XX`bqld`N?)6RME}#47&dyxEa~r?v_6t=qamub2r>Mf{S zsL+1|-9mY*vMfI)+sND*HFt4NR@Yy=<#dy7v(i1c9(ntd^L#cnxbHWU!*XY7eBQ3l zSR?z#%R^iI8YLz8pkb<~xlv|vtqqHh(2!z;5$8gD>nd6AvXQL^>*D;Op7BCIQt)|n zl}0;R{7Es=RMpR3XrNNv%{j}n=-m@$BS2JK2vyu_u?Jb$%_j7L!&KKG=MhQjU4*++ z@J#%HDBzFj?y1V6w6vE#pF(`Noq_qV`Jp|NE@Rqd%A}B0`jnkhs$JPdC5#a3ZJ`id zEF~H8&mcrP_aw?9>)Q0fwqHAVy^l^un@Ib4RU&p$9Y3gHU6*cxW$t|Sme&+!BTF~# z;aOf0z|I3EY(R0F)9shY^>nzVLoQJ(>c2=-EEF)fawE8#MuYF<)9(+U|Px*b? z><>!wDY4k_ckW$y3D0wjj?P6X>0Ar~7k5 zy4x>y{x?0t8%UmO)wi?a`|qxzZ~H-4G9el9Uk#eKtpy{(GA>81cWp9goa6360r#~S zts~vvwoD~QKRUhj@KxOVcQ^UDl0)9e|K)(Y{#wen|1cZaYRQj{%>^%&=PiOD_k#;y z>-j?Om#dEbp`LKv{jqnFs?)|Q6?ah|`<`+B)@DV@V_1CUzCg?FY?Kba~P%Gju@7uWUq;<6xCydy;feWm9hi84) zllo=JK^6ATf4uyuvE{?~W!HnyfA;zRbjiAh=Nx5n*~bG9>4IK_e65+PqJzjirfjsa zmc)Nmlf7h~IC#Z#E#7XblSOvR$Jk8Mp-%l<2A_8)brK;qJmb1#5R)!DB;z3e%-PX` z?su#=g@?<1L|p%0)#%(6ejWcECkz}F{7W@eqYEj!MQHFlfqtI;@BKUcY=Fd(;KN_~ z7w^F@&~sOmec*oYqk)ap**wIAw1hl_q%DF%C@%}d6VvR^$vdE*BqPwmP|rbcPO!-G z&q#UQ4N9I1@~A&42~)iO91-AF!N6ke@cVcuaL`hcl6us!u+;D%a#;+`+>l1jqzJG$TB!BIIq; zW1v7ymbAW7ZLkc`OiLEG=U{U(NYF z($q53(8R zoJBy{6BCHZ;h7J>`YsYc*vvU%g~X)U1Y!=)f<=i!2}6lR2|3qS?* z`-*@>f?2#nm@|S>Ll{9RVU3`bz@QwLm(Z8G zG{-8ajMKlfzxXdeD7bq@t+7z_jM>p(YU=!N{Iy9@cJz>~)&#j94Sk0|Bca1U(_qly z(V-!*R%M~mz*sPV*fS1-8QPpxfF8x8y>E}o`4PBp=ZgRThB$nrkaCcr zpF6?dr$D>TK)juOy19MmIRB4?%;Ce!`7a{i@R8>9`FsAq&Q-C{nUnkz~{~Pw0K=qiAUxa#2A^A_Vr=LRt*s+M(F_p<-FSgs?l*)lG zh6UsxQ_KW^JMM?Q--UR!1@q+a<;wsN$owzxuk3%pJ&XTQ6m;kx<&NvW9S^`d{DDZ} z20fI?`9{&R32W;Nwvh?q{7=@-391ZX9}a?VPYbbxOAaw74PXsu^x5z=2X_To1TBMC zLN4Kx1J4NpMEdafVuQ1Sf^FME5jZ4)5)m@wm;fey3HdmhA*CF3GA!{5#Eww>3zxOQtf3d%=s&brFXWD6wH! z-|)yZR}r_7LJd;IcXHSvryOJrs_Pq$#&PsRvwV&*?g)vd?$m(x6gh_JtHBG%L`Voz zEFQ{81sUnsrGV8XFmDEX+bJu=WtjHh4+J$t?ZL#JSPnslQ!EajobQtli4H;!s8rh|n1#sW|DNE-PS&mZGI28yFG?2gC%t1h08!4(`$zap zLZ4cI9xVtqAu3j2pN~Mjpz8N4eR;X?1!IBg&v>~Oc1+_#nTcQ&NJtX;St(aVlu(sh zckrvg_gMnrZw`J5L5gx!=Oz~d`Rd&S2QujKVXd-d*Km1Zt;uq{Kpi}lU$-6U?Uxzc zA%|r8@?fC`UcJXTNEuE2@qV-q!tM%}c3rySh&ZpV7oUU5v+qW zsg(-5&e#k7J|`ueZZqM^e|IDGzB!q56D|r08Lo&RLDgu034g?q&wxfvT=<60y8vw$ zYCS{{#&I$)L_94ig+ztpdsso57hoPrX$gbZee?SezOzI|SOLWW(_zqqvbGU~fPU5m zAy=p%!ya5$mp%+_ZMO)=Hm^?+CUD1=owZoJrBZ~FH}f@q)6voP=cL^W^C%J|seY^h z40VB?5dXI+ykQqOq(D|7T#ad}q9l!hK1DE|gv8kt7uqmn6)fyuLKpGWHDIhl?caZJ zuk7aN05e@9SI~crBf|%*Q3p2DnLiGI_ix(G`_Dl((J6LP{~nueIq~I#0;U~ye4}d~ zyb}35DKIPTkuDd37;2qAO%ws=a2($bA#9c&%fK^~!!oUQ!R@4b-G3{?QXu_@u@&mu zb1AWo0B*A0N*guP^MEqp`=wMv#NWwyEDNd4rlz+s#eLre0;kf21<;UG9_t0iGQpnW zsCHsH-X(T!@iw;w7e+I2bdMf?paraUFuh7Ktn zGMWfWiZS`D%A?yJ#~eLqoVByVNlTI+Cdo`pYwFe48wrNJZ`m|`%DnVE1t-lvd5iQooCpS{y~<|PmBoD$CX96`Ef)t95y#awrSzM1 z&~UDZ6JmV;0cT@;cT`LIPOVqr2>0jGq(ab;VkZ1%LuwH7p9}#ewF$8tHX*B$8d`p%%?wsCUti%!QmCKIa}DAxo(~DXEzggl zoZeY8HJRGHrUS75`?tfq^+7|H>+1g8k551h&1tXhkV;Qioa+aZw;9@}0wzr2 zkdrBS$ay!;T7!P3OY~=aqxZ~JeFx3B$pJiObxHI0=1O;F#Eh zQEjmg9h;BoRidLtWu!)D=k>QY_Z?2rez7V>!7&;9=<;WH@LqG*hbP`(K`i;%p%s7k z0+vQa2Ind%SA)7D6{)_Hm;p|*wbx@5+9RD>_VIN$_08!?P~HO(0Vvb~Iw8S*KY)@c zgM`s?(uZvl=&OSD#HDxh8|u+`EB&{mnGNX#!E)`?iTOMV=#?*P2Q{= zi%$r;B<4$NW^{Lthx6^jO`~(oXv#iV1LYz}Q!ZYLX{?Sc7m>$&E8oXh{;=gp>Ma)g zrZ}1&j;64ymQF}2$P5b525<|TO-5OD+@OCMi^Qew2tK%21V)wEDBcH=Tel|FGx z)+V6_H)z#)o;DVl1NeD_8A;qE<@y)4gWvKJiW}qVHi$El72^*X;vvI2aeV^?;{i0` zY@4B{bqv zxz!@#j1t9a-$lNhL>g$Cd57(`nrK63B&dUN_X4ds{o2wzh;(Zsj>uBiqOTBfi`9F6^Bcxz;4hy;jE9i7o_;v zD6&5G#NB1U)*gbyi8YA!;A35+TQph`+*ZPAELwFB3d_3oPSLeJw!Ej!&-E|Wi1r{f zQ)WajeRtaELENU1PX+M%xsL+x^{ok{jT3*Td@3w87Rh5O9LI=^rxg&)F!xl0kVZ2r zz#a*rHQP^n$UJpr%yP(?ZK-8<{c$Ev6zVDTu#?VQI+mrWn5848ZlS0jq zT^<*|n+vnArbSrp;JCyj zFWr9+O80|2i-rciB~ng^J{MiG>MS$>?a;O zg5P)b@}*&^EIRoGuP@=7&SvOXnufW(d4s**9n2u2l@>s5$Shzq60Knu*0;BZBb49{ z%PQ5+_4Di+k$&3XQGmOv-bgqDJzAJJ`nXOM5dtl6G2S+GY)^+aOydt`J>&2UOGY1J zyg9FUFa8;8eBom`xh+rU-*(~r zqPnAdN_K|3$_6bJXvW_AvQiqOSD#U>{axK<5?mZ~s3n5$CX2XpOG`D_h_78zVfFp4 zelsbq_ubXbPEScF%`+5AUA)f(uxzY*wKM}!iZ5Nkv;?1?KF7k=n4=ydT8nF8Difj5 zgy#iA#5mIlq}(dl#yJdvYA;jfp^I8&HUBVLqs0EEW2-og_v;h1n;$SkL|eJa1TM0< zpoT+vc3H=u-_4$Y(TpTu2Ta7=0WoPdkX^jn<^2(0XlbZ_xNno`JuS4cr zZ)~jfwA!4(2SIWE-ZPJ>jja7`YeGrd%93e)V)bV|ezD@mvTTtlDeNgv>*7Y=maN`m z6g^L460e^vO!1{lgvuqwRbc6RQ`X4`Yv+p2e1{a$x)At$tI+NhQHpz|-y=tqsQf4x zc$J&L#R$*KwSnc8@|~065rdvZli0>O!@|Y;!Cu6Pyru;+SzkKxx9PcW#Bw1@dYj)( zKvgd|(&lDE)d3lao5|v7cvV7qvAF{Xu+$Ei0TW>`8s!Ip4i+`NFr&KBU7GFkwB=a6 z(=z=ut)?o=5G-L?u3UVGkxxU?!Nc$y4Qo&G?R zdgON`Rlv@|E(N(d} zOPPTFy*|pQbLC=Ie9o5(cst+WU6M%nq1N@f_PO|->PEo>LT z^`ElS#|4kyWn_BogyYXu((j;DG?#`2eY>T-!fddl_; zAlKb&LZ!~1({Xrj%Jl9g-d}LmcGsQEo6{w+`-6r>iMS=T1X|?Gb2_2P@a&qgv@d-$ zd?cP8uPxWM^x9`@{kv1gqCB=$fEAo*qHZi`T3;!TkeN@B?K6HFX^U`~Xqd|oB_*Dw z5uxck&yq-k&~Rg{)(c|VkgFY%nOnz9dd1{|Jy~?ygU!`S(!;pj{{5LwX$tZm7Ur$s z!?GjiVNOdx?cb0x;{s0UGnG6_Zkn|NH)WcpnNw1Hm7-g?pi)*j}q|I%?LAWP#g zIOqtb0wfYPtW8V$hC{yf>^SIf^>piqXIhBuPM(9Y8F*`w;-lp2UbzceU@MPBQO&6{ zh&329NoBanc%=^_&xEh0y{RE#_p~}#I{WZ3)Qvjt z3a+I8x>azzOU;r4$Gt^nh{v>5Xu{%muX5?nNktW4pVfCR7LDh3q+wZ-+ye5>$#@j` z;L4Iks>=0QxX$z@*&?)O$BjW%lO=n3(ojQY-y<;nQf=mtpGbo}Nd@;sMD9E`#ht8L z>zibWjmmOj0-;HZ%PogzyFF4t@#Hmi`Ivtdeo~D7*S|TTrE|~Us15JAR)rY#{+*T+ zoVC}v+wTS=3`HoCDXm5h!EDnkip@<^qecot0sjs`9>(CIPB0?4{4B0S4Fi`jvLKLP zx_Ws+vlULT)*H%>@z@(kGyqdbmXVP96Bg@e-TU|wQF=nTt?A%vR$;mVQ9DzO4%l9eD%{T^Q5lVvhRB_ z@0HsRN4XSg!j<;}#IXXff%=d{1bHn$vg8t;E@Zzv@@*Y^+&l6BXA8lzyvg}m80L}- zTYRg&qeD4j3Cn8i+;QbX(Lrd2{+K+uizv@aB#)nvOTZ-*D@!k3?;Du@f#^4RJ zq*v3q))cs9Y}ifGLZ(*hpN2ytlsl0k+MIP+tSF0a0tc`!V){RR&V6ZEIN;lu3eAgw z-AfN=koO@C=wW`Vl#=Q1U@%>VyPIdVI5h4GZI(E8js=QlJ?$2+=z7M!-#uIsghzu2 za&>y;++nczIqX7k`uR`b=AeI@IN2E{aQ5vI5s+(MZ2Z{J7lzty$%|Sx>8EkRKa4~d zW_NzwyP$QbLB^yTfwpg}^vvvRI{p;ENhuY1AA)jd;A2Y}%s5KDRmxH;O`o!=HwK0$ zo)U^$_;)l!B&10G3~754O+B}DIA+7N!llv$HsO4a=G=6_q+~HKdkW+$eQI6Ng@oWW zXn)y^-m)FlA5Za=m=60^1huNQ=dqO5CHChKnT}J$Ae!)EIriZ`{xg1b{c&RWszHn6a~EsZFD%ZnKqoC5!rno^1!qgtp*;4A zlm+^AUQ2IUV5Fobg&o3w_^F#ttcN=d(Qk+vh{jX$80jQu3T zG^<0uThMJ4+M_aB@e5>+;_Fmzq^P$06Py*T6sXhK-!4%S6xPUgfe$GtbXK<5vCSwr z=)-uW7K=fyJ!Es~1x;09rR8A|Toc}V(p|uuF=Y=~WOzxkDDpAG=hwx47^_fKrExv4 zj1x(ZMpvU<;qYmDrL{xHg^Syj!HWJ;u29G?198#mwUPVPK1f-Pw=`GRs%GR`$5yP( z?Y}tU2*~i4T-Rl7cn4%C#WQRRtvw z7L|8^VjT$3+8MX=B>;+rT${Oxvpje7^1;L*{;*8+4rxMTZ$zLl?J%@qb3~xl9_G?7 z&T=G#JSN~9e&s?F+;3ZME>1|z;qx~w?dJXBv#JRW=HBC6D2Q9DQc%&D_nU|*JJjkJ z9_A|1gN}ZX!G+iIyCQJ@1$U!U*`96cqi$R+2+7<=BLywE8bk{g$>JLK~_K3=47K;3^v`e-;uWT3d5PX zr|N?zmRwfpPKfr?g(?E2uwpxVX$OkEmxCg#^K<0pS2#GK?<@|kuYf@XcfCqv+*Usc zVxl-b&v>3*Tj=;_*2~!l(fJl#^1K+YDf{&==9+bV)VMQCmc%3`JYeLic-sn8NeXH> zFjs408KrS_$cIy4jU5gXoR;^*lMTg^P0`BR6McRrKnS#bOBCQE4u;3hs;I+V5U*^@ z;8%=Yw_^@x?^NjBv*EZd3#Bmyo%a5IvbT?nbe`eWnS`n!rctDeOE5Mh0hc&PU{kxO zreB_$nyeuV=1Um`U`oHn@s*&X$y1SHW(l-5M*r@x^bR-zL&{)j+gbs~HlzIg=QCRX zq^+Fa8Ci0{M~3}lCW43B!_yWE{kBMZ7+PkLK}YB5axtfo6`jGPmIlmQzVn2)#+^M9 zZE-%8XzJXoDA$h{e=48H7nwFJ6mAG3bH$lHnD78{EyT zR-c;nnN=1h(3MzGG;=RHk4BXUFDDH9sbuF2m^GvA*#jgqNs)w)rOHp~gL zAb+|OJ>!0#C%$wi7b1DNqiG!5Y>nRMEVQ_yXMznYis-?~y?JGP`Ut#^kNYF`nOk<= zX@$5blEHU$WXw1}bm5U)_rqP9(<*TzEpE&fo8$y&0Lt0u$mP0>M^&U2V%C<|oxYBA z2W_S{IKFZ1_}9MhT;vl!rdqkCkdi=pSwY+=G;P-CQB2kk(wtJB{V0@xww5avxXOR!>^{C2;l71QkH zq4R>jlZdKiDhpQP^gDvfPf46T`X8P_=QxK&XW9t*jCVVi{hD911c~e3b}#YZMRw$3 zVZYGfUqhy<=))Z<@TIOmYH@ZW2tpT@M*PWgF$-V`2Jj{Tq`_yv zZk;tk<$#W7+jIl)Kx;y?ga@!ot`rIeTmObMv_@ucIv2r{Sy5*l<0PEZx9s_83}mPE zsfRGo#7f8P_j#=y*XkmCMnzu2KiqDq2g~zUcV2O4e7wq{F^BKU%Jl}Rb{Fz|`$@}* zG+U0i<+*(Z*7oC-(3hO$=F>!{Li21l(OkV{`8%b$lgt@pqAH+V?c1Wo@E$yp1Xn3!H~H zZ>@$hPH%mH_YTi5EfTrSv;FNyTow22@R4uY-d3X-xrXhZMP5~|>1uhkM=;l8)PrY* z2}kFM(Tl9k(G5wIYRPv!hT$eMvu}Zi-3S)X3s)syx=FE z(9O-$)gmhL$4l;_c$rh_ZmP6b;5q z@H~6L^-h;JIb->B_rB~nUT*D z#gcKvX4t-w2tT$*!Y&=6j2?A_BH=+29d-DQjF^|;i%`VG3x}X8!dX%K&3bZt!zy}f zS@GIJkC|eDXEa&D}1MtJ9#|+QEG(0s*U*&XSbAW zfMpd@^a0QN0jk-t^3NR62EWXbjR(1OkNt)us8`w7CN)D8mJ**)sKX2N_s9 zqh(6mn(Rl3^pC+_6^50f>*nYpWO$upO25SG8{0byak`2a{WW?S4OmHwV4CX}n?QzL z>`UILM!)0w8g;M-R;xWC1|b~D*R98DCDd=Fqq=U&SXrG+Nzc@{g;)ZM3tM->$CqBU z4>Vjo*|^FUi>oRg@vZyUB(MuNF*40FE909e;3` z3*X3gt0tB%;V$}g>Xc)qPUC%}bF^FEBxMFcA|Vmmp5J&zT}BZob9IbV%$n$aT}Z~M z<#S=I)!oIKe)a_knQ8$sB;qHbWGm^hg=k!@t(+Sl!T9thHuCKMuK!u)!x?L>yqCpz zVB^-z`$q-@>+aCcP)!FaSY62J=H1znUcoJ6C^oW%5M>V56Y8#n8cw*+4>TcL4@t+b zEfY$A)07GypX(yb>+Uk53 zV=Jg!bx`4GZz_85ODn1))(Cz5SMkZUl%p6ml;sIv;R(_Ck>$$TH1Y%d-Ku zve$KOB@NS6V)=GT^Cn?%H-JyXx)R?n zOA>U~Jj7b)5~h-QtZ#6aUI0j^w^=6Z9(YKC8?DK{ zCk#IO&WU&npIZ34I?<4Q!|$M-W)c*2c!4=tlenlNjXQVP#uvnpm6*Z19}v!UicWg> zbF_7c++042y*qX@j9g_jky#mkq%a`{)vnAKyC!RBAnk7+G}0liN4h*^VM4Ao<2)>P zwov!AMat6%WN&LKn@{xaBDh+;p;r~8?g8Z4E{nHMl_OQOGTP~PeP&!h5BtWuT|o}t zxT8NyW(-Z0EFI#t50V7BGDd%&q$ds(zU`GjDxOtPEb7;Jr0f3yW6dXwP`t|oJB-hg z9#N~}BG&LJ;}uO_VE8$pTn5S6SO&TO)^2WA9>mvC=X3YMqv@5y@b*+m->i8a5@mRQ z1k@FA>_XWjFi5Da*xFwb`&?z8b~P?r!=?+v*L`@Qhd=PgX|X65TP_8!LLi-q2_3ZL zX%UpT-5_w3g1cG4ok0G>zM(Jz=Tk<>FSlCWEOaXE{k<*2flby>*jLK=jFWz|MY4G> zMx7Z919n&kqsBGoXf$M(G$Q|o6k(pr4|m*WY4?dWU)8Z8H+RV);4T;Nd8l_)#3!q! znCwpR?FgwTbL>YZciC3gi)!|UyPmP?(Gt$P!bMLTlBSFL;CZ8H>0Ygz@n>MW1oN5* z*yTTV-r3PY$0WeKv;v=<`Voo*H5Knk6dX{7il()}Aw7vSk>@~9L zovr&sZ#9$h7{5l><_wN6q=Lxbe3AX>0<&^YPruh0j(Lr8etOOf^=-rCM6GC@8B3qq znwL9>8*oJGFA%AZh)1krRev}ebn3*xae!al^kvU;MS9s_a<0cMnKb}cS)p|2iH(BT zlg(Y$@ac`WIxklM3OzbukfD-B)veC7CxcV1^WVD=r~#~N+=)ul^h-ed#KuQEGp-r# zLZ%E2&d9`tH1^oYj`8cOm3eEJ0or>U88bl{SFqMANFF&q9lGzP*^pvq?)X;^SX2?3D~KHjm9;*f>%P3f?taPy8<{R{{MQjv z{!kS;D2gzh8R3k|Xl|U5y6*%3>Ztn}7Ftt?s{#b5F6~4xhRaD(kN(a~}d)f>gb&xV6O(#Sssi(`}AlHz?&6 zykgvE3%aJ;8%q;M$lTIZLp4Lue4#(}oT*bPLFjl*TXYJ;iPr|&jg3_}PDd~M`2_MU zV+oex=&{9v#0%j^tw^1!pCaD#EsQmz_h!%g2pDoV{w`Daokz#@8YYRsb5U>q@w#<3 zG397sUDDqP{ia@9ruimFKVnwXc-i{uCM#OnoH*0DOlp-z8)>%MUL}%g(S5sT3|&xD z_Ce_7Y`=5XhtAEk3W!%-(jJ$f8YtJ)lL>}lnVC_O4d;CI{KH!zfqDA=aZYH3Ay=K7 z2wxd(YKQ;WqgbJO@w~ds*zfflba3+zpZ#{A83bu2FQVJq0UttY6T#eJxC*=*&v9{e znkx5K$Tz6F>0rN3iT=&aCZ_*8gevaGiOzD}#M~0r|J8`~F(GWF@7oQ20Z&KO0 zvDtv_DV+n0XfWXdjtWY7vGD$_ipI90J8~AYac96*a#=_exf$2eBpF420s@5~(J$wB zTkGs%9oF(S_j<#Px3v}RQm%ZFHszc6RT1`SVU!*Pu(!kE(q_uEtN8cV@*g-Yx92r* z98VtPOuhF_fv_jl5;SaMCZm*%$pZT39aBeAqp>Vpno5w#B9) zG@Oj10pon>EQ5nfWq)`1Q$;jVYn;dd<6xD_X&H6BaqO|r%ruI+ zi?eOk^F1dK$z+*Tr%vet-Im| zrJMUI>6j?~PcND|dmOd1K;$2h`O45|E_~8u%KG3V;ldM$tKx!qKy`){J1F~C4$4}u-9wIyZWKZ-e3mRY zP|`c2VC&e;uOSx2XAz5dW4*;&E$A!NS5SYo`8sqa_yS|q0t zz^XnFGQ#@gj@7>1x(e&*NwnQ9nkA8STR#I@xTxTu^5FVu+LVOs#98fF{JZ*yZiSrr zuM%z|K(4wVIRY@S-WMPX4u{CG&{hhkP9wYHw-m$sjfJM08n6IJHoYAGDvSTr zDdq)cB5&`ovB2pELIUWy**~a zs#MBilUjb|tE8nYs29Hx4-@&t-ez*n_4}hA_mVSpI5BMV`c4TNW$I#Eol2%j2ju2Y z6WRH+S%k{`S~K<4XCI-m)F|KidqAvV2{A(NL+77qlWg# za!e+K1HH+l`NWA_$kOuDXcDqU?M$QVDI`ru;=EIfXgLogg6hKn;5Bd~o7$hX*xOvr ztB>pwq;JT2vq8Kfov-QD9ZcgssgH9`4` zMT|*mEbT%AWqy2<5-T!6KA3V{oEk%;SDK7t=b`J716|_Gw7N3d;@qEDq}hxDd4>R|@^%Z05Rwu;PrQ@HmkRQoHzQmrRxoozXR1Sih5s)}Afp z(Gia3{!y<=Vq#+xho*^BxO2gAE^+r_Uy4v{kH0`oHLx&ybOP|0!nV(x#a|Pt$Z(HF zl$HF}^r}+Y?4LD(Q1`|*S|mzB5L`E8yb=sgnwZj+D@a!BIO`z5%teA47OHkEj zsm0!L(J#Ad%_k8zHC<_3k-w?7wtUA>3VrZ!aAotW2Z`qn`RHfH5UxhYOoV%@ztGv0 z5lx`&8(hWbV-Sw{$dyenQ{Dv$a^+kt8!N#4?~bY5zZno79u-lCFg{qDQ=>@FtItbP zI(gadVJP8#bZtWQ9y7TP{K6^>>pG??Wx25&(sKnV&IYS!v>@pyFGHelJsL9Oq+TFi)RV! z-~6&oXrcdtKb1sLGL;R)WueLZ%IZl1m@ALEoX;J$GOh6(P9_{PL=QzAd~4an@kL_^ z*|8Ifftzbyt-^|q!I6Z~QpDdMl%`Jl!a&+|;xjZbI$IB>LYN{03M{778Wu4AX2g!4 zATY#^JWn4~gy-!gjm}cs+SO0g61vsLLh2$Aq^))O@-E7!3xDbe#<7_zHSF^xt|LBC z7yMo8t&8Fb^zK5zyikJom1m{_una|XR-`e-;uH=T#9l7!2!s~XHLyw9p|!g*UJDa| zVShIai$9+Hp898CV>7;yRoQNCyO9Cn(aU}<>MCNNAbhRZnetG1jx7`Qp$UPb2&(E=x{*f0rJ8mZ4y^K$T=|&sdMkF^@f|_W%z0x`{QqJ;i~w2n`Vz3 z^WMJ?y=Y~>uw;Pt@W;IPhnRln_Ee!9wFEhM^hJ|KkZh$-AX_<6Q`zADo93tvxw$aa zQ=@lxKK0=hzFt>fq(?MuUY1K>F8_YPL#yqbWg$xr8Pb!Uk?=Nbl&(*c;ozpfjk%Bs zt|3rIQO&`Hx$6&2^{MA8v-U_S#+YrJauSjzdu-`Q(A`Ry-Gsy`)v&x&xfW%gGss-McWN)ihBTgF z|6GGX+|HRU1ZWo1M4Xs^#>1IoE_oK@L^q^L+p`;@%+t2^YTOU7F3&-FtO}bPoRd`o zuj@)>4k8Xz-dH@_w5f}-)5E%?4m=D@cTz27tq5D{`i#jMvc=9v+8oFr)!uDtNSbSg zqqi8rzdkC-R9sZK>8fx|?lNAC`bAXh(b??Weo6~bZ8Vyu(v!sI{?1se|7`!w5|70M zm%?U-sF5~p^yM8+7KwPVIjmY;a4W)sQNDT+DURTx&D+}baHIFdpZuqde;H$Eh6q+x z=hVj#uz<3`akT+RH4qo`+6}!OOy5U(| zLRK^psiIkT{E@;^Jghc$QgU34d#o10-Qs6P&<#e$bimp*xYe)GZ1A=3^o0g#o%)qh ztQ0y!YPHu3of&Cih^on;*zt`0swDO26Iy5Da zFa}ooANvqV`c@ub9xcGh(%|GgAe1{EMQff7C%i=^N4pf-BD34mICSsKsq&ZZjY1@Z zG@cskez&-3LR{LSb>0w8%lLLeeTrfQJid;XIgr}h*c*>u8vpd{r(jOBOun1DG?99P zAJyGtNadXNbG>XDx-&U_9h_@aeln5&8ByV3v!2yPOR9>SfKc2elbbvGJQyYylua2{ zePSD<0ULNOAc%nKh#r?3;J!3jKqTto^KD)!HLq_li z@*sZy1P>c2`P2xE3$raNHmv(-Xpt3vx!=KqfBiH6qF=3Xi9fJ)Sk_wNNVS7@*^N1G zwm?+ayM{?yGM#h|zPZGxl*ny=U>i3hmDO^eS%-l(m1BUB7_T%)IDhoWz_=nZTv~xs z=->u%W5CCuZ$y*SSSY8->w@hMN?@^>;^ZPo#m>9+$`Xgf0U>?imG!&Y%dHQCdA+;;Di+=sEIyKjlN2=4uQT9)+zC*O&ORAe%fN_wT-1)uY ztdOb0J$Q2J%`-LU(>J=aB*00%XMk0hYpcFqYxhL#`({Nb$=j#jeyEo2)7gL;fmag~it@v@EJ2pw*i9g$8AS&<<2|G6nS*-Mh*&4@Dh8Y2k0`x6%J#1Jmgcv>p6R-<$; zwnuJ%qs-nRyJv2tOh?Dl(+Cf9)F41@LGQCQDWL`Oi_K^X@`2mRGE2iUX zI*0t^Bv$xAdoG3ZdS;>)%1|YjPPSl5#MYxzk^Fgys-6epRHp=A^axQs+S|x!{RJo@ zK2Q5idVMrAq8{_9hF9@eW^WTZ6pgzujXcCiGcxH4+D=|3O>LU~om9)h4G^r&9V+3(%5BPMP6**e;b#VqW`!mVD6GPEL(bhLd=Q()y)xh|W_zR7yij#NZxP@UkfJ zP{ZI%y4npoq}E#Q)j5L>-PUm9`f;f~yQm=|Md2Z7p7X1)T}%&7V)!g_IK}3C%x4nc z+G2zG)-ssC?960L+ zB~lnCuXZV#4`Y~@p%vF9un=50&u9^CB3fQY_UchXekXCqP-Nq4YRGGn$aO`ouJbgn zVWSQ^Tw2%ZI*F7pY#Uo=)7Tka(mOkF7ba?pPPE2$^=)AvhM)J)`dqlNN?CoaJ|4M# z%jjC%U8UO1pLfHnUIW!hQN8LrcUCM56 zFPHKOOqa5E;x6tuhP$|pKKw-n?Igc>Nxau;Sp$s)+vjUTMy-Y$DxOBH1;hdtB5T(^ zre!s(q1Wq8PmENDju%fTd4)MQ!|LQGwRIT3ICcMKaW*-*m0>0_P>x5kra`mTq1$B| zE4Gp+h$y57I^gR+@qsnM#exE$c(5`sBKr=r)!99> zJ4UIy4&#Ng?cnz>FBgef|A2WjO^MF&NdD$@p;>(gtB&eO>+`WWGZ)&jR-o#2CDl4k zSF(=uv~3cGfH)3yAZk6ku;X}{mc`aV!{%md%n!EJ<`mKkOEL7s@d|*>njh%w3$UV` zTAf$9tDZqPI4Nno>hTVnT_5UlI~t8eZ{Kh3_ENPQ+_hxfgji0;qSc;mi9{WVS`mTN z%cz@;*7uLq!u;}YTvo-?`i#vNVvTM%?zQSe$qXca|JY^gr6PJ=f3i0fcX5V93j*5Z z*C(Y6r;e4fjZ5o1WgJ^>hhe*>T^-8v5C7Cc-d4jAwk~tZc*#0K=4gpl7imVoqZCN0 z^Kci}lL^8b78nvzY6(~PmMx&h=P=|avn_MSNucKcv*(>y z&MeF6iVlupy26%ZH~Zo)(oGUH=`JVr2@Q+Xx?Ym~A{5$_`0RO_biZI#yw03!#fJM` zb0SK;(L!v-TDPiRS7Wf15P#?kL}(d}Y;nrtPLZ z>rue)%cstDf8|#*0v&GcEv=;jh)pef67*_uv z;xqMwSE-WBLy73>S*4c9)G%qm{4j@0bUvR%(e_!^0GVktFDo;=1>_8won1^seAe8D zLo{A@q%)%FO*!DmX}%&(OYzr5ST+w+x^FGUrVu zp05|eW_XhGQ6U}`Ky?2uZqvD+OQ6#vWRjtmpz9ZJ8xEg#_`nx*{qGkQOz~~-AFV_E zeTRm%@D2BErSJ{%CpVPne5fAmP82h9iz{rVftMYcemR3PiF8ZW6eJm zYaiCk*KeqQ@PPLIK_QBfI<1(8!XJGH`Rm-!yi2IRU%0}TbW*nm#X=P!e*R%$0l|Yh zMfgW>b-`ZeYu$Yn0lxhM{liFWy5F}{f1iBAjQpj0tAN08e}$(s%s&tWlPdfH-td5s zU}<I#xxUt&+v^C%7>*QL8!Yt=mIl=t#r!`YSlZj8Vf`9XkA{uO7Pz<+>L zjkZ*-QV+gyZE0Aze`r_%_o?flb)QC5I0f`(7~Ic4C^USeGz?pgTR}IdLvEmB=|TSC hz5~D=w9QD*l+{M5`m>d;6X@FP{{f=cjDRGB0sz{GFv0); literal 0 HcmV?d00001 diff --git a/tests/testthat.R b/tests/testthat.R index 948ef9b74..8508877e9 100644 --- a/tests/testthat.R +++ b/tests/testthat.R @@ -64,6 +64,20 @@ pest_dda <- findChromPeaks(pest_dda, param = cwp) ## sciex_data <- readMSData(fl, mode = "onDisk") ## sciex_data <- pickPeaks(sciex_data) +## IM data +if(system.file("data/IM_test_data.rds", package = "xcms") != ""){ + library(MsExperiment) + im_spec <- readRDS(system.file("data/IM_test_data.rds", package = "xcms")) + im_exp <- MsExperiment() + spectra(im_exp) <- im_spec + sampleData(im_exp) <- DataFrame( + raw_file = "testfile" + ) + im_exp <- linkSampleData(im_exp, + with = "sampleData.raw_file = spectra.dataOrigin") + +} + test_check("xcms") bpstop(prm) diff --git a/tests/testthat/test_Param_classes.R b/tests/testthat/test_Param_classes.R index a20ff3554..1a133ee00 100644 --- a/tests/testthat/test_Param_classes.R +++ b/tests/testthat/test_Param_classes.R @@ -482,6 +482,121 @@ test_that("CentWavePredIsoParam works", { expect_equal(L$snthresh, 123) }) +test_that("IMCentWaveParam works", { + skip_on_os(os = "windows", arch = "i386") + + ## Check getter/setter methods: + p <- new("IMCentWaveParam", ppm = 14) + expect_equal(ppm(p), 14) + ppm(p) <- 13 + expect_equal(ppm(p), 13) + p <- IMCentWaveParam(ppm = 21) + expect_equal(ppm(p), 21) + expect_error(IMCentWaveParam(ppm = -4)) + + p <- new("IMCentWaveParam", peakwidth = c(1, 2)) + expect_equal(peakwidth(p), c(1, 2)) + peakwidth(p) <- c(2, 3) + expect_equal(peakwidth(p), c(2, 3)) + p <- IMCentWaveParam(peakwidth = c(3, 4)) + expect_equal(peakwidth(p), c(3, 4)) + expect_error(IMCentWaveParam(peakwidth = 1:3)) + + p <- new("IMCentWaveParam", snthresh = 5) + expect_equal(snthresh(p), 5) + snthresh(p) <- 6 + expect_equal(snthresh(p), 6) + p <- IMCentWaveParam(snthresh = 7) + expect_equal(snthresh(p), 7) + expect_error(IMCentWaveParam(snthresh = 1:3)) + + p <- new("IMCentWaveParam", prefilter = c(1, 2)) + expect_equal(prefilter(p), c(1, 2)) + prefilter(p) <- c(2, 3) + expect_equal(prefilter(p), c(2, 3)) + p <- IMCentWaveParam(prefilter = c(3, 4)) + expect_equal(prefilter(p), c(3, 4)) + expect_error(IMCentWaveParam(prefilter = 1:3)) + + p <- new("IMCentWaveParam", mzCenterFun = "meanApex3") + expect_equal(mzCenterFun(p), "meanApex3") + mzCenterFun(p) <- "mean" + expect_equal(mzCenterFun(p), "mean") + p <- IMCentWaveParam(mzCenterFun = "mean") + expect_equal(mzCenterFun(p), "mean") + expect_error(IMCentWaveParam(mzCenterFun = "median")) + + p <- new("IMCentWaveParam", integrate = 2L) + expect_equal(integrate(p), 2L) + integrate(p) <- 1L + expect_equal(integrate(p), 1L) + p <- IMCentWaveParam(integrate = 2L) + expect_equal(integrate(p), 2L) + expect_error(IMCentWaveParam(integrate = 3L)) + + p <- new("IMCentWaveParam", mzdiff = 1) + expect_equal(mzdiff(p), 1) + mzdiff(p) <- 2 + expect_equal(mzdiff(p), 2) + p <- IMCentWaveParam(mzdiff = 3) + expect_equal(mzdiff(p), 3) + expect_error(IMCentWaveParam(mzdiff = 2:3)) + + p <- new("IMCentWaveParam", fitgauss = TRUE) + expect_equal(fitgauss(p), TRUE) + fitgauss(p) <- FALSE + expect_equal(fitgauss(p), FALSE) + p <- IMCentWaveParam(fitgauss = TRUE) + expect_equal(fitgauss(p), TRUE) + expect_error(IMCentWaveParam(fitgauss = c(TRUE, FALSE))) + + p <- new("IMCentWaveParam", noise = 3) + expect_equal(noise(p), 3) + noise(p) <- 6 + expect_equal(noise(p), 6) + p <- IMCentWaveParam(noise = 8) + expect_equal(noise(p), 8) + expect_error(IMCentWaveParam(noise = c(3, 5))) + + p <- new("IMCentWaveParam", verboseColumns = TRUE) + expect_equal(verboseColumns(p), TRUE) + verboseColumns(p) <- FALSE + expect_equal(verboseColumns(p), FALSE) + p <- IMCentWaveParam(verboseColumns = TRUE) + expect_equal(verboseColumns(p), TRUE) + expect_error(IMCentWaveParam(verboseColumns = c(TRUE, FALSE))) + + p <- new("IMCentWaveParam", firstBaselineCheck = TRUE) + expect_equal(firstBaselineCheck(p), TRUE) + firstBaselineCheck(p) <- FALSE + expect_equal(firstBaselineCheck(p), FALSE) + p <- IMCentWaveParam(firstBaselineCheck = FALSE) + expect_equal(firstBaselineCheck(p), FALSE) + expect_error(IMCentWaveParam(firstBaselineCheck = c(TRUE, FALSE))) + + p <- new("IMCentWaveParam", ppmMerging = 10) + expect_equal(ppmMerging(p), 10) + ppmMerging(p) <- 20 + expect_equal(ppmMerging(p), 20) + p <- IMCentWaveParam(ppmMerging = 30) + expect_equal(ppmMerging(p), 30) + expect_error(IMCentWaveParam(ppmMerging = c(10, 60))) + + p <- new("IMCentWaveParam", binWidthIM = 0.03) + expect_equal(binWidthIM(p), 0.03) + binWidthIM(p) <- 0.05 + expect_equal(binWidthIM(p), 0.05) + p <- IMCentWaveParam(binWidthIM = 0.01) + expect_equal(binWidthIM(p), 0.01) + expect_error(IMCentWaveParam(binWidthIM = c(0.01, 0.08))) + + ## Check the .param2list method: + p <- new("IMCentWaveParam", snthresh = 123) + L <- .param2list(p) + expect_equal(L$snthresh, 123) +}) + + test_that("PeakDensityParam works", { skip_on_os(os = "windows", arch = "i386") diff --git a/tests/testthat/test_do_findChromPeaks-functions.R b/tests/testthat/test_do_findChromPeaks-functions.R index 94aebf335..0de569b1c 100644 --- a/tests/testthat/test_do_findChromPeaks-functions.R +++ b/tests/testthat/test_do_findChromPeaks-functions.R @@ -88,6 +88,48 @@ test_that("do_findChromPeaks_massifquant works", { expect_true(nrow(res_3) < nrow(res_2)) }) +test_that("do_findChromPeaks_IM_centWave works", { + skip_on_os(os = "windows", arch = "i386") + skip_if(system.file("data/IM_test_data.rds", package = "xcms") == "") + res_1 <- do_findChromPeaks_IM_centWave(im_spec, + peakwidth = c(1,5), + snthresh = 0, + ppmMerging = 50) + res_2 <- do_findChromPeaks_IM_centWave(im_spec, + peakwidth = c(1,5), + snthresh = 0, + ppmMerging = 10) + expect_true(nrow(res_2) > nrow(res_1)) # Using a higher ppmMerging yields less peaks +}) + +test_that("extract_mobilogram works", { + skip_on_os(os = "windows", arch = "i386") + skip_if(system.file("data/IM_test_data.rds", package = "xcms") == "") + mobilogram <- .extract_mobilogram(peaksData(im_spec), + peak = data.frame(mzmin = 286.1982, mzmax = 286.2051, + rtmin = 357.0347, rtmax = 359.838), + rt = rtime(im_spec), + im = im_spec$inv_ion_mobility) + expect_length(mobilogram, 2) + expect_equal(range(mobilogram$x), c(0, 12360727)) + expect_equal(range(mobilogram$mids), c(0.005, 0.995)) +}) + + +test_that("split_mobilogram_CWT works", { + skip_on_os(os = "windows", arch = "i386") + skip_if(system.file("data/IM_test_data.rds", package = "xcms") == "") + mobilogram <- .extract_mobilogram(peaksData(im_spec), + peak = data.frame(mzmin = 286.1982, mzmax = 286.2051, + rtmin = 357.0347, rtmax = 359.838), + rt = rtime(im_spec), + im = im_spec$inv_ion_mobility) + split <- .split_mobilogram_CWT(mobilogram) + expect_length(split, 1) + expect_equal(split[[1]], c(0.795, 0.895)) +}) + + test_that("do_findChromPeaks_matchedFilter works", { skip_on_os(os = "windows", arch = "i386") From 120627a39d64e13bb95ff68f98052fb58e0d68e5 Mon Sep 17 00:00:00 2001 From: RogerGinBer <45530914+RogerGinBer@users.noreply.github.com> Date: Wed, 11 Jan 2023 15:50:33 +0100 Subject: [PATCH 8/8] refactor: simplify split_mobilogram_CWT Simplify algorithm for split_mobilgram_CWT: instead of a parametric CWT, use a much simpler/robust `localMaxima` approach to detect apexes and `descentMinTol` to determine the regions. Add start-apex-end information in output vector. Change tests accordingly. - Also, fix bug in `extract_mobilogram`, where bin required the y values to be ordered (mentioned in https://github.com/rformassspectrometry/MsCoreUtils/issues/108) --- R/do_findChromPeaks-functions.R | 33 +++++++------------ .../test_do_findChromPeaks-functions.R | 7 ++-- 2 files changed, 16 insertions(+), 24 deletions(-) diff --git a/R/do_findChromPeaks-functions.R b/R/do_findChromPeaks-functions.R index dd3dc2783..39357f5b9 100644 --- a/R/do_findChromPeaks-functions.R +++ b/R/do_findChromPeaks-functions.R @@ -2316,7 +2316,6 @@ do_findChromPeaks_IM_centWave <- function(spec, rt <- rtime(spec) im <- spec$inv_ion_mobility - ## Resolving peaks across IM dimension resolved_peaks <- vector("list", nrow(peaks)) for (i in seq_len(nrow(peaks))) { @@ -2327,7 +2326,7 @@ do_findChromPeaks_IM_centWave <- function(spec, next } - bounds <- .split_mobilogram_CWT(mobilogram) + bounds <- .split_mobilogram(mobilogram) new_peaks <- data.frame( mz = current_peak["mz"], mzmin = current_peak["mzmin"], @@ -2335,9 +2334,9 @@ do_findChromPeaks_IM_centWave <- function(spec, rt = current_peak["rt"], rtmin = current_peak["rtmin"], rtmax = current_peak["rtmax"], - im = vapply(bounds, mean, numeric(1)), - immin = vapply(bounds, min, numeric(1)), - immax = vapply(bounds, max, numeric(1)), + im = vapply(bounds, function(x) x[[2]], numeric(1)), + immin = vapply(bounds, function(x) x[[1]], numeric(1)), + immax = vapply(bounds, function(x) x[[3]], numeric(1)), row.names = NULL ) resolved_peaks[[i]] <- new_peaks @@ -2399,34 +2398,26 @@ do_findChromPeaks_IM_centWave <- function(spec, ints <- vapply(pdata[keep], xcms:::.aggregate_intensities, mzr = mzr, INTFUN = sum, na.rm = TRUE, numeric(1)) if(all(ints == 0)) return() - mob <- MsCoreUtils::bin(x = ints, y = ims, size = binWidthIM, FUN = sum) + mob <- MsCoreUtils::bin(x = ints[order(ims)], y = sort(ims), + size = binWidthIM, FUN = sum) mob } -#' @importFrom MassSpecWavelet peakDetectionCWT -.split_mobilogram_CWT <- function(mob){ +.split_mobilogram <- function(mob){ if(length(mob$x) == 0){return()} vec <- mob$x - #Add some padding, which will be removed after - padding_size <- 5 - vec <- c(rep(0, padding_size), vec, rep(0, padding_size)) - pks <- MassSpecWavelet::peakDetectionCWT(vec, scales = c(1:7)) - left <- sapply(pks$majorPeakInfo$peakCenterIndex - pks$majorPeakInfo$peakScale, function(x) max(1, x)) - right <- sapply(pks$majorPeakInfo$peakCenterIndex + pks$majorPeakInfo$peakScale, function(x) min(x, length(vec) - padding_size)) + apex <- which(MsCoreUtils::localMaxima(vec, hws = 4)) limits <- list() - for (i in seq_along(pks$majorPeakInfo$peakCenterIndex)){ - ranges <- xcms:::descendMinTol(vec, startpos = c(left[i], right[i]), maxDescOutlier = 1) - padding_size - ranges[1] <- min(max(1, ranges[1]), length(mob$mids)) - ranges[2] <- min(ranges[2], length(mob$mids)) - limits[[i]] <- mob$mids[ranges] + for (i in seq_along(apex)){ + ranges <- descendMinTol(vec, startpos = c(apex[i], apex[i]), maxDescOutlier = 2) + limits[[i]] <- mob$mids[c(ranges[1], apex[i], ranges[2])] } limits <- limits[vapply(limits, function(x){!any(is.na(x))}, logical(1))] - return(limits) + limits } - #' @importFrom dplyr between .extract_EIC_IM <- function(peak, pdata, rt, im){ rtr <- c(peak["rtmin"], peak["rtmax"]) diff --git a/tests/testthat/test_do_findChromPeaks-functions.R b/tests/testthat/test_do_findChromPeaks-functions.R index 0de569b1c..e708ab2a2 100644 --- a/tests/testthat/test_do_findChromPeaks-functions.R +++ b/tests/testthat/test_do_findChromPeaks-functions.R @@ -116,7 +116,7 @@ test_that("extract_mobilogram works", { }) -test_that("split_mobilogram_CWT works", { +test_that("split_mobilogram works", { skip_on_os(os = "windows", arch = "i386") skip_if(system.file("data/IM_test_data.rds", package = "xcms") == "") mobilogram <- .extract_mobilogram(peaksData(im_spec), @@ -124,9 +124,10 @@ test_that("split_mobilogram_CWT works", { rtmin = 357.0347, rtmax = 359.838), rt = rtime(im_spec), im = im_spec$inv_ion_mobility) - split <- .split_mobilogram_CWT(mobilogram) + split <- .split_mobilogram(mobilogram) expect_length(split, 1) - expect_equal(split[[1]], c(0.795, 0.895)) + expect_length(split[[1]], 3) + expect_equal(split[[1]], c(0.795, 0.825, 0.865)) })