Merge pull request #761 from sneumann/issue_755

Some updates and bug fixes

Author: sneumann

.github/workflows/check-bioc.yml

@@ -124,32 +124,7 @@ jobs:
       - name: Install macOS system dependencies
         if: matrix.config.os == 'macOS-latest'
         run: |
-          ## Enable installing XML from source if needed
-          brew install libxml2
-          echo "XML_CONFIG=/usr/local/opt/libxml2/bin/xml2-config" >> $GITHUB_ENV
-
-          ## Required to install magick as noted at
-          ## https://github.com/r-lib/usethis/commit/f1f1e0d10c1ebc75fd4c18fa7e2de4551fd9978f#diff-9bfee71065492f63457918efcd912cf2
-          brew install imagemagick@6
-
-          ## For textshaping, required by ragg, and required by pkgdown
-          brew install harfbuzz fribidi
-
-          ## For installing usethis's dependency gert
-          brew install libgit2
-
-          ## required for ncdf4
-          ## brew install netcdf ## Does not work as it is compiled with gcc
-          ## Use pre-compiled libraries from https://mac.r-project.org/libs-4/
-          # curl -O https://mac.r-project.org/libs-4/netcdf-4.7.4-darwin.17-x86_64.tar.gz
-          # tar fvxzm netcdf-4.7.4-darwin.17-x86_64.tar.gz -C /
-          # rm netcdf-4.7.4-darwin.17-x86_64.tar.gz
-          # curl -O https://mac.r-project.org/libs-4/hdf5-1.12.0-darwin.17-x86_64.tar.gz
-          # tar fvxzm hdf5-1.12.0-darwin.17-x86_64.tar.gz -C /
-          # rm hdf5-1.12.0-darwin.17-x86_64.tar.gz
-          # curl -O https://mac.r-project.org/libs-4/szip-2.1.1-darwin.17-x86_64.tar.gz
-          # tar fvxzm szip-2.1.1-darwin.17-x86_64.tar.gz -C /
-          # rm szip-2.1.1-darwin.17-x86_64.tar.gz
+        shell: Rscript {0}
 
       - name: Install Windows system dependencies
         if: runner.os == 'Windows'

DESCRIPTION

@@ -1,5 +1,5 @@
 Package: xcms
-Version: 4.3.2
+Version: 4.3.3
 Title: LC-MS and GC-MS Data Analysis
 Description: Framework for processing and visualization of chromatographically
     separated and single-spectra mass spectral data. Imports from AIA/ANDI NetCDF,
@@ -54,7 +54,7 @@ Imports:
     methods,
     Biobase,
     BiocGenerics,
-    ProtGenerics (>= 1.35.4),
+    ProtGenerics (>= 1.37.1),
     lattice,
     MassSpecWavelet (>= 1.66.0),
     S4Vectors,
@@ -63,7 +63,7 @@ Imports:
     MsCoreUtils (>= 1.15.5),
     MsFeatures,
     MsExperiment (>= 1.5.4),
-    Spectra (>= 1.13.7),
+    Spectra (>= 1.15.7),
     progress,
     RColorBrewer,
     MetaboCoreUtils (>= 1.11.2)

NAMESPACE

@@ -4,7 +4,8 @@ import("methods")
 importMethodsFrom("ProtGenerics", "peaks", "chromatogram", "writeMSData",
     "polarity<-", "centroided", "isCentroided", "peaks<-",
     "isolationWindowTargetMz", "quantify", "bin", "spectrapply",
-    "filterFeatures", "filterMzRange", "filterRt", "filterMz", "filterMsLevel")
+    "filterFeatures", "filterMzRange", "filterRt", "filterMz", "filterMsLevel",
+    "estimatePrecursorIntensity")
 importClassesFrom("ProtGenerics", "Param")
 importFrom("BiocGenerics", "updateObject", "fileName", "subset",
     "dirname", "dirname<-")
@@ -250,7 +251,6 @@ export(
     "hasFilledChromPeaks",
     "plotAdjustedRtime",
     "groupOverlaps",
-    "estimatePrecursorIntensity",
     "featureArea",
     "loadXcmsData",
     "matchLamasChromPeaks",
@@ -263,6 +263,7 @@ export(
 exportMethods(
     "showError",
     "findChromPeaks",
+    "estimatePrecursorIntensity",
     "groupChromPeaks",
     "adjustRtime",
     "findChromPeaksIsolationWindow",
@@ -568,6 +569,7 @@ importMethodsFrom("Spectra", "precursorMz")
 importMethodsFrom("Spectra", "$")
 importMethodsFrom("Spectra", "uniqueMsLevels")
 importMethodsFrom("Spectra", "backendBpparam")
+importMethodsFrom("Spectra", "estimatePrecursorIntensity")
 importFrom("Spectra", "MsBackendMemory", "filterEmptySpectra")
 
 ## MsExperiment things

NEWS.md

@@ -1,5 +1,15 @@
 # xcms 4.3
 
+## Changes in version 4.3.3
+
+- Fix issue #755: `chromatogram()` with `msLevel = 2` fails to extract
+  chromatographic data if `isolationWindowTargetMz` is not specified or
+  available (e.g. for MSe data).
+- Support coercing from `XcmsExperiment` to `XCMSnExp` with
+  `as(object, "XCMSnExp")`.
+- Change `estimatePrecursorIntensity()` to a method and add an implementation
+  for `MsExperiment` objects.
+
 ## Changes in version 4.3.2
 
 - Remove data/results import/export functionality as it is being developed in

R/MsExperiment-functions.R

@@ -462,7 +462,7 @@
     if (length(msLevel) != npks)
         msLevel <- rep(msLevel[1L], npks)
     if (!length(isolationWindow))
-        isolationWindow <- rep(1L, npks)
+        isolationWindow <- rep(NA_real_, npks)
     if (length(isolationWindow) && length(isolationWindow) != npks)
         stop("Length of 'isolationWindow' (if provided) should match the ",
              "number of chromatograms to extract.")

R/MsExperiment.R

@@ -114,3 +114,16 @@ setMethod(
             msLevel = msLevel, isolationWindow = isolationWindowTargetMz,
             chunkSize = chunkSize, BPPARAM = BPPARAM)
     })
+
+#' @rdname estimatePrecursorIntensity
+setMethod(
+    "estimatePrecursorIntensity", "MsExperiment",
+    function(object, ppm = 10, tolerance = 0,
+             method = c("previous", "interpolation"),
+             BPPARAM = bpparam()) {
+        method <- match.arg(method)
+        estimatePrecursorIntensity(spectra(object), ppm = ppm,
+                                   tolerance = tolerance, method = method,
+                                   f = spectraSampleIndex(object),
+                                   BPPARAM = BPPARAM)
+})

R/XcmsExperiment-functions.R

@@ -1136,6 +1136,78 @@ featureArea <- function(object, mzmin = min, mzmax = max, rtmin = min,
 }
 
 #' function to create an empty `XcmsExperiment` object
+#'
+#' @noRd
 XcmsExperiment <- function() {
     as(MsExperiment(), "XcmsExperiment")
 }
+
+#' function to convert an XcmsExperiment into an XCMSnExp.
+#'
+#' @author Johannes Rainer
+#'
+#' @noRd
+.xcms_experiment_to_xcms_n_exp <- function(from) {
+    requireNamespace("MSnbase", quietly = TRUE)
+    n <- new("XCMSnExp")
+    if (!length(from))
+        return(n)
+    m <- new("MsFeatureData")
+    ## works only if we have a MsBackendMzR
+    if (!inherits(spectra(from)@backend, "MsBackendMzR"))
+        stop("Can not convert 'from' to a 'XCMSnExp' object: 'spectra(x)' uses",
+             "an MS backend other than 'MsBackendMzR'. Currently, only",
+             " coercion of an object with a 'MsBackendMzR' is supported.")
+    ## works only if we have an empty processing queue
+    if (length(spectra(from)@processingQueue) > 0)
+        stop("Can not convert 'from' to a 'XCMSnExp' object: processing queue",
+             " of the Spectra object is not empty.")
+    ## -> OnDiskMSnExp
+    n@processingData <- new("MSnProcess",
+                            processing = paste0("Data converted [", date(), "]"),
+                            files = fileNames(from),
+                            smoothed = NA)
+    n@phenoData <- new("NAnnotatedDataFrame", as.data.frame(sampleData(from)))
+    fd <- as.data.frame(from@spectra@backend@spectraData)
+    fnames <- unique(fd$dataStorage)
+    fd$fileIdx <- match(fd$dataStorage, fnames)
+    fd <- fd[, !colnames(fd) %in% c("dataStorage", "dataOrigin")]
+    colnames(fd) <- sub("scanIndex", "spIdx", colnames(fd))
+    colnames(fd) <- sub("rtime", "retentionTime", colnames(fd))
+    colnames(fd) <- sub("precScanNum", "precursorScanNum", colnames(fd))
+    colnames(fd) <- sub("precursorMz", "precursorMZ", colnames(fd))
+    fd$spectrum <- seq_len(nrow(fd))
+    rownames(fd) <- MSnbase:::formatFileSpectrumNames(
+                                  fd$fileIdx, fd$spIdx,
+                                  max(fd$fileIdx), max(fd$spIdx))
+    n@featureData <- new("AnnotatedDataFrame", fd)
+    nf <- length(fnames)
+    n@experimentData <- new("MIAPE",
+                            instrumentManufacturer = rep(NA_character_, nf),
+                            instrumentModel = rep(NA_character_, nf),
+                            ionSource = rep(NA_character_, nf),
+                            analyser = rep(NA_character_, nf),
+                            detectorType = rep(NA_character_, nf))
+    n@processingData <- new("MSnProcess",
+                            processing = paste0("Coercion from ",
+                                                class(from)[1L],
+                                                " [", date(), "]"),
+                            files = fnames)
+    ## -> XCMSnExp
+    if (hasChromPeaks(from)) {
+        chromPeaks(m) <- chromPeaks(from)
+        chromPeakData(m) <- chromPeakData(from)
+    }
+    if (hasAdjustedRtime(from)) {
+        art <- fd$retentionTime_adjusted
+        names(art) <- rownames(fd)
+        adjustedRtime(m) <- split(art, fd$fileIdx)
+    }
+    if (hasFeatures(from))
+        featureDefinitions(m) <- DataFrame(featureDefinitions(from))
+    lockEnvironment(m, bindings = TRUE)
+    n@msFeatureData <- m
+    n@.processHistory <- from@processHistory
+    validObject(n)
+    n
+}

R/XcmsExperiment.R

@@ -124,14 +124,16 @@
 #'   associated feature definitions will be included in the returned
 #'   `XChromatograms`. By default the function returns chromatograms from MS1
 #'   data, but by setting parameter `msLevel = 2L` it is possible to e.g.
-#'   extract also MS2 chromatograms. For `msLevel` other than 1 it is in
-#'   addition important to also specify the `isolationWindowTargetMz` for which
-#'   MS2 data should be extracted (e.g. for SWATH data MS2 spectra are created
-#'   for different m/z isolation windows and the `isolationWindowTargetMz`
-#'   parameter allows to define from which of these the MS2 chromatogram
-#'   should be extracted.
-#'   Note that in future more efficient data structures for chromatographic
-#'   data will be available as well.
+#'   extract also MS2 chromatograms. By default, with parameter
+#'   `isolationWindowTargetMz = NULL` or `isolationWindowTargetMz = NA_real_`,
+#'   data from **all** MS2 spectra will be considered in the chromatogram
+#'   extraction. If MS2 data was generated within different m/z isolation
+#'   windows (such as e.g. with Scies SWATH data), the parameter
+#'   `isolationWindowTargetMz` should be used to ensure signal is only extracted
+#'   from the respective isolation window. The `isolationWindowTargetMz()`
+#'   function on the `Spectra` object can be used to inspect/list available
+#'   isolation windows of a data set. See also the xcms *LC-MS/MS vignette* for
+#'   examples and details.
 #'
 #' - `chromPeaks`: returns a `numeric` matrix with the identified
 #'   chromatographic peaks. Each row represents a chromatographic peak

R/functions-OnDiskMSnExp.R

@@ -722,7 +722,7 @@ setReplaceMethod("dirname", "OnDiskMSnExp", function(path, value) {
 #' par(mfrow = c(1, 2))
 #' plot(res, precursorIntensity(tmt))
 #' plot(res_2, precursorIntensity(tmt))
-.estimate_prec_intensity <- function(x, ppm = 10,
+.estimate_prec_intensity <- function(x, ppm = 10, tolerance = 0,
                                      method = c("previous", "interpolation")) {
     method <- match.arg(method)
     pmz <- precursorMz(x)
@@ -739,7 +739,8 @@ setReplaceMethod("dirname", "OnDiskMSnExp", function(path, value) {
             before_int <- numeric()
             if (length(before_idx)) {
                 sp <- sps[[before_idx[length(before_idx)]]]
-                before_idx <- closest(pmz[i], sp@mz, ppm = ppm, tolerance = 0,
+                before_idx <- closest(pmz[i], sp@mz, ppm = ppm,
+                                      tolerance = tolerance,
                                       duplicates = "closest")
                 if (!is.na(before_idx)) {
                     before_rt <- sp@rt
@@ -794,48 +795,6 @@ setReplaceMethod("dirname", "OnDiskMSnExp", function(path, value) {
     pmi
 }
 
-#' @title Estimate precursor intensity for MS level 2 spectra
-#'
-#' @description
-#'
-#' `estimatePrecursorIntensity` determines the precursor intensity for a MS 2
-#' spectrum based on the intensity of the respective signal from the
-#' neighboring MS 1 spectra (i.e. based on the peak with the m/z matching the
-#' precursor m/z of the MS 2 spectrum). Based on parameter `method` either the
-#' intensity of the peak from the previous MS 1 scan is used
-#' (`method = "previous"`) or an interpolation between the intensity from the
-#' previous and subsequent MS1 scan is used (`method = "interpolation"`, which
-#' considers also the retention times of the two MS1 scans and the retention
-#' time of the MS2 spectrum).
-#'
-#' @param x `OnDiskMSnExp` or `XCMSnExp` object.
-#'
-#' @param ppm `numeric(1)` defining the maximal acceptable difference (in ppm)
-#'     of the precursor m/z and the m/z of the corresponding peak in the MS 1
-#'     scan.
-#'
-#' @param method `character(1)` defining the method how the precursor intensity
-#'     should be determined (see description above for details). Defaults to
-#'     `method = "previous"`.
-#'
-#' @param BPPARAM parallel processing setup. See [bpparam()] for details.
-#'
-#' @return `numeric` with length equal to the number of spectra in `x`. `NA` is
-#'     returned for MS 1 spectra or if no matching peak in a MS 1 scan can be
-#'     found for an MS 2 spectrum
-#'
-#' @author Johannes Rainer
-#'
-#' @md
-estimatePrecursorIntensity <- function(x, ppm = 10,
-                                       method = c("previous", "interpolation"),
-                                       BPPARAM = bpparam()) {
-    method <- match.arg(method)
-    unlist(bplapply(.split_by_file2(x, subsetFeatureData = FALSE),
-                    .estimate_prec_intensity, ppm = ppm, method = method,
-                    BPPARAM = BPPARAM), use.names = FALSE)
-}
-
 #' Helper function to convert an OnDiskMSnExp to a Spectra object. This will
 #' only convert the spectra data, but no sample information.
 #'

R/functions-utils.R

@@ -791,8 +791,11 @@ groupOverlaps <- function(xmin, xmax) {
 #'     chromatograms should be extracted.
 #'
 #' @param pks_tmz `numeric` with the isolation window target m/z in which
-#'     the (MS2) chromatographic peak was detected. Does not need to be
-#'     provided for MS1 data.
+#'     the (MS2) chromatographic peak was detected. For `pks_msl > 1L` only
+#'     spectra with their `isolationWindowTargetMz` being equal to this value
+#'     are considered for the chromatogram extraction. Set to
+#'     `pks_tmz = NA_real_` to use **all** spectra with matching MS level and
+#'     ignore the isolation window.
 #'
 #' @param file_idx `integer(1)` allowing to optionally set the index of the
 #'     file the EIC is from (parameter `fromFile`).
@@ -803,8 +806,9 @@ groupOverlaps <- function(xmin, xmax) {
 #'
 #' @noRd
 .chromatograms_for_peaks <- function(pd, rt, msl, file_idx = 1L,
-                                     tmz = rep(1L, length(pd)), pks, pks_msl,
-                                     pks_tmz = rep(1L, nrow(pks)),
+                                     tmz = rep(NA_real_, length(pd)), pks,
+                                     pks_msl,
+                                     pks_tmz = rep(NA_real_, nrow(pks)),
                                      aggregationFun = "sum") {
     nr <- nrow(pks)
     pks_msl <- as.integer(pks_msl)
@@ -826,9 +830,9 @@ groupOverlaps <- function(xmin, xmax) {
         slot(res[[i]], "msLevel", check = FALSE) <- pks_msl[i]
         ## if pks_msl > 1: precursor m/z has to match!
         keep <- between(rt, pks[i, rtc]) & msl == pks_msl[i]
-        if (pks_msl[i] > 1L) {
+        if (pks_msl[i] > 1L && !is.na(pks_tmz[i])) {
             ## for DIA MS2: spectra have to match the isolation window.
-            keep <- keep & tmz == pks_tmz[i]
+            keep <- keep & tmz %in% pks_tmz[i]
         }
         keep <- which(keep)             # the get rid of `NA`.
         if (length(keep)) {