Mass decomposition with the Rdisop package

Chemical background

Atoms are composed of electrons, protons, and neutrons. The number of protons (the atomic number) is fixed and defines what element the atom is. The number of neutrons, on the other hand, can vary: Atoms with the same number of protons but different numbers of neutrons are called isotopes of the element. Each of these isotopes occurs in nature with a certain abundance. The nominal mass of a molecule is the sum of protons and neutrons of the constituting atoms. The mass of the molecule is the sum of masses of these atoms. The monoisotopic (nominal) mass of a molecule is the sum of (nominal) masses of the constituting atoms where for every element its most abundant natural isotope is chosen. Clearly, nominal mass and mass depend on the isotopes the molecule consists of, thus on the isotope species of the molecule.

No present-day analysis technique is capable of resolving isotope species with identical nominal mass. Instead, these isotope species appear as one single peak in the mass spectrometry output. For this reason, we merge isotope species with identical nominal mass and refer to the resulting distribution as the molecule’s isotope pattern.

Identification schema

Obtaining an accurate isotope pattern from a high resolution mass spectrometer, we apply this information to identify the elemental composition of the sample molecule. Our input is a list of masses with normalized abundances that corresponds to the isotope pattern of the sample molecule. We want to find that molecule’s elemental composition whose isotope pattern best matches the input.

Solving this task is divided into the following parts: First, all elemental compositions are calculated that share some property, for example monoisotopic mass, with the input spectrum. Second, to remove those compositions that do not exist in nature, chemical bonding rules are applied, discarding formulas that have negative or non-integer degree of unsaturation. And third, for every remaining composition, its theoretical isotope pattern is calculated and compared to the measured isotope pattern. Candidate patterns are ranked using Bayesian Statistics, and the one with the highest score is chosen.

Functions decomposeMass and decomposeIsotopes

The function decomposeMass returns a list of molecules which have a given exact mass (within an error window in ppm):

molecules <- decomposeMass(46.042, ppm=20, maxisotopes = 4)
molecules
#> $formula
#> [1] "C2H6O"
#> 
#> $score
#> [1] 1
#> 
#> $exactmass
#> [1] 46.04186
#> 
#> $charge
#> [1] 0
#> 
#> $parity
#> [1] "e"
#> 
#> $valid
#> [1] "Valid"
#> 
#> $DBE
#> [1] 0
#> 
#> $isotopes
#> $isotopes[[1]]
#>            [,1]        [,2]        [,3]         [,4]
#> [1,] 46.0418648 47.04534542 48.04631711 4.904960e+01
#> [2,]  0.9749152  0.02293559  0.00210353 4.540559e-05

This call produces a list of potential molecules (with a single element in this case). The larger the masses, the allowed ppm deviation and the allowed elements list, the larger the result list will grow. For each hypothesis there is its formula and weight and score. The parity, validity (using the nitrogen rule) and double bond equivalents (DBE) are simple, yet commonly used hints for the plausibility of a solution and can be used for filtering the results list. For an amino acid this simple method guesses already eight hypotheses:

length(decomposeMass(147.053))
#> [1] 8

On modern mass spectrometers a full isotope pattern can be obtained for a molecule, and the masses and intensities improve the accuracy of the sum formula prediction. Accessor functions return only subsets of the molecule data structure:

# glutamic acid (C5H9NO4)
masses <- c(147.053, 148.056)
intensities <- c(93, 5.8)

molecules <- decomposeIsotopes(masses, intensities)
data.frame(getFormula(molecules), getScore(molecules), getValid(molecules))
#>   getFormula.molecules. getScore.molecules. getValid.molecules.
#> 1               C5H9NO4        1.000000e+00               Valid
#> 2              C3H17P2S        1.935336e-09             Invalid

The first ranked solution already has a score close to one, and if using an N-rule filter, only one solution would remain. These cases are not removed by default, because a few compound classes do not obey the N-rule, which after all is just a simple heuristic.

If the masses were obtained by an LC-ESI-MS, it is likely that the measured mass signal actually resembles an adduct ion, such as [M+H]⁺. The sum formula obtained through decomposeIsotopes will have one H too much, and will not be found in PubChem or other libraries, unless the adduct has been removed:

getFormula(subMolecules("C5H10NO4", "H"))
#> [1] "C5H9NO4"

Similarly, if during ionization an in-source fragmentation occurred, the lost fragment can be added before querying using addMolecules.

Interaction with other BioConductor packages

This section will give some suggestions how the Rdisop functionality can be combined with other BioConductor packages.

Usually the masses and intensities will be obtained from a high-resolution mass spectrometer such as an FT-ICR-MS or QTOF-MS. BioConductor currently has two packages dealing with peak picking on raw machine data, MassSpecWavelet and XCMS. The latter contains a wrapper for MassSpecWavelet, so we need to deal with XCMS peak lists only. The ESI package¹ can extract a set of isotope clusters from peak lists.

After Rdisop has created a set of candidate molecular formulae, the open-access compound databases PubChem or ChEBI can be queried whether any information about this compound exists. Nota bene: a hit or non-hit does not indicate a correct or incorrect formula, but merely helps in further verification or structure elucidation steps. For other cheminformatics functionality in BioConductor see e.g. RCDK.