Multilayer regulation and control network key node identification method based on global interactive transitivity and feature fusion

Abstract

The invention provides a multi-layer regulatory network key node identification method based on global interactive transitivity and feature fusion, and belongs to the field of bioinformatics and computational biology, and the method comprises the following steps: constructing an lncRNA-miRNA-mRNA-TF ceRNA regulatory network of a tumor; selecting a classic network topology attribute and a biological attribute (expressing a variable coefficient (Corr) and phenotypic correlation (Coef)) as node features; node features are fused by using a Bayes algorithm, a Bayes value obtained by each node is used as an importance weight of the node, a PageRank algorithm is used for transmission under the interaction of a global network, and a BPR importance index of a plurality of feature combinations is further constructed; a tumor-related node in the RNADisease database is adopted as a positive sample to verify the prediction effect of the BPR index; comparing the prediction capabilities of the single-feature index, the multi-feature index and the optimal BPR index in RF and Bayes; a result shows that the optimal BPR index (including Degree, EPC, BottleNeck and Coef) has the strongest capability in key node identification. Taking lung adenocarcinoma as an example, the key nodes related to lung adenocarcinoma progress are successfully identified, and the potential important value of the method in biomedical research and clinical application is shown.

Description

Technical field:

[0001] The present invention relates to the fields of bioinformatics and computational biology, and in particular to a method for identifying key nodes in a multi-level regulatory network by combining network topological features, biological features and global interactive transitivity. The method is optimized and implemented based on a Bayes-PageRank algorithm. Background technology:

[0002] The occurrence and development of biological processes depend on gene regulatory networks, which are scale-free networks in which a few nodes play a core role. Identifying these key nodes is of great significance for a deep understanding of disease mechanisms, building precise disease models, and discovering drug targets. Especially in the study of network disorder diseases such as cancer, identifying key molecules is crucial for early screening, diagnosis, and treatment of diseases.

[0003] Traditional key node identification methods are mostly focused on single-level networks (such as protein-protein interaction networks), while multi-level gene regulatory networks (such as gene-miRNA-lncRNA-TF) can more accurately reflect the complex process of disease occurrence and development. Traditional methods usually rely on the topological characteristics of the network to identify key nodes, but often ignore the influence of other biological properties. Network topological characteristics mainly evaluate the importance of nodes based on their neighbor information, such as degree, and judge their importance by the number of nodes' neighbors. Other methods focus on the shortest path between nodes, including betweenness, closeness, stress, bottleneck, and radiality. In recent years, some methods based on a combination of multiple topological attributes have also been proposed, such as the TOPSIS method, which evaluates the importance of nodes by combining multiple centrality indicators. At the same time, some researchers have proposed identifying key nodes in the network from a multi-level perspective. However, these methods generally ignore the importance of global interaction transitivity and low-connectivity nodes. Existing studies have shown that there are low-connectivity key nodes in multi-level regulatory networks that independently regulate hub genes. Some researchers have proposed the edge permeability component (EPC) method, which evaluates the connectivity of nodes by randomly deleting a small number of edges and found that some low-connectivity nodes are crucial to the integrity of the network structure. In addition, methods based on network node clustering, such as the IPCA algorithm, have been used to analyze key genes in gastric cancer, but these methods mainly rely on network topological properties and ignore the biological characteristics of nodes.

[0004] Based on this, we proposed a comprehensive importance index, Bayes-PageRank (BPR), which can be applied to multi-level lncRNA-miRNA-mRNA-TF ceRNA regulatory networks. As an effective tool for evaluating the importance of network nodes, this method can more comprehensively evaluate the importance of nodes in the network by combining the topological structure of the network, the biological characteristics of the nodes (such as the coefficient of variation of gene expression and the correlation between gene expression and phenotype), and global interaction transitivity. Summary of the invention:

[0005] The purpose of the present invention is to overcome the shortcomings of the existing background technology and provide a method for identifying key nodes in a multi-layer regulatory network based on global interactive transitivity and feature fusion. The method comprises the following steps:

[0006] Step 1: From TCGA (The Cancer Genome Atlas Program, https: / / portal.gdc.cancer.gov / ) database to download cancer RNA-seq (including mRNA, lncRNA and TF) and mature miRNA Isoform-Seq gene expression data and clinical sample data, and perform logCPM normalization on the expression data of all nodes.

[0007] Step 2: The R software package "Deseq2" was used to perform differential analysis of mRNA, miRNA, and lncRNA between adjacent normal samples and stage I cancer samples, and differential TFs were screened from the differential mRNAs.

[0008] Step 3: Construct miRNA-mRNA regulatory relationship pairs.

[0009] Step 4: Construct miRNA-lncRNA regulatory relationship pairs.

[0010] Step 5: Construct mRNA-TF interaction pairs.

[0011] Step 6: Use CytoscapeV 3.7.1 to construct the lncRNA-miRNA-mRNA-TF ceRNA regulatory network of four types of nodes and visualize it.

[0012] Step 7: Calculation of eigenvalues ​​and feature selection of network nodes. Seven topological features including degree, edge permeability component, bottleneck centrality, closeness centrality, radial degree centrality, betweenness centrality and stress centrality of nodes in the network were selected for calculation, as well as two biological features including coefficient of variation and phenotypic correlation of nodes. The above features were screened for subsequent analysis.

[0013] Step 8: Construct the Bayes-PageRank (BPR) importance index. The Bayesian network structure is learned from the lncRNA-miRNA-mRNA-TF regulatory network constructed above by the score search algorithm. The features screened in step 8 are placed individually or randomly combined into the unsupervised Bayesian model for parameter learning. The conditional probability value of the node is obtained by using the maximum likelihood estimation for parameter estimation and the variable elimination algorithm for inference. It is input into the PageRank algorithm as a "weight". Various BPR indicators with single features or combined features are thus constructed.

[0014] Step 9: Obtain the optimal BPR importance index. http: / / www.rnadisease.org / ) Screen miRNA, mRNA (including TF) and lncRNA molecules related to cancer and marked them as positive samples; mark other molecules as negative samples. Calculate the AUC value of each BPR indicator based on the ROC curve. Evaluate the predictive ability of all BPR indicators for positive samples in step 8, compare the AUC values ​​of different indicators, and select the optimal BPR indicator (BPR-best).

[0015] Step 10: Evaluate the correlation between the optimal BPR index and the key nodes of the disease. Based on the optimal BPR model, calculate the BPR value of each node in the ceRNA regulatory network and analyze the proportion of positive samples in different BPR value ranges, that is, power law distribution.

[0016] Step 11: Comparison of the prediction ability of the optimal BPR indicator and other single-feature or multi-feature indicators in the Bayesian model and the random forest model.

[0017] Step 12: Key node screening and literature verification: Select the top 5% nodes with the best BPR index as key nodes, and use literature to verify their relevance to the occurrence and development of cancer.

[0018] Furthermore, when performing differential analysis on mRNA, miRNA and lncRNA in step 2, the threshold |log2(FC)|≥1, adj.p≤0.05. According to the experimental database Hocomocov11( https: / / hocomoco12.autosome.org ) and Animal TFDB( http: / / bioinfo.life.hust.edu.cn / AnimalTFDB4 / ), manually curated human TFs, and screened differential TFs from differential mRNAs.

[0019] Furthermore, in step 3, four biological experiment databases mir2disease ( http: / / www.mir2disease.org / )、miRecords( http: / / mirecords.biolead.org / )、TarBase( http: / / www.microrna.gr / tarbase )、miR2Disease( http: / / mirtarbase.mbc.nctu.edu.tw / php / index.php ) and a prediction database Starbase (starbase.sysu.edu.cn), and the interaction relationship pairs obtained in the above databases were crossed to obtain miRNA-mRNA regulatory pairs.

[0020] Furthermore, in step 4, the biological experiment database LncBase ( https: / / diana.e- ce.uth.gr / lncbasev3 / interactions ) and the interaction relationship pairs obtained from the prediction database StarBase (starbase.sysu.edu.cn) were crossed to obtain the miRNA-lncRNA regulatory factor pairs.

[0021] Furthermore, in step 5, the biological experiment database TRRUST ( www.grnpedia.org / trrust ) and the prediction database hTFtarget( http: / / bioinfo.life.hust.edu.cn / hTFtarget ) The two databases were crossed to obtain mRNA-TF interaction pairs.

[0022] Furthermore, in step 6, in order to reduce the false positives of the obtained regulatory relationship pairs, we used the Pearson correlation coefficient to calculate the correlation between the above relationship pairs and the gene expression in the paired samples. Since TFs with weak binding ability show stronger synergistic effects, different thresholds were selected for the above relationship pairs. In miRNA-mRNA and miRNA-lncRNA, all negative correlation pairs with p-value <0.05 (negative regulation of lncRNA and mRNA by miRNA) were retained. In TFs-mRNA, interaction pairs with |r|≥0.2 and p-value <0.05 were retained. In addition, the phenotypic correlation (0 or 1) was calculated, and nodes with correlation r≥0.1 and p-value <0.05 were retained.

[0023] Furthermore, the step 7 includes:

[0024] (1) Degree (Deg):

[0025] Deg(v)=|N(v)|

[0026] The higher the degree of a node, the more nodes it is connected to, which means it is more important.

[0027] (2) Edge Percolated Component (EPC):

[0028]

[0029] The edge permeation component is a method to identify network communities based on the location and connectivity of nodes in the network. It is not based on the degree of a single node, but on the connectivity of a collection of nodes. The higher the value of EPC, the more likely the other nodes connected to the node are to be connected to each other, and the more important the node is.

[0030] (3) Bottleneck Centrality (BN):

[0031] BN(v)=∑ s∈V p s (v)

[0032] Bottleneck centrality evaluates the importance of a node by calculating how often a node acts as a bottleneck in the shortest paths between different pairs of nodes in the network. The more times a node appears as a bottleneck in more shortest paths, the higher its bottleneck centrality value is, and it is considered to be more important in the network.

[0033] (4) Closeness Centrality (Clo):

[0034]

[0035] The higher the Clo value, the shorter the average distance from the node to other nodes in the network, and the higher the centrality of the node in the network.

[0036] (5) Radial centrality (Radiality, Rad):

[0037]

[0038] Radial centrality evaluates the centrality of a node based on its distance from the center of the network (usually minimizing the average distance from all nodes to the center). The radial score is the inverse of the distance from the node to the center of the network. The higher the Rad, the shorter the path from the node to other nodes in the network, and the better the connectivity and accessibility.

[0039] (6) Betweenness (BC):

[0040]

[0041] Betweenness centrality evaluates the importance of a node based on how often it appears on the shortest paths between all pairs of nodes in the network. The higher the BC score, the more shortest paths pass through the node, and the more important the node is.

[0042] (7) Stress centrality (Stress, Str):

[0043] Str(v)=∑ s≠t≠v∈C(v) σst (v)

[0044] Stress centrality takes into account the sum of the shortest paths from a node to other nodes, so the higher the stress centrality value of a node, the more important the connectivity of the node in the network.

[0045] (8) Coefficient of variation (Corr):

[0046]

[0047] The coefficient of variation can be used to indicate the degree of dispersion of gene expression. A larger CV value means that the expression level of the gene varies greatly in different tissues, disease states, or treatment conditions, and may be associated with a specific biological process or disease state.

[0048] (9) Phenotypic correlation (Coef):

[0049]

[0050] Phenotypic correlation is used to describe the degree of association between gene expression levels and phenotypes. A large Corr value indicates a strong association between a gene (or biological feature) and a specific phenotype, which may indicate that these genes play an important role in causing or regulating the phenotype.

[0051] The coexistence of highly correlated features will cause the information contained in these features to be overweighted, which will lead to an increasing loss of information contained in other features. To avoid this, we use the Spearman correlation analysis method to calculate the correlation between these features. If there are features with high correlation, only one of them will be selected to build the model.

[0052] Furthermore, the calculation formula of the conditional probability value of the node in the Bayes algorithm in step 8 is as follows:

[0053] P(Inf)=∏P(Inf|Attr)

[0054] Where P(Inf|Attr) represents the conditional probability value of a node given a certain feature; P(Inf) represents the joint probability product of all features of a node, that is, the gene importance value. P(Inf|Attr) represents the conditional probability value of a node given a certain feature.

[0055] The calculation formula of PageRank algorithm is as follows:

[0056]

[0057] Where P i (trf) represents the importance of node i. M(i) represents the set of genes connected to other genes i, Rij represents the Pearson correlation coefficient between nodes i and j, and degree(j) represents the degree of gene j. n represents the number of nodes that interact with node i, and d represents the probability of node i interacting with other adjacent nodes, which is generally set to 0.85.

[0058] In this formula, the importance of node i is obtained by the importance of all its neighboring nodes through their interactions. The importance of neighboring node j also receives all the interactions around it, and then passes the importance to node i through their interactions.

[0059] Furthermore, the prediction ability of different models was evaluated by the area under the receiver operating characteristic curve (AUC), sensitivity, specificity, Youden index and accuracy in step 11. A P value < 0.05 was statistically significant. Description of the drawings:

[0060] Figure 1 This is a flow chart of a method for identifying key nodes in a multi-layer regulatory network based on global interactive transitivity and feature fusion according to an embodiment of the present invention;

[0061] Figure 2 It is a lncRNA-miRNA-mRNA-TF ceRNA regulatory network diagram of lung adenocarcinoma in an embodiment of the present invention;

[0062] Figure 3 is the correlation between BPR and other features in the embodiment of the present invention;

[0063] Figure 4 It is the intersection of the lung adenocarcinoma ceRNA regulatory network nodes in the embodiment of the present invention and the LUAD-related nodes in the RNADisaese database;

[0064] Figure 5 The predictive capability of the BPR index constructed by a single feature or a combination of features in the embodiment of the present invention;

[0065] Figure 6 The BPR value distribution of the lung adenocarcinoma ceRNA regulatory network nodes in the optimal BPR index in the embodiment of the present invention;

[0066] Figure 7 It is the Wilcoxon signed rank test result of the BPR value between the positive sample and the negative sample in the embodiment of the present invention. Specific implementation method:

[0067] The following will take lung adenocarcinoma (LUAD) as an example and illustrate the embodiments of the present invention in conjunction with the accompanying drawings. Those skilled in the art can easily understand other advantages and effects of the present invention from the contents disclosed in this specification. The present invention can also be implemented or applied through other different specific examples, and the details in this specification can also be modified and changed in various ways based on different viewpoints and applications without departing from the spirit of the present invention.

[0068] like Figure 1 FIG. 1 is a flowchart of a method for identifying key nodes in a multi-layer regulatory network based on global interactive transitivity and feature fusion according to an embodiment of the present invention, comprising the following steps:

[0069] Step 1: Data download and preprocessing. Download RNA-seq (including mRNA, lncRNA and TF) and Isoform-Seq gene expression data of mature miRNA and clinical sample data of LUAD from the TCGA database, and perform logCPM normalization on the expression data of all nodes.

[0070] Step 2: Differential analysis. The R software package "Deseq2" was used to perform differential analysis of mRNA, miRNA, and lncRNA between adjacent normal samples and stage I cancer samples, and differential TFs were screened from the differential mRNAs.

[0071] Step 3: Construct miRNA-mRNA regulatory relationship pairs.

[0072] Step 4: Construct miRNA-lncRNA regulatory relationship pairs.

[0073] Step 5: Construct mRNA-TF interaction pairs.

[0074] Step 6: Network visualization: Cytoscape V 3.7.1 was used to construct and visualize the lncRNA-miRNA-mRNA-TF ceRNA regulatory network of four types of nodes.

[0075] Step 7: Calculation of eigenvalues ​​and feature selection of network nodes. Seven topological features, including degree, edge permeability component, bottleneck centrality, closeness centrality, radial centrality, betweenness centrality, and stress centrality, and two biological features, including coefficient of variation and phenotypic correlation, were selected and used for subsequent analysis.

[0076] Step 8: Construct the Bayes-PageRank (BPR) importance index. The Bayesian network structure is learned from the lncRNA-miRNA-mRNA-TF regulatory network constructed above by the score search algorithm. The features screened in step 7 are placed individually or randomly combined into the unsupervised Bayesian model for parameter learning. The conditional probability value of the node is obtained by using the maximum likelihood estimation for parameter estimation and the variable elimination algorithm for inference. It is input into the PageRank algorithm as a "weight". Various BPR indicators with single features or combined features are thus constructed.

[0077] Step 9: Obtain the optimal BPR importance index. Screen the experimentally verified miRNA, mRNA (including TF) and lncRNA molecules related to cancer through the RNADisease 4.0 database, and mark these molecules as positive samples; mark other molecules as negative samples. Calculate the AUC value of each BPR index based on the ROC curve. Evaluate the predictive ability of all BPR indicators in step 8 for positive samples, compare the AUC values ​​of different indicators, and select the optimal BPR indicator (BPR-best).

[0078] Step 10: Evaluate the correlation between the optimal BPR index and the key nodes of the disease. Calculate the BPR value of each node in the ceRNA regulatory network based on the optimal BPR index, and analyze the proportion of positive samples in different BPR value ranges, that is, power law distribution.

[0079] Step 11: Comparison of the prediction ability of the optimal BPR indicator and other single-feature or multi-feature indicators in the Bayesian model and the random forest model.

[0080] Step 12: Key node screening and literature verification: Based on the BPR value obtained by the optimal BPR indicator, the top 5% of nodes are selected as key nodes, and their relevance to the occurrence and development of cancer is verified using literature.

[0081] In an embodiment, a method for identifying key nodes in a multi-layer regulatory network based on global interactive transitivity and feature fusion. After preliminary filtering of the raw data in step 1, mRNA (including TFs) and lncRNA of 296 stage I cancer samples and 59 adjacent normal samples were retained; miRNA of 283 stage I cancer samples and 46 adjacent normal samples were retained. Each sample contained 18,749 mRNAs, 16,433 lncRNAs, and 2,213 mature miRNAs.

[0082] In an embodiment, a method for identifying key nodes in a multi-layer regulatory network based on global interactive transitivity and feature fusion. When performing differential analysis in step 2, the threshold |log2(FC)|≥1, adj.p≤0.05. According to the experimental database Hocomocov11 and the Animal TFDB database, 2092 human TFs were manually sorted. Based on these TFs, different TFs were screened from different mRNAs. There were 5105 differentially expressed mRNAs (including 512 TFs), of which 3362 were upregulated and 1741 were downregulated; 3879 of the 4879 differentially expressed lncRNAs were upregulated and 1000 were downregulated; 227 of the 411 differentially expressed miRNAs were upregulated and 184 were downregulated.

[0083] In an embodiment, a method for identifying key nodes in a multi-layer regulatory network based on global interactive transitivity and feature fusion is provided. In steps 3, 4, and 5, mRNA-miRNA, lncRNA-miRNA, and TF-mRNA regulatory pairs are found in multiple databases according to differentially expressed mRNA, miRNA, lncRNA, and TF, and the results are shown in Table 1.

[0084] Table 1 Statistics of regulatory relationships screened from the database

[0085]

[0086] In an embodiment, a method for identifying key nodes in a multi-layer regulatory network based on global interactive transitivity and feature fusion. In step 6, miRNA-mRNA and miRNA-lncRNA interaction pairs with p-value <0.05 and r <-0.3, as well as TFs-mRNA interaction pairs with |r| ≥ 0.3 and p-value <0.05 are retained. In addition, the phenotypic correlation (0 or 1) is calculated, and nodes with correlation r ≥ 0.1 and p-value <0.05 are retained. Finally, a total of 563 pairs of regulatory pairs were screened out, including 402 pairs of miRNA-mRNA, 56 pairs of lncRNA-miRNA, and 105 pairs of TF-mRNA, with a total of 319 nodes, including 89 miRNAs, 173 mRNAs, 5 lncRNAs, and 52 TFs.

[0087] like Figure 2 As shown, it is the lncRNA-miRNA-mRNA-TF ceRNA regulatory network diagram of lung adenocarcinoma in the embodiment of the present invention and the ceRNA regulatory network diagram after screening.

[0088] In an embodiment, a method for identifying key nodes in a multi-layer regulatory network based on global interactive transitivity and feature fusion is provided. In step 7, seven topological properties and two biological properties of nodes in the network are calculated, and the results are shown in Table 2. After all eigenvalues ​​are obtained, the correlation between the eigenvalues ​​is analyzed.

[0089] Table 2 Summary of attribute values ​​of each node in the network (partial)

[0090]

[0091] like Figure 3 As shown, the correlation between BPR and other features in the embodiment of the present invention is shown.

[0092] In the LUAD network, the topological features of the nodes, EPC, Radiality, and Closeness, are highly correlated (r>0.6), and Degree, Betweenness, and Stress are strongly correlated (r>0.89). Therefore, one of the features with higher correlation is selected for subsequent model construction research (EPC, Degree).

[0093] like Figure 4 As shown, it is the intersection of the lung adenocarcinoma ceRNA regulatory network nodes in the embodiment of the present invention and the LUAD related nodes in the RNADisaese database.

[0094] In an embodiment, a method for identifying key nodes of a multi-layer regulatory network based on global interactive transitivity and feature fusion is provided. In step 9, a total of 1518 molecules in the RNADisease database were confirmed to be associated with LUAD. In the constructed LUAD-specific ceRNA regulatory network, there were 100 positive samples and 219 negative samples.

[0095] like Figure 5 As shown, the prediction ability of the BPR index constructed by a single feature or a combination of features in the embodiment of the present invention is

[0096] The BPR index constructed by the single feature or any combination of features obtained in step 7 was evaluated for its ability to screen positive samples / negative samples. By comparing the AUC values ​​of all BPR indicators, it was found that the AUC values ​​of the BPR indicators constructed by the combined features were mostly higher than the AUC values ​​of the BPR indicators constructed by the single features. Among them, the BPR index obtained by the combination of four features had the best effect. The BPR index obtained by combining the four features of Degree+EPC+BottleNeck+Coef had the highest AUC value of 0.8553, which was the optimal BPR indicator.

[0097] This indicator not only reflects the connectivity of nodes, the higher the network connectivity, the more it reflects the importance of nodes, but also the key hubs between modules, which also reflects the importance of nodes. It also reflects that the greater the fluctuation of the biological expression of nodes, the more important it is for the occurrence and development of diseases. The model also reflects the transitivity of a node linking global network nodes, and in the network transmission, each node has a different "weight".

[0098] like Figure 6 As shown, the power-law distribution statistics of the BPR values ​​of the lung adenocarcinoma ceRNA regulatory network nodes in the optimal BPR index in an embodiment of the present invention.

[0099] In an embodiment, a method for identifying key nodes in a multi-layer regulatory network based on global interactive transitivity and feature fusion is provided. Evaluation of the optimal BPR index in step 10. The BPR values ​​of each node in the network obtained by using the optimal BPR index are statistically analyzed by power law distribution, and we only have a few nodes with high BPR values.

[0100] Among the 319 nodes, 250 (78.37%) have a BPR value less than 1, 69 (21.63%) have a BPR value greater than 1, 25 (7.84%) have a BPR value greater than 3, and only 10 (3.13%) have a BPR value greater than 5.

[0101] In addition, we found that the higher the BPR value of the node, the higher the proportion of key nodes related to LUAD.

[0102] Eight out of 10 nodes with BPR values ​​greater than 5 were known to be associated with LUAD (80%); 19 out of 25 nodes with BPR values ​​greater than 3 were known to be associated with LUAD (76.00%), and 109 out of 294 nodes with BPR values ​​less than 3 were known to be associated with LUAD (37.07%).

[0103] The distribution of BPR values ​​of nodes in the network also conforms to the power-law distribution characteristics of biological regulatory networks.

[0104] From the above, we can see that in LUAD, the node with a larger BPR value has a higher probability of being related to the disease and is more likely to be a key node.

[0105] like Figure 7 , which is the Wilcoxon signed rank test result of the BPR value between the positive sample and the negative sample in the embodiment of the present invention.

[0106] Through the statistical method of Wilcoxon signed rank test, the BPR value of the positive sample node is significantly higher than that of the negative sample node (p value < 0.0001), indicating that the BPR feature can effectively distinguish positive samples from negative samples.

[0107] In an embodiment, a method for identifying key nodes in a multi-layer regulatory network based on global interactive transitivity and feature fusion is provided. In step 11, the optimal BPR index is compared with other indexes in different models.

[0108] Table 3 shows the predictive ability of single feature index, multi-feature index and optimal BPR index in RF and Bayes models. It can be seen that in the random forest (RF) model, the predictive ability of the RF model constructed based on the optimal BPR index (AUC = 0.77, Specificity = 0.83, Youden index = 0.44) is higher than the best predictive ability of other single features and multi-feature combinations (AUC = 0.75, Specificity = 0.50-0.64, Youden index = 0.30-0.38). In the naive Bayesian (Bayes) model, the predictive ability of the Bayesian model constructed based on the optimal BPR index (AUC = 0.84, Sensitivity = 0.88, Specificity = 0.73, Youden index = 0.50) is higher than the best predictive ability of other single features and multi-feature combinations (AUC = 0.78-0.82, Sensitivity = 0.83-0.86, Specificity = 0.70-0.72, Youden index = 0.48).

[0109] Table 3 Comparison of the prediction capabilities of single-feature indicators, multi-feature indicators and the optimal BPR indicator in RF and Bayes models.

[0110]

[0111] Note: The black bold font represents the optimal value of the model evaluation index

[0112] In an embodiment, a method for identifying key nodes in a multi-layer regulatory network based on global interactive transitivity and feature fusion is provided. In step 12, key nodes are screened and verified by literature. The top 5% of nodes are selected according to the node BPR value obtained by the optimal BPR indicator. As shown in Table 4, 15 key nodes were screened in the LUAD network, of which 12 nodes have been included in the RNADisease database, and 3 have not been included in the database, but all have experimental literature to verify that they are related to the occurrence and development of LUAD.

[0113] Table 4 Nodes with top 5% BPR value

[0114]

[0115]

[0116] Note: Bold indicates key nodes known in the database

[0117] From the results in Table 5, it can be seen that most nodes are ranked in the top 15% under the three traditional topological methods (Degree, EPC, BottleNeck). Although some nodes perform poorly in phenotypic correlation (Coef) (hsa-miR-205-5p, 76.49%), they can still be effectively identified by the BPR method. Taking MALAT1 as an example, the node ranks at the top in Degree, EPC, and BottleNeck in the network.

[0118] In addition, there are 3 nodes (has-miR-93-5p, has-miR-326, has-miR-183-5p) that rank relatively low in traditional topological features, but perform well in Coef. These nodes are also identified as the top 5% important nodes through the optimal BPR indicator. For example, has-mir-326 is not prominent in Degree (17.24%), EPC (33.86%), and BottleNeck (19.75%), but performs well in Coef (5.64%). This shows that in addition to the topological features based on network structure, phenotypic correlation (Coef) also plays an important role in the identification of key nodes, which can make up for the shortcomings of traditional topological features.

[0119] Table 5. Ranking of traditional topological features (Degree, EPC, BottleNeck) and phenotypic correlation (Coef) corresponding to the top 5% of nodes in BPR value

[0120]

[0121] Note: Bold indicates key nodes known in the database

[0122] The above technical solution has at least the following advantages and beneficial effects compared with the prior art:

[0123] (1) Integrate multi-source data to construct a ceRNA regulatory network with four nodes of lncRNA-miRNA-mRNA-TF:

[0124] The present invention integrates information from different databases, such as miR2Disease, TarBase, miRTarbase, lncBase, StarBase, etc., when constructing the network, and can construct a more comprehensive and biologically rich ceRNA regulatory network. This integration helps to reveal the complex interactions and regulatory relationships between different types of RNA molecules (including lncRNA, miRNA, mRNA) and transcription factors (TF), thereby providing more comprehensive regulatory relationship information; and the quadruple node network has a larger amount of information and is more tolerant, and can recruit more disease-related factors and regulatory relationships into the regulatory network. Based on the quadruple node network to screen key molecules, non-single-type key molecules will be obtained based on a state close to the real regulation of the human body.

[0125] (2) Optimize the identification method of key nodes:

[0126] The present invention proposes a new node screening index, BPR, which comprehensively considers topological features, biological features and global interactive transmission. This method is not only based on the classic network topological properties, but also introduces the coefficient of variation (Corr) and phenotypic correlation (Coef) of molecular expression as important features. The single or multiple features of the node are combined by the Bayes algorithm, and the obtained Bayes value is used as the importance weight of the node. The PageRank algorithm is used to transmit under the interaction of the global network, and finally the importance value of all nodes, BPR, is obtained. This innovative weight distribution mechanism effectively overcomes the limitation of the traditional method that focuses on highly connected nodes, and can identify nodes with low connectivity but key roles, making up for the traditional method's neglect of the importance assessment of low-connectivity key nodes in multi-level gene regulatory networks. In addition, by comparing and verifying the prediction results of the optimal BPR index and other single-feature and multi-feature indicators in the Bayesian model and the random forest model, the significant advantages of this method in key node identification are demonstrated.

[0127] (3) Method evaluation and verification of key nodes:

[0128] The BPR index in the present invention is mainly used to identify key nodes in the regulatory network. Among the BPR values ​​obtained by the optimal BPR index, there are 15 nodes ranked in the top 5%, of which 12 nodes have been included in the RNADisease database, and although the other 3 nodes have not been included in the database, there are experimental literatures proving that they are closely related to the occurrence and development of lung adenocarcinoma (LUAD). This result shows that the optimal BPR index can effectively identify key nodes related to the occurrence and development of the disease in the molecular regulatory network of complex diseases, especially lung adenocarcinoma. The wide applicability of the BPR index provides strong support for disease network analysis, and provides a theoretical basis and technical support for further exploring the molecular mechanisms of complex diseases and screening potential disease markers and drug targets.

[0129] The description herein is combined with the drawings and specific embodiments to help understand the method and core idea of ​​the present invention. The method described in the present invention is not limited to the embodiments described in the specific embodiments. Other embodiments derived by those skilled in the art based on the method and ideas of the present invention also belong to the technical innovation scope of the present invention. The content of this specification should not be understood as limiting the present invention.

Claims

1. A method for identifying key nodes in a multi-layer regulatory network based on global interactive transitivity and feature fusion, characterized in that: The following steps are involved: Step 1: Data download and preprocessing. Download cancer RNA-seq (including mRNA, lncRNA and TF) and mature miRNA Isoform-Seq gene expression data and clinical sample data from the TCGA database, and perform logCPM normalization on the expression data of all nodes. Step 2: Differential analysis. The R package "Deseq2" was used to perform differential analysis of mRNA, miRNA, and lncRNA between adjacent normal samples and stage I cancer samples. Human TFs were sorted according to the experimental databases Hocomocov11 and Animal TFDB, and differential TFs were screened from differential mRNAs. Step 3: Construct miRNA-mRNA regulatory relationship pairs. Based on four biological experiment databases, mir2disease, miRecords, TarBase, miR2Disease, and one prediction database, StarBase, the interaction relationship pairs obtained in the above databases were crossed to obtain miRNA-mRNA regulatory pairs. Step 4: Construct miRNA-lncRNA regulatory relationship pairs. The interaction relationship pairs obtained from the biological experiment database LncBase and the prediction database StarBase were cross-referenced to obtain miRNA-lncRNA regulatory factor pairs. Step 5: Construct mRNA-TF interaction pairs. The biological experiment database TRRUST and the prediction database hTFtarget are crossed to obtain mRNA-TF interaction pairs. Step 6: Network visualization: Cytoscape V 3.7.1 was used to construct and visualize the lncRNA-miRNA-mRNA-TFceRNA regulatory network of four types of nodes. Step 7: Calculation of eigenvalues ​​and feature selection of network nodes. Seven topological features, including degree, edge permeability component, bottleneck centrality, closeness centrality, radial centrality, betweenness centrality, and stress centrality, and two biological features, including coefficient of variation and phenotypic correlation, were selected and used for subsequent analysis. Step 8: Construct the Bayes-PageRank (BPR) importance index. The Bayesian network structure is learned from the lncRNA-miRNA-mRNA-TF regulatory network constructed above by the score search algorithm. The features screened in step 8 are placed individually or randomly combined into the unsupervised Bayesian model for parameter learning. The conditional probability value of the node is obtained by using the maximum likelihood estimation for parameter estimation and the variable elimination algorithm for inference. It is input into the PageRank algorithm as a "weight". Various BPR indicators with single features or combined features are thus constructed. The calculation formula for the conditional probability value of a node in the Bayes algorithm is as follows: Where P(Inf|Attr) represents the conditional probability value of a node given a certain feature; P(Inf) represents the joint probability product of all features of a node, that is, the gene importance value. P(Inf|Attr) represents the conditional probability value of a node given a certain feature. The calculation formula of PageRank algorithm is as follows: Where P i (trf) represents the importance of node i. M(i) represents the set of genes connected to other genes i, R ij represents the Pearson correlation coefficient between nodes i and j, and degree(j) represents the degree of gene j. n represents the number of nodes that interact with node i, and d represents the probability of node i interacting with other adjacent nodes, which is generally set to 0.

85. Step 9: Obtain the optimal BPR importance index. Screen the experimentally verified miRNA, mRNA (including TF) and lncRNA molecules related to cancer through the RNADisease 4.0 database, and mark these molecules as positive samples; Mark other molecules as negative samples. Calculate the AUC value of each BPR indicator based on the ROC curve. Evaluate the predictive ability of all BPR indicators in step 8 for positive samples, compare the AUC values ​​of different indicators, and select the optimal BPR indicator (BPR-best). Step 10: Evaluate the correlation between the optimal BPR index and the key nodes of the disease. Calculate the BPR value of each node in the ceRNA regulatory network based on the optimal BPR index, and analyze the proportion of positive samples in different BPR value ranges, that is, power law distribution. Step 11: Comparison of the prediction ability of the optimal BPR indicator and other single-feature or multi-feature indicators in the Bayesian model and the random forest model. Step 12: Key node screening and literature verification: Select the top 5% nodes with the best BPR index as key nodes, and use literature to verify their relevance to the occurrence and development of cancer.

2. According to claim 1, a method for identifying key nodes in a multi-layer regulatory network based on global interactive transitivity and feature fusion, characterized in that In step 2, when performing differential analysis on mRNA, miRNA, and lncRNA, the threshold |log2(FC)|≥1 and adj.p≤0.

05.

3. According to claim 1, a method for identifying key nodes in a multi-layer regulatory network based on global interactive transitivity and feature fusion, characterized in that In steps 3, 4, and 5, data on the interactions between lncRNA, miRNA, mRNA, and TF were collected from various databases. mir2Disease, miRecords, Tar Base, miRTarBase, and LncBase mainly collect data from biological experiments. StarBase is mainly predicted by computer algorithms. This study screened out the interaction relationship pairs predicted by at least three algorithms from the StarBase database. Combining multiple experimental verification databases and prediction result databases, more comprehensive information on the regulatory relationship between miRNA-mRNA, TF-mRNA, and miRNA-lncRNA was collected and integrated.

4. According to claim 1, a method for identifying key nodes in a multi-layer regulatory network based on global interactive transitivity and feature fusion, characterized in that In step 6, in order to reduce the false positives of the obtained regulatory relationship pairs, the Pearson correlation coefficient was used to calculate the correlation between the above relationship pairs and the gene expression in the paired samples. Since TFs with weak binding ability showed stronger synergistic effects, different thresholds were selected for the above relationship pairs. In miRNA-mRNA and miRNA-lncRNA, all negative correlation pairs with p-value < 0.05 were retained (negative regulation of lncRNA and mRNA by miRNA). In TFs-mRNA, interaction pairs with |r| ≥ 0.2 and p-value < 0.05 were retained. In addition, the phenotypic correlation (0 or 1) was calculated, and the nodes with correlation r ≥ 0.1 and p-value < 0.05 were retained.

5. According to claim 1, a method for identifying key nodes in a multi-layer regulatory network based on global interactive transitivity and feature fusion, characterized in that The calculation formula in step 7 is as follows: (1) Degree (Deg): Deg(v)=|N(v)| The higher the degree of a node, the more nodes it is connected to, which means it is more important; (2) Edge Percolated Component (EPC): The edge permeation component is a method to identify network communities based on the location and connectivity of nodes in the network. It is not based on the degree of a single node, but on the connectivity of a collection of nodes; the higher the value of EPC, the more likely the other nodes connected to the node are to be connected to each other, and the more important the node is; (3) Bottleneck Centrality (BN): Bottleneck centrality evaluates the importance of a node by calculating the frequency with which a node acts as a bottleneck in the shortest paths between different pairs of nodes in the network; the more times a node appears as a bottleneck in more shortest paths, the higher its bottleneck centrality value is, and it is considered to be more important in the network; (4) Closeness Centrality (Clo): The higher the Clo value, the shorter the average distance from the node to other nodes in the network, and the higher the centrality of the node in the network; (5) Radial centrality (Radiality, Rad): The radial degree centrality evaluates the centrality of a node based on the distance from the node to the center of the network (usually the average distance from all nodes to the center is minimized); the radial degree score is the inverse of the distance from the node to the center of the network; the higher the RD, the shorter the path from the node to other nodes in the network, and the better the connectivity and accessibility. (6) Betweenness (BC): Betweenness centrality evaluates the importance of a node based on how often it appears on the shortest paths between all pairs of nodes in the network; The higher the BC score, the more shortest paths pass through the node, and the more important the node is. (7) Stress centrality (Stress, Str): Stress centrality considers the sum of the shortest paths from a node to other nodes, so the higher the stress centrality value of a node, the more important the connectivity of the node in the network; (8) Coefficient of Variation (CV): The coefficient of variation can be used to indicate the degree of discreteness of gene expression; a larger CV value means that the expression level of the gene varies greatly in different tissues, disease states, or treatment conditions, and may be associated with a specific biological process or disease state; (9) Phenotypic Correlation (Corr): Phenotypic correlation is used to describe the degree of association between gene expression levels and phenotypes; a large Corr value indicates a strong association between genes (or biological traits) and specific phenotypes, which may indicate that these genes play an important role in causing or regulating the phenotype. The coexistence of highly correlated features will cause the information contained in these features to be overweighted, which will lead to an increasing loss of information contained in other features. To avoid this, we use the Spearman correlation analysis method to calculate the correlation between these features. If there are features with high correlation, only one of them will be selected to build the model.

6. According to claim 1, a method for identifying key nodes in a multi-layer regulatory network based on global interactive transitivity and feature fusion, characterized in that: In step 8, various network topology features and biological features are used as Bayesian prior knowledge to constrain and regularize the data, effectively reducing the risk of overfitting. In the PageRank formula, the importance of node i is obtained by the importance of all its neighboring nodes through their interactions. The importance of neighboring node j also receives all the interactions around it, and then passes the importance to node i through their interactions. By adding the conditional probability value obtained by the Bayes algorithm as "weight" to the PageRank algorithm, the rules for transferring the importance between nodes are optimized, and the accuracy and stability of node importance evaluation are improved.

7. The method for identifying key nodes in a multi-layer regulatory network based on global interactive transitivity and feature fusion according to claim 1, characterized in that In step 11, the prediction effects of various indicators in the model were compared, mainly using the area under the receiver operating characteristic curve (AUC), sensitivity, specificity, Youden index and accuracy to evaluate the prediction performance of different models. A P value less than 0.05 was considered statistically significant.

Images