Multi-scale Sinusoidal Embeddings Enable Learning on High Resolution Mass Spectrometry Data

Small molecules in biological samples are studied to provide information about disease states, environmental toxins, natural product drug discovery, and many other applications. The primary window into the composition of small molecule mixtures is tandem mass spectrometry (MS2), which produces data that are of high sensitivity and part per million resolution. We adopt multi-scale sinusoidal embeddings of the mass data in MS2 designed to meet the challenge of learning from the full resolution of MS2 data. Using these embeddings, we provide a new state of the art model for spectral library search, the standard task for initial evaluation of MS2 data. We also introduce a new task, chemical property prediction from MS2 data, that has natural applications in high-throughput MS2 experiments and show that an average $R^2$ of 80\% for novel compounds can be achieved across 10 chemical properties prioritized by medicinal chemists. We use dimensionality reduction techniques and experiments with different floating point resolutions to show the essential role multi-scale sinusoidal embeddings play in learning from MS2 data.


Introduction
Tandem Mass Spectrometry (MS/MS) is one of the primary tools used in metabolomics, the characterization of the small molecule ( 1,000 Daltons) contents of complex biological samples. Tandem Mass Spectrometry works by measuring with very high precision the masses of molecules and their constituent fragments. While MS techniques are highly sensitive and extremely precise, inferring the identity of the molecules and their properties from the resulting mass spectra is one of metabolomics' primary bottlenecks [7]. Unlocking the information contained in mass spectra will lead to an improvement in a number of applications including characterization of disease pathways, development of new agrochemicals, improved forensics analysis, and the discovery of new drugs [48]. same manner. The resulting spectrum contains the m/z's and intensities (together, the "peaks") of all resulting fragments, and is called the MS2 spectrum. See [9].
In recent years, several machine learning methods have been developed to identify the structures and properties of small molecules from their mass spectra. These approaches [13,14,19,26,36,41] historically discretize m/z (via tokenization or binning). However, the m/z values obtained in modern mass spectrometry experiments are collected with parts per million levels of precision. There are three critical reasons to expect that modeling m/z values as numeric quantities is the appropriate technique. First, we know that relevant chemical information is present at the millidalton level [15,32] and discretization schemes typically strip this information away. Second, mass differences between peaks represent fragmentation patterns, and are therefore relevant to understanding the molecule. Learning to utilize this information from tokens is exceedingly difficult. Finally, any laboratory measurement always comes with some degree of experimental error, and discretization techniques are particularly susceptible to measurement errors. In light of this, we proceed by modeling m/z as a numerical value using sinusoidal embeddings, which we hypothesize will enable us to capture information across many orders of magnitude.
In this work, we apply a numerical representation of m/z that uses sinusoidal embeddings across multiple scales to retain the information content of MS2 data across its entire mass resolution. We use these embeddings to achieve state of the art accuracy on spectral library search [39], a task which identifies compounds by comparing mass spectra from experiments to spectra from previously characterized compounds. This task is the primary method used in standard metabolomics analyses, but is challenging because spectra from a given compound vary widely according to experimental conditions. We also introduce a novel task of predicting chemical properties relevant to drug discovery from MS2 data and apply the same modeling approach to this new task. We achieve 80% average R 2 for out of sample molecules across 10 properties on this task, which is high enough to enable first-pass filtering and selection of candidate drug molecules in high-throughput experiments based solely on spectral data. Finally, we qualitatively demonstrate that the multi-scale sinusoidal m/z embeddings capture high resolution structure by studying UMAP [28] projections of the embeddings.

Related work
Modeling numerical data in terms of sinusoidal functions has a long history in many scientific fields. In machine learning, sinusoidal embeddings are most commonly used to encode the discrete positions of natural language token inputs in transformer models, which are otherwise position-agnostic [42]. Other work has used sinusoidal embeddings for multi-dimensional positional encoding in image recognition [23]. In mass spectrometry, sinusoidal embeddings have also been described in the context of modeling protein sequences from mass values, either with [33] or without [46] an initial mass binning step.
For modeling m/z values in the field of metabolomics, previous machine learning models have relied primarily on discretization of the continuous mass inputs. This is usually accomplished by binning m/z values into fixed length vectors with peak intensity as the value for each element. Various authors have used binned representations of spectra for spectral library search [14], unsupervised topic modeling [41], and molecule identification [19,26]. Alternatively, masses have been tokenized via rounding for tasks such as unsupervised spectral similarity [13] and molecule prediction from synthetic data [36].
Other discretization approaches rely on tokenization of m/z values by direct molecular formula assignment. Even when this method is able to uniquely identify molecular formulas for every fragment ion, it still faces a problem endogenous to tokenization: it's very difficult to reason about the relationships between discrete tokens except by pattern recognition over a very large number of training examples. Böcker et al. [4] addressed this difficulty by modeling the fragmentation process with a tree structure, and Dührkop et al. [6] used these fragmentation trees as inputs to a molecular fingerprint prediction model that is used to assign molecules to spectra from molecular structure databases.
is composed of a precursor m/z and a set of N fragment peaks at various m/z's and intensities (I).
Various other data are typically collected, including precursor abundance, charge, collision energy, etc, but these are not used in this work. Transformer encoders [42] without the positional encoding are explicitly fully-symmetric functions and hence are ideally suited to model a set of fragmentation peaks. We therefore take our base model, SpectrumEncoder (S), to be a transformer encoder, whose inputs are a set of (m/z, I) pairs that includes the precursor along with all of the fragment peaks. We normalize the intensities to a maximum of 1 for fragments and assign an intensity of 2 to the precursor. Finally, we take as output the embedding vector from the final transformer layer corresponding to the precursor input. See Figure 2. Note, the use of transformers to study synthetic MS2 spectra has been explored in [36].
This flexible approach enables us to experiment with and compare various representations of MS2 data. We pursue two approaches: tokenization of m/z and modeling m/z as numerical values via sinusoidal embeddings. Before describing these, we need the following definition for a simple two layer feed forward MLP, Note that each occurrence of FF we employ below is a separate instance with distinct weights.

Tokenized m/z peak embedding
Our first approach is to discretize m/z by rounding to 0.1 precision and treating these objects as tokens (as in an NLP context). These tokens are embedded in a dense vector space as in [13,29], via an embedding function TE. We then construct the token peak embedding, where denotes concatenation. A diagram of PE token is shown in Figure 1a.

Sinusoidal m/z peak embedding
We now consider a numerical peak embedding that utilizes the full m/z precision. We employ a sinusoidal embedding as follows: The frequencies are chosen so that the wavelengths are log-spaced from λ min = 10 −2.5 Daltons to λ max = 10 3. 3 Daltons, corresponding to the mass scales we wish to resolve. This embedding is inspired by [23,42]. Similarly to Equation 3, we construct the sinusoidal peak embedding, PE sin (m/z, I) = FF (FF (SE (m/z)) I) .
We also experiment with casting the m/z inputs as half, single, and double precision floating point values. Equations 4-5 are computed at the precision specified by m/z and the results are then cast to the precision of the model weights used for training. A diagram of PE sin is shown in Figure 1b.   Figure 1b (right) is produced by a sinusoidal embedding followed by a two layer MLP. In both cases, an intensity value I is appended to make a vector of dimension d + 1. This is then fed into a simple feed-forward network to produce a peak embedding of dimension d. We use blue to highlight the network components involved in producing an m/z embedding and purple to highlight the full peak embedding after intensity is incorporated.

Transformer
The set of fragment peaks and the precursor are passed through either of the two peak embeddings specified above. The result is a sequence of embedding vectors which is then passed through a transformer encoder: Here PE stands in for either PE token or PE sin . TransformerEncoder has embedding dimension d = 512 and six layers, each with 32 attention heads and an inner hidden dimension of d. Note that all hidden and embedding layers in our peak embeddings above have dimension d as well. As mass spectra have no intrinsic ordering, we opt to not include a positional encoding.
In order to get a single embedding vector as an output, the final transformer layer query only attends to the first embedding, corresponding to the position of the precursor m/z. A diagram of our full model architecture is shown in Figure 2.

Model training
We train the base models discussed here via two independent tasks aimed at applications within metabolomics and medicinal chemistry. These tasks are described in detail below. We build these models using PyTorch [8,30] and train them with the Adam [16] optimization algorithm with parameters β 1 = 0.9, β 2 = 0.999, and learning rate α = 5.0 × 10 −5 . We also use a weight decay [18] parameter of 0.1, a gradient clipping a parameter of 0.5, and a dropout [38] parameter of 0.1. All models were trained using between 25 and 50 epochs. Finally, all model weights were trained using half precision floating point values, regardless of the precision of the m/z inputs.

Spectral similarity with Siamese networks
A common workflow for the analysis of mass spectrometry data is to compute some notion of spectral similarity for pairs of spectra that is intended to correlate with the molecular similarity of the underlying compounds [47]. This can then be used for tasks such as spectral library search [39] and molecular networking [34,45]. Historically, most spectral similarity metrics used in the field have been heuristic-based [24,47]. More recently, spectral similarity methods that make use of techniques from deep learning have been developed [13,14]. The latter of these, MS2Deepscore, is a Siamese network that trains a simple MLP (by binning m/z) to predict a measure of molecular similarity [3] from pairs of mass spectra. The loss function is the M SE between the ground truth molecular similarity and the cosine similarity evaluated on the dense embedding vectors generated by the base model.
Molecular similarity scores range from 0 to 1. The scores of randomly sampled pairs from our labeled dataset are not uniformly distributed and strong similarity matches are rare. Huber et al., in [14], describe a procedure to sample pairs of spectra such that the molecular similarities will be uniformly distributed. We use this procedure to sample pairs for training and evaluation tasks.
Following Huber et al. [14] we also train a Siamese network, but we do so using the base model transformer architectures outlined in Section 2.1 instead of with the binned mass vectors used by Huber. We report performance results for models trained using PE token and PE sin below in Section 3.2.

Property prediction
We also train the base model architecture to predict 10 chemical properties relevant to medicinal chemistry directly from an MS2 spectrum. We do so by training the following network, where the final layer has dimension equal to the number of properties we wish to predict. We also scale all of our training labels to zero mean and unit variance. For inference, we run predicted properties through an inverse-scaler. As in Section 2.2.1, we train multiple versions of this property prediction model. Due to the lack of existing methods for comparison, we also train a version of Equation 8 on binned spectrum representations for a simple baseline. We replace SpectrumEncoder with the MS2Deepscore feed forward architecture found in [14] for this baseline.
The properties we predict are standard indicators of druglikeness and bioavailability in the field of medicinal chemistry. These include the properties that make up the Quantitative Estimate of Druglikeness (QED) along with several others [5,25,43]. For the complete list of properties, see Table 2. All properties are computed deterministically from chemical structure using RDKit [20], so there is no additional dependency on experimental data.

Data
We construct a set of labeled MS2 spectra by combining a number of publicly available datasets [1, 2, 10, 11,12,17,21,22,27,31,35,37,40,44] with our internal proprietary dataset. We allow spectra collected in both positive and negative ion mode, and across a variety of instrument types, collision energies, and other important instrument parameters. We constrain that spectra in our dataset have at least 5 peaks and at least 3 decimal places of m/z resolution. We strip all stereochemistry from our molecular structure labels, which is a common step taken in MS2 modeling and allows for better molecule-disjoint splitting. Chemical properties and fingerprints are then computed from the cleaned molecules using RDKit [20]. Our resulting dataset has 1,251,830 spectra corresponding to 45,351 distinct molecules.
In metabolomics experiments, biological samples contain both compounds that have previously been profiled in MS2 libraries and compounds that have not been profiled. Therefore, we separately characterize performance on "known" molecules that have been previously profiled with MS2, and "novel" molecules for which no labeled MS2 spectra are available. We partition our set of labeled spectra into a training, development and test set of size 1,214,812, 18,750, and 18,268 respectively. All spectra associated with 1,002 randomly selected molecules are partitioned together into each of the development and test sets. In this manner, a subset of both the development and test sets is molecule-disjoint from the training set. Henceforth, we refer to any spectrum in the development or test set as corresponding to a known molecule if that molecule is present in the training set, and a novel molecule if it is not. Our resulting development and test sets each have 998 known spectra and 17,752 and 17,270 novel spectra, respectively. While a subset of molecules in our development and test sets overlap with our training set, the spectra in all three sets are fully disjoint. We use the development set for hyperparameter optimization and only access the test set to compute reported metrics.

Spectral similarity
A key application of any spectral similarity method is spectral library search [39] (or dereplication), whereby an unknown MS2 spectrum is identified by searching through a database of curated spectra obtained from known molecules. To benchmark our various spectral similarity models, we query spectra from our test set against a lookup reference (which is the same as our training set). Every query spectrum is matched, by maximum similarity, to a reference spectrum. We then assess accuracy on a set of spectra queries by determining the fraction of matches that exceed a particular molecular similarity threshold. We use two such thresholds, 0.60 for approximate match and 0.95 for nearly exact match. Finally, we typically have a varying number of recorded spectra per distinct molecule so all reported accuracies are macro averaged over query molecules.
Because novel molecules do not occur in the reference set, we estimate the maximum theoretical performance on the prediction task for novel molecules by evaluating the best molecular similarity matches between the novel test set and the reference set. A spectral similarity model that perfectly approximates molecular similarity would score 0.31 and 0.89 on nearly exact and approximate matching, respectively. We also evaluate the performance of a random baseline by uniformly sampling pairs of spectra and computing their corresponding molecular similarity. Only 0.713% and 0.03% of random spectra pairs from our data have a corresponding molecular similarity above the threshold for an approximate and nearly exact match, respectively.
We train a number of spectral similarity models, as described in Section 2.2.1, and will refer to the sinusoidal Siamese transformer and tokenized Siamese transformer models as the sinusoidal and tokenized models respectively. For baselines, we train both a Spec2Vec [13] and an MS2Deepscore [14] model on our data. In Table 1 we report spectral library search accuracies (based on molecular similarity thresholds) for both known and novel compounds obtained from a number of spectral similarity models. In the same table we also report the M SE loss between predicted and actual molecular similarity on spectrum pairs drawn from train, known, and novel spectra sets. No M SE can be reported for Spec2Vec as this method does not predict a molecular similarity.
We find that MS2Deepscore outperforms Spec2Vec and is competitive with our tokenized model on novel spectral library search, but underperforms on known spectral library search. MS2Deepscore is substantially less accurate at predicting the specific pair molecular similarity on known and novel spectra when compared to all transformer approaches. The double precision sinusoidal model produced the best performance on all evaluation tasks, beating all baselines. This model produces state of the art known spectral library search accuracies of 0.943 and 0.97 for nearly exact and approximate matching respectively. Numerical representation of MS2 m/z via sinusoidal embeddings consistently outperforms all the discretization based approaches we benchmarked.   Table 1, drops to be on par with the tokenized model (worse on known and better on novel). We do not report results with single precision m/z since these were indistinguishable from double precision m/z. On all evaluation metrics, using m/z inputs of precision greater than 16 bits dramatically improves performance of the sinusoidal model. From this, we draw the implication that sinusoidal embedding of m/z is able to effectively learn from with the wide range of scales present in high precision mass spectrometry measurements and this partially accounts for the superior performance of this method.

Qualitative embedding analysis
To further characterize the respective properties of token and multi-scale sinusoidal embeddings, we inspect UMAP [28] projections of our m/z embeddings in Figure 3. For this analysis, we use siamese transformer models described in Section 2.2.1. For our sinusoidal models, we embed 50,000 m/z values between 0 and 1,000 Daltons using FF (SE (m/z)). We do not use the full peak embedding function as we are not interested in the embedding of intensity information. Because our tokenization procedure involves rounding m/z values to 1 decimal place, for the tokenized model we only embed 10,000 m/z values between 0 and 1,000 Daltons using the embedding function TE (m/z).
While m/z measurements are continuous, the space they represent is inherently discrete. Two molecular fragments differing in atomic composition (let alone structure) could nevertheless have a very similar m/z's. However, the large separation between molecular mass scales and the mass scale at which nuclear forces manifest themselves means that the relative atomic composition may be deduced from the fractional m/z [15,32], defined by, Therefore, we expect a quality embedding of m/z to embed fragments similar in both m/z and {m/z} close to one another, while paying greater attention to the latter.
As is seen in Figure 3, embeddings from tokenized and low precision sinusoidal models are able to capture general trends in m/z. However, they fail to preserve distances in m/z and show little to no structure in {m/z}. In contrast, the high precision sinusoidal embeddings allow our models to represent important information in {m/z}, while preserving distance in m/z.

Property prediction
Many applications in metabolomics can be unblocked with knowledge of just a limited set of chemical properties, without need for identification of molecular structure. Towards this end, we introduce a novel task of chemical property prediction from MS2 spectra. We evaluate sinusoidal vs. tokenized m/z embedding performance on a list of properties that have compelling applications in drug discovery and are easily computable from our molecule labels. For this experiment, we train a sinusoidal model and several baseline models as described in Section 2.2.2. In Table 2, property prediction models are evaluated on known and novel structures using R 2 . All transformer models outperform our feed forward baseline, illustrating the difficulty of the task and the significance of our results. Half precision sinusoidal models produced an R 2 of 0.746 (0.936) on known (novel) molecules averaged across all 10 properties. These results are competitive with tokenization but substantially underperform when compared to higher precision sinusoidal embeddings. When m/z values are specified at double precision, sinusoidal embeddings show improvements over tokenization on each property, resulting in an 11% improvement over tokenization on R 2 for novel molecules averaged across all properties.
This is further evidence that the multi-scale information captured by sinusoidal embeddings allows models to better understand the nuances of MS2 and chemical space. In addition, the substantial performance of our sinusoidal property prediction model is strong enough to enable prioritization over a large number of unknown metabolites. This has major implications in the drug discovery context, as inspecting the properties of unknown molecules in complex mixtures was previously not possible without isolation and individual experimentation.

Conclusions
Through a number of experiments, we demonstrate that MS2 modeling methods that use a numerical representation of m/z via multi-scale sinusoidal embeddings are superior to methods which make use of discretized m/z. In order to obtain insight into the performance attributed to sinusoidal embeddings, we present evidence that indicates that sinusoidal embeddings enable learning from the high-precision data that is produced from MS2 experiments.
The modeling work represented here produces state of the art results on spectral dereplication tasks. We also produce, as far as we know, the first results on ML chemical property prediction from MS2 spectra. The property prediction results of are sufficient quality to functionally inform and advance drug-discovery efforts. We expect that sinusoidal embeddings of m/z will be a critical component of further research that incorporates MS2 spectra modeling, and we look forward to exploring their use in other contexts.