Method and system for constructing genetic metabolic disease screening model based on metabolic network enhancement

Abstract

This invention discloses a method and system for constructing a genetic metabolic disease screening model based on metabolic network enhancement. The method performs liquid chromatography-mass spectrometry analysis on neonatal dried blood spot samples in the training set to obtain initial features of target metabolites, and stores high-resolution full-scan and data-independent acquisition spectra as a digital sample archive. On a metabolic reaction network with edge weights based on the product of carbon atom conservation, reaction competitiveness, and reaction reliability, forward and reverse band restart random walks are performed using differentially expressed metabolites as seeds. A comprehensive recommendation score is obtained after global baseline correction, and recommended nodes are selected. After correcting the predicted retention time with the measured retention time of the isotope internal standard, the chromatographic peaks of the recommended metabolites are extracted from the digital sample archive as complementary features, which are then used to train the classification model together with the initial features. This invention expands the feature dimensions of dried blood spot metabolic detection without requiring re-collection of blood or additional sample injections.

Description

Technical Field

[0001] This invention relates to the field of bioinformatics technology, specifically to a method and system for constructing a screening model for genetic metabolic diseases based on metabolic network enhancement. Background Technology

[0002] Newborn disease screening is a crucial part of the three-tiered prevention of birth defects. Tandem mass spectrometry can simultaneously determine dozens of metabolites, such as amino acids and acylcarnitines, in a single dried blood spot test, and has been used as a first-line screening method for various inherited metabolic diseases.

[0003] Current initial screening methods often rely on whether the concentration of one or a few biomarkers exceeds the cutoff value as the basis for judgment. However, an increase in the same biomarker can be caused by a variety of non-disease factors. For example, propionylcarnitine (C3) is elevated in methylmalonic acidemia and propionic acidemia, and it can also be elevated in maternal vitamin B12 deficiency, preterm birth, parenteral nutrition support, and carnitine supplementation. It is difficult to distinguish based on its concentration alone, thus the recall rate is too high.

[0004] Enzyme deficiencies not only alter the concentrations of their direct substrates and products but also affect the concentrations of upstream and downstream metabolites in the metabolic network. Existing cutoff value methods and multivariate statistical methods treat each metabolite as an independent variable, failing to utilize the correlation information within the metabolic reaction network. Existing methods use disease-related metabolites as seed nodes, performing a restarted random walk on the metabolic reaction network to rank candidate metabolites according to steady-state probability, but do not determine whether the candidate metabolites can be measured in dried blood spots. The top-ranked candidates are mostly intracellular intermediates such as coenzyme A thioesters; these molecules do not enter peripheral blood and have poor stability, making them undetectable in dried blood spots.

[0005] Furthermore, triple quadrupole mass spectrometry, in multiple reaction monitoring mode, only records signals for preset ion pairs, does not retain full-scan spectra, and operates at unit mass resolution, making it impossible to extract chromatographic peaks of unpreset metabolites at precise mass. For newly recommended candidate biomarkers, newborns need to be recalled for re-blood collection or a separate detection method needs to be established; when including new diseases in screening, biomarkers for these diseases also need to be validated separately and cutoff values ​​established.

[0006] This invention relates to eight types of inherited metabolic diseases, covering four categories: organic acidemia, fatty acid oxidation disorders, urea cycle disorders, and lysosomal storage diseases. Type I glutamate hyperacidity is caused by a deficiency of glutaryl-CoA dehydrogenase, resulting in elevated serum glutarylcarnitine (C5DC) and urinary glutaric acid and 3-hydroxyglutaric acid. Affected children often develop acute encephalopathy in infancy due to infection or fasting, leaving residual basal ganglia damage. Methylmalonic acidemia is caused by a deficiency of methylmalonic acid-CoA mutase or cobalamin metabolism disorders, leading to the accumulation of methylmalonic acid, propionylcarnitine, and methylcitric acid. Approximately 85% of affected children develop the disease before the age of one. Ornithine carbamoyltransferase deficiency is an X-linked urea cycle disorder; obstruction of the ammonia detoxification pathway leads to hyperammonemia. Neonatal onset can rapidly progress to cerebral edema. Very long chain acyl-CoA dehydrogenase deficiency is caused by mutations in the ACADVL gene, with the main marker being tetradecenoylcarnitine (C14:1). The cardiopathic form can lead to fatal multiple organ failure even in the neonatal period. Medium chain acyl-CoA dehydrogenase deficiency is caused by mutations in the ACADM gene, with elevated capryloylcarnitine (C8) as the main marker. Carriers are at risk of metabolic crisis under infection or fasting stress. Pompe disease is caused by mutations in the GAA gene leading to acid α-glucosidase deficiency. The infantile form presents with myocardial hypertrophy, hypotonia, and respiratory failure; the effectiveness of enzyme replacement therapy depends on early diagnosis. Fabry disease is caused by mutations in the GLA gene leading to α-galactosidase A deficiency, with progressive accumulation of trihexosylceramide and its deacylated product, trihexosphosamine. Mucopolysaccharidosis type I is caused by a mutation in the IDUA gene leading to a deficiency in α-L-iduronidase, resulting in the accumulation of dermatin sulfate and heparan sulfate in lysosomes. Once irreversible neurological and skeletal damage occurs, the effectiveness of treatment is significantly reduced. All eight diseases mentioned above have a narrow treatment window, and early diagnosis can significantly alter the prognosis. Summary of the Invention

[0007] The purpose of this invention is to provide a method and system for constructing a genetic metabolic disease screening model based on metabolic network enhancement, so that the extended metabolites recommended by the metabolic network can be obtained from the collected test data without the need for re-collection of blood or additional injections.

[0008] To achieve the above objectives, this invention provides a method for constructing a genetic metabolic disease screening model based on metabolic network enhancement, comprising the following steps: After adding isotope internal standards to dried blood spot samples from newborns in the training set, liquid chromatography-mass spectrometry analysis was performed. The content of each metabolite in the preset target metabolite group was determined by the internal standard method as the initial feature. High-resolution full-scan spectra covering the entire chromatographic retention time range and data-independent acquisition spectra were collected and stored as digital sample files for the sample. A metabolic reaction network is constructed with metabolites as nodes and directed edges pointing from substrates to products to represent enzyme-catalyzed reactions. The probability transition matrix is ​​obtained by normalizing the edges using the product of carbon atom conservation factor, reaction competition factor, and reaction confidence factor. For each target disease, the differential metabolites between positive samples and control samples of the disease are used as seed nodes. Random walks with restart are performed on the metabolic reaction network along the substrate-to-product direction and along the product-to-substrate direction to obtain the forward propagation score and backward propagation score of each node. The same calculation is performed with all nodes having an equal probability initial distribution to obtain the global baseline forward score and global baseline backward score. The product of the scores in the two directions after correction by the global baseline score is used as the comprehensive recommendation score. Several non-seed nodes with the highest comprehensive recommendation score are selected as recommendation nodes. The predicted retention time of the metabolite corresponding to the recommended node is corrected by the measured retention time of the isotope internal standard in the same sample. The extracted ion chromatographic peak of the metabolite is extracted from the digital sample archive and integrated within the corrected retention time window to obtain the recovery feature. The initial features and the supplementary features are concatenated into an enhanced feature matrix, and the classification model for each target disease is trained.

[0009] Furthermore, the liquid chromatography-mass spectrometry analysis employs high-resolution mass spectrometry, alternately performing targeted quantitative acquisition, full scan, and data-independent acquisition within the same chromatographic run, with the mass resolution of full scan and data-independent acquisition being no less than 15,000.

[0010] Furthermore, the reaction competitiveness factor is the reciprocal of the number of organic product types in the reaction, and the reaction credibility factor is determined according to the reaction annotation evidence level; the restart probability of the random walk is determined by grid search with the area under the average receiver operating characteristic curve of each target disease on the validation set as the objective function.

[0011] Furthermore, after selecting recommended nodes, the measurability of dried blood spots of each recommended node is determined based on the biological sample attribution annotations and molecular physicochemical properties in the metabolite database. Unmeasurable recommended nodes are mapped to measurable proxies in body fluids along transport reaction edges or binding reaction edges, and the proxies are used as the metabolites corresponding to the recommended nodes.

[0012] Furthermore, a piecewise linear mapping is established using the predicted and measured retention times of the internal standards of each isotope within the same sample as anchor points to correct the predicted retention times of the metabolites corresponding to the recommended nodes.

[0013] Furthermore, the forward propagation score and the backward propagation score determine the upstream and downstream positions of each feature in the enhanced feature matrix relative to the seed node and the expected direction of concentration change. The directional consistency score of the sample is calculated by combining the reliability weight of each feature, which is used as the input feature of the classification model, and a consistency threshold is determined. The supplementary features with reliability weights lower than the preset threshold are removed.

[0014] Furthermore, the classification model is a binary classifier trained on a gradient boosting decision tree as the base learner and using a one-to-many strategy for each target disease.

[0015] This invention also provides a system for constructing a genetic metabolic disease screening model based on metabolic network enhancement, including a dual-channel acquisition and archiving module, a metabolic network management module, a bidirectional propagation module, a complementation module, and a joint judgment and reporting module. The joint judgment and reporting module outputs the risk probability of the newborn to be assessed for each target disease using a classification model constructed according to the above method. The system may also include a measurability mapping module and a consistency scoring module, which jointly determine the risk level using risk probability and directional consistency score.

[0016] Compared with the prior art, the beneficial effects of the present invention are: 1. This invention provides a method and system for constructing a genetic metabolic disease screening model based on metabolic network enhancement. During liquid chromatography-mass spectrometry analysis, full scan and data-independent acquisition spectra are simultaneously acquired and stored. The extended metabolites recommended by the metabolic network can be directly extracted from the digital sample archive without recalling newborns for blood collection or establishing separate detection methods for the extended metabolites. When adding new screening diseases, the constructed metabolic response network and the stored digital sample archive can be reused.

[0017] 2. This invention provides a method and system for constructing a genetic metabolic disease screening model based on metabolic network enhancement. It uses the product of carbon atom conservation factor, response competition factor and response reliability factor as edge weights, and corrects the forward and backward propagation scores with global baseline scores, thereby reducing the proportion of highly connected cofactors such as ATP and coenzyme A in the recommendation results.

[0018] 3. This invention provides a method and system for constructing a genetic metabolic disease screening model based on metabolic network enhancement. It determines the measurability of recommended nodes in dried blood spots, mapping intracellular intermediates such as coenzyme A thioesters to surrogates measurable in dried blood spots. The predicted retention time is corrected by the measured retention time of the isotope internal standard, reducing missed extractions and peak misalignment within the extraction window. The expected direction of changes in the concentration of each characteristic is determined by forward and backward propagation scores, and a direction consistency score is calculated, which can be used to distinguish samples with consistent changes in upstream and downstream metabolites in the pathway from samples with only a few elevated markers. Attached Figure Description

[0019] To more clearly illustrate the specific embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the specific embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained from these drawings without creative effort.

[0020] Figure 1 This is a schematic diagram of the overall process of the method of the present invention.

[0021] Figure 2 This is a timing diagram for dual-channel acquisition.

[0022] Figure 3 This is a schematic diagram of the construction of the three-factor weights of a directed edge.

[0023] Figure 4 This is a schematic diagram of bidirectional propagation and global baseline correction.

[0024] Figure 5 This is a schematic diagram of measurability filtering and proxy mapping.

[0025] Figure 6 A schematic diagram of the internal standard anchorage retention time mapping correction.

[0026] Figure 7 This is a schematic diagram comparing the expected direction vector with the measured deviation direction. The left image shows true positive samples, and the right image shows confirmed false positive samples.

[0027] Figure 8 This is a comparison chart of the individual orientation consistency score distribution between true positive samples and confirmed false positive samples.

[0028] Figure 9 This is a flowchart illustrating the process of determining risk under two conditions.

[0029] Figure 10 A comparison chart of the average AUC of eight diseases under different extended feature selection strategies.

[0030] Figure 11 This is a schematic diagram of the module structure of the system of the present invention.

[0031] The diagram shows: Dual-channel acquisition and archiving module 110, metabolic network management module 120, bidirectional propagation module 130, measurability mapping module 140, backfilling module 150, consistency scoring module 160, and joint judgment and reporting module 170. Detailed Implementation

[0032] To make the objectives, technical solutions, and advantages of the present invention clearer, the technical solutions are described in detail below with reference to the accompanying drawings and specific embodiments. The described embodiments are some, but not all, of the embodiments of the present invention. Where specific conditions are not specified in the embodiments, they are performed under conventional conditions in the art or conditions recommended by the manufacturer; reagents, consumables, chromatographic columns, mass spectrometers, etc., whose manufacturers are not specified, are all conventional products that can be obtained through commercial channels.

[0033] The screening model construction method and screening assessment system of the present invention can be used according to... Figure 1 The process is as shown. Model training in steps S1 to S8 and step S9 is performed on the training set to build the screening model. For newborns to be evaluated, the screening assessment system obtains the consistency score of their enhancement features and orientation according to steps S1, S6, and S8, and determines the risk level according to step S9.

[0034] Step S1: Take a perforated slide of dried blood spot sample, extract it with an extractant containing an isotopic internal standard, and then perform liquid chromatography-mass spectrometry (LC-MS) analysis. The content of each metabolite in the targeted metabolite group is determined using the internal standard method as the initial characteristic. Simultaneously, high-resolution full-scan chromatograms and data-independent acquisition chromatograms covering the entire chromatographic retention time range are acquired and stored in an open format as a digital sample file for the sample. When using high-resolution mass spectrometry, targeted quantitative acquisition, full scan, and data-independent acquisition are performed within the same chromatographic run. Figure 2 The timing sequence is executed alternately, with a full scan and data-independent acquisition quality resolution of no less than 15,000. The acquisition cycle is set to ensure that each chromatographic peak is acquired at least 12 times in both types of acquisition. For samples already quantified using triple quadrupole mass spectrometry, a second perforated slide of the same dried blood spot can be acquired on a high-resolution mass spectrometer and stored as a digital sample file.

[0035] Step S2: Construct a directed weighted metabolic reaction network G=(V,E,W), where V is the set of metabolite nodes, E is the set of directed edges, and W is the set of edge weights; a directed edge e=(u,v) indicates that metabolite u can be converted into metabolite v via an enzymatic reaction. Before construction, remove inorganic small molecule nodes with a molecular weight below 50 Da. Figure 3 The calculation of the edge weights for each directed edge is shown below: ; In the formula, It is the carbon atom conservation factor. As a competitive factor for the response, The reaction confidence factor is used. The carbon atom conservation factor and the reaction competition factor are calculated using the following formulas: ; ; In the formula, The total number of carbon atoms in metabolite u. For the reaction from u to v The number of carbon atoms retained from u to v. The carbon skeleton of u is fully preserved at v. The value is 1.0; when u participates in the reaction only as a cofactor and does not contribute carbon atoms to v (including the case where u does not contain carbon atoms). We set it to 0.01 to preserve network connectivity. For the reaction A collection of organic products, To disregard the number of organic products resulting from the reaction with protons and water, reversible reactions are counted in both directions. Reaction reliability factor. The values ​​are assigned based on the level of evidence for the response annotation: 1.0 for responses that have been manually verified, 0.75 for responses that have been cross-validated across multiple databases, and 0.5 for responses that are calculated and predicted from only a single database.

[0036] The three factors are combined in a product. Taking a bypass reaction involving a cofactor as an example, its... It is 0.01. and All three factors have a value of 1.0, resulting in a product of 0.01. The arithmetic mean of these three factors is 0.67, indicating that the cofactor edge still holds a significant weight in the probability transition matrix. Random walks can transfer between unrelated pathways via highly connected nodes such as ATP and coenzyme A. When using a product combination, a lower value for any factor will reduce the weight of that edge. A comparison of the two combination methods is shown in Table 6.

[0037] Normalize the edge weights by row to obtain the probability transition matrix P, whose elements... for Divide by the sum of the edge weights of all outgoing edges of node u; when a node has no outgoing edges, its self-loop transition probability is set to 1.

[0038] Step S3: For each target disease d, compare each feature of the positive samples and control samples for that disease in the training set. The compared features include initial features and ratio-type features calculated from the initial features. Differential metabolites that simultaneously satisfy both the fold change threshold and the statistical significance threshold are used as seed nodes, forming a seed node set S(d). Ratio-type differential features are mapped to network nodes according to the metabolites corresponding to their molecules. After mapping, duplicates with existing seed nodes are merged. Ratio-type features can use metabolite concentration ratios commonly used in neonatal tandem mass spectrometry screening, such as C3 / C0, C3 / C2, C4 / C0, and C4 / C3. Seed nodes are determined from samples with known diagnostic results in the training set and do not depend on the samples to be evaluated.

[0039] Step S4: Construct the initial distribution vector using S(d). The element of its corresponding seed node is taken All other elements are 0. Number of seed nodes. (By...) Figure 4 The process involves bidirectional propagation, with forward and backward propagation iterated using the following formulas: ; ; In the formula, t is the iteration round. and Let be the forward propagation score vector and the backward propagation score vector for the t-th iteration, respectively. This is the transpose of P. For the probability of restarting, .by During iteration, the score shifts along the substrate to product direction, and during iteration P, it backtracks along the product to substrate direction. The iteration termination condition for forward propagation is: ; In the formula, The Euclidean norm is given; the iteration termination condition for backpropagation is given. Since P is a row random matrix and , The spectral radii are all The iterations of both equations converge to a unique fixed point. Restart probability. The objective function is determined by grid search, using the average AUC of each target disease on the validation set as the objective function (if the validation set does not contain positive samples of a certain disease, that disease is not included in the average).

[0040] Furthermore, using an initial distribution with equal probability for all nodes, the same forward and reverse calculations are performed to obtain the global baseline forward score for each node. Inverse scoring from global baseline .node The overall recommendation score is calculated using the following formula: ; In the formula, and These are the corresponding nodes in f and b, respectively. elements, To prevent constants with zero denominators, cofactors such as ATP, coenzyme A, and NADH participate in numerous reactions and have high scores in both f and b without correction, but are unrelated to specific diseases; the global baseline scores of such nodes are also high, but their overall recommendation scores decrease after ratio correction. When using the product form, only nodes with strong correlations to the seed node in both the forward and reverse directions receive high overall recommendation scores. The top K non-seed nodes are selected as recommendation nodes in descending order of their overall recommendation scores.

[0041] Step S5, press Figure 5As shown, based on the biological sample attribution annotations in the metabolite database, the presence of blood or plasma detection records, and molecular physicochemical properties, including whether the sample is in a coenzyme-bound state, whether it carries a permanent charge, and the stability of ester or thioester bonds, the measurability of dried blood spots at each recommended node is determined. Recommended nodes meeting any of the following conditions are determined to be unmeasurable dried blood spots: they are in a coenzyme-bound state such as coenzyme A thioester; the molecule contains chemical bonds such as thioester bonds that are easily hydrolyzed during the drying and storage of dried blood spots; there are no blood or plasma detection records in the metabolite database, and the spots contain permanently charged groups or phosphate groups that hinder their transmembrane entry into peripheral blood at physiological pH. All other recommended nodes are determined to be measurable.

[0042] Nodes identified as intracellularly confined, unstable, or not entering peripheral blood are mapped along transport or conjugation reaction edges to fluid-measurable surrogates, which are then used as the metabolites corresponding to the recommended node. Acyl-CoA thioesters are mapped to their hydrolyzed free acid, homoacylglycine, or homoacylcarnitine. For dicarboxylic acid acyl-CoA thioesters, the free acid is used; for the others, acylglycine is preferred because it is more stable than short-chain acylcarnitine in dried blood spots and is unaffected by carnitine supplementation. Phosphorylated sugar metabolic intermediates are mapped to their dephosphorylated products or downstream oligosaccharides. If multiple candidate surrogates still exist according to the above rules, the edge with the largest weight connecting the original node and the candidate surrogate is selected. If the surrogate is already in the targeted metabolite group, no additional feature is generated.

[0043] Step S6: For the metabolites corresponding to the recommended nodes, predict the retention time based on their chemical structures, and then... Figure 6 The method shown uses in-sample isotope internal standard correction: For each sample, using the predicted and measured retention times of all isotope internal standards as anchor points, a piecewise linear mapping is established in ascending order of predicted values. The predicted retention time of the metabolite is transformed into a correction value. When the predicted retention time is outside the range of the first and last anchor points, it is extrapolated using the straight line determined by the two nearest anchor points. The extraction window is centered on the correction value, and its half-width is three times the standard deviation of the residuals of the two adjacent anchor points on both sides of the correction value (when the correction value is outside the anchor point range, the two nearest anchor points are used). The anchor point residual is the difference between the measured retention time and the predicted retention time of the internal standard. Within the extraction window, the extracted ion chromatographic peak of the metabolite is extracted from the digital sample archive according to the exact mass number and integrated. The area of ​​the isotope internal standard peak with the closest retention time is normalized to obtain the recovery feature. Recovery is completed on the stored digital sample archive without re-collecting or re-injecting samples.

[0044] The recovered features are semi-quantitative results, and their reliability is characterized by reliability weights: ; In the formula, The reliability weights for compensating feature i; For the signal-to-noise ratio normalization term, take the features The ratio of the median signal-to-noise ratio in all samples to 10 is taken as 1 if it is greater than 1; To retain the time deviation decay term, take , The median of the ratio of the absolute value of the difference between the measured peak retention time and the correction value in each sample to the half-width of the extraction window; As a penalty term for the coefficient of variation of quality control, the minimum value of the coefficient of variation of all compensated features in the mixed quality control samples is taken, along with the feature value. The ratio of coefficients of variation. The reliability weight of the initial feature is set to 1. Complementary features with a reliability weight lower than a preset threshold, such as 0.3, are discarded, and the remaining complementary features are concatenated with the initial features to form an enhanced feature matrix.

[0045] Step S7, for each target disease d and each feature i in the enhanced feature matrix, calculate the topological location quantity using the following formula: ; In the formula, , , and The value of the node corresponding to feature i is taken. The node corresponding to the initial feature is the node for that metabolite, and the node corresponding to the supplementary feature is the recommended node from which it originated. When enzyme deficiency blocks the metabolic pathway, metabolites accumulate upstream of the blockage point, while downstream metabolites decrease. A value greater than 0 indicates a feature It is primarily accessed via reverse propagation, located upstream of the seed node, and its expected sign is +1; Less than 0 indicates a feature Located downstream of the seed node, the expected sign is -1. The expected signs of each feature constitute the directional expected vector s(d) of disease d.

[0046] Step S8: Using the reference population (normal control sample, stratified by blood collection date and birth weight) as a baseline, calculate the robust z-values ​​for each characteristic of the sample. ,in For sample features The measured value, and These are the median and interquartile range of feature i in the respective stratum. The directional consistency score is calculated using the following formula: ; In the formula, the summation covers all features of the enhanced feature matrix. For symbolic functions, The weighting coefficients are determined by the following formula:

[0047] The larger the absolute value of the topological location quantity, the more clearly the upstream and downstream positions of the feature are defined, and the larger its weight coefficient. For complementary features with low reliability weights, the weight coefficients decrease accordingly. In enzyme-deficient samples, upstream features tend to increase while downstream features tend to decrease, and the measured deviation direction is mostly consistent with the expected sign. The values ​​are relatively large positive values. Elevated C3 levels caused by non-disease factors such as maternal vitamin B12 deficiency only involve a few markers; upstream markers such as propionylglycine and isobutyrylglycine do not change accordingly with downstream markers such as succinic acid, and the positive and negative values ​​cancel each other out. Approaching 0, comparison method shows opinion Figure 7 .

[0048] Training samples It can be used as one of the input features for a classification model. This involves comparing positive and negative samples of the disease in the training set. The distribution determines the consistency threshold, taking the cut point with the largest Youden index; negative samples are prioritized from those that initially screened positive but were subsequently excluded from the disease, as indicated in the opinion. Figure 8 .

[0049] Step S9: Using the enhanced feature matrix as input, or the enhanced feature matrix and the orientation consistency score as input, train a classification model for each target disease using a one-to-many strategy. For the newborn to be evaluated, the classification model outputs the risk probability p of each target disease; p is not lower than the risk probability threshold and A risk level is defined as high if the value is not lower than the consistency threshold, and p is defined as high if the value is not lower than the risk probability threshold. If the value is below the consistency threshold, it is determined that a review is required; if p is below the risk probability threshold, it is determined that the risk is low. See the determination process below. Figure 9 .

[0050] The modular structure of the screening and assessment system is as follows: Figure 11 As shown, the system includes a dual-channel acquisition and archiving module 110, a metabolic network management module 120, a bidirectional propagation module 130, a measurability mapping module 140, a backfilling module 150, a consistency scoring module 160, and a joint judgment and reporting module 170. Module 110 executes step S1; module 120 stores and maintains the network and probability transition matrix constructed in step S2; module 130 executes steps S3, S4, and S7; module 140 executes step S5; module 150 executes step S6; module 160 executes step S8; and module 170 executes model training and risk assessment in step S9 and generates a screening report.

[0051] As a specific embodiment 1, this embodiment constructs and validates screening models for eight inherited metabolic diseases. Samples were obtained from routine screening samples and confirmatory test archives from three provincial neonatal screening centers from January 2019 to October 2024. All dried blood spot samples were collected within 24 to 72 hours after birth, dried at room temperature for at least 4 hours, and then sealed and stored at 4°C. The research protocol was reviewed and approved by the ethics committees of all participating institutions, and the legal guardians of all participating newborns signed written informed consent forms.

[0052] Inclusion criteria were newborns with complete initial screening records and confirmatory test results. Exclusion criteria were samples with obvious hemolysis of blood spots, blood spot diameter less than 6 mm, or damaged and damp packaging. A total of 856 cases were ultimately included and randomly assigned in a 7:1.5:1.5 ratio to form a training set of 600 cases, a validation set of 128 cases, and an independent test set of 128 cases. The cohort composition is shown in Table 1.

[0053] Table 1

[0054] Among them, non-target disease IEM positive refers to samples diagnosed with other inherited metabolic diseases such as phenylketonuria, maple syrup urine disease, propionic acidemia, and tyrosinemia, but not the eight target diseases of this invention. These samples are treated as negative classes in the training of each binary classifier to enhance the model's ability to distinguish the similarity of metabolic characteristics between diseases. Normal controls are defined as samples with normal initial screening results and no clinical symptoms related to inherited metabolic diseases during at least 12 months of follow-up.

[0055] Process the samples according to step S1. Using a 3.2 mm diameter punch, take a circular filter paper from the center of the dried blood spot collection circle and place it in a 2 mL polypropylene 96-well deep plate. Place a blank filter paper at each of the four corners and the center of each plate as a blank control for extraction. Add 100 μL of extraction working solution to each well. This working solution is prepared by mixing 90 mL of methanol, 10 mL of ultrapure water, and 0.1 mL of formic acid, and then adding the stock solutions of each isotope. The final concentration of each internal standard is set in the range of 0.2 to 2.0 μmol / L. Incubate the deep plate at 600 rpm for 30 minutes at room temperature, centrifuge at 2500×g for 5 minutes, transfer 80 μL of supernatant to the sample plate, dry it under a nitrogen stream at 45°C, and reconstitute it with 100 μL of reconstitution solution (water:methanol = 80:20, containing 0.1% formic acid). Take equal volumes of each reconstitution solution from all the test samples and mix them as a pooled quality control sample.

[0056] The chromatographic column was a C18 reversed-phase column (ACQUITY UPLC HSS T3, 2.1 × 50 mm, 1.8 μm, Waters Corporation). Mobile phase A was an aqueous solution containing 0.1% formic acid, and mobile phase B was an acetonitrile solution containing 0.1% formic acid. The gradient elution program was as follows: 0 to 0.5 min, maintain 0% B; 0.5 to 2.5 min, B phase linearly increases from 0% to 95%; 2.5 to 3.0 min, maintain 95% B; 3.0 to 3.5 min, B phase linearly decreases from 95% to 0%; 3.5 to 5.0 min, maintain 0% B equilibration. The flow rate was 0.5 mL / min, the column temperature was 40℃, and the injection volume was 5 μL.

[0057] Mass spectrometry was performed using a high-resolution mass spectrometer with full scan and data-independent acquisition capabilities, an electrospray ionization source, and positive / negative ion switching. Within the same chromatographic run, [the following was done]: Figure 2 The timing sequence shown alternates between the two types of data acquisition.

[0058] The quantitative channel collected data on 82 targeted metabolites according to the preset precursor ion mass and retention time window using parallel reaction monitoring or targeted selected ion monitoring, with a resolution of 17500. Each metabolite was quantified using the internal standard method.

[0059] The archiving channel uses full scan, m / z 50 to 1000, resolution 60000, combined with data-independent acquisition, isolation window width 20 Da, covering m / z 50 to 1000, resolution 15000, without a target list.

[0060] The two types of acquisitions are performed alternately within the same acquisition cycle, with the acquisition cycle set to ensure that each chromatographic peak is acquired at least 12 times in both types of acquisitions. Archived data is stored in the mzML open format and saved together with the sample number, acquisition batch, and instrument parameters to form a digital sample archive.

[0061] The targeted metabolites in this embodiment were selected based on the following criteria: they possess known biomarker value for the target disease; they can be directly quantified using the aforementioned reversed-phase chromatography method without derivatization, except for glycosaminoglycan-derived disaccharides; they are stable in dried blood spots; and commercially available isotope internal standards are available. A total of 82 metabolites were selected, belonging to five classes. The 22 amino acids are alanine, arginine, aspartic acid, asparagine, citrulline, glutamic acid, glutamine, glycine, histidine, leucine, isoleucine, alloleucine, lysine, methionine, ornithine, phenylalanine, proline, serine, threonine, tryptophan, tyrosine, and valine. Leucine, isoleucine, and alloleucine were quantified separately after being separated using the aforementioned gradient baseline. There are 28 types of acylcarnitines, namely C0, C2, C3, C4, C5, C5:1, C5DC, C5OH, C6, C8, C8:1, C10, C10:1, C12, C12:1, C14, C14:1, C14:2, C14OH, C16, C16:1, C16OH, C18, C18:1, C18:2, C18OH, C18:1OH and C4DC.

[0062] Sixteen organic acids were identified: methylmalonic acid, succinic acid, glutaric acid, 3-hydroxyglutaric acid, 2-hydroxyglutaric acid, ethylmalonic acid, methylsuccinic acid, adipic acid, octanoic acid, sebacic acid, 3-hydroxypropionic acid, orotic acid, 4-hydroxyphenyllactic acid, methylcitric acid, acetoacetic acid, and 2-methyl-3-hydroxybutyric acid. Eight sphingolipids and polysaccharides were identified: lysphingomyelin, trihexosphosin, glucosphosin, galactosphosin, dermatan sulfate-derived disaccharide, heparan sulfate-derived disaccharide HS-NS and HS-OS, and glucose tetrasaccharide. Eight auxiliary biomarkers were identified: creatine, guanidinoacetic acid, dimethylglycine, homocysteine, homocysteine, uridine, kynurenine, and serotonin. Dermatin sulfate-derived disaccharides and heparan sulfate-derived disaccharides, namely HS-NS and HS-OS, are highly polar and show almost no retention on a reversed-phase C18 column. Their determination is performed as follows: The dried blood spot extract is pretreated with methanol or enzymatic digestion to break the glycosaminoglycan polysaccharide chains into detectable disaccharide units, which are then injected separately and separated on a HILIC column (2.1 × 100 mm, 1.7 μm). The mobile phase is 10 mM ammonium acetate buffer (pH 5.0) and acetonitrile. Alternatively, reversed-phase ion-pair chromatography using dibutylamine as the ion-pairing reagent can be used. Mass spectrometry monitors [MH] in negative ion mode. - Or [M=-2H] 2- Ions were quantified using the internal standard method.

[0063] The precursor ion mass, characteristic fragment ions, collision energy, and retention time of each metabolite were pre-determined through standard injection. Quantification was performed using the internal standard method. For the 16 metabolites without specific isotopic internal standards, the average of the two internal standards with the closest retention times was normalized. Quantitative results were uniformly expressed in μmol / L, yielding an 82-dimensional initial metabolic feature vector.

[0064] The network was constructed according to step S2. Human metabolic reaction data were extracted from four public databases: KEGG, HMDB, MetaCyc, and Recon3D. Metabolites were deduplicated and merged using the first-level hash code of the InChI Key. Reactions were represented as reactant-product pairs, and reversible reactions were split into two directed edges, one forward and one reverse, and weighted accordingly. Nodes of inorganic small molecules with a molecular weight below 50 Da and nodes recorded only in a single database without cross-validation evidence were removed. The SMILES structures of each metabolite were obtained through the PubChem and ChEBI interfaces, and nodes for which no clear structure could be obtained were removed. The total number of carbon atoms was statistically analyzed using RDKit to resolve the SMILES structures, and the number of carbon atoms retained in the product was determined by the RXNMapper chemical reaction atom mapping algorithm. The reaction confidence factor was assigned according to the database evidence code: 1.0 for those from KEGG RPAIR or Reactome with manually verified or literature evidence, 0.75 for those from Rhea with cross-validation by two or more databases, and 0.5 for those calculated and predicted only from a single database. The resulting network contains thousands of metabolite nodes and tens of thousands of directed edges, with non-zero elements accounting for approximately 0.04%. It is stored in a sparse matrix format and occupies approximately 28 MB of memory.

[0065] Determine the seed nodes according to step S3. For each target disease, perform the Mann-Whitney U test on each of the 82 initial features and ratio features of the confirmed positive samples and normal control samples in the training set, and calculate the concentration fold change. Select the features that simultaneously satisfy |log2(FC)|>0.585, i.e., FC>1.5 or FC<0.667 and p<0.01, and select the top 10 in descending order of |log2(FC)|. Map them to network nodes and merge duplicates to use as seed nodes.

[0066] Taking methylmalonic acidemia as an example, the differential characteristics included: propionylcarnitine C3, FC=9.8, p=3.2×10 -12 Methylmalonic acid, FC=15.1, p=8.9×10 -16 Glycine, FC=2.9, p=7.1×10 -7 Free carnitine C0, FC=0.35, p=5.3×10 -6 Acetylcarnitine C2, FC=1.9, p=2.8×10 -4Methylcitric acid, FC=4.6, p=4.3×10 -8 3-Hydroxypropionic acid, FC=3.0, p=1.2×10 -5 The two ratios, C3 / C0 and C3 / C2, were also observed, both of which were mapped to and merged with the propionylcarnitine C3 node. Due to the limited number of positive samples for Pompeii etiology, the first five differential features were selected: glucose tetrasaccharide, FC=12.3, p=1.8×10⁻⁶. -9 Creatine, FC=2.3, p=2.5×10 -4 Aspartic acid, FC=1.8, p=3.6×10 -3 Free carnitine C0, FC=0.68, p=1.5×10 -2 And the C4 / C0 ratio characteristic, which is mapped to the butyrylcarnitine C4 node on a molecular basis.

[0067] Perform bidirectional propagation according to step S4. Perform a grid search on the restart probability ρ among 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.40, and 0.50, with the average AUC of each disease on the validation set as the objective function. Since there are no MPS I positive samples in the validation set, this disease is not included in the average, and the optimal value of 0.20 is taken. Take 10 -8 The iterations converged after 50 to 120 iterations. Taking methylmalonic acidemia as an example, before baseline correction, cofactors and inorganic molecules accounted for 9 of the top 20 recommended items and 45 of the top 100, or 45.0%. After correction, these nodes were all ranked outside the top 100. In the recommended list of eight diseases after baseline correction, the average proportion of cofactors was 3.1%, which is consistent with the values ​​of the present invention listed in Table 6.

[0068] For each disease, the top 20 non-seed nodes with the comprehensive recommendation score were selected, K=20. The union of the recommended nodes for each disease was then taken, resulting in a total of 83 nodes. Taking methylmalonic acidemia as an example, the top three in the comprehensive recommendation score were succinyl-CoA, L-methylmalonyl-CoA, and isobutyryl-CoA, with adipic acid monoacyl-CoA at the seventh position. All of these are intracellular coenzyme A thioesters, which do not enter the plasma and whose thioester bonds are easily hydrolyzed, making them undetectable in dried blood spots. The detectability of the recommended nodes and surrogate mapping were performed according to step S5. The results of the top 9 recommended nodes in the comprehensive recommendation score for this disease are shown in Table 2.

[0069] Table 2

[0070] Among them, the recommended nodes and their rankings are the network propagation results of this embodiment. The rankings are obtained by sorting according to the comprehensive recommendation scores after removing seed nodes. The testability determination and proxy mapping are performed according to the rules described in step S5.

[0071] Perform the backfilling in step S6. For the metabolites corresponding to each recommended node, obtain their precise monoisotope masses and major adduct forms from HMDB and PubChem, [M+H]. + or[M+Na] + In negative ion mode, it is [MH]. - The predicted retention time was obtained using the Retip retention time prediction model based on the SMILES structure. Ion chromatographic peaks were extracted within the extraction window at an exact mass value ± 5 ppm. After Gaussian smoothing and baseline correction, the peaks were integrated, and the average of the peak areas of the two internal standard isotopes whose retention time most closely matched the corrected value of the metabolite was used as a normalized reference to obtain a semi-quantitative relative concentration.

[0072] In this embodiment, after mapping in step S5, the corresponding metabolites of the 83 recommended nodes that are already in the target metabolite group do not generate additional complementary features; those mapped to the same metabolite are merged; and the remaining metabolites are complemented in the digital sample archive. A retention condition is a signal-to-noise ratio (SNR) of at least 3 for all chromatographic peaks; those with an SNR below this standard are discarded, resulting in 68 effective complementary features. The 82 initial features and 68 complementary features are concatenated into a 150-dimensional enhanced feature matrix.

[0073] Train the classification model according to step S9. A one-to-many strategy is used, training an independent gradient boosting decision tree binary classifier for each target disease. The positive class represents confirmed positive samples for that disease, and the negative class is the union of positive samples from other diseases and normal controls. The input features are a 150-dimensional augmented feature matrix. The hyperparameter search domain is: number of decision trees 100 to 800, stride 50, learning rate 0.005 to 0.2, maximum depth 3 to 10, row subsampling ratio 0.6 to 1.0, and column subsampling ratio 0.5 to 1.0. Seven-fold hierarchical cross-validation combined with Bayesian optimization is used to determine the hyperparameters. The surrogate model is a Gaussian process, and the acquisition function is the expected improvement criterion. A total of 50 iterations are performed, with the first 10 iterations using random sampling to initialize the search space.

[0074] The 7-fold stratified cross-validation uses "disease category + positive / negative" as the stratification unit, namely, six units: normal control, positive IEM for non-target diseases, positive for organic acidemia, positive for fatty acid oxidation disorders, positive for urea cycle disorders, and positive for lysosomal storage diseases, rather than stratifying by individual disease labels. When the number of positive samples for a certain disease in the training set is less than the number of folds, the positive samples for that disease are not included in the proportional split, but are allocated to different folds one by one according to the order determined by a fixed random seed, so that the validation part of each fold contains at most one positive sample for that disease; when the training part or validation part of a fold does not contain a positive sample for that disease, that fold is not included in the calculation of the AUC for that disease, and the cross-validation target value for that disease is the average of the AUCs of the remaining folds. For diseases with only one positive sample in the training set, such as MPS I, there is no fold from which AUC can be calculated. The hyperparameters are not optimized separately. Instead, the optimal hyperparameters determined by cross-validation from Pompe disease or Fabry disease with a larger sample size in the same disease category are used, or the preset default hyperparameters are adopted, namely, 200 decision trees, learning rate 0.05, maximum depth 4, row subsampling ratio 0.8, and column subsampling ratio 0.8. The performance of this disease is evaluated by the results of the independent test set combined with the bootstrap confidence interval.

[0075] The incidence rates of the eight target diseases varied by tens of times, resulting in severe class imbalance in the training set. Therefore, the class weights of each binary classifier were set according to the inverse of the class frequency; hyperparameter optimization used AUC rather than accuracy as the objective function, as AUC is insensitive to the prior distribution of classes; model evaluation also primarily used AUC as the metric. In a classifier with a positive-to-negative ratio of approximately 1:600, predicting all samples as negative yielded an accuracy of over 99.8% with zero sensitivity; therefore, accuracy was not used as the evaluation metric.

[0076] The risk probability threshold was set at the point on the ROC curve of the training set with the largest Youden exponent, with an initial reference value of 0.50. After hyperparameter optimization and training on the training set, performance was evaluated on an independent test set, and the results are shown in Table 3.

[0077] Table 3

[0078] 1000 bootstrap resampling cycles with replacement were performed on the independent test set. The 95% confidence intervals of the AUC for each disease were calculated using the percentile method. The interval widths for all eight diseases did not exceed 0.08. For example, the AUC for MMA was 0.982, and the 95% confidence intervals ranged from 0.953 to 0.998. Twenty independent 7-fold stratified cross-validation cycles were performed on the training set. Except for MPS I, the standard deviation of the AUC for each disease was less than 0.008.

[0079] In this embodiment, ablation experiments were conducted under the same conditions. Table 4 shows the comparison between the AUC of each disease and the enhancement feature matrix when all reinforcing features were removed and only 82 initial features were retained.

[0080] Table 4

[0081] The diseases with the largest increases in AUC were VLCADD, Pompe, and MPS I, for which there are few or no initial biomarkers available.

[0082] The 68 recommended metabolites were replaced with 68 randomly selected metabolites from HMDB and 68 metabolites selected based on their frequency of appearance in literature, respectively. Training and evaluation were performed under the same conditions. The randomized group was repeated 10 times, and the average was taken. The results are shown in Table 5. Figure 10 .

[0083] Table 5

[0084] The three groups of extended features had the same number of features. Random selection and selection based on literature frequency increased the average AUC by 0.002 and 0.007, respectively, while network recommendation increased the average AUC by 0.024. This indicates that the increase in AUC mainly depends on the selection method of extended metabolites and is not significantly related to the increase in the number of features.

[0085] Using the average AUC of eight diseases and the proportion of cofactors in the recommended list as indicators, cofactors include ATP, coenzyme A, NADH, NADPH, FAD, FMN, etc., the comparison of the five edge weight schemes is shown in Table 6.

[0086] Table 6

[0087] The proportion of auxiliary factors in the unweighted network recommendation list is 23.8%, decreasing to 12.5% ​​when only the carbon atom conservation factor is used, and further decreasing to 3.1% in the three-factor product scheme, a reduction of 87.0% compared to the unweighted network. The proportion of auxiliary factors in the three-factor addition scheme is 11.2%, higher than the 8.9% in the two-factor multiplication scheme. This is because when the carbon atom conservation factor in the additive combination is 0.01, the other two factors still maintain a high weight for the auxiliary factor edge. The three-factor addition scheme and the product scheme use the same factors, yet the proportion of auxiliary factors differs by approximately 3.6 times, indicating that the reduction in the proportion of auxiliary factors mainly depends on the product combination method.

[0088] Under the same conditions, the comparison between the baseline correction scheme and the two-way scoring fusion method is shown in Tables 7 and 8, respectively.

[0089] Table 7

[0090] Normalization by node out-degree is less effective than global baseline correction because out-degree only reflects the number of local connections of a node, while global baseline score also reflects the impact of the node's location in the network and reachable paths on the steady-state distribution of random walks.

[0091] Table 8

[0092] When using only unidirectional propagation scoring, the average AUC is lower than that of bidirectional fusion. Forward propagation reflects downstream products of differentially metabolites, while backward propagation reflects their upstream precursors. Additive fusion cannot distinguish between nodes with high scores in one direction and low scores in the other, and nodes with moderate scores in both directions; maximum value fusion ignores the score in the other direction; geometric mean fusion does not correct for the influence of node connectivity.

[0093] As a specific embodiment 2, this embodiment uses short-chain acyl-CoA dehydrogenase deficiency as an example to illustrate the process of adding a new target disease.

[0094] Eighteen confirmed positive samples of this disease were extracted from the training set, and differential metabolite analysis was performed between them and normal controls. The selected differentially expressed features included butyrylcarnitine C4, FC=5.8, p=4.2×10⁻⁶. -8 Ethylmalonic acid, FC=9.1, p=7.5×10 -11 The C4 / C3 ratio feature was also used, which was mapped to the butyrylcarnitine C4 node and merged with it. Using the metabolic response network and probability transition matrix constructed in Example 1, bidirectional propagation was performed with this seed set, taking approximately 3 minutes. Features for the metabolites corresponding to newly recommended nodes were retrieved from the stored digital sample archive, merged with the existing enhanced feature matrix, and then used to train a binary classifier for the disease.

[0095] The classifier achieved an AUC of 0.963 on the independent test set. The entire process, from differential metabolite statistics to classifier training and evaluation, took 1 to 2 hours. The network construction was completed in Example 1, during which no new samples were collected or additional samples were injected.

[0096] As a specific embodiment 3, laboratories that have already adopted the triple quadrupole multiple reaction monitoring method for quantification can complete the quantification with triple quadrupole and the archiving with high-resolution mass spectrometry according to this embodiment.

[0097] Two 3.2 mm perforated slides were taken from the same dried blood spot. The first slide was used for quantification of 82 targeted metabolites on a triple quadrupole using conventional multiple reaction monitoring (MRM). The second slide was acquired and archived within the same analytical batch using a full scan combined with data-independent acquisition on a high-resolution mass spectrometry platform. Both slides were taken from the same blood spot, without increasing the number of blood samples collected. Subsequent steps S2 to S9 were the same as in Example 1, and were performed on the digital sample archive obtained from the second slide.

[0098] In this embodiment, the chromatographic and mass spectrometric conditions for the triple quadrupole quantitative channel are as follows: the column, mobile phase, and gradient are the same as in Example 1; the mass spectrometer uses an electrospray ionization source with switching between positive and negative ion modes, multiple reaction monitoring (MRM), capillary voltage of 3.0 kV, ion source temperature of 150 °C, desolvation gas temperature of 500 °C, desolvation gas flow rate of 900 L / h, cone gas flow rate of 50 L / h, and residence time for each channel of 5 to 20 ms. Method validation data for representative metabolites are shown in Table 9. Method validation was performed according to CLSI C62-A guidelines and WS / T 492-2016: the linear range was determined by calibration curves with at least 8 concentration levels, with each level repeated 3 times, using weighted least squares regression, requiring a correlation coefficient of not less than 0.990; the limit of detection and limit of quantitation were determined with a signal-to-noise ratio of not less than 3 and not less than 10, respectively; the recovery rate was calculated by adding three concentration levels of standards (low, medium, and high) to a blank dried blood spot matrix, with 6 replicates per level; intra-batch precision was determined by 20 consecutive injections of the same quality control sample within a single batch, and inter-batch precision was determined by the results of measurements over 20 consecutive working days.

[0099] Table 9

[0100] Where CE represents collision energy, RT represents retention time, LOD represents limit of detection, LOQ represents limit of quantitation, and CV represents coefficient of variation. LOD and LOQ units are within the same linear range. Recovery rates are expressed as mean ± SD, representing the combined statistical results for low, medium, and high concentration levels. Ion-pair parameters for organic acids were obtained in negative ion mode, and the capillary voltage was optimized to -2.5 kV for negative ion mode. Data for GAbiose were obtained on a HILIC column (2.1 × 100 mm, 1.7 μm) in negative ion mode, with a mobile phase of 10 mM ammonium acetate buffer (pH 5.0) and acetonitrile.

[0101] When implementing the high-resolution dual-channel method described in Example 1, the above-mentioned methodological parameters of the quantitative channel need to be re-verified on this platform, and the verification items and acceptance criteria are the same as above.

[0102] In the above embodiments, the system suitability was assessed for each analytical batch according to the following procedure: At the beginning of the batch, blank solvent and extraction blank were injected twice each to confirm the absence of interference residues; standard quality control samples were injected to confirm that the retention time drift of each channel was less than 0.05 min and the peak area variation was less than 15%; mixed quality control samples were continuously injected at least 5 times to balance the column and ion source, and then the batch was run at a frequency of one mixed quality control sample inserted for every 12 analytes; at the end of the batch, standard quality control samples were injected again to assess instrument performance drift. A peak area variation coefficient of each metabolite in the mixed quality control sample not exceeding 15% was considered acceptable; for those exceeding 15% but not exceeding 25%, local weighted regression smoothing correction was performed using the response values ​​of the three mixed quality control samples before and after that metabolite; for those exceeding 25%, the metabolite was not reported quantitatively in that batch.

[0103] The screening report output by the screening assessment system includes: basic information of the newborn; risk probability, directional consistency score and risk level of each target disease; for diseases determined to be high-risk or requiring re-examination, the top 5 initial features and top 5 supplementary features with the greatest contribution to the classification and their contribution, and indicate whether the measured deviation direction of each feature is consistent with the expected direction; for high-risk diseases, confirmatory testing suggestions are given, such as plasma amino acid analysis, urine organic acid analysis, lysosomal enzyme activity measurement or targeted gene sequencing.

[0104] The above-described specific embodiments are merely preferred embodiments of the present invention and are not intended to limit the scope of protection of the present invention. Various modifications, substitutions, and improvements made by those skilled in the art to the technical solutions of the present invention based on the provided textual description and drawings, without departing from the design concept and spirit of the present invention, should all fall within the scope of protection of the present invention. The scope of protection of the present invention is determined by the claims.

Claims

1. A method for constructing a genetic metabolic disease screening model based on metabolic network enhancement, characterized in that, Includes the following steps: S1. After adding isotope internal standards to dried blood spot samples from each newborn in the training set, liquid chromatography-mass spectrometry analysis was performed. The content of each metabolite in the preset target metabolite group was determined by the internal standard method as the initial feature. High-resolution full-scan spectrum and data-independent acquisition spectrum covering the entire chromatographic retention time range were collected and stored as the digital sample file of the sample. S2. Construct a metabolic reaction network with metabolites as nodes and directed edges from substrate to product to represent enzyme-catalyzed reactions. Use the product of carbon atom conservation factor, reaction competition factor and reaction confidence factor as edge weights. The carbon atom conservation factor is the proportion of the number of carbon atoms in the substrate that are retained in the product to the total number of carbon atoms in the substrate. Normalize the edge weights by row to obtain the probability transition matrix P. S3. For each target disease, using the differential metabolites between positive samples and control samples of the disease as seed nodes, perform a random walk with restart along the substrate to product direction on the metabolic reaction network to obtain the forward propagation score of each node, and perform a random walk with restart along the product to substrate direction to obtain the backward propagation score of each node. S4. Perform the same calculation with all nodes having an equal initial distribution to obtain the global baseline positive score and global baseline negative score for each node. S5. The comprehensive recommendation score is the product of the ratio of the forward propagation score to the global baseline forward score and the ratio of the backward propagation score to the global baseline backward score. The K non-seed nodes with the highest comprehensive recommendation scores are selected as recommendation nodes, where K is a preset positive integer. S6. For each recommended node, the predicted retention time of the metabolite is corrected by the measured retention time of the isotope internal standard in the same sample. The extracted ion chromatographic peak of the metabolite is extracted from the digital sample file within the corrected retention time window and integrated to obtain the recovery feature. S7. The initial features and the supplementary features are concatenated into an enhanced feature matrix, and the enhanced feature matrix is ​​used to train a classification model for each target disease.

2. The method according to claim 1, characterized in that, The liquid chromatography-mass spectrometry analysis employs high-resolution mass spectrometry, alternately performing targeted quantitative acquisition, full scan, and data-independent acquisition within the same chromatographic run. The mass resolution of the high-resolution full scan spectrum and the data-independent acquisition spectrum is not less than 15,000.

3. The method according to claim 1, characterized in that, Edge weight of the directed edge from metabolite u to metabolite v Calculate using the following formula: ； ； ； In the formula, , and These are the carbon atom conservation factor, the reaction competition factor, and the reaction credibility factor, respectively. Values ​​are assigned based on the level of evidence indicated in the response annotations. The total number of carbon atoms in metabolite u. For the reaction from u to v The number of carbon atoms retained by u up to v, when u does not contribute carbon atoms to v. Take 0.01; For the reaction A collection of organic products, This refers to the number of organic product types that exclude protons and water.

4. The method according to claim 1, characterized in that, The forward propagation scores of each node constitute the forward propagation score vector f, and the backward propagation scores of each node constitute the backward propagation score vector b. f and b are calculated iteratively according to the following formula: ； ； In the formula, t is the iteration round. and Let f and b be the values ​​obtained in the t-th iteration. This is the transpose of P. The restart probability is determined by grid search; This is the initial distribution vector, where the element corresponding to the seed node is the reciprocal of the number of seed nodes, and all other elements are 0. The iterations of f and b have Euclidean norms less than 10 for the difference between two adjacent iterations. -9 The process will end at that time.

5. The method according to claim 4, characterized in that, node Overall recommendation rating Calculate using the following formula: ； In the formula, and These are the elements of node i in f and b, respectively. and These represent the global baseline positive score and the global baseline negative score for node i, respectively. To prevent constants with a denominator of zero.

6. The method according to claim 1, characterized in that, After selecting recommended nodes, the following also includes: The measurability of dried blood spots at each recommended node is determined based on the biological sample attribution annotations and molecular physicochemical properties in the metabolite database. Unmeasurable recommended nodes are mapped to measurable proxies in body fluids along transport reaction edges or binding reaction edges, and these proxies are used as the metabolites corresponding to the recommended nodes. Among them, acyl-CoA thioester is mapped to its free acid obtained from hydrolysis, acyl-glycine with the same acyl group, or acyl-carnitine with the same acyl group, and phosphorylated sugar metabolic intermediate is mapped to its dephosphorylated product or downstream oligosaccharide. The targeted metabolite group includes amino acids, acylcarnitines, organic acids, sphingolipids and polysaccharides, and cofactor marker metabolites; the target diseases include one or more of the following: glutaric acidemia type I, methylmalonic acidemia, ornithine carbamoyltransferase deficiency, very long chain acyl-CoA dehydrogenase deficiency, medium chain acyl-CoA dehydrogenase deficiency, Pompe disease, Fabry disease, and mucopolysaccharidosis type I; the classification model is a binary classifier trained for each target disease using a gradient boosting decision tree as the base learner and a one-to-many strategy.

7. The method according to claim 1, characterized in that, When correcting the predicted retention time, a piecewise linear mapping is established using the predicted retention time and the measured retention time of each isotope internal standard in the same sample as anchor points, and the predicted retention time is transformed using this piecewise linear mapping.

8. The method according to claim 5, characterized in that, It also includes: calculating the directional consistency score for each training sample for each target disease d. The features are used as input features to the classification model, and the positive and negative samples of the disease in the training set are used as input features. The distribution determines the consistency threshold; Calculate using the following formula: ； ； ； In the formula, n is the number of features of the enhanced feature matrix, and the summation covers all features of the enhanced feature matrix. ; For samples in features The robust z-value relative to the reference population; It is a symbolic function; For topological location quantities, , , and Take features The value of the corresponding node, the node corresponding to the feature is the recommendation node from which it originated; For weighting coefficients; For reliability weights, the initial feature is set to 1, and the complemented feature is the product of the signal-to-noise ratio normalization term, the retention time deviation attenuation term, and the quality control variation coefficient penalty term; Eliminate features whose reliability weight is below a preset threshold.

9. A system for constructing a genetic metabolic disease screening model based on metabolic network enhancement, used to implement the method for constructing a genetic metabolic disease screening model based on metabolic network enhancement as described in any one of claims 1-8, characterized in that, It includes a dual-channel acquisition and archiving module (110), a metabolic network management module (120), a bidirectional propagation module (130), a recovery module (150), and a joint judgment and reporting module (170). The dual-channel acquisition and archiving module (110) is used to acquire initial features and store digital sample archives; The metabolic network management module (120) is used to store the metabolic reaction network and the probability transition matrix; The bidirectional propagation module (130) is used to select recommended nodes; The backfilling module (150) is used to extract backfilling features and splice enhanced features; The joint determination and reporting module (170) is used to output the risk probability of the newborn to be assessed having each target disease using the classification model constructed according to the method and to generate a screening report.

10. The system according to claim 9, characterized in that, It also includes a measurability mapping module (140) for mapping unmeasurable recommendation nodes to measurable proxies of bodily fluids, and a method for calculating the orientation consistency score of the newborn to be evaluated according to claim 8. Consistency scoring module (160); The joint judgment and reporting module (170) is in operation when the risk probability is not lower than a preset risk probability threshold and A risk level is determined to be high if the risk probability is not lower than the consistency threshold described in claim 8, and if the risk probability is not lower than the risk probability threshold. If the value is below the consistency threshold, it is determined that a review is required.

Images