Method and system for identifying potential target spot in complex disease based on network

Abstract

The invention provides a method ComplexDnet for identifying a disease target based on a network for the first time. According to the method, on the basis of a WGCNA result, a protein interaction network is constructed, and structural key analysis, robustness analysis and functional consistency analysis are further performed on the network to obtain three scores respectively. And constructing a comprehensive score algorithm according to the three scores, and sorting and screening the target spots according to the comprehensive score so as to obtain key treatment target spots. According to the method, the discovery period from data to candidate drugs can be shortened, the prediction accuracy is high, and the result is accurate.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of bioinformatics and network analysis, and in particular, to a method and system for identifying potential targets in complex diseases based on networks. BACKGROUND

[0002] Complex diseases such as metabolic-associated steatohepatitis (MASH), diabetes, cancer, Alzheimer's disease, etc. have multiple etiologies and complex molecular mechanisms. Traditional drug development for these diseases often fails due to insufficient understanding of the pathogenesis of the disease, inaccurate selection of drug targets, and inefficient evaluation of drug activity, etc.

[0003] In recent years, with the development of high-throughput sequencing technology and artificial intelligence (AI), various computational-assisted drug discovery methods have been proposed, mainly including differential expression gene (DEG) analysis, weighted gene co-expression network analysis (WGCNA), network-based inference (NBI), network topology analysis, and deep learning methods, etc. At present, the most commonly used and mainstream target screening process in bioinformatics analysis is to construct a gene co-expression network using WGCNA. By identifying the gene module with the strongest correlation with the disease phenotype, the key module is analyzed for functional enrichment (such as GO or KEGG analysis) to understand its main biological function. Within the key module, the gene with the highest connectivity (degree) is selected as the "hub gene". These hub genes are used as potential key targets for subsequent literature research or experimental validation. WGCNA can reduce the research scope from tens of thousands of genes to a few hundred genes, with high efficiency and good accuracy.

[0004] However, the method of using WGCNA to identify key targets of diseases has the following shortcomings: (1) the evaluation standard of the target is too simple, mainly relying on "connectivity" to define the hub gene, which may miss the real "key" target such as the node that "bridges" two sub-networks; (2) there is a lack of systematic evaluation of network robustness, and the importance of nodes is not evaluated from the perspective of "network disintegration" or "maintaining network integrity", which cannot quantify the degree of damage to the entire disease-related network after a target is attacked (such as being inhibited by a drug); (3) the biological function of the target itself is not quantitatively compared with the core function of the entire disease module, resulting in an inaccurate evaluation of the matching degree of the function; (4) there is a lack of a unified quantitative scoring system, and there is a lack of a comprehensive and quantifiable score to integrate various important attributes of the target (such as structural importance and functional consistency), making the priority ranking of the targets not scientific and reliable enough.

[0005] Therefore, there is an urgent need in the art for a multi-dimensional and comprehensive target importance evaluation system, which can more accurately identify the real core target in the disease network by integrating the keyness of network topology, network robustness and functional consistency. SUMMARY

[0006] The purpose of the present application is to provide a network-based method and system for identifying potential targets in complex diseases.

[0007] In a first aspect, the present application provides a network-based method for identifying disease targets, comprising the steps of:

[0008] (1) providing transcriptome data, protein-protein interaction network data and functional annotation data of a disease to be tested; the transcriptome data includes gene expression data and phenotype data of the disease;

[0009] (2) performing weighted gene co-expression network analysis, comprising the steps of:

[0010] (2a) calculating the correlation between each pair of genes according to the transcriptome data using a correlation algorithm to obtain a correlation value; taking each gene as a node of the network and taking the value obtained by performing soft thresholding on the correlation value as an edge of the network, constructing a co-expression network;

[0011] (2b) dividing the co-expression network into multiple gene modules using a clustering algorithm;

[0012] (2c) calculating the correlation between the module characteristic genes in each gene module and the phenotype data using a correlation algorithm to obtain a correlation value; selecting a module with a strong correlation value as a key disease module;

[0013] (3) identification and evaluation of key targets, comprising the steps of:

[0014] (3a) mapping all genes in the key disease module obtained in step (2) to the protein-protein interaction network data, extracting a subnetwork composed of the genes and their interaction relationships, and obtaining a disease-specific protein-protein interaction network;

[0015] (3b) identifying potential targets in the disease-specific protein-protein interaction network using a network centrality algorithm and / or a network disassembly algorithm, calculating the change in the largest connected component (LCC) after removing the potential targets, and obtaining a structural keyness score;

[0016] (3c) removing several nodes in the disease-specific protein-protein interaction network using a machine learning algorithm to minimize the cumulative normalized connectivity (R), and then gradually introducing the removed nodes to calculate the change in the cumulative normalized connectivity, and normalizing the change to obtain a robustness score;

[0017] (3d) performing functional enrichment analysis on the key disease module to obtain a function list of the key disease module; querying the function list in which each of the potential target points participates, calculating the overlap degree of the function list in which the potential target point participates and the function list of the key disease module, and obtaining a functional consistency score;

[0018] Wherein, steps (3b) and (3c) are performed in sequence, interchangeably or simultaneously.

[0019] (4) For each of the potential target points, multiply the structural key score, the robustness score and the functional consistency score to obtain a comprehensive score; select the potential target points with high ranking / comprehensive score higher than a threshold value as therapeutic target points.

[0020] In another preferred embodiment, the disease to be tested is a complex disease.

[0021] In another preferred embodiment, the complex disease includes metabolic-associated steatohepatitis (MASH), cancer, diabetes, Alzheimer's disease.

[0022] In another preferred embodiment, before step (2a), it further includes the step of: preprocessing the transcriptome data.

[0023] In another preferred embodiment, the preprocessing includes: processing outliers, removing batch effects, and processing null values.

[0024] In another preferred embodiment, the data form of the transcriptome data includes: TPM, RKPM.

[0025] In another preferred embodiment, the correlation algorithm includes: Pearson correlation coefficient, Spearman correlation coefficient, Kendall correlation coefficient.

[0026] In another preferred embodiment, the soft threshold is 10.

[0027] In another preferred embodiment, the clustering algorithm includes: dynamic pruning tree algorithm, k-means clustering algorithm, hierarchical clustering algorithm.

[0028] In another preferred embodiment, the clustering algorithm is a dynamic pruning tree algorithm.

[0029] In another preferred embodiment, in step (2c), the correlation includes positive correlation and negative correlation.

[0030] In another preferred embodiment, in step (2c), the "correlation value is strong" means that the absolute value of the correlation value is >0.5, such as >0.6, >0.7, >0.8, >0.9.

[0031] In another preferred embodiment, the potential targets in the disease-specific protein-protein interaction network are identified by a network disassembly algorithm.

[0032] In another preferred embodiment, the network centrality algorithm comprises: degree centrality, betweenness centrality, closeness centrality.

[0033] In another preferred embodiment, the network disassembly algorithm comprises: GAPR.

[0034] In another preferred embodiment, the structural criticality score is calculated as follows:

[0035]

[0036] wherein, LCCi represents the LCC of each potential target, LCC original LCCoriginal represents the LCC of the original network, LCC new LCCnew represents the LCC of the new network after removing the potential target.

[0037] In another preferred embodiment, the machine learning algorithm comprises: a reinforcement learning algorithm, a deep learning algorithm, or a combination thereof.

[0038] In another preferred embodiment, the machine learning algorithm is a deep reinforcement learning algorithm.

[0039] In another preferred embodiment, the deep reinforcement learning algorithm comprises: FINDER.

[0040] In another preferred embodiment, the robustness score is calculated as follows:

[0041]

[0042] wherein, (n 1, n2, …, n N ) is the sequence of nodes in the disease-specific protein-protein interaction network G; wherein, C i is the ith connectivity in network G, which corresponds to the influence on network connectivity after removing the key node i

[0043] In another preferred embodiment, the functional enrichment analysis is KEGG pathway enrichment analysis.

[0044] In another preferred embodiment, the functional list comprises key pathway information.

[0045] In another preferred embodiment, the functional consistency score is calculated as follows:

[0046]

[0047] wherein, represents the functional list of key disease module n i , which comprises the top 20 key pathways, P ni represents the number of all the key pathways involved by the potential target n i .

[0048] In another preferred embodiment, the comprehensive score is calculated as follows:

[0049]

[0050] In another preferred embodiment, in step (4), the "top ranking" comprises top 5, top 10, top 50%.

[0051] In another preferred embodiment, in step (4), the top 5 potential targets are selected as therapeutic targets.

[0052] In another preferred embodiment, the method further comprises step: (5) applying validation, comprising steps:

[0053] (5a) for the therapeutic targets obtained in step (4), screening compounds interacting with the therapeutic targets from a compound library by network inference method;

[0054] (5b) experimentally validating the compounds, thereby confirming the compounds as lead compounds.

[0055] In another preferred embodiment, the network inference method comprises: wSDTNBI.

[0056] In another preferred embodiment, the experimental validation comprises: in vitro activity test, structural biology study, in vivo animal model pharmacodynamics validation, or a combination thereof.

[0057] In a second aspect of the present application, an electronic device is provided, comprising a processor and a memory, the memory having a plurality of executable instructions, the processor being configured to read the instructions and perform the steps in the identification method as described in the first aspect of the present application.

[0058] In a third aspect of the present application, a computer readable storage medium is provided, which stores computer executable instructions, which, when read and executed by a processor, implement the steps of the identification method according to the first aspect of the present application.

[0059] In a fourth aspect of the present application, a computer program product is provided, which comprises computer executable instructions, which, when executed by a processor, implement the steps of the identification method according to the first aspect of the present application.

[0060] It should be understood that, within the scope of the present application, each of the technical features described above and in the following (e.g. in the examples) can be combined with each other to form new or preferred technical solutions. Due to the limited space, they are not listed one by one here. BRIEF DESCRIPTION OF DRAWINGS

[0061] Figure 1 A schematic diagram of the ComplexDnet technical solution is shown. The diagram shows an overview of the overall process of the present application, which includes four main parts: (A) data collection, (B) key module identification, (C) key target identification, and (D) screening and experimental verification.

[0062] Figure 2 A WGCNA module analysis result diagram is shown. The diagram shows the gene modules identified in the MASH case and their correlation with the disease, as well as the gene correlation and pathway enrichment results of the key modules.

[0063] Figure 3 A network target identification and pathway enrichment analysis result diagram is shown. The diagram shows the robustness of the GAPR and FINDER algorithms, as well as the distribution of the four scoring methods, and the GO enrichment results of the target points finally screened by Disease-Score.

[0064] Figure 4 An immune cell infiltration and single cell analysis diagram is shown. The diagram shows the relationship between the key target points and the 28 immune cell types, as well as the expression distribution of the key target points at the single cell level, which is used to verify the biological function of the target points.

[0065] Figure 5 A RORγt inverse agonist screening and crystal structure diagram is shown. The diagram shows the in vitro activity data (IC50 value) of the compound PXT and its binding mode with the target point RORγt.

[0066] Figure 6 A PXT pharmacodynamic diagram in a BDL-induced liver fibrosis model is shown. The diagram shows experimental data of PXT in improving liver pathology and inhibiting fibrosis-related indicators in an animal model. DETAILED DESCRIPTION

[0067] The inventors, through extensive and in-depth tests and experiments, first proposed a network-based method for identifying disease targets, ComplexDnet. Based on the results of WGCNA, the method constructs a protein interaction network, and further uses the GAPR algorithm, the FINDER algorithm and biological pathway enrichment analysis to analyze the network structure, robustness and functional consistency, respectively, to obtain three scores. According to the above three scores, a comprehensive scoring algorithm is constructed, and the target points are sorted and screened according to the comprehensive score, so as to obtain the key treatment target. The method of the application is verified in 8 kinds of cancers, and the recall rate of target discovery is high, which significantly exceeds the traditional differential expression gene method. In addition, further verification is carried out on the liver data of MASH patients, and a novel treatment target RORγt is successfully discovered, and an effective drug targeting the RORγt target point is further determined. On this basis, the present application is completed.

[0068] It should be understood that the following description of specific methods and experimental conditions of the application in various details is provided to give a substantial understanding of the application. The definitions of certain terms used in this specification are provided below. Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which the application belongs.

[0069] The term

[0070] Where a numerical range is provided, unless the context clearly dictates otherwise, it is intended to include every delimited integer within the range, every sub-range between the recited upper and lower limits, and any other intervening range or point within the recited range, unless otherwise indicated. The upper and lower limits of these smaller ranges can independently be included in the smaller ranges, and are also encompassed, but only within the context of the larger range, unless otherwise indicated. For example, "1 to 50" includes "2 to 25," "5 to 20," "25 to 50," "1 to 10," etc.

[0071] As used herein, the term "containing" or "including" can be open, semi-closed and closed. In other words, the term also includes "consisting essentially of" or "consisting of."

[0072] As used herein, the term "and / or" relates to and encompasses any and all possible combinations of one or more of the associated listed items.

[0073] As used herein, the term "significant" refers to the observed effect (e.g., the difference between the experimental group and the control group) in a hypothesis test is unlikely to have arisen from random error alone. A hypothesis test includes: a null hypothesis (H0), which assumes that the observed effect does not exist (e.g., there is no difference between the experimental group and the control group); a p-value, which is the probability of observing the current or more extreme effect when H0 is true; and a significance threshold (a). The hypothesis test is usually judged to be significant by the significance threshold, which is typically 0.05. If the p-value is < a, then H0 is rejected, i.e., the observed effect exists, and the result is said to be "significant".

[0074] As used herein, the term "robustness" refers to the ability of a system to maintain stable core functions under internal parameter variations or external perturbations.

[0075] As used herein, the "RORyt-GAL4 luciferase reporter assay", "HTRF", and other assays employed are routine compound validation methods, and "bile duct ligation", "CCl4-induced mouse model of liver fibrosis", and other assays employed are routine in vivo validation methods, all of which are well known to those skilled in the art.

[0076] WGCNA

[0077] Provided herein is a method for identifying disease targets, which further combines network centrality algorithm and / or network disassembly algorithm, machine learning algorithm, and functional analysis on the basis of WGCNA, thereby forming ComplexDnet.

[0078] As used herein, the term "weighted gene co-expression network analysis" and "WGCNA" can be used interchangeably, which is a systems biology method for constructing co-expression networks based on gene expression data, narrowing down the range of potential targets by identifying gene modules highly correlated with disease phenotypes.

[0079] As used herein, the term "dynamic pruning tree algorithm" is a method of cluster analysis, which can be dynamically cut according to the set parameters (e.g., the minimum number of module characteristic genes), thereby obtaining multiple modules, and merging similar modules.

[0080] As used herein, the term "minimum module characteristic gene" is defined as the first principal component of a given module, which can be regarded as a representative of the gene expression profile in the module.

[0081] The commonly used methods and procedures of weighted gene co-expression network analysis are known to those skilled in the art, which can be analyzed by WGCNA package in R language. The procedure of WGCNA usually includes: quality control and preprocessing of input data; constructing co-expression network according to the transcriptome data; dividing the co-expression network into multiple gene modules according to the clustering algorithm; calculating the correlation of each gene module with the phenotype data. Finally, the correlation degree of different gene modules with the phenotype can be obtained.

[0082] The method for quality control and preprocessing of input data is known to those skilled in the art, which includes processing outliers, removing batch effects, and processing null values. For outliers, box plots, cluster plots, or principal component analysis can be used to identify outliers deviating from the overall level. Generally, outliers are deleted before analysis, or can be retained according to the analysis situation. For batch effects, various methods can be used to remove batch effects, such as ComBat method, proxy variable method, ratio-based method, average center method, and distance weighted discrimination method. For null values, deletion and imputation are usually used for processing, wherein the imputation methods include special value filling, average value filling, hot card filling, K-means filling, using regression algorithm filling, using expectation maximization method filling, and multiple imputation.

[0083] The finally obtained correlation degree can be positive correlation or negative correlation.

[0084] The identification method of the present application

[0085] The present application provides an identification method of disease target based on network, comprising the steps of:

[0086] (1) providing transcriptome data, protein-protein interaction network data and functional annotation data of a disease to be tested; the transcriptome data includes gene expression data and phenotype data of the disease.

[0087] (2) performing weighted gene co-expression network analysis, comprising the steps of:

[0088] (2a) calculating the correlation between each pair of genes according to the transcriptome data using correlation algorithm to obtain correlation value; taking each gene as the node of the network, and taking the value obtained by performing soft thresholding operation on the correlation value as the edge of the network to construct the co-expression network; preferably, the soft threshold is 10;

[0089] (2b) dividing the co-expression network into multiple gene modules by using clustering algorithm; preferably, the clustering algorithm is dynamic pruning tree algorithm; other gene clustering or co-expression analysis algorithms such as k-means clustering and hierarchical clustering can also be used to identify gene modules;

[0090] (2c) calculating the correlation between the module eigengene of each of the gene modules and the phenotype data using a correlation algorithm to obtain a correlation value; and selecting a module with a strong correlation value as a key disease module.

[0091] Preferably, the WGCNA analysis package is selected for step (2).

[0092] (3) identification and evaluation of key targets, comprising steps of:

[0093] (3a) mapping all genes in the key disease module obtained in step (2) to the protein-protein interaction network data to extract a sub-network composed of the genes and their interaction relationships to obtain a disease-specific protein-protein interaction network.

[0094] (3b) identifying potential targets in the disease-specific protein-protein interaction network using a network centrality algorithm and / or a network disassembly algorithm, and calculating the change in the largest connected component (LCC) after removing the potential targets to obtain a structural key score. Preferably, the network disassembly algorithm is used to identify the potential targets in the disease-specific protein-protein interaction network. Preferably, the network disassembly algorithm comprises GAPR. Other network centrality or network disassembly algorithms can also be used, for example, the Betweenness Centrality can be used to identify the "bridge" nodes in the network instead of GAPR. Other network control theory or key node identification algorithms can be used instead of FINDER.

[0095] (3c) removing several nodes in the disease-specific protein-protein interaction network using a machine learning algorithm to minimize the cumulative normalized connectivity (R), and then gradually introducing the removed nodes to calculate the change in the cumulative normalized connectivity and normalize the change to obtain a robustness score.

[0096] (3d) performing functional enrichment analysis on the key disease module to obtain a function list of the key disease module; querying the function list involved in each of the potential targets obtained in step (3b) to calculate the overlap degree between the function list involved in the potential target and the function list of the key disease module to obtain a functional consistency score.

[0097] Steps (3b) and (3c) can be performed sequentially, interchangeably or simultaneously.

[0098] (4) For each of the potential target points, multiply the structure criticality score, the robustness score and the functional consistency score to obtain a comprehensive score; select the potential target points with a high ranking / the comprehensive score higher than a threshold value as therapeutic target points.

[0099] The method further comprises step (5) applying verification, comprising:

[0100] (5a) For the therapeutic target points obtained in step (4), screen compounds interacting with the therapeutic target points from a compound library by using a network inference method.

[0101] (5b) Experimentally verify the compounds, thereby confirming the compounds as lead compounds.

[0102] In step (5), other mature drug-target interaction prediction methods such as molecular docking, deep learning-based DTI prediction models, etc. can be used to screen compounds binding to the target points.

[0103] These alternatives can replace the original modules in specific scenarios, but the core innovation of the present application lies in its integrated design idea and integrated process framework in three dimensions of structure, robustness and function.

[0104] The main advantages of the present application include:

[0105] (1) The method of the present application is more accurate and more reliable in prediction, and can identify target points that are not necessarily the highest in connectivity but are crucial to the structure and function stability of the disease network. The recall rate of target point discovery in the present application in 8 kinds of cancers reaches an average of 77.63%, which is about 20% higher than that of traditional DEG analysis, and is also superior to methods relying only on centrality, effectively reducing the false positive rate.

[0106] (2) The method of the present application first develops a target comprehensive score algorithm Disease-Score, which quantitatively integrates the structural importance, robustness influence and functional correlation of the target points, so that the basis for ranking the target points is more scientific and the physical meaning is more explicit. This not only can screen target points, but also can reveal the important reasons of the target points in a specific disease from a system level, providing a clearer mechanism hypothesis for subsequent research.

[0107] (3) The present application provides a complete integrated process, which only needs to input omics data to obtain a target point list evaluated and ranked in multiple dimensions, and can directly link to the next step of virtual screening, greatly shortening the discovery period from data to candidate drugs, and has accessibility and application efficiency.

[0108] The application will be further described in conjunction with specific examples. It should be understood that these examples are only used to illustrate but not limit the scope of the application. The experimental methods in the following examples, if not otherwise specified, are generally carried out according to the conventional conditions, for example, the conditions described in Sambrook et al., Molecular Cloning: A Laboratory Manual (New York: Cold Spring Harbor Laboratory Press, 1989), or the conditions recommended by the manufacturer. Unless otherwise specified, percentages and parts are weight percentages and weight parts.

[0109] Materials and methods:

[0110] 1. Data collection:

[0111] Transcriptomic data (FPKM format) of 8 cancer types were obtained from GDC PANCAN cohort. Quality control was performed on the data, including removing samples with incomplete clinical annotation or ambiguous diagnosis, resulting in 955 sample data. MASH transcriptomic data (GSE126848) was obtained from NCBI GEO database, including 57 liver biopsies (31 MASH patients and 26 age-matched healthy controls). Liver scRNA-seq data of hepatocytes under high sucrose / fat diet (HSD) condition (GSE182365) was obtained from GEO database to study diet-induced transcriptional changes. PPI network data was obtained from publicly available human interactome resource. Cancer driver genes (IntOGen database) were mapped to corresponding cell lines using standardized nomenclature in Cancer Cell Line Encyclopedia (CCLE). Clock genes (n=2,510) were identified by GeneCards (v4.14) with stringent relevance score (>7) to ensure specificity of biological pathways.

[0112] 2. WGCNA analysis:

[0113] WGCNA package in R was used for analysis. The basic process is as follows: (1) Calculate the normalized TPM or FPKM values of transcriptomic data, then perform duplicate data deletion and averaging to ensure the robustness of network analysis. (2) Screen and select appropriate soft threshold to optimize network construction, thereby improving specificity and sensitivity. (3) Use dynamic cutting tree method to cluster and identify key modules, then calculate the correlation between modules, diseases and individual genes. Screen out the modules with the highest correlation comprehensively, and determine the key genes in these important modules.

[0114] 3. Construction of protein interaction network to identify key proteins:

[0115] Disease-specific PPI networks were drawn based on PPI network data. Then two network-based methods were used to determine key targets: GAPR and FINDER. The articulation points (APs) identified in the network by GAPR were defined as key nodes, and removing the node would disrupt the network connection or increase the number of connected components. First, all articulation points in the network were identified, and then they were systematically removed in an iterative peeling process until there were no articulation points. Each articulation point identified in each layer was considered a key protein in the PPI network. FINDER is a deep reinforcement learning-based method that first trains on small synthetic graphs tailored for specific applications, and then applies it to larger real-world networks, effectively identifying key nodes and dismantling network structures.

[0116] 4. Single-sample gene set enrichment analysis (ssGSEA):

[0117] ssGSEA is a commonly used method for analyzing immune cell infiltration, and the GSVA package in R was used to validate the key targets determined by the above method. This method can assess changes in pathways and biological processes in the sample population, providing insight into the heterogeneity of the immune response in patients with monocytosis. Specifically, ssGSEA was used to assess the infiltration levels of 28 different immune cell types and determine the changes in key genes between the MASH group and the control group. Wilcoxon rank-sum test was used to compare the results between groups. In addition, Spearman correlation analysis was performed to explore the relationship between different immune cell types, providing a comprehensive understanding of the interactions between immune cells.

[0118] 5. Single-cell RNA-seq data analysis:

[0119] The single-cell transcriptome dataset (GSE182365) included hepatocyte and hepatic stellate cell data from high-fat and high-sucrose diets. These data were pre-processed using the Seurat R package (v4.3.0). Key indicators such as the number of molecules (nCount RNA) and the number of genes detected per cell (nFeature RNA) were compared with sequencing reads to ensure data integrity. To address the issue of mitochondrial genome contamination (commonly found in low-quality cells or dead cells), quality control was adopted to filter cells with 200 < nFeature_RNA < 2500 and mitochondrial count (%.mt) < 5%.

[0120] Cell clustering was performed using PCA, and then UMAP was used to show clear cell clustering results. By comparing the marker genes in the DISCO database and the specific genes of each cell type, statistically significant cell marker genes (corrected p value less than 0.05) were determined and used for the classification of cell groups, thereby gaining a deep understanding of the distribution and abundance of key genes in various cell subpopulations.

[0121] 6. Network-based RORyt ligand virtual screening:

[0122] RORyt-related ligands were identified by wSDTNBI according to the substructure-drug-target network. DTI data were collected from published literature. Subsequently, the optimal wSDTNBI model Global-wSDTNBI-FCFP_4 was selected for network-based virtual screening (VS) of RORyt (UniProtion number: P51449, gene symbol: RORC). First, a chemical library consisting of 198 triterpenoids was obtained from Topscience. After generating FCFP_4 fingerprints for compounds outside the global DTI network, the newly obtained compound-substructure associations were integrated into the model, enhancing the model's ability to predict the potential targets of these compounds. Then, the model was used to predict 20 new targets for drugs within and outside the global DTI network. Compounds with RORyt as a newly predicted target were designated as potential RORyt ligands. Finally, considering the availability of raw materials and maximizing the demand for triterpenoids, compounds that appeared more than twice were retained. In addition, raw materials with potential hepatoprotective effects were also studied. Therefore, 28 compounds were screened from 7 plants for further mechanism and activity studies. Subsequently, the 28 predicted RORyt ligands were purchased from Topscience for experimental validation.

[0123] Example 1: Construction of a network-based process for identifying disease targets ComplexDnet

[0124] The present invention proposes a network-based strategy named ComplexDnet, whose technical solution comprises the following steps: Figure 1 ):

[0125] S1: Data collection and preprocessing Figure 1 A). This process requires obtaining multi-source data related to the target complex disease. Mainly including:

[0126] (1) Transcriptome data: Obtain gene expression profile data (such as FPKM or TPM format) of disease samples and control samples from public databases (such as GEO) or self-experiments.

[0127] (2) Protein-protein interaction (PPI) network data: Obtain a comprehensive human PPI background network from public databases (e.g., Human Interactome).

[0128] (3) Functional annotation data: Obtain gene functional annotation information such as KEGG pathways.

[0129] S2: Identification of key disease modules Figure 1 B). This step uses the weighted gene co-expression network analysis (WGCNA) method, and the specific steps are as follows:

[0130] (1) Construction of co-expression network: Preprocess the transcriptome data, calculate the expression correlation between genes, and select an appropriate soft threshold (soft threshold, referred to as power in this step) to construct a scale-free gene co-expression network.

[0131] Specifically, the input data is the transcriptome data obtained in S1, which is usually an expression matrix, with genes as rows and samples as columns. Quality control can be performed according to the quality of the transcriptome data, such as checking for outliers by hierarchical clustering, removing batch effects, removing null values, etc.

[0132] Calculate the correlation (such as Pearson correlation coefficient) between gene pairs (such as genex and geney) to obtain a correlation value, where each gene constitutes a point in the co-expression network; based on the soft threshold, the correlation value is processed to form the edge of the co-expression network. The type of network constructed is different, and the calculation method of the edge attribute is also different. For an undirected network, the calculation method of the edge attribute is abs(cor(genex, geney)) power , where abs is the absolute value, and cor(genex, geney) is the correlation value of the gene pair; for a directed network, the calculation method of the edge attribute is (1+cor(genex, geney) / 2) power . Usually, an undirected network is used to construct a co-expression network.

[0133] The selection of the soft threshold power follows the principle of making the constructed network more consistent with the characteristics of a scale-free network, i.e., making the structure R 2 of the network spectrum reach 0.8, and more preferably reach 0.85. When the present application selects power = 10, it can reach 0.85.

[0134] A preferred weighted co-expression network is constructed according to the above method.

[0135] (2) Module identification: clustering algorithm is used to divide the co-expression network into multiple gene modules (modules); each module is composed of a group of genes with highly coordinated expression patterns. Among them, the clustering algorithm includes dynamic tree-cutting algorithm, k-means clustering algorithm, hierarchical clustering algorithm, etc. The dynamic tree-cutting algorithm is selected by the present application for module identification.

[0136] (3) Screening of key modules: the correlation between the module eigengene (ME) of each module and the clinical phenotype of the disease (such as MASH state) is calculated. The module with the strongest correlation and the most statistically significant is selected as the key disease module for subsequent analysis.

[0137] Specifically, a numerical feature matrix related to the expression matrix needs to be provided. For disease traits, only the state corresponding to the sample in the expression matrix needs to be provided, such as:

[0138] TCGA-AU-3779-01A Disease

[0139] TCGA-A6-5662-11A Health.,

[0140] Wherein the first column is the sample and the second column is the state.

[0141] The module can be screened by functional analysis. The functional analysis includes gene ontology (GO) analysis, KEGG enrichment analysis, etc.

[0142] S3: Identification and multidimensional evaluation of key targets. This step is to construct a module-specific PPI network and apply multiple algorithms for comprehensive evaluation.

[0143] S3a: Construction of disease-specific PPI network. All genes in the key module screened in S2 are mapped to the background PPI network, and a subnetwork composed of these genes and their interactions is extracted to obtain a "disease-specific PPI network".

[0144] S3b: Network structure key analysis. The greedy articulation points removal (GAPR) algorithm is applied to iteratively identify and remove articulation points in the network. Articulation points are nodes that, when removed, will cause the network to split into more connected components or be directly disconnected.

[0145] For each gene (node ​​ni) identified as a potential target, calculate the change in the largest connected component (LCC) resulting from its removal, denoted as the LCC change value. ni This value quantifies the structural importance of a target to network connectivity. The calculation method is as follows:

[0146]

[0147] Among them, LCC ni LCC represents the LCC of each potential target gene. original LCC represents the original network. new This indicates the LCC of the new network after the node is removed.

[0148] S3c: Network robustness analysis. The FINDER (Finding key players in complex networks) algorithm, based on deep reinforcement learning, is applied. This algorithm aims to find an optimal sequence of node removals to dismantle the entire network as quickly as possible.

[0149] The algorithm calculates a robustness score for each node in the network. ni This score reflects the importance of the node in maintaining the overall robustness of the network. A higher score indicates that removing the node would be more destructive to the network. Specifically, for a protein-protein interaction network G=(n,E) (including the node set n... N Given the edge set E and a predefined connectivity index θ, the FINDER method removes the node sequence (n1, n2, ..., nn). N This minimizes the cumulative normalized connectivity (R). To assess node importance, the change in R is calculated by iteratively reintroducing the removed nodes, and the results are normalized to obtain the node's robustness. ni .

[0150] The calculation method is as follows:

[0151] Cumulative normalized connectivity R:

[0152]

[0153] in, For the predefined connectivity measure σ, similar to the LCC above (G\{n1,n2,…,n N} represents the connectivity of the network remaining after removing nodes from the set in sequence from the network G. Among them, Ci For the ith connectivity in the network G, it corresponds to the impact on the network connectivity after the deletion of the key node i.

[0154] S3d: Functional consistency analysis. KEGG pathway enrichment analysis was performed on the key disease module to obtain a list of pathways most relevant to the core mechanism of the disease. For each potential target, all KEGG pathways in which it participates were queried. At the same time, the degree of overlap between the pathways in which the target participates and the list of core pathways of the module was calculated to obtain a pathway similarity score (Similarity ni ). This score quantifies the degree to which the target matches the core pathology of the disease in terms of function. The calculation method is as follows:

[0155]

[0156] wherein, denotes the functional list of the key disease module n i , which includes the top 20 key pathways, P ni denotes the number of all the key pathways in which the potential target n i participates.

[0157] S4: Comprehensive ranking of key targets.

[0158] S4a: Calculation of comprehensive score (Disease-Score). The present application innovatively multiplies the evaluation indicators of the above three dimensions to obtain the final score of each potential target. The calculation method is as follows:

[0159]

[0160] S4b: Ranking of targets. According to the Disease-Score from high to low, all potential targets are ranked, and the target with the highest score is the most critical therapeutic target predicted by the present application.

[0161] In order to evaluate the universality across diseases, two indicators were used to quantify the stability of the prediction:

[0162] (1) Recall rate of top-ranked candidate targets (Recall j10 denotes the recall rate of the top 10 candidate targets of disease j; Rank 10 denotes the top 10 candidate targets; T j denotes the driver gene of disease j):

[0163]

[0164] (2) Recall rate of total candidate targets (Recall j50Ranking top 50% of candidate targets for disease j; Rank 50 Ranking top 50% of candidate targets; T j Ranking top 50% of candidate targets for disease j; Rank

[0165]

[0166] S5: Downstream applications and validation.

[0167] S5a: Virtual screening: for the top ranked target (e.g. RORyt), use a network inference method (e.g. wSDTNBI) to virtually screen a compound library (e.g. natural product library) to predict active compounds that can interact with the target.

[0168] S5b: Experimental validation: test the virtually screened compounds for in vitro activity (e.g. HTRF, fluorescence quenching), structural biology studies (e.g. crystal structure analysis) and pharmacodynamics validation in animal models (e.g. BDL or CCl4-induced liver fibrosis model) to finally confirm active lead compounds (e.g. Panaxatriol PXT).

[0169] Example 2: Performance evaluation of ComplexDnet

[0170] To demonstrate the accuracy and robustness of ComplexDnet, the method was applied in 8 cancer datasets to observe the success rate of discovered targets. With the IntOGen-DriverGenes database as the benchmark, different performance patterns were discovered: In bladder cancer (BLCA), cholangiocarcinoma (CHOL), glioblastoma (GBM) and hepatocellular carcinoma (LIHC), ComplexDnet showed statistically higher driver gene priority compared to differential expression gene (DEG) analysis (P < 0.01, Mann-Whitney U test). The performance of esophageal adenocarcinoma (ESCA) and colorectal adenocarcinoma (COAD) was comparable. Although DEG showed nominal ranking advantage in breast cancer (BRCA) and lung adenocarcinoma (LUAD), key analysis found that ComplexDnet discovered more effective driver genes (BRCA: 33 vs 17; LUAD: 16 vs 3) despite lower position ranking.

[0171] To objectively evaluate the biological relevance, precision-sensitive indicators (Recall 50 , Recall 10 ) were adopted. In all cancer types, the Recall 50 of ComplexDnet was ≥ 60%, with the highest Recall 50The highest, reaching 87.09%. It is worth noting that there are four malignant tumors of Recall 50 The critical value exceeds 80%. This method shows very high precision (> 50%) in BRCA, GBM and LIHC, and more than half of the top 10 predicted genes are cataloged driver genes. Most importantly, the system benchmark test confirms that ComplexDnet improves the performance of 10-40% compared to traditional DEG methods, thus establishing its ability to solve the complexity of multi-genome signaling in pan-cancer target identification while minimizing false discovery rates.

[0172] Example 3: Analysis results of liver transcriptome data of MASH patients

[0173] Liver transcriptome data (GSE126848) of MASH patients and healthy controls were used. Through WGCNA analysis, two key modules highly negatively correlated with MASH pathological state were successfully identified: the "greenyellow" module (correlation coefficient -0.81) and the "blue" module (correlation coefficient -0.76). KEGG enrichment analysis showed that both modules were significantly enriched in the "Th17 cell differentiation" pathway, revealing the core role of this pathway in MASH.

[0174] Disease-specific PPI networks were constructed for these two modules. GAPR and FINDER algorithms were applied for analysis, and combined with pathway similarity calculation, the top five key targets were finally selected by Disease-Score scoring, which were A2M, CD274, HCK, PTPRD and RORyt. Among them, RORyt is a key regulatory factor for Th17 cell differentiation, but its direct role in MASH has been reported less, which is a novel potential target.

[0175] Example 4: Target verification

[0176] Immune infiltration analysis (ssGSEA) showed that the expression of RORyt was significantly positively correlated with the infiltration level of Th17 cells in the liver of MASH patients (r = 0.48), which verified the relationship between the target and the core pathway at the level of clinical samples.

[0177] Single-cell RNA sequencing (scRNA-seq) analysis further confirmed that the expression of RORyt was significantly in the T cell subpopulation of MASH liver.

[0178] Example 5: Active compound screening and verification

[0179] Using wSDTNBI method to screen a natural product library, 89 potential RORγt ligands were predicted, and finally 28 compounds were screened for experimental verification.

[0180] First, in vitro verification experiments were carried out. In the RORγt-GAL4 luciferase reporter gene experiment, Panaxatriol (PXT) was found to have very strong inhibitory activity, with an IC50 value of 0.01 μM.

[0181] In the HTRF coactivator recruitment experiment, the IC50 value of PXT was 0.11 μM.

[0182] Through X-ray crystal diffraction, the crystal structure of the complex of PXT and RORγt-LBD (PDB ID: 7E3M) was successfully resolved, which confirmed that PXT directly binds and inhibits RORγt as an inverse agonist at the atomic level.

[0183] Second, in vivo pharmacodynamic experiments were carried out. In the bile duct ligation (BDL) and carbon tetrachloride (CCl4) induced mouse liver fibrosis model, PXT treatment significantly reduced liver damage, inflammation and collagen deposition, and its effect was comparable or superior to that of the positive drug obeticholic acid (OCA).

[0184] PXT significantly down-regulated the expression of IL-17 and TGF-β in the liver, and verified its anti-fibrosis effect by inhibiting the RORγt-IL-17 axis from the molecular mechanism.

[0185] This example fully demonstrates the whole process of the present application from data analysis to target discovery, to active compound confirmation and in vivo verification, which proves the feasibility, advancement and practical value of the technical scheme of the present application.

[0186] All the documents mentioned in the present application are cited as references in the present application, just as each document is cited as a reference. In addition, it should be understood that those skilled in the art can make various modifications or modifications to the present application after reading the above teachings of the present application, and these equivalent forms also fall within the scope defined by the claims attached to the present application.

Claims

1. A network-based method for identifying disease targets, characterized in that, The method comprises the steps of: (1) providing transcriptome data, protein-protein interaction network data and functional annotation data of a disease to be tested; the transcriptome data comprises gene expression data and phenotype data of the disease; (2) performing weighted gene co-expression network analysis, comprising the steps of: (2a) calculating the correlation between each pair of genes according to the transcriptome data by using a correlation algorithm to obtain a correlation value; taking each gene as a node of a network and taking a value obtained by performing soft thresholding on the correlation value as an edge of the network to construct a co-expression network; (2b) dividing the co-expression network into a plurality of gene modules by using a clustering algorithm; (2c) calculating the correlation between a module characteristic gene in each gene module and the phenotype data by using a correlation algorithm to obtain a correlation value; and selecting a module with a strong correlation value as a key disease module; (3) identifying and evaluating key targets, comprising the steps of: (3a) mapping all genes in the key disease module obtained in step (2) to the protein-protein interaction network data to extract a subnetwork composed of the genes and their interaction relationships to obtain a disease-specific protein-protein interaction network; (3b) identifying potential targets in the disease-specific protein-protein interaction network by using a network centrality algorithm and / or a network disassembly algorithm, calculating the change in the largest connected component (LCC) after removing the potential targets to obtain a structural key score; (3c) removing a plurality of nodes in the disease-specific protein-protein interaction network by using a machine learning algorithm to minimize the cumulative normalized connectivity (R), and then gradually introducing the removed nodes to calculate the change in the cumulative normalized connectivity and normalizing the change to obtain a robustness score; (3d) performing functional enrichment analysis on the key disease module to obtain a functional list of the key disease module; querying the functional list in which each potential target participates, calculating the overlap degree between the functional list in which the potential target participates and the functional list of the key disease module to obtain a functional consistency score; wherein steps (3b) and (3c) are performed in sequence, interchangeably or simultaneously; (4) for each potential target, multiplying the structural key score, the robustness score and the functional consistency score to obtain a comprehensive score; and selecting a potential target with a high comprehensive score and / or a high ranking as a therapeutic target.

2. The method of authentication of claim 1, wherein, The network disassembly algorithm comprises GAPR.

3. The identification method as described in claim 1, characterized in that, The structural criticality score is calculated as follows: wherein, LCC represents the LCC of each potential target, LCC original LCC represents the LCC of the original network, LCC new LCC represents the LCC of the new network after removal of the potential target.

4. The method of claim 1, wherein the step of identifying comprises the step of: 5 determining the presence of a mutation in the nucleic acid sequence of the patient. The deep reinforcement learning algorithm comprises FINDER.

5. The method of authentication of claim 1, wherein, The robustness score is calculated as follows: wherein (n1, n2,..., n N ) is a sequence of nodes in the disease-specific protein-protein interaction network, is a predetermined connectivity.

6. The method of authentication of claim 1, wherein, the functional consistency score is calculated as follows: wherein, represents the key disease module n i a list of functions, including the top 20 key pathways, P ni represents the number of all the key pathways involved in the potential target n i .

7. The method of authentication of claim 1, wherein, The comprehensive score is calculated as follows:

8. An electronic device comprising a processor and a memory, characterized in that The memory has a plurality of executable instructions, and the processor is configured to read the instructions and perform the steps in the identification method of claim 1.

9. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores computer-executable instructions, and the computer-executable instructions are read and executed by the processor to implement the steps in the identification method of claim 1.

10. A computer program product comprising computer executable instructions, characterised in that, The computer-executable instructions are executed by the processor to implement the steps in the identification method of claim 1.

Images