Key node identification method of disease marker expression regulation and control network

Abstract

The invention provides a key node identification method for a disease marker expression regulation network, and belongs to the technical field of disease markers, and the method comprises the steps: firstly carrying out the preprocessing and quality control of original data, including batch effect removal, abnormal sample identification and the like; then identifying differential expression genes through multiple difference analysis and a pre-training model, and constructing a gene expression correlation network; and integrating multi-source regulation and control data to construct a multi-level weighted network, calculating network node features, and carrying out representation learning and module division. And based on multi-dimensional features such as network topology features, module contribution degree and biological importance, a neural network model is trained to carry out key node identification. And finally, optimizing the model through multi-layer verification such as pathway enrichment, disease gene overlapping, expression stability, time sequence change and network disturbance, and finally obtaining a verified key node set. The problem that in the prior art, the interaction relation between molecules is ignored, and consequently some key regulation and control nodes are possibly missed is solved.

Description

Technical Field

[0001] The present invention belongs to the technical field of disease markers, and more particularly, relates to a method for identifying key nodes of a disease marker expression regulation network. Background Art

[0002] In recent years, with the continuous development of biomedical technologies, obtaining multi-dimensional biological data related to diseases using high-throughput omics technologies has become a common means for disease diagnosis and discovery of therapeutic targets. These omics data include gene expression profiles, protein expression profiles, metabolomics data, epigenomic data, etc., which can comprehensively reflect the molecular changes during the occurrence and development of diseases. By performing computational analysis on these omics data, key biomarkers related to diseases can be identified, providing a basis for early disease diagnosis and potential treatment.

[0003] Existing methods for identifying disease biomarkers mainly focus on differential expression analysis, functional enrichment analysis, and network analysis. Differential expression analysis methods screen candidate biomarkers with significant changes by comparing the gene or protein expression levels between disease and normal samples. Functional enrichment analysis further evaluates the biological processes and pathways involved by these differentially expressed genes or proteins to find biological functions closely related to the occurrence of diseases. Network analysis methods, on the other hand, depict the regulatory relationships between disease-related molecules at the system level to identify core regulatory nodes that play key roles during the occurrence of diseases. These methods have improved the discovery rate and reliability of biomarkers to a certain extent.

[0004] However, existing disease biomarker expression regulation networks, due to focusing on the functional annotation of differential molecules and ignoring the interaction relationships between molecules, may miss some key regulatory nodes, affecting the accuracy of key node identification. Summary of the Invention

[0005] In view of this, the present invention provides a method for identifying key nodes of a disease marker expression regulation network, which can solve the technical problem that the prior art ignores the interaction relationships between molecules, resulting in possible omission of some key regulatory nodes.

[0006] The present invention is implemented as follows:

[0007] The present invention provides a method for identifying key nodes of a disease biomarker expression regulation network, including performing data preprocessing on the original omics data related to disease tissue and biomarker samples to obtain expression stable components and expression variable components; calculating the inter-sample correlation of the expression components, constructing a sample quality score matrix, identifying and removing abnormal samples to obtain quality-controlled data; identifying significantly changed genes under disease states, obtaining differential gene importance scores, and screening a differential expression gene set; constructing an expression correlation coefficient matrix, screening significantly correlated gene pairs, and constructing a gene expression correlation network; integrating multiple types of data, calculating confidence coefficients, and constructing a multi-level weighted regulation network; constructing a network node feature matrix, calculating the contribution degree of network topological features; using a graph neural network model for representation learning, performing module division, and calculating the node time series contribution degree score; integrating multiple features to construct a node comprehensive scoring model to identify key nodes; and performing hierarchical verification on the key nodes, optimizing the model, and obtaining a finally verified key node set. Specifically, it includes the following steps:

[0008] S10. Perform data preprocessing on the original omics data related to disease tissue and its biomarker samples, including batch effect removal, missing value filling, data transformation, and data normalization, and obtain expression stable components and expression variable components through decomposition calculation;

[0009] S20. Calculate the inter-sample correlation of the expression stable components and the expression variable components, perform principal component analysis, construct a sample quality score matrix, identify and remove abnormal samples based on the sample quality score matrix to obtain quality-controlled data;

[0010] S30. Use a multiple differential expression analysis method to identify significantly changed genes under disease states, input the quality-controlled data into a pre-trained differential gene importance score model to obtain differential gene importance scores, and screen a differential expression gene set based on the differential gene importance scores and a differential fold change threshold;

[0011] S40. Based on the differential expression genes, construct an expression correlation coefficient matrix, calculate the stability coefficient of the expression correlation coefficient matrix, screen significantly correlated gene pairs according to the stability coefficient and a correlation threshold, and construct a gene expression correlation network;

[0012] S50. Integrate transcription factor binding site data, protein-protein interaction data, metabolite-protein interaction data, calculate the confidence coefficient of each type of data, and construct a multi-level weighted regulation network based on the confidence coefficient;

[0013] S60. Construct a network node feature matrix, where the network node feature matrix includes degree centrality, betweenness centrality, closeness centrality, and eigenvector centrality, and use the principal component analysis method to calculate the contribution degree of network topological features;

[0014] S70. Use a graph neural network model to perform representation learning on the multi-level weighted regulation network. Based on the results of the representation learning, use a community detection algorithm for module division and calculate the temporal contribution score of nodes in the module.

[0015] S80. Integrate the network topology feature contribution, the temporal contribution score of nodes in the module, the biological importance of nodes, and the relevance between nodes and diseases to construct a multi-layer neural network model, train to obtain a node comprehensive scoring model, and identify key nodes.

[0016] S90. Conduct hierarchical verification on the key nodes, including pathway enrichment analysis, disease gene overlap analysis, expression stability analysis, temporal variation analysis, and network perturbation analysis. Optimize the node comprehensive scoring model based on the results of the hierarchical verification to obtain a set of finally verified key nodes.

[0017] Among them, the original omics data is specifically multi-omics integrated data, including gene expression profile data, protein expression profile data, metabolomics data, epigenomics data, and clinical phenotype data.

[0018] Further, the sample - to - sample correlation between the expression stable component and the expression variable component specifically refers to calculating the expression correlation degree between samples using the Pearson correlation coefficient.

[0019] Further, the method for identifying significantly changed genes under disease states in the multiple differential expression analysis specifically combines three differential analysis methods, namely DESeq2, limma, and SAM, and integrates the results of multiple methods through Meta - analysis.

[0020] Further, the differential gene importance scoring model adopts an embedded deep - learning structure, including an embedded mathematical model layer, a feature extraction layer, an attention mechanism layer, a deep neural network layer, and an output layer.

[0021] Further, the differential gene importance scoring model adopts an embedded deep - learning structure, including an embedded mathematical model layer, a feature extraction layer, an attention mechanism layer, a deep neural network layer, and an output layer;

[0022] The embedded mathematical model layer is a system of mathematical equations, including a feature selection equation, a sample weight equation, and a network structure equation;

[0023] The feature selection equation is used for feature screening based on L1 regularization. The inputs include a gene expression matrix and phenotype data, and the output is a feature importance score;

[0024] The sample weight equation is used for sample weighting and integration. The inputs include sample quality scores and batch information, and the output is a sample weighting coefficient;

[0025] The network structure equation is used for network structure feature extraction. The inputs include a gene interaction network and pathway annotations, and the output is network topology features;

[0026] The feature extraction layer is used for data dimensionality reduction and feature learning. The input includes original expression data, and the output is a low-dimensional feature representation;

[0027] The attention mechanism layer is used for feature importance evaluation. The input includes feature representations, and the output is feature weights;

[0028] The deep neural network layer is used for non-linear feature transformation. The input includes weighted features, and the output is a high-level feature representation;

[0029] The output layer is used to output the importance scores of differentially expressed genes.

[0030] Furthermore, the differentially expressed gene set is screened based on the importance scores of the differentially expressed genes and the fold change threshold. Specifically, differentially expressed genes are selected according to the criteria that the absolute value of the fold change is greater than 1 and the significance is less than 0.05, combined with the sorting of importance scores.

[0031] Furthermore, significant related gene pairs are screened according to the stability coefficient and correlation threshold. Specifically, gene pairs with an absolute value of the correlation coefficient greater than 0.7 and a stability coefficient greater than 0.6 are set as significant related gene pairs.

[0032] Furthermore, the gene expression correlation network is constructed. Specifically, an adjacency matrix is constructed based on the screened significant related gene pairs, and the correlation coefficient is assigned as the edge weight.

[0033] Furthermore, the confidence coefficient of each type of data is calculated. Specifically, a comprehensive confidence coefficient is calculated based on data quality scores, data integrity, and experimental repeatability;

[0034] The multi-level weighted regulatory network is constructed based on the confidence coefficient. Specifically, different types of regulatory relationships are weighted and integrated into a heterogeneous network according to the confidence coefficient;

[0035] The graph neural network model is used for representation learning of the multi-level weighted regulatory network. Specifically, a graph convolutional neural network is used for node feature encoding and structural feature learning;

[0036] Module division is performed based on the representation learning results using a community detection algorithm. Specifically, the Louvain algorithm is used for community detection and module stability evaluation.

[0037] Compared with the prior art, a method for identifying key nodes of a disease biomarker expression regulation network provided by the present invention adopts a strategy of multi-omics data integration. First, preprocessing techniques such as batch effect correction and missing value imputation are used to convert the original data into a high-quality expression matrix. Then, the sample quality is evaluated by principal component analysis, and abnormal samples that may have problems are excluded. Next, a variety of differential expression analysis methods such as DESeq2, limma, and SAM are comprehensively used to identify genes that significantly change under disease states. An importance scoring model based on deep learning is further trained for these differential genes, and on the basis of fully considering factors such as sample quality and network structure, the importance scores of each differential gene are calculated. Based on the screened differential genes, the present invention constructs a gene expression correlation network, and conducts multi-level integration of the network, including transcriptional regulation, protein-protein interaction, and metabolic relationships, etc., and finally obtains a weighted network that comprehensively reflects disease molecular regulation. On this basis, the topological features of network nodes are extracted, and combined with information such as the temporal contribution degree, biological importance of nodes in network modules, and correlation with diseases, a multi-layer neural network model is trained to evaluate the overall importance of network nodes. By conducting hierarchical verification on key nodes, such as pathway enrichment analysis, disease gene overlap analysis, expression stability analysis, etc., the biological significance is further verified, and finally a set of key nodes with high credibility is screened out.

[0038] Compared with existing methods, the present invention has the following advantages: First, the strategy of multi-omics data integration is fully utilized, and the comprehensive application of multi-dimensional information is reflected in key steps such as data preprocessing, differential gene identification, and network construction, greatly improving the accuracy and reliability of the analysis results; Second, factors in multiple dimensions such as topological features, module contribution degree, and biological significance are introduced in the node importance evaluation, and a more comprehensive node scoring model is constructed, which can more accurately identify key regulatory nodes; It solves the technical problem that the prior art ignores the interaction relationship between molecules, resulting in possible omission of some key regulatory nodes. Brief Description of the Drawings

[0039] Figure 1 It is a flowchart of the method provided by the present invention;

[0040] Figure 2 It is a comparison diagram of batch effect correction effects in Example 2, including two subgraphs. Subgraph (a) is the batch effect before correction, and subgraph (b) is the batch effect after correction;

[0041] Figure 3 It is a scatter plot of principal component analysis in Example 2;

[0042] Figure 4 It is a heat map of differential gene expression in Example 2. Specific Embodiments

[0043] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings in the embodiments of the present invention.

[0044] As Figure 1 shown, it is a flowchart of a method for identifying key nodes of a disease biomarker expression regulation network provided by the present invention. This method includes the following steps:

[0045] S10. Perform data preprocessing on the original omics data related to the samples of disease tissues and their biomarkers, including batch effect removal, missing value filling, data transformation, and data normalization. After decomposition calculation, obtain the expression stable component and the expression variable component;

[0046] S20. Calculate the sample - to - sample correlation between the expression stable component and the expression variable component, perform principal component analysis, construct a sample quality score matrix, identify and remove abnormal samples based on the sample quality score matrix, and obtain the data after quality control;

[0047] S30. Use a multiple differential expression analysis method to identify genes with significant changes in the disease state. Input the data after quality control into a pre - trained differential gene importance score model to obtain the differential gene importance score, and screen the differential expression gene set based on the differential gene importance score and the fold - change threshold;

[0048] S40. Based on the differentially expressed genes, construct an expression correlation coefficient matrix, calculate the stability coefficient of the expression correlation coefficient matrix, screen significant correlation gene pairs according to the stability coefficient and the correlation threshold, and construct a gene expression correlation network;

[0049] S50. Integrate transcription factor binding site data, protein - protein interaction data, metabolite - protein interaction data, calculate the confidence coefficient of each type of data, and construct a multi - level weighted regulation network based on the confidence coefficient;

[0050] S60. Construct a network node feature matrix, where the network node feature matrix includes degree centrality, betweenness centrality, closeness centrality, and eigenvector centrality, and use the principal component analysis method to calculate the contribution degree of network topological features;

[0051] S70. Use a graph neural network model to perform representation learning on the multi - level weighted regulation network, perform module division using a community detection algorithm based on the representation learning results, and calculate the temporal contribution degree score of the node in the module;

[0052] S80, integrating the contribution of network topology features, the temporal contribution score of nodes in the module, the biological importance of nodes, and the correlation between nodes and diseases, building a multi-layer neural network model, training to obtain a node comprehensive scoring model, and identifying key nodes;

[0053] S90. Perform hierarchical verification on key nodes, including pathway enrichment analysis, disease gene overlap analysis, expression stability analysis, temporal change analysis, and network perturbation analysis. Optimize the node comprehensive scoring model based on the hierarchical verification results to obtain the final verified key node set.

[0054] The specific implementation methods of the above steps are described in detail below:

[0055] The specific implementation method of step S10 is: first, the original omics data is corrected for batch effects. A linear mixed model is used to characterize the batch effect, in which the sample expression value is used as the dependent variable, and the gene ID and batch ID are used as independent variables for modeling, so as to separate the influence of the batch effect and obtain the corrected expression matrix. Secondly, the missing values ​​in the original data are filled by the iterative singular value decomposition method. Specifically, the expression matrix is ​​decomposed into a left singular matrix, a singular value diagonal matrix, and a right singular matrix, and the missing values ​​are reconstructed using these components, and iterative optimization is performed until convergence. Next, the expression data is logarithmically transformed to alleviate the highly skewed distribution characteristics of the data. At the same time, the data is normalized by the Z-score standardization method to make the expression values ​​of each gene comparable. Finally, the total expression matrix is ​​decomposed into an expression stability component and an expression change component using the principal component analysis method. The sample correlation characteristics of the two types of components are characterized by calculating the Pearson correlation coefficient between samples.

[0056] The specific implementation of step S20 is: first, calculate the Pearson correlation coefficient matrix between samples to quantify the expression correlation between samples. Based on the correlation coefficient matrix, the principal component analysis method is used to extract the principal component score matrix, and the quality score of each sample is calculated according to the principal component score. The calculation formula of the sample quality score is: sample quality score = weighted sum of the principal component scores, where the weight coefficient is determined by the variance contribution rate of the principal component explanation. Then, a threshold is set according to the sample quality score, and abnormal samples are identified and eliminated to obtain quality-controlled data. The value of the threshold can be set empirically according to the distribution of the sample quality score and actual needs, and is usually 0.5 or 0.6.

[0057] The specific implementation of step S30 is as follows: First, multiple differential expression analysis methods such as DESeq2, limma, and SAM are used to analyze the quality-controlled data to identify genes with significant changes in the disease state. DESeq2 uses a negative binomial distribution model to characterize gene expression count data and uses a generalized linear model for parameter estimation and hypothesis testing; limma uses an empirical Bayesian model to fit the expression data, improving the detection ability in the case of small sample sizes; the SAM method evaluates gene differences based on the t-test statistic. The results of the three methods are integrated through Meta-analysis to obtain a more reliable set of differential genes. Second, the set of differential genes is input into a pre-trained deep learning model, which includes an embedded mathematics layer, a feature extraction layer, an attention mechanism layer, and a deep neural network layer, to calculate the differential importance score for each gene. The design of the differential gene importance score model fully considers various factors such as gene expression characteristics, network topology structure, and sample quality. Finally, according to the thresholds of fold change (|log2FC|>1) and FDR significance (FDR<0.05), the differential genes are screened and determined in combination with the differential gene importance score.

[0058] The specific implementation of step S40 is as follows: First, a gene expression correlation network is constructed based on the screened differential genes. Specifically, the Pearson correlation coefficient matrix between differential genes is calculated, and significant gene pairs with a correlation coefficient absolute value greater than 0.7 and a stability coefficient greater than 0.6 are selected as network edges. The calculation method of the stability coefficient is: by repeatedly sampling to calculate the correlation coefficient multiple times, and the frequency of the correlation coefficient higher than the threshold is statistically used as the stability coefficient. Second, the transcription factor binding site data, protein-protein interaction data, and metabolite-protein interaction data are integrated, and the confidence coefficient of each type of data is calculated respectively. The calculation of the confidence coefficient takes into account factors such as data quality score, data integrity, and experimental repeatability. Finally, according to the confidence coefficients of various types of data, different types of regulatory relationships are weighted and fused into a multi-level weighted regulatory network.

[0059] The specific implementation of step S50 is as follows: First, for the constructed multi-level weighted regulatory network, topological features such as degree centrality, betweenness centrality, closeness centrality, and eigenvector centrality of network nodes are extracted to construct a network node feature matrix. Then, the method of principal component analysis is used to calculate the relative contribution degrees of these network topological features, and the importance scores of each node in the network topological structure are obtained. Secondly, graph neural network models such as graph convolutional neural network are used for representation learning of the multi-level weighted regulatory network, so as to obtain the low-dimensional feature representation of each node. Based on these representation learning results, the Louvain algorithm is applied to partition the network into communities, and the temporal contribution degree scores of each node in its affiliated community are calculated. Finally, information from multiple dimensions such as the contribution degree of network topological features, the temporal contribution degree of nodes in modules, the biological importance of nodes, and the correlation between nodes and diseases are integrated into a multi-layer neural network model, and a node comprehensive scoring model is trained to identify key nodes.

[0060] The specific implementation of step S60 is as follows: First, based on the constructed multi-level weighted regulatory network, topological features such as degree centrality, betweenness centrality, closeness centrality, and eigenvector centrality of each node are extracted to construct a network node feature matrix. Degree centrality reflects the connectivity of nodes; betweenness centrality describes the hub role of nodes in the network; closeness centrality depicts the average distance from a node to other nodes; eigenvector centrality takes into account the importance of a node's neighbors. Then, the method of principal component analysis is used to calculate the relative contribution degrees of these network topological features to the importance of nodes, and the importance scores of each node in the network topological structure are obtained. This process can reveal the structural features of nodes in the network.

[0061] The specific implementation of step S70 is as follows: First, graph neural network models such as graph convolutional neural network are used for representation learning of the constructed multi-level weighted regulatory network. Graph neural networks can effectively learn the low-dimensional feature representations of graph-structured data and capture the structural and attribute information of nodes and their neighbors. Specifically, each layer of the graph neural network updates the representation of nodes according to the adjacency relationship and feature information of nodes, and finally outputs the low-dimensional feature vectors of each node. Secondly, based on these node feature representations, the Louvain algorithm is applied to partition the network into communities. The Louvain algorithm is a community detection algorithm based on modularity optimization that can quickly identify tightly connected modules in the network. At the same time, the temporal contribution degree scores of each node in its affiliated community can also be calculated to reflect the importance of nodes in the module. These steps aim to explore the biological significance of nodes from the perspectives of network topology and module organization.

[0062] The specific implementation of step S80 is as follows: First, construct a multi-layer neural network model. The input includes information in multiple dimensions such as the contribution degree of network topology features, the temporal contribution degree of nodes in the module, the biological importance of nodes, and the correlation between nodes and diseases. The output is the comprehensive score of the nodes. The multi-layer neural network can effectively fuse these heterogeneous features through non-linear transformation and learn the comprehensive evaluation function of node importance. Second, the obtained node comprehensive score model can be used to identify key nodes. Specifically, calculate the comprehensive scores for all nodes, and select the nodes with higher scores as key nodes after sorting. The calculation formula for the comprehensive score is: Comprehensive score = Contribution degree of network topology features × Weight 1 + Temporal contribution degree of nodes in the module × Weight 2 + Biological importance of nodes × Weight 3 + Correlation between nodes and diseases × Weight 4. These weight coefficients can be optimized through cross-validation.

[0063] The specific implementation of step S90 is as follows: First, conduct hierarchical verification on the key nodes identified in step S80. Specifically, it includes the following aspects: 1) Pathway enrichment analysis to evaluate whether the biological pathways where the key nodes are located are related to diseases; 2) Disease gene overlap analysis to check whether the key nodes overlap with known disease-related genes; 3) Expression stability analysis to evaluate the stability of the expression changes of the key nodes in different samples; 4) Temporal change analysis to explore the dynamic change trends of the key nodes during the disease development process; 5) Network perturbation analysis to observe the impact on the entire network topology structure after the key nodes are removed. Second, according to the results of the above hierarchical verification, optimize and adjust the node comprehensive score model, such as adjusting feature weights, adding new features, etc., and finally obtain a set of verified key nodes. These hierarchical verifications aim to evaluate the biological significance and reliability of the key nodes from multiple perspectives and ensure the strong credibility of the identification of key nodes.

[0064] The following provides a specific description of the calculation process or mathematical model involved in the present invention.

[0065] 1. Mathematical expression in the data preprocessing stage (S10):

[0066] The mathematical model for removing batch effects is as follows:

[0067] Y ijk =μ i +β j +γ ijk +∈ ijk ;

[0068] In the formula, Y ijk represents the expression value of the i-th gene in the k-th sample of the j-th batch; μ i represents the average expression level of the i-th gene; β jrepresents the effect of the j-th batch; γ ijk represents biological variation; ∈ ijk represents random error, following the normal distribution N(0, σ 2 ).

[0069] Missing value filling adopts the iterative SVD decomposition method:

[0070] X = U∑V T + E;

[0071] where X is the expression matrix; U is the left singular matrix; ∑ is the diagonal matrix of singular values; V is the right singular matrix; E is the residual matrix.

[0072] Data transformation adopts logarithmic transformation:

[0073] X transformed = log2(X raw + α);

[0074] where X raw is the original expression value; α is the smoothing factor, and its value range is 0.1 - 1.

[0075] Data normalization adopts Z-score standardization:

[0076]

[0077] where Z ij is the standardized expression value; X ij is the original expression value; μ i is the mean of the i-th feature; σ i is the standard deviation of the i-th feature.

[0078] Expression component decomposition model:

[0079] X = X stable + X variable + ∈;

[0080] where X is the total expression matrix; X stable is the stable expression component; X variable is the variable expression component; ∈ is the noise term.

[0081] 2. Mathematical expression of sample quality assessment (S20):

[0082] Calculation of sample - to - sample correlation:

[0083]

[0084] where r ij is the Pearson correlation coefficient between samples i and j; X kiis the expression value of the k-th gene in sample i; μ i is the average expression value of sample i.

[0085] Principal component analysis model:

[0086] T = XW;

[0087] In the formula, T is the principal component score matrix; X is the centered data matrix; W is the eigenvector matrix.

[0088] Calculation of sample quality score:

[0089]

[0090] In the formula, Q i is the quality score of the i-th sample; w j is the weight of the j-th principal component; PC ij is the score of the i-th sample on the j-th principal component.

[0091] 3. Mathematical expression of differential gene identification (S30):

[0092] DESeq2 model:

[0093] K ij ~NB(μ ij ,α i );

[0094] log2(μ ij ) = β i0 +β i1 x j +log2(s j );

[0095] In the formula, K ij is the count value; μ ij is the mean parameter; α i is the dispersion parameter; β i0 is the baseline expression level; β i1 is the treatment effect; s j is the scale factor.

[0096] Limma model:

[0097] E[y g = Xα g ;

[0098]

[0099] In the formula, y g is the expression vector of gene g; X is the design matrix; α g is the coefficient vector; is the variance parameter; w gi is the weight.

[0100] SAM analysis model:

[0101]

[0102] In the formula, d i is the statistic of gene i; is the mean of the treatment group; is the mean of the control group; s i is the standard deviation; s0 is the adjustment factor.

[0103] Meta-analysis integration model:

[0104]

[0105] In the formula, Z combined is the comprehensive statistic; w k is the weight of the k-th method; Z k is the standardized statistic of the k-th method.

[0106] 4. Differential gene importance scoring model (feature selection equation):

[0107] S i =β1||θ i ||1 + β2∑ j∈N(i) w ij S j +β3f(X i ) + ∈ i ;

[0108] In the formula, S i is the importance score of gene i; θ i is the L1 regularization coefficient vector; w ij is the network edge weight; N(i) is the neighbor set of gene i; f(X i ) is the expression feature function; β1, β2, β3 are weight parameters; ∈ i is the error term.

[0109] Sample weight equation:

[0110] W s =α1Q s +α2B s +α3T s ;

[0111] In the formula, W s is the weight of sample s; Q s is the quality score; B s is the batch effect score; T sis the temporal correlation score; α1, α2, α3 are adjustment parameters.

[0112] Network structure equation:

[0113] F i = γ1C i + γ2B i + γ3E i + γ4P i ;

[0114] In the formula, F i is the network feature of node i; C i is the centrality index; B i is the bottleneck coefficient; E i is the modularity index; P i is the pathway enrichment score; γ1, γ2, γ3, γ4 are weight coefficients.

[0115] 5. Mathematical expression for constructing the expression correlation network (S40):

[0116] Calculation of the correlation coefficient matrix:

[0117]

[0118] In the formula, R ij is the correlation coefficient between genes i and j; cov(X i , X j ) is the covariance; σ i , σ j are the standard deviations.

[0119] Calculation of the stability coefficient:

[0120]

[0121] In the formula, S ij is the stability coefficient of the gene pair (i, j); N is the number of resampling times; is the correlation coefficient for the k-th sampling; θ is the correlation threshold; I() is the indicator function.

[0122] 6. Mathematical expression for constructing the multi-level weighted regulatory network (S50):

[0123] Calculation of the confidence coefficient:

[0124] C k = λ1Q k + λ2P k + λ3R k ;

[0125] In the formula, C k is the confidence of the k-th type of data; Q kis the quality score; P k is the integrity score; R k is the repeatability score; λ1, λ2, λ3 are the weight parameters.

[0126] Multi-layer network integration:

[0127]

[0128] where M ij is the integrated edge weight; C k is the confidence of the k-th layer network; is the element of the adjacency matrix of the k-th layer network.

[0129] 7. Mathematical expression for evaluating the importance of network nodes (S60 - S70):

[0130] Degree centrality:

[0131]

[0132] where DC i is the degree centrality of node i; A ij is the element of the adjacency matrix.

[0133] Betweenness centrality:

[0134]

[0135] where BC i is the betweenness centrality of node i; σ st is the number of shortest paths from node s to t; σ st (i) is the number of shortest paths passing through node i.

[0136] Closeness centrality:

[0137]

[0138] where CC i is the closeness centrality of node i; d ij is the shortest distance from node i to j; n is the total number of nodes.

[0139] Eigenvector centrality:

[0140]

[0141] where EC i is the eigenvector centrality of node i; λ is the eigenvalue; A ij is the element of the adjacency matrix.

[0142] Graph neural network representation learning:

[0143]

[0144] In the formula, H l is the representation of the l-th layer node; A is the adjacency matrix; D is the degree matrix; W l is the weight matrix; σ is the activation function.

[0145] 8. Mathematical expression of the node comprehensive scoring model (S80):

[0146] Comprehensive scoring calculation:

[0147] Score i = ω1T i + ω2M i + ω3B i + ω4D i ;

[0148] In the formula, Score i is the comprehensive score of node i; T i is the topological feature score; M i is the module contribution degree; B i is the biological importance; D i is the disease relevance; ω1, ω2, ω3, ω4 are weight parameters.

[0149] The parameters of all the above equations can be obtained in the following ways:

[0150] 1. Network parameters (w ij , A ij etc.) are obtained by constructing a network with experimental data;

[0151] 2. Weight parameters (α, β, γ, λ, ω, etc.) are obtained by optimizing through cross-validation;

[0152] 3. Threshold parameters (θ, etc.) are determined by empirical setting or data distribution;

[0153] 4. Quality scores, etc. (Q s , Q k etc.) are obtained by calculating through the data quality control process;

[0154] 5. Biological parameters are obtained through database annotation or experimental verification.

[0155] The principles and meanings of these equations are as follows:

[0156] 1. The data preprocessing equations consider batch effects, technical noise, and biological variations, and improve data quality through decomposition and standardization;

[0157] 2. The quality assessment equation is based on sample correlation and principal component analysis to accurately identify abnormal samples;

[0158] 3. The differential analysis equation integrates multiple statistical methods to improve the reliability of differential gene identification;

[0159] 4. The network construction equation takes into account multi-level biological relationships and provides a more complete view of the regulatory network;

[0160] 5. The node importance evaluation equation characterizes node features from multiple dimensions to accurately identify key nodes.

[0161] Specifically, the principle of the present invention is as follows: First, in the data preprocessing stage, the present invention uses technical means such as batch effect correction, missing value imputation, data transformation, and normalization for the original omics data, and converts it into a high-quality expression matrix. Among them, batch effect correction uses a linear mixed model to fit batch effect parameters, thereby eliminating the influence of technical noise; missing value imputation uses the method of singular value decomposition, which can retain the inherent structure of the data and improve the accuracy of imputation; the purpose of data transformation and normalization is to alleviate the highly skewed distribution and incomparability of expression data. Through these preprocessing steps, a high-quality expression matrix is obtained, laying a foundation for subsequent analysis.

[0162] Secondly, in the sample quality assessment stage, the present invention uses the method of principal component analysis to calculate the quality score of each sample based on the expression correlation between samples. This scoring method can effectively identify abnormal samples that may have problems and provide high-quality data input for subsequent differential analysis.

[0163] Next, in the differential gene analysis stage, the present invention comprehensively uses multiple statistical methods such as DESeq2, limma, and SAM to identify genes that change significantly under disease conditions. These methods are respectively based on the negative binomial distribution model, empirical Bayesian model, and t-test statistic, and can effectively capture differential features in different types of data. By integrating these results through Meta-analysis, the reliability of differential gene identification can be further improved. In addition, a differential gene importance scoring model based on deep learning is also trained. This model takes into account multiple factors such as gene expression characteristics, network topology, and sample quality, and can more accurately evaluate the key role of each differential gene in the disease.

[0164] In the network construction stage, the present invention first constructs a gene co-expression correlation network based on the screened differential genes. By calculating the Pearson correlation coefficient between genes and setting a correlation threshold, a network reflecting gene co-expression relationships is obtained. On this basis, multi-level biological regulatory relationships such as transcriptional regulation, protein-protein interaction, and metabolic relationships are integrated to construct a comprehensive weighted regulatory network. This multi-level network integration method can more comprehensively characterize the complex regulatory mechanisms between disease-related molecules.

[0165] Finally, in the node importance evaluation stage, topological features such as degree centrality, betweenness centrality, closeness centrality, and eigenvector centrality of network nodes were extracted, and combined with information such as the temporal contribution degree, biological importance, and disease relevance of nodes in network modules, a multi-layer neural network model was trained. This model can comprehensively evaluate the overall importance of network nodes and provide a basis for the identification of key nodes. By performing hierarchical verification on key nodes, such as pathway enrichment analysis and disease gene overlap analysis, the key role of these nodes in the disease occurrence mechanism was further verified, ensuring the reliability of the results.

[0166] A specific embodiment 1 of the present invention is provided below. The specific implementation manners of each step in this embodiment 1 are described in detail as follows: The specific implementation manner of step S10 is as follows:

[0167] First, batch effect correction was performed on the original omics data. The following linear mixed model was used to characterize the batch effect:

[0168] Y ijk = μ i + β j + γ ijk + ∈ ijk ;

[0169] where Y ijk represents the expression value of the i-th gene in the k-th sample of the j-th batch; μ i represents the average expression level of the i-th gene; β j represents the effect of the j-th batch; γ ijk represents biological variation; ∈ ijk represents random error, which follows the normal distribution N(0, σ 2 ). By fitting this model, the influence of the batch effect can be separated and the corrected expression matrix can be obtained.

[0170] Secondly, for the missing values in the original data, the method of iterative singular value decomposition was used for filling. Specifically, the expression matrix X was decomposed into the left singular matrix U, the singular value diagonal matrix ∑, and the right singular matrix V T , as follows:

[0171] X = U∑V T + E;

[0172] where E is the residual matrix. The missing values were reconstructed using these components and iteratively optimized until convergence.

[0173] Next, logarithmic transformation was performed on the expression data, and the formula is as follows:

[0174] X transformed = log2(Xraw +α);

[0175] Among them, X raw is the original expression value, α is the smoothing factor, and its value range is 0.1 to 1. At the same time, the Z-score normalization method is used to normalize the data:

[0176]

[0177] Among them, Z ij is the normalized expression value; X ij is the original expression value; μ i is the mean of the i-th feature; σ i is the standard deviation of the i-th feature.

[0178] Finally, the total expression matrix X is decomposed into a stable expression component X stable and a variable expression component X variable , as follows:

[0179] X = X stable + X variable + ∈;

[0180] Among them, ∈ is the noise term. The sample correlation characteristics of the two types of components are characterized by calculating the Pearson correlation coefficient between samples:

[0181]

[0182] Among them, r ij is the Pearson correlation coefficient between samples i and j; X ki is the expression value of the k-th gene in sample i; μ i is the average expression value of sample i.

[0183] The specific implementation manner of step S20 is as follows:

[0184] First, calculate the Pearson correlation coefficient matrix between samples to quantify the expression correlation between samples. Based on the correlation coefficient matrix, the principal component score matrix T is extracted by using the principal component analysis method, as follows:

[0185] T = XW;

[0186] Among them, X is the centered data matrix; W is the eigenvector matrix. Then, calculate the quality score Q of each sample according to the principal component scores i :

[0187]

[0188] Among them, w j is the weight of the j-th principal component; PCij is the score of the i-th sample on the j-th principal component.

[0189] Next, a threshold is set according to the sample quality score to identify and remove abnormal samples, thereby obtaining the quality-controlled data. The value of the threshold can be empirically set according to the distribution of the sample quality score and actual requirements, usually taking values of 0.5 or 0.6.

[0190] The specific implementation manner of step S30 is as follows:

[0191] First, multiple differential expression analysis methods such as DESeq2, limma, and SAM are used to analyze the quality-controlled data to identify genes that change significantly under disease states. DESeq2 uses a negative binomial distribution model to characterize gene expression count data, and uses a generalized linear model for parameter estimation and hypothesis testing. The formula is as follows:

[0192] K ij ~NB(μ ij ,α i );

[0193] log2(μ ij ) = β i0 +β i1 x j +log2(s j );

[0194] where K ij is the count value; μ ij is the mean parameter; α i is the dispersion parameter; β i0 is the baseline expression level; β i1 is the treatment effect; s j is the scale factor.

[0195] limma uses an empirical Bayes model to fit the expression data. The formula is as follows:

[0196] E[y g = Xα g ;

[0197]

[0198] where y g is the expression vector of gene g; X is the design matrix; α g is the coefficient vector; is the variance parameter; w gi is the weight.

[0199] The SAM method evaluates gene differences based on the t-test statistic. The formula is as follows:

[0200]

[0201] Among them, d i is the statistic of gene i; is the mean of the treatment group; is the mean of the control group; s i is the standard deviation; s0 is the adjustment factor.

[0202] The results of the three methods are integrated by Meta-analysis, and the formula is as follows:

[0203]

[0204] Among them, Z combined is the combined statistic; w k is the weight of the k-th method; Z k is the standardized statistic of the k-th method.

[0205] Secondly, the differentially expressed gene set is input into a pre-trained deep learning model, which includes an embedded mathematics layer, a feature extraction layer, an attention mechanism layer, and a deep neural network layer, for calculating the differential importance score S i of each gene. The design formula of the differentially expressed gene importance score model is as follows:

[0206] S i =β1||θ i ||1 + β2∑ j∈N(i) w ij S j +β3f(X i ) + ∈ i ;

[0207] Among them, θ i is the L1 regularization coefficient vector; w ij is the network edge weight; N(i) is the neighbor set of gene i; f(X i ) is the expression feature function; β1, β2, β3 are weight parameters; ∈ i is the error term.

[0208] Finally, according to the fold change (|log2FC|>1) and FDR significance (FDR<0.05) thresholds, combined with the differentially expressed gene importance score, the differentially expressed genes are screened and determined.

[0209] The specific implementation manner of step S40 is as follows:

[0210] First, a gene expression correlation network is constructed based on the screened differentially expressed genes. Calculate the Pearson correlation coefficient matrix R between the differentially expressed genes:

[0211]

[0212] where R ij is the correlation coefficient between genes i and j; cov(X i , X j ) is the covariance; σ i , σ j are the standard deviations. According to the criterion that the absolute value of the correlation coefficient is greater than 0.7 and the stability coefficient S ij is greater than 0.6, significant relevant gene pairs are screened out as network edges:

[0213]

[0214] where N is the number of resampling times; is the correlation coefficient of the k-th sampling; θ is the correlation threshold; I() is the indicator function.

[0215] Secondly, integrate the transcription factor binding site data, protein-protein interaction data, and metabolite-protein interaction data, and calculate the confidence coefficient C k for each type of data respectively:

[0216] C k = λ1Q k + λ2P k + λ3R k ;

[0217] where Q k is the quality score; P k is the integrity score; R k is the repeatability score; λ1, λ2, λ3 are weight parameters. According to the confidence coefficients of various types of data, different types of regulatory relationships are weighted and fused into a multi-level weighted regulatory network, and its edge weight M ij is calculated as follows:

[0218]

[0219] where is the element of the adjacency matrix of the k-th layer network.

[0220] The specific implementation manner of step S50 is as follows:

[0221] First, for the constructed multi-level weighted regulatory network, extract topological features such as the degree centrality DC i , betweenness centrality BC i , closeness centrality CC i and eigenvector centrality EC i of the network nodes, and construct a network node feature matrix:

[0222]

[0223] Among them, A ij is an adjacency matrix element; σ st is the number of shortest paths from node s to t; σ st (i) is the number of shortest paths passing through node i; d ij is the shortest distance from node i to j; λ is the eigenvalue. Then, the relative contribution degrees of these network topology features are calculated by the method of principal component analysis to obtain the importance scores of each node in the network topology structure.

[0224] Secondly, use the graph neural network model to perform representation learning on the multi-level weighted regulation network:

[0225]

[0226] Among them, H l is the node representation of the l-th layer; A is the adjacency matrix; D is the degree matrix; W l is the weight matrix; σ is the activation function. Based on these representation learning results, the Louvain algorithm is applied to partition the network, and the temporal contribution degree scores of each node in its affiliated community are calculated.

[0227] The specific implementation manner of step S60 is as follows:

[0228] Construct a multi-layer neural network model, the input includes the contribution degree T of network topology features i , the temporal contribution degree M of the node in the module i , the biological importance B of the node i and the disease correlation D of the node i and other multi-dimensional information, and the output is the comprehensive score Score of the node i , and the formula is as follows:

[0229] Score i = ω1T i + ω2M i + ω3B i + ω4D i ;

[0230] Among them, ω1, ω2, ω3, ω4 are weight parameters. The multi-layer neural network can effectively fuse these heterogeneous features through non-linear transformation and learn the comprehensive evaluation function of node importance. These weight parameters can be optimized by cross-validation.

[0231] Finally, calculate the comprehensive scores for all nodes, and select the nodes with higher scores as key nodes after sorting.

[0232] The specific implementation manner of step S90 is as follows:

[0233] Perform hierarchical verification on the key nodes identified in step S80. First, evaluate whether the biological pathways where the key nodes are located are related to the disease through pathway enrichment analysis. Second, check whether the key nodes overlap with known disease-related genes. Third, evaluate the stability of the expression changes of the key nodes in different samples. Fourth, explore the dynamic change trends of the key nodes during the disease development process. Finally, observe the impact on the entire network topology after removing the key nodes.

[0234] According to the results of the above hierarchical verification, optimize and adjust the node comprehensive scoring model, such as adjusting feature weights, adding new features, etc., and finally obtain a set of verified key nodes.

[0235] Generally speaking, the key node identification method proposed by the present invention involves key steps such as data preprocessing, sample quality assessment, differential gene analysis, network construction, and node importance assessment. Through technical means such as mathematical modeling, machine learning, and bioinformatics analysis, the importance of network nodes is comprehensively evaluated, providing a systematic and effective method for identifying key nodes in the disease biomarker expression regulation network. This method utilizes multi-omics data, fully considering various factors such as batch effects, sample quality, differential expression, and network topology. While improving the accuracy of key node identification, it also enhances the reliability and biological significance of the results through hierarchical verification.

[0236] To better understand and implement the present invention, Example 2 of a specific application scenario of the present invention is provided below: A research team intends to systematically identify key biomarkers of Alzheimer’s Disease (AD) based on the method proposed by the present invention. First, they collected gene expression data, protein expression data, metabolomics data, and epigenomics data of AD patients and healthy control groups from public databases, totaling 4 groups of omics data.

[0237] In the data preprocessing stage, the research team first corrected for batch effects. Taking gene expression data as an example, they established the following linear mixed model:

[0238] Y ijk =μ i +β j +γ ijk +∈ ijk ;

[0239] where, Y ijk represents the expression value of the i-th gene in the k-th sample of the j-th batch; μ i represents the average expression level of the i-th gene; β j represents the effect of the j-th batch; γ ijk represents biological variation; ∈ ijkrepresents the random error, which follows the normal distribution N(0,σ 2 ). By fitting this model, the research team successfully eliminated the influence of batch effects and obtained the corrected gene expression matrix. As Figure 2 (Comparison of batch effect correction effects) shows, it demonstrates the effect of batch effect correction in data preprocessing. The left figure shows the gene expression distributions of two different batches before correction, and obvious batch differences can be seen; the right figure shows the distribution after correction, and the expression distributions of the two batches basically overlap, indicating that the batch effects have been effectively eliminated.

[0240] Next, for the missing values in the original data, the research team used the method of iterative singular value decomposition for filling. Specifically, they decomposed the expression matrix X into the left singular matrix U, the singular value diagonal matrix ∑, and the right singular matrix V T , as follows:

[0241] X = U∑V T + E;

[0242] where E is the residual matrix. Through iterative optimization, the research team finally obtained the complete filled expression matrix.

[0243] To mitigate the highly skewed distribution of the expression data, the research team performed logarithmic transformation on the data:

[0244] X transformed = log2(X raw + 0.5);

[0245] At the same time, they also used the Z-score standardization method to normalize the data:

[0246]

[0247] where Z ij is the standardized expression value; X ij is the original expression value; μ i is the mean of the i-th feature; σ i is the standard deviation of the i-th feature.

[0248] Finally, the research team used the method of principal component analysis to decompose the total expression matrix X into the expression stable component X stable and the expression variable component X variable :

[0249] X = X stable + X variable + ∈;

[0250] And by calculating the Pearson correlation coefficient between samples, the sample correlation characteristics of these two types of components were analyzed.

[0251] In the sample quality assessment stage, the research team first calculated the Pearson correlation coefficient matrix R between samples:

[0252]

[0253] where R ij is the correlation coefficient between samples i and j; cov(X i , X j ) is the covariance; σ i , σ j are the standard deviations. Based on the correlation coefficient matrix, the research team used principal component analysis to extract the principal component score matrix T:

[0254] T = XW;

[0255] where X is the centered data matrix; W is the eigenvector matrix. According to the principal component scores, they calculated the quality score Q of each sample i :

[0256]

[0257] where w j is the weight of the j-th principal component; PC ij is the score of the i-th sample on the j-th principal component. The research team set the quality score threshold at 0.6 and excluded 4 abnormal samples. As Figure 3 (Principal component analysis scatter plot) shows, it presents the distribution of samples in the AD patient group and the control group in the space of the first two principal components. The red dots in the figure represent AD patient samples, the blue dots represent control group samples, and several key genes that contribute the most to the sample distribution are labeled. This

[0258] In the differential gene analysis stage, the research team first used three methods, DESeq2, limma, and SAM, to perform differential analysis on the gene expression data of the AD patient group and the healthy control group respectively. DESeq2 uses a negative binomial distribution model to fit the gene expression count data:

[0259] K ij ~NB(μ ij , α i );

[0260] log2(μ ij ) = β i0 +β i1 x j +log2(s j );

[0261] where K ij is the count value of the i-th gene in the j-th sample; μ ijis the mean parameter; α i is the dispersion parameter; β i0 is the baseline expression level; β i1 is the treatment effect; s j is the scale factor.

[0262] Limma then uses an empirical Bayes model to fit gene expression data:

[0263] E[y g = Xα g ;

[0264]

[0265] where, y g is the expression vector of gene g; X is the design matrix; α g is the coefficient vector; is the variance parameter; w gi is the weight.

[0266] The SAM method then evaluates gene differences based on the t-test statistic:

[0267]

[0268] where, d i is the statistic of gene i; is the mean of the AD group; is the mean of the control group; s i is the standard deviation; s0 is the adjustment factor.

[0269] Through meta-analysis, the research team integrated the results of these three methods to obtain a relatively reliable set of differentially expressed genes. Next, they input these differentially expressed genes into a pre-trained deep learning model to calculate the importance score S i :

[0270] S i = β1||θ i ||1 + β2∑ j∈N(i) w ij S j + β3f(X i ) + ∈ i ;

[0271] where, θ i is the L1 regularization coefficient vector; w ij is the network edge weight; N(i) is the neighbor set of gene i; f(X i ) is the expression feature function; β1, β2, β3 are weight parameters; ∈ iis the error term. According to the fold change (|log2FC|>1) and FDR significance (FDR<0.05) thresholds, combined with the importance scores of differentially expressed genes, the research team finally screened out 471 differentially expressed genes. As Figure 3 (Heatmap of differentially expressed genes) shows, it presents the expression patterns of 20 important differentially expressed genes in AD patient and control samples. Red indicates high expression and blue indicates low expression. The heatmap clearly shows the expression differences of these genes in the two groups of samples.

[0272] In the network construction stage, the research team first constructed a gene co-expression correlation network based on these 471 differentially expressed genes. Specifically, they calculated the Pearson correlation coefficient matrix R between genes, and according to the criterion that the absolute value of the correlation coefficient is greater than 0.7 and the stability coefficient S ij is greater than 0.6, 709 significantly correlated edges were selected:

[0273]

[0274] where, N is the number of resampling times; is the correlation coefficient of the k-th sampling; I() is the indicator function.

[0275] In addition, the research team also integrated transcription factor binding site data, protein-protein interaction data, and metabolite-protein interaction data, and calculated the confidence coefficient C of each type of data respectively k :

[0276] C k = 0.4Q k + 0.3P k + 0.3P k ;

[0277] where, Q k is the quality score; P k is the integrity score; R k is the repeatability score. According to the confidence of each type of data, they constructed a comprehensive weighted regulatory network, and its edge weight M ij is calculated as follows:

[0278]

[0279] where, is the element of the adjacency matrix of the k-th layer network.

[0280] In the node importance evaluation stage, the research team first extracted the degree centrality DC i , betweenness centrality BC i , closeness centrality CC i and eigenvector centrality EC iTopological features such as:

[0281]

[0282]

[0283] where A ij is an element of the adjacency matrix; σ st is the number of shortest paths from node s to t; σ st (i) is the number of shortest paths passing through node i; d ij is the shortest distance from node i to j; n is the total number of nodes; λ is the eigenvalue. Then, they used a graph neural network model to perform representation learning on the network and calculated the temporal contribution score of each node in the network module based on this. Finally, the research team trained a multi-layer neural network model, comprehensively considering network topological features, module contribution, biological importance, and disease relevance, to obtain the comprehensive score Score i :

[0284] Score i = 0.4T i + 0.3M i + 0.2B i + 0.1D i ;

[0285] where T i is the contribution of the network topological feature of node i; M i is the temporal contribution of node i in the module; B i is the biological importance of node i; D i is the correlation between node i and AD. Based on this, the research team identified the top 100 key nodes with the highest comprehensive scores.

[0286] To further verify the biological significance of these key nodes, the research team conducted the following stratified analyses respectively:

[0287] 1) Pathway enrichment analysis: The research team mapped the genes corresponding to these 100 key nodes to the KEGG database and found that they were mainly involved in pathways closely related to the pathogenesis of AD, such as neurotransmitter regulation, immune response, and apoptosis. This indicates that these key nodes play an important regulatory role in the AD pathogenesis process.

[0288] 2) Disease gene overlap analysis: By querying the AD-related gene database, the research team found that 38 key node genes have been identified as pathogenic genes or susceptibility genes for AD. This further proves the core position of these nodes in the AD pathogenesis.

[0289] 3) Expression stability analysis: The research team examined the expression changes of these 100 key nodes in different AD patient samples and found that most of them showed high expression stability, which is beneficial for their use as reliable biomarkers.

[0290] 4) Temporal change analysis: By integrating gene expression data at different stages of AD progression, the research team found that the expression changes of these key nodes often reflect the development process of AD. Some nodes showed significant upregulation in the early stage, while others showed downregulation in the late stage, providing a basis for using these nodes to monitor the dynamic changes of AD.

[0291] 5) Network perturbation analysis: By simulating the removal of these key nodes, the research team found that their deletion significantly changed the topological structure of the entire AD-related network, thereby affecting network function. This indicates that these nodes occupy a key hub position in the AD regulatory network.

[0292] Based on the above analysis results, the research team finally confirmed 50 high-confidence AD key biomarker nodes. These nodes cover multiple key biological processes closely related to the pathogenesis of AD, such as neurotransmitter regulation, immune response, and apoptosis, providing valuable clues for further understanding the pathogenesis of AD and developing new diagnostic and treatment strategies.

[0293] The variables involved in the present invention are explained as shown in Table 1 below.

[0294] Table 1 Variable Explanation Table

[0295]

[0296]

[0297] The above is only the specific implementation manner of the present invention, but the protection scope of the present invention is not limited thereto. Any person skilled in the art can easily think of changes or substitutions within the technical scope disclosed by the present invention, and all should be covered by the protection scope of the present invention.

Claims

1. A method for identifying key nodes of a disease biomarker expression regulation network, characterized in that, Including the following steps: Perform data preprocessing on the original omics data related to disease tissue and biomarker samples to obtain expression stable components and expression variable components; Calculate the inter-sample correlation of the expression components, construct a sample quality score matrix, identify and remove abnormal samples to obtain quality-controlled data; Identify genes with significant changes in the disease state, obtain the importance scores of differential genes, and screen the differentially expressed gene set; construct an expression correlation coefficient matrix, screen significantly correlated gene pairs, and construct a gene expression correlation network; Integrate multiple types of data, calculate the confidence coefficient, and construct a multi-level weighted regulatory network; construct a network node feature matrix and calculate the contribution degree of network topological features; Use a graph neural network model for representation learning, perform module division, and calculate the node time series contribution degree score; Integrate multiple features to construct a node comprehensive scoring model and identify key nodes; Perform hierarchical verification on the key nodes, optimize the model, and obtain the final set of verified key nodes.

2. The key node identification method for a disease biomarker expression regulation network according to claim 1, characterized in that The original omics data is multi-omics integrated data, including gene expression profile data, protein expression profile data, metabolomics data, epigenomic data, and clinical phenotype data.

3. The method for identifying key nodes of a disease marker expression regulation network according to claim 2, wherein, The inter-sample correlation of the expression stable components and expression variable components is calculated using the Pearson correlation coefficient to measure the degree of expression correlation between samples.

4. The key node identification method for a disease biomarker expression regulation network according to claim 3, wherein Identifying genes with significant changes in the disease state combines three differential analysis methods, DESeq2, limma, and SAM, and integrates the results of multiple methods through Meta-analysis.

5. The key node identification method for a disease biomarker expression regulation network according to claim 4, wherein The differential gene importance score model adopts an embedded deep learning structure, including an embedded mathematical model layer, a feature extraction layer, an attention mechanism layer, a deep neural network layer, and an output layer.

6. The key node identification method for a disease biomarker expression regulation network according to claim 5, wherein, The embedded mathematical model layer is a system of mathematical equations, including a feature selection equation, a sample weight equation, and a network structure equation; the feature selection equation is used for feature screening based on L1 regularization; the sample weight equation is used for sample weighting and integration; the network structure equation is used for network structure feature extraction; the feature extraction layer is used for data dimensionality reduction and feature learning; the attention mechanism layer is used for feature importance evaluation; the deep neural network layer is used for non-linear feature transformation; the output layer is used to output the differential gene importance score.

7. The key node recognition method for a disease biomarker expression regulation network according to claim 6, characterized in that Screening the differentially expressed gene set selects differentially expressed genes according to the criteria that the absolute value of the fold change is greater than 1 and the significance is less than 0.05, combined with the importance score ranking.

8. The method for identifying key nodes of a disease biomarker expression regulation network according to claim 7, characterized in that Screening significantly correlated gene pairs sets gene pairs with an absolute value of the correlation coefficient greater than 0.7 and a stability coefficient greater than 0.6 as significantly correlated gene pairs.

9. The key node identification method for a disease marker expression regulation network according to claim 8, characterized in that Constructing a gene expression correlation network constructs an adjacency matrix based on the screened significantly correlated gene pairs and assigns the correlation coefficient as the edge weight.

10. The method for identifying key nodes of a disease biomarker expression regulatory network according to claim 9, characterized in that, Calculating the confidence coefficient for each type of data is based on data quality scores, data integrity, and experimental repeatability to calculate a comprehensive confidence coefficient; constructing a multi-level weighted regulatory network integrates different types of regulatory relationships into a heterogeneous network according to the confidence coefficient; performing representation learning uses a graph convolutional neural network for node feature encoding and structural feature learning; Module division is to use the Louvain algorithm for community detection and evaluate the module stability.

Images