feat: cubic calibration support; fixes #3

DESCRIPTION

@@ -1,11 +1,21 @@
 Package: readBrukerFlexData
-Version: 1.8.5
-Date: 2017-04-22
+Version: 1.8.5.99
+Date: 2022-04-15
 Title: Reads Mass Spectrometry Data in Bruker *flex Format
-Authors@R: person("Sebastian", "Gibb", role=c("aut", "cre"),
-    email="mail@sebastiangibb.de")
+Authors@R:
+    c(
+        person(
+            "Sebastian", "Gibb", role = c("aut", "cre"),
+            email = "mail@sebastiangibb.de",
+            comment = c(ORCID = "0000-0001-7406-4443")
+        ),
+        person("Samuel", "Granjeaud", role = "ctb"),
+        person("Alan", "Race", role = "ctb")
+    )
 Depends:
     R (>= 3.3.0)
+Suggests:
+    testthat
 Description: Reads data files acquired by Bruker Daltonics' matrix-assisted
     laser desorption/ionization-time-of-flight mass spectrometer of the *flex
     series.
@@ -14,4 +24,4 @@ URL: http://strimmerlab.org/software/maldiquant/
     https://github.com/sgibb/readBrukerFlexData/
 BugReports: https://github.com/sgibb/readBrukerFlexData/issues/
 LazyLoad: yes
-RoxygenNote: 6.0.1
+RoxygenNote: 7.1.2

NEWS

@@ -1,6 +1,12 @@
 RELEASE HISTORY OF THE "readBrukerFlexData" PACKAGE
 ===================================================
 
+Version 1.8.5.99 [2022-04-15]:
+- Add basic support for cubic calibration as encoded in
+  "V1.0TOF2CalibrationConstants" (fixes
+  https://github.com/sgibb/readBrukerFlexData/issues/3).
+  Contributed by Samuel Granjeaud (@SamGG) and Alan Race (@AlanRace).
+
 Version 1.8.5 [2017-04-22]:
 - Don't errorously convert tof to mass if "V1.0CTOF2CalibrationConstants" was
   used and throw a warning.

R/ntbcal-functions.R

@@ -0,0 +1,78 @@
+#' Extracts NTBCal coefficients.
+#'
+#' This is a helper function to extract coefficients and constants from
+#' metaData$NtbCalStr
+#'
+#' @param x \code{character}, NTBCal line
+#'
+#' \preformatted{
+#' "##$NTBCal= <V3.0CCalibrator 12 1 0 0 -2147483648 2147483647  V1.0CHPCData  Order 0 vCoeff V1.0VectorDouble 0  c2 0 c0 0 minMass 0 maxMass 0 bUse 0 endCHPCData V1.0CTOF2CalibrationConstants 48 0.25 3884.868239229128 1164156.9998801197 6.1682474689786355 -0.063083671657452864 -34.362308485323794 2 V1.0CTOF2CalibrationConstants 0 0 0 0 0 0 0 -1  >"
+#' }
+#'
+#' @return a \code{double}, coefficients
+#'
+#' @seealso \code{\link[readBrukerFlexData]{.ctof2calibration}}
+#' @noRd
+#' @keywords internal
+#' @return \code{double}, vector of the calibration constants.
+#'
+#' @keywords internal
+#' @noRd
+.extractV10CTOF2CalibrationConstants <- function(x) {
+    x <- gsub("^.*V1.0CTOF2CalibrationConstants ([0-9. -]*) V1.0CTOF2CalibrationConstants.*$", "\\1", x)
+    as.numeric(unlist(strsplit(x, " ", fixed = TRUE)))
+}
+
+#' Cubic calibration
+#'
+#' Only basic support (not 100\% identical results) for Bruker Daltonics' cubic
+#' calibration (encoded in "##\$NTBCal" as "V1.0CTOF2CalibrationConstants").
+#' This is an internal function and should normally not used by the user.
+#'
+#' @param tof \code{double}, time of flight.
+#' @param d \code{double}, calibration constants, return value of
+#' \code{.extractV10CTOF2CalibrationConstants}.
+#' @return \code{double}, mz values.
+#'
+#' @note
+#' Bruker Daltonics doesn't explain how the cubic calibration works.
+#' All formula are results of \dQuote{trial and error}.
+#' That is why mass calculated by \code{.ctof2calibration}
+#' differs little from original mass.
+#' See also \url{https://github.com/sgibb/readBrukerFlexData/issues/3}.
+#'
+#' @seealso
+#' \url{https://github.com/sgibb/readBrukerFlexData/issues/3}
+#'
+#' @author
+#' Alan Race, Samuel Granjeaud, Sebastian Gibb
+#'
+#' @keywords internal
+#' @rdname ctof2calibration
+.ctof2calibration <- function(tof, d) {
+    ## e.g.  "V1.0CTOF2CalibrationConstants 48 0.25 3884.868239229128 1164156.9998801197 6.1682474689786355 -0.063083671657452864 -34.362308485323794 2 V1.0CTOF2CalibrationConstants"
+    ## DELAY DW ML2 ML1 ML3 A E 2
+    mz <- .tof2mass(tof, d[4], d[3], d[5])
+
+    A <- d[6]
+    B <- d[5]
+    C <- sqrt(1e12 / d[4])
+    D <- d[3]
+    E <- d[7]
+
+    ## quadratic: 0 = B * (sqrt(m/z))^2 + C * sqrt(m/z) + D(times)
+    ## same formula as tof2mass but instead of D(times) = ML3 - tof
+    ## it seems to be necessary to use D(times) = tof - ML3
+    m <- (-C + sqrt(abs((C * C) - (4 * B * (tof - D))))) / (2 * B)
+
+    s <- sign(tof - D)
+    belowZero <- s < 0
+
+    m[!belowZero] <- vapply(tof[!belowZero], function(x) {
+        ## Assuming cubic with A*x^3 + B*x^2 + C*x + (D - tof)
+        ## where x = sqrt(mz)
+        Re(polyroot(c(D - x, C, B, A))[1L])
+    }, NA_real_)
+
+    m^2 * s - E
+}

R/package.R

@@ -26,20 +26,6 @@
 #' \cr \cr
 #' All trademarks are owned by or licensed to Bruker Daltonics.
 #'
-#' \tabular{ll}{
-#'
-#' Package: \tab readBrukerFlexData \cr
-#'
-#' Version: \tab 1.8.5\cr
-#'
-#' Date: \tab 2017-04-22\cr
-#'
-#' License: \tab GPL (>= 3)\cr
-#'
-#' URL: \tab \url{https://github.com/sgibb/readBrukerFlexData}\cr
-#'
-#' }
-#'
 #' @docType package
 #' @name readBrukerFlexData-package
 #' @author Sebastian Gibb \email{mail@@sebastiangibb.de}

R/readAcquFile-functions.R

@@ -187,9 +187,14 @@
   }
 
   # https://github.com/sgibb/MALDIquantForeign/issues/19
+  # https://github.com/sgibb/readBrukerFlexData/issues/3
   metaData$v1tofCalibration <-
     grepl("V1.0CTOF2CalibrationConstants",
           .grepAcquValue("##\\$NTBCal=", acquLines))
+  metaData$v1tofCalibrationConstants <- 
+    .extractV10CTOF2CalibrationConstants(
+      .grepAcquValue("##\\$NTBCal=", acquLines)
+    )
 
   ## obligate LIFT
   metaData$lift <- c(.grepAcquDoubleValue("##\\$Lift1=", acquLines),

R/readBrukerFlexFile-functions.R

@@ -231,11 +231,24 @@ readBrukerFlexFile <- function(fidFile, removeMetaData=FALSE, useHpc=TRUE,
     tof <- tof[notNull]
   }
 
-  ## calculate mass of TOFs
-  mass <- .tof2mass(tof,
-                    metaData$calibrationConstants[[1L]],
-                    metaData$calibrationConstants[[2L]],
-                    metaData$calibrationConstants[[3L]])
+  ## TODO: fix V1.0CTOF2 calibration
+  ## https://github.com/sgibb/MALDIquantForeign/issues/19
+  ## https://github.com/sgibb/readBrukerFlexData/issues/3
+  if (isTRUE(metaData$v1tofCalibration)) {
+    warning("The spectrum file ", sQuote(fidFile), " uses ",
+            "V1.0CTOF2CalibrationConstants.\n",
+            "V1.0CTOF2CalibrationConstants aren't fully supported by ",
+            "readBrukerFlexFile. See ",
+            "https://github.com/sgibb/readBrukerFlexData/issues/3 ",
+            "for details.")
+    mass <- .ctof2calibration(tof, metaData$v1tofCalibrationConstants)
+  } else {
+    ## calculate mass of TOFs
+    mass <- .tof2mass(tof,
+                      metaData$calibrationConstants[[1L]],
+                      metaData$calibrationConstants[[2L]],
+                      metaData$calibrationConstants[[3L]])
+  }
 
   ## TODO: fix equations in .hpc
 
@@ -254,20 +267,6 @@ readBrukerFlexFile <- function(fidFile, removeMetaData=FALSE, useHpc=TRUE,
                  hpcCoefficients=metaData$hpcCoefficients)
   }
 
-  ## TODO: fix V1.0CTOF2 calibration
-  ## https://github.com/sgibb/MALDIquantForeign/issues/19
-  if (isTRUE(metaData$v1tofCalibration)) {
-    warning("The spectrum file ", sQuote(fidFile), " uses ",
-            "V1.0CTOF2CalibrationConstants.\n",
-            "V1.0CTOF2CalibrationConstants aren't supported by ",
-            "readBrukerFlexFile. See ",
-            "https://github.com/sgibb/MALDIquantForeign/issues/19 ",
-            "for details.\n",
-            "Could not convert time-of-flight values into mass!")
-    mass <- tof
-
-  }
-
   ## TODO: add LIFT support
 
   ## was LIFT involved?

tests/testthat.R

@@ -0,0 +1,2 @@
+library("testthat")
+test_check("readBrukerFlexData")

tests/testthat/test_ntbcal-functions.R

@@ -0,0 +1,5 @@
+test_that(".extractV10CTOF2CalibrationConstants", {
+    x <- "##$NTBCal= <V3.0CCalibrator 12 1 0 0 -2147483648 2147483647  V1.0CHPCData  Order 0 vCoeff V1.0VectorDouble 0  c2 0 c0 0 minMass 0 maxMass 0 bUse 0 endCHPCData V1.0CTOF2CalibrationConstants 48 0.25 3884.868239229128 1164156.9998801197 6.1682474689786355 -0.063083671657452864 -34.362308485323794 2 V1.0CTOF2CalibrationConstants 0 0 0 0 0 0 0 -1  >"
+    r <- c(48, 0.25, 3884.868239229128, 1164156.9998801197, 6.1682474689786355, -0.063083671657452864, -34.362308485323794, 2)
+    expect_identical(.extractV10CTOF2CalibrationConstants(x), r)
+})

tests/testthat/test_readBrukerFlexFile.R

@@ -0,0 +1,5 @@
+test_that("general", {
+    expect_error(readBrukerFlexFile("./"))
+    expect_error(readBrukerFlexFile("test_readBrukerFlexFile.R"))
+    expect_error(readBrukerFlexFile("ape"))
+})

tests/testthat/test_sampleName.R

@@ -0,0 +1,11 @@
+test_that("sampleName", {
+    files <- c(
+        "./Run1/2-100kDa/0_A1/1/1SLin/fid",
+        "foo/Pankreas_HB_L_061019_A10/0_a19/1/1SLin/fid",
+        "LT_1631_05_smear/0_H1_1SLin/fid"
+    )
+    name <- c("2_100kDa", "Pankreas_HB_L_061019_A10", "LT_1631_05_smear")
+
+    for (i in seq(along=files))
+        expect_equal(.sampleName(files[i]), name[i])
+})