A method for constructing a genome-scale metabolic network model

Abstract

This invention belongs to the field of systems biology and discloses a method for constructing a genome-scale metabolic network model. The method includes: constructing a biochemical reaction dataset, a metabolite dataset, and a gene-protein-reaction rule dataset; performing gene function annotation and screening; mapping biochemical reactions, constructing a biochemical reaction set, and adding relevant metabolite dataset information; mapping local biomass equations corresponding to similar species and adding them to the biochemical reaction set to generate an initial genome-scale metabolic network model; and sequentially performing gap filling, redundant reaction removal, and annotation information addition on the initial genome-scale metabolic network model to obtain the final genome-scale metabolic network model. This invention addresses the shortcomings of existing GEMs construction tools in terms of model coverage and prediction accuracy, develops an automated de novo construction process from genome sequences to GEMs, reduces time-consuming manual correction processes, and improves construction efficiency.

Description

Technical Field

[0001] This invention relates to the field of systems biology, and in particular to a method for constructing a genome-scale metabolic network model. Background Technology

[0002] Genome-scale metabolic network models (GEMs) are mathematical frameworks for describing the metabolic capacity of organisms. They are built based on genome annotation, metabolic databases, and experimental data. They describe all metabolic responses of a target organism and their gene-protein-reaction (GPR) rules through stoichiometry matrices and use linear programming methods such as flux balance analysis (FBA) to simulate cellular metabolic growth.

[0003] GEMs contain complete metabolic information about an organism, mainly including metabolites, genes, and reactions. Metabolites include substrates, intermediates, and products in the cell. For ease of representation, they are usually represented using metabolite IDs, with different formats for each database. Metabolites in GEMs need to include compartmentalization information to indicate the cellular location where the reaction occurs. For example, "_c" represents the cytoplasm, "_e" represents the extracellular environment, "_m" represents mitochondria, and "glc__D_c" represents D-glucose in the cytoplasm. GEMs also include name information, i.e., the common name of the metabolite. For example, the common name of the metabolite corresponding to D-glucose is D-Glucose. GEMs record the molecular formula information of the metabolites to represent their molecular composition. For example, the molecular formula corresponding to D-glucose is "...". GEMs also contain charge information, representing the net charge of the metabolite under physiological conditions; typically, "glc__D_c" has a charge of 0. Furthermore, GEM annotations provide cross-reference information for the metabolite in multiple databases such as BiGG, KEGG, and MetaCyc. For example, D-glucose is identified as "glc__D" in BiGG and has the ID "C00031" in KEGG.

[0004] Currently, GEMs have become an important tool in metabolic engineering, providing guidance and improving research efficiency. Due to the large amount of data contained in metabolic models, manually constructing GEMs is time-consuming and labor-intensive, leading to a growing demand for automated methods and tools for constructing genome-scale metabolic network models. In recent years, genome annotation and metabolic databases have developed rapidly, and various automated reconstruction tools have been developed. However, traditional tools support a very limited number of species for model construction; a few tools support prokaryotic models, and only a handful can construct eukaryotic GEMs, but the quality of these models needs improvement. This is mainly because the tools rely on single databases, have limited species coverage, and lack response data for specific species. Eukaryotic metabolism is complex, with numerous organelles; therefore, model construction requires extensive cellular sublocalization information to determine reaction compartments, as well as sufficient transport reactions and specialized construction strategies.

[0005] Furthermore, research on the construction tools revealed that GEMs built by these tools typically do not support accurate quantitative simulations. The primary reason is the presence of errors in the database, such as unlimited energy generation and net metabolite production. These errors lead to non-conservation of mass, ultimately affecting quantitative calculations. Secondly, issues with the tool's workflow and algorithms necessitate manual intervention to adjust the constructed models for accurate quantitative calculations.

[0006] Therefore, the main technical challenges and difficulties in constructing eukaryotic models are the limited availability of databases, the complexity of eukaryotic model construction, and the need for accurate quantitative calculations. How to provide a method for constructing genome-scale metabolic network models is an urgent problem to be solved. Summary of the Invention

[0007] This invention provides a method for constructing a genome-scale metabolic network model to address the problems mentioned above in the prior art.

[0008] According to a first aspect of the present invention, a method for constructing a genome-scale metabolic network model is provided.

[0009] In one embodiment, the method for constructing the genome-scale metabolic network model includes: Biochemical reaction data from biological metabolism databases are integrated and preprocessed. Based on the preprocessed biochemical reaction data, biochemical reaction datasets, metabolite datasets, and gene-protein-reaction rule datasets are constructed. A local alignment library is constructed based on the protein sequences of the target species. The similarity between the input sequence and the sequences in the local alignment library is compared, and gene function annotation is performed based on the alignment results. The alignment score threshold of the input sequence is determined based on the dynamic threshold method, and the gene function annotation is screened in combination with the alignment score threshold. Based on the selected gene function annotation results, biochemical reaction mapping is performed by combining the biochemical reaction dataset and the gene-protein-reaction rule dataset. Based on the mapping results, a biochemical reaction set is constructed, and relevant metabolite dataset information is added to the biochemical reaction set. Species with the most similar metabolism to the input sequence are selected, the local biomass equations corresponding to the similar species are mapped and added to the biochemical reaction set, and the added biochemical reaction set is integrated into the metabolite file. An initial genome-scale metabolic network model is generated based on the metabolite file. The initial genome-scale metabolic network model was sequentially processed by gap filling, redundant response removal, and annotation information addition to obtain the final genome-scale metabolic network model.

[0010] In one embodiment, the integration of biochemical reaction data from a biological metabolism database, followed by preprocessing of the biochemical reaction data, and the construction of a biochemical reaction dataset, a metabolite dataset, and a gene-protein-reaction rule dataset based on the preprocessed biochemical reaction data, includes: Biochemical reaction data from biological metabolism databases are integrated, and the biochemical reaction data are sequentially cleaned and standardized, error corrected, and processed for respiratory chain reactions. Based on the processed biochemical reaction data, biochemical reaction datasets for prokaryotes, eukaryotes, and archaea are constructed. Extract metabolite information for each biochemical reaction from the biochemical reaction dataset and remove duplicates. Construct a metabolite dataset based on the deduplicated metabolite information. Extract gene-protein-reaction rule data from the biochemical reaction dataset and obtain the correspondence between reactions and gene IDs. Based on the correspondence between reactions and gene IDs, construct a gene-protein-reaction rule dataset.

[0011] In one embodiment, the metabolite information includes metabolite ID, name, molecular formula, charge, compartment information, and annotation information.

[0012] In one embodiment, the steps of constructing a local alignment library based on the protein sequences of the target species, performing similarity comparisons between the input sequence and sequences in the local alignment library, and annotating gene functions based on the alignment results; determining the alignment score threshold of the input sequence based on a dynamic threshold method, and filtering gene function annotations in conjunction with the alignment score threshold include: The protein sequences of the target species were obtained, and text-based bioinformatics standard format files for prokaryotes, eukaryotes and archaea were constructed respectively. The text-based bioinformatics standard format files were then converted into the special format files of sequence alignment software to build a local alignment library. The system receives input sequences uploaded by users, automatically identifies the type of input sequences using built-in functions and scripts of the sequence alignment software, calls the amino acid sequence or nucleotide sequence model in the sequence alignment software to perform sequence similarity comparison with the local alignment library, outputs the sequence alignment result file, and completes gene function annotation. The sequence alignment result file is checked and ensured to exist and not be empty. The annotation results are sorted in descending order according to the alignment score of the sequence alignment result file, and duplicate gene matches are removed. Based on the sorting distribution of alignment scores, a dynamic threshold method is used to determine the alignment score threshold of the input sequence, and gene function annotation results with alignment scores lower than the alignment score threshold are filtered out.

[0013] In one embodiment, determining the alignment score threshold of the input sequence using a dynamic thresholding method based on the ranking distribution of alignment scores includes: The total number of genes after deduplication is counted, a preset percentage of the total number of genes is selected as the threshold number of genes, and the threshold number of genes is stored in the critical value variable; The total number of genes falling into the preset similarity score interval is counted and the sum is calculated. Based on the calculation results, the range of intervals into which the critical value falls is determined. Determine whether the range into which the critical value falls belongs to the preset similarity score range. If so, use the lower bound of the range as the comparison score threshold; otherwise, select the preset threshold as the comparison score threshold.

[0014] In one embodiment, the step of mapping biochemical reactions based on the screened gene function annotation results, combined with biochemical reaction datasets and gene-protein-reaction rule datasets, constructing a biochemical reaction set based on the mapping results, and adding relevant metabolite dataset information to the biochemical reaction set includes: Gene IDs are extracted from the sequence alignment results, and gene IDs whose scores meet the threshold are searched in the gene-protein-reaction rule dataset to obtain the corresponding reaction IDs and protein combinations. Based on the reaction ID, match rows with the same ID column in the biochemical reaction dataset and extract the entire row of data; repeat the search and extraction operation for all gene IDs that meet the threshold, and summarize to obtain the biochemical reaction set, which includes reaction ID, reaction name, metric matrix, upper and lower boundaries of reaction, gene-protein-reaction rules, annotation information, target reaction and reaction equation; Based on the result with the highest alignment score, find all subunits of catalytic reactions in the gene-protein-reaction rule dataset, search for the compartment where the subunit is located in the Uniprot database and replace it with the correct compartment to complete the compartment mapping of the effective subunit; By employing or linking isozymes and connecting them to complexes composed of multiple subunits, the Boolean connection logic can be reconstructed. By combining compartment information, biochemical reaction datasets, and gene-protein-reaction rule datasets, we extract the catalytic metabolic reactions of gene-annotated enzymes and complete the adaptation of rules to reactions. Add the effective subunit and complex related dataset information after compartment mapping to the biochemical reaction set.

[0015] In one embodiment, after adding relevant metabolite dataset information to the biochemical reaction set, the method further includes adding spontaneous reactions, cellular sublocalization information of the reactions, and culture medium reactions to the biochemical reaction set.

[0016] In one embodiment, gap-filling of an initial genome-scale metabolic network model includes: The sequence alignment scores are divided into several confidence levels, and a weighting factor is assigned to the response corresponding to each confidence level. The gaps in the initial genome-scale metabolic network model are filled based on the weighting factors.

[0017] In one embodiment, redundant response removal processing of the initial genome-scale metabolic network model includes: Traverse all reactions in the initial genome-scale metabolic network model, remove circulating metabolites in the reactions, and identify redundant reaction sets with the same substrates and products. Remove nested subsets from redundant reaction groups and retain the optimal reaction in each group based on mass balance, reversibility, and metabolite quantity indices; Eliminate reactions containing adenosine triphosphate and adenosine diphosphate, as well as empty reactions, from redundant reaction groups; Remove reactions that were not retained from the redundant reaction group and merge the gene rules of the unretained reactions into the retained reactions; Identify and remove reversible reactions that contain only adenosine triphosphate and adenosine diphosphate, as well as completely repeatable reaction formulas; Based on the results of flux variability analysis, unnecessary redundant reactions associated with individual genes are cyclically removed, and metabolites that do not participate in the reaction are cleaned up, so as to achieve redundant reaction removal processing of the initial genome-scale metabolic network model.

[0018] In one embodiment, adding annotation information to the initial genome-scale metabolic network model includes: At the gene level, add the correspondence between gene and protein sequence and function databases; At the reaction level, information related to metabolic reactions and enzyme catalysis is added.

[0019] According to a second aspect of the present invention, a system for constructing a genome-scale metabolic network model is provided.

[0020] In one embodiment, the genome-scale metabolic network model construction system includes: The dataset construction module is used to integrate biochemical reaction data from the biological metabolism database, preprocess the biochemical reaction data, and construct biochemical reaction datasets, metabolite datasets, and gene-protein-reaction rule datasets based on the preprocessed biochemical reaction data. The gene function annotation module is used to construct a local alignment library based on the protein sequences of the target species, perform similarity comparison between the input sequence and the sequences in the local alignment library, and perform gene function annotation based on the alignment results; based on the dynamic threshold method, the alignment score threshold of the input sequence is determined, and the gene function annotation is screened in combination with the alignment score threshold. The biochemical reaction mapping module is used to perform biochemical reaction mapping based on the screened gene function annotation results, combined with the biochemical reaction dataset and the gene-protein-reaction rule dataset. Based on the mapping results, a biochemical reaction set is constructed, and relevant metabolite dataset information is added to the biochemical reaction set. The initial model building module is used to screen out species that are most metabolically similar to the input sequence, map the local biomass equations corresponding to the similar species and add them to the biochemical reaction set, and integrate the added biochemical reaction set into the metabolite file to generate an initial genome-scale metabolic network model based on the metabolite file. The model optimization module is used to sequentially perform gap filling, redundant response removal, and annotation information addition on the initial genome-scale metabolic network model to obtain the final genome-scale metabolic network model.

[0021] The technical solutions provided by the embodiments of the present invention may include the following beneficial effects: This invention integrates biochemical reactions from the BiGG, KEGG, and MetaCyc databases, performs quality control and error correction on the reaction datasets, and organizes them into three biochemical reaction datasets: prokaryotes, eukaryotes, and archaea. Gene function annotation is performed using the DIAMOND alignment tool, and a dynamic threshold method is developed to filter more reasonable alignment results, solving the value selection problem caused by differences in genome size between species and improving the accuracy of gene-protein-reaction rules. Gap-filling and optimization of the model are performed using methods such as simplified throughput balance analysis, weighting strategies, and confidence scoring, ultimately generating GEMs in multiple formats. This effectively fills the gaps in model coverage and prediction accuracy of existing GEMs construction tools. A de novo automated construction process from genome sequences to GEMs is developed, reducing time-consuming manual correction processes and improving construction efficiency.

[0022] It should be understood that the above general description and the following detailed description are exemplary and explanatory only, and are not intended to limit the invention. Attached Figure Description

[0023] The accompanying drawings, which are incorporated in and form part of this specification, illustrate embodiments consistent with the invention and, together with the description, serve to explain the principles of the invention.

[0024] Figure 1 This is a flowchart illustrating the construction of a genome-scale metabolic network model according to an exemplary embodiment; Figure 2 This is a diagram illustrating a metabolite data structure according to an exemplary embodiment; Figure 3 This is a schematic diagram illustrating GPR rule extraction according to an exemplary embodiment; Figure 4 This is a diagram illustrating the data structure of the comparison results according to an exemplary embodiment; Figure 5 This is a diagram illustrating the reaction set data structure according to an exemplary embodiment; Figure 6 This is a diagram illustrating the GPR rule data structure according to an exemplary embodiment; Figure 7 This is a flowchart illustrating Biomass assignment according to an exemplary embodiment. Detailed Implementation

[0025] The following description and accompanying drawings fully illustrate specific embodiments described herein to enable those skilled in the art to practice them. Some portions and features of certain embodiments may be included in or replace portions and features of other embodiments. The scope of the embodiments herein includes the entire scope of the claims and all available equivalents thereof. The various embodiments described herein are presented in a progressive manner, with each embodiment focusing on its differences from other embodiments; similar or identical parts between embodiments can be referred to interchangeably.

[0026] The modules in the apparatus or system of this application can be implemented entirely or partially through software, hardware, or a combination thereof. These modules can be embedded in or independent of the processor in a computer device in hardware form, or stored in the memory of a computer device in software form, so that the processor can call and execute the operations corresponding to each module.

[0027] Where there is no conflict, the embodiments and features in the embodiments of the present invention can be combined with each other.

[0028] Figure 1 An embodiment of a genome-scale metabolic network model construction method of the present invention is shown.

[0029] In this optional embodiment, the method for constructing the genome-scale metabolic network model includes: Step S101: Integrate biochemical reaction data from the biological metabolism database, preprocess the biochemical reaction data, and construct a biochemical reaction dataset, a metabolite dataset, and a gene-protein-reaction rule dataset based on the preprocessed biochemical reaction data. Step S102: Construct a local alignment library based on the protein sequences of the target species, perform similarity comparison between the input sequence and the sequences in the local alignment library, and perform gene function annotation based on the alignment results; determine the alignment score threshold of the input sequence based on the dynamic threshold method, and filter gene function annotations in combination with the alignment score threshold. Step S103: Based on the filtered gene function annotation results, perform biochemical reaction mapping by combining the biochemical reaction dataset and the gene-protein-reaction rule dataset. Based on the mapping results, construct a biochemical reaction set and add relevant metabolite dataset information to the biochemical reaction set. Step S104: Select the species that are most metabolically similar to the input sequence, map the local biomass equation corresponding to the similar species and add it to the biochemical reaction set, and integrate the added biochemical reaction set into the metabolite file. Generate an initial genome-scale metabolic network model based on the metabolite file. Step S105: Perform gap filling, redundant reaction removal and annotation information addition on the initial genome-scale metabolic network model in sequence to obtain the final genome-scale metabolic network model.

[0030] In this optional embodiment, the integration of biochemical reaction data from the biological metabolism database, the preprocessing of the biochemical reaction data, and the construction of a biochemical reaction dataset, a metabolite dataset, and a gene-protein-reaction rule dataset based on the preprocessed biochemical reaction data include: Biochemical reaction data from biological metabolism databases are integrated, and the biochemical reaction data are sequentially cleaned and standardized, error corrected, and processed for respiratory chain reactions. Based on the processed biochemical reaction data, biochemical reaction datasets for prokaryotes, eukaryotes, and archaea are constructed. Extract metabolite information for each biochemical reaction from the biochemical reaction dataset and remove duplicates. Construct a metabolite dataset based on the deduplicated metabolite information. Extract gene-protein-reaction rule data from the biochemical reaction dataset and obtain the correspondence between reactions and gene IDs. Based on the correspondence between reactions and gene IDs, construct a gene-protein-reaction rule dataset.

[0031] In this optional embodiment, the metabolite information includes metabolite ID, name, molecular formula, charge, compartment information, and annotation information.

[0032] In this optional embodiment, the steps of constructing a local alignment library based on the protein sequences of the target species, performing similarity comparisons between the input sequence and sequences in the local alignment library, and performing gene function annotation based on the alignment results; determining the alignment score threshold of the input sequence based on the dynamic threshold method, and filtering gene function annotations in conjunction with the alignment score threshold include: The protein sequences of the target species were obtained, and text-based bioinformatics standard format files for prokaryotes, eukaryotes and archaea were constructed respectively. The text-based bioinformatics standard format files were then converted into the special format files of sequence alignment software to build a local alignment library. The system receives input sequences uploaded by users, automatically identifies the type of input sequences using built-in functions and scripts of the sequence alignment software, calls the amino acid sequence or nucleotide sequence model in the sequence alignment software to perform sequence similarity comparison with the local alignment library, outputs the sequence alignment result file, and completes gene function annotation. The sequence alignment result file is checked and ensured to exist and not be empty. The annotation results are sorted in descending order according to the alignment score of the sequence alignment result file, and duplicate gene matches are removed. Based on the sorting distribution of alignment scores, a dynamic threshold method is used to determine the alignment score threshold of the input sequence, and gene function annotation results with alignment scores lower than the alignment score threshold are filtered out.

[0033] In this optional embodiment, determining the alignment score threshold of the input sequence using a dynamic thresholding method based on the ranking distribution of alignment scores includes: The total number of genes after deduplication is counted, a preset percentage of the total number of genes is selected as the threshold number of genes, and the threshold number of genes is stored in the critical value variable; The total number of genes falling into the preset similarity score interval is counted and the sum is calculated. Based on the calculation results, the range of intervals into which the critical value falls is determined. Determine whether the range into which the critical value falls belongs to the preset similarity score range. If so, use the lower bound of the range as the comparison score threshold; otherwise, select the preset threshold as the comparison score threshold.

[0034] In this optional embodiment, the step of performing biochemical reaction mapping based on the screened gene function annotation results, combined with the biochemical reaction dataset and the gene-protein-reaction rule dataset, constructing a biochemical reaction set based on the mapping results, and adding relevant metabolite dataset information to the biochemical reaction set includes: Gene IDs are extracted from the sequence alignment results, and gene IDs whose scores meet the threshold (i.e., the alignment result threshold when the total number of genes is greater than 75% of the total number of annotated genes) are searched in the gene-protein-reaction rule dataset to obtain the corresponding reaction IDs and protein combinations. Based on the reaction ID, match rows with the same ID column in the biochemical reaction dataset and extract the entire row of data; repeat the search and extraction operation for all gene IDs that meet the threshold, and summarize to obtain the biochemical reaction set, which includes reaction ID, reaction name, metric matrix, upper and lower boundaries of reaction, gene-protein-reaction rules, annotation information, target reaction and reaction equation; Based on the result with the highest alignment score, find all subunits of catalytic reactions in the gene-protein-reaction rule dataset, search for the compartment where the subunit is located in the Uniprot database and replace it with the correct compartment to complete the compartment mapping of the effective subunit; By employing or linking isozymes and connecting them to complexes composed of multiple subunits, the Boolean connection logic can be reconstructed. By combining compartment information, biochemical reaction datasets, and gene-protein-reaction rule datasets, we extract the catalytic metabolic reactions of gene-annotated enzymes and complete the adaptation of rules to reactions. Add the effective subunit and complex related dataset information after compartment mapping to the biochemical reaction set.

[0035] In this optional embodiment, after adding relevant metabolite dataset information to the biochemical reaction set, the method further includes adding spontaneous reactions, cellular sublocalization information of the reactions, and culture medium reactions to the biochemical reaction set.

[0036] In this optional embodiment, gap-filling processing of the initial genome-scale metabolic network model includes: The sequence alignment scores are divided into several confidence levels, and a weighting factor is assigned to the response corresponding to each confidence level. The gaps in the initial genome-scale metabolic network model are filled based on the weighting factors.

[0037] In this optional embodiment, the redundant response removal process for the initial genome-scale metabolic network model includes: Traverse all reactions in the initial genome-scale metabolic network model, remove circulating metabolites in the reactions, and identify redundant reaction sets with the same substrates and products. Remove nested subsets from redundant reaction groups and retain the optimal reaction in each group based on mass balance, reversibility, and metabolite quantity indices; Eliminate reactions containing adenosine triphosphate and adenosine diphosphate, as well as empty reactions, from redundant reaction groups; Remove reactions that were not retained from the redundant reaction group and merge the gene rules of the unretained reactions into the retained reactions; Identify and remove reversible reactions that contain only adenosine triphosphate and adenosine diphosphate, as well as completely repeatable reaction formulas; Based on the results of flux variability analysis, unnecessary redundant reactions associated with individual genes are cyclically removed, and metabolites that do not participate in the reaction are cleaned up, so as to achieve redundant reaction removal processing of the initial genome-scale metabolic network model.

[0038] In this optional embodiment, the annotation information addition process for the initial genome-scale metabolic network model includes: At the gene level, add the correspondence between gene and protein sequence and function databases; At the reaction level, information related to metabolic reactions and enzyme catalysis is added.

[0039] To facilitate understanding of the above technical solutions of the present invention, the following further explains the above technical solutions of the present invention from the perspective of architecture and principle, as follows: I. Dataset Construction 1. Construction of biochemical reaction dataset Data Integration and Collection: Addressing the issues of existing tools' heavy reliance on single databases and insufficient data comprehensiveness, this study integrates biochemical reactions from three databases—BiGG, KEGG, and MetaCyc—based on the rigorously quality-controlled BKMD database. Furthermore, the data was expanded using the latest biochemical reaction data from the KEGG and MetaCyc databases. This comprehensive and up-to-date data improves the accuracy of the final GEM construction.

[0040] The BiGG data were primarily extracted from Universalmodel v1.6 on the website (http: / / bigg.ucsd.edu / data_access), containing a total of 28,301 reactions from 108 published high-quality GEMs representing 35 species. These included GEMs from model organisms such as *Escherichia coli*, *Bacillus subtilis*, and *Saccharomyces cerevisiae*. Table 1 lists the major species (excluding subspecies) among these 35 species.

[0041] Table 1 Species in the BiGG model

[0042] KEGG reaction data are from its website. MetaCyc reaction data are extracted from SBML models published on its website (https: / / metacyc.org / download.shtml), version 28.0. Among them, 11,443 prokaryotic models have UniProt IDs, meaning the reactions have corresponding protein sequences; therefore, these were downloaded and integrated, resulting in 42,988 reactions. There are 222 archaea models, of which 170 have UniProt ID annotations; these were integrated, resulting in 28,183 reactions. There are 39 eukaryotic models, of which 17 have UniProt annotations; these were integrated, resulting in 10,822 reactions. Table 2 shows examples of 5 prokaryotic, 5 archaea, and 5 eukaryotic species from MetaCyc.

[0043] Table 2 Species representation in the MetaCyc model

[0044] Due to significant metabolic differences among species, prokaryotes have relatively more compact metabolic pathways, while archaea possess more unique metabolic reactions, such as methanogenesis. Eukaryotes contain complex organelles, such as mitochondria and chloroplasts, resulting in even more complex metabolic networks. Therefore, to improve the accuracy of model construction and facilitate subsequent dataset maintenance and updates, this invention constructs three BKMD datasets for prokaryotes, eukaryotes, and archaea, respectively.

[0045] During the quality control of these datasets, errors were primarily corrected regarding reaction mass balance, charge balance, reaction direction, and net ATP and reducing power generation. Furthermore, duplicate reactions in the BKMD dataset were removed, and metabolite IDs and reaction IDs were standardized, prioritizing BiGG database identifiers, followed by KEGG and then MetaCyc. This standardized labeling makes the metabolic pathways in the model more clearly defined and easier to represent. In the process of removing redundant reactions, reaction IDs were retained according to the priority order of BiGG > KEGG > MetaCyc, and redundant reactions were removed without considering circulating metabolites.

[0046] 2. Metabolite Dataset The BKMD dataset primarily stores reaction data, encompassing 45,734 metabolites. Since building GEMs requires not only reaction data but also information on all metabolites included in the reactions within the model, extracting this information from the reaction dataset for each build would consume significant computational resources and time. Therefore, a dedicated metabolite dataset was created to store metabolite data. Independently storing metabolite information avoids redundant calculations, reduces I / O accesses, and accelerates modeling. It also facilitates subsequent updates and expansions, aids in ID standardization, and reduces data redundancy.

[0047] Based on the cleaned and standardized BKMD dataset of prokaryotes, archaea, and eukaryotes, a Python script was used to extract metabolites from each reaction in the dataset, obtaining information for each metabolite. Then, duplicate metabolites were removed. The final metabolite set was obtained and stored in an Excel worksheet named "metabolites," facilitating the subsequent generation of an Excel spreadsheet used to build the initial model.

[0048] The metabolite dataset contains comprehensive metabolite information, including metabolite ID, name, molecular formula, charge, compartmentalization information, and annotations. BiGG IDs are concise and readable; for example, D-glucose is typically represented as "glc__D_c", where "glc" is an abbreviation for glucose, "__D" represents D-type glucose, "_c" represents cell compartmentalization information, and "c" is an abbreviation for cytosol. BiGG IDs are highly readable, which is the main reason for their priority use in research. In contrast to BiGG, MetaCyc reaction IDs lack standardization. For example, "CPD__45__16942_c" is the ID for 4-oxovalerate, using internal data labels. KEGG has similar issues; for example, rotenone is represented as "C07593" in KEGG IDs. This internal numbering format requires searching the dataset, which is inconvenient for researchers. The processed metabolite data is mainly divided into 6 columns, with the data structure as follows: Figure 2 As shown, the specific content is as follows: (1) The id column is the metabolite ID.

[0049] (2) The name column contains the English names of the metabolites.

[0050] (3) The "compartment" column is an abbreviation of the cell compartment name.

[0051] (4) The values ​​in the charge column represent the amount of charge carried by the metabolite.

[0052] (5) The contents of the formula column are molecular formulas.

[0053] (6) The annotation column contains annotation information. It is currently empty in the metabolite dataset. An annotation information will be added separately in subsequent steps.

[0054] 3. Gene-protein-reaction rule dataset The GPR rule dataset is a crucial step in the construction of this research tool, serving as a bridge between sequence alignment and metabolic response set extraction. The accuracy and quality of the GPR rule data directly affect the genetic basis of the metabolic response set and the accuracy of the constructed model in predicting gene knockout targets. By analyzing the BiGG, KEGG, and MetaCyc models, the GPR rule data contained within the models was extracted. Then, Pandas data processing was used, and a Python script was employed to extract the correspondence between responses and gene IDs, saving the data to the file bigg_MetaCYC_gprs.csv. The data extracted from the BiGG and MetaCyc databases included 267,900 prokaryotic GPR rules, 364,007 eukaryotic GPR rules, and 722,510 archaeal GPR rules. For KEGG, the ko_reaction.list and ko_genes.list files from the provided ko.tar file were used to extract species-indiscriminate GPR rules by ko number, yielding 15,300,235 rules. For the detailed GPR rule dataset construction process, see [link to relevant documentation]. Figure 3 As shown.

[0055] II. GEMs Construction 1. Gene function annotation This invention constructs GEMs from scratch based on protein sequences. First, a local alignment library is built by collecting protein sequences of the target species. Then, the DIAMOND tool is used to score sequence similarity, enabling efficient gene function annotation.

[0056] This invention downloads the protein sequence of a bacterial species from the UniProt database using the species information, then extracts the amino acid sequences of all these genes and constructs FASTA files for prokaryotes, eukaryotes, and archaea, respectively. The dataset contains 31,260 prokaryotic sequences, 58,898 eukaryotic sequences, and 77,823 archaea sequences. These sequences are then converted into dmnd format files, which serve as the sequence alignment database for this tool. The dmnd format is a dedicated alignment database format for DIAMOND, offering higher alignment efficiency compared to traditional FASTA files. The dmnd format is built upon FASTA files, creating an index during construction for rapid target sequence location during DIAMOND alignment. Traditional FASTA formats store sequences in text format; dmnd uses a binary format, reducing file size. The dmnd format can be directly used for alignment with DIAMOND's BLASTP and BLASTX models. FASTA files can be easily converted to dmnd files using built-in DIAMOND functions.

[0057] This tool supports inputting FASTA format amino acid sequences (.faa) or nucleotide sequences (.fna). It uses the DIAMOND alignment tool and the pre-constructed dmnd sequence database for functional annotation. Through DIAMOND built-in functions and Python scripts, it automatically identifies whether the uploaded file is an amino acid or nucleotide sequence, and then automatically performs BLASTP or BLASTX functionality.

[0058] Following the above process, a TSV format file is obtained, containing the alignment results between the user-uploaded FASTA file and the local sequence database. The result file contains 12 columns, with the data structure using the alignment results of E. coli K-12 strain as an example. Figure 4 As shown, the specific content is as follows: (1) Query accession, the sequence identifier in the FASTA file entered by the user.

[0059] (2) Target accession represents the sequence ID in the local database. The first element is the local model name, the second element is the gene ID of the sequence, and the third column is the UniProt ID, separated by ".".

[0060] (3) The third column is Sequence identity, which is the proportion of amino acids that are completely matched in the alignment region.

[0061] (4) Alignment length, the total length of local alignment, including the positions of matches, mismatches and insertions / deletions.

[0062] (5) Mismatches: The number of mismatched amino acids in the comparison region.

[0063] (6) Gap openings: The number of gaps (insertions / deletions) that appear during the alignment.

[0064] (7) Query start, the starting position of the comparison area in the query sequence (Query).

[0065] (8) Query end: The end position of the comparison area in the query sequence.

[0066] (9) Target start, the starting position of the alignment region in the target sequence (Target).

[0067] (10) Target end, the end position of the alignment region in the target sequence.

[0068] (11) E-value, expected value (E-value), represents the probability that the alignment result appears in a random database. The smaller the value, the more significant the alignment (usually, expected value < 1e). -5 (Considered a significant match).

[0069] (12) The comparison score is calculated based on the Smith-Waterman algorithm. The higher the value, the higher the degree of matching between the Query and the Target.

[0070] 2. Dynamic threshold method for filtering annotation results After obtaining the DIAMOND alignment results, the results need to be filtered and sorted to remove low-scoring results and ensure accuracy for subsequent biochemical reaction mapping. Next, Pandas is used to read the data and check the contents of the result file to ensure that the file exists and is not empty; otherwise, an error will be thrown. Then, the data is sorted in descending order based on the alignment scores in the third column, duplicate genes are removed, and higher-scoring matches are retained to ensure high-scoring results in the final dataset.

[0071] To accommodate the alignment distribution of genomic data of varying sizes, this invention employs a dynamic thresholding method, which sets screening criteria based on the distribution of gene alignment scores. This step effectively removes erroneous and redundant responses introduced by sequences with low scores.

[0072] 3. Biochemical reaction mapping based on annotation results The annotation results allow mapping biochemical reactions from the GPR rule dataset to the BKMD dataset. First, in the TSV of the DIAMOND alignment results, the second column contains the gene ID. The GPR rule dataset has four columns: gene ID, protein combination, response ID, and model ID. Gene IDs in the TSV results that meet certain conditions are searched in the GPR rule dataset to obtain the corresponding response ID and protein combination. Then, the response ID is searched in the ID column of the BKMD response dataset; if a match is found, the entire row is extracted. This process is repeated for all gene IDs that meet the conditions, ultimately yielding the response set.

[0073] This reaction set contains comprehensive information, including reaction ID, reaction name, econometric matrix, reaction upper and lower boundaries, GPR rules, annotation information, target reaction, and reaction equation. We store this reaction set in a DataFrame for later retrieval. Using a dynamic thresholding procedure, we extract the rows from the reaction dataset that match the 'id' column, using the second element (reaction ID) of the second column. This result is the reaction set stored in `reactions.csv`. The data structure is shown below. Figure 5 The specific information is as follows: (1) The contents of the id column are the reaction IDs.

[0074] (2) The name column is the English name of the reaction.

[0075] (3) The metabolites column is a stoichiometric matrix generated based on the stoichiometric coefficients in the chemical reaction equation of the reaction, which is used to construct the SBML model and perform FBA calculations.

[0076] (4) The lower_bound and upper_bound columns are the lower and upper boundaries of the reaction, respectively, and are used to constrain the reaction rate.

[0077] (5) The gene_reaction_rule column contains GPR rules. When obtaining the reaction set, there is no GPR rule information. During the initial model building step, GPR rule data will be added through a Python script.

[0078] (6) The annotation column contains annotation information, and there is a dedicated step to add annotations later.

[0079] (7) The objective_coefficient column is the target reaction. Enter the number 1 in the target reaction row and leave the rest of the rows blank.

[0080] (8) The reaction_eq column contains chemical reaction equations.

[0081] 4. Add GPR rules to the biochemical reaction center Based on the GPR rule dataset constructed above, the data structure is as follows: Figure 6 As shown. The specific rules are as follows: The user sequence and the local sequence are aligned using DIAMOND. The score obtained is the alignment score between the local gene ID in the `gene` column and the gene ID of the input sequence. The corresponding `protein` column contains enzyme catalysis information. For each reaction, it may be catalyzed by a single enzyme or by multiple enzymes simultaneously. Based on the result with the highest DIAMOND alignment score, all subunits catalyzing that reaction are found in the GPR rule dataset and replaced. For isoenzymes, "or" is used for connection, while for complexes composed of multiple subunits, "and" is used.

[0082] By combining compartmentalization information and the BKMD GPR rule dataset, this invention further extracts the metabolic reactions catalyzed by the enzyme and constructs an initial metabolic model accordingly, ensuring the model's biological rationality and consistency with cell structure. The specific details are as follows: (1) The content of the gene column is the gene ID of the local dataset.

[0083] (2) The protein column follows the GPR rule.

[0084] (3) The reaction column is the reaction ID in the corresponding reaction dataset.

[0085] (4) The model column is the source of the corresponding model.

[0086] 5. Construct the Biomass equation dataset The biomass models used in this invention are primarily derived from the BiGG database. Specifically, the prokaryotic biomass models mainly reference models such as *Escherichia coli* iML1515, *Bacillus subtilis* iYO844, *Klebsiella pneumoniae* iYL1228, *Corynebacterium glutamicum* iCW773, and *Bacillus putida* iJN1463; the eukaryotic biomass models include *Plasmodium berghei* iAM_Pb448, *Trypanosoma cruzi* iIS312, diatoms iLB1027_lipid, and *Saccharomyces cerevisiae* iMM904.

[0087] During model construction, this invention generates dmnd format alignment libraries for the sequences of all species in Table 3, which are used as local libraries for DIAMOND alignment. Subsequently, based on the FASTA sequence files uploaded by users as input, the input sequences are compared with the local libraries. This allows for the identification of species with the most metabolic similarity to the user-uploaded sequences. The corresponding local Biomass mappings are then added to the response set to ensure the model's rationality and biological relevance.

[0088] Table 3 Biomass reaction

[0089] 6. Generate the initial model The Python script automates the process of integrating the added culture medium, spontaneous reactions, Biomass equations, and reaction sets into a metabolite Excel file, creating a separate worksheet named "reactions". This Excel file now contains all metabolites and the constructed reaction set. The Python script then directly generates an initial model from the Excel file in SBML format.

[0090] 7. Fill in the gaps in the model The initial models constructed above suffer from metabolic gaps during the construction process due to incomplete gene annotation, missing metabolic databases, and incomplete experimental data. This results in the models failing to grow normally or making unreasonable predictions.

[0091] This invention supports customizing the culture medium substrate according to gap-filling requirements, allowing users to choose and add specific amino acids, vitamins, and other nutrients. During gap-filling, this invention differentiates the DIAMOND alignment results into confidence levels 1, 2, and 3, assigning weighting factors to prioritize responses with higher similarity scores during the gap-filling process. This makes GEMs more accurate and interpretable. It also increases the probability of metabolic pathways with stronger genetic evidence being filled into GEMs, ensuring the biological validity of the model.

[0092] 8. Remove redundant reactions The GEMs constructed using the above steps possess all the functions and metrics of high-quality GEMs. However, since the extraction reactions are based on gene annotation, they may introduce similar reactions catalyzed by the same enzyme or similar reactions from the database. To make the GEMs more concise and accurate, redundant reactions need to be removed. First, obtain all metabolic reaction IDs and construct a reaction dictionary. Then, remove the influence of circulating metabolites in the reactions, such as non-specific metabolites like water, inorganic phosphorus, and hydrogen ions. Retain unique reactions. Next, score the correctness of the reactions, for example, using mass balance detection. Calculate the reaction with the highest mass imbalance score and prioritize its deletion. If all these reactions satisfy mass balance, prioritize retaining the reaction with the highest reversibility. It is necessary to check whether deleting a reaction affects the growth of the metabolic network; if not, delete it.

[0093] Finally, an FVA analysis was performed on the model to calculate the minimum and maximum fluxes for each reaction. If both were zero, the reaction could be deleted; otherwise, the reaction with the least impact on flux was removed. After optimization, the model's target growth rate and ATP production rate remained stable, avoiding overfitting.

[0094] 9. Annotation Information and Model Evaluation To further improve the biological accuracy and resolvability of the model, this invention integrates rich and comprehensive annotation information. After generating the model, a script is used to add annotation information to the constructed model. This includes database information such as KEGG, ChEBI, BiGG, MetaNetX, and PubChem; at the gene level, it provides correspondences with databases such as UniProt, KEGG, RefSeq, and NCBI; and in the reaction annotation, it covers EC codes, Rhea, KEGG, and BiGG database information.

[0095] In the model evaluation, the MEMOTE tool (https: / / github.com / opencobra / memote) was used. Since its results are usually output in HTML format, and our research aims to develop a website, this invention uses the constructed model as input and then embeds the generated MEMOTE results page into the page developed by this invention for display.

[0096] III. Platform Implementation Technologies and Development Tools 1. Front-end development framework This website's front-end development framework is based on the FastAPI framework, leveraging the resource advantages of the Taihang-1 high-performance data computing platform to provide users with efficient and scalable data computing and task management services. The website achieves front-end and back-end linkage, task allocation, and resource optimization through distributed computing resources and a modern scheduling system. This invention uses the pFBA method to calculate the flux distribution of the optimal path, converts the reaction list in the path into a metabolite map, and then returns the generated map file to the front end for visualization. D3flux (https: / / github.com / pstjohn / d3flux) is used to automatically generate the layout of the path map and present it intuitively on the webpage. With PEAR, users can easily obtain GEMs models, calculate the optimal synthesis rate and metabolic flux distribution of specific products, and intuitively visualize and analyze synthetic pathways, providing important guidance and support for metabolic engineering and biomanufacturing.

[0097] 2. Programming tools and software packages The entire process of the program and data processing in this invention is based on Python and its third-party libraries. For details and version information, please refer to Table 4.

[0098] Table 4 Content and Version Information

[0099] IV. Results and Discussion 1. Development of automated build processes for GEMs This invention develops PEAR (Platform for Efficient Automated Reconstruction), an online automated GEMs construction tool that integrates functions such as model building, computational simulation, model evaluation, and pathway visualization.

[0100] First, metabolic reaction and metabolite data were collected from the BiGG, KEGG, and MetaCyc databases. Metabolite IDs were standardized, and duplicate reaction data were removed. Metabolite names from the three databases were standardized, prioritizing the use of BiGG IDs as unique identifiers. If a metabolite was not found in the BiGG database, then the KEGG or MetaCyc databases were used to represent it, ensuring each metabolite has a unique identifier and maintaining the connectivity of the metabolic network. Then, errors in energy, reducing power, and net metabolite production cycles were corrected, and three BKMD datasets (prokaryotes, archaea, and eukaryotes) were constructed to store biochemical reaction data.

[0101] Next, based on the GPR rule dataset constructed in the above steps, the corresponding amino acid sequences are obtained from the UniProt database based on gene IDs, and a dmnd sequence alignment database is constructed. After inputting FASTA sequence files, DIAMOND is used for alignment, and the similarity score results are filtered. Enzymatic reactions are mapped to the BKMD dataset, and spontaneous reactions and culture medium reactions are added to construct a biochemical reaction set for the target organism. Subsequently, through sequence similarity comparison, the bacterial species most similar to the user's bacterial metabolic network is found. Based on the Biomass equation and bacterial species information from the above steps, the most suitable Biomass equation is assigned to the model, completing the initial model construction.

[0102] However, the initial model may contain gaps, mainly due to insufficient completeness of the reaction database and gene annotations. These gaps typically require gap-filling algorithms. During the gap-filling stage, PEAR uses the pFBA method and a centralized strategy to prioritize the introduction of enzymatic reactions with high similarity scores. This ensures that all reactions added during gap-filling are genetically based, making GEMs more accurate and interpretable.

[0103] For eukaryotic model construction, PEAR adds a step to obtain cellular sublocalization information before generating the initial model. PEAR parses the input FASTA sequence, obtains cellular sublocalization information of genes from the UniProt database via the Python API, adds compartmentalization information for specific reactions in the model, and then performs the same gap-filling process as for prokaryotes. Finally, it adds annotation information for metabolites, reactions, and genes to the model and generates the final model. The construction process will be described in detail below.

[0104] 2. Construction of Biochemical Reaction Dataset 1) Data cleaning and standardization In the preceding steps of constructing the biochemical dataset, the BKMD dataset was obtained. Since the biochemical reaction dataset integrates biochemical reactions from three databases—BiGG, KEGG, and MetaCyc—the metabolite IDs differ between these databases. This leads to duplicates where the same reaction is represented by different IDs, or redundancy due to different circulating metabolites in the reaction. Therefore, the BKMD dataset underwent data cleaning and standardization.

[0105] To eliminate redundant reactions, the IDs of each metabolite first need to be standardized and uniformly represented. This effectively distinguishes duplicate reactions, making the metabolic pathways in the model more hierarchical and easier to express. For example, the MetaCyc and KEGG reaction IDs "RXN-19380" and "R12505_2" respectively, after standardized identification of metabolites, clearly show that they represent the same reaction: "aspsa_c + 2.0 fdxrd_c + h2s_c + 2.0 h_c --> 2.0fdxox_c + h2o_c + hcys__L_c". Among the metabolite IDs in the BiGG, KEGG, and MetaCyc databases, BiGG IDs are highly readable, while MetaCyc IDs can also display some reaction information. KEGG IDs use internal labels, making them less readable. Therefore, a Python script is used to uniformly represent the metabolite IDs in the reactions, retaining the reaction IDs according to the priority order of BiGG > KEGG > MetaCyc. Then, based on the database annotation information of metabolites in the model, a Python script was used to process the equations, searching for identical and similar reactions. Only one identical reaction was retained, with the reaction ID retained in the same order as the metabolite. Rules were established to remove redundancy from similar reactions: reactions with identical equations after removing circulating metabolites were considered similar. The script checked the mass and charge balance of the reactions, prioritizing the removal of unbalanced reactions. Finally, data cleaning and standardization of the BKMD biochemical reaction datasets for the three species were completed. The prokaryotic, archaea, and eukaryotic BKMD biochemical reaction datasets were stored in BKMD_Prokaryotes.csv, BKMD_Archaea.csv, and BKMD_Eukaryotes.csv, respectively, for subsequent biochemical reaction mapping.

[0106] 2) Dataset error correction After data cleaning, the three processed biochemical reaction datasets were converted into SBML models, and error checking was performed on the three models. When using D-glucose as the substrate as input, calculations revealed that the actual yield of each carbon-containing metabolite in the model exceeded the theoretical yield, which is inconsistent with biological reality and indicates that the current model contains errors and cannot accurately quantify product yields. The main reasons for these errors are issues such as non-conservation of matter, unlimited generation of reducing power, and unlimited energy generation. Examples of corrections are shown in Table 5. For example, in SUCptspp_1: pep_c => pyr_c + suc6p_c, an error occurs in the conservation of matter, leading to the unlimited generation of D-gluconic acid. Other errors due to non-conservation of matter, such as the unlimited generation of reducing power (nadh) and net ATP energy generation, also prevent the model from accurately quantifying these issues. Furthermore, the mass imbalance reaction RE2477: amet_c + CE3087_c => ahcys_c + h_c + C09209_c causes the error of unlimited nadh generation. The reaction atp_c + coa_c + ddca_c <=> amp_c + ddcacoa_c + ppi_c is conserved in mass, but it will cause the error of infinite ATP production. The actual direction of the reaction should be irreversible.

[0107] Table 5 Examples of reaction correction results

[0108] To address the aforementioned errors, a Python script was used for automated detection and correction. First, exchange and transport reactions were excluded. Reactions with mass or charge imbalances were penalized 3 points, and reactions lacking annotation information were penalized 1 point. For reactions with mass balance, the reaction direction was corrected. If a reaction was mass-imbalanced but had annotation information, it was rebalanced; reactions without annotation information were deleted.

[0109] 3) Respiratory chain reaction treatment The respiratory chain reaction generates energy by forming a proton gradient, which could lead to unlimited energy generation. Therefore, 21 respiratory chain reactions were deleted, and a merged respiratory chain and proton exchange reaction ADD_h_c: h_c <=> was added to make the proton gradient of ATP generation irrelevant and avoid errors that could lead to unlimited energy generation. The specific reactions that were processed are shown in Table 6.

[0110] Ultimately, after data cleaning, standardization, and error correction, the prokaryotic BKMD dataset contained 62,488 reactions; the eukaryotic BKMD dataset contained 70,736 reactions; and the archaea BKMD dataset contained 79,647 reactions. These three BKMD datasets will be used for biochemical reaction mapping.

[0111] Table 6 Respiratory chain reaction treatment

[0112] 3. Design of the genome functional annotation and result screening workflow PEAR employs a top-down strategy, constructing GEMs from scratch. First, it uses DIAMOND for similarity comparison, scoring the similarity between the input sequence and the alignment library, and then annotating gene functions based on the alignment scores. In reality, due to the significant differences in genome size and sensitivity to alignment libraries among species, this invention develops a dynamic threshold strategy that automatically calculates an appropriate alignment score threshold for each user-uploaded sequence.

[0113] First, the number of duplicate genes is counted and stored in the variable `gene_num`. 75% of the total number of genes is taken as the threshold gene count and stored in the variable `threshold`. Predefined similarity score intervals are used: [100, 80], [79.99, 60], [59.99, 40], [39.99, 35], [34.99, 30], [29.99, 25], and [24.99, 20]. The total number of genes falling into each interval is counted and summed to determine the range within which the threshold falls. The lower bound of this interval is used as the final screening threshold. If no gene falls into any interval, the default threshold of 15 is used to ensure proper model construction. Sequences scoring below 15 are automatically considered to have low genetic evidence and are discarded. This threshold serves as the minimum score for screening annotation results, addressing the limitations and errors of fixed parameter settings and effectively resolving errors such as under-prediction and over-prediction.

[0114] 4. Mapping biochemical reactions based on annotation results 1) Mapping enzyme reaction Based on the above-selected annotation results, PEAR performs biochemical reaction mapping using the BKMD dataset to obtain enzymatic and transport reactions. Then, based on the gene annotations for each enzymatic reaction, it maps them using a GPR rule dataset to assign GPR rules to the biochemical reactions. This ultimately generates a set of enzymatic reactions.

[0115] 2) Add spontaneous reaction The enzyme-catalyzed reaction set was obtained through gene annotation, while spontaneous reactions typically occur without enzyme catalysis or any external factors. Therefore, these reactions cannot be added to the reaction set via gene annotation. These reactions are crucial; their absence will prevent the model from growing properly. Spontaneous reactions do not require enzyme catalysis and are primarily determined by chemical kinetics and environmental conditions, including processes such as redox reactions, chemical degradation, or metal ion coordination. Referring to the iML1515 model, the necessary spontaneous reactions listed in Table 7 were added to the model to ensure proper growth simulation. Finally, the model was integrated with the enzyme-catalyzed reaction set.

[0116] Table 7 List of Spontaneous Reactions

[0117] 3) Add cellular sublocalization information of the reaction This step is only used for constructing eukaryotic metabolic network models. Constructing eukaryotic metabolic network models is typically quite complex, as eukaryotes contain more intricate metabolic systems, including mitochondria, the nucleus, the Golgi apparatus, and many other organelles. Each organelle not only has its own specific reactions, but there is also material exchange between the organelle and the cytoplasm; these reactions are indispensable. Although most eukaryotic-specific reactions can be obtained through gene annotation and sequence alignment, it is still not comprehensive enough. To address this issue, this invention uses a Python script via the UniProt API to obtain cell sub-localization information and add it to the model. First, we parse the user-uploaded sequence file, iterating through each sequence in the Python script to obtain cell sub-localization information one by one. The obtained sub-cell localization information is saved as a DataFrame. Then, we parse and identify the cell sub-localization information string; specific compartmentalization information is shown in Table 8. Finally, the corresponding compartmentalization information is added to the reaction set.

[0118] Table 8 Cell compartment information

[0119] 4) Add culture medium for reaction The aforementioned culture media provide the necessary nutrient environment for microbial or cell growth, offering essential carbon sources, nitrogen sources, inorganic salts, vitamins, and other growth factors to maintain normal cell growth. The composition of the culture medium directly affects the organism's growth rate, metabolic flux, product synthesis, and metabolic regulation. Culture media are typically categorized as complex media (LB) and chemically synthesized media (M9). Different media can be used to optimize growth conditions, regulate metabolic flux, and simulate physiological environments, such as increasing cell density, enhancing the synthesis of target products, and simulating intestinal and soil environments. When studying auxotrophic strains, which lack the ability to synthesize certain substances, thus hindering their normal growth, specific culture media must be used to artificially supplement these substances, such as certain amino acids or vitamins.

[0120] Culture media are equally indispensable in metabolic network models, and are typically set up to simulate actual growth conditions, such as LB and M9 media. Currently, almost all published metabolic network model building tools do not allow users to customize the culture medium composition; they can only simulate the culture medium through built-in functions or pre-set parameters. This reduces the applicability and accuracy of GEM metabolic simulations. Therefore, this tool addresses these issues by adding a user-defined culture medium function. It also supports adding essential substances that auxotrophic bacteria cannot synthesize, such as amino acids or vitamins, thus broadening the tool's applicability.

[0121] In addition, referencing the exchange reactions in the *E. coli* model iML1515, the ion exchange reactions listed in Table 9 were added to each culture medium. The upper and lower boundaries of these reactions were set to 1000 and -1000, respectively. Other exchange reactions were only allowed to proceed in the forward direction, with upper and lower boundaries of 1000 and 0, respectively. The carbon sources and vitamins that can be added to the culture medium are shown in Table 10, and all 20 amino acids can be added to the model.

[0122] Table 9 Reversible ion exchange reactions

[0123] Table 10 Carbon Sources and Vitamins

[0124] 5) Add Biomass equations In the Biomass equation addition step, DIAMOND was used to compare the user-input sequence with the sequence of the species corresponding to the Biomass equation. The highest alignment score indicates greater metabolic similarity, thus the Biomass equation for that species is introduced. Compared to current construction tools, PEAR has 10 Biomass equations, covering species such as bacteria, archaea, protozoa, algae, and fungi, making it more accurate. For example, CarveMe indicates it has 2 Biomass equations. See [link to Biomass equation documentation] for the detailed Biomass process. Figure 7 As shown.

[0125] 5. Initial Model Generation The constructed biochemical reaction dataset is saved in the "reactions" worksheet of an Excel spreadsheet. Next, the gene information of the reaction alignment is output to the "gene" worksheet of an Excel spreadsheet, and the constructed metabolite information is output to the "metabolites" worksheet of an Excel spreadsheet. All of the above content is saved in an Excel spreadsheet. Finally, an SBML format model is generated using a Python program.

[0126] The script first reads the Excel file, identifies the information in the `metabolites` and `reactions` worksheets, parses the data using pandas to remove null values ​​and correct formatting, converts the data into metabolite objects recognizable by COBRApy, and adds the metabolites to the model. Then, it parses the matrix information, annotation information, lower_bound and upper_bound constraints, GPR rules, and other information in the `OrderedDict` column of the reaction data, converts it into an SBML language model, and outputs it in XML format.

[0127] 6. Development of an automated model refinement process 1) Develop weighting algorithms for model gap filling The initial model may contain metabolic gaps, requiring gap-filling. To minimize the introduction of metabolic pathways and avoid redundancy, while ensuring that the introduced responses have strong genetic evidence, this study categorizes responses into three confidence levels based on the selected DIAMOND scores: the top 50% are classified as confidence level 1, 50% to 80% as confidence level 2, and the bottom 20% as confidence level 3. Subsequently, the coefficients of the responses corresponding to confidence levels 1, 2, and 3 are multiplied by weights of 1000, 600, and 300, respectively, and added to the gap-filling comparison library. This prioritizes the inclusion of responses with stronger genetic evidence during gap-filling, resulting in a more accurate and reliable model.

[0128] 2) Redundant reaction removal process To improve the accuracy of the constructed GEMs, this invention incorporates a redundant reaction removal step for automated model refinement. The process begins by loading the model and calculating the initial objective function value. The program then iterates through all reactions, removing transient metabolites and identifying redundant reaction groups with identical substrates and products. Nested subsets within these redundant groups are then removed, and the optimal reaction in each group is retained based on indicators such as mass balance, reversibility, and metabolite quantity. Reactions involving ATP and ADP, and empty reactions (containing only water or similar substances), are categorized and treated specially. Next, all unretained redundant reactions are attempted to be deleted without affecting the model's objective function value, and their gene rules are merged into the retained reactions. The system further identifies and removes reversible reactions involving only ATP and ADP, as well as completely repetitive reactions. Finally, based on FVA analysis results, unnecessary redundant reactions associated with single genes are iteratively removed, and metabolites not involved in the reaction are cleaned from the model, resulting in a streamlined metabolic model with complete information retention.

[0129] 3) Add annotation information and generate the final model Annotation information is added to the refined model, and the cross-referenced database annotation information for reactions, metabolites, and genes in the model is improved. Database names and URLs are shown in Table 11. Finally, the model is saved in SBML, JSON, MAT, or EXCEL formats for users to download.

[0130] Table 11 Database Names and URLs

[0131] This invention is not limited to the structures described above and shown in the accompanying drawings, and various modifications and changes can be made without departing from its scope. The scope of this invention is limited only by the appended claims.

Claims

1. A method for constructing a genome-scale metabolic network model, characterized in that, include: Biochemical reaction data from biological metabolism databases are integrated and preprocessed. Based on the preprocessed biochemical reaction data, biochemical reaction datasets, metabolite datasets, and gene-protein-reaction rule datasets are constructed. A local alignment library is constructed based on the protein sequences of the target species. The similarity between the input sequence and the sequences in the local alignment library is compared, and gene function annotation is performed based on the alignment results. Based on the dynamic threshold method, the alignment score threshold of the input sequence is determined, and the gene function annotation is screened in combination with the alignment score threshold. Based on the selected gene function annotation results, biochemical reaction mapping is performed by combining the biochemical reaction dataset and the gene-protein-reaction rule dataset. Based on the mapping results, a biochemical reaction set is constructed, and relevant metabolite dataset information is added to the biochemical reaction set. Species with the most similar metabolism to the input sequence are selected, the local biomass equations corresponding to the similar species are mapped and added to the biochemical reaction set, and the added biochemical reaction set is integrated into the metabolite file. An initial genome-scale metabolic network model is generated based on the metabolite file. The initial genome-scale metabolic network model was sequentially processed by gap filling, redundant response removal, and annotation information addition to obtain the final genome-scale metabolic network model. Specifically, the process of mapping biochemical reactions based on the selected gene function annotation results, combined with biochemical reaction datasets and gene-protein-reaction rule datasets, constructing a biochemical reaction set based on the mapping results, and adding relevant metabolite dataset information to the biochemical reaction set includes: Gene IDs are extracted from the sequence alignment results, and gene IDs whose scores meet the threshold are searched in the gene-protein-reaction rule dataset to obtain the corresponding reaction IDs and protein combinations. Based on the reaction ID, match rows with the same ID column in the biochemical reaction dataset and extract the entire row of data; repeat the search and extraction operation for all gene IDs that meet the threshold, and summarize to obtain the biochemical reaction set, which includes reaction ID, reaction name, metric matrix, upper and lower boundaries of reaction, gene-protein-reaction rules, annotation information, target reaction and reaction equation; Based on the result with the highest alignment score, find all subunits of catalytic reactions in the gene-protein-reaction rule dataset, search for the compartment where the subunit is located in the Uniprot database and replace it with the correct compartment to complete the compartment mapping of the effective subunit; By employing or linking isozymes and connecting them to complexes composed of multiple subunits, the Boolean connection logic can be reconstructed. By combining compartment information, biochemical reaction datasets, and gene-protein-reaction rule datasets, we extract the catalytic metabolic reactions of gene-annotated enzymes and complete the adaptation of rules to reactions. Add the effective subunit and complex related dataset information after compartment mapping to the biochemical reaction set.

2. The method for constructing a genome-scale metabolic network model according to claim 1, characterized in that, The integrated biochemical reaction data from the biological metabolism database is preprocessed, and based on the preprocessed biochemical reaction data, a biochemical reaction dataset, a metabolite dataset, and a gene-protein-reaction rule dataset are constructed, including: Biochemical reaction data from biological metabolism databases are integrated, and the biochemical reaction data are sequentially cleaned and standardized, error corrected, and processed for respiratory chain reactions. Based on the processed biochemical reaction data, biochemical reaction datasets for prokaryotes, eukaryotes, and archaea are constructed. Extract metabolite information for each biochemical reaction from the biochemical reaction dataset and remove duplicates. Construct a metabolite dataset based on the deduplicated metabolite information. Extract gene-protein-reaction rule data from the biochemical reaction dataset and obtain the correspondence between reactions and gene IDs. Based on the correspondence between reactions and gene IDs, construct a gene-protein-reaction rule dataset.

3. The method for constructing a genome-scale metabolic network model according to claim 2, characterized in that, The metabolite information includes metabolite ID, name, molecular formula, charge, compartment information, and annotation information.

4. The method for constructing a genome-scale metabolic network model according to claim 1, characterized in that, The process involves constructing a local alignment library based on the protein sequences of the target species, performing similarity comparisons between the input sequence and the sequences in the local alignment library, and annotating gene functions based on the alignment results. Based on the dynamic thresholding method, the alignment score threshold of the input sequence is determined, and the gene functional annotation is screened in combination with the alignment score threshold, including: The protein sequences of the target species were obtained, and text-based bioinformatics standard format files for prokaryotes, eukaryotes and archaea were constructed respectively. The text-based bioinformatics standard format files were then converted into the special format files of sequence alignment software to build a local alignment library. The system receives input sequences uploaded by users, automatically identifies the type of input sequences using built-in functions and scripts of the sequence alignment software, calls the amino acid sequence or nucleotide sequence model in the sequence alignment software to perform sequence similarity comparison with the local alignment library, outputs the sequence alignment result file, and completes gene function annotation. The sequence alignment result file is checked and ensured to exist and not be empty. The annotation results are sorted in descending order according to the alignment score of the sequence alignment result file, and duplicate gene matches are removed. Based on the sorting distribution of alignment scores, a dynamic threshold method is used to determine the alignment score threshold of the input sequence, and gene function annotation results with alignment scores lower than the alignment score threshold are filtered out.

5. The method for constructing a genome-scale metabolic network model according to claim 4, characterized in that, The ranking distribution based on alignment scores, and the determination of the alignment score threshold for the input sequence using a dynamic threshold method, include: The total number of genes after deduplication is counted, a preset percentage of the total number of genes is selected as the threshold number of genes, and the threshold number of genes is stored in the critical value variable; The total number of genes falling into the preset similarity score interval is counted and the sum is calculated. Based on the calculation results, the range of intervals into which the critical value falls is determined. Determine whether the range into which the critical value falls belongs to the preset similarity score range. If so, use the lower bound of the range as the comparison score threshold; otherwise, select the preset threshold as the comparison score threshold.

6. The method for constructing a genome-scale metabolic network model according to claim 1, characterized in that, Adding relevant metabolite dataset information to the biochemical reaction set also includes adding spontaneous reactions, cellular sublocalization information of the reactions, and culture medium reactions to the biochemical reaction set.

7. The method for constructing a genome-scale metabolic network model according to claim 1, characterized in that, Gap-filling in the initial genome-scale metabolic network model includes: The sequence alignment scores are divided into several confidence levels, and a weighting factor is assigned to the response corresponding to each confidence level. The gaps in the initial genome-scale metabolic network model are filled based on the weighting factors.

8. The method for constructing a genome-scale metabolic network model according to claim 1, characterized in that, Redundancy removal from the initial genome-scale metabolic network model includes: Traverse all reactions in the initial genome-scale metabolic network model, remove circulating metabolites in the reactions, and identify redundant reaction sets with the same substrates and products. Remove nested subsets from redundant reaction groups and retain the optimal reaction in each group based on mass balance, reversibility, and metabolite quantity indices; Eliminate reactions containing adenosine triphosphate and adenosine diphosphate, as well as empty reactions, from redundant reaction groups; Remove reactions that were not retained from the redundant reaction group and merge the gene rules of the unretained reactions into the retained reactions; Identify and remove reversible reactions that contain only adenosine triphosphate and adenosine diphosphate, as well as completely repeatable reaction formulas; Based on the results of flux variability analysis, unnecessary redundant reactions associated with individual genes are cyclically removed, and metabolites that do not participate in the reaction are cleaned up, so as to achieve redundant reaction removal processing of the initial genome-scale metabolic network model.

9. The method for constructing a genome-scale metabolic network model according to claim 1, characterized in that, Adding annotation information to the initial genome-scale metabolic network model includes: At the gene level, add the correspondence between gene and protein sequence and function databases; At the reaction level, information related to metabolic reactions and enzyme catalysis is added.

Images