Adding new peak quality metrics to findChromPeaks for CentWave

NAMESPACE

@@ -51,7 +51,7 @@ importFrom("stats", "aov", "approx", "convolve", "cor", "deriv3",
  "dist", "fft", "fitted", "lm", "loess", "lsfit", "median",
  "na.omit", "nextn", "nls", "predict", "pt", "quantile",
  "runmed", "sd", "stepfun", "weighted.mean", "density", "approxfun",
-  "rnorm", "runif")
+  "rnorm", "runif", "dbeta")
 importFrom("utils", "flush.console", "head", "object.size",
  "packageVersion", "read.csv", "tail", "write.csv",
  "write.table")

R/DataClasses.R

@@ -441,6 +441,14 @@ setClass("XProcessHistory",
 #' 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.
+#' 
+#' @param verboseBetaColumns Option to calculate two additional metrics of peak
+#' quality via comparison to an idealized bell curve. Adds \code{beta_cor} and
+#' \code{beta_snr} to the \code{chromPeaks} output, corresponding to a Pearson
+#' correlation coefficient to a bell curve with several degrees of skew as well
+#' as an estimate of signal-to-noise using the residuals from the best-fitting
+#' bell curve. See https://github.com/sneumann/xcms/pull/685 and 
+#' https://doi.org/10.1186/s12859-023-05533-4 for more information.
 #'
 #' @details
 #'
@@ -491,7 +499,7 @@ NULL
 #' for a chromatographic peak detection using the centWave method. Instances
 #' should be created with the \code{CentWaveParam} constructor.
 #'
-#' @slot ppm,peakwidth,snthresh,prefilter,mzCenterFun,integrate,mzdiff,fitgauss,noise,verboseColumns,roiList,firstBaselineCheck,roiScales,extendLengthMSW See corresponding parameter above. Slots values should exclusively be accessed
+#' @slot ppm,peakwidth,snthresh,prefilter,mzCenterFun,integrate,mzdiff,fitgauss,noise,verboseColumns,roiList,firstBaselineCheck,roiScales,extendLengthMSW,verboseBetaColumns See corresponding parameter above. Slots values should exclusively be accessed
 #' \emph{via} the corresponding getter and setter methods listed above.
 #'
 #' @rdname findChromPeaks-centWave
@@ -533,7 +541,8 @@ setClass("CentWaveParam",
  roiList = "list",
  firstBaselineCheck = "logical",
  roiScales = "numeric",
- extendLengthMSW = "logical"
+ extendLengthMSW = "logical",
+ verboseBetaColumns = "logical"
  ),
  contains = c("Param"),
  prototype = prototype(
@@ -550,7 +559,8 @@ setClass("CentWaveParam",
  roiList = list(),
  firstBaselineCheck = TRUE,
  roiScales = numeric(),
- extendLengthMSW = FALSE
+ extendLengthMSW = FALSE,
+ verboseBetaColumns = FALSE
  ),
  validity = function(object) {
  msg <- character()
@@ -1242,7 +1252,7 @@ NULL
 #' \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,snthreshIsoROIs,maxCharge,maxIso,mzIntervalExtension,polarity
+#' @slot ppm,peakwidth,snthresh,prefilter,mzCenterFun,integrate,mzdiff,fitgauss,noise,verboseColumns,roiList,firstBaselineCheck,roiScales,extendLengthMSW,verboseBetaColumns,snthreshIsoROIs,maxCharge,maxIso,mzIntervalExtension,polarity
 #' See corresponding parameter above.
 #'
 #' @rdname findChromPeaks-centWaveWithPredIsoROIs

R/do_findChromPeaks-functions.R

@@ -104,6 +104,12 @@
 #' \item{scmin}{Left peak limit found by wavelet analysis (scan number).}
 #' \item{scmax}{Right peak limit found by wavelet analysis (scan numer).}
 #' }
+#' Additional columns for \code{verboseBetaColumns = TRUE}:
+#' \describe{
+#'
+#' \item{beta_cor}{Correlation between an "ideal" bell curve and the raw data}
+#' \item{beta_snr}{Signal-to-noise residuals calculated from the beta_cor fit}
+#' }
 #'
 #' @author Ralf Tautenhahn, Johannes Rainer
 #'
@@ -145,7 +151,8 @@ do_findChromPeaks_centWave <- function(mz, int, scantime, valsPerSpect,
  firstBaselineCheck = TRUE,
  roiScales = NULL,
  sleep = 0,
- extendLengthMSW = FALSE) {
+ extendLengthMSW = FALSE,
+ verboseBetaColumns = FALSE) {
  if (getOption("originalCentWave", default = TRUE)) {
  ## message("DEBUG: using original centWave.")
  .centWave_orig(mz = mz, int = int, scantime = scantime,
@@ -156,7 +163,8 @@ do_findChromPeaks_centWave <- function(mz, int, scantime, valsPerSpect,
  verboseColumns = verboseColumns, roiList = roiList,
  firstBaselineCheck = firstBaselineCheck,
  roiScales = roiScales, sleep = sleep,
- extendLengthMSW = extendLengthMSW)
+ extendLengthMSW = extendLengthMSW, 
+ verboseBetaColumns = verboseBetaColumns)
  } else {
  ## message("DEBUG: using modified centWave.")
  .centWave_new(mz = mz, int = int, scantime = scantime,
@@ -178,7 +186,7 @@ do_findChromPeaks_centWave <- function(mz, int, scantime, valsPerSpect,
  noise = 0, ## noise.local=TRUE,
  sleep = 0, verboseColumns = FALSE, roiList = list(),
  firstBaselineCheck = TRUE, roiScales = NULL,
- extendLengthMSW = FALSE) {
+ extendLengthMSW = FALSE, verboseBetaColumns = FALSE) {
  ## Input argument checking.
  if (missing(mz) | missing(int) | missing(scantime) | missing(valsPerSpect))
  stop("Arguments 'mz', 'int', 'scantime' and 'valsPerSpect'",
@@ -218,6 +226,7 @@ do_findChromPeaks_centWave <- function(mz, int, scantime, valsPerSpect,
  "into", "intb", "maxo", "sn")
  verbosenames <- c("egauss", "mu", "sigma", "h", "f", "dppm", "scale",
  "scpos", "scmin", "scmax", "lmin", "lmax")
+ betanames <- c("beta_cor", "beta_snr")
 
  ## Peak width: seconds to scales
  scalerange <- round((peakwidth / mean(diff(scantime))) / 2)
@@ -226,15 +235,19 @@ do_findChromPeaks_centWave <- function(mz, int, scantime, valsPerSpect,
  scalerange <- scalerange[-z]
  if (length(scalerange) < 1) {
  warning("No scales? Please check peak width!")
- if (verboseColumns) {
- nopeaks <- matrix(nrow = 0, ncol = length(basenames) +
- length(verbosenames))
- colnames(nopeaks) <- c(basenames, verbosenames)
- } else {
- nopeaks <- matrix(nrow = 0, ncol = length(basenames))
- colnames(nopeaks) <- c(basenames)
- }
- return(invisible(nopeaks))
+ matrix_length <- length(basenames)
+ matrix_names <- basenames
+ if (verboseColumns) {
+ matrix_length <- matrix_length + length(verbosenames)
+ matrix_names <- c(matrix_names, verbosenames)
+ }
+ if (verboseBetaColumns) {
+ matrix_length <- matrix_length + length(betanames)
+ matrix_names <- c(matrix_names, betanames)
+ }
+ nopeaks <- matrix(nrow = 0, ncol = matrix_length)
+ colnames(nopeaks) <- matrix_names
+ return(invisible(nopeaks))
  }
 
  if (length(scalerange) > 1)
@@ -319,15 +332,19 @@ do_findChromPeaks_centWave <- function(mz, int, scantime, valsPerSpect,
  ## ROI.list <- findmzROI(object,scanrange=scanrange,dev=ppm * 1e-6,minCentroids=minCentroids, prefilter=prefilter, noise=noise)
  if (length(roiList) == 0) {
  warning("No ROIs found! \n")
- if (verboseColumns) {
- nopeaks <- matrix(nrow = 0, ncol = length(basenames) +
- length(verbosenames))
- colnames(nopeaks) <- c(basenames, verbosenames)
- } else {
- nopeaks <- matrix(nrow = 0, ncol = length(basenames))
- colnames(nopeaks) <- c(basenames)
- }
- return(invisible(nopeaks))
+ matrix_length <- length(basenames)
+ matrix_names <- basenames
+ if (verboseColumns) {
+ matrix_length <- matrix_length + length(verbosenames)
+ matrix_names <- c(matrix_names, verbosenames)
+ }
+ if (verboseBetaColumns) {
+ matrix_length <- matrix_length + length(betanames)
+ matrix_names <- c(matrix_names, betanames)
+ }
+ nopeaks <- matrix(nrow = 0, ncol = matrix_length)
+ colnames(nopeaks) <- matrix_names
+ return(invisible(nopeaks))
  }
  }
 
@@ -525,7 +542,8 @@ do_findChromPeaks_centWave <- function(mz, int, scantime, valsPerSpect,
  td[best.scale.pos],
  td[lwpos],
  td[rwpos], ## Peak positions guessed from the wavelet's (scan nr)
- NA, NA)) ## Peak limits (scan nr)
+ NA, NA, ## Peak limits (scan nr)
+ NA, NA)) ## Beta fitting values
  peakinfo <- rbind(peakinfo,
  c(best.scale, best.scale.nr,
  best.scale.pos, lwpos, rwpos))
@@ -538,7 +556,7 @@ do_findChromPeaks_centWave <- function(mz, int, scantime, valsPerSpect,
 
  ## postprocessing
  if (!is.null(peaks)) {
- colnames(peaks) <- c(basenames, verbosenames)
+ colnames(peaks) <- c(basenames, verbosenames, betanames)
  colnames(peakinfo) <- c("scale", "scaleNr", "scpos",
  "scmin", "scmax")
  for (p in 1:dim(peaks)[1]) {
@@ -561,6 +579,16 @@ do_findChromPeaks_centWave <- function(mz, int, scantime, valsPerSpect,
  lm <- .narrow_rt_boundaries(lm, d)
  lm_seq <- lm[1]:lm[2]
  pd <- d[lm_seq]
+
+ # Implement a fit of a skewed gaussian (beta distribution)
+ # for peak shape and within-peak signal-to-noise ratio
+ # See https://doi.org/10.1186/s12859-023-05533-4 and
+ # https://github.com/sneumann/xcms/pull/685
+ if(verboseBetaColumns){
+ beta_vals <- .get_beta_values(pd)
+ peaks[p, "beta_cor"] <- beta_vals$best_cor
+ peaks[p, "beta_snr"] <- beta_vals$beta_snr
+ }
 
  peakrange <- td[lm]
  peaks[p, "rtmin"] <- scantime[peakrange[1]]
@@ -675,22 +703,30 @@ do_findChromPeaks_centWave <- function(mz, int, scantime, valsPerSpect,
 
  if (length(peaklist) == 0) {
  warning("No peaks found!")
-
- if (verboseColumns) {
- nopeaks <- matrix(nrow = 0, ncol = length(basenames) +
- length(verbosenames))
- colnames(nopeaks) <- c(basenames, verbosenames)
- } else {
- nopeaks <- matrix(nrow = 0, ncol = length(basenames))
- colnames(nopeaks) <- c(basenames)
- }
- message(" FAIL: none found!")
- return(nopeaks)
+ matrix_length <- length(basenames)
+ matrix_names <- basenames
+ if (verboseColumns) {
+ matrix_length <- matrix_length + length(verbosenames)
+ matrix_names <- c(matrix_names, verbosenames)
+ }
+ if (verboseBetaColumns) {
+ matrix_length <- matrix_length + length(betanames)
+ matrix_names <- c(matrix_names, betanames)
+ }
+ nopeaks <- matrix(nrow = 0, ncol = matrix_length)
+ colnames(nopeaks) <- matrix_names
+ message(" FAIL: none found!")
+ return(nopeaks)
  }
  p <- do.call(rbind, peaklist)
- if (!verboseColumns)
- p <- p[, basenames, drop = FALSE]
-
+ keepcols <- basenames
+ if (verboseColumns){
+ keepcols <- c(keepcols, verbosenames)
+ }
+ if(verboseBetaColumns){
+ keepcols <- c(keepcols, betanames)
+ }
+ p <- p[, keepcols, drop = FALSE]
  uorder <- order(p[, "into"], decreasing = TRUE)
  pm <- as.matrix(p[,c("mzmin", "mzmax", "rtmin", "rtmax"), drop = FALSE])
  uindex <- rectUnique(pm, uorder, mzdiff, ydiff = -0.00001) ## allow adjacent peaks
@@ -2623,6 +2659,11 @@ do_findKalmanROI <- function(mz, int, scantime, valsPerSpect,
 #' \item{scmin}{Left peak limit found by wavelet analysis (scan number).}
 #' \item{scmax}{Right peak limit found by wavelet analysis (scan numer).}
 #' }
+#' Additional columns for \code{verboseBetaColumns = TRUE}:
+#' \describe{
+#' \item{beta_cor}{Correlation between an "ideal" bell curve and the raw data}
+#' \item{beta_snr}{Signal-to-noise residuals calculated from the beta_cor fit}
+#' }
 #'
 #' @rdname do_findChromPeaks_centWaveWithPredIsoROIs
 #'
@@ -2634,7 +2675,8 @@ do_findChromPeaks_centWaveWithPredIsoROIs <-
  verboseColumns = FALSE, roiList = list(),
  firstBaselineCheck = TRUE, roiScales = NULL, snthreshIsoROIs = 6.25,
  maxCharge = 3, maxIso = 5, mzIntervalExtension = TRUE,
- polarity = "unknown", extendLengthMSW = FALSE) {
+ polarity = "unknown", extendLengthMSW = FALSE, 
+ verboseBetaColumns = FALSE) {
  ## Input argument checking: most of it will be done in
  ## do_findChromPeaks_centWave
  polarity <- match.arg(polarity, c("positive", "negative", "unknown"))
@@ -2655,7 +2697,8 @@ do_findChromPeaks_centWaveWithPredIsoROIs <-
  roiList = roiList,
  firstBaselineCheck = firstBaselineCheck,
  roiScales = roiScales,
- extendLengthMSW = extendLengthMSW)
+ extendLengthMSW = extendLengthMSW,
+ verboseBetaColumns = verboseBetaColumns)
  return(do_findChromPeaks_addPredIsoROIs(mz = mz, int = int,
  scantime = scantime,
  valsPerSpect = valsPerSpect,
@@ -3668,3 +3711,49 @@ peaksWithCentWave <- function(int, rt,
  }
  lm
 }
+
+
+#' @description
+#'
+#' Calculate beta parameters for a chromatographic peak, both its similarity
+#' to a bell curve of varying degrees of skew and the standard deviation of the
+#' residuals after the best-fit bell is normalized and subtracted.
+#'
+#' @param eic_ints A numeric vector corresponding to the peak intensities
+#' @param skews A numeric vector of the skews to try, corresponding to the
+#' shape1 of dbeta with a shape2 of 5. Values less than 5 will be increasingly
+#' right-skewed, while values greater than 5 will be left-skewed.
+#'
+#' @author William Kumler
+#'
+#' @noRd
+.get_beta_values <- function(eic_ints, skews=c(3, 3.5, 4, 4.5, 5)){
+ if(length(eic_ints)<5){
+ best_cor <- NA
+ beta_snr <- NA
+ } else {
+ beta_sequence <- rep(seq(0, 1, length.out=length(eic_ints)), each=5)
+ beta_vals <- t(matrix(dbeta(beta_sequence, shape1 = skews, shape2 = 5), nrow = 5))
+ # matplot(beta_vals)
+ beta_cors <- cor(eic_ints, beta_vals)
+ best_cor <- max(beta_cors)
+ best_curve <- beta_vals[,which.max(beta_cors)]
+ noise_level <- sd(diff(.scale_zero_one(best_curve)-.scale_zero_one(eic_ints)))
+ beta_snr <- log10(max(eic_ints)/noise_level)
+ }
+ list(best_cor=best_cor, beta_snr=beta_snr)
+}
+
+#' @description
+#'
+#' Simple helper function to scale a numeric vector between zero and one by
+#' subtracting the lowest value and dividing by the maximum.
+#'
+#' @param num_vec `numeric` vector, typically chromatographic intensities.
+#'
+#' @author William Kumler
+#'
+#' @noRd
+.scale_zero_one <- function(num_vec){
+ (num_vec-min(num_vec))/(max(num_vec)-min(num_vec))
+}

R/functions-Params.R

@@ -86,14 +86,15 @@ CentWaveParam <- function(ppm = 25, peakwidth = c(20, 50), snthresh = 10,
  integrate = 1L, mzdiff = -0.001, fitgauss = FALSE,
  noise = 0, verboseColumns = FALSE, roiList = list(),
  firstBaselineCheck = TRUE, roiScales = numeric(),
- extendLengthMSW = FALSE) {
+ extendLengthMSW = FALSE, verboseBetaColumns = FALSE) {
  return(new("CentWaveParam", 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))
+ extendLengthMSW = extendLengthMSW, 
+ verboseBetaColumns=verboseBetaColumns))
 }
 
 #' @return The \code{MatchedFilterParam} function returns a
@@ -213,6 +214,7 @@ CentWavePredIsoParam <- function(ppm = 25, peakwidth = c(20, 50), snthresh = 10,
  integrate = 1L, mzdiff = -0.001, fitgauss = FALSE,
  noise = 0, verboseColumns = FALSE, roiList = list(),
  firstBaselineCheck = TRUE, roiScales = numeric(),
+ extendLengthMSW = FALSE, verboseBetaColumns = FALSE,
  snthreshIsoROIs = 6.25, maxCharge = 3, maxIso = 5,
  mzIntervalExtension = TRUE, polarity = "unknown") {
  return(new("CentWavePredIsoParam", ppm = ppm, peakwidth = peakwidth,
@@ -221,6 +223,8 @@ CentWavePredIsoParam <- function(ppm = 25, peakwidth = c(20, 50), snthresh = 10,
  mzdiff = mzdiff, fitgauss = fitgauss, noise = noise,
  verboseColumns = verboseColumns, roiList = roiList,
  firstBaselineCheck = firstBaselineCheck, roiScales = roiScales,
+ extendLengthMSW = extendLengthMSW, 
+ verboseBetaColumns = verboseBetaColumns,
  snthreshIsoROIs = snthreshIsoROIs, maxIso = as.integer(maxIso),
  maxCharge = as.integer(maxCharge),
  mzIntervalExtension = mzIntervalExtension, polarity = polarity))

R/functions-XCMSnExp.R

@@ -295,6 +295,12 @@ dropGenericProcessHistory <- function(x, fun) {
  meanMz <- do.call(mzCenterFun, list(mtx[, 2], mtx[, 3]))
  if (is.na(meanMz)) meanMz <- mtx[maxi[1], 2]
  res[i, "mz"] <- meanMz
+
+ if("beta_cor"%in%cn){
+ beta_vals <- .get_beta_values(mtx[,3])
+ res[i, "beta_cor"] <- beta_vals$best_cor
+ res[i, "beta_snr"] <- beta_vals$beta_snr
+ }
  } else {
  res[i, ] <- rep(NA_real_, ncols)
  }