A method for gas chromatography-mass spectrometry data deconvolution

Abstract

The application discloses a gas chromatography-mass spectrometry data deconvolution method, and relates to the technical field of gas chromatography-mass spectrometry data analysis. The method preliminarily decomposes a time-mass-to-charge ratio two-dimensional intensity matrix by using a non-negative matrix decomposition function, obtains an initial elution curve matrix and an initial mass spectrum matrix, and adopts an exponential modified Gaussian peak shape modeling function to fit an initial peak shape parameter set. For each component number: according to the initial peak shape parameter set, the exponential modified Gaussian peak shape modeling function is used to generate an original elution curve matrix; the non-negative least square optimization and the nonlinear least square optimization method are used to cyclically optimize the mass spectrum matrix, the peak shape parameter set and the elution curve matrix; and thus the optimal component number is selected by using a comprehensive score function, and the optimal peak shape parameter set, the mass spectrum matrix and the elution curve matrix under the optimal component number are sequentially output. The application can realize accurate separation of co-flow-out components and improve the accuracy and reliability of complex mixture analysis.

Description

Technical Field

[0001] This application relates to the field of gas chromatography-mass spectrometry (GC-MS) data analysis technology, and in particular to a method for deconvolution of GC-MS data. Background Technology

[0002] Gas chromatography-mass spectrometry (GC-MS) technology, with its high separation efficiency and high detection sensitivity, has become a core method for the qualitative and quantitative analysis of complex mixtures. It is widely used in environmental pollutant detection (such as the analysis of polycyclic aromatic hydrocarbons in water), drug component quality control (such as the screening of impurities in active pharmaceutical ingredients), biosample metabolomics research (such as the identification of endogenous metabolites in blood), detection of pesticide and veterinary drug residues in food, and characterization of petrochemical product components.

[0003] In practical GC-MS analysis, limitations imposed by the separation capabilities of the stationary phase of the chromatographic column, the complexity of the sample matrix (such as a large number of matrix interfering substances in biological samples), and the differences in component elution behavior can lead to partial or complete co-elution of different components during column elution. This results in the superposition of acquired two-dimensional mass spectrometry signals (retention time-mass-charge ratio), forming complex mixed peaks. This peak overlap problem causes two major technical bottlenecks: firstly, the characteristic mass spectrometry signals of trace components are masked by the signals of high-concentration components, leading to the missed detection of trace components; secondly, the superimposed mixed peak signals cause distortion in mass spectral library similarity matching, and at the same time, it is impossible to accurately divide the peak areas of each component during integration quantification. Ultimately, this leads to a decrease in the accuracy of component identification and significant deviations in quantitative results, severely restricting the improvement of the efficiency of GC-MS technology in the analysis of complex systems.

[0004] To address the issue of overlapping chromatographic peaks, existing mainstream deconvolution methods can be broadly categorized into two types:

[0005] Peak shape model-based fitting methods. These methods describe chromatographic peak shapes using pre-defined mathematical models (such as Gaussian models or Exponentially Modified Gaussian (EMG) models) and employ optimization techniques such as nonlinear least squares to estimate peak positions, peak widths, and other parameters, thereby achieving the separation of mixed signals. While these methods have clear physical meaning, their performance is highly dependent on the reasonableness of the peak shape assumptions. When multiple co-elution peaks exist in the sample or matrix interference causes peak shape distortion, the fitting accuracy decreases significantly, resulting in insufficient anti-interference capability and making it difficult to meet the analytical needs of complex systems.

[0006] Blind source separation methods based on matrix factorization. These methods include Principal Component Analysis (PCA), Independent Component Analysis (ICA), and Non-negative Matrix Factorization (NMF). They do not rely on pre-defined peak shapes but extract characteristic components from mixed signals through mathematical decomposition, offering data-driven flexibility. NMF, due to the introduction of non-negativity constraints, better reflects the non-negativity of chromatographic intensity and mass spectrometry signals compared to PCA and ICA. However, its decomposition results may still deviate from actual chromatographic elution patterns or lack corresponding mass spectrometry characteristics of real compounds, thus limiting chemical interpretability. Furthermore, these methods are sensitive to noise and prone to generating meaningless spurious components. It is important to note that in a typical NMF model, the parameter scale is approximately k×(t+m), where t is the number of time points, m is the mass-to-charge ratio channel number, and k is the number of decomposed components. This high degree of freedom, while explaining more variance, also makes it prone to overfitting noise, leading to inaccurate signals in the extraction results and affecting the reliability of qualitative and quantitative analysis.

[0007] From the perspective of variance decomposition, total variance can be decomposed into explained variance and residual variance. For the same input data... X A higher explained variance usually means a better model fit, but the risk of overfitting must be taken into account.

[0008] In terms of model degrees of freedom, a typical NMF model includes (Chromatography curve matrix) and (Mass spectrum matrix), the total number of parameters is approximately The higher degree of freedom inherent in this approach carries a risk of overfitting. In contrast, peak modeling methods based on EMG require fitting 3-4 peak shape parameters (such as mean μ, standard deviation σ, and attenuation factor) for each peak. and optional amplitude A), plus There are several mass spectrometry intensity parameters, with a total of approximately [number] parameters. Its degrees of freedom are relatively low, resulting in a lower risk of overfitting. The difference in degrees of freedom directly determines the fitting flexibility of the two types of methods: NMF can explain more variance but is more prone to fitting noise; EMG has a more robust fit but lacks flexibility.

[0009] In summary, existing deconvolution methods struggle to balance the fit of physical models with the data-driven advantages of matrix factorization, and their ability to handle severe co-elution in complex systems is limited. Therefore, there is an urgent need to develop a novel GC-MS data deconvolution strategy that can achieve precise separation of co-eluting components while preserving the physical characteristics of chromatographic peaks and the chemical meaning of mass spectrometry signals, thereby improving the accuracy and reliability of complex mixture analysis. Summary of the Invention

[0010] The purpose of this application is to provide a method for deconvolution of gas chromatography-mass spectrometry data, which can achieve accurate separation of co-eluting components and improve the accuracy and reliability of complex mixture analysis.

[0011] To achieve the above objectives, this application provides the following solution:

[0012] This application provides a method for deconvolution of gas chromatography-mass spectrometry data, including:

[0013] Based on the raw gas chromatography-mass spectrometry data, a two-dimensional intensity matrix of time-mass-charge ratio was constructed.

[0014] Multiple groups are determined, and based on these groups, the time-mass-charge ratio two-dimensional intensity matrix is ​​initially decomposed using a non-negative matrix factorization function to obtain the initial elution curve matrix and the initial mass spectrum matrix.

[0015] Based on the initial elution curve matrix, an exponentially modified Gaussian peak modeling function is used to fit the initial peak parameter set;

[0016] Based on the group score and the initial peak shape parameter set, the original elution curve matrix is ​​generated by using the exponentially modified Gaussian peak shape modeling function.

[0017] Based on the original elution curve matrix, the initial mass spectrometry matrix is ​​optimized using non-negative least squares to obtain the mass spectrometry matrix for the current iteration in each group.

[0018] Based on the mass spectrometry matrix of the current iteration for each group score, the initial peak shape parameter set is optimized using nonlinear least squares to obtain the optimized peak shape parameter set for the current iteration. The optimized peak shape parameter set replaces the initial peak shape parameter set, and the process returns to the step "Based on each group score and the initial peak shape parameter set, the original elution curve matrix is ​​generated using the exponentially modified Gaussian peak shape modeling function". This process continues until the iteration termination condition is met, thus obtaining the optimal peak shape parameter set, mass spectrometry matrix, and elution curve matrix for each group score.

[0019] Based on the optimal set of peak shape parameters, mass spectrometry matrix, and elution curve matrix for each group score, a comprehensive scoring function is used to determine the score value for each group score; the comprehensive scoring function includes reconstruction error and complexity penalty terms.

[0020] The optimal component score is selected based on the minimum score value. The optimal set of peak shape parameters, mass spectrometry matrix, and elution curve matrix under the optimal component score are then arranged and output according to the peak position of the elution curve.

[0021] According to the specific embodiments provided in this application, this application has the following technical effects:

[0022] This application provides a deconvolution method for gas chromatography-mass spectrometry (GC-MS) data. In each iteration, the mass spectrometry matrix is ​​first updated using a non-negative least squares method, followed by an update of the peak shape parameter set using nonlinear least squares optimization, thereby generating a new elution curve matrix. This process, while preserving the physical meaning of the exponentially modified Gaussian model, fully leverages the data adaptability of matrix decomposition to achieve accurate separation of co-eluting components. To prevent model overfitting, the comprehensive scoring function includes reconstruction error and complexity penalty terms, introducing a penalty term based on reconstruction error and complexity into the algorithm. The value selection mechanism seeks the optimal balance between goodness of fit and model complexity, thereby automatically determining the number of components. This ensures that while pursuing high explained variance, the algorithm effectively suppresses the risk of overfitting, providing accurate and reliable deconvolution results for the analysis of complex systems. Attached Figure Description

[0023] To more clearly illustrate the technical solutions in the embodiments of this application or related technologies, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0024] Figure 1 A schematic flowchart of a gas chromatography-mass spectrometry data deconvolution method provided in this application embodiment;

[0025] Figure 2 A simplified flowchart illustrating a gas chromatography-mass spectrometry data deconvolution method provided in this application embodiment;

[0026] Figure 3 This is a schematic diagram illustrating the principle of the gas chromatography-mass spectrometry data deconvolution method provided in the embodiments of this application. Detailed Implementation

[0027] The technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, and not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0028] To make the above-mentioned objectives, features and advantages of this application more apparent and understandable, the application will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0029] Existing deconvolution methods for GC-MS data have the following shortcomings:

[0030] 1. Incomplete peak separation: Conventional Gaussian models are unable to accurately describe chromatographic peak tailing, making it difficult to effectively distinguish co-eluted components;

[0031] 2. Trace components are difficult to detect: weak signals are easily masked by strong peaks, resulting in missed detection or mismatch.

[0032] 3. Lack of standards for model selection: Traditional methods struggle to determine the optimal number of groups, making them prone to overfitting or underfitting.

[0033] 4. Insufficient interpretability of results: It relies solely on data-driven decomposition, ignores the physical laws of chromatographic peaks, and lacks chemical significance.

[0034] Therefore, there is a need for a deconvolution method that can automatically select components in the retention time-mass-charge ratio two-dimensional space, taking into account both physical constraints and the advantages of data-driven approaches, in order to improve the accuracy and reliability of GC-MS analysis of complex systems.

[0035] In view of this, in an exemplary embodiment, such as Figure 1 As shown, a method for deconvolution of gas chromatography-mass spectrometry (GC-MS) data is provided. This method is executed by a computer device, specifically by a terminal or server alone, or by both a terminal and a server. In this embodiment, it includes steps 101 to 108. Wherein:

[0036] Step 101: Construct a two-dimensional intensity matrix of time-mass ratio based on the original gas chromatography-mass spectrometry data.

[0037] Step 102: Determine multiple groups and, based on the multiple groups, use a non-negative matrix decomposition function to initially decompose the time-mass-charge ratio two-dimensional intensity matrix to obtain the initial elution curve matrix and the initial mass spectrum matrix.

[0038] Step 103: Based on the initial elution curve matrix, use the exponentially modified Gaussian peak modeling function to fit the initial peak parameter set.

[0039] Step 104: Based on the scores of each group and the initial set of peak shape parameters, use the exponentially modified Gaussian peak shape modeling function to generate the original elution curve matrix.

[0040] Step 105: Based on the original elution curve matrix, optimize the initial mass spectrometry matrix using non-negative least squares to obtain the mass spectrometry matrix for the current iteration in each group.

[0041] Step 106: Based on the mass spectrometry matrix of the current iteration for each group number, use nonlinear least squares to optimize the initial peak shape parameter set, obtain the optimized peak shape parameter set for the current iteration, replace the initial peak shape parameter set with the optimized peak shape parameter set, return to step 104, until the iteration termination condition is met, and obtain the optimal peak shape parameter set, mass spectrometry matrix and elution curve matrix for each group number.

[0042] Step 107: Based on the optimal set of peak shape parameters, mass spectrometry matrix, and elution curve matrix for each group score, determine the score value for each group score using a comprehensive scoring function; the comprehensive scoring function includes reconstruction error and complexity penalty terms.

[0043] Step 108: Select the group score corresponding to the minimum score as the optimal group score, and then arrange and output the optimal set of peak shape parameters, mass spectrometry matrix and elution curve matrix under the optimal group score according to the peak position of the elution curve.

[0044] Implementing steps 101 to 108 above provides a GC-MS data deconvolution method based on joint iterative optimization of chromatographic physical models and matrix decomposition, with automatic selection of component numbers. This method achieves accurate separation of co-eluting components in the two-dimensional space of retention time and mass-to-charge ratio. It fully utilizes the physical constraints of chromatographic peak shapes to ensure the authenticity of results, while leveraging the data-driven advantages of matrix decomposition to improve separation flexibility. Ultimately, it solves the problems of trace component signal masking, mass spectrometry similarity matching distortion, and integral quantification deviation, providing a reliable technical solution for GC-MS analysis of complex systems.

[0045] like Figure 2 As shown, the method of this application includes two stages:

[0046] Phase 1: Data Preprocessing and Initial Decomposition

[0047] The raw GC-MS data is preprocessed to construct a time-to-mass-to-charge ratio two-dimensional intensity matrix. A non-negative matrix factorization function is called to perform preliminary matrix decomposition, outputting the initial elution profile matrix and initial mass spectrum matrix of each component, thus determining the benchmark for iterative optimization. Simultaneously, the data dimensionality range is defined (number of retention time points t, number of mass-to-charge ratio (m / z) intervals m, and optimal component number). ).

[0048] Phase Two: Iterative Optimization and Optimal Model Selection

[0049] Based on the initial elution profile matrix, initial mass spectrometry matrix, and initial peak shape parameters obtained in the first stage, the process involves: updating the elution curve matrix C with θ → updating the mass spectrometry matrix S using non-negative least squares based on C → fixing S and updating the peak shape parameters using least squares fitting. → The iterative process of “calculating the relative reconstruction error and determining the convergence of the iteration” ensures the non-negativity of the signal and the physical morphology of the chromatographic peaks, while suppressing overfitting and avoiding parameter redundancy through the penalty term, thus achieving synergistic optimization.

[0050] The method of this application can automate data processing through computer programs. The following details the technical process, clarifying the operational logic and technical details of each step to ensure that the method is repeatable and scalable.

[0051] In another exemplary embodiment of this application, step 101 above involves data input and preprocessing, and step 101 can be replaced by the following steps 201 to 207:

[0052] Step 201: Extract mass spectrometry data from the colorimetric-mass linkage data file.

[0053] Data extraction: The mass spectrometry data parsing module (supporting mainstream formats such as mzXML) is called by a computer program to read the GC-MS data file and extract key information from the GC-MS raw data file - MS1 (first-level mass spectrometry) scan data.

[0054] Step 202: Determine the retention time points corresponding to the preset number of scan points.

[0055] Retention Time Series: Extract the retention time corresponding to a preset number (e.g., 200) of scan points to form a sequence. That is, the length of the time dimension.

[0056] Step 203: Obtain all original mass-to-charge ratios and the corresponding net intensities for each original mass-to-charge ratio from the mass spectrometry data.

[0057] Original mass-to-charge ratio and intensity: for each scan point This function extracts all original mass-to-charge ratios (m / z) and corresponding signal intensities at that time point. Therefore, it returns two arrays: mass-to-charge ratio (m / z) and intensity (Intensity). Each peak consists of both m / z and intensity components.

[0058] Step 204: Preset multiple mass-to-load ratio bin intervals.

[0059] Binning Parameter Settings: Preset m / z binning parameters according to analysis requirements: m / z range: Set the m / z interval for the target analysis (e.g., 35.0~400.0, unit: Da); Binning interval: Set the binning width (e.g., 0.5 or 0.1), dividing the continuous m / z range into m equal-width intervals (m is the mass-to-charge ratio dimension length), forming binning intervals. .

[0060] Step 205: For each scanning point, accumulate the net intensity corresponding to all original mass-to-charge ratios within each mass-to-charge ratio bin interval to obtain the intensity of each scanning point within each mass-to-charge ratio bin interval.

[0061] Box-by-box strength calculation: For each scanning point i, the corrected net strength is accumulated over the box-by-box interval in m / z: For the box-by-box interval The net intensity corresponding to all original m / z values ​​within the interval is summed to obtain the binned intensity; if there is no original m / z value within the interval, the intensity is set to 0 (to ensure matrix integrity).

[0062] The mass spectrometry data (mz and inten) are binned, and the cumulative intensity inten of the mass spectrometry data in each bin is returned. The center value of each bin is calculated as the mz value of that bin. This is to discretize the data within a predetermined interval to form a histogram.

[0063] Step 206: Construct an initial time-mass-charge ratio two-dimensional intensity matrix with the retained time points as rows, the mass-charge ratio bin intervals as columns, and the intensity of each scan point in each mass-charge ratio bin interval as elements.

[0064] By filtering out the low-quality (<35) noise region and the high-quality (>400) no-signal region, a two-dimensional intensity matrix of time-mass-charge ratio with dimension t×m is formed (where m is the number of m / z bin intervals), which reduces the data dimension while retaining the effective signal.

[0065] Initial time-mass-to-charge ratio two-dimensional intensity matrix :

[0066] ;

[0067] in, , representing the signal strength at the i-th retention time point and the j-th m / z sub-bin interval. Complete the data structuring process.

[0068] Step 207: Denoise the initial time-mass-charge ratio two-dimensional intensity matrix to obtain the final time-mass-charge ratio two-dimensional intensity matrix.

[0069] Savitzky-Golay smoothing, SNIP baseline correction, and weak signal filtering (such as intensity below 1% of the maximum value) are performed on each scan (row) to obtain the denoising matrix.

[0070] In another exemplary embodiment of this application, step 102 above involves NMF initial decomposition and parameter initialization, and step 102 can be replaced by the following steps 301 to 302:

[0071] Step 301: Determine the range of values ​​for the grouped fractions, and select multiple grouped fractions from the range of values ​​for the grouped fractions.

[0072] NMF decomposition parameter settings: Define the initial grouping range (e.g., ... to ), group score for each attempt Perform NMF decomposition. Determine the possible range of the number of components k to avoid overfitting or underfitting in subsequent optimization, and provide a starting baseline that conforms to the data characteristics for iterative optimization.

[0073] Step 302: Based on the multiple component numbers, call the non-negative matrix decomposition function to perform a preliminary decomposition of the time-mass-charge ratio two-dimensional intensity matrix with the goal of minimizing the decomposition error, obtaining the initial elution curve matrix and the initial mass spectrum matrix; the expression for the decomposition error is: ;in, The time-mass-charge ratio two-dimensional intensity matrix, This is the initial elution curve matrix (also known as the initial elution profile matrix, which reflects the elution trend of the components over time). This is the initial mass spectrum matrix (reflecting the intensity distribution of components in the m / z dimension). It is the Frobenius norm.

[0074] Decompose the observation matrix X into the elution curve matrix. With mass spectrum matrix The product of , k is the number of components, and the peak parameters of each component are fitted simultaneously. This allows it to interpret each peak in the data.

[0075] The initial matrices C and S are generated by a nonnegative matrix factorization multiplication update algorithm:

[0076] .

[0077] in This indicates element-wise multiplication. For a small constant, to prevent division by zero, X≈WH, and As the initial C, H0 is the initial S.

[0078] Set NMF iteration parameters: maximum number of iterations, convergence tolerance, to avoid excessive time consumption or insufficient accuracy in decomposition.

[0079] Initial matrix output: Two core initial matrices are obtained after decomposition. Initial elution profile matrix. Each column Represents the elution intensity change trend of the i-th component with retention time; initial mass spectrometry matrix Each line This represents the signal intensity distribution of the i-th component in each m / z bin interval (all satisfying the non-negativity constraint). NMF is iteratively optimized using a multiplicative update rule, satisfying the nonnegativity constraint. The outputs W0 and H0 are used as initial values ​​for subsequent fine-tuning.

[0080] In another exemplary embodiment of this application, an EMG function is selected to describe the dynamic behavior of chromatographic peaks, specifically to describe asymmetric chromatographic peaks with steep leading edges and tailing characteristics. This function uses peak position (μ), peak width (σ), and tailing factor (σ) to describe these peaks. Three parameters are used to accurately characterize the asymmetric features of actual chromatographic peaks, characterized by a steep leading edge and a trailing tail (solving the problem that Gaussian models cannot describe tailed peaks). Initial peak shape parameters for each component are fitted based on the initial elution profile matrix, and the maximum number of EMG iterations (e.g., 30 times) and convergence threshold (e.g., 1e-5) are set to ensure parameter convergence. The exponentially modified Gaussian (EMG) function models the chromatographic peaks, with parameters... ,formula:

[0081] ;

[0082] Where μ is the peak position parameter, used to represent the central retention time of the component elution peak; σ is the peak width parameter, used to represent the degree of peak expansion; – Tail parameter, used to indicate the degree of asymmetric tailing of the peak shape.

[0083] .

[0084] It is the elution curve matrix generated by the EMG function.

[0085] The EMG function is used to simulate the asymmetric shape of chromatographic peaks, through parameters. To adjust the tailing degree, the exponential term in the formula needs to be clipped to [-700, 700] to avoid numerical overflow and accurately reflect the tailing phenomenon of the peak elution curve over retention time. Step 103 above involves peak shape function and deconvolution: initial peak shape parameter fitting: based on the initial elution profile matrix. Call the exponentially modified Gaussian function for each column. (The elution curves of individual components) were initially fitted.

[0086] Extracting the initial peak position: The retention time corresponding to the maximum value is used as the initial value of the peak.

[0087] Peak Time : The retention time corresponding to the maximum value. .

[0088] Optional initial guess ( ), It is adapted to the elution characteristics of capillary columns (peak width is usually 1 second to 10 seconds) and can be modified according to the actual situation of different types of columns.

[0089] Fitting peak shape parameters: The parameters of the Gaussian function (such as peak center, peak width, tailing coefficient, etc.) are optimized using the least squares method to make the fitted curve closely resemble the peak shape. Minimum error;

[0090] Fitting EMG to... using least squares method : ,in For scaling factor (matching) (Intensity range), optional .

[0091] Forming an initial set of peak shape parameters This provides a starting point for iterative optimization.

[0092] In another exemplary embodiment of this application, steps 104 to 106 constitute an alternating iterative optimization strategy, which involves cyclic iteration of the groupings and cyclic iteration under each grouping. Alternating iterative optimization strategy: sequentially taking k from... arrive Different component sets are optimized using a three-step iterative process: generating elution curves, optimizing the mass spectrometry matrix, and optimizing peak shape parameters. This process achieves convergence of the objective function while simultaneously considering physical model constraints and data-driven optimization. The specific workflow is as follows:

[0093] Parameter basis: EMG peak shape parameters of the current iteration round For input;

[0094] Peak shape generation: Based on the current peak shape parameters, generate the elution curve matrix C:

[0095] Calculate the original elution curves of each component using the EMG function. During the calculation, the exponent term needs to be clipped to [-700, 700]; the peak area of ​​the original elution curve of each component is calculated by trapezoidal integration. Then normalize the original curve: To ensure that the peak areas of the elution curves of each component are consistent, eliminate the interference of intensity redundancy on subsequent optimization, and avoid zero-area peaks. A small positive number is added to the denominator to prevent division by zero. For the first The original elution curves of each component, For the first The original elution peak shapes of each component, For the first Peak area of ​​the original elution peak shape of each component.

[0096] fixed The mass spectrum matrix S is optimized using Non-Negative Least Squares (NNLS) to ensure that the mass spectrum intensity is non-negative.

[0097] Solving column by column: For each mass-to-charge ratio bin j (j=1,2,...,m), solve for the mass spectral intensity vector S[:,j] under non-negativity constraints. The mathematical expression is:

[0098] .

[0099] In the formula, For the first The mass spectrum intensity vector of each mass-to-charge ratio bin. The second time-mass-charge ratio two-dimensional intensity matrix is ​​the first... The column data of the mass-to-charge ratio bins, This is the original elution curve matrix. For the first Optimization results of mass spectrometry intensity for each mass-charge ratio bin The mass-to-load ratio is the number of boxes.

[0100] Principle: The strength of each m / z sub-bin is the sum of the contributions of each component. NNLS ensures that the contribution of each component is non-negative (consistent with the quantitative logic of chemistry).

[0101] Matrix integration: Optimization results of all mass-load bins Integrate to form the complete mass spectrometry matrix S for the current iteration round.

[0102] With S fixed, nonlinear least squares optimization of peak shape parameters is performed. :

[0103] Residual focusing: With the mass spectrometry matrix S fixed, calculate the "residual signal of target component i after removing contributions from other components" to eliminate interference from signal superposition between components. The formula is:

[0104] .

[0105] Parameter optimization: with the goal of "minimizing" and Optimize the EMG peak shape parameters of component i with the residual F-norm as the objective. The mathematical expression is:

[0106] .

[0107] in, For the first Peak shape parameters of each component, For parameter boundaries, Indicates the first The peak shape parameters of each component are The corresponding elution curve at time, The first in the mass spectrum matrix Signal intensity distribution of each component in each mass-to-charge ratio bin interval For the first The residual signal of each component , The time-mass-charge ratio two-dimensional intensity matrix, For the first The original elution curves of each component, The first in the mass spectrum matrix Signal intensity distribution of each component in each mass-to-charge ratio bin interval. The parameter boundary bounds must satisfy... >0、 For values ​​>0, the optimization tool can be a nonlinear algorithm that supports boundary constraints.

[0108] Convergence criterion and loop: Calculate the relative reconstruction error of the current iteration:

[0109] .

[0110] If satisfied ( For the convergence threshold, such as If the iteration stops when the number of iterations reaches max_iter (e.g., 20-40 times), the iteration of the current k is terminated; otherwise, the process returns to the "Generate elution curve matrix C based on the current peak shape parameters" step and starts the next round of iteration. In addition to convergence determination, an absolute threshold for reconstruction error can be set. If the error is too large, even if convergence is achieved, it will not be considered the final result.

[0111] Repeat the above three steps until the relative error change of reconstruction in two consecutive iterations reaches the threshold or the number of iterations reaches the upper limit. Save the current joint objective function value, complete the parameter optimization under the current component number, and realize the two-way collaboration of "physical model constraining chromatographic behavior and matrix decomposition optimizing signal allocation".

[0112] As can be seen from the above, step 104 can be replaced by steps 401 to 404, for example:

[0113] Step 401: Based on the fractions of each group and the initial peak shape parameter set, use the exponentially modified Gaussian peak shape modeling function to generate the original elution peak shape of each component;

[0114] Step 402: Calculate the peak area of ​​the original elution peak shape for each component;

[0115] Step 403: Based on the peak area, use the formula The original elution peak shapes of each component were normalized to obtain the original elution curves of each component.

[0116] Step 404: Construct a matrix of original elution curves for all components.

[0117] For example, step 105 above can be replaced by steps 501 to 502:

[0118] Step 501: For each mass-to-charge ratio bin, based on the original elution curve matrix, optimize the initial mass spectrum matrix using non-negative least squares, and apply the formula... The mass spectrometry intensity optimization results for each mass-to-charge ratio bin under non-negative constraints are obtained.

[0119] Step 502: Integrate the mass spectrometry intensity optimization results of all mass-charge ratio bins to form the mass spectrometry matrix for the current iteration round.

[0120] For example, step 106 above can be replaced by steps 601 to 602:

[0121] Step 601: The optimization function for establishing the peak shape parameter set is: ;

[0122] Step 602: Based on the optimization function of the peak shape parameter set and the mass spectrum matrix of the current iteration under each group, the initial peak shape parameter set is optimized using nonlinear least squares to obtain the optimized peak shape parameter set for the current iteration.

[0123] For example, the iteration termination conditions include: , or the number of iterations equals the maximum number of iterations. Where, To reconstruct the relative error, The time-mass-charge ratio two-dimensional intensity matrix, For the set of peak shape parameters The corresponding elution curve matrix, For mass spectrometry matrix, For the number of iterations, This represents the relative reconstruction error from the previous iteration. This represents the relative reconstruction error for the current iteration. The convergence threshold, As a precision requirement / convergence criterion, 1e-4 can be taken. For minimum iteration threshold control, Used for numerical protection to prevent division by zero, such as taking 1e-9.

[0124] In another exemplary embodiment of this application, a multivariate collaborative optimization framework is constructed to achieve global convergence of the objective function through alternating iterations.

[0125] Let W0 be the initial C and H0 be the initial S. That is, W0 corresponds to the initial elution curve matrix (which becomes C in subsequent iterations) and H0 corresponds to the initial mass spectrometry matrix (which becomes S in subsequent iterations).

[0126] Objective function definition: With minimizing the difference between the observed matrix and the reconstructed matrix as the core objective, the mathematical expression of the comprehensive scoring function is:

[0127] ;

[0128] To avoid overfitting due to excessively large group numbers, a penalty term is introduced that increases with the group number k. Wherein: This is the penalty coefficient (which can be adjusted according to data complexity to balance fitting accuracy and model simplicity; for example, 0.1 can be selected). m +3) represents the number of independent parameters for a single component, where 3 corresponds to peak position, peak width, and tailing parameters; log(size(X)) is the logarithm of the data size ( (i.e., the total number of elements in matrix X), allowing the penalty term to adapt to data of different dimensions. Constraints: . Based on peak shape parameter matrix Constructed elution profile matrix; This is the optimized mass spectrometry matrix. It's worth noting that the comprehensive scoring function is calculated based on the time-mass-charge ratio two-dimensional intensity matrix after global intensity normalization.

[0129] For the component range set in the first stage Compare the objective function values ​​corresponding to the scores of all candidate groups. ;

[0130] ;

[0131] In the formula, For the first The score value under each group score The time-mass-charge ratio two-dimensional intensity matrix, For the first The set of peak shape parameters under each group score for The corresponding elution curve matrix, For the first Mass spectrum matrix for each group fraction It is the Frobenius norm. The penalty coefficient is... To retain time points.

[0132] The component with the minimum objective function value is selected as the optimal solution, i.e.:

[0133] ;

[0134] Finally, the model with the lowest score is selected as the output, thus filtering out the components that achieve the optimal balance between fitting accuracy and model simplicity. This achieves a globally optimal balance between "fitting accuracy and model complexity".

[0135] As can be seen, the purpose of the iterative updates and joint optimization (EMG parameter fitting, S matrix NNLS optimization) in this application is to combine the initial value of NMF with EMG peak shape constraints, and alternately optimize C (elution curve), S (mass spectrometry) and (EMG parameters), minimize the objective function value .

[0136] In another exemplary embodiment of this application, the components are sorted in ascending order based on the peak position (apex) of the elution curve, and the C, S, and The order of columns / elements.

[0137] First, calculate the peak of the elution peak for each component (i.e., the retention time corresponding to the point with the strongest intensity in the elution curve).

[0138] Then, the components are sorted according to the retention time of these vertices to obtain the sorting index "order".

[0139] Finally, the elution curve matrix C, mass spectrometry matrix S, and peak parameter list theta are rearranged according to this sorting index so that the components are arranged in elution order (retention time from smallest to largest).

[0140] Then step 108 above can be replaced by the following steps 701 to 703:

[0141] Step 701: In the optimal elution curve matrix, determine the retention time point corresponding to the point with the greatest intensity on each elution curve, and use it as the retention time point of the elution peak of each component; one elution curve corresponds to one component.

[0142] Step 702: Sort the components according to the retention time of the elution peak of each component in ascending order of retention time to obtain the sorting index.

[0143] Step 703: According to the sorting index, rearrange the optimal set of peak shape parameters, mass spectrometry matrix and elution curve matrix so that the components are arranged in the elution order.

[0144] In another exemplary embodiment of this application, in order to verify the effectiveness of the method of this application, after step 108 above, the method may further include steps 801 to 802. Wherein:

[0145] Step 801: Obtain the reconstructed time-mass-charge ratio two-dimensional intensity matrix based on the optimal mass spectrum matrix and elution curve matrix under the optimal component.

[0146] Step 802: Using the explained variance, signal-to-noise ratio, and spectral similarity as evaluation metrics, evaluate the optimal mass spectrometry matrix and elution curve matrix under the optimal grouping.

[0147] Explanation of variance proportions:

[0148] ;

[0149] ;

[0150] in, This represents the signal strength at the i-th retention time point and the j-th mass-to-charge ratio bin interval.

[0151] Signal-to-noise ratio:

[0152] ;

[0153] Wherein, the original data variance Var(X): calculates the global variance based on all elements of the two-dimensional intensity matrix X; the residual variance Var(X-CS): calculates the global variance based on all elements of the residual matrix.

[0154] Spectral similarity:

[0155] ;

[0156] in, Represents the true spectrum (such as the mass spectrum of a known standard substance, the true signal matrix of chromatographic peaks, etc.); This represents the estimated / reconstructed spectrum (such as the target spectrum obtained through deconvolution, model fitting, etc.); the similarity is mapped to the [-1,1] interval by "dot product divided by the product of their respective norms" (usually in [0,1] if the spectrum elements are non-negative), and the closer the value is to 1, the more similar the two spectra are.

[0157] Optimal Results Processing and Output: Components are sorted according to the retention time of the elution peak apex. Peak parameters (component + mu / sigma / tau), mass spectrometry matrix (component + mass-to-charge ratio intensity), elution curves (retention time + relative intensity of each component), reconstruction matrix, and model indices (k-value + score + explained variance + signal-to-noise ratio) are saved as a CSV file for convenient subsequent analysis. Total Ion Current (TIC), elution curves, mass spectra, and reconstruction comparison plots are plotted to visually display the deconvolution results. Key deconvolution results, such as data source, optimal component count, explained variance ratio, and signal-to-noise ratio, are output to the console.

[0158] Output a CSV file and visualization charts, including:

[0159] 1. Peak shape parameter matrix ; ;

[0160] 2. Pure component mass spectrometry matrix (Can be directly used for mass spectrometry library matching);

[0161] 3. Elution curve matrix ;

[0162] 4. Reconstructing the matrix This is used to verify the separation effect;

[0163] 5. Model evaluation metrics (explained variance EV, signal-to-noise ratio SNR, similarity).

[0164] Generate PNG charts: elution curve overlay, mass spectrum (first 20 m / z), preprocessing comparison chart (raw data vs. smoothed data vs. baseline-corrected data), k-value and index relationship chart, and reconstructed comparison chart, providing direct basis for subsequent qualitative identification and quantitative analysis.

[0165] The overall process of this application is as follows: First, the raw GC-MS data is preprocessed. A two-dimensional intensity matrix of time-mass-charge ratio is constructed using a mass-charge ratio bin. A non-negative matrix factorization function is called to initially decompose the matrix, outputting the initial elution profile matrix and the initial mass spectrum matrix. This determines the starting baseline for iterative optimization and clarifies the initial data dimension range (including the number of retention time points, mass-charge ratio interval division, and the initial component trial range). Second, the exponentially modified Gaussian peak shape modeling function is called to fit the peak shape parameters (peak position μ, peak width σ, and tailing factor) of each component based on the initial elution profile. This invention constructs a chromatographic physical constraint model; then, with peak shape parameters fixed, it updates the pure component mass spectrometry matrix channel by channel using a non-negative least squares algorithm; next, with the mass spectrometry matrix fixed, it optimizes the peak shape parameters of each component using nonlinear least squares, outputting the peak shape parameter matrix and mass spectrometry intensity matrix under the current iteration, while calculating the reconstruction error; then, it moves to the next iteration, repeating the steps of "fitting peak shape parameters - updating mass spectrometry matrix - calculating error" until the reconstruction error is less than a set threshold or the maximum number of iterations is reached; finally, it selects the optimal value from the attempted component scores using a comprehensive scoring function (reconstruction error + complexity penalty term); after sorting by elution peak vertex time, it outputs the final peak shape parameters, pure mass spectrometry matrix, elution curve matrix, and reconstructed data matrix. This invention combines the realistic fitting of physical models with the data-driven advantages of matrix factorization, enabling accurate separation of co-eluting components, significantly improving the accuracy and quantitative reliability of trace component identification in complex systems, and is suitable for GC-MS applications in environmental monitoring, drug analysis, metabolomics, food safety, and petrochemicals.

[0166] The working principle of the method in this application is as follows: Figure 3 As shown, this application constructs a fusion framework of chromatographic physical model and matrix decomposition through a closed-loop process of "data input → data preprocessing → peak shape function and deconvolution → initialization → iterative optimization → model selection → component sorting → result output".

[0167] The beneficial effects of this application are as follows:

[0168] 1. Overcoming the shortcomings of existing methods: By introducing joint optimization of EMG physical model and NMF data decomposition, both physical realism and data-driven flexibility are taken into account;

[0169] 2. High separation accuracy and reliability: NNLS is used to ensure non-negativity, and the number of components is screened by combining reconstruction error and complexity penalty terms to avoid the generation of false components, which can accurately separate severely co-efferent components;

[0170] 3. Strong detection capability for trace components: Weak signals are preserved under physical constraints, solving the problems of missed detection of trace components caused by peak overlap and mass spectrometry matching distortion;

[0171] 4. Wide applicability: The method has a clear process, adjustable parameters, and the output results (elution curve, pure mass spectrometry matrix, peak shape parameters) have clear chemical meaning. It is suitable for core GC-MS applications such as environmental monitoring, drug analysis, metabolomics, food safety, and petrochemicals.

[0172] The technical features of the above embodiments can be combined in any way. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0173] This document uses specific examples to illustrate the principles and implementation methods of this application. The descriptions of the above embodiments are only for the purpose of helping to understand the methods and core ideas of this application. Furthermore, those skilled in the art will recognize that, based on the ideas of this application, there will be changes in the specific implementation methods and application scope. Therefore, the content of this specification should not be construed as a limitation of this application.

Claims

1. A method for deconvolution of gas chromatography-mass spectrometry data, characterized in that, include: Based on the raw gas chromatography-mass spectrometry data, a two-dimensional intensity matrix of time-mass-charge ratio was constructed. Multiple groups are determined, and based on these groups, the time-mass-charge ratio two-dimensional intensity matrix is ​​initially decomposed using a non-negative matrix factorization function to obtain the initial elution curve matrix and the initial mass spectrum matrix. Based on the initial elution curve matrix, an exponentially modified Gaussian peak modeling function is used to fit the initial peak parameter set; Based on the group score and the initial peak shape parameter set, the original elution curve matrix is ​​generated by using the exponentially modified Gaussian peak shape modeling function. Based on the original elution curve matrix, the initial mass spectrometry matrix is ​​optimized using non-negative least squares to obtain the mass spectrometry matrix for the current iteration in each group. Based on the mass spectrometry matrix of the current iteration for each group score, the initial peak shape parameter set is optimized using nonlinear least squares to obtain the optimized peak shape parameter set for the current iteration. The optimized peak shape parameter set replaces the initial peak shape parameter set, and the process returns to the step "Based on each group score and the initial peak shape parameter set, the original elution curve matrix is ​​generated using the exponentially modified Gaussian peak shape modeling function". This process continues until the iteration termination condition is met, thus obtaining the optimal peak shape parameter set, mass spectrometry matrix, and elution curve matrix for each group score. Based on the optimal set of peak shape parameters, mass spectrometry matrix, and elution curve matrix for each group score, a comprehensive scoring function is used to determine the score value for each group score; the comprehensive scoring function includes reconstruction error and complexity penalty terms. The optimal component score is selected based on the minimum score value. The optimal set of peak shape parameters, mass spectrometry matrix, and elution curve matrix under the optimal component score are then arranged and output according to the peak position of the elution curve.

2. The method for deconvolution of gas chromatography-mass spectrometry data according to claim 1, characterized in that, Based on the raw gas chromatography-mass spectrometry data, a two-dimensional intensity matrix of time-mass-charge ratio was constructed, specifically including: Extract mass spectrometry data from a colorimetric-mass-linked data file; Determine the retention time points corresponding to a preset number of scan points; Obtain all original mass-to-charge ratios and the corresponding net intensities at each scan point from the mass spectrometry data; Multiple mass-to-load ratio bin intervals are preset; For each scanning point, the net intensity corresponding to all original mass-to-charge ratios within each mass-to-charge ratio bin interval is accumulated to obtain the intensity of each scanning point within each mass-to-charge ratio bin interval; Using the retained time points as rows and the mass-to-charge ratio bin intervals as columns, the intensity of each scan point in each mass-to-charge ratio bin interval is used as an element to construct an initial time-mass-to-charge ratio two-dimensional intensity matrix; The initial time-mass-charge ratio two-dimensional intensity matrix is ​​denoised to obtain the final time-mass-charge ratio two-dimensional intensity matrix.

3. The method for deconvolution of gas chromatography-mass spectrometry data according to claim 1, characterized in that, Multiple component numbers are determined, and based on these component numbers, a nonnegative matrix factorization function is used to initially decompose the time-mass-charge ratio two-dimensional intensity matrix to obtain the initial elution curve matrix and the initial mass spectrum matrix. Specifically, this includes: Determine the range of values ​​for the grouped numbers, and select multiple grouped numbers from the range of values ​​for the grouped numbers; Based on multiple component numbers, a nonnegative matrix factorization function is invoked to perform a preliminary decomposition of the time-mass-charge ratio two-dimensional intensity matrix with the goal of minimizing the decomposition error, thereby obtaining the initial elution curve matrix and the initial mass spectrum matrix; the expression for the decomposition error is: ;in, The time-mass-charge ratio two-dimensional intensity matrix, This is the initial elution curve matrix. The initial mass spectrum matrix, It is the Frobenius norm.

4. The method for deconvolution of gas chromatography-mass spectrometry data according to claim 1, characterized in that, Based on the scores of each group and the initial set of peak shape parameters, an exponentially modified Gaussian peak shape modeling function is used to generate the original elution curve matrix, which specifically includes: Based on the fractions of each group and the initial set of peak shape parameters, the original elution peak shape of each component is generated using an exponentially modified Gaussian peak shape modeling function. Calculate the peak area of ​​the original elution peak shape for each component; Based on the peak area, using the formula The original elution peak shapes of each component are normalized to obtain the original elution curves of each component; where, For the first The original elution curves of each component, For the first The original elution peak shapes of each component, For the first Peak area of ​​the original elution peak shape of each component; The original elution curves of all components are used to construct an original elution curve matrix.

5. The method for deconvolution of gas chromatography-mass spectrometry data according to claim 1, characterized in that, Based on the original elution curve matrix, the initial mass spectrometry matrix is ​​optimized using non-negative least squares to obtain the mass spectrometry matrix for the current iteration in each group number, specifically including: For each mass-to-charge ratio bin, based on the original elution curve matrix, the initial mass spectrum matrix is ​​optimized using non-negative least squares, and the formula is applied. Solve for the mass spectral intensity optimization results for each mass-to-charge ratio bin under non-negative constraints; where, For the first The mass spectrum intensity vector of each mass-to-charge ratio bin. The second time-mass-charge ratio two-dimensional intensity matrix is ​​the first... The column data of the mass-to-charge ratio bins, This is the original elution curve matrix. For the first Optimization results of mass spectrometry intensity for each mass-charge ratio bin , The mass-to-load ratio is the number of boxes; The mass spectrometry intensity optimization results of all mass-charge ratio bins are integrated to form the mass spectrometry matrix for the current iteration.

6. The method for deconvolution of gas chromatography-mass spectrometry data according to claim 1, characterized in that, Based on the mass spectrum matrix of the current iteration for each group, the initial peak shape parameter set is optimized using nonlinear least squares to obtain the optimized peak shape parameter set for the current iteration, which specifically includes: The optimization function for establishing the peak shape parameter set is: ;in, For the first Peak shape parameters of each component, s is the parameter boundary, Indicates the first The peak shape parameters of each component are The corresponding elution curve at time, The first in the mass spectrum matrix Signal intensity distribution of each component in each mass-to-charge ratio bin interval It is the Frobenius norm. For the first The residual signal of each component , The time-mass-charge ratio two-dimensional intensity matrix, For the first The original elution curves of each component, The first in the mass spectrum matrix Signal intensity distribution of each component in each mass-to-charge ratio bin interval; Based on the optimization function of the peak shape parameter set and the mass spectrum matrix of the current iteration under each group, the initial peak shape parameter set is optimized by nonlinear least squares to obtain the optimized peak shape parameter set for the current iteration.

7. The method for deconvolution of gas chromatography-mass spectrometry data according to claim 1, characterized in that, The iteration termination conditions include: Or the number of iterations equals the maximum number of iterations; in, To reconstruct the relative error, , The time-mass-charge ratio two-dimensional intensity matrix, For the set of peak shape parameters The corresponding elution curve matrix, For mass spectrometry matrix, It is the Frobenius norm; For the number of iterations, This represents the relative reconstruction error from the previous iteration. This represents the relative reconstruction error for the current iteration. The convergence threshold, This is for minimum iteration threshold control.

8. The method for deconvolution of gas chromatography-mass spectrometry data according to claim 1, characterized in that, The comprehensive scoring function is: ； In the formula, For the first The score value under each group score The time-mass-charge ratio two-dimensional intensity matrix, For the first The set of peak shape parameters under each group score for The corresponding elution curve matrix, For the first Mass spectrum matrix for each group fraction It is the Frobenius norm. The penalty coefficient is... To preserve time points, The mass-to-load ratio is the number of boxes. The number of independent parameters for a single component.

9. The method for deconvolution of gas chromatography-mass spectrometry data according to claim 1, characterized in that, Based on the peak positions of the elution curves, the optimal set of peak shape parameters, mass spectrometry matrix, and elution curve matrix for the optimal component groups are arranged, specifically including: In the optimal elution curve matrix, the retention time point corresponding to the point with the greatest intensity on each elution curve is determined as the retention time point of the elution peak of each component; one elution curve corresponds to one component. Based on the retention time of the elution peak of each component, the components are sorted in ascending order of retention time to obtain a sorting index; According to the sorting index, the optimal set of peak shape parameters, mass spectrometry matrix, and elution curve matrix are rearranged so that the components are arranged in the elution order.

10. The method for deconvolution of gas chromatography-mass spectrometry data according to claim 1, characterized in that, The component score corresponding to the minimum score is selected as the optimal component score. Then, based on the peak position of the elution curve, the optimal set of peak shape parameters, mass spectrometry matrix, and elution curve matrix under the optimal component score are arranged and output. This is followed by: Based on the optimal mass spectrum matrix and elution curve matrix under the optimal grouping, the reconstructed time-mass-charge ratio two-dimensional intensity matrix is ​​obtained; The optimal mass spectrometry matrix and elution curve matrix under the optimal grouping are evaluated using explained variance, signal-to-noise ratio, and spectral similarity as evaluation metrics.

Images