Most abundant isotope peaks and efficient selection on $Y=X_1+X_2+\cdots + X_m$

TL;DR

This paper presents a novel, efficient method for computing the most abundant isotope peaks in chemical compounds by reducing the problem to sorting sums of element isotopes, with practical applications demonstrated on complex molecules.

Contribution

The paper introduces a new computational approach that simplifies and accelerates the calculation of isotope distributions in chemistry, linking it to sorting algorithms.

Findings

Efficient algorithm for isotope peak computation

Successful application to large chemical compounds

Improved accuracy in isotope abundance estimation

Abstract

The isotope masses and relative abundances for each element are fundamental chemical knowledge. Computing the isotope masses of a compound and their relative abundances is an important and difficult analytical chemistry problem. We demonstrate that this problem is equivalent to sorting $Y=X_1+X_2+\cdots+X_m$. We introduce a novel, practically efficient method for computing the top values in $Y$. then demonstrate the applicability of this method by computing the most abundant isotope masses (and their abundances) from compounds of nontrivial size.

Equations

$$Y=$$

$$=$$

$$(X_{\frac{n}{2}+1}+X_{\frac{n}{2}+2}+\cdots +X_{\frac{3n}{4}})+(X_{\frac{3n}{4}+1}+X_{\frac{3n}{4}+2}+\cdots +X_n$$

$$=$$

Full text

1 Introduction

Atoms of the same element may have a variable numbers of neutrons in their nulcei. The number of neutrons in the nucleus of an atom influences its mass. The relative abundances with which different isotopes naturally occur is well established.[5]

In compounds composed of several elements, finding the relative abundances of the most prevalent isotopes of the compound (and the respective masses at which these isotopes occur) is a difficult combinatorial problem.

For small problems, this can be solved via brute force: it is possible to compute all isotope masses and their respective abundances, sort them in descending order of abundance, and then retrieving the isotope peaks (we will refer to the mass and relative abundance of a particular isotope as an isotope peak, because that is how they are observed in mass spectrometry) with with greatest abundance; however, the runtime of this brute-force approach is far too inefficient, because it grows $\in\Omega(2^{n})$.

Łącki et al.[5] recently introduced a statistical approach to this problem, by which the top isotope peaks of a compound may be efficiently approximated. For each element, the method generates data from the isotopologue (i.e., from the contributions all possible isotopes from that element). For each element, its possible contributions follow a multinomial distribution, which Łącki et al. approximate using multivarite Gaußians and generate in descending order or probability. The isotopologue contributions are combined over the relevant element sets. At this point, they generate the isotope configurations of the compound in descending order of probability, which is equivalent to finding the $k$ most probable isotopologues of the molecule.

In this paper, we demonstrate that finding the top $k$ isotope peaks of a compound composed of $m$ elements is equivalent to finding the top $k$ values (and respective indices) in $Y=X_{1}+X_{2}+\cdots+X_{m}$, where $X_{i}$ are vectors of length $n$ and $Y$ is the Cartesian product of these vectors under the operator $+$. This problem, which is important to other problems such as max-convolution[2] and max-product Bayesian inference[8, 6], is the generalization of the pairwise problem in which we compute the top $k$ values in $C=A+B$.

Finding the top $k$ values in $C=A+B$ is nontrivial. In fact, there is no known approach that sorts all $n^{2}$ values of $C$ faster than naively computing and sorting them in $O(n^{2}\log(n^{2}))=O(n^{2}\log(n))$[1]; however, Frederickson & Johnson demonstrated that the top $n$ values in $C$ can be computed $\in O(n\log(n))$[4]. Frederickson & Johnson generalized the problem to obtaining the $k^{\text{th}}$ top value in a matrix which is sorted by both columns and rows[3]. A sorted matrix may be built with $X_{1}+X_{2}$, but the method presented assumes the matrix is already of this form and does not take into account the work to produce the matrix from the vectors.

This can be observed by sorting both $A$ and $B$ in descending order, and then sparsely building a matrix of $A_{i}+B_{j}$. If $A^{\prime}$ and $B^{\prime}$ represent sorted vectors such that $A_{1}^{\prime}\geq A_{2}^{\prime}\geq\cdots$ and $B_{1}^{\prime}\geq B_{2}^{\prime}\geq\cdots$, then the maximal value of $C$ is $A_{1}^{\prime}+B_{1}^{\prime}$. The second largest value in $C$ is $\max(A_{1}^{\prime}+B_{2}^{\prime},A_{2}^{\prime}+B_{1}^{\prime})$.

W.l.o.g., we know that we will never insert $A_{i+1}^{\prime}+B_{j}^{\prime}$ into the sorted result of top values in $C$ without first inserting the larger (or equal) value $A_{i}^{\prime}+B_{j}^{\prime}$. Thus, each time a value from the matrix $A_{i}^{\prime}+B_{j}^{\prime}$ is found to be the next largest in $C$, the subsequent next largest value in $C$ may be in the unvisited existing neighbors (row or column, but not diagonal) of previously unvisited values. These considered values form a “fringe” around the indices whose values are already inserted into the result. We find the next value to insert into the result by finding the minimum value currently in this fringe. Because $\leq 2$ values will be inserted into the fringe in each iteration (i.e., if $A_{i}^{\prime}+B_{j}$ was just inserted into the result, then $A_{i+1}^{\prime}+B_{j}^{\prime}$ and $A_{i}^{\prime}+B_{j+1}^{\prime}$ will be added into the fringe if they are in bounds). Thus, the minimum value in the fringe can be updated sparsely by using a binary heap. Note that only the indices and values comprising the fringe are the only values stored, and that the full matrix, which would have space $n^{2}$ and thus runtime $\in\Omega(n^{2})$, is never realized. The fringe can never have size $\in\omega(n)$ (because in the worst-case scenario, it moves $n$ steps up and $n$ steps to the right), and so each next largest value in $C$ will have cost $\in O(\log(n))$. Computing the top $k$ values in $C$ is thus $\in O(k\log(n))$.

In this paper, we first construct a direct generalization of Frederickson & Johnson’s method to the problem of finding the top $k$ values in $Y=X_{1}+X_{2}+\cdots+X_{m}$. We then create a more efficient method by generalizing Frederickson & Johnson’s $C=A+B$ method to the case where $A$ and $B$ are arbitrary, heap data structures, and compute the largest $k$ values in $Y$ by constructing a balanced binary tree whose nodes each are one of these data structures. This method is then applied to finding the most abundant isotope peaks for a given molecular formula.

2 Methods

2.1 A direct $m$-dimensional generalization of Frederickson & Johnson

The direct generalization of Frederickson & Johnson’s method, which closely resembles Łącki et al.’s approach to generating the top isotope peaks[5], is straightforward: Instead of a matrix, we have an $\mathbb{R}^{m}$ tensor. As before, we store only the current fringe, which is stored in a heap. In each iteration, we remove the minimal value in the fringe and append that to the result vector. Let this minimal value come from index $(i_{1},i_{2},\ldots i_{m})$. Now we insert the $m$ values from the neighbors of index $(i_{1},i_{2},\ldots i_{m})$ into the heap holding the fringe: $(i_{1}+1,i_{2},\ldots,i_{m}),(i_{1},i_{2}+1,\ldots,i_{m}),\ldots(i_{1},i_{2},\ldots,i_{m}+1)$. As with the two-dimensional method, it is possible to store not only the $X_{1,i_{1}}+X_{2,i_{2}}+X_{3,i_{3}}+\cdots X_{m,i_{m}}$ in the heap, but also store the index tuple from which it came. This is shown in Listing 2.

This $m$-dimensional method will be substantially slower than the two-dimensional method: the fringe in this version of the $m$-dimensional method is a plane of width $\leq n$ and dimension $m-1$, and thus can have size up to $O(n^{m-1})$.

2.2 A hierarchical $m$-dimensional method for finding the top $k$ values of $Y$

First, observe that in the two-dimensional case $C=A+B$, it is not necessary for $A$ and $B$ to be sorted vectors; instead, it is sufficient that $A$ and $B$ simply be max-heap data structures, from which we can repeatedly request the next largest value. The method thus runs similarly to the two-dimensional case: the top-left corner of $C$ is computed via the maximal value of $A$ and the maximal value of $B$. This inserts two values into the fringe: either the sum of the largest value in $A$ and the second-largest value in $B$ or the sum of the second-largest value in $A$ and the largest value in $B$. Neither the full, sorted contents of $A$ nor the full, sorted contents of $B$ are needed.

We thus construct a balanced binary tree of these heap-like structures:

[TABLE]

Each heap-like structure is of the form $C=A+B$, where $A$ and $B$ are heap-like structures (Figure 1). The base case (at the leaves) is achieved by simply using a binary heap of an input vector (Listing 4).

2.3 Computing the most abundant isotope peaks

Let carbon be represented as the vector $C=(\log(0.9892),\log(0.0108))=(-0.0108,-4.5282)$ and hydrogen as the vector $H=(\log(0.9999),\log(0.0001))=(-0.0001,-9.2103)$.[5] Propane, \ceC3H8, has abundances composed of $Y=C+C+C+H+H+H+H+H+H+H+H$. To reduce the number of vectors from 11 to two, a multinomial for each element is computed to find the contributions from all isotopes. This multinomial is represented as a vector which takes the place of all vectors for that element. Thus, propane could be estimated as the sum of two multinomials, which are encoded as vectors: one for $C$ and one for $H$. Molecules that consist of any amounts of hydrogen, carbon, nitrogen, oxygen, sulfur, etc. can solved using vectors: $H,C,N,O,S\ldots$. The most abundant isotope peaks are found via the top values in $Y$.

The specific isotopes to which these abundances correspond (and from which we can compute the masses that correspond to each abundance) can be computed easily from the tuple indices. Python code implementing the hierarchical method to calculate the most abundant isotope peaks is shown in Listing 5.

3 Results

3.1 Efficient isotope peak computation

On a fake compound composed of

\ceCl800V800He800C800H800N800O100S6Cu800Ga800Ag800Tl800Ne800 with the cost of computing the multinomial isotopologues not included (this cost was identical for both methods and required $\approx 1$ second), computing the top 512 isotope beaks via a TensorCartesianSumHeap.py-based implementation took 0.002984046 seconds, while the TreeCartesianSumHeap.py-based implementation took 0.00146198272 seconds.

3.2 Time and memory use on arbitrary $Y=X_{1}+X_{2}+\cdots+X_{m}$

Problems of various sizes $m$ were run with vector length $n=m$, retrieving the top $m$ values in $Y$. The brute-force approach was not considered, for efficiency reasons. Time and memory use of TensorCartesianSumHeap.py and TreeCartesianSumHeap.py on a dual Xeon with 128GB of RAM are shown in Figure 2.

4 Discussion

As $m$ increases, the method we introduce here is far more time efficient, but more importantly, far more space efficient than direct $m$-dimensional generalization of Frederickson & Johnson.

Although the approach we propose here does have benefit to computation of intense isotope peaks, the limited number of elements (currently at $m=118$) benefits only slightly from our approach. Furthermore, it is rare for many elements to be combined in a single compound. Our demonstration implementation generates multinomials naively, unlike Łącki et al.; the longer runtime from this unnecessary naïvete in preprocessing mutes the speedup of the hiercharical method for elements with several isotopes that are probable; however, it is possible to use this hierarchical approach but with multinomials generated non-naively as Łącki et al. do, which would likely achieve a modest speedup over their current approach.

However, there are cases outside of computation of intense isotope peaks, in which the hierarchical method we propose would yield large practical benefit. These include maximum a posteriori Bayesian inference on dependencies of the form $Y=X_{1}+X_{2}+\cdots+X_{m}$[7, 6]. Operations research applications include financial markets, e.g., retrieving the $k$ lowest overall bids in each sector of supply lines for a product.

Using a language like javascript, this method can be parallelized easily by parallelizing heap pop operations in nodes at the same layer of the tree. If enough hardware parallelism is available, the runtime of a full propagation through the tree would thus be the height of the tree, which is $\in\Theta(\log(m))$ pop operations from fringe heaps. Each of these fringe heap pop operations is $\in O(\log(n))$, and thus the runtime would be $\in O(\log(m)\log(n))$.

5 Availability

Python source code from this method is available at https://bitbucket.org/seranglab/theodolite/ (MIT license, free for both academic and commercial use).