Fix for issue #22

ChangeLog

@@ -1,3 +1,12 @@
+2025-01-22 Jan Lisec <[email protected]>
+
+    * Rdisop-1.67.5
+    * fix for issue #22; implementing atomic masses with higher precision
+    * the new set of elemental definitions will be used as the default which can
+      be considered a breaking change as it might lead to minor changes in 
+      previous analyses; previous results can be obtained by setting `method = 'IUPAC'`
+      in function initialize
+
 2024-11-08 Jan Lisec <[email protected]>
 
     * Rdisop-1.67.4

DESCRIPTION

@@ -8,7 +8,7 @@ Authors@R:
         person("Steffen", "Neumann", , "[email protected]", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-7899-7192")),
         person("Jan", "Lisec", , "[email protected]", role = c("ctb"), comment = c(ORCID = "0000-0003-1220-2286")),
         person("Miao", "Yu", , "[email protected]", role = c("ctb")),
-	person("Roberto", "Canteri", , , role = c("ctb"))
+	    person("Roberto", "Canteri", , , role = c("ctb"))
     )
 Description: In high resolution mass spectrometry (HR-MS), the measured masses
     can be decomposed into potential element combinations (chemical sum formulas). 

R/initializeElements.R

@@ -11,8 +11,8 @@
 #'     reduce the number of decomposition hypotheses, subsets of elements can be 
 #'     created.
 #'     
-#' @param method Use isotope mass and abundance data from either "IUPAC"
-#(default) or "NIST".
+#' @param method Use isotope mass and abundance data from either "NIST"
+#'     (default) or "IUPAC".
 #'
 #' @param names Vector of element names within PSE.
 #'
@@ -32,8 +32,9 @@
 #' @author Steffen Neumann <[email protected]>
 #' @references For a description of the underlying IMS see citation("Rdisop").
 #'     Isotope patterns were obtained through wikipedia.org
-initializeElements <- function(names) {
-    elements <- initializePSE()
+initializeElements <- function(names, method=c("NIST","IUPAC")) {
+    method <- match.arg(method)
+    elements <- initializePSE(method = method)
     lapply(names, function (name) {.getElement(name, elements)})
 }
 
@@ -72,7 +73,8 @@ initializeCHNOPSNaK <- function() {
 
 #' @rdname initializeElements
 #' @export
-initializePSE <- function(method="IUPAC") {
+initializePSE <- function(method=c("NIST","IUPAC")) {
+    method <- match.arg(method)
     if (method=="IUPAC")
         return(.initializePSE_IUPAC())
     else if (method=="NIST")
@@ -81,8 +83,8 @@ initializePSE <- function(method="IUPAC") {
         stop("Unknown table requested")
 }
 
-#' @rdname initializeElements
-#' @export
+#' @noRd
+#' @keywords internal
 .initializePSE_IUPAC <- function() {	
     D <- list(name="D", mass=2,  isotope = list(mass=c(0.014102), abundance=c(1))) #Heavy Water
     Ac <- list(name= 'Ac', mass=227, isotope=list(mass=c(0.02775),abundance=c(1)))
@@ -192,10 +194,10 @@ initializePSE <- function(method="IUPAC") {
     list(D, Ac, Ag, Al, Am, Ar, As, At, Au, B, Ba, Be, Bi, Bk, Br, C, Ca, Cd, Ce, Cf, Cl, Cm, Co, Cr, Cs, Cu, Dy, Er, Es, Eu, F, Fe, Fm, Fr, Ga, Gd, Ge, H, He, Hf, Hg, Ho, I, In, Ir, K, Kr, La, Li, Lr, Lu, Md, Mg, Mn, Mo, N, Na, Nb, Nd, Ne, Ni, No, Np, O, Os, P, Pa, Pb, Pd, Pm, Po, Pr, Pt, Pu, Ra, Rb, Re, Rh, Rn, Ru, S, Sb, Sc, Se, Si, Sm, Sn, Sr, Ta, Tb, Tc, Te, Th, Ti, Tl, Tm, U, V, W, Xe, Y, Yb, Zn, Zr)
 }
 
-#' @rdname initializeElements_NIST
-#' @export
+#' @noRd
+#' @keywords internal
 .initializePSE_NIST <- function() {
-    D <- list(name="D", mass=2,  isotope = list(mass=c(0.014102), abundance=c(1))) #Heavy Water
+    D <- list(name="D", mass=2,  isotope = list(mass=c(0.01410178), abundance=c(1))) #Heavy Water
     Ac <- list(name='Ac', mass=227, isotope=list(mass=c(0.02775230), abundance=c(1.00000000)))
     Ag <- list(name='Ag', mass=107, isotope=list(mass=c(-0.09490840, 0.00000000, -0.09524470), abundance=c(0.51839000, 0.00000000, 0.48161000)))
     Al <- list(name='Al', mass=27, isotope=list(mass=c(-0.01846147), abundance=c(1.00000000)))
@@ -263,7 +265,7 @@ initializePSE <- function(method="IUPAC") {
     P <- list(name='P', mass=31, isotope=list(mass=c(-0.02623800), abundance=c(1.00000000)))
     Pa <- list(name='Pa', mass=231, isotope=list(mass=c(0.03588420), abundance=c(1.00000000)))
     Pb <- list(name='Pb', mass=204, isotope=list(mass=c(-0.02695600, 0.00000000, -0.02553430, -0.02410270, -0.02334750), abundance=c(0.01400000, 0.00000000, 0.24100000, 0.22100000, 0.52400000)))
-    Pd <- list(name='Pd', mass=102, isotope=list(mass=c(-0.09439780, 0.00000000, -0.09596950, -0.09492040, -0.09651960, 0.00000000, -0.09610840, 0.00000000, -0.09484700), abundance=c(0.01020000, 0.00000000, 0.11140000, 0.22330000, 0.27330000, 0.00000000, 0.26460000, 0.00000000, 0.11720000)))
+    Pd <- list(name='Pd', mass=102, isotope=list(mass=c(-0.09439780, 0.00000000, -0.09596950, -0.09492040, -0.09651960, 0.00000000, -0.09610840, 0.00000000, -0.09482780), abundance=c(0.01020000, 0.00000000, 0.11140000, 0.22330000, 0.27330000, 0.00000000, 0.26460000, 0.00000000, 0.11720000)))
     Pm <- list(name='Pm', mass=145, isotope=list(mass=c(-0.08724410), abundance=c(1.00000000)))
     Po <- list(name='Po', mass=209, isotope=list(mass=c(-0.01756920), abundance=c(1.00000000)))
     Pr <- list(name='Pr', mass=141, isotope=list(mass=c(-0.09234240), abundance=c(1.00000000)))

data-raw/isotopes.R

@@ -1,7 +1,7 @@
 # to be more comparable with package `enviPat` and to allow easier access to the
 # isotope information the original functions defining the PSE as a list is used
 # to set up a data.frame which can be attended by Rdisop::isotopes
-ele <- initializePSE()
+ele <- initializePSE(method = "NIST")
 isotopes <- plyr::ldply(ele, function(x) {
     data.frame(
         "element" = rep(x$name, length(x$isotope$mass)),

data-raw/mono_masses.R

@@ -1,7 +1,7 @@
 # to provide a fast way calculating the mono-isotopic mass of an element, here a
 # named vector containing the main isotopes of all PSE elements is prepared to 
 # be accessed by Rdisop::mono_masses
-ele <- initializePSE()
+ele <- initializePSE(method = "NIST")
 mono_masses <- sort(sapply(ele, function(x) {
     m <- x$mass + x$isotope$mass + (1:length(x$isotope$mass))-1
     m <- m[which.max(x$isotope$abundance)]

tests/testthat/test-getMolecule.R

@@ -23,7 +23,9 @@ testthat::test_that(
         out <- getMolecule(fml)
 
         # the exact mass is the second isotope
-        testthat::expect_equal(getMass(out), getIsotope(out, 2)[[1]][1,])
+        # this test is no longer valid since v.1.67.5 due to defaulting to more precise NIST elemental definitions
+        # before v1.67.5 the rounding errors led to the 'exactmass' being the second isotope
+        # testthat::expect_equal(getMass(out), getIsotope(out, 2)[[1]][1,])
 
         # the monoisotopic mass is the first isotope
         testthat::expect_equal(unname(getMonoisotopic(out)), getIsotope(out, 1)[[1]][1,])
@@ -48,5 +50,24 @@ testthat::test_that(
         # elements contained in the formula need to exist in the provided elements list
         testthat::expect_error(getMolecule(formula = "CH4", elements = initializeElements("H")))
 
+    }
+)
+
+testthat::test_that(
+    desc = "getMolecule improves precision with NIST elemental definitions compared to IUPAC", 
+    code = {
+
+        # this was an error reported in issue #22
+        # example fml returns a wrong monoisotopic peak due to rounding differences
+        fml <- "C89H166O17P2"
+        res_iupac <- Rdisop::getMolecule(formula = fml, elements = initializePSE(method = "IUPAC"), maxisotopes = 2)$isotopes[[1]]
+        res_nist <- Rdisop::getMolecule(formula = fml, elements = initializePSE(method = "NIST"), maxisotopes = 2)$isotopes[[1]]
+
+        testthat::expect_equal(res_iupac, structure(c(1569.16002852, 0.494738964752399, 1570.16346006151, 0.505261035247601), dim = c(2L, 2L)))
+        testthat::expect_equal(res_nist, structure(c(1569.16002752, 0.502975666994795, 1570.16344446454, 0.497024333005205), dim = c(2L, 2L)))
+
+        # the monoisotopic peak has changed
+        testthat::expect_true(res_iupac[2,1]<0.5 & res_nist[2,1]>0.5)
+
     }
 )