Genome scale metabolic network model construction method

Abstract

The invention belongs to the technical field of system biology, and discloses a genome scale metabolic network model construction method, which comprises the following steps: constructing a biochemical reaction data set, a metabolite data set and a gene-protein-reaction rule data set; performing gene function annotation and screening; performing biochemical reaction mapping, constructing a biochemical reaction set, and adding related metabolite data set information; mapping and adding local biomass equations corresponding to similar species into a biochemical reaction set to generate an initial genome scale metabolic network model; and sequentially carrying out vacancy filling, redundant reaction removal and annotation information addition processing on the initial genome scale metabolic network model to obtain a final genome scale metabolic network model. According to the method, the defects of an existing GEMs construction tool in the aspects of model coverage and prediction precision are overcome, the automatic de novo construction process from the genome sequence to the GEMs is developed, the time-consuming manual correction process is reduced, and the construction efficiency is improved.

Description

A method for constructing a genome-scale metabolic network model 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, the annotations in GEMs 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 a genome-scale metabolic network model includes: integrating biochemical reaction data from a biological metabolic database and preprocessing the biochemical reaction data; constructing a biochemical reaction dataset, a metabolite dataset, and a gene-protein-reaction rule dataset based on the preprocessed biochemical reaction data; 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 comparison results; determining the alignment score threshold of the input sequence based on a dynamic threshold method; and filtering gene function annotations based on the alignment score threshold; and performing gene function annotation based on the filtered gene function annotation results. The results were analyzed by combining the biochemical reaction dataset and the gene-protein-reaction rule dataset to perform biochemical reaction mapping. Based on the mapping results, a biochemical reaction set was constructed, and relevant metabolite dataset information was added to the biochemical reaction set. Species with the most similar metabolism to the input sequence were selected, and the local biomass equations corresponding to these similar species were mapped and added to the biochemical reaction set. The added biochemical reaction set was then integrated into the metabolite file, and an initial genome-scale metabolic network model was generated based on the metabolite file. The initial genome-scale metabolic network model was then subjected to gap filling, redundant reaction removal, and annotation information addition processing to obtain the final genome-scale metabolic network model.

[0010] In one embodiment, the process of integrating biochemical reaction data from a biological metabolism database, preprocessing the biochemical reaction data, and constructing a biochemical reaction dataset, a metabolite dataset, and a gene-protein-reaction rule dataset based on the preprocessed biochemical reaction data includes: integrating biochemical reaction data from a biological metabolism database, and sequentially performing data cleaning and standardization, error correction, and respiratory chain reaction processing on the biochemical reaction data; constructing biochemical reaction datasets for prokaryotes, eukaryotes, and archaea based on the processed biochemical reaction data; extracting metabolite information for each biochemical reaction in the biochemical reaction dataset and performing deduplication; constructing a metabolite dataset based on the deduplicated metabolite information; and extracting gene-protein-reaction rule data from the biochemical reaction dataset and obtaining the correspondence between reactions and gene IDs; constructing a gene-protein-reaction rule dataset based on the correspondence between reactions and gene IDs.

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

[0012] In one embodiment, the step 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 based on the alignment score threshold includes: obtaining the protein sequences of the target species, constructing text-based bioinformatics standard format files for prokaryotes, eukaryotes, and archaea respectively, converting the text-based bioinformatics standard format files into a dedicated format file for the sequence alignment software, and building a local alignment library; receiving user-uploaded... The input sequence is automatically identified using built-in functions and scripts of the sequence alignment software. The software's amino acid or nucleotide sequence model is then used to perform sequence similarity comparisons with a local alignment library, outputting a sequence alignment result file to complete gene function annotation. The sequence alignment result file is checked to ensure it exists and is not empty. Annotation results are sorted in descending order based on alignment scores, and duplicate gene matches are removed. Based on the alignment score distribution, a dynamic threshold method is used to determine the alignment score threshold for the input sequence, and gene function annotation results with alignment scores below the threshold are filtered out.

[0013] In one embodiment, determining the alignment score threshold of the input sequence using a dynamic threshold method based on the sorting distribution of alignment scores includes: counting the total number of genes after deduplication, selecting a preset percentage of the total number of genes as the threshold gene number, and storing the threshold gene number in a critical value variable; counting the total number of genes falling into a preset similarity score interval and calculating the sum, determining the interval range into which the critical value falls based on the calculation result; determining whether the interval range into which the critical value falls belongs to the preset similarity score interval, if so, using the lower bound of the interval range as the alignment score threshold, if not, selecting the preset threshold as the alignment score threshold.

[0014] In one embodiment, the step of mapping biochemical reactions 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: extracting gene IDs from the sequence alignment results, and selecting gene IDs whose scores meet the threshold to search in the gene-protein-reaction rule dataset to obtain the corresponding reaction IDs and protein combinations; matching rows with the same ID column in the biochemical reaction dataset according to the reaction IDs, and extracting the entire row of data; repeating the search and extraction operation for all gene IDs that meet the threshold, and summarizing to obtain the biochemical reaction set, wherein the biochemical reaction set includes reaction IDs and reaction names. The dataset includes: a metric matrix, reaction upper and lower boundaries, gene-protein-reaction rules, annotation information, target reactions, and reaction equations. Based on the highest alignment score, all subunits of catalytic reactions are located in the gene-protein-reaction rule dataset. The subunit compartments are searched in the Uniprot database, and the correct compartments are replaced to complete the effective subunit compartment mapping. Isoenzymes are linked using OR, and complexes composed of multiple subunits are linked to achieve Boolean connection logic reconstruction. Combining compartment information, biochemical reaction datasets, and gene-protein-reaction rule datasets, the catalytic metabolic reactions of gene-annotated enzymes are extracted to adapt the rules to the reactions. The effective subunits and complex-related dataset information after compartment mapping are added 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 the initial genome-scale metabolic network model includes: dividing the initial model into several confidence levels based on sequence alignment scores, assigning a weighting factor to the response corresponding to each confidence level, and filling gaps in the initial genome-scale metabolic network model based on the weighting factor.

[0017] In one embodiment, redundant reaction removal processing of the initial genome-scale metabolic network model includes: traversing all reactions in the initial genome-scale metabolic network model, removing circulating metabolites in the reactions to identify redundant reaction groups with the same substrate and product; removing nested subsets in the redundant reaction groups and retaining the optimal reaction in each group based on mass balance, reversibility, and metabolite quantity indices; removing reactions containing only adenosine triphosphate (ATP) and adenosine diphosphate (ADP) and empty reactions from the redundant reaction groups; removing reactions not retained from the redundant reaction groups and merging the gene rules of the not retained reactions into the retained reactions; identifying and deleting reversible reactions containing only ATP and ADP and completely duplicated reaction formulas; and cyclically deleting unnecessary redundant reactions associated with single genes based on flux variability analysis results and cleaning up metabolites that do not participate in the reactions, thereby achieving 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, adding the correspondence between genes and protein sequences and functional databases; and at the reaction level, adding metabolic reaction and enzyme catalysis association information.

[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: a dataset construction module for integrating biochemical reaction data from a biological metabolic database, preprocessing the biochemical reaction data, and constructing a biochemical reaction dataset, a metabolite dataset, and a gene-protein-reaction rule dataset based on the preprocessed biochemical reaction data; a gene function annotation module for 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 performing gene function annotation based on the comparison results; determining the alignment score threshold of the input sequence based on a dynamic threshold method, and filtering gene function annotations based on the alignment score threshold; and a biochemical reaction mapping module for mapping the gene function annotations based on the filtered gene function annotations. Based on the functional annotation results, biochemical reaction mapping is performed using 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 select the species most similar to the input sequence in metabolism, 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. An initial genome-scale metabolic network model is generated based on the metabolite file. The model optimization module is used to sequentially perform 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.

[0021] The technical solutions provided by the embodiments of this invention can 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 screen 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, and develops an automated de novo construction process from genome sequences to GEMs, 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 is a flowchart illustrating the construction of a genome-scale metabolic network model according to an exemplary embodiment; Figure 2 is a diagram illustrating the metabolite data structure according to an exemplary embodiment; Figure 3 is a schematic diagram illustrating GPR rule extraction according to an exemplary embodiment; Figure 4 is a diagram illustrating the alignment result data structure according to an exemplary embodiment; Figure 5 is a diagram illustrating the reaction set data structure according to an exemplary embodiment; Figure 6 is a diagram illustrating the GPR rule data structure according to an exemplary embodiment; Figure 7 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 illustrates an embodiment of a genome-scale metabolic network model construction method according to the present invention.

[0029] In this optional embodiment, the method for constructing the genome-scale metabolic network model includes: Step S101, integrating biochemical reaction data from a biological metabolic database, preprocessing the biochemical reaction data, and constructing a biochemical reaction dataset, a metabolite dataset, and a gene-protein-reaction rule dataset based on the preprocessed biochemical reaction data; Step S102, constructing a local alignment library based on the protein sequences of the target species, performing similarity comparison between the input sequence and the 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 the dynamic threshold method, and filtering gene function annotations based on the alignment score threshold; Step S103, based on the filtered gene functions... The annotation results are combined with the biochemical reaction dataset and the gene-protein-reaction rule dataset to perform biochemical reaction mapping. Based on the mapping results, a biochemical reaction set is constructed, and relevant metabolite dataset information is added to the biochemical reaction set. Step S104: Select the species that are most similar to the input sequence in metabolism, 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. Based on the metabolite file, generate an initial genome-scale metabolic network model. Step S105: Perform gap filling, redundant reaction removal and annotation information addition processing 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, 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: integrating biochemical reaction data from the biological metabolism database, and sequentially performing data cleaning and standardization, error correction, and respiratory chain reaction processing on the biochemical reaction data; constructing biochemical reaction datasets for prokaryotes, eukaryotes, and archaea based on the processed biochemical reaction data; extracting metabolite information for each biochemical reaction in the biochemical reaction dataset and performing deduplication; constructing a metabolite dataset based on the deduplicated metabolite information; and extracting gene-protein-reaction rule data from the biochemical reaction dataset and obtaining the correspondence between reactions and gene IDs; constructing a gene-protein-reaction rule dataset based on the correspondence between reactions and gene IDs.

[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 step 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 the dynamic threshold method, and filtering gene function annotations in conjunction with the alignment score threshold includes: obtaining the protein sequences of the target species, constructing text-based bioinformatics standard format files for prokaryotes, eukaryotes, and archaea respectively, converting the text-based bioinformatics standard format files into a dedicated format file for the sequence alignment software, and building a local alignment library; receiving user-uploaded... The input sequence is automatically identified using built-in functions and scripts of the sequence alignment software. The software's amino acid or nucleotide sequence model is then used to perform sequence similarity comparisons with a local alignment library, outputting a sequence alignment result file to complete gene function annotation. The sequence alignment result file is checked to ensure it exists and is not empty. Annotation results are sorted in descending order based on alignment scores, and duplicate gene matches are removed. Based on the alignment score distribution, a dynamic threshold method is used to determine the alignment score threshold for the input sequence, and gene function annotation results with alignment scores below the threshold are filtered out.

[0033] In this optional embodiment, the step of determining the alignment score threshold of the input sequence using a dynamic threshold method based on the sorting distribution of the alignment scores includes: counting the total number of genes after deduplication, selecting a preset percentage of the total number of genes as the threshold gene number, and storing the threshold gene number in a critical value variable; counting the total number of genes falling into a preset similarity score interval and calculating the sum, and determining the interval range into which the critical value falls based on the calculation result; determining whether the interval range into which the critical value falls belongs to the preset similarity score interval, and if so, using the lower bound of the interval range as the alignment score threshold, and if not, selecting the preset threshold as the alignment 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: extracting gene IDs from the sequence alignment results, and selecting gene IDs whose scores meet a threshold (i.e., the alignment result threshold when the total number of genes is greater than 75% of the total number of annotated genes) to search in the gene-protein-reaction rule dataset to obtain the corresponding reaction IDs and protein combinations; matching rows with the same ID column in the biochemical reaction dataset according to the reaction IDs, and extracting the entire row of data; repeating the search and extraction operation for all gene IDs that meet the threshold, and summarizing to obtain the biochemical reaction set, wherein the biochemical reaction... The dataset should include reaction ID, reaction name, metric matrix, upper and lower boundaries of the reaction, gene-protein-reaction rules, annotation information, target reaction, and reaction equation. Based on the result with the highest alignment score, all subunits of catalytic reactions are searched in the gene-protein-reaction rule dataset. The compartments of the subunits are searched in the Uniprot database, and the correct compartments are replaced to complete the compartment mapping of effective subunits. Isoenzymes are linked using OR, and complexes composed of multiple subunits are linked to achieve the reconstruction of Boolean connection logic. Combining compartment information, biochemical reaction datasets, and gene-protein-reaction rule datasets, the catalytic metabolic reactions of gene-annotated enzymes are extracted to complete the adaptation of rules and reactions. The relevant dataset information of effective subunits and complexes after compartment mapping is added to the biochemical reaction dataset.

[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 in the initial genome-scale metabolic network model includes: dividing the sequence alignment score into several confidence levels, assigning a weighting factor to the response corresponding to each confidence level, and filling gaps in the initial genome-scale metabolic network model according to the weighting factor.

[0037] In this optional embodiment, redundant reaction removal processing of the initial genome-scale metabolic network model includes: traversing all reactions in the initial genome-scale metabolic network model, removing circulating metabolites in the reactions to identify redundant reaction groups with the same substrate and product; removing nested subsets in the redundant reaction groups and retaining the optimal reaction in each group based on mass balance, reversibility, and metabolite quantity indices; eliminating reactions containing only adenosine triphosphate (ATP) and adenosine diphosphate (ADP) and empty reactions in the redundant reaction groups; eliminating reactions not retained in the redundant reaction groups and merging the gene rules of the not retained reactions into the retained reactions; identifying and deleting reversible reactions containing only ATP and ADP and completely duplicated reaction formulas; and cyclically deleting unnecessary redundant reactions associated with a single gene based on flux variability analysis results and cleaning up metabolites that do not participate in the reaction, thereby achieving 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, adding the correspondence between gene and protein sequence and functional databases; at the reaction level, adding metabolic reaction and enzyme catalysis association information.

[0039] To facilitate understanding of the above technical solutions of this invention, the following further explains the technical solutions from the perspectives of architecture and principles, as follows: I. Dataset Construction 1. Construction of Biochemical Reaction Datasets Data Integration and Collection: Addressing the issues of existing tools' strong dependence 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, it expands upon the latest biochemical reaction data from the KEGG and MetaCyc databases, ensuring comprehensive and up-to-date data to improve 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. The BKMD metabolite 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 data 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, compartment information and annotation information. BiGG ID is concise and readable. For example, D-glucose is usually represented as "glc__D_c", where "glc" is the abbreviation for glucose, "__D" represents D-type glucose, "_c" represents cell compartment information, and c is the abbreviation for cytosol. BiGG ID is highly readable, which is the main reason why it is preferred in research. Compared with BiGG, MetaCyc reaction ID lacks standardization. For example, "CPD__45__16942_c" is the ID of 4-oxovalerate, which uses internal data labeling. KEGG also has similar problems. For example, rotenone in KEGG ID is "C07593". This internal numbering format requires querying the dataset, which is not conducive to researchers' identification. The processed metabolite data is mainly divided into 6 columns. The data structure is shown in Figure 2. The specific contents are 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. The Gene-Protein-Response Rule Dataset (GPR) 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 models from BiGG, KEGG, and MetaCyc, 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 includes 267,900 prokaryotic GPR rules, 364,007 eukaryotic GPR rules, and 722,510 archaea GPR rules. KEGG uses the ko_reaction.list and ko_genes.list files in the ko.tar file provided in the database to extract species-indiscriminate GPR rules by ko number, totaling 15,300,235 rules. The specific GPR rule dataset construction process is shown in Figure 3.

[0055] II. GEMs Construction 1. Gene Function Annotation This invention constructs GEMs from scratch based on protein sequences. First, a local alignment library is constructed 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] After the above process, a TSV format file can be obtained, which contains the comparison results between the FASTA file uploaded by the user and the local sequence data database. The result file has 12 columns. The data structure is shown in Figure 4, taking the comparison result of Escherichia coli K-12 strain as an example. The specific contents are 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 alignment score data in the third column is sorted in descending order to remove duplicate genes, and then higher-scoring matches are retained to ensure high-scoring results in the final product.

[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 Response Mapping Based on Annotation Results: Using the annotation results, biochemical responses can be mapped 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 they match, the entire row is extracted. This process is repeated for all gene IDs that meet the conditions, ultimately yielding the response set.

[0073] The 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 the reaction set in a DataFrame for later retrieval. By filtering with a dynamic thresholding program, the second element of the second column of the result, reaction ID, is used to extract the row in the reaction dataset that is the same as the id column. This result is the reaction set stored in reactions.csv. The data structure is shown in Figure 5. The specific information is as follows: (1) The id column contains the reaction ID.

[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 set. Based on the GPR rule dataset constructed above, the data structure is shown in Figure 6. The specific rules are as follows: The user sequence and the local sequence are aligned using DIAMOND. The score of the alignment between the local gene ID and the gene ID of the input sequence is obtained in the `gene` column. 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] Combining compartment 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 biological rationality of the model and the consistency of cell structure. The specific contents are as follows: (1) The gene column contains the gene IDs 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. Constructing the Biomass Equation Dataset The Biomass data in this invention primarily originates from the BiGG database. Specifically, the prokaryotic Biomass data mainly references models such as *Escherichia coli* iML1515, *Bacillus subtilis* iYO844, *Klebsiella pneumoniae* iYL1228, *Corynebacterium glutamicum* iCW773, and *Bacillus putida* iJN1463; the eukaryotic Biomass data references models including *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. Generating the Initial Model: A Python script automatically integrates 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 the initial model from the Excel file in SBML format.

[0090] 7. Filling Model Gaps: Due to incomplete gene annotation, missing metabolic databases, and incomplete experimental data, metabolic gaps may occur in the initial models constructed above during the GEMs construction process, resulting in the models failing to grow normally or making unreasonable simulations or 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. Removing Redundant Reactions: The GEMs constructed through the above steps possess all the functions and metrics of high-quality GEMs. However, since the extraction reactions are based on gene annotation, similar reactions catalyzed by the same enzyme or similar reactions from the database may be introduced. 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. Then, score the correctness of the reactions, for example, through 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 there is no impact, 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 Technology and Development Tools 1. Front-end Development Framework This website is developed 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 the help of PEAR, users can easily obtain GEMs models, calculate the optimal synthesis rate and metabolic flux distribution of specific products, and intuitively visualize and analyze the synthesis path, providing important guidance and support for metabolic engineering and biomanufacturing.

[0097] 2. Programming Implementation Tools and Software Packages The entire program and data processing of this invention are based on Python and its third-party libraries. See Table 4 for details and version information.

[0098] Table 4 Content and Version Information

[0099] IV. Results and Discussion 1. Development of Automated GEMs Construction Process This invention developed 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. Biochemical Reaction Dataset Construction 1) Data Cleaning and Standardization In the previous 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 checks were 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 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 occurred in the conservation of matter, leading to the unlimited generation of D-gluconic acid. Other errors due to non-conservation of matter include the unlimited generation of reducing power such as nadh and net ATP energy generation, all of which prevent the model from accurately quantifying the results. Furthermore, the mass-imbalanced reaction RE2477: amet_c + CE3087_c => ahcys_c + h_c + C09209_c will cause an infinite error in the production of nadh. The reaction atp_c + coa_c + ddca_c <=> amp_c + ddcacoa_c + ppi_c is mass-conserved, but it will cause an infinite error in the production of ATP; its actual reaction direction 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) Treatment of respiratory chain reactions: Respiratory chain reactions generate energy by forming a proton gradient, which may 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 lead to unlimited energy generation. The specific reactions treated 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. Genome Functional Annotation and Result Screening Workflow Design: PEAR employs a top-down strategy, constructing GEMs from scratch. First, similarity comparisons are performed using DIAMOND to score the sequence similarity between the input sequence and the alignment library, thereby performing gene function annotation based on the alignment score. In reality, due to the significant differences in genome size and sensitivity to alignment libraries among species, this invention develops a dynamic threshold strategy, automatically calculating 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 Reactions PEAR: Based on the above-screened annotation results, biochemical reactions are mapped using the BKMD dataset to obtain enzyme reactions and transport reactions. Then, based on the gene annotations for each enzyme reaction, GPR rules are assigned to the biochemical reactions using a GPR rule dataset. Finally, a set of enzyme reactions is generated.

[0115] 2) The spontaneous reaction set was obtained through gene annotation, and spontaneous reactions typically occur without enzyme catalysis or any external factors. Therefore, these reactions cannot be added to the reaction set through 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 normal growth simulation. Finally, the model was integrated with the enzyme reaction set.

[0116] Table 7 List of Spontaneous Reactions

[0117] 3) 3.1.3.3 Adding Cellular Sublocalization Information for Reactions This step is only used for constructing eukaryotic metabolic network models. The construction of eukaryotic metabolic network models is usually quite complex, as eukaryotes contain more complex 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 sublocalization information and add it to the model. First, we parse the user-uploaded sequence file, iterate through each sequence in the Python script, and obtain the cell sublocalization information one by one. The obtained subcellular localization information is saved as a DataFrame. Then, we parse and identify the cell sublocalization information string; the 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) Adding Culture Media: All the above-mentioned 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. They are generally classified into complex media (LB) and chemically synthesized media (M9). Different culture 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, because they lack the ability to synthesize certain substances, preventing normal growth, it is necessary to use specific culture media 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) Adding Biomass Equations: In this step, DIAMOND is 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. The specific Biomass workflow is shown in Figure 7.

[0125] 5. The initial model generation and constructed biochemical reaction dataset are 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, the 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 Automated Model Refinement Process 1) Development of a Weighted Algorithm for Gap Filling: The initial model may contain metabolic gaps, requiring gap-filling. To minimize the introduction of metabolic pathways and avoid redundancy, and to ensure 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 considered confidence level 1, 50% to 80% are considered confidence level 2, and the bottom 20% are considered 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 allows for the priority 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 make the constructed GEMs more accurate, this invention adds a redundant reaction removal step for automated model refinement. The process first loads the model and calculates the initial objective function value. The program then iterates through all reactions, removing circulating metabolites and identifying redundant reaction groups with the same substrate and product. Next, nested subsets within these redundant groups are 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. Then, 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 containing only ATP and ADP, as well as completely repetitive reactions. Finally, based on FVA analysis results, unnecessary redundant reactions associated with single genes are cyclically removed, and metabolites not involved in the reaction are cleaned from the model, resulting in a simplified metabolic model with complete information retention.

[0129] 3) Adding Annotations and Generating the Final Model: Annotations are added to the refined model to improve the cross-references between the reaction, metabolite, and gene databases. Database names and URLs are listed in Table 11. Finally, the model is saved in SBML, JSON, MAT, or EXCEL formats for user 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. Similarity comparisons are performed between the input sequences and the sequences in the local alignment library, and gene function annotations are 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 was performed by combining the biochemical reaction dataset and the gene-protein-reaction rule dataset. Based on the mapping results, a biochemical reaction set was constructed, and relevant metabolite dataset information was added to the biochemical reaction set. Species with the most similar metabolism to the input sequence were selected, and the local biomass equation corresponding to the similar species was mapped and added to the biochemical reaction set. The added biochemical reaction set was then integrated into the metabolite file, and an initial genome-scale metabolic network model was generated based on the metabolite file. The initial genome-scale metabolic network model was then subjected to gap filling, redundant reaction removal, and annotation information addition processing to obtain the final genome-scale metabolic network model.

2. The method for constructing a genome-scale metabolic network model according to claim 1, characterized in that, The process of integrating biochemical reaction data from a biological metabolism database and preprocessing the data to construct biochemical reaction datasets, metabolite datasets, and gene-protein-reaction rule datasets includes: integrating biochemical reaction data from a biological metabolism database and sequentially performing data cleaning and standardization, error correction, and respiratory chain reaction processing on the biochemical reaction data; constructing biochemical reaction datasets for prokaryotes, eukaryotes, and archaea based on the processed data; extracting metabolite information for each biochemical reaction in the biochemical reaction datasets and performing deduplication; constructing a metabolite dataset based on the deduplicated metabolite information; and extracting gene-protein-reaction rule data from the biochemical reaction datasets and obtaining the correspondence between reactions and gene IDs; constructing a gene-protein-reaction rule dataset based on the correspondence between reactions and gene IDs.

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 threshold method, the alignment score threshold of the input sequence is determined, and gene function annotation is screened in conjunction with the alignment score threshold. This includes: obtaining the protein sequence of the target species, constructing text-based bioinformatics standard format files for prokaryotes, eukaryotes, and archaea respectively, converting the text-based bioinformatics standard format files into the special format file of the sequence alignment software, and building a local alignment library; receiving the input sequence uploaded by the user, automatically identifying the type of the input sequence using the built-in functions and scripts of the sequence alignment software, and calling the amino acid sequence or nucleotide sequence model in the sequence alignment software to perform sequence similarity comparison with the local alignment library, outputting the sequence alignment result file, and completing the gene function annotation; detecting the sequence alignment result file, ensuring that the file exists and is not empty, sorting the annotation results in descending order according to the alignment score of the sequence alignment result file, and removing duplicate gene matches; based on the sorting distribution of the alignment scores, the 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 screened out.

5. The method for constructing a genome-scale metabolic network model according to claim 4, characterized in that, The method of determining the alignment score threshold of the input sequence based on the sorting distribution of the alignment score includes: counting the total number of genes after deduplication, selecting a preset percentage of the total number of genes as the threshold gene number, and storing the threshold gene number in the critical value variable; counting the total number of genes falling into the preset similarity score interval and calculating the sum, and determining the interval range into which the critical value falls based on the calculation result; determining whether the interval range into which the critical value falls belongs to the preset similarity score interval, and if so, using the lower bound of the interval range as the alignment score threshold, and if not, selecting the preset threshold as the alignment score threshold.

6. The method for constructing a genome-scale metabolic network model according to claim 1, characterized in that, The process involves mapping biochemical reactions based on the filtered 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. This includes: extracting gene IDs from sequence alignment results and selecting gene IDs that meet a threshold for searching in the gene-protein-reaction rule dataset to obtain the corresponding reaction ID and protein combination; matching rows with the same ID column in the biochemical reaction dataset based on the reaction ID and extracting the entire row; repeating the search and extraction operation for all gene IDs that meet the threshold, and summarizing to obtain the biochemical reaction set. The biochemical reaction set includes reaction ID, reaction name, and metric moments. The dataset includes: array, reaction upper and lower boundaries, gene-protein-reaction rules, annotation information, target reaction, and reaction equation. Based on the highest alignment score, all subunits of catalytic reactions are searched in the gene-protein-reaction rule dataset. The subunit compartments are searched in the Uniprot database, and the correct compartments are replaced to complete the effective subunit compartment mapping. Isoenzymes are linked using OR, and complexes composed of multiple subunits are linked to achieve Boolean connection logic reconstruction. Combining compartment information, biochemical reaction datasets, and gene-protein-reaction rule datasets, the catalytic metabolic reactions of gene-annotated enzymes are extracted to complete the adaptation of rules and reactions. The effective subunits and complex-related dataset information after compartment mapping are added to the biochemical reaction set.

7. 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.

8. The method for constructing a genome-scale metabolic network model according to claim 1, characterized in that, Gap-filling processing of the initial genome-scale metabolic network model includes: dividing the confidence level into several levels based on sequence alignment scores, assigning a weighting factor to the response corresponding to each confidence level, and filling gaps in the initial genome-scale metabolic network model based on the weighting factor.

9. The method for constructing a genome-scale metabolic network model according to claim 1, characterized in that, Redundant reaction removal from the initial genome-scale metabolic network model includes: traversing all reactions in the initial genome-scale metabolic network model, removing circulating metabolites in the reactions to identify redundant reaction groups with the same substrate and product; removing nested subsets from the redundant reaction groups and retaining the optimal reaction in each group based on mass balance, reversibility, and metabolite quantity indices; eliminating reactions containing only adenosine triphosphate (ATP) and adenosine diphosphate (ADP) and empty reactions from the redundant reaction groups; eliminating reactions not retained from the redundant reaction groups and merging the gene rules of the unretained reactions into the retained reactions; identifying and deleting reversible reactions containing only ATP and ADP and completely duplicated reaction formulas; and cyclically deleting unnecessary redundant reactions associated with single genes based on flux variability analysis results and cleaning up metabolites that do not participate in the reaction, thereby achieving redundant reaction removal from the initial genome-scale metabolic network model.

10. 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, adding the correspondence between gene and protein sequences and functional databases; at the reaction level, adding information on metabolic reactions and enzyme catalysis associations.

Images