A method for identifying key nodes of a multi-layer regulatory network based on global interaction transitivity and feature fusion

Abstract

The application provides a kind of based on global interaction transmissibility and feature fusion multi-layer regulation network key node identification method, belongs to the field of bioinformatics and computational biology, method includes: by constructing lncRNA-miRNA-mRNA-TF ceRNA regulation network of tumor;Select classic network topological properties and biological properties (expression variation coefficient (Corr) and phenotype correlation (Coef)) as node features;Bayes algorithm is used to fuse node features, the Bayes value obtained by each node is used as the importance weight of the node, and the PageRank algorithm is used to transmit under the interaction of the global network, and further construct the BPR importance index of various feature combinations;The nodes related to tumors in the RNADisease database are used as positive samples to verify the prediction effect of the BPR index;Compare the prediction ability of single feature index, multi-feature index and optimal BPR index in RF and Bayes;The results show that the optimal BPR index (including Degree, EPC, BottleNeck and Coef) has the strongest ability in key node identification. Taking lung adenocarcinoma as an example, the application successfully identifies the key nodes related to lung adenocarcinoma progression, and shows its potential important value in biomedical research and clinical application.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of bioinformatics and computational biology, and particularly relates to a method for identifying key nodes in multi-level regulatory networks by combining network topological features, biological features and global interaction transitivity. BACKGROUND

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

[0003] Traditional key node identification methods mostly focus 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 network topological features to identify key nodes, but often ignore the influence of other biological attributes. Network topological features mainly rely on the neighbor information of nodes to evaluate their importance, such as degree (Degree), which assesses the importance of a node by the number of its neighbors. Other methods focus on the shortest path between nodes, including betweenness centrality (Betweenness), closeness centrality (Closeness), stress centrality (Stress), bottle neck centrality (BottleNeck) and radiality centrality (Radiality). In recent years, some methods based on the combination of multiple topological attributes have been proposed, such as the TOPSIS method, which evaluates the importance of nodes by combining multiple centrality indicators. Meanwhile, some researchers have proposed identifying key nodes in networks from a multi-level perspective. However, these methods generally ignore the global interaction transitivity and the importance of 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 permeation component (EPC) method, which assesses the connectivity of nodes by randomly deleting a small number of edges, and finds 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 attributes and ignore the biological features of nodes.

[0004] Based on this, we propose a comprehensive importance index-Bayes-PageRank (BPR) that can be applied to multi-level lncRNA-miRNA-mRNA-TF ceRNA regulatory network, as an effective tool to evaluate 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 the global interaction transitivity. SUMMARY

[0005] The purpose of the present application is to overcome the deficiencies in the prior art, and to provide a multi-level regulatory network key node identification method based on global interaction transitivity and feature fusion. The method comprises the following steps:

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

[0007] Step 2: Use the R package "Deseq2" to perform differential analysis on mRNA, miRNA and lncRNA between adjacent normal samples and I stage cancer samples, and select differential TFs from differential mRNA.

[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 lncRNA-miRNA-mRNA-TF ceRNA regulatory network of four types of nodes and visualize it.

[0012] Step 7: Network node feature value calculation and feature selection, select 7 topological features of nodes in the network, including degree, edge penetration component, bottleneck centrality, tightness centrality, radial degree centrality, betweenness centrality and stress centrality, and 2 biological features of nodes, including coefficient of variation and phenotype correlation, and select the above features for subsequent analysis.

[0013] Step 8: Constructing Bayes-PageRank (BPR) importance index. Learn the structure of the Bayesian network from the lncRNA-miRNA-mRNA-TF regulatory network constructed above by the scoring search algorithm. Put the features screened out in step 8 into the unsupervised Bayesian model for parameter learning, alone or randomly combined. Use maximum likelihood estimation for parameter estimation and variable elimination algorithm for inference to obtain the conditional probability value of the node. Input it as "weight" into the PageRank algorithm. Thus, various BPR indexes with single features or combined features are constructed.

[0014] Step 9: Obtain the optimal BPR importance index. Filter the experimentally verified miRNA, mRNA (including TF) and lncRNA molecules related to cancer from the RNADisease 4.0 database ( http: / / www.rnadisease.org / ) 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 prediction ability of all BPR indexes in step 8 for positive samples, compare the AUC values of different indexes, and select the optimal BPR index (BPR-best).

[0015] Step 10: Evaluate the correlation of 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 model, and analyze the proportion of positive samples in different BPR value ranges, i.e. power-law distribution.

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

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

[0018] Further, 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 collate human TFs, and screen differential TFs from differential mRNAs.

[0019] Further, 4 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 1 prediction database Starbase (starbase.sysu.edu.cn), the interaction pairs obtained from the above databases are crossed to obtain miRNA-mRNA regulatory pairs.

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

[0021] Further, the biological experimental database TRRUST www.grnpedia.org / trrust ) and the prediction database hTFtarget http: / / bioinfo.life.hust.edu.cn / hTFtarget ) are crossed to obtain mRNA-TF interaction pairs in step 5.

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

[0023] Further, step 7 includes:

[0024] (1) Degree (Deg):

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

[0026] The higher the degree of the node, the more nodes connected to the node, representing the higher importance of the node.

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

[0028]

[0029] EPC 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 set of nodes. The higher the EPC value, the more the other nodes connected to this node tend to be connected to each other, the more important the node is.

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

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

[0032] Bottleneck Centrality assesses the importance of a node by calculating the frequency of the node 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 the shortest paths, the higher its bottleneck centrality value, and the more important it is considered to be in the network.

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

[0034]

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

[0036] (5) Radiality Centrality (Rad):

[0037]

[0038] Radiality Centrality assesses the centrality of a node based on its distance to the center of the network (usually the node that minimizes the average distance of all nodes to the center). The radiality score is the inverse of the distance of the node to the center of the network. The higher the Rad, the shorter the path of the node to other nodes in the network, and the better the connectivity and accessibility.

[0039] (6) Betweenness Centrality (BC):

[0040]

[0041] Betweenness Centrality assesses the importance of a node based on the frequency of its appearance in 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 (Str):

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

[0044] 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.

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

[0046]

[0047] The coefficient of variation can be used to represent the degree of dispersion of gene expression. The larger the CV value, the greater the variation in the expression level of the gene under different tissues, disease states or treatment conditions, which may be associated with a specific biological process or disease state.

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

[0049]

[0050] Phenotype correlation is used to describe the degree of association between gene expression level and phenotype. A large Corr value indicates a strong association between the 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] Highly correlated features coexist, which can lead to an overemphasis on the information contained in these features, resulting in a loss of information contained in other features. To avoid this situation, we use the Spearman correlation analysis method to calculate the correlation between these features. If there are highly correlated features, only one of them is selected to build the model.

[0052] Further, the conditional probability value of the node in the Bayes algorithm in step 8 is calculated as follows:

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

[0054] Where P(Inf|Attr) represents the conditional probability value of the node given a certain feature; P(Inf) represents the joint probability product of all features of the node, i.e. the gene importance value. P(Inf|Attr) represents the conditional probability value of the 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 Pierson correlation coefficient between nodes i and j, degree(j) represents the degree of gene j. n represents the number of nodes interacting with node i, d represents the probability of node i interacting with other adjacent nodes, generally set to 0.85.

[0058] In the formula, the importance of node i is obtained by the importance of all adjacent nodes of node i through their interaction. The importance of adjacent node j also receives all the interactions around it, and then transmits the importance to node i through the interaction between them.

[0059] Further, the prediction ability of different models in step 11 is evaluated by area under the curve (AUC), sensitivity, specificity, Youden index and accuracy of the receiver operating characteristic curve. P<0.05 has statistical significance. BRIEF DESCRIPTION OF DRAWINGS

[0060] Figure 1 A flow chart of a key node identification method of a multi-layer regulatory network based on global interaction transitivity and feature fusion according to an embodiment of the present application;

[0061] Figure 2 A lncRNA-miRNA-mRNA-TF ceRNA regulatory network diagram of lung adenocarcinoma in an embodiment of the present application;

[0062] Figure 3 The correlation of BPR and other features in an embodiment of the present application;

[0063] Figure 4 The intersection of lung adenocarcinoma ceRNA regulatory network nodes and LUAD related nodes in RNADisaese database in an embodiment of the present application;

[0064] Figure 5 The prediction ability of the BPR index constructed by single feature or combined feature in an embodiment of the present application;

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

[0066] Figure 7 The Wilcoxon signed rank test result of BPR value between positive samples and negative samples in an embodiment of the present application. DETAILED DESCRIPTION

[0067] The present application will be described below with reference to lung adenocarcinoma (LUAD) as an example and in conjunction with the accompanying drawings, and other advantages and effects of the present application can be easily understood by those skilled in the art from the content disclosed in the present specification. The present application can also be implemented or applied through other different specific examples, and various modifications and changes can be made to the details in the present specification based on different views and applications without departing from the spirit of the present application.

[0068] As shown in FIG. 1, it is a flow chart of a key node identification method of a multi-layer regulatory network based on global interaction transitivity and feature fusion according to an embodiment of the present application, which comprises the following steps: Figure 1

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

[0070] Step 2: Differential analysis. The mRNA, miRNA and lncRNA between adjacent normal samples and I stage cancer samples are subjected to differential analysis using the R software package "Deseq2", and the differential TFs are screened from the differential mRNA.

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

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

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

[0074] Step 6: Network visualization. The lncRNA-miRNA-mRNA-TF ceRNA regulatory network of the four types of nodes is constructed and visualized by using CytoscapeV 3.7.1.

[0075] Step 7: Calculation of characteristic values of network nodes and feature selection. Seven topological features of the nodes in the network, including degree, edge penetration component, bottleneck centrality, tightness centrality, radial degree centrality, betweenness centrality and stress centrality, and two biological features of the nodes, including coefficient of variation and phenotype correlation, are selected, and the above features are screened for subsequent analysis.

[0076] ​Step 8: Constructing Bayes-PageRank (BPR) importance index. Learn the structure of the Bayesian network from the lncRNA-miRNA-mRNA-TF regulatory network constructed above by the scoring search algorithm. Put the features screened in step 7 into the unsupervised Bayesian model for parameter learning, either alone or randomly combined. Use maximum likelihood estimation for parameter estimation and variable elimination algorithm for inference to obtain the conditional probability value of the node. Input it as "weight" into the PageRank algorithm. Thus, various BPR indexes with single features or combined features are constructed.

[0077] Step 9: Obtain the optimal BPR importance index. Filter the experimentally verified miRNA, mRNA (including TF) and lncRNA molecules related to cancer from 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 prediction ability of all BPR indexes in step 8 for positive samples, compare the AUC values of different indexes, and select the optimal BPR index (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, i.e. power-law distribution.

[0079] Step 11: Comparison of the prediction ability of the optimal BPR index and other single feature or multi-feature indexes 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 index, select the top 5% nodes as key nodes, and verify their correlation with the occurrence and development of cancer by literature.

[0081] In an embodiment, a multi-layer regulatory network key node identification method based on global interaction transitivity and feature fusion. After preliminary filtering of the original data in step 1, 296 I-phase cancer samples and 59 adjacent normal samples of mRNA (including TFs) and lncRNA are retained; 283 I-phase cancer samples and 46 adjacent normal samples of miRNA are retained. Each sample contains 18,749 mRNA, 16,433 lncRNA and 2,213 mature miRNA.

[0082] In the embodiment, a multi-layer regulatory network key node identification method based on global interaction transitivity and feature fusion is provided. In step 2, when performing difference analysis, the threshold |log2(FC)|≥1, and adj.p≤0.05. According to the experimental database Hocomocov11 and the Animal TFDB database, 2092 human TFs are manually sorted out. According to these TFs, different TFs are screened from different mRNAs. The differentially expressed mRNAs have 5105 (including 512 TFs), of which 3362 are up-regulated, and 1741 are down-regulated; of the 4879 differentially expressed lncRNAs, 3879 are up-regulated, and 1000 are down-regulated; of the 411 differentially expressed miRNAs, 227 are up-regulated, and 184 are down-regulated.

[0083] In the embodiment, a multi-layer regulatory network key node identification method based on global interaction transitivity and feature fusion is provided. In steps 3, 4, and 5, according to the differentially expressed mRNAs, miRNAs, lncRNAs, and TFs, mRNA-miRNA, lncRNA-miRNA, and TF-mRNA regulatory pairs are found in multiple databases, and the results are shown in Table 1.

[0084] Table 1 Statistics of regulatory relationship pairs screened in the database

[0085]

[0086] In the embodiment, a multi-layer regulatory network key node identification method based on global interaction transitivity and feature fusion is provided. In step 6, miRNA-mRNA and miRNA-lncRNA interaction pairs with p-value<0.05 and r<-0.3, and TFs-mRNA interaction pairs with |r|≥0.3 and p-value<0.05 are retained. In addition, the phenotype 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 regulatory pairs are screened, including 402 miRNA-mRNA pairs, 56 lncRNA-miRNA pairs, and 105 TF-mRNA pairs, and a total of 319 nodes, including 89 miRNAs, 173 mRNAs, 5 lncRNAs, and 52 TFs.

[0087] As shown in Figure 2 , it is an lncRNA-miRNA-mRNA-TF ceRNA regulatory network diagram of lung adenocarcinoma in the embodiment of the application and a screened ceRNA regulatory network diagram.

[0088] In the embodiment, a multi-layer regulatory network key node identification method based on global interaction transitivity and feature fusion. In step 7, 7 topological properties and 2 biological properties of nodes in the network are calculated, and the results are shown in Table 2. After obtaining all the characteristic values, the correlation between the characteristic values is analyzed.

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

[0090]

[0091] As Figure 3 shown, the correlation of BPR and other features in the embodiment of the application

[0092] In the LUAD network, the correlation between EPC, Radiality and Closeness among the topological features of the nodes is high (r>0.6), and the correlation between Degree, Betweenness and Stress is strong (r>0.89). Therefore, in the subsequent model construction research, one of the highly correlated features is selected (EPC, Degree).

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

[0094] In the embodiment, a multi-layer regulatory network key node identification method based on global interaction transitivity and feature fusion. In step 9, there are 1518 molecules in the RNADisease database that are confirmed to be related to LUAD. In the constructed LUAD-specific ceRNA regulatory network, the positive samples are 100 and the negative samples are 219.

[0095] As Figure 5 shown, the prediction ability of the BPR index constructed by the single feature or the combined feature in the embodiment of the application

[0096] The BPR index constructed by the single feature or any combination of features obtained in step 7 is evaluated for positive sample / negative sample screening ability. By comparing the AUC values of all BPR indexes, it is found that the AUC values of the BPR indexes constructed by the combined features are higher than those of the BPR indexes constructed by the single features. Among them, the four feature combinations have the best effect, and the BPR index constructed by the combination of Degree, EPC, BottleNeck and Coef has the highest AUC value, which is 0.8553. This index is the optimal BPR index.

[0097] The index not only reflects the importance of the node, but also reflects the importance of the node in the module, and the greater the fluctuation of the biological expression amount of the node, the more important the node is to the occurrence and development of the disease. The model also reflects the transferability of a node linking global network nodes, and the different "weights" of each node in network transmission.

[0098] As shown in Figure 6 , it is the BPR value power-law distribution statistical situation of the lung adenocarcinoma ceRNA regulation network node in the optimal BPR index in the embodiment of the application.

[0099] In the embodiment, a multi-layer regulation network key node identification method based on global interaction transferability and feature fusion. In step 10, the evaluation of the optimal BPR index. The BPR values of each node in the network obtained by using the optimal BPR index are subjected to power-law distribution statistics, and only a few nodes have high BPR values.

[0100] Among the 319 nodes, 250 (78.37%) have BPR values less than 1, 69 (21.63%) have BPR values greater than 1, 25 (7.84%) have BPR values greater than 3, and only 10 (3.13%) have BPR values 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] Among the 10 nodes with BPR values greater than 5, 8 are known to be related to LUAD (80%); among the 25 nodes with BPR values greater than 3, 19 are known to be related to LUAD (76.00%), and among the 294 nodes with BPR values less than 3, 109 are related to LUAD nodes (37.07%).

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

[0104] As can be seen from the above, in LUAD, the nodes with higher BPR values are more likely to be related to the disease and are more likely to be key nodes.

[0105] As shown in Figure 7 , it is the Wilcoxon signed rank test result of the BPR value between the positive samples and the negative samples in the embodiment of the application.

[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 and negative samples.

[0107] In an embodiment, a method for identifying key nodes in a multi-layer regulatory network based on global interaction transitivity and feature fusion. In step 11, comparison of the optimal BPR indicator and other indicators in different models.

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

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

[0110]

[0111] Note: black bold is the optimal value in the model evaluation indicator

[0112] In an embodiment, a method for identifying key nodes in a multi-layer regulatory network based on global interaction transitivity and feature fusion. In step 12, key node screening and literature verification. According to the node BPR value obtained by the optimal BPR indicator, the top 5% nodes are selected. As shown in Table 4, 15 key nodes are screened in the LUAD network, of which 12 nodes have been included in the RNADisease database, and 3 nodes have not been included in the database, but all have experimental literature to verify their correlation with the development of LUAD.

[0113] Table 4 shows the top 5% nodes in BPR value

[0114]

[0115]

[0116] Note: Bold is the key node known in the database

[0117] From the results of Table 5, it can be seen that most nodes are in the top 15% under the three traditional topology methods (Degree, EPC, BottleNeck), although some nodes perform poorly in phenotype correlation (Coef) (hsa-miR-205-5p, 76.49%), they can still be effectively identified by the BPR method. Take MALAT1 as an example, the node ranks in 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 topology features, but perform well in Coef, and these nodes are also identified as the top 5% important nodes by the optimal BPR index. For example, has-mir-326 does not rank outstandingly 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 topology features based on network structure, the phenotype correlation (Coef) also plays an important role in the identification of key nodes, which can make up for the shortcomings of traditional topology features.

[0119] Table 5 The top 5% nodes in BPR value correspond to the ranking of traditional topology features (Degree, EPC, BottleNeck) and phenotype correlation (Coef)

[0120]

[0121] Note: Bold is the key node known in the database

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

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

[0124] The application integrates information from different databases such as miR2Disease, TarBase, miRTarbase, lncBase, StarBase and the like when constructing a network, so that a more comprehensive and biologically meaningful ceRNA regulation network can be constructed. Such 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 information amount of the four-node network is larger and more tolerant, and more disease-related factors and regulatory relationships can be recruited into the regulatory network. Based on the four-node network, key molecules are screened, and non-single type of key molecules is obtained based on the proximity to the real human regulation state.

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

[0126] The application proposes a new node screening index BPR which comprehensively considers the topological characteristics, biological characteristics and global interaction transmissibility. The method is not only based on the classic network topological properties, but also introduces the coefficient of variation (Corr) of molecular expression and the phenotype correlation (Coef) as important features. Through the Bayes algorithm, the single or multiple features of the node are combined, and the obtained Bayes value is used as the importance weight of the node. The PageRank algorithm is used for transmission under the interaction of the global network, and finally the importance value BPR of all nodes is obtained. This innovative weight distribution mechanism effectively overcomes the limitation of traditional methods which give priority to high connectivity nodes, and can identify low connectivity but key nodes, which makes up for the importance evaluation problem of ignoring low connectivity key nodes in traditional methods in multi-level gene regulation network. In addition, by comparing and verifying the prediction results of the optimal BPR index and other single feature and multi-feature indexes in the Bayes model and the random forest model, the significant advantage of the method in key node identification is shown.

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

[0128] The BPR index in the application is mainly used for identifying key nodes in the regulation network. In the BPR value obtained by the optimal BPR index, the top 5% nodes are a total of 15, of which 12 nodes have been included in the RNADisease database, and the other 3 nodes have not been included in the database, but experimental literature has proved 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 diseases in the molecular regulation network of complex diseases, especially in 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 mechanism of complex diseases and screening potential disease markers and drug targets.

[0129] The method of the present application is illustrated by the accompanying drawings and specific embodiments. The method of the present application is not limited to the embodiments described in the specific embodiments. Other embodiments derived from the method and idea of the present application by those skilled in the art also belong to the technical innovation range of the present application. The content of the specification should not be understood as a limitation of the present application.

Claims

1. A method for identifying key nodes in a multi-layer regulatory network based on global interactive transitivity and feature fusion, characterized by: The following steps are involved: Step 1: Data download and preprocessing: Download cancer RNA-seq data, including mRNA, lncRNA, and TF, as well as 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 organized based on the experimental databases Hocomocov11 and Animal TFDB, and differential TFs were screened from the 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; use the interaction relationship pairs obtained from the biological experiment database LncBase and the prediction database StarBase to cross-reference and obtain miRNA-lncRNA regulatory factor pairs; Step 5: Construct mRNA-TF interaction pairs; cross-reference the biological experiment database TRRUST and the prediction database hTFtarget 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: Calculate the eigenvalues ​​and select features of network nodes. Seven topological features, including degree, edge permeability component, bottleneck centrality, closeness centrality, radial degree centrality, betweenness centrality, and stress centrality, as well as two biological features, including coefficient of variation and phenotypic correlation, were selected and used for subsequent analysis. Step 8: Construct the Bayesian-PageRank importance index (BPR). Use the scoring search algorithm to learn the Bayesian network structure from the lncRNA-miRNA-mRNA-TF regulatory network constructed above. The features selected in step 8 are placed individually or randomly combined into an unsupervised Bayesian model for parameter learning. Maximum likelihood estimation is used for parameter estimation and variable elimination algorithm is used for inference to obtain the conditional probability value of the node. This is input as the "weight" into the PageRank algorithm to construct various BPR indicators with single features or combined features. The calculation formula for the conditional probability value of a node in the Bayes algorithm is as follows: P(Inf)=∏P(Inf|Attr) Where P(Inf|Attr) represents the conditional probability value of the node given a certain feature; P(Inf) represents the joint probability product of all features of the node, that is, the gene importance value; P(Inf|Attr) represents the conditional probability value of the node given a certain feature; The calculation formula of the PageRank algorithm is as follows: Among them, P i (trf) represents the importance of node i, M(i) represents the set of genes connected to other genes i, and R ij represents the Pearson correlation coefficient between nodes i and j, degree(j) represents the degree of gene j, n represents the number of nodes interacting 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 RNADisease 4.0 database for experimentally validated cancer-related miRNAs, mRNA-containing TFs, and lncRNA molecules, label these molecules as positive samples, and label 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 indexes in step 8 for positive samples, compare the AUC values ​​of different indexes, and select the optimal BPR index, namely BPR-best; Step 10: Evaluate the correlation between the optimal BPR index and the key nodes of the disease; Based on the optimal BPR index, the BPR value of each node in the ceRNA regulatory network was calculated, and the proportion of positive samples in different BPR value ranges was analyzed, that is, the power law distribution; Step 11: Comparison of the prediction ability of the optimal BPR index with 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% of nodes with the best BPR index as key nodes, and use literature to verify their correlation with the occurrence and development of cancer.

2. A 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 2, when performing differential analysis on mRNA, miRNA, and lncRNA, the threshold |log2(FC)| ≥ 1 and adj.p ≤ 0.

05.

3. A 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 steps 3, 4, and 5, data on the interactions between lncRNA, miRNA, mRNA, and TF were collected from various databases. mir2Disease, miRecords, Tar Base, and miRTarBase, while LncBase mainly collected data from biological experiments, were used, and StarBase was mainly used for predictions through computer algorithms. This study screened out interaction pairs predicted by at least three algorithms from the StarBase database, combined with multiple experimental verification databases and prediction result databases, to more comprehensively collect and integrate information on the regulatory relationships between miRNA-mRNA, TF-mRNA, and miRNA-lncRNA.

4. A 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 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; because 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, that is, the negative regulation of miRNA on lncRNA and mRNA; in TFs-mRNA, interaction pairs with |r| ≥ 0.2 and p-value < 0.05 were retained; in addition, phenotypic correlation was calculated, with a value of 0 or 1, and nodes with correlation r ≥ 0.1 and p-value < 0.05 were retained.

5. A 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 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 for identifying network communities based on the position and connectivity of nodes in the network. It is not based on the degree of a single node, but on the connectivity of a set of nodes. The higher the EPC value, the more likely the nodes connected to the node are to be connected to each other, and the more important the node is. (3) Bottleneck Centrality—Bottleneck, BN: BN(v)=∑ s∈V p s (in) Bottleneck centrality evaluates the importance of a node by calculating how often it 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—Closeness, Clo: The higher the Clo value, the shorter the average distance between the node and other nodes in the network, and the higher the centrality of the node in the network; (5) Radial centrality—Radiality, Rad: Radial degree centrality evaluates the centrality of a node based on its distance from the network center, usually minimizing the average distance from all nodes to the center. The radial degree score is the inverse of the distance from a node to the network center. The higher the RD, the shorter the path from the node to other nodes in the network, indicating better 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 passing through the node, and the more important the node is. (7) Stress centrality—Stress, Str: Str(v)=∑ s≠t≠v∈C(v) σ st (v) 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 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. (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 a gene or biological characteristic and a specific phenotype, 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 weight contained in these features to be too high, which in turn leads 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. 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 8, multiple network topology features and biological characteristics 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 transmits its importance to node i through their interactions. By adding the conditional probability value obtained by the Bayes algorithm as the "weight" into 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. A 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. The area under the receiver operating characteristic curve (AUC), sensitivity, specificity, Youden index, and accuracy were mainly used to evaluate the prediction performance of different models. A P value less than 0.05 was considered statistically significant.

Images