ORIGINAL RESEARCH article

Front. Anal. Sci., 08 January 2024

Sec. Chemometrics

Volume 3 - 2023 | https://doi.org/10.3389/frans.2023.1301602

A Bayesian approach for constituent estimation in nucleic acid mixture models

  • 1. Department of Reasoning for Intelligence, The Institute of Scientific and Industrial Research, Osaka University, Ōsaka, Japan

  • 2. AI Solution Unit, Technology Research Laboratory, Shimadzu Corporation, Kyoto, Japan

  • 3. Shimadzu Analytical Innovation Research Laboratories, Osaka University, Ōsaka, Japan

  • 4. Life Science Business Department, Analytical and Measuring Instruments Division, Shimadzu Corporation, Kyoto, Japan

Abstract

Mass spectrometry (MS) is a powerful analytical method used for various purposes such as drug development, quality assurance, food inspection, and monitoring of pollutants in the environment. In recent years, with the active development of antibodies and nucleic acid-based drugs, impurities with various modifications are produced. These can lead to a decrease in drug stability, pharmacokinetics, and efficacy, making it crucial to differentiate these impurities. Previously, attempts have been made to estimate the monoisotopic mass and ion amounts in the spectrum generated by electrospray ionization (ESI). However, conventional methods could not explicitly estimate the number of constituents, and discrete state evaluations, such as the probability that the number of constituents is k or k+1, were not possible. We propose a method where, for each possible number of constituents in the sample, mass spectrometry is modeled using parameters like monoisotopic mass and ion counts. Using Simulated Annealing, NUTS, and stochastic variational inference, we determine the parameters for each constituent number model and the maximum posterior probability. Finally, by comparing the maximum posterior probabilities between models, we select the optimal number of constituents and estimate the monoisotopic mass and ion counts under that scenario.

1 Introduction

Mass Spectrometry (MS) is a powerful analytical technique used for various purposes such as drug development and quality assurance, food inspection, and monitoring of pollutants in the environment. In recent years, with the active development of antibodies and nucleic acid drugs, impurities with different modifications are produced. These can cause a decrease in the stability of the drug, its pharmacokinetics, and its efficacy (; ; ; ). Therefore, it is crucial in drug development and quality assurance to distinguish these multiple impurities in pharmaceuticals and take measures against them. Moreover, if we know the monoisotopic mass of the constituents, it can provide valuable information for considering the cause of impurity generation, and if the ion amounts of the constituents are known, it helps estimate the effect of the impurity.

However, in current mass spectrometry, it is generally difficult to directly distinguish impurities in targets of middle molecules or higher that have slight modifications, and to estimate the correct number of constituents and their monoisotopic mass and ion amounts. This challenge arises because separating such impurities using conventional chromatography methods is problematic. Furthermore, the MS spectra become complex due to the isotopes in the target constituent. Especially with the most used ionization method, Electrospray Ionization (ESI), more complex spectra are produced because it creates multivalent ions, leading to many interpretations and degrees of freedom, making analysis even more challenging.

While increasing hardware resolution can differentiate slight differences between isotopes and modifications, methods like FT-ICR (Fourier Transform Ion Cyclotron Resonance), which have high resolution, require massive equipment and substantial costs, making it inconvenient to handle. Thus, it is preferable to analyze with devices that can be managed in general labs, such as Triple-Quadrupole-MS and Quadrupole -Time-of-Flight-MS(Q-TOF-MS).

Therefore, there is active research in approaching signal analysis through software. Various attempts have been made to estimate mass from mass spectrometer data. Simple methods to derive m/z lists from spectra include wavelet transformation (). Recently, peak detection algorithms have been developed that combine continuous wavelet transformation (CWT) and image processing (). This method applies image segmentation to CWT coefficients, creating masks, and combining them with CWT maximum values to improve peak detection ROC (receiver operating characteristic) curve. Such techniques are useful for resolving close mass constituents in spectra measured with low molecules, which have a relatively narrow isotope distribution, or when ionized with methods producing simple charge distributions like EI (electron ionization) or MALDI (matrix-assisted laser desorption ionization). However, for spectra of middle to high molecules with a wide isotope distribution, especially those generated by ESI, which produce multivalent ions with a charge distribution, distinguishing the monoisotopic mass of interest becomes challenging.

For charge deconvolution and deisotoping from multivalent ion spectra, a lot of algorithms such as heuristic gaussian fitting using nonlinear least squares minimization () have been proposed. The ReSpect algorithm using the Max Entropy method () has been long used (; ; ). This algorithm integrates m/z lists based on charge distribution constraints, enabling the determination of monoisotopic mass. However, ReSpect cannot explicitly estimate the number of constituents in the spectrum, and it cannot evaluate discrete states, like the probability that there are constituents or constituents. Also, the entropy term of the objective function increases as the number of peaks in the deconvolution spectrum increases, leading to the selection of spectra with many peaks.

Recently, new methods like UniDec using Bayesian deconvolution have emerged (; ). UniDec adopts a unique algorithm similar to the Richardson-Lucy method (; ) and is faster than ReSpect. However, its iterative method to approximate the observed data by a convoluted spectrum does not resolve the issue of not being able to evaluate the probability of a certain number of constituents.

In this paper, we newly propose a method to select the optimal number of constituents by comparing the probability of each constituent count, and to estimate the monoisotopic mass and ion counts under that condition. This can suggest the presence of impurities in pharmaceuticals, assist in the search for better synthesis conditions for middle to high molecular pharmaceuticals, and be useful for quality assurance in factories. For this study, we target Time-of-Flight mass spectrometers, which are frequently used in drug development due to their high sensitivity and resolution.

2 Proposed method

2.1 Analytical method framework

First, we model the mass spectrometry system based on parameters like the mass and charge of each constituent, assuming a certain number of constituents in the sample. Here, a constituent is defined as a substance with a specific monoisotopic mass. We then perform a MAP (Maximum A Posteriori) estimation of these parameters from the observed spectrum. By comparing the maximum posterior probability in models with different numbers of constituents, we determine the model with the most appropriate number of constituents.

However, this model has a large dimensionality of the number of constituents multiplied by 6. Moreover, the posterior probability for one of the parameters, the monoisotopic mass, is flat over a large portion of the search space and has several sharp peaks locally. Hence, gradient-based methods are not suitable for this case due to anticipated gradient vanishing. Therefore, to estimate the parameters, we combine the No-U-Turn Sampler (NUTS (), a type of Markov Chain Monte Carlo (MCMC), with Simulated Annealing ().

The purpose of using Simulated Annealing is to introduce a temperature parameter. By selecting a high-temperature exploration parameter distribution, we can actively explore parameters even in areas where the posterior probability is flat or has sharp peaks. This ensures a broader search across the parameter space, reducing the chance of overlooking the global solution and getting trapped in local minima.

Furthermore, NUTS can explore parameters sparsely in areas with small gradients and can explore parameters in detail in areas with large gradients. Thus, introducing NUTS allows efficient exploration of the vast, high-dimensional parameter space.

On the other hand, while MCMC is good at searching for global solutions, it does not always reach the optimal solution within a certain number of search steps. Therefore, we use the parameters with the highest posterior probabilities obtained from NUTS and Simulated Annealing as initial values and apply stochastic variational inference. By doing this, we search for the optimal parameter where the posterior probability is maximized in the vicinity of that initial value, aiming to improve the accuracy of parameter estimation.

However, simultaneously searching for parameters for all possible numbers of constituents leads to a curse of dimensionality, where the search space explosively expands as the number of constituents increases, potentially reducing search efficiency and accuracy. To avoid this problem, we sequentially increase the number of constituents from to the maximum conceivable number , and for constituents calculate the optimal parameters and their posterior probabilities, and efficiently focus the parameter search areas for the k+1 constituents by these posterior probabilities.

To balance the complexity of the model (number of constituents) and its fit (loss against the data), in addition to the prior distribution of each parameter, we introduce a prior distribution for the number of constituents. We also incorporate a prior distribution on the differences between the monoisotopic masses of multiple constituents. For analytical purposes, we have defined a single constituent as a substance with a distinct monoisotopic mass, thereby ensuring that their masses do not mutually take the same value. When seeking to separate isomers, it is essential to integrate other techniques such as fragmentation, ion mobility spectrometry, and chromatography, in addition to the proposed method. We first construct a model with constituent, obtain the optimal parameters and the maximum posterior probability based on the above prior distributions and observed data.

Next, we construct a model with constituents. For one of the two constituents, we use a prior distribution centered on the optimal parameters already estimated for , narrowing its range. This suppresses the significant increase in the parameter search space. Based on this new prior distribution, we estimate the optimal parameters and obtain the maximum posterior probability.

Subsequently, we seek the maximum posterior probability for each model with constituent numbers up to the upper limit by efficiently exploring the optimal parameters in the same manner.

Finally, we compare the maximum posterior probabilities corresponding to each model with different numbers of constituents. We select the model with the highest probability and obtain the estimates for the monoisotopic masses and ion counts.

2.2 Physical model of mass spectrometers

In a time-of-flight mass spectrometer, the relationship between time of filght , mass , and charge can ideally be represented as in Eq. 1 (). Here, represents the acceleration voltage to eject ions, and is the flight distance of ions.

The actual signal obtained is a convolution of the delta function with the detector’s response waveform . Preliminary experiments have shown that has a stochastic ToF jitter following the logarithmic normal distribution (). Additionally, the height (intensity) of the response waveform also has a stochastic response factor following the logarithmic normal distribution .

Reflecting these variances in time of flight and detector response intensity, the response when a single ion with mass m and charge z enters the detector according to Eq. 1 can be written as:

The mass and charge of constituent follow a distribution as described below. Since the natural isotopic abundances differ by element, and each functional group has different charge rates, the probability mass functions of mass and charge are represented by a multinomial distribution as:

However, this model becomes too complex for Bayesian inference due to its large number of parameters . Therefore, we introduce an approximation that assumes equal isotopic abundances for all elements, an increase of one neutron in isotopes, and equal charge probabilities for all functional groups.

Consequently, the mass distribution of constituent j can be represented by a binomial distribution based on the monoisotopic mass as a reference, where atoms contained in a molecule are replaced by isotopes with an increase of +1Da at a probability of . Similarly, for the charge distribution, functional groups can be represented by a binomial distribution where they acquire a charge of +1 (in the case the mass spectrometry system is in positive mode) at a rate of . This reduces the number of parameters down to six dimensions.

Considering that constituent j contains up to ions, its spectrum is represented as:

If the ion counts are sufficiently large, it can be approximated as:

Noise in the time-of-flight mass spectrometer is known to be stationary and follows a normal distribution based on preliminary experiments. It is known that the thermal noise of detection circuits in such as mass spectrometers follows a normal distribution (), suggesting that in this case, thermal noise is the dominant factor in overall noise. Therefore, the noise to be added to the entire spectrum is represented as .

Considering these, the conclusive spectrum combined with the multiple spectra of single constituent is represented as:

2.3 Bayesian estimation of number of constituents and parameters

When the observation data from the mass spectrometer is obtained, assuming the number of constituents as , the posterior probability distribution for parameters is defined as per Bayes’ theorem. Note that represents the likelihood of parameters θ when is provided, and denotes the prior distribution.

Here, in addition to the prior distribution of each parameter (uniform distribution), we incorporate a regularization term, to achieve a suitable balance between model complexity (number of constituents) and model fit (loss with respect to data). We also introduce a regularization term, , to prevent multiple constituents within the same model from assuming the same monoisotopic mass. Hence, we introduce the following logarithmic prior distribution:

To determine the appropriate number of constituents , we define the regularization term representing the complexity of the model with based on the Bayesian Information Criterion (BIC). The BIC is a statistical measure that balances the fit to the data and model complexity (; ). Here, represents the dimension of the observation data , which in this study is the number of data points in the time direction.

Furthermore, we define a constituent by its unique monoisotopic mass. Therefore, if the estimated values of the monoisotopic mass parameters of multiple constituents are the same in the algorithm, the count of constituents will not be accurate. Here, we define the logarithmic prior distribution (regularization term) as follows, using a penalty that increases exponentially according to the difference in estimated monoisotopic mass values, as shown in Figure 1. The integral of the spectrum is also multiplied as a coefficient to ensure that the impact of the penalty does not change depending on the scale of the observed data.

FIGURE 1

Here, by substituting the parameter generated from MCMC into model (9), we obtain the spectrum as . We assume a normal distribution for the noise. The standard deviation of the noise denoted as is set to 2,000. Consequently, the logarithm of the posterior probability distribution is as follows.

Here, represents the inverse temperature. As described in 2.1, we use Simulated Annealing to ensure active parameter exploration in flat areas or sharp peaks of posterior probability. This is achieved by multiplying the inverse temperature to the posterior probability. Initially starting from a low inverse temperature value (i.e., high temperature) and gradually increasing to a higher value (i.e., low temperature). At low inverse temperatures (high temperatures), the system explores a wide parameter space. Conversely, at high inverse temperatures (low temperatures), the system converges to the optimal solution. This time, we set the temperature change in three stages: .

To obtain the maximum posterior probability and parameters that maximize the posterior probability (formula (13)), we conduct sampling from this posterior probability distribution using MCMC.

2.4 Parameter exploration and optimization

From the posterior probability distribution , we sample the parameter to select the one that maximizes the posterior probability. We use the No-U-Turn Sampler (NUTS) for sampling, a recent and popular variant of the Markov Chain Monte Carlo (MCMC) method. NUTS is a type of MCMC, especially a derivative of the Hamiltonian Monte Carlo method (HMC) (; ).

After executing MCMC, the parameters of the maximum posterior probability obtained are inherited as initial values, and optimization of the parameters is performed using Stochastic Variational Inference [SVI (; ; )]. For more details, please refer to the Supporting Material.

2.5 Workflow for inferring constituents in a sample

The overall picture of the workflow to determine the optimal parameters and posterior probability for each assumed number of constituents from the observational data of the mass spectrometer is as shown in Figure 2.

FIGURE 2

First, as described in 2.3.1, (i) input the observational data of the mass spectrometer with dimensions of flight time and ion counts. Then (ii) assume that the number of constituents, , contained in the sample is 1. (iii) Set the inverse temperature to .

Next, as described in 2.3.2, (iv) sample times from the posterior probability distribution, (v) set the MAP solution obtained by MCMC as the initial value for the next MCMC. Then (vi) divide the inverse temperature by 0.2. Repeat steps (iv) to (vi) three times.

As described in 2.3.3, (vii) set the parameter of the maximum posterior probability obtained by MCMC as the initial value for SVI. Then (viii) optimize the parameters with SVI, and (ix) set the parameter of the maximum posterior probability obtained by SVI as the initial value for the next MCMC.

Following 2.3.1, increase the number of constituents, , by 1. Repeat steps (iii) to (vi). (x) Continue this until the maximum possible number of constituents, . Finally, (xi) compare the maximum posterior probabilities of models from constituents to , and (xii) select the model with the largest posterior probability. Also, obtain the optimal parameters at that time.

3 Results

3.1 Validation environment

The specifications of the PC used for verifying the proposed method, as well as the software versions, are detailed in Table 1. The proposed method handles data with 1 million dimensions along the time axis, requiring a large memory size. Additionally, to rapidly explore a wide 6-dimensional parameter space using MCMC, the high-speed probabilistic programming library, NumPyro, along with its compatible CUDA and GPU, were used.

TABLE 1

CPUIntel(R) Xeon(R) CPU E5-2698 v4 @ 2.20 GHz
GPUTesla V100-DGXS-16GB
RAM264 GB
OSUbuntu 20.04.6 LTS
SoftwarePython 3.8.10
Numpyro 0.11.0
jax 0.4.7
CUDA 11.8

Validation environment.

3.2 Creation of simulation data for validation

Based on the nucleic acid drug Fomivirsen () (ID: A), four impurity constituents with modified base sequences were added, and spectra for a total of five constituents were generated via simulation. Specific values are as per Table 2. This enables the replication of a system where the principal constituent’s isotopic distribution and the impurity spectra are mixed.

TABLE 2

IDSequenceMolecular FormulaMonoisotopic Mass [Da]Representative Functional Group Number Representative Charge Rate Ion Counts
Agcg​ttt​gct​ctt​ctt​ctt​gcgC204H263N63O134P206361.0882240.035200,000
Bgcg​ttt​gut​ctt​ctt​ctt​gcgC204H262N62O135P206362.0722240.035200,000
Cgug​ttt​gut​ctt​ctt​ctt​gcgC204H261N61O136P206363.0572240.035200,000
Dgug​ttt​gut​ctt​ctt​ctt​gugC204H260N60O137P206364.0422240.035200,000
Egug​ttt​gut​ctt​utt​ctt​gugC204H259N59O138P206365.0272240.035200,000

Settings for constituent spectrum generation.

Ion counts for each constituent were set at 20,000. To facilitate the interpretation of results and to ensure that the algorithm treats each constituent fairly, we will conduct evaluations using a 1:1 concentration ratio for each component in the proposed method. The atomic counts for H, C, N, O, and S were obtained from the molecular formula of each constituent. Natural isotopic abundance ratios followed the NIST Atomic Weights and Isotopic Compositions for All Elements. The representative functional group number and the representative charge rate were set to 224 and 0.035, respectively, to ensure that the generated spectra resembled real data.

The procedure involved sampling from the multinomial distribution represented by Eqs 3 and 4 20,000 times (total incoming ion counts) for each constituent. Subsequently, spectra were formed following the procedures in Eqs 2 and 7.

The mutation from C (Cytosine) to U (Uracil) is called deamination and is generated in the synthesis process due to solvent conditions and thermal stress (; ).

The spectra of the generated single constituents A to E were combined according to Eq. 9 in the 15 combinations listed in Table 3. This allows for a comprehensive combination of 2-3 constituents based on constituent A, as well as an evaluation of each individual constituent. We use these as test data.

TABLE 3

Mixture No.Constituents
1A,B,C
2A,B,D
3A,B,E
4A,C,D
5A,C,E
6A,D,E
7A,B
8A,C
9A,D
10A,E
11A
12B
13C
14D
15E

Combinations of constituents when generating spectra.

3.3 Evaluation of constituent count estimation accuracy

The results of estimating the number of constituents in the spectra of the test data (Mixture No.1∼15) using our proposed method are as shown in Table 4. The values within the table represent the negative logarithm of the maximum posterior probability in the model of constituent count k. Therefore, the smallest value should be selected.

TABLE 4

Mixture No.k = 1k = 2k = 3k = 4k = 5
14373749.7673756983.6913752996.6593758705.4213763612.302
24278474.5703765457.0603753649.1503756752.5293762496.339
34194715.3363771155.2813759533.9793765154.9703763711.852
44219672.2983748972.0913754868.3303761667.3203763951.346
54319572.8183773245.5363757747.3073758337.6773763773.03
63824787.1153750044.8063752258.2073757899.0923763075.106
73798372.6833746176.0703747192.8053752543.8883758004.000
83795441.1963744561.0273748947.3323756355.2903759389.714
93824787.1153750044.8063752258.2073757899.0923763075.106
103825137.5933769113.7743758565.4093769515.0123771413.809
113733728.2033738454.0983743333.7703748732.1243754347.052
123736353.9103739259.2613744921.2133750141.5793755513.083
133734850.6883738732.1373743820.6673751713.9793754223.045
143735192.1973740628.6773745980.7213751377.4483755751.922
153734867.2463738788.0103744300.3273749555.6343755907.017

Negative logarithm of the maximum posterior probability assuming each constituent count (Orange background indicates the true number of constituents, blue text indicates the minimum value across models).

By choosing the most suitable number of constituents based on this criterion, the success rate for estimating the true number of constituents was 80% (12/15). Additionally, the presence or absence of impurities (distinguishing between k = 1 and k≧2) could be determined with 100% accuracy. We believe this is sufficient as a standard for recognizing the presence and number of impurities in pharmaceuticals and taking appropriate measures.

3.4 Accuracy of parameter estimation with maximum posterior

The optimal monoisotopic masses and ion counts estimated in the model where the posterior probability is maximum for each test data are shown in Table 5. The monoisotopic mass had an average error of 1.348 Da and a maximum error of 4.931 Da. This is insufficient to determine how many mutations have occurred, making it unsuitable for examining the cause of impurity generation with a difference of 1 Da. Regarding the ion counts, there was an average error of 4% and a maximum error of 82%. For instance, the standards for total desamido impurity and total impurities in injectable glucagon are 14% or less and 31% or less, respectively (). Therefore, the accuracy of the ion count estimation in the proposed method is insufficient to estimate the impact of impurities.

TABLE 5

Mixture No.ConstituentsMass [Da] (Infer)Mass [Da] (True)Absolute Error [Da]Ion counts [ions](Infer)Ion counts [ions](True)Relative Error [%]
1A,B,C6358.0736361.088−3.015138,290200,000−31%
6361.0886362.072−0.984299,930200,00050%
6363.0476363.057−0.010172,510200,000−14%
2A,B,D6360.0886361.088−1.000207,760200,0004%
6361.0436362.072−1.029270,470200,00035%
6361.0816364.042−2.961132,170200,000−34%
3A,B,E6359.0476361.088−2.041299,970200,00050%
6360.1036362.072−1.969239,990200,00020%
6366.0086365.0270.98174,160200,000−63%
4A,C,D6360.0886361.088−1.000298,940200,00049%
6363.0436363.057−0.014299,980200,00050%
-6364.042--200,000-
5A,C,E6360.0246361.088−1.064238,440200,00019%
6361.076363.057−1.987296,510200,00048%
6361.1166365.027−3.91180,810200,000−60%
6A,D,E6360.0796361.088−1.009297,500200,00049%
6362.0276364.042−2.015299,940200,00050%
-6365.027--200,000-
7A,B6357.0886361.088−4.000191,670200,000−4%
6362.0736362.0720.001220,850200,00010%
8A,C6361.0436361.088−0.045113,870200,000−43%
6361.0806363.057−1.977283,890200,00042%
9A,D6359.0446361.088−2.044280,530200,00040%
6359.1116364.042−4.931138,700200,000−31%
10A,E6357.088--227,540--
6361.0296361.088−0.05935,400200,000−82%
6364.0106365.027−1.017157,560200,000−21%
11A6361.0886361.0880.000191,840200,000−4%
12B6361.0726362.072−1.000207,340200,0004%
13C6363.0586363.0570.001189,640200,000−5%
14D6363.0426364.042−1.000205,290200,0003%
15E6365.0276365.0270.000190,240200,000−5%

Optimal monoisotopic masses and ion counts of the model with the maximum posterior probability.

For reference, a comparison between the spectra reconstructed from the estimated parameters and the original signal is shown in Figure 3. The overall view in (a) represents the charge distribution, and the enlarged view in (b) represents the isotopic distribution. From these results, it is clear that the spectrum we generated closely matches the observed data. Despite the spectra matching, errors in parameter estimation occurred because of the high degree of freedom in isotopic parameters that trade-off with monoisotopic mass. Even if the monoisotopic mass was lower than the true value, by increasing the representative atomic number or the representative isotopic natural abundance , it is possible to make it fit the observed data to some extent.

FIGURE 3

Also, the estimated ion counts of each constituent showed errors of up to 82% from the true values. This is presumed to be due to the trade-off relationship between the ion counts of each constituent, with a decrease in the ion count of one constituent being compensated by an increase in another. This is further supported by the fact that the average error in ion counts settles at 4%.

3.5 Comparison with UniDec

Deconvolution of the test data was performed using the existing method, UniDec as well. The results of deconvolution for each observed spectrum by UniDec are shown in Table 6. According to these results, the accuracy for the correct number of constituents was 13% (2/15). This is presumed to be because the UniDec algorithm, which obtains the number of constituents after multiple iterations of deconvolution, does not necessarily guarantee the number of constituents. Please note that this use of UniDec to determine the number of constituents is not its intended application.

TABLE 6

Mixture No.ConstituentsMass [Da] (Infer)Mass [Da] (True)Absolute Error [Da]Intensity [a.u.](Infer)Intensity [a.u.](True)Relative Error [%]
1A,B,C6359.9006361.088−1.188100.000100.000100%
6360.9006362.072−1.17254.614100.00055%
-6363.057--100.000-
2A,B,D6359.9006361.088−1.188100.000100.000100%
6360.9006362.072−1.17268.122100.00068%
6361.8006364.042−2.24223.490100.00023%
3A,B,E6359.900--100.000--
6360.900--47.326--
6361.8006361.0880.71222.533100.00023%
6362.8006362.0720.72813.473100.00013%
6363.8006365.027−1.22713.496100.00013%
4A,C,D6359.900--100.000--
6360.9006361.088−0.18894.673100.00095%
6361.8006363.057−1.25764.641100.00065%
6362.8006364.042−1.24219.369100.00019%
5A,C,E6359.900--100.000--
6360.900--63.684--
6361.8006361.0880.71256.992100.00057%
6362.8006363.057−0.25733.851100.00034%
6363.8006365.027−1.22719.330100.00019%
6A,D,E6359.900--100.000--
6360.900--53.209--
6361.8006361.0880.71261.898100.00062%
6362.8006364.042−1.24270.845100.00071%
6363.8006365.027-39.057100.00039%
7A,B6359.9006361.088−1.188100.000100.000100%
6361.0006362.072−1.07211.538100.00012%
8A,C6359.900--100.000--
6361.0006361.088−0.08840.696100.00041%
6361.8006363.057−1.25710.199100.00010%
9A,D6359.900--100.000--
6361.0006361.088−0.08841.897100.00042%
6361.800--26.937--
6362.8006364.042−1.24216.351100.00016%
10A,E6359.000--19.045--
6359.900--100.000--
6360.900--27.472--
6361.8006361.0880.71219.107100.00019%
6362.800--27.075--
6363.900--33.325--
6364.8006365.027−0.22713.659100.00014%
11A6358.900--48.595--
6359.8006361.088−1.288100.000100.000100%
12B6359.900--40.161--
6360.8006362.072−1.272100.000100.000100%
13C6360.900--41.609--
6361.8006363.057−1.257100.000100.000100%
14D6361.800--52.753--
6362.8006364.042−1.242100.000100.000100%
15E6362.800--54.440--
6363.9006365.027−1.127100.000100.000100%

Deconvolution results for each observed spectrum by UniDec.

For the verification above, we used UniDec (Version 6.0.2). The particularly set parameters during this verification are shown in Table 7. The Mass Range was set to the same range as the proposed method, and Sample Mass Every (Da) was set to 0.1 to sufficiently detect impurities with a difference of 1 Da. For parameters not mentioned, default values were used.

TABLE 7

ParameterSetting value
UniDec ParametersCharge Range1 to 50
Mass Range6300–6400 Da
Sample Mass Every (Da)0.1
Additional Deconvolution ParametersIsotopesMono
Peak Selection and PlottingPeak Detection Range (Da)0.1
Peak Detection Threshold0.01

UniDec setting parameters.

*The other settings are using default values.

4 Discussion

Using NUTS, Simulated Annealing, and stochastic variational inference, we estimated parameters such as monoisotopic masses from observed data, and were able to choose the correct number of constituents with a higher probability than existing methods. This is thought to be due to the fact that we created models for each number of constituents, allowing for the comparative evaluation and selection of models for each number of constituents. This made it possible to suggest the presence of impurities in pharmaceuticals, which is useful for searching for better synthesis conditions for middle to high molecular weight pharmaceuticals, and for quality assurance in factories.

On the other hand, as shown in Table 5, the estimated monoisotopic mass had a maximum error of 4.931Da from the true value. This is thought to be due to the trade-off relationship between the monoisotopic mass and the parameters and that determine the isotopic distribution. Additionally, there was a relative error of several tens of percent from the true value in the ion counts of each estimated constituent. This is speculated to be because the ion counts of each constituent trade off with each other, with a decrease in one ion being compensated for by an increase in another ion. A potential solution to these problems is to represent monoisotopic masses and ion counts as probability distributions. By considering the uncertainty in monoisotopic masses and ion counts of constituents in the sample, improvements in estimation satisfaction can be expected.

Furthermore, it took about 50 h for deconvolution assuming 5 constituents per data. This is long compared to the few seconds to a few minutes processing time of UniDec. Also, this processing time is expected to increase almost linearly with the assumed number of constituents. Therefore, it is expected to take a long time when analyzing samples with many constituents, such as serum or environmental samples. A possible countermeasure to this problem is to divide the monoisotopic mass space into mini-batches and perform parallel calculations.

5 Conclusion

We assumed multiple numbers of constituents in the sample and created a mass spectrometry model from parameters such as monoisotopic masses and ion counts. We then sought the maximum posterior probability in the model of each number of constituents against observed data using NUTS, Simulated Annealing, and stochastic variational inference. As a result, we were able to estimate the number of constituents with high accuracy. We were also able to estimate parameters such as monoisotopic masses and ion counts at the same time.

Future challenges include reducing computation time, improving mass accuracy, and improving ion count accuracy. Incorporating chromatography or ion mobility information, addressing more stringent concentration ratios between constituents (e.g., greater than 10:1), and adapting to complex samples with more constituents will be pursued to expand the applicability of this method.

Statements

Data availability statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Funding

The author(s) declare that no financial support was received for the research, authorship, and/or publication of this article.

Acknowledgments

We extend our deepest gratitude to Yoshihiro Ueno, Yusuke Tagawa, Daisuke Okumura, and Daisuke Hiramaru for their advice and coordination on the project as a whole. We also thank Akira Noda and Yusuke Tamai for their insights on Bayesian estimation. Atsuhiko Toyama, Natsuyo Asano, Hiroaki Waki, Kiyoshi Ogawa, Hideaki Izumi, Masahiro Takebe, Takashi Kawabe, and Yusuke Tateishi for their advice on the needs and trends of mass spectrometers. Our appreciation extends to Makoto Yamada and Tomoyuki Oshiro for obtaining the single response waveforms of the detectors, to Masaru Nishiguchi and Hiroyuki Miura for providing samples, to Momoka Hayashida and Noriko Kato for offering insights on nucleic acid analysis, to Yuta Miyazaki for advice on noise characteristics, and to Tomoya Kudo for guidance on ion optics simulations.

Conflict of interest

Authors TT, YN, KT, and JI were employed by Shimadzu Corporation.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/frans.2023.1301602/full#supplementary-material

References

Summary

Keywords

LC-MS, ESI, chemometrics, Bayesian inference, deconvolution, signal processing, nucleic-acid-drugs

Citation

Tomono T, Hara S, Nakai Y, Takahara K, Iida J and Washio T (2024) A Bayesian approach for constituent estimation in nucleic acid mixture models. Front. Anal. Sci. 3:1301602. doi: 10.3389/frans.2023.1301602

Received

05 September 2023

Accepted

14 December 2023

Published

08 January 2024

Volume

3 - 2023

Edited by

Krzysztof Bernard Bec, University of Innsbruck, Austria

Reviewed by

Rui Vitorino, University of Aveiro, Portugal

Marina De Gea Neves, University of Duisburg-Essen, Germany

Updates

Copyright

*Correspondence: Taichi Tomono,

Disclaimer

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.

Outline

Figures

Cite article

Copy to clipboard


Export citation file


Share article

Article metrics