Skip to content

Commit

Permalink
refactor: little fixes
Browse files Browse the repository at this point in the history
  • Loading branch information
jorainer committed Sep 24, 2024
1 parent 13cfaf3 commit 156788a
Show file tree
Hide file tree
Showing 7 changed files with 96 additions and 42 deletions.
2 changes: 2 additions & 0 deletions NEWS.md
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,8 @@
an idealized bell shape (*beta*) during gap filling for centWave-based
chromatographic peak detection with parameter `verboseBetaColumns = TRUE`.
- Add `chromPeakSummary` generic (issue #705).
- Add `chromPeakSummary()` method to calculate the *beta* quality metrics.


## Changes in version 4.3.3

Expand Down
20 changes: 20 additions & 0 deletions R/AllGenerics.R
Original file line number Diff line number Diff line change
Expand Up @@ -551,6 +551,18 @@ setGeneric("chromPeakSpectra", function(object, ...)
#'
#' Currently implemented methods/parameter classes are:
#'
#' - `BetaDistributionParam`: calculates the *beta_cor* and *beta_snr* quality
#' metrics as described in (Kumler 2023) representing the result from a
#' (correlation) test of similarity to a bell curve and the signal-to-noise
#' ratio calculated on the residulas of this test.
#'
#' @param BPPARAM Parallel processing setup. See [bpparam()] for details.
#'
#' @param chunkSize `integer(1)` defining the number of samples from which data
#' should be loaded and processed at a time.
#'
#' @param msLevel `integer(1)` with the MS level of the chromatographic peaks
#' on which the metric should be calculated.
#'
#' @param object an *xcms* result object containing information on
#' identified chromatographic peaks.
Expand All @@ -569,7 +581,15 @@ setGeneric("chromPeakSpectra", function(object, ...)
#'
#' @author Pablo Vangeenderhuysen, Johannes Rainer
#'
#' @md
#'
#' @references
#'
#' Kumler W, Hazelton B J and Ingalls A E (2023) "Picky with peakpicking:
#' assessing chromatographic peak quality with simple metrics in metabolomics"
#' *BMC Bioinformatics* 24(1):404. doi: 10.1186/s12859-023-05533-4
#'
#' @export
setGeneric("chromPeakSummary", function(object, param, ...)
standardGeneric("chromPeakSummary"))

Expand Down
11 changes: 6 additions & 5 deletions R/XcmsExperiment.R
Original file line number Diff line number Diff line change
Expand Up @@ -2024,11 +2024,11 @@ setMethod(
if (!hasChromPeaks(object, msLevel = msLevel))
stop("No ChromPeaks definitions for MS level ", msLevel, " present.")
## Define region to calculate metrics from for each file
cp <- chromPeaks(object)
f <- factor(cp[,"sample"],seq_along(object))
pal <- split.data.frame(cp[,c("mzmin","mzmax","rtmin","rtmax")],f)
cp <- chromPeaks(object, msLevel = msLevel)
f <- factor(cp[,"sample"], seq_along(object))
pal <- split.data.frame(cp[, c("mzmin", "mzmax", "rtmin", "rtmax")], f)
names(pal) <- seq_along(pal)
## Manual chunk processing because we have to split `object` and `pal`
## Manual chunk processi ng because we have to split `object` and `pal`
idx <- seq_along(object)
chunks <- split(idx, ceiling(idx / chunkSize))
pb <- progress_bar$new(format = paste0("[:bar] :current/:",
Expand All @@ -2042,7 +2042,8 @@ setMethod(
.xmse_integrate_chrom_peaks(
.subset_xcms_experiment(object, i = z, keepAdjustedRtime = TRUE,
ignoreHistory = TRUE),
pal = pal[z], intFun = .chrom_peak_beta_metrics, BPPARAM = BPPARAM)
pal = pal[z], intFun = .chrom_peak_beta_metrics,
msLevel = msLevel, BPPARAM = BPPARAM)
})
res <- do.call(rbind, res)
pb$tick()
Expand Down
50 changes: 41 additions & 9 deletions man/chromPeakSummary.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

5 changes: 1 addition & 4 deletions tests/testthat/test_Param_classes.R
Original file line number Diff line number Diff line change
Expand Up @@ -1005,9 +1005,6 @@ test_that("FilterIntensityParam works", {


test_that("BetaDistributionParam works", {
skip_on_os(os = "windows", arch = "i386")

res <- BetaDistributionParam()
expect_true(is(res, "BetaDistributionParam "))

expect_true(is(res, "BetaDistributionParam"))
})
8 changes: 3 additions & 5 deletions tests/testthat/test_XcmsExperiment.R
Original file line number Diff line number Diff line change
Expand Up @@ -831,11 +831,11 @@ test_that(".chrom_peak_beta_metrics works", {
x <- Spectra::peaksData(spectra(xmse[2L]))
rt <- rtime(spectra(xmse[2L]))
pks <- chromPeaks(xmse)[chromPeaks(xmse)[, "sample"] == 2L, ]

res <- .chrom_peak_beta_metrics(x, rt, pks, sampleIndex = 2L,
cn = colnames(pks))
expect_equal(nrow(res), nrow(pks))

})

## That's from XcmsExperiment-functions.R
Expand Down Expand Up @@ -1458,8 +1458,6 @@ test_that("chromPeakSummary,XcmsExperiment works", {
verboseBetaColumns = FALSE)
xmse <- findChromPeaks(mse, param = p)
mat <- chromPeakSummary(xmse,BetaDistributionParam())
expect_true(all(c("beta_cor", "beta_snr") %in% colnames(res))).
expect_true(all(c("beta_cor", "beta_snr") %in% colnames(mat)))
expect_true(is.numeric(mat))

})

42 changes: 23 additions & 19 deletions vignettes/xcms.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -433,42 +433,46 @@ chromPeakData(faahko)

### Chromatographic peak quality

Based on the publication by Kumler et al. published in 2023 [@Kumler2023], new
quality metrics (`beta_cor` and `beta_snr`) were added to *xcms*. `beta_cor`
indicates how bell-shaped the peak is and `beta_snr` is estimating the
signal-to-noise ratio using the residuals. These metrics can be calculated
during peak picking by setting `verboseBetaColumns = TRUE` in the
`CentWaveParam` object, or they can becalculated afterwards by using the
`chromPeakSummary` function with the `XcmsExperiment` object and
`BetaDistributionParam` parameter object as input.
Based on the publication by Kumler et al. published in 2023 [@Kumler2023], new
quality metrics (*beta_cor* and *beta_snr*) were added to *xcms*. *beta_cor*
indicates how bell-shaped the chromatographic peak is and *beta_snr* is
estimating the signal-to-noise ratio using the residuals from this fit. These
metrics can be calculated during peak picking by setting `verboseBetaColumns =
TRUE` in the `CentWaveParam` object, or they can be calculated afterwards by
using the `chromPeakSummary()` function with the `XcmsExperiment` object and
the `BetaDistributionParam` parameter object as input:

```{r peak-detection-chromPeakSummary}
beta_metrics <- chromPeakSummary(faahko,BetaDistributionParam())
beta_metrics <- chromPeakSummary(faahko, BetaDistributionParam())
head(beta_metrics)
```

Using summary statistics, one can explore the distribution of these metrics in
The result returned by `chromPeakSummary()` is thus a numeric matrix with the
values for these quality estimates, one row for each chromatographic peak.
Using summary statistics, one can explore the distribution of these metrics in
the data.

```{r beta-metrics}
summary(beta_metrics)
```

Visual inspection gives a better idea of what these metrics represent in terms
of peak quality in the data at hand. This information can be used to e.g.
filter out peaks that don't meet a chosen quality metric threshold. In order to
plot the detected peaks, their EIC can be extracted with the function
`chromPeakChromatograms`. An example of a peak with a high `beta_cor` a peak
with a low `beta_cor` score is given below.
Visual inspection gives a better idea of what these metrics represent in terms
of peak quality in the data at hand. This information can be used to e.g.
filter out peaks that don't meet a chosen quality metric threshold. In order to
plot the detected peaks, their EIC can be extracted with the function
`chromPeakChromatograms`. An example of a peak with a high *beta_cor* and for
a peak with a low *beta_cor* score is given below.

```{r chromPeakChromatograms, message=FALSE}
eics <- chromPeakChromatograms(faahko)
beta_metrics[c(4, 6), ]
eics <- chromPeakChromatograms(
faahko, peaks = rownames(chromPeaks(faahko))[c(4, 6)])
```

```{r peak-quality-metrics, fig.widht = 10, fig.height = 5, fig.cap = "Plots of high and low quality peaks. Left: peak CP0004 with a beta_cor = 0.98, right: peak CP0006 with a beta_cor = 0.13."}
peak_1 <- eics[4]
peak_2 <- eics[6]
peak_1 <- eics[1]
peak_2 <- eics[2]
par(mfrow = c(1, 2))
plot(peak_1)
plot(peak_2)
Expand Down

0 comments on commit 156788a

Please sign in to comment.