Diagnostic model construction method based on novel network toxicology research normal form

Abstract

The invention discloses a diagnostic model construction method based on a novel network toxicology research normal form, and belongs to the field of health risk assess.The diagnostic model construction method specifically comprises the following steps that firstly, a potential poison and disease target library is constructed, and common potential targets between ATBC and fracture healing damage are screened and recognized; iI, constructing a PPI network to identify a core target, and analyzing functions and pathway enrichment of the core target; iII, identifying and screening the biomarkers, and constructing a diagnostic model to evaluate the diagnostic value of the biomarkers; according to the method, the evaluation precision of ATBC on the potential risk of fracture healing is remarkably improved, the conformation limitation of traditional static docking is broken through, stability evaluation in the time dimension is provided, the limitation of a rigid docking method is effectively overcome, meanwhile, visual representation of the dynamic behavior of a compound is provided, and the accuracy of the dynamic behavior of the compound is improved. The method provides a quantitative basis for mechanism analysis of molecular interaction, provides an effective technical means for subsequent high-risk substance screening, and is simple and convenient to operate and wide in application range.

Description

Technical Field

[0001] The present invention relates to the field of health risk assessment, and in particular to a method for constructing a diagnostic model based on a novel network toxicology research paradigm. Background Art

[0002] With the mass production and widespread use of plastic products, the concentration of microplastics in the environment has skyrocketed. Consequently, there is growing concern about exposure to microplastics and their potential risks to the environment and human health. Recent studies have shown that the concentration of microplastics in the human body, particularly in the brain, is increasing significantly. Plasticizers, widely used to soften plastics, have a strong positive correlation between their usage and plastic production. Acetyl tributyl citrate (ATBC), a plasticizer, has been widely used as a superior alternative to traditional phthalate plasticizers due to its excellent thermal and chemical stability, insolubility in water, and solubility in most organic solvents. It is widely used in food packaging materials (such as dairy products and beverage containers), medical PVC products (such as infusion tubing), and children's toys. Although ATBC is widely considered a safe plasticizer, reports have increasingly emerged on the potential risks of ATBC exposure to humans.

[0003] In recent decades, the amount of plastic waste has increased exponentially. As plastic waste is released into the environment, it may be deposited in the human body through various pathways and have deleterious effects on human health. However, the toxicological characteristics of this type of plastic waste remain largely understudied. The emerging field of network toxicology overcomes the limitations of traditional toxicology, such as limited technical tools. By using multiple databases to construct networks of relationships between compounds, toxicities, and targets, it can comprehensively and effectively assess the potential health risks of the growing number of environmental pollutants in the environment. With its advantages of accurate results, ease of use, wide applicability, and strong potential for in-depth exploration, network toxicology provides important technical support for predicting the potential toxicity of chemicals, elucidating their molecular mechanisms of action, and screening toxicants. Although ATBC is considered a safer alternative plasticizer, the potential harm of ATBC aggregation to the human body remains a concern due to prolonged exposure to high concentrations of ATBC. Currently, our understanding of the toxicological mechanisms of ATBC on fracture healing remains incomplete. To this end, we propose a method for constructing a diagnostic model based on this novel network toxicology research paradigm. Summary of the Invention

[0004] The purpose of the present invention is to solve the defects in the existing technology and propose a diagnostic model construction method based on a new network toxicology research paradigm.

[0005] In order to achieve the above object, the present invention adopts the following technical solutions:

[0006] A diagnostic model construction method based on a novel network toxicology research paradigm, the specific steps of the diagnostic model construction method are as follows:

[0007] I. Construct a library of potential toxicants and disease targets and screen to identify common potential targets between ATBC and impaired fracture healing;

[0008] II. Construct a PPI network to identify core targets and analyze their functions and pathway enrichment;

[0009] III. Identify and screen biomarkers and build diagnostic models to evaluate the diagnostic value of biomarkers;

[0010] IV. Analyze the infiltration level of each immune cell type and analyze the gene set enrichment of biomarkers, and then construct a regulatory network based on the biomarkers;

[0011] V. Online docking of the biomarker protein structure with the ATBC structure and dynamic simulation of the ligand-protein complex generated by docking.

[0012] As a further embodiment of the present invention, the specific steps of screening and identifying common potential targets between ATBC and impaired fracture healing in step I are as follows:

[0013] S1.1: Determine the standard structure and SMILES symbol of ATBC by searching for "acetyl tributyl citrate" in the PubChem database. Based on the determined standard structure and SMILES symbol of ATBC, set the keyword "acetyltributyl citrate" and limit the species to "Homo sapiens";

[0014] S1.2: Based on the set keywords and restrictions, potential targets of ATBC were obtained from the ChEMBL database. At the same time, the SMILE symbols of ATBC were uploaded to the Swiss Target Prediction database and the STITCH database to identify potential targets that have been overlooked.

[0015] S1.3: Integrate targets from different sources, remove duplicate data, and verify structural consistency. Use the UniProt database to standardize the obtained target names to establish a comprehensive ATBC target library to store all merged and verified targets;

[0016] S1.4: Search and extract targets corresponding to impaired fracture healing in the GeneCards and OMIM databases, and obtain datasets corresponding to impaired fracture healing by searching the keyword “fracture healing” in the GEO database;

[0017] S1.5: Download the GSE93215 dataset from the GPL6244 platform. The GSE93215 dataset contains peripheral blood microarray data from 19 patients with acute fractures and 9 healthy volunteers. Then, use the R language "limma" package to analyze the differentially expressed genes in the GSE93215 dataset, and screen out differentially expressed genes with a P value less than 0.05.

[0018] S1.6: Merge data from GeneCards and OMIM databases, remove duplicate data, and merge the merged results with the data from the GSE93215 dataset, retaining duplicate data, to establish a fracture healing impaired target library;

[0019] S1.7: Use the “VennDiagram” package to identify common potential targets between the ATBC target library and the impaired fracture healing target library. Then use the “ggplot2” package to visualize the results and design the intersection parts as potential targets of ATBC-induced impaired fracture healing.

[0020] As a further solution of the present invention, the specific steps of constructing the PPI network to identify core targets in step I are as follows:

[0021] S2.1: Potential target genes of ATBC-induced impaired fracture healing were entered into the STRING database, the species was restricted to "Homo sapiens", and the "minimum required interactions core" was set to "medium confidence > 0.4" for analysis;

[0022] S2.2: Import each set of data generated by STRING database analysis into Cytoscape software for visualization. Then, use the cytoHubba plug-in to identify core targets based on six algorithms: MCC, MNC, Degree, Closeness, Betweenness, and Stress.

[0023] S2.3: Upset plots were used to identify the top ten core targets in the intersection of the six algorithms, and the “circlize” package was used to perform chromosome localization analysis on the core targets.

[0024] As a further embodiment of the present invention, the specific steps of analyzing the function and pathway enrichment of the core target in step II are as follows:

[0025] S3.1: Use the “clusterProfiler” package to analyze core targets for Gene Ontology and Kyoto Encyclopedia of Genes and Genomes enrichment, and perform comprehensive GO analysis based on biological processes, cellular components, and molecular functions;

[0026] S3.2: KEGG analysis was used to identify the corresponding pathways of ATBC-induced fracture healing, and ATBC-induced fracture healing pathways with adjusted P values ​​lower than 0.05 were screened out. The screening results were visualized in the form of bar charts and bubble charts using the “ggplot2” package.

[0027] As a further embodiment of the present invention, the specific steps of constructing a diagnostic model to evaluate the diagnostic value of biomarkers in step III are as follows:

[0028] S4.1: Use the least absolute shrinkage and selection operator method LASSO and the random forest method RF to identify biomarkers, and build corresponding regression models using the "glmnet" package, and then determine the penalty coefficient through 10-fold cross-validation;

[0029] S4.2: The RF model was constructed using the “randomForest” package. The top four molecules were then screened based on the relative importance scores of the RF model. The intersection between the LASSO and RF algorithm results was then obtained using the “VennDiagram” package to identify key biomarkers for ATBC-induced impaired fracture healing.

[0030] S4.3: Use the GLM function to construct a binary logistic model, and then use the "rms" package to construct a nomogram diagnostic model to evaluate the diagnostic value of the biomarker. Then, based on the expression of the biomarker in the GSE93215 dataset, the diagnostic efficacy of the biomarker was verified by group comparison plots and receiver operating characteristic (ROC) curves;

[0031] S4.4: Based on the validation results, the nomogram diagnostic model was verified using the “rmda” and “pROC” packages, and the decision curve (DCA), precision-recall curve (PR), and ROC curve were visualized using the “ggplot2” package.

[0032] As a further embodiment of the present invention, the specific steps of analyzing the infiltration degree of each immune cell type and analyzing the gene set enrichment of biomarkers in step IV are as follows:

[0033] S5.1: Calculate the infiltration degree of 24 immune cell types based on the gene expression levels of 24 published immune cell gene sets using the single-sample gene set enrichment analysis method ssGSEA in the "GSVA" package, and evaluate the relationship between biomarkers and each immune cell type;

[0034] S5.2: Gene set enrichment analysis was performed using the “clusterProfiler” package to obtain the biological functions of each biomarker, and the “ggplot2” package was used for data visualization.

[0035] As a further embodiment of the present invention, the specific steps of constructing a regulatory network based on biomarkers in step IV are as follows:

[0036] S6.1: Obtain potential genes for biomarkers using the GeneMANIA database and construct a gene-gene interaction network;

[0037] S6.2: Obtain potential transcription factors (TFs) for biomarkers using the ChIPBase database and construct a gene-transcription factor interaction network;

[0038] S6.3: Obtain potential miRNAs and lncRNAs for biomarkers through the starBase database, and construct a competitive endogenous RNA regulatory network of mRNA-miRNA-lncRNA, namely the ceRNA regulatory network.

[0039] As a further embodiment of the present invention, the specific steps of performing dynamic simulation on the ligand-protein complex generated by docking in step V are as follows:

[0040] S7.1: Download the biomarker protein structure from the RCSB PDB database and obtain the ATBC structure from PubChem. Use PyMOL software to remove water molecules and native ligands to preprocess the protein and save the PyMOL processing results as a PDB file.

[0041] S7.2: Use the molecular docking online website CB-Dock2 to dock the biomarker protein structure with the ATBC structure online. Then use PyMOL software to compare and visualize the binding patterns between the biomarker and ATBC.

[0042] S7.3: Molecular dynamics simulations were performed using Gromacs software. AmberTools22 was used to add the GAFF force field to the small molecules. Gaussian16W was used to perform hydrogenation and calculate the RESP potential to complete the pretreatment of the small molecules.

[0043] S7.4: Add the generated potential data to the molecular dynamics system topology file. Set the simulation conditions to a static temperature of 300 K, a pressure of 1 bar, the amber99sb-ildn force field, and water as the solvent. Run the simulation. During the simulation, add the required amount of Na+ to neutralize the total charge of the simulated system.

[0044] S7.5: Minimize the energy using the steepest descent method, followed by 100,000 steps of isothermal and isochoric ensemble equilibrium NVT and isothermal and isobaric ensemble equilibrium NPT, with a steepest descent coupling constant of 0.1 ps and a duration of 100 ps, ​​followed by a free molecular dynamics simulation with a duration of 5,000,000 steps, a step size of 2 fs, and a total duration of 100 ns.

[0045] S7.6: After the simulation is completed, the trajectory is analyzed using Gromacs software to calculate the root mean square deviation (RMSD), root mean square fluctuation (RMSF), solvent accessible surface area (SASA), protein radius of gyration (Rg), binding free energy, and free energy map (FEL) of each amino acid trajectory.

[0046] Compared with the prior art, the present invention has the following beneficial effects:

[0047] This diagnostic model construction method based on the new network toxicology research paradigm established a library of ATBC and fracture healing-related targets through multiple databases, and screened the intersection of the two as common potential targets. Then, it used protein interaction networks and multiple algorithms to identify core targets, and performed functional and pathway enrichment analysis to explore their biological mechanisms. Subsequently, diagnostic biomarkers were screened through LASSO and random forest models, and a logistic regression diagnostic model was constructed and verified. After that, the immune infiltration characteristics were analyzed, the regulatory network was established, and the binding stability of ATBC and key targets was verified through molecular docking and molecular dynamics simulation. This significantly improved the accuracy of assessing the potential risk of ATBC to fracture healing, broke through the conformational limitations of traditional static docking, provided stability assessment in the time dimension, and effectively overcame the limitations of rigid docking methods. At the same time, it not only provided a visual representation of the dynamic behavior of the complex, but also provided a quantitative basis for the analysis of the mechanism of molecular interactions, providing an effective technical means for the subsequent screening of high-risk substances. It is easy to operate and has a wide range of applications. BRIEF DESCRIPTION OF THE DRAWINGS

[0048] The accompanying drawings are used to provide further understanding of the present invention and constitute a part of the specification. They are used to explain the present invention together with the embodiments of the present invention and do not constitute a limitation of the present invention.

[0049] Figure 1 This is a flowchart of a diagnostic model construction method based on a novel network toxicology research paradigm proposed in the present invention;

[0050] Figure 2 This is a design flow chart of a diagnostic model construction method based on a novel network toxicology research paradigm proposed in the present invention;

[0051] Figure 3This is a Venn diagram of the intersection of ATBC and impaired fracture healing targets for a diagnostic model construction method based on a novel network toxicology research paradigm proposed in the present invention;

[0052] Figure 4 This is a schematic diagram of the core target screening of a diagnostic model construction method based on a novel network toxicology research paradigm proposed in the present invention;

[0053] Figure 5 This is a schematic diagram of GO and KEGG analysis of core targets of a diagnostic model construction method based on a new network toxicology research paradigm proposed in the present invention;

[0054] Figure 6 This is a schematic diagram of diagnostic biomarker screening based on a diagnostic model construction method proposed in the present invention based on a novel network toxicology research paradigm;

[0055] Figure 7 This is a schematic diagram of the construction and verification results of a diagnostic model based on a novel network toxicology research paradigm proposed in the present invention;

[0056] Figure 8 This is a schematic diagram of immune infiltration and GSEA analysis results for a diagnostic model construction method based on a novel network toxicology research paradigm proposed in the present invention;

[0057] Figure 9 This is a schematic diagram of the regulatory network of a diagnostic model construction method based on a novel network toxicology research paradigm proposed in the present invention;

[0058] Figure 10 This is a schematic diagram of the molecular docking validation results of a diagnostic model construction method based on a novel network toxicology research paradigm proposed in the present invention;

[0059] Figure 11 This is a schematic diagram of the HDAC2 molecular dynamics simulation validation results of a diagnostic model construction method based on a novel network toxicology research paradigm proposed in the present invention;

[0060] Figure 12 This is a schematic diagram of the TP53 molecular dynamics simulation validation results of a diagnostic model construction method based on a novel network toxicology research paradigm proposed in the present invention;

[0061] Figure 13 This is a free energy profile diagram of a diagnostic model construction method based on a new network toxicology research paradigm proposed in the present invention. DETAILED DESCRIPTION

[0062] Example 1, with reference to Figure 1-13A diagnostic model construction method based on a new network toxicology research paradigm is described. The specific steps of the diagnostic model construction method are as follows:

[0063] A library of potential toxicants and disease targets was constructed and screened to identify common potential targets between ATBC and impaired fracture healing.

[0064] Specifically, the standard structure and SMILES symbol of ATBC were determined by querying "acetyl tributyl citrate" in the PubChem database. Based on the determined standard structure and SMILES symbol of ATBC, the keyword "acetyltributyl citrate" was set, and the species was restricted to "Homo sapiens". According to the set keywords and restrictions, potential targets of ATBC were obtained from the ChEMBL database. At the same time, the SMILE symbol of ATBC was uploaded to the Swiss TargetPrediction database and the STITCH database to identify potential targets that were overlooked. Targets from different sources were integrated, duplicate data were deleted, and structural consistency was verified. The UniProt database was used to standardize the obtained target names to establish a comprehensive ATBC target library, which stored all merged and verified targets. Targets corresponding to impaired fracture healing were retrieved and extracted in the GeneCards and OMIM databases, and the keywords "fracture" were searched in the GEO database. "healing", to obtain the corresponding dataset of impaired fracture healing. The GSE93215 dataset was downloaded from the GPL6244 platform. The GSE93215 dataset contains peripheral blood microarray data of 19 patients with acute fractures and 9 healthy volunteers. The "limma" package of R language was then used to analyze the differentially expressed genes in the GSE93215 dataset, and the differentially expressed genes with a P value less than 0.05 were screened out. The data from the GeneCards and OMIM databases were merged, and the duplicate data were removed. The merged results were merged with the data of the GSE93215 dataset, and the duplicate data were retained to establish a target library for impaired fracture healing. The "VennDiagram" package was used to identify common potential targets between the ATBC target library and the impaired fracture healing target library. The "ggplot2" package was then used to visualize the results, and the intersection was designed as a potential target of ATBC-induced impaired fracture healing.

[0065] It should be further explained that the reference Figure 3It can be seen that 362 ATBC-related targets were obtained by integrating the ChEMBL, SwissTargetPrediction and STITCH databases, and 3,730 fracture healing-related targets were obtained through the GeneCards and OMIM databases. 4,903 DEGs were screened out through the GSE93215 dataset. The results of the integrated database and dataset obtained 834 fracture healing-related targets. Subsequently, 132 cross-targets shared by ATBC and OA were identified. Each group of overlapping targets is a potential candidate for the effect of ATBC on impaired fracture healing.

[0066] A PPI network was constructed to identify core targets, and the functions and pathway enrichment of core targets were analyzed.

[0067] Specifically, the potential target genes of ATBC-induced impaired fracture healing were entered into the STRING database, the species was limited to "Homo sapiens", and the "minimum required interactions core" was set to "medium confidence > 0.4" for analysis. The data sets generated by the STRING database analysis were imported into the Cytoscape software for visualization. The cytoHubba plug-in was then used to identify core targets based on six algorithms, namely MCC, MNC, Degree, Closeness, Betweenness, and Stress. The Upset graph was used to identify the intersection core targets ranked in the top ten among the six algorithms, and the "circlize" package was used to perform chromosomal localization analysis of the core targets.

[0068] Specifically, the "clusterProfiler" package was used to analyze core targets for gene ontology and Kyoto Encyclopedia of Genes and Genomes enrichment, and a comprehensive GO analysis was performed based on biological processes, cellular components, and molecular functions. The corresponding pathways of ATBC-induced fracture healing were identified through KEGG analysis, and ATBC-induced fracture healing pathways with adjusted P values ​​below 0.05 were screened out. The screening results were visualized in the form of bar charts and bubble charts using the "ggplot2" package.

[0069] In addition, reference Figure 4-5It can be seen that it should be noted that a PPI interaction network with 47 nodes and 382 edges was generated through the STRING database, and the top 10 potential targets of MCC, MNC, Degree, Closeness, Betweenness and Stress were generated through the cytoHubba plug-in. Finally, based on 6 algorithms, a total of 6 intersection targets were obtained, namely HDAC2, KAT2B, TP53, HDAC3, SMARCA4 and EP300. The 6 core targets were then enriched. A total of 573 statistically significant GO entries were generated through set analysis, including 479 BP entries, 25 CC entries and 69 MF entries. Among them, in the BP category, targets were mainly enriched in protein deacetylation, regulation of glucose metabolism, and androgen receptor signaling pathway; in the CC category, targets were mainly enriched in transcriptional regulatory complexes, histone deacetylase complexes, ATPase complexes, etc.; in the MF category, targets were mainly enriched in DNA binding transcription factor binding, histone deacetylase binding, NF-κB binding, etc.

[0070] In addition, a total of 19 signaling pathways with statistical differences were enriched, including the core targets of ATBC and fracture healing, which are mainly related to the thyroid hormone signaling pathway, Notch signaling pathway, cell cycle, miRNA in cancer, Wnt signaling pathway and neutrophil extracellular trap formation.

[0071] Identify and screen biomarkers, and build diagnostic models to evaluate the diagnostic value of biomarkers.

[0072] Specifically, the least absolute shrinkage and selection operator method LASSO and the random forest method RF were used to identify biomarkers, and the corresponding regression models were established using the "glmnet" package. The penalty coefficient was then determined through 10 cross-validations, and the RF model was constructed using the "randomForest" package. The top four molecules were then screened based on the relative importance score of the RF model. The intersection between the LASSO and RF algorithm results was obtained using the "VennDiagram" package to identify key biomarkers for ATBC-induced impaired fracture healing. A binary logistic model was constructed using the GLM function, and a nomogram diagnostic model was constructed using the "rms" package to evaluate the diagnostic value of the biomarkers. The diagnostic efficacy of the biomarkers was then verified using group comparison plots and receiver operating characteristic (ROC) curves based on their expression in the GSE93215 dataset. Based on the verification results, the nomogram diagnostic model was validated using the "rmda" and "pROC" packages, and the decision curve DCA, precision-recall curve PR, and ROC curve were visualized using the "ggplot2" package.

[0073] It should be further explained that the reference Figure 6-7 As shown, the LASSO algorithm identified two biomarkers, while the RF algorithm identified the top four genes within the core target. The intersection of the results from these two algorithms then identified two diagnostic biomarkers. Furthermore, a nomogram model incorporating the gene scores of the two biomarkers was constructed and validated for HDAC2 and TP53. HDAC2 showed significant differential expression between the disease and control groups, with upregulation in the disease group and an area under the curve (AUC) of 0.930. TP53 also showed significant differential expression between the disease and control groups, with downregulation in the disease group and an AUC of 0.906. The constructed diagnostic model was then validated. The decision curve demonstrated a good net benefit, while the PR and ROC curves demonstrated good diagnostic efficacy.

[0074] Analyze the infiltration degree of each immune cell type and the gene set enrichment of biomarkers, and then construct a regulatory network based on the biomarkers.

[0075] Specifically, the single-sample gene set enrichment analysis method ssGSEA in the "GSVA" package was used to calculate the infiltration degree of 24 immune cell types based on the gene expression levels of 24 published immune cell gene sets, and the relationship between biomarkers and each immune cell was evaluated. The "clusterProfiler" package was used for gene set enrichment analysis to obtain the biological function of each biomarker, and the "ggplot2" package was used for data visualization.

[0076] In addition, reference Figure 8 The ssGSEA algorithm was used to evaluate immune cell infiltration between the disease and control groups in the impaired fracture healing dataset. The results showed that the number of neutrophils was increased in the disease group, while the number of aDCs, cytotoxic cells, NKCD56bright cells, NKCD56dim cells, T cells, TFH, Tgd, Th17 cells, and TReg cells was decreased. In addition, two diagnostic biomarkers were closely associated with various immune cells. GSEA analysis of the two diagnostic biomarkers was performed to evaluate the potential impact of HDAC2 and TP53 on the occurrence and progression of impaired fracture healing. The results showed that the HDAC2 gene was mainly enriched in RNA polymerase II transcription termination and cornified envelope formation, while the TP53 gene was mainly enriched in T cell receptor pathways and secretory factors.

[0077] Specifically, the potential genes of biomarkers were obtained through the GeneMANIA database, and a gene-gene interaction network was constructed. The potential transcription factors TFs of biomarkers were obtained through the ChIPBase database, and a gene-transcription factor interaction network was constructed. The potential miRNAs and lncRNAs of biomarkers were obtained through the starBase database, and a competitive endogenous RNA regulatory network of mRNA-miRNA-lncRNA, namely, the ceRNA regulatory network, was constructed.

[0078] It should be further explained that the reference Figure 9 As can be seen, the gene-gene interaction network of biomarkers has a central hub node representing these genes, surrounded by 20 nodes (FAU, MDM2, MCFD2, MYB, EXOC6B, SIN3A, RCOR1, CIR1, PHF21A, COP1, TP53BP1, HMG20B, TP53BP2, KDM1A, PRDM4, PPARD, HDAC1, MDM4, TP63 and HIF1A), representing genes significantly associated with them. In addition, the network highlights the top seven related functions associated with these genes, such as the intrinsic apoptosis signaling pathway, protein deacetylase activity and inhibitory transcription factor binding;

[0079] The gene-transcription factor interaction network of biomarkers, with the central hub node representing these genes, surrounded by 19 nodes, representing TFs significantly associated with them, of which 11 TFs related to HDAC2 and 12 TFs related to TP53 were predicted;

[0080] The ceRNA regulatory network of biomarker mRNA, miRNA, and lncRNA was constructed. Through miRNA prediction, three miRNAs with potential associations with the two biomarkers were discovered, including hsa-miR-494-3p, hsa-miR-218-5p, and hsa-miR-3163. Furthermore, lncRNA prediction was performed on the discovered miRNAs, and 21 lncRNAs with potential associations with the miRNAs were identified.

[0081] The biomarker protein structure was docked with the ATBC structure online, and the ligand-protein complex generated by docking was dynamically simulated.

[0082] Specifically, the biomarker protein structure was downloaded from the RCSB PDB database, and the ATBC structure was obtained from PubChem. The water molecules and original ligands were removed using PyMOL software to pre-process the protein, and the PyMOL software processing results were saved as PDB files. The biomarker protein structure was docked with the ATBC structure online using the molecular docking online website CB-Dock2. The binding mode of the biomarker and ATBC was compared and visualized using PyMOL software. Molecular dynamics simulation was performed using Gromacs software. The GAFF force field was added to the small molecule through AmberTools22. At the same time, Gaussian16W was used to perform hydrogenation and calculate the RESP potential of the small molecule to complete the pre-processing of the small molecule. The generated potential data was added to the molecular dynamics system topology file, and the simulation conditions were set to a static temperature of 300 K, 1 bar atmospheric pressure, amber99sb-ildn force field, and water molecules as solvent, and simulations were carried out. During the simulation, the total charge of the simulation system was neutralized by adding the required amount of Na+, and the steepest descent method was used for energy minimization. Then, 100,000 steps of isothermal and isovolumetric ensemble equilibrium NVT and isothermal and isobaric ensemble equilibrium NPT were performed respectively. Among them, the coupling constant of the steepest descent method was 0.1 ps, the duration was 100 ps, ​​and then the free molecular dynamics simulation was run with a simulation process of 5,000,000 steps, a step size of 2 fs, and a total duration of 100 ns. After the simulation was completed, Gromacs software was used to analyze the trajectory, and the root mean square deviation RMSD, root mean square fluctuation RMSF, solvent accessible surface area SASA, protein gyration radius, Rg, binding free energy, and free energy map FEL of each amino acid motion trajectory were calculated.

[0083] As a further explanation of this embodiment, refer to Figure 10-12 Molecular docking technology was used to explore the interaction between ATBC and two diagnostic proteins. The results showed that the binding energy of HDAC2 was -4.8kcal / mol and that of TP53 was -5.2kcal / mol. Both diagnostic biomarkers showed strong binding to ATBC, with binding energies less than -4.25kcal / mol. This indicates that ATBC and each marker can spontaneously bind, indicating that they play an important role in the molecular mechanism of ATBC-induced fracture healing.

[0084] Meanwhile, in the molecular dynamics simulation, the RMSD curve of the HDAC2 / TP53-ATBC complex showed equilibrium after 20 ns, with average RMSD values ​​of 0.18 nm and 0.17 nm, respectively. RMSF analysis identified key flexible regions in HDAC2 and TP53 proteins, especially around amino acid residues 0-50 and 300-350 of the HDAC2-ATBC complex and around amino acid residues 1500-1600 of the TP53-ATBC complex, especially around 1530 and 1570, with corresponding SASA values ​​of 165 nm. 2 and 72nm 2 The Rg values ​​of the two complexes remained stable, with average values ​​of 2.02 nm and 1.44 nm, respectively. Among the dissociation free energies of the HDAC2-ATBC complex, amino acid A:PHE:206 had the highest affinity for the molecule, with a specific value of -3.083 kcal / mol. Among the dissociation free energies of the TP53-ATBC complex, amino acid A:GLN:1577 had the highest affinity for the molecule, with a specific value of -1.468 kcal / mol. After the complex system stabilized, we used the MM / GBSA method to calculate the average binding free energies of the two groups of proteins and small molecules. As shown in the table below, the average free energies of binding between HDAC2 and TP53 and ATBC were -24.89 and -19.86 kcal / mol, respectively. The results showed that HDAC2 and TP53 had good binding energies and high binding affinities with ATBC, indicating that the binding of HDAC2, TP53 and ATBC is thermodynamically favorable, mainly driven by hydrophobic effects and electrostatic interactions, further validating the results of molecular docking.

[0085] Binding free energy calculation results

[0086]

[0087] In the table, ΔVDWAALS: van der Waals energy; ΔE elec : electrostatic energy; ΔE GB : polar solvation energy; ΔE surf : nonpolar solvation energy; ΔG gas :Molecular mechanics energy (meteorological energy) = ΔVDWAALS + ΔE elec ;ΔG solvation : Solvation energy = ΔEGB + ΔEsurf; ΔTOTAL: Total energy = ΔG gas +ΔG solvation .

[0088] refer to Figure 13It can be seen that the PC1 and PC2 obtained by principal component analysis of the RMSD value and Rg value of the small molecule and protein receptor docking complex are used as coordinate axes, and the Gibbs relative free energy is used as the Z axis to draw a three-dimensional free energy topography. This graphical representation method not only provides rich visual information, but also helps to gain a deeper understanding of the dynamic behavior and stability of the complex.

[0089] Figure 13 The stable conformations of the complex are prominently highlighted in the simulation. These conformations cluster within the minimum free energy region at low energies. The smoothness and depth of this region reflect the stability of the complex. If the protein-ligand interaction is weak or unstable, the free energy landscape will exhibit multiple, rough-surfaced minimum energy clusters, typically implying multiple possible stable states, with frequent or facile transitions between these states. However, in both systems, the free energy landscapes of the small molecule-protein receptor docking complex exhibited a nearly single, sharp minimum energy cluster. This indicates that the complex maintained a highly stable state during the simulation. These results, combined with the RMSD curves of the complex, further confirm the exceptional stability of the HDAC2 / TP53-ATBC complex. This stability is crucial for studying the potential impact of ATBC on fracture healing, as it implies that the complex maintains its structural and functional integrity under physiological conditions, allowing it to more effectively exert its biological activity.

Claims

1. A diagnostic model construction method based on a new network toxicology research paradigm, characterized in that: The specific steps of the diagnostic model construction method are as follows: I. Construct a library of potential toxicants and disease targets and screen to identify common potential targets between ATBC and impaired fracture healing; II. Construct a PPI network to identify core targets and analyze their functions and pathway enrichment; III. Identify and screen biomarkers and build diagnostic models to evaluate the diagnostic value of biomarkers; IV. Analyze the infiltration level of each immune cell type and analyze the gene set enrichment of biomarkers, and then construct a regulatory network based on the biomarkers; V. Online docking of the biomarker protein structure with the ATBC structure and dynamic simulation of the ligand-protein complex generated by docking.

2. A diagnostic model construction method based on a novel network toxicology research paradigm according to claim 1, characterized in that: The specific steps for screening and identifying common potential targets between ATBC and impaired fracture healing described in Step I are as follows: S1.1: Determine the standard structure and SMILES symbol of ATBC by searching for "acetyl tributyl citrate" in the PubChem database. Based on the determined standard structure and SMILES symbol of ATBC, set the keyword "acetyltributyl citrate" and limit the species to "Homo sapiens"; S1.2: Based on the set keywords and restrictions, potential targets of ATBC were obtained from the ChEMBL database. At the same time, the SMILE symbols of ATBC were uploaded to the Swiss Target Prediction database and the STITCH database to identify potential targets that have been overlooked. S1.3: Integrate targets from different sources, remove duplicate data, and verify structural consistency. Use the UniProt database to standardize the obtained target names to establish a comprehensive ATBC target library to store all merged and verified targets; S1.4: Search and extract targets corresponding to impaired fracture healing in the GeneCards and OMIM databases, and obtain datasets corresponding to impaired fracture healing by searching the keyword "fracture healing" in the GEO database; S1.5: Download the GSE93215 dataset from the GPL6244 platform. The GSE93215 dataset contains peripheral blood microarray data from 19 patients with acute fractures and 9 healthy volunteers. Then, use the "limma" package in R to analyze the differentially expressed genes in the GSE93215 dataset, and screen out differentially expressed genes with a P value less than 0.

05. S1.6: Merge data from GeneCards and OMIM databases, remove duplicate data, and merge the merged results with the data from the GSE93215 dataset, retaining duplicate data, to establish a fracture healing impaired target library; S1.7: Use the "VennDiagram" package to identify common potential targets between the ATBC target library and the impaired fracture healing target library. Use the "ggplot2" package to visualize the results and designate the intersections as potential targets of ATBC-induced impaired fracture healing.

3. A diagnostic model construction method based on a novel network toxicology research paradigm according to claim 2, characterized in that: The specific steps for constructing the PPI network to identify core targets described in Step I are as follows: S2.1: Potential target genes of ATBC-induced impaired fracture healing were entered into the STRING database, the species was restricted to "Homo sapiens", and the "minimum required interactions core" was set to "medium confidence > 0.4" for analysis; S2.2: Import each set of data generated by STRING database analysis into Cytoscape software for visualization. Then, use the cytoHubba plug-in to identify core targets based on six algorithms: MCC, MNC, Degree, Closeness, Betweenness, and Stress. S2.3: Upset plots were used to identify the top ten core targets in the intersection of the six algorithms, and the "circlize" package was used to analyze the chromosome localization of the core targets.

4. A diagnostic model construction method based on a novel network toxicology research paradigm according to claim 3, characterized in that: The specific steps for analyzing the function and pathway enrichment of core targets described in Step II are as follows: S3.1: Use the "clusterProfiler" package to analyze core targets for Gene Ontology and Kyoto Encyclopedia of Genes and Genomes enrichment, and perform comprehensive GO analysis based on biological processes, cellular components, and molecular functions; S3.2: KEGG analysis was used to identify the corresponding pathways of ATBC-induced fracture healing, and ATBC-induced fracture healing pathways with adjusted P values ​​less than 0.05 were screened out. The screening results were visualized in the form of bar charts and bubble charts using the "ggplot2" package.

5. The method for constructing a diagnostic model based on a novel network toxicology research paradigm according to claim 4, characterized in that: The specific steps for constructing a diagnostic model to evaluate the diagnostic value of biomarkers in step III are as follows: S4.1: Use the least absolute shrinkage and selection operator method LASSO and the random forest method RF to identify biomarkers, and build corresponding regression models using the "glmnet" package. Then, determine the penalty coefficient through 10-fold cross-validation; S4.2: The RF model was constructed using the "randomForest" package. The top four molecules were then screened based on the relative importance scores of the RF model. The intersection of the LASSO and RF algorithm results was then obtained using the "VennDiagram" package to identify key biomarkers for ATBC-induced impaired fracture healing. S4.3: Use the GLM function to construct a binary logistic model, and then use the "rms" package to construct a nomogram diagnostic model to evaluate the diagnostic value of the biomarker. Then, based on the expression of the biomarker in the GSE93215 dataset, the diagnostic efficacy of the biomarker was verified using group comparison plots and receiver operating characteristic (ROC) curves. S4.4: Based on the validation results, verify the nomogram diagnostic model using the "rmda" and "pROC" packages, and use the "ggplot2" package to visualize the decision curve (DCA), precision-recall curve (PR), and receiver operating characteristic (ROC) curve.

6. A diagnostic model construction method based on a novel network toxicology research paradigm according to claim 5, characterized in that: The specific steps for analyzing the infiltration level of each immune cell type and analyzing the gene set enrichment of biomarkers described in step IV are as follows: S5.1: Calculate the infiltration degree of 24 immune cell types based on the gene expression levels of 24 published immune cell gene sets using the single-sample gene set enrichment analysis method ssGSEA in the "GSVA" package, and evaluate the relationship between biomarkers and each immune cell type; S5.2: Gene set enrichment analysis was performed using the "clusterProfiler" package to obtain the biological functions of each biomarker, and the "ggplot2" package was used for data visualization.

7. A diagnostic model construction method based on a novel network toxicology research paradigm according to claim 6, characterized in that: The specific steps for constructing a regulatory network based on biomarkers described in step IV are as follows: S6.1: Obtain potential genes for biomarkers using the GeneMANIA database and construct a gene-gene interaction network; S6.2: Obtain potential transcription factors (TFs) for biomarkers using the ChIPBase database and construct a gene-transcription factor interaction network; S6.3: Obtain potential miRNAs and lncRNAs for biomarkers through the starBase database, and construct a competitive endogenous RNA regulatory network of mRNA-miRNA-lncRNA, namely the ceRNA regulatory network.

8. The method for constructing a diagnostic model based on a novel network toxicology research paradigm according to claim 7, characterized in that: The specific steps of performing dynamic simulation on the ligand-protein complex generated by docking in step V are as follows: S7.1: Download the biomarker protein structure from the RCSB PDB database and obtain the ATBC structure from PubChem. Use PyMOL software to remove water molecules and native ligands to preprocess the protein and save the PyMOL processing results as a PDB file. S7.2: Use the molecular docking online website CB-Dock2 to dock the biomarker protein structure with the ATBC structure online. Then use PyMOL software to compare and visualize the binding patterns between the biomarker and ATBC. S7.3: Molecular dynamics simulations were performed using Gromacs software. AmberTools22 was used to add the GAFF force field to the small molecules. Gaussian16W was used to perform hydrogenation and calculate the RESP potential to complete the pretreatment of the small molecules. S7.4: Add the generated potential data to the molecular dynamics system topology file. Set the simulation conditions to a static temperature of 300 K, a pressure of 1 bar, the amber99sb-ildn force field, and water as the solvent. Run the simulation. During the simulation, add the required amount of Na+ to neutralize the total charge of the simulated system. S7.5: Minimize the energy using the steepest descent method, followed by 100,000 steps of isothermal and isochoric ensemble equilibrium NVT and isothermal and isobaric ensemble equilibrium NPT, with a steepest descent coupling constant of 0.1 ps and a duration of 100 ps, ​​followed by a free molecular dynamics simulation with a duration of 5,000,000 steps, a step size of 2 fs, and a total duration of 100 ns. S7.6: After the simulation is completed, the trajectory is analyzed using Gromacs software to calculate the root mean square deviation (RMSD), root mean square fluctuation (RMSF), solvent accessible surface area (SASA), protein radius of gyration (Rg), binding free energy, and free energy map (FEL) of each amino acid trajectory.

Images