Metabonomics data batch correction method based on multi-kernel learning

Abstract

The invention discloses a metabonomics data batch correction method based on multi-kernel learning, and belongs to the cross technical field of bioinformatics and analytical chemistry. According to the method, the multi-kernel learning technology is utilized, the advantages of different kernel functions are fused in a self-adaptive mode, a model conforming to data reality is constructed, complex drift characteristics of metabolite signals are accurately captured, and efficient and accurate normalization processing of metabonomics data is achieved. Compared with traditional data standardization methods such as SVR and LOESS, the method has the advantages that the performance is excellent in the aspect of reducing the metabolite peak intensity variability, and the data stability is remarkably improved. In the subsequent multivariate statistical analysis, the classification accuracy is greatly improved, the comparability among different batches of data is also remarkably enhanced, reliable data support can be provided for discovery of disease biomarkers, and the method plays a key role in large-scale metabonomics research.

Description

Technical Field

[0001] This application relates to the interdisciplinary fields of bioinformatics and analytical chemistry, and in particular to a batch correction method for metabolomics data based on multi-kernel learning. Background Technology

[0002] Metabolomics, as an important branch of systems biology, focuses on the comprehensive analysis of small-molecule metabolites in organisms, aiming to reveal the metabolic changes in biological systems under different physiological and pathological states. Its research findings have shown enormous application potential in many fields such as disease prediction, drug development, and environmental toxicology. High-resolution mass spectrometry techniques, such as liquid chromatography-mass spectrometry (LC-MS) and gas chromatography-mass spectrometry (GC-MS), have become the most commonly used data acquisition methods in metabolomics research due to their advantages of high sensitivity, high resolution, and wide dynamic range.

[0003] However, in large-scale metabolomics studies, batch effects have become a key issue affecting data quality and the reliability of research results. Due to the large number of samples, the analysis process often needs to be performed in batches. It is difficult to ensure that the experimental conditions between different batches are completely consistent, such as instrument status, mobile phase composition, ambient temperature, and humidity. These small differences can cause metabolite signals to drift and fluctuate between batches, manifesting as changes in characteristics such as retention time, peak intensity, and mass-to-charge ratio of metabolite peaks. The existence of batch effects makes data from different batches of samples incomparable, seriously interfering with subsequent data analysis and biological interpretation.

[0004] To eliminate the batch effect, various batch correction methods for metabolomics data have been developed. Early methods were mainly based on internal standards, which involve adding a known concentration of an internal standard to the sample and using the signal change of the internal standard to correct for the peak intensity of metabolites. However, this method has limitations. It is difficult to find a suitable internal standard to accurately normalize all metabolites because the chemical properties and chromatographic behaviors of different metabolites vary greatly, and the internal standard may not fully reflect their true changes. In addition, internal standards are susceptible to chromatographic peak overlap and ion suppression effects, further reducing the accuracy of correction.

[0005] Sample scalar correction methods rely on the self-averaging properties of the data, eliminating batch-to-batch differences by scaling the peak intensities of metabolites across all samples. This method is simple and easy to implement, but it assumes that all metabolites exhibit consistent trends across different batches, which is often not the case in reality. Therefore, sample scalar correction methods are ineffective in handling complex batch effects and may introduce new errors.

[0006] Regression methods based on quality control (QC) samples are currently a widely used batch calibration method. These methods insert a certain number of QC samples into the sample sequence and use the signal changes of the QC samples to build a classification model, which is then used to correct the data of other samples. Linear regression methods assume that the drift of metabolite signals is linear, but in actual metabolomics data, signal drift often exhibits complex nonlinear characteristics, thus limiting the calibration effect of linear regression methods. While locally weighted regression with scatter smoothing (LOESS) can fit nonlinear signal drift to some extent, its predictive ability is relatively weak and it is prone to overfitting, especially when the data volume is small or the signal drift is complex.

[0007] Support Vector Regression (SVR) has also been applied to batch correction of metabolomics data. SVR can handle nonlinear data by introducing a kernel function to map the data to a high-dimensional space for linear regression. However, traditional SVR methods typically use a single kernel function, which struggles to fully capture the complex characteristics of metabolomics data. Different kernel functions are suitable for different types of data distributions, and a single kernel function may not fully realize its advantages, leading to inaccurate fitting of nonlinear signal drift and thus affecting the effectiveness of batch correction. Methods such as multinomial regression and neural networks can eliminate systematic errors between batches, resulting in more stable and comparable metabolite peak intensity data after correction.

[0008] In summary, existing batch correction methods for metabolomics data have certain limitations in handling complex nonlinear signal drift and batch effects, and cannot meet the requirements of large-scale metabolomics research for data accuracy and reliability. Summary of the Invention

[0009] This invention discloses a batch correction method for metabolomics data based on multi-kernel learning, aiming to address the limitations of existing methods in handling nonlinear signal drift and complex batch effects. This method utilizes multi-kernel learning technology to adaptively fuse the advantages of different kernel functions, constructing a model that fits the actual data, accurately capturing the complex drift characteristics of metabolite signals, and achieving efficient and accurate normalization of metabolomics data. Compared with traditional data normalization methods such as SVR and LOESS, this invention's method demonstrates superior performance in reducing metabolite peak intensity variability and significantly improves data stability. In subsequent multivariate statistical analyses, classification accuracy is greatly improved, and the comparability between different batches of data is also significantly enhanced, providing reliable data support for the discovery of disease biomarkers and playing a crucial role in large-scale metabolomics research.

[0010] The method of the present invention includes the following steps:

[0011] S1: Data Acquisition and Preprocessing

[0012] S11: Based on the purpose of metabolomics research, collect biological samples covering different physiological states, disease stages, or experimental conditions. The biological samples include human blood, urine, saliva, tissue samples, or cell and tissue extracts from animals and plants. Blood samples are centrifuged at 4°C and 3,000 rpm for 10 minutes to separate serum or plasma. Tissue samples are homogenized using a liquid nitrogen grinding combined with ultrasonic disruption. Metabolites are extracted using a methanol-water mixture. The processed samples are arranged in a predetermined order, and one quality control (QC) sample is inserted every 5-10 samples to ensure that the QC samples are evenly distributed in the sample sequence.

[0013] S12: Use a high-resolution mass spectrometer to analyze the sample and obtain raw mass spectrometry data. If using liquid chromatography-mass spectrometry (LC-MS), set a C18 column, acetonitrile-water mobile phase system with the addition of an appropriate amount of formic acid, flow rate of 0.3 mL / min, and mass spectrometry scan range of m / z 50-1500. If using gas chromatography-mass spectrometry (GC-MS), use an HP-5MS column, set an initial temperature of 40℃ and hold for 2 minutes, increase to 300℃ at 10℃ / min and hold for 5 minutes, and use full scan mass spectrometry acquisition mode.

[0014] S13: Use Meexplorer Ultimate software to convert the raw mass spectrometry data files, i.e., .d, .raw, etc., to mzXML or mzML format. Use Meexplorer Ultimate software and set appropriate peak width range, signal-to-noise ratio threshold, etc. to perform peak identification (clarify peak identification parameters), peak filtering (clarify peak filtering conditions), alignment, and removal of noise peaks, isotope peaks, adduct ion peaks, and metabolite peaks that were not detected in more than 80% of the samples.

[0015] S2: Construction of Multi-Kernel Learning Support Vector Machine (MKL-SVM) Model:

[0016] S21: Choose a linear kernel function Polynomial kernel function Where, x i and x j Let K represent the sample feature vectors in the metabolomics data, γ be the kernel coefficient, r be the offset, d be the polynomial order, T be the basic matrix / vector operation in linear algebra, and K be the Gaussian radial basis function (RBF). rbf (x i ,x j )=exp(-γ||x i -x j || 2), where exp is the natural exponential function, the kernel coefficient γ of the Gaussian radial basis function RBF is preferably in the range of 0.01-1.0, and the order d of the polynomial kernel function is preferably 2-4. The goodness of fit of the model is determined by principal component analysis (PCA) or linear regression (R²). 2 Determine whether the data has a linear trend. If R... 2 If the cumulative contribution rate of the first two principal components of PCA is ≥0.7 or ≥70%, the data is determined to have a certain linear trend, and both linear kernel function and RBF kernel function are selected; otherwise, the data is determined to be complex nonlinear data, and the use of polynomial kernel function is added.

[0017] S22: Construct a multi-kernel learning framework, using gradient descent / particle swarm optimization to adaptively learn the weights of different kernel functions. Using QC sample data as a reference, and incorporating information such as injection order and metabolite peak intensity as input features, construct an MKL-SVM model; define the multi-kernel function. Where K k (x i ,x j ) is the k-th kernel function, α k Let it be its corresponding weight, and satisfy the following conditions: By optimizing the objective function To determine the optimal weight combination, where y i f(x) is the true value. i ) represents the predicted value, and λ represents the regularization coefficient;

[0018] S23: The parameters of the MKL-SVM model are comprehensively optimized using cross-validation. These parameters include the penalty parameter C of the support vector machine, the order d of the polynomial kernel function, the kernel coefficient γ of the RBF kernel function, and the weight parameter α in multi-kernel learning. k The dataset is divided into training and validation sets in a 7:3 or 8:2 ratio. The model is trained on the training set and evaluated on the validation set. The optimal parameter combination for model performance is found by using grid search or random search, with the goal of minimizing the mean squared error or maximizing the classification accuracy. The parameter range for grid search is as follows: penalty parameter C∈[0.01,0.1,1,10,100], RBF kernel γ∈[0.0001,0.001,0.01,0.1,1].

[0019] S3: Data Batch Correction

[0020] S31: The MKL-SVM model is trained in depth using preprocessed QC sample data. The injection order and metabolite peak intensity of the QC samples are used as input features, and the true values ​​of metabolite peak intensities are used as output labels. During training, stochastic gradient descent or particle swarm optimization is used to optimize the model. The number of training rounds is set to 100-500 and the learning rate is set to 0.01-0.1. The model parameters are continuously adjusted so that the model learns the drift pattern of metabolite signals over time.

[0021] S32: Input the sample data to be corrected into the trained MKL-SVM model. Based on the drift pattern learned from the QC samples, the model predicts the theoretical values ​​of metabolite peak intensities in the sample to be corrected. Compare the actual measured values ​​with the theoretical values ​​and use a linear transformation y. corrected =α×y measured +b, where y measured y represents the actual measured peak intensity of the metabolite. corrected The corrected peak intensity is represented by a, where a is the slope coefficient and b is the intercept coefficient, both obtained through model fitting.

[0022] S4: Evaluation of Correction Effectiveness

[0023] S41: Calculate the relative standard deviation (RSD) of the metabolite peak intensities before and after correction. The formula is as follows: Where SD is the standard deviation. The average value is used. The RSD of metabolite peak intensities in QC samples before and after correction and in all samples are calculated and an RSD distribution histogram is plotted. By calculating and comparing the RSD values ​​before and after correction, and combining the trend analysis of the RSD values ​​with their changing trends and distribution characteristics, the effect of the correction method on improving data stability can be evaluated.

[0024] S42: Perform principal component analysis (PCA) and partial least squares discriminant analysis (PLS-DA) on the corrected data, using multivariate statistical analysis. In PCA, calculate the covariance matrix of the data, and use the singular value decomposition (SVD) algorithm to solve for eigenvalues ​​and eigenvectors, projecting the high-dimensional data into a two-dimensional / three-dimensional space to observe the clustering and distribution patterns among samples. In PLS-DA, use the SIMPLS or NIPALS algorithm to build a predictive model to identify differential metabolites among samples from different groups. By comparing the clustering density of the PCA score maps before and after correction, the classification accuracy of the PLS-DA model, and the variable importance projection (VIP) value, evaluate the impact of the correction method on sample classification accuracy and differential metabolite screening.

[0025] S5: Discovery and application of biomarkers:

[0026] S51: Based on the corrected data, statistical methods such as t-test and ANOVA are used to screen differential metabolites among samples from different groups with a p-value less than 0.05 and a fold change greater than 2 or less than 0.5 as thresholds. Metabolic pathway analysis is then performed using databases and tools such as MeexplorerUltimate and KEGG to identify potential biomarkers. S52: A classification model or support vector machine classification model based on biomarkers is established, and the screened biomarkers are applied to disease prediction and prognostic assessment.

[0027] During model parameter optimization, either stochastic gradient descent or batch gradient descent is used for parameter updates. Stochastic gradient descent updates parameters using only one sample at a time, making it suitable for rapid iteration on large-scale datasets; batch gradient descent updates parameters using all samples, resulting in slower convergence but more stable parameter updates, and can be flexibly chosen based on the data size.

[0028] Particle swarm optimization algorithm implementation: The formulas for updating particle velocity and position are as follows: and in Let be the velocity of particle i in generation t. Let be the position of particle i in generation t, ω be the inertial weight (ranging from 0.4 to 0.9), C1 and C2 be the acceleration constants, and r1 and r2 be random numbers between 0 and 1. i Let be the individual optimal position of the particle, and g be the global optimal position. This algorithm can efficiently search for the optimal weight combination of a multi-kernel function.

[0029] Parameter search strategies: The grid search method sets a discrete value sequence for each parameter to be optimized and traverses all combinations to find the optimal solution, which is suitable for scenarios with fewer parameter dimensions; the random search method randomly samples and evaluates within the parameter value range, which has lower computational cost and is more efficient in high-dimensional parameter spaces, and the appropriate strategy can be selected according to the parameter size.

[0030] Metabolic pathway analysis workflow: Based on the screened differentially enriched metabolites, pathway enrichment analysis is performed using Meexplorer Ultimate and the KEGG database. The number and proportion of enriched metabolites in each pathway are statistically analyzed, enrichment scores are calculated, and significantly enriched pathways (p<0.05) are screened to identify the core biological processes involved by the differentially enriched metabolites. In Meexplorer Ultimate, the "Pathway Enrichment" module can be selected, and the species database and p-value threshold parameters can be set to complete the analysis.

[0031] Classification model construction method: When building a biomarker classification model, the maximum likelihood estimation method is used to estimate the model parameters, and the optimal parameters are solved by maximizing the sample likelihood function; when building a support vector machine classification model, a linear kernel or RBF kernel function is selected, and the penalty parameter C is optimized by combining grid search. The optimal model parameters are determined by 5-fold cross-validation. The final model can be used for disease risk prediction and prognosis assessment.

[0032] Technical effect comparison

[0033] Compared with traditional SVR methods, the method of this invention improves the RSD reduction rate of more than 90% of metabolites in QC samples by 20%-35%; compared with the LOESS method, the classification accuracy of the PLS-DA model is improved by an average of 12.6% (specific data: accuracy before correction 76.3%, after correction by this method 88.9%, after correction by LOESS 81.5%). Through multi-kernel function adaptive fusion and parameter optimization strategies, this method can more accurately capture the nonlinear drift characteristics of metabolite signals, significantly improving the stability and comparability of cross-batch data.

[0034] The specific implementation process of the technical solution of this invention is as follows: First, a multi-kernel learning support vector machine (MKL-SVM) model is constructed using the injection order and metabolite peak intensity data of quality control (QC) samples. The selection criteria for the kernel function are as follows: for linear trend data, both linear kernel functions and Gaussian radial basis function (RBF) kernel functions are used simultaneously; for complex nonlinear data, a polynomial kernel function is added. Second, the kernel function weights are optimized using gradient descent or particle swarm optimization (PSO) algorithm. The inertia weight ω of the PSO algorithm ranges from 0.4 to 0.9, and the acceleration constants c1 and c2 are both set to 1.5. Then, the sample data to be corrected is input into the trained model, and batch correction is performed using the linear transformation formula y_corrected=a×y_measured+b, where the coefficients a and b are determined through model optimization. Finally, the relative standard deviation (RSD) calculation formula is used to verify the correction effect, ensuring that the RSD value of the QC samples meets the threshold standard required by the technical requirements.

[0035] Compared with the traditional SVR method, the method of this invention improves the RSD reduction rate of more than 90% of metabolites in QC samples by 20%-35%; compared with the LOESS method, the classification accuracy of the PLS-DA model is improved by an average of 12.6% (specific data: accuracy before correction 76.3%, after correction by this method 88.9%, after correction by LOESS 81.5%).

[0036] The specific implementation process of the technical solution of this invention is as follows:

[0037] 1. Sample preprocessing operation specifications:

[0038] After blood sample collection, EDTA anticoagulant (final concentration 1.8 mg / mL) should be added immediately, centrifuged at 3000 rpm for 10 minutes at 4°C, and the separated serum / plasma should be aliquoted into 0.5 mL centrifuge tubes and stored at -80°C for no more than 6 months.

[0039] After the tissue samples were collected, they were quickly frozen in liquid nitrogen and then stored at -80°C. During the experiment, 50 mg of tissue block was taken and 1 mL of methanol-water mixture (volume ratio 7:3) was added. The mixture was homogenized in an ice bath for 3 minutes using a tissue homogenizer (power 300W). The tissue was then sonicated (power 200W, working for 3 seconds and stopping for 2 seconds, for a total of 10 minutes). The mixture was then centrifuged at 12000 rpm at 4°C for 15 minutes and the supernatant was collected.

[0040] QC sample preparation: Mix 10 μL of each of the samples, insert one QC sample for every 5 samples, and insert one QC sample for every 10 samples when the total number of samples exceeds 100, to ensure that the QC samples are evenly distributed in the sequence.

[0041] In this step, the blood sample centrifugation parameters (4℃, 3,000 rpm × 10 minutes) and the QC sample insertion standard (1 insertion for every 5 samples) correspond to S11 of claim 1.

[0042] 2. Mass spectrometry parameter settings details:

[0043] LC-MS: An Agilent ZORBAX SB-C18 column (250 mm × 4.6 mm, 5 μm) was used. Mobile phase A was 0.1% formic acid aqueous solution, and mobile phase B was 0.1% formic acid acetonitrile solution. Gradient elution program: 0-5 min B phase 5%-30%, 5-15 min B phase 30%-80%, 15-20 min B phase 80%-95%, flow rate 0.3 mL / min, column temperature 30℃, injection volume 5 μL.

[0044] GC-MS: HP-5MS column (30m×0.25mm, 0.25μm), carrier gas helium (purity 99.999%), flow rate 1mL / min, injection port temperature 280℃, split ratio 10:1, mass spectrometer ion source temperature 230℃, quadrupole temperature 150℃.

[0045] 3. Key parameters for data processing:

[0046] Use Mexplorer Ultimate 3.0 software to convert the data format. Peak identification parameters: peak width 0.5-2 min, signal-to-noise ratio threshold 5, minimum peak height 10000.

[0047] Peak alignment employs a dynamic time warping algorithm, retaining a time deviation tolerance of 0.2 min and an m / z deviation tolerance of 5 ppm.

[0048] 4. Specific steps for model construction:

[0049] Kernel function selection: When the cumulative contribution rate of the first two principal components of PCA is ≥70%, a linear kernel (γ=1.0) and an RBF kernel (γ=0.1) are selected; otherwise, a polynomial kernel (γ=0.1, r=0.5, d=3) is added; in this embodiment, a linear kernel and an RBF kernel (γ=0.1) are selected, which meets the kernel function selection criteria of claim 1S21 when the cumulative contribution rate of the first two principal components of PCA is ≥70%.

[0050] Particle swarm optimization algorithm parameters: population size 30, maximum number of iterations 100, inertia weight ω decreases linearly from 0.9 to 0.4, c1=c2=1.5;

[0051] The cross-validation uses 5-fold cross-validation, and the grid search parameter range is: C∈[0.1,1,10,100], γ∈[0.001,0.01,0.1,1].

[0052] In a multi-kernel learning framework, the initial weights of different kernel functions are set to be evenly distributed (i.e., α). k =1 / K, where K is the total number of kernel functions), during the iteration of the particle swarm optimization algorithm, the weight sum must be forcibly verified to be 1 after each iteration (i.e., ∑α). k =1), to avoid weight drift that could cause the model to fail.

[0053] 5. Correction and evaluation implementation methods:

[0054] The linear transformation coefficients a and b are calculated using the least squares fitting method:

[0055] a=cov(y_measured,y_predicted) / var(y_measured),

[0056] b=mean(y_predicted)-a×mean(y_measured);

[0057] RSD calculations should include all QC samples (no less than 20) and 20% of the samples randomly selected from each batch. When plotting the histogram, the horizontal axis should be spaced 0.5% apart, and the vertical axis should be the frequency.

[0058] PCA analysis was implemented using the scikit-learn library, with 50 iterations for the SVD algorithm. PLS-DA used the SIMPLS algorithm with -components=5 and 100 iterations. Attached Figure Description

[0059] Figure 1 This application provides a flowchart of the multi-kernel correction process in a batch correction method for metabolomics data based on multi-kernel learning.

[0060] Figure 2 The RT-m / z plot of peak table data after peak extraction is provided in the batch correction method for metabolomics data based on multi-kernel learning provided in this application.

[0061] Figure 3 The QC quality control samples and heatmaps of actual samples after batch correction are provided in this application for a batch correction method for metabolomics data based on multi-kernel learning.

[0062] Figure 4 The results of PCA multivariate correction analysis of peak table data after batch correction in a batch correction method for metabolomics data based on multi-kernel learning provided in this application.

[0063] Figure 5 The PLS-DA score results of the peak table data after batch correction in the batch correction method for metabolomics data based on multi-kernel learning provided in this application. Detailed Implementation

[0064] The specific technical solution is as follows:

[0065] 1. Sample processing

[0066] S11: Based on the purpose of metabolomics research, collect biological samples covering different physiological states, disease stages, or experimental conditions. The biological samples include human blood, urine, saliva, tissue samples, or cell and tissue extracts from animals and plants. Blood samples are centrifuged at 4°C and 3,000 rpm for 10 minutes to separate serum or plasma. Tissue samples are homogenized using a liquid nitrogen grinding combined with ultrasonic disruption. Metabolites are extracted using a methanol-water mixture. The processed samples are arranged in a predetermined order, and one quality control (QC) sample is inserted every 5-10 samples to ensure that the QC samples are evenly distributed in the sample sequence.

[0067] S12: Use a high-resolution mass spectrometer to analyze the sample and obtain raw mass spectrometry data. If using liquid chromatography-mass spectrometry (LC-MS), set a C18 column, acetonitrile-water mobile phase system with the addition of an appropriate amount of formic acid, flow rate of 0.3 mL / min, and mass spectrometry scan range of m / z 50-1500. If using gas chromatography-mass spectrometry (GC-MS), use an HP-5MS column, set an initial temperature of 40℃ and hold for 2 minutes, increase to 300℃ at 10℃ / min and hold for 5 minutes, and use full scan mass spectrometry acquisition mode.

[0068] Samples: 10 μL of serum from 60 healthy individuals (control group) and 60 diabetic patients (case group) were taken from each group and mixed to prepare QC samples. One QC sample was inserted for every 5 samples (a total of 24 QC samples) in accordance with the requirement of S11: "One quality control sample is inserted for every 5-10 samples to ensure that the QC sample distribution interval is uniform".

[0069] Pretreatment: Serum was centrifuged at 3000 rpm for 10 minutes at 4℃ and stored at -80℃. Operating parameters conformed to S11, which specifies "blood samples are centrifuged at 300 rpm for 10 minutes to separate serum." LC-MS analysis was performed on an Agilent ZORBAX SB-C18 column. Mobile phase A: 0.1% formic acid in water, mobile phase B: 0.1% formic acid in acetonitrile, 0-5 min, B = 5%-30%, conforming to S12.

[0070] The technical features include "LC-MS setup with a C18 column and formic acid added to the acetonitrile-water mobile phase system".

[0071] 2. Data Preprocessing

[0072] S13: Use Meexplorer Ultimate software to convert the raw mass spectrometry data files, i.e., .d, .raw, etc., to mzXML or mzML format. Use Meexplorer Ultimate software and set appropriate peak width range, signal-to-noise ratio threshold, etc. to perform peak identification (clarify peak identification parameters), peak filtering (clarify peak filtering conditions), alignment, and removal of noise peaks, isotope peaks, adduct ion peaks, and metabolite peaks that were not detected in more than 80% of the samples.

[0073] Format conversion: Mexplorer Ultimate 3.0 converts .raw to mzXML. Peak recognition parameters: peak width 0.8min, signal-to-noise ratio = 5, minimum peak height = 10000, which corresponds to the requirement of S13 "Use Mexplorer Ultimate software to convert the format and set peak recognition parameters such as peak width and signal-to-noise ratio".

[0074] Peak alignment: Dynamic time warping algorithm (retaining time deviation ≤ 0.2 min) removes undetected peaks in 80% of samples, resulting in a matrix of 1500 metabolite features, thus fulfilling the processing requirement of S13 "peak alignment and removal of more than 80% of undetected metabolite peaks".

[0075] 3. Model Building

[0076] S21: Choose a linear kernel function Polynomial kernel function Where, x i and x jLet K represent the sample feature vectors in the metabolomics data, γ be the kernel coefficient, r be the offset, d be the polynomial order, T be the basic matrix / vector operation in linear algebra, and K be the Gaussian radial basis function (RBF). rbf (x i ,x j )=exp(-γ||x i -x j || 2 ), where exp is the natural exponential function, the kernel coefficient γ of the Gaussian radial basis function RBF is preferably in the range of 0.01-1.0, and the order d of the polynomial kernel function is preferably 2-4. The goodness of fit of the model is determined by principal component analysis (PCA) or linear regression (R²). 2 Determine whether the data has a linear trend. If R... 2 If the cumulative contribution rate of the first two principal components of PCA is ≥0.7 or ≥70%, the data is determined to have a certain linear trend, and both linear kernel function and RBF kernel function are selected; otherwise, the data is determined to be complex nonlinear data, and the use of polynomial kernel function is added.

[0077] S22: Construct a multi-kernel learning framework, using gradient descent / particle swarm optimization to adaptively learn the weights of different kernel functions. Using QC sample data as a reference, and incorporating information such as injection order and metabolite peak intensity as input features, construct an MKL-SVM model; define the multi-kernel function. Where K k (x i ,x j ) is the k-th kernel function, α k Let it be its corresponding weight, and satisfy the following conditions: By optimizing the objective function To determine the optimal weight combination, where y i f(x) is the true value. i ) represents the predicted value, and λ represents the regularization coefficient;

[0078] S23: The parameters of the MKL-SVM model are comprehensively optimized using cross-validation. These parameters include the penalty parameter C of the support vector machine, the order d of the multinomial kernel function, the kernel coefficient γ of the RBF kernel function, and the weight parameter a in multi-kernel learning. k The dataset is divided into training and validation sets in a 7:3 or 8:2 ratio. The model is trained on the training set and evaluated on the validation set. The optimal parameter combination for model performance is found by using grid search or random search, with the goal of minimizing the mean squared error or maximizing the classification accuracy. The parameter range for grid search is as follows: penalty parameter C∈[0.01,0.1,1,10,100], RBF kernel γ∈[0.0001,0.001,0.01,0.1,1].

[0079] S31: The MKL-SVM model is trained in depth using preprocessed QC sample data. The injection order and metabolite peak intensity of the QC samples are used as input features, and the true values ​​of metabolite peak intensities are used as output labels. During training, stochastic gradient descent or particle swarm optimization is used to optimize the model. The number of training rounds is set to 100-500 and the learning rate is set to 0.01-0.1. The model parameters are continuously adjusted so that the model learns the drift pattern of metabolite signals over time.

[0080] Kernel function selection: The cumulative contribution rate of the first two principal components of PCA is 73% (≥70%). Linear kernel (γ=1.0) and RBF kernel (γ=0.1) are selected, which fully meets the judgment criterion S21 "when the cumulative contribution rate of the first two principal components of PCA is ≥70%, both linear function and RBF kernel function are selected".

[0081] Parameter optimization: 5-fold cross-validation and grid search yielded the optimal parameters C=10 and γ.

[0082] =0.1, mean square error =0.028, parameter search range covers S23 "penalty parameter C∈

[0083] The technical constraints of "[0.01,100], RBF kernel γ∈[0.0001,1]" are met, and the method characteristics of "network search traversing parameter combinations to find the optimal solution" are also met.

[0084] Training: Stochastic gradient descent (learning rate 0.05, 300 rounds), input QC sample injection sequence (1-24) and peak intensity, output true values, fit drift patterns (such as nonlinear decay of peak intensity with the number of injections) in accordance with the parameter requirements of S31 "use stochastic gradient descent optimization method, set 100-500 training rounds and 0.01-0.1 learning rate", and meet the technical feature of "stochastic gradient descent updates parameters with one sample each time".

[0085] 4. Batch calibration

[0086] S32: Input the sample data to be corrected into the trained MKL-SVM model. Based on the drift pattern learned from the QC samples, the model predicts the theoretical values ​​of metabolite peak intensities in the sample to be corrected. Compare the actual measured values ​​with the theoretical values ​​and use a linear transformation y. corrected =a×y measured , where y measured y represents the actual measured peak intensity of the metabolite. corrected The corrected peak intensity is represented by a, where a is the slope coefficient and b is the intercept coefficient, both obtained through model fitting.

[0087] Correction calculation: The measured value of a certain metabolite is y_measured = 1500, and the model prediction is y_predicted = 1420. a = cov(1500, 1420) / var(1500) = 0.89, b = 1420 - 0.89 × 1500 = 45. After correction, y_corrested = 0.89 × 1500 + 48 = 1380. The correction formula of S32 "linear transformation y_corrested = a × y_measured + b" is strictly adopted, where a and b are determined by model fitting.

[0088] 3. Effectiveness Evaluation

[0089] S41: Calculate the relative standard deviation (RSD) of the metabolite peak intensities before and after correction. The formula is as follows: Where SD is the standard deviation and x is the mean; RSD is calculated for the peak intensity of metabolites in QC samples before and after correction and in all samples, and RSD distribution histograms are plotted. By calculating and comparing the RSD values ​​before and after correction, and combining the trend analysis of the RSD values ​​with their changing trends and distribution characteristics, the effect of the correction method on improving data stability can be evaluated.

[0090] S42: Perform principal component analysis (PCA) and partial least squares discriminant analysis (PLS-DA) on the corrected data, using multivariate statistical analysis. In PCA, calculate the covariance matrix of the data, and use the singular value decomposition (SVD) algorithm to solve for eigenvalues ​​and eigenvectors, projecting the high-dimensional data into a two-dimensional / three-dimensional space to observe the clustering and distribution patterns among samples. In PLS-DA, use the SIMPLS or NIPALS algorithm to build a predictive model to identify differential metabolites among samples from different groups. By comparing the clustering density of the PCA score maps before and after correction, the classification accuracy of the PLS-DA model, and the variable importance projection (VIP) value, evaluate the impact of the correction method on sample classification accuracy and differential metabolite screening.

[0091] RSD: The average RSD of the QC samples before correction was 19.2%, and after correction it was 5.7%. The evaluation was completed using the calculation formula "RSD = SD / average value × 100%" in S41, which realizes the quantitative verification of data stability.

[0092] The PLS-DA classification accuracy was improved from 75.6% to 89.3%, using the method of S42 "evaluating the correction effect through the classification accuracy of the PLS-DA model", and the results confirmed the advantages of the present invention compared with traditional methods.

[0093] 6. Biomarker discovery

[0094] S51: Based on the corrected data, statistical methods such as t-test and ANOVA are used to screen differential metabolites among samples from different groups with a p-value less than 0.05 and a fold change greater than 2 or less than 0.5 as thresholds. Metabolic pathway analysis is then performed using databases and tools such as MeexplorerUltimate and KEGG to identify potential biomarkers. S52: A classification model or support vector machine classification model based on biomarkers is established, and the screened biomarkers are applied to disease prediction and prognostic assessment.

[0095] The gradient descent optimization method should be explained according to its technical characteristics: During the model parameter optimization process, either stochastic gradient descent or batch gradient descent is used for parameter updates. Stochastic gradient descent updates parameters using only one sample at a time, making it suitable for rapid iteration on large-scale datasets; batch gradient descent updates parameters using all samples, although its convergence speed is slower, its parameter updates are more stable, and the choice can be made flexibly according to the data size.

[0096] Particle swarm optimization algorithm implementation: The formulas for updating particle velocity and position are as follows: and in Let be the velocity of particle i in generation t. Let be the position of particle i in generation t, ω be the inertial weight (ranging from 0.4 to 0.9), C1 and C2 be the acceleration constants, and r1 and r2 be random numbers between 0 and 1. i Let be the individual optimal position of the particle, and g be the global optimal position. This algorithm can efficiently search for the optimal weight combination of a multi-kernel function.

[0097] Parameter search strategies: The grid search method sets a discrete value sequence for each parameter to be optimized and traverses all combinations to find the optimal solution, which is suitable for scenarios with fewer parameter dimensions; the random search method randomly samples and evaluates within the parameter value range, which has lower computational cost and is more efficient in high-dimensional parameter spaces, and the appropriate strategy can be selected according to the parameter size.

[0098] Metabolic pathway analysis workflow: Based on the screened differentially enriched metabolites, pathway enrichment analysis is performed using Meexplorer Ultimate and the KEGG database. The number and proportion of enriched metabolites in each pathway are statistically analyzed, enrichment scores are calculated, and significantly enriched pathways (p<0.05) are screened to identify the core biological processes involved by the differentially enriched metabolites. In Meexplorer Ultimate, the "Pathway Enrichment" module can be selected, and the species database and p-value threshold parameters can be set to complete the analysis.

[0099] Classification model construction method: When building a biomarker classification model, the maximum likelihood estimation method is used to estimate the model parameters, and the optimal parameters are solved by maximizing the sample likelihood function; when building a support vector machine classification model, a linear kernel or RBF kernel function is selected, and the penalty parameter C is optimized by combining grid search. The optimal model parameters are determined by 5-fold cross-validation. The final model can be used for disease risk prediction and prognosis assessment.

[0100] Thirty differentially expressed metabolites (p<0.05, fold change>2) were screened and enriched in the "insulin signaling pathway" (KEGGID:map04910), which meets the technical requirements of S51: "differentially expressed metabolites are screened with p-values ​​less than 0.05 and fold changes greater than 2 as thresholds, and metabolic pathway analysis is performed in conjunction with the KEGG database".

[0101] An SVM model (RBF kernel, C=10) was constructed, achieving a disease prediction accuracy of 88.5%, following the construction specification of "selecting the RBF kernel function and optimizing the penalty parameter C when building a support vector machine classification model".

[0102] This invention provides a batch correction method for metabolomics data based on multi-kernel learning, which overcomes the limitations of existing batch correction methods for metabolomics data when dealing with complex nonlinear signal drift and batch effects, thereby improving the accuracy and reliability of the data. After introducing the basic principles of this invention, various non-limiting embodiments will be described in detail below with reference to the accompanying drawings.

[0103] like Figure 1 As shown, the specific implementation process is as follows:

[0104] 1. Data Input:

[0105] Input the peak table data (including retention time, m / z, peak intensity, e.g., a matrix of 1000 samples × 1200 metabolites) for batches 1 to m.

[0106] 2. Data preprocessing:

[0107] Missing value handling: k-nearest neighbor imputation (k=5, based on the mean peak intensity of similar samples) is used.

[0108] Standardization: z-score standardization (x' = (x - μ) / σ, where μ is the mean and σ is the standard deviation).

[0109] Peak alignment: Dynamic time warping algorithm (allowing retention time deviation ±0.2 min).

[0110] 3. Kernel function selection:

[0111] Calculate the cumulative contribution rate of the first two principal components of PCA (e.g., 72%), and determine it as a linear trend.

[0112] Linear kernels (γ = 1.0) and RBF kernels (γ = 0.1) were selected.

[0113] 4. Construction of a multi-core learning framework:

[0114] Define a multi-kernel function K = a1K linear +a2K rβp The initial a1 = a2 = 0.5 is used to optimize the objective function.

[0115] (λ=0.01), the optimal values ​​of a1=0.6 and a2=0.4 are obtained by particle swarm optimization.

[0116] 5. Model parameter optimization:

[0117] A grid search was performed to traverse C = 0.01-100 and γ = 0.0001-1, and the optimal parameters (C = 10, γ = 0.1) were obtained through 5-fold verification.

[0118] 6. Data correction:

[0119] Model training: Input the injection sequence (1-20) and peak intensity of the QC samples to learn the drift law of the signal over time (such as the nonlinear decay trend of peak intensity with the increase of injection number).

[0120] Correction calculation: For each metabolite of each sample, adjust using y_corrected = a × y_measured + b (a and b are obtained by fitting the model to the QC sample).

[0121] 7. Effectiveness Evaluation and Application:

[0122] RSD calculation: Plot the RSD histogram of QC samples before and after correction (horizontal axis interval 0.5%, such as 80% of metabolites RSD < 8% after correction).

[0123] Multivariate statistics: PCA score plots showed that the corrected samples were more tightly clustered (within-group distance reduced by 40%), and the number of metabolites with PLS-DAVIP values ​​>1 increased by 30%.

[0124] like Figure 2As shown, in the batch correction method for metabolomics data based on multi-core learning and support vector machine (MKL-SVM) of this invention, the distribution of metabolites before correction is illustrated by a scatter plot of the distribution of metabolites before correction in terms of mass-to-charge ratio (m / z) and retention time (tR). Each scatter point in the plot corresponds to a metabolite detection point, with the horizontal axis representing retention time and the vertical axis representing mass-to-charge ratio. By observing the distribution of these scatter points, the distribution of metabolites before correction in these two key dimensions can be intuitively presented, reflecting the distribution characteristics of metabolite signals in the original data and providing an original state reference for subsequent evaluation of the batch correction effect.

[0125] like Figure 3 As shown, this invention, based on multi-core learning and support vector machine (MKL-SVM) batch correction method for metabolomics data, includes a heatmap used to display data correlation or clustering relationships. Different colors represent the degree of correlation or clustering between data, with the intensity of the color reflecting the strength of the correlation or the tightness of the clustering. This heatmap allows observation of similarities or differences between different samples or metabolites, aiding in the analysis of sample grouping characteristics and metabolite association patterns in metabolomics research. It can also be used to compare structural changes in data before and after correction.

[0126] like Figure 4 As shown, in the batch correction method for metabolomics data based on multi-core learning and support vector machine (MKL-SVM) of this invention, the principal component analysis (PCA) results plot shows the contribution rates of the first five principal components (PC1-PC5), which are 48.5%, 13.3%, 6%, 3.6%, and 2.8%, respectively. Each principal component's corresponding subplot includes both a scatter plot showing the distribution of principal component scores, reflecting the distribution of samples in that principal component space, and a histogram showing the frequency distribution of principal components. Through this information, we can understand the dispersion and distribution characteristics of samples across different principal component dimensions, assess the overall structure of the data and batch-to-batch differences, and assist in judging the batch correction effect.

[0127] like Figure 5 As shown, in the batch correction method for metabolomics data based on Multi-Kernel Learning Support Vector Machine (MKL-SVM) of this invention, the scatter plot used to display the sample classification or clustering situation shows that different colored scatter points represent samples of different categories, and the distribution of scatter points in two-dimensional space shows the clustering relationship or classification boundary between samples. By observing this plot, one can intuitively judge the changes in clustering or classification of samples before and after batch correction, and evaluate the impact of the correction method on the accuracy and stability of sample grouping.

Claims

1. A batch correction method for metabolomics data based on multi-kernel learning, characterized in that, The method includes: S1: Data Acquisition and Preprocessing S11: Based on the purpose of metabolomics research, collect biological samples covering different physiological states, disease stages, or experimental conditions. The biological samples include human blood, urine, saliva, tissue samples, or cell and tissue extracts from animals and plants. Blood samples are centrifuged at 4°C and 3,000 rpm for 10 minutes to separate serum or plasma. Tissue samples are homogenized using a liquid nitrogen grinding combined with ultrasonic disruption. Metabolites are extracted using a methanol-water mixture. The processed samples are arranged in a predetermined order, and one quality control (QC) sample is inserted every 5-10 samples to ensure that the QC samples are evenly distributed in the sample sequence. S12: Use a high-resolution mass spectrometer to analyze the sample and obtain raw mass spectrometry data. If using liquid chromatography-mass spectrometry (LC-MS), you can set a C18 column, acetonitrile-water mobile phase system with the addition of an appropriate amount of formic acid, flow rate of 0.3 mL / min, and mass spectrometry scan range of m / z 50-1500. If using gas chromatography-mass spectrometry (GC-MS), you can use an HP-5MS column, set an initial temperature of 40℃ and hold for 2 minutes, increase to 300℃ at 10℃ / min and hold for 5 minutes, and use full scan mass spectrometry acquisition mode. S13: Use Meexplorer Ultimate software to convert the raw mass spectrometry data files, i.e., .d, .raw, etc., to mzXML or mzML format. Use Meexplorer Ultimate software and set appropriate peak width range, signal-to-noise ratio threshold, etc. to perform peak identification (clarify peak identification parameters), peak filtering (clarify peak filtering conditions), alignment, and removal of noise peaks, isotope peaks, adduct ion peaks, and metabolite peaks that were not detected in more than 80% of the samples. S2: Construction of Multi-Kernel Learning Support Vector Machine (MKL-SVM) Model: S21: Choose a linear kernel function Polynomial kernel function Where, x i and x j Let K represent the sample feature vectors in the metabolomics data, γ be the kernel coefficient, r be the offset, d be the polynomial order, T be the basic matrix / vector operation in linear algebra, and K be the Gaussian radial basis function (RBF). rbf (x i ,x j )=exp(-γ||x i -x j || 2 ), where exp is the natural exponential function, the kernel coefficient γ of the Gaussian radial basis function RBF ranges from 0.01 to 1.0, and the order d of the polynomial kernel function is 2 to 4. The goodness of fit of the model is determined by principal component analysis (PCA) or linear regression (R²). 2 Determine whether the data has a linear trend. If R... 2 If the cumulative contribution rate of the first two principal components of PCA is ≥0.7 or ≥70%, the data is determined to have a certain linear trend, and both linear kernel function and RBF kernel function are selected; otherwise, the data is determined to be complex nonlinear data, and the use of polynomial kernel function is added. S22: Construct a multi-kernel learning framework, using gradient descent / particle swarm optimization to adaptively learn the weights of different kernel functions. Using QC sample data as a reference, and incorporating information such as injection order and metabolite peak intensity as input features, construct an MKL-SVM model; define the multi-kernel function. Where K k (x i ,x j ) is the k-th kernel function, α k Let it be its corresponding weight, and satisfy the following conditions: By optimizing the objective function To determine the optimal weight combination, where y i f(x) is the true value. i ) represents the predicted value, and λ represents the regularization coefficient; S23: The parameters of the MKL-SVM model are comprehensively optimized using cross-validation. These parameters include the penalty parameter C of the support vector machine, the order d of the polynomial kernel function, the kernel coefficient γ of the RBF kernel function, and the weight parameter α in multi-kernel learning. k The dataset is divided into training and validation sets in a 7:3 or 8:2 ratio. The model is trained on the training set and evaluated on the validation set. The optimal parameter combination for model performance is found by using grid search or random search, with the goal of minimizing the mean squared error or maximizing the classification accuracy. The parameter range for grid search is as follows: penalty parameter C∈[0.01,0.1,1,10,100], RBF kernel γ∈[0.0001,0.001,0.01,0.1,1]. S3: Data Batch Correction S31: The MKL-SVM model is trained in depth using preprocessed QC sample data. The injection order and metabolite peak intensity of the QC samples are used as input features, and the true values ​​of metabolite peak intensities are used as output labels. During training, stochastic gradient descent or particle swarm optimization is used to optimize the model. The number of training rounds is set to 100-500 and the learning rate is set to 0.01-0.

1. The model parameters are continuously adjusted so that the model learns the drift pattern of metabolite signals over time. S32: Input the sample data to be corrected into the trained MKL-SVM model. Based on the drift pattern learned from the QC samples, the model predicts the theoretical values ​​of metabolite peak intensities in the sample to be corrected. Compare the actual measured values ​​with the theoretical values ​​and use a linear transformation y. corrected =a×y measured +b, where y measured y represents the actual measured peak intensity of the metabolite. corrected The corrected peak intensity is represented by a, where a is the slope coefficient and b is the intercept coefficient, both obtained through model fitting. S4: Evaluation of Correction Effectiveness S41: Calculate the relative standard deviation (RSD) of the metabolite peak intensities before and after correction. The formula is as follows: Where SD is the standard deviation. The average value is used. The RSD of metabolite peak intensities in QC samples before and after correction and in all samples are calculated and an RSD distribution histogram is plotted. By calculating and comparing the RSD values ​​before and after correction, and combining the trend analysis of the RSD values ​​with their changing trends and distribution characteristics, the effect of the correction method on improving data stability can be evaluated. S42: Perform principal component analysis (PCA) and partial least squares discriminant analysis (PLS-DA) on the corrected data, using multivariate statistical analysis. In PCA, calculate the covariance matrix of the data, and use the singular value decomposition (SVD) algorithm to solve for eigenvalues ​​and eigenvectors, projecting the high-dimensional data into a two-dimensional / three-dimensional space to observe the clustering and distribution patterns among samples. In PLS-DA, use the SIMPLS or NIPALS algorithm to build a predictive model to identify differential metabolites among samples from different groups. By comparing the clustering density of the PCA score maps before and after correction, the classification accuracy of the PLS-DA model, and the variable importance projection (VIP) value, evaluate the impact of the correction method on sample classification accuracy and differential metabolite screening. S5: Discovery and application of biomarkers: S51: Based on the corrected data, statistical methods such as t-test and ANOVA are used to screen for differential metabolites among samples from different groups with a p-value less than 0.05 and a fold change greater than 2 or less than 0.5 as thresholds. Metabolic pathway analysis is then performed using databases and tools such as MeexplorerUltimate and KEGG to uncover potential biomarkers. S52: Establish a biomarker-based classification model or support vector machine classification model, and apply the selected biomarkers to disease prediction and prognosis assessment.

2. The batch correction method for metabolomics data based on multi-kernel learning according to claim 1, characterized in that, The gradient descent method includes: Parameter updates are performed using either stochastic gradient descent or batch gradient descent. Stochastic gradient descent updates parameters using only one sample at a time, while batch gradient descent updates parameters using all samples.

3. The batch correction method for metabolomics data based on multi-kernel learning according to claim 1, characterized in that, The particle swarm optimization algorithm includes: The formulas for updating the particle's velocity and position are as follows: and in Let be the velocity of particle i in generation t. Let ω be the position of particle i in generation t, ω be the inertial weight, c1 and c2 be the acceleration constants, t1 and r2 be random numbers between 0 and 1, and p be the position of particle i in generation t. i ω represents the individual optimal position of the particle, and g represents the global optimal position; the inertial weight ω ranges from 0.4 to 0.

9.

4. The batch correction method for metabolomics data based on multi-kernel learning according to claim 1, characterized in that, The method of grid search or random search aims to find the optimal parameter combination for model performance, with the goal of minimizing mean squared error or maximizing classification accuracy. This includes: The grid search method sets a series of discrete values ​​for each parameter to be optimized and finds the optimal solution by traversing all possible combinations of parameters; The random search method evaluates a certain number of parameter combinations randomly sampled within the range of parameter values.

5. The batch correction method for metabolomics data based on multi-kernel learning according to claim 1, characterized in that, The metabolic pathway analysis includes: Based on the enrichment of differential metabolites in the Meexplorer Ultimate and KEGG pathways, the main biological metabolic pathways in which the metabolites participate were identified. In the Meexplorer Ultimate tool, different analysis modules and parameters can be selected for metabolic pathway analysis.

6. The batch correction method for metabolomics data based on multi-kernel learning according to claim 1, characterized in that, Establishing biomarker-based classification models or support vector machine classification models includes: When establishing a classification model based on biomarkers, the maximum likelihood estimation method is used to estimate the parameters of the model; When establishing a support vector machine classification model, a kernel function and penalty parameter are selected for model training. The kernel function includes linear kernel and RBF kernel.

Images