A pathway enrichment analysis method and apparatus based on Moran's I
By constructing a spatial weight matrix based on the Moran's index and utilizing Moran's I statistic, the problems of topological structure and multi-omics data integration in pathway enrichment analysis were solved, key regulatory nodes were identified, and the accuracy and relevance of the analysis were improved.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- TSINGHUA UNIVERSITY
- Filing Date
- 2026-02-12
- Publication Date
- 2026-06-05
AI Technical Summary
Existing pathway enrichment analysis methods ignore pathway topology, cannot integrate multi-omics data, fail to properly handle missing data, and are difficult to identify key regulatory nodes, resulting in insufficient accuracy and relevance of the analysis results.
Using a masked Moran index-based approach, a spatial weight matrix is constructed, and global and local Moran's I statistics are utilized in conjunction with Monte Carlo permutation tests to identify the spatial autocorrelation and key regulatory nodes of pathways.
It effectively integrates pathway topology and multi-omics data, identifies key regulatory nodes, improves the accuracy and relevance of the analysis, and provides more precise information for subsequent experimental verification.
Smart Images

Figure CN122157773A_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the fields of bioinformatics and systems biology, and in particular to a pathway enrichment analysis method and apparatus based on Moran's I. Background Technology
[0002] In biomedical research, high-throughput omics technologies (such as proteomics and metabolomics) can detect changes in molecular expression in biological samples at a global level. To extract biological significance from a large number of differentially expressed molecules, researchers typically use pathway enrichment analysis to map differentially expressed molecules to known biological pathways and assess which pathways have changed significantly under specific biological conditions.
[0003] Currently, commonly used pathway enrichment analysis methods mainly include over-representation analysis (ORA) based on hypergeometric tests and gene set enrichment analysis (GSEA) based on ordination. However, these methods have the following limitations:
[0004] 1. Ignoring pathway topology: Traditional methods treat pathways as simple collections of genes / proteins, neglecting the interactions between molecules within the pathway and upstream and downstream regulatory information. Biological pathways are actually complex regulatory networks with various relationships between molecules, such as activation and inhibition. This topological information is crucial for understanding the functional state of pathways.
[0005] 2. Inability to integrate multi-omics data: Proteins and metabolites are closely linked in cell signaling, for example, enzymes (proteins) catalyze the conversion of substrates (metabolites) into products (metabolites). Traditional GO enrichment analysis is only applicable to genes / proteins, and although KEGG enrichment analysis includes metabolic pathways, in practical applications, proteins and metabolites are often analyzed independently, which cannot reflect their synergistic regulatory patterns in the same pathway.
[0006] 3. Improper handling of missing data: In actual experiments, due to limitations in detection technology, only some molecules in a pathway are often identified. Traditional methods typically ignore undetected molecules, which may lead to the loss of pathway connectivity information and affect the accuracy of the analytical results.
[0007] 4. Inability to identify key regulatory nodes: Traditional methods mainly focus on the overall enrichment of the entire pathway or gene set, making it difficult to further identify which specific proteins or metabolites in the pathway are the core nodes that play a key regulatory role. This limits the specificity of subsequent experimental verification.
[0008] Therefore, there is an urgent need to develop a new pathway enrichment analysis method that can comprehensively consider the network topology of the pathway, integrate multi-omics data of proteins and metabolites, properly handle missing data, and identify key regulatory nodes in the pathway, so as to provide more accurate information for biological mechanism research and experimental verification. Summary of the Invention
[0009] The present invention aims to at least partially solve one of the technical problems in the related art.
[0010] Therefore, the first objective of this invention is to propose a path enrichment analysis method based on the mask Moran index, comprising: S1, acquire proteomics and / or metabolomics data of biological samples, and extract molecular identifiers and their corresponding expression values; S2, map the molecular identifier to the nodes of the pathway network in the preset pathway database, and construct a spatial weight matrix based on the topology of the pathway network. In the pathway network, the nodes represent biomolecules, the edges between the nodes represent the regulatory or interaction relationships between biomolecules, and the spatial weight matrix represents the proximity relationship between the nodes in the pathway network. S3, masking is performed on nodes that are not mapped to the expressed data, and the observed nodes and missing data nodes are marked by the mask vector; S4, extract sub-graphs based on mask vectors, and calculate the global spatial autocorrelation of pathways using the global Moran's I statistic; S5. The significance was assessed by the Monte Carlo permutation test, and the p-value corresponding to the global Moran's I statistic was calculated. S6. Repeat steps S2-S5 to obtain the global Moran's I statistic and its corresponding p-value for all pathways. Sort the pathways according to the p-value and select pathways with p-values less than the significance threshold as enriched pathways.
[0011] In one embodiment of the present invention, S2 further includes: S21, Obtain pathway network data from a preset pathway database and convert molecular identifiers of proteins and metabolites; When the molecular identifier is a protein identifier, use a gene annotation database to convert the identifier type and convert the input protein identifier into a standard identifier used by the pathway database; When the molecular identifier is a metabolite identifier, use a metabolite mapping table to convert the identifier type and the metabolite mapping table contains the correspondence between metabolite identifiers in different databases. S22, Match the converted identifier with the identifier of the path network node to establish the correspondence between the expressed data and the network node; S23, determine whether the transformed identifier matches the node successfully. For nodes that match successfully, assign the corresponding expression difference value to the node; for nodes that do not match successfully, mark them as missing data nodes. S24, obtain the adjacency matrix of the path network, the adjacency matrix indicates whether there is a direct connection relationship between each node, and calculate the spatial weight between nodes based on the network topology distance; The method for calculating the spatial weights between nodes based on network topology distance is as follows: Calculate the spatial weights between nodes using the shortest path distance weighting method, the inverse distance weighting method, or the Gaussian kernel weighting method; Before calculating the spatial weights, the spatial weight matrix is self-connected.
[0012] In one embodiment of the present invention, S4 further includes: S41, using the observation nodes in the mask vector, extract the expression values of the observation nodes from the complete expression value vector to form expression value sub-vectors, extract the weight relationships between the observation nodes from the complete spatial weight matrix to form observation node weight sub-matrix, and use the observation nodes, expression value sub-vectors and observation node weight sub-matrix to form a subgraph; S42, based on the extracted subgraphs, calculate the sum of subgraph weights and the sum of squares of subgraph representation values, and determine the subgraph degradation based on the sum of subgraph weights and the sum of squares of subgraph representation values. When the sum of subgraph weights is greater than 0, calculate Moran's I statistic; otherwise, output an invalid value. S43. Based on the results of Moran's I statistic, the global spatial autocorrelation of the pathway is determined. A Moran's I value close to 1 indicates strong positive spatial autocorrelation, where adjacent nodes in the pathway tend to be upregulated or downregulated simultaneously. A Moran's I value close to -1 indicates strong negative spatial autocorrelation, where adjacent nodes in the pathway tend to have opposite expression changes. A Moran's I value close to 0 indicates no spatial autocorrelation, where expression changes are randomly distributed in the pathway.
[0013] In one embodiment of the present invention, the formula for calculating the Moran's I statistic is as follows: ; Where, n eff To observe the number of rows or columns in the node weight submatrix, This represents the sum of the subgraph weights. The subgraph represents the sum of squared values. To express the value sub-vector, This is the weight submatrix for the observation nodes.
[0014] In one embodiment of the present invention, S5 further includes: S51, construct a permutation vector of length n, and randomly sample n from all expression values. eff The extracted expression values are then used to fill the first n values of the permutation vector. eff There are n positions, and the remaining positions are filled with 0, where n eff To observe the number of rows or columns in the weight submatrix of the nodes; S52, randomly sort the elements in the permutation vector to obtain the reordered permutation vector, and construct the permutation mask based on the non-zero positions of the reordered permutation vector. S53, calculate the permutation mask statistic. Based on the positive and negative values of the permutation mask statistic and the relationship between the permutation mask statistic and the global Moran's I statistic, perform extreme value statistics. When the permutation mask statistic is ≥0 and the global Moran's I statistic is ≥ the permutation mask statistic, increment the first counter by one. When the permutation mask statistic is <0 and the global Moran's I statistic is ≤ the permutation mask statistic, increment the second counter by one. S54, calculate the one-sided p-value based on the sign of the permutation mask statistic. When the permutation mask statistic is ≥ 0, the formula for calculating the p-value is: p = (ge + 1) / (N + 1); Where ge is the count of the first counter, and N is the set total number of replacements; When the permutation mask statistic is less than 0, the formula for calculating the p-value is: p = (le + 1) / (N + 1); Where le is the number of times the second counter counts; When the one-sided p-value is less than the preset significance threshold, the spatial autocorrelation of the pathway is considered significant.
[0015] To achieve the above objectives, a second aspect of the present invention proposes a method for analyzing key regulatory nodes in enriched pathways based on local Moran's I, comprising: A1. Obtain proteomics and / or metabolomics data from biological samples, and extract molecular identifiers and their corresponding expression values. A2, map the molecular identifier to the nodes of the pathway network in the preset pathway database, and construct a spatial weight matrix based on the topology of the pathway network. In the pathway network, the nodes represent biomolecules, the edges between the nodes represent the regulatory or interaction relationships between biomolecules, and the spatial weight matrix represents the proximity relationship between the nodes in the pathway network. A3, performs masking on nodes that are not mapped to the expressed data, marks the observation nodes and missing data nodes with mask vectors, and uses the observation nodes as candidate key nodes; A4. Calculate the local Moran's I statistic for each key node and evaluate the significance of the local Moran's I for each candidate key node through permutation test. A5. Based on the expression value, local Moran's I statistic, and significance of candidate key nodes, determine whether they are key nodes and classify them. The key node classification includes high-high clustering pattern, low-low clustering pattern, high-low clustering pattern, and low-high clustering pattern.
[0016] In one embodiment of the present invention, the formula for calculating the local Moran's I statistic for each observation node is as follows: ; Where, n eff S0 is the number of observation nodes; S0 is the spatial weight matrix W. sub The sum of weights; z i Let z be the centered expression variance value of node i, which is the expression variance value of that node minus the mean of the expression values of the observed nodes; m T This represents the transpose of a vector.
[0017] In one embodiment of the present invention, A4 further includes: A41, construct a permutation vector of length n, and randomly sample n values from all representation values. eff The extracted expression values are then used to fill the first n values of the permutation vector. eff There are n positions, and the remaining positions are filled with 0, where n eff To observe the number of rows or columns in the weight submatrix of the nodes; A42, randomly sort the elements in the permutation vector to obtain the reordered permutation vector, and construct the permutation mask based on the non-zero positions of the reordered permutation vector. A43. Calculate the local Moran's I statistic and local permutation mask statistic for each node. When the local Moran's I statistic of a node is a finite value, increment the local first counter by one. Perform extreme value statistics based on the positive and negative values of the local permutation mask statistic and the relationship between the local permutation mask statistic and the local Moran's I statistic. When the local permutation mask statistic is ≥0 and the local Moran's I statistic is ≥ the local permutation mask statistic, increment the local second counter by one. When the local permutation mask statistic is <0 and the local Moran's I statistic is ≤ the local permutation mask statistic, increment the local third counter by one. A44, calculate one-sided p based on the sign of the local permutation mask statistic. local When the local permutation mask statistic is ≥0, p local The formula for calculating the value is: plocal [i] = (cnt ge [i] + 1) / (cnt included [i] + 1); Among them, cnt included The local first counter counts the number of times, cnt ge This is the number of times the local second counter counts. When the local permutation mask statistic is < 0, p local The formula for calculating the value is: p local [i] = (cnt le [i] + 1) / (cnt included [i] + 1); Among them, cnt le This is the number of times the local third counter counts. When p local When the value is less than the preset significance threshold, the node is considered significant.
[0018] In one embodiment of the present invention, the method for determining the classification of key nodes in step A5 is as follows: When the expression value of the key node is higher than the mean, the local Moran's I statistic is greater than 0, and the key node is significant, the key node is a high-high clustering pattern node. When the expression value of the key node is lower than the mean, the local Moran's I statistic is greater than 0, and the key node is significant, the key node is a low-low clustering pattern node. When the expression value of the key node is higher than the mean, the local Moran's I statistic is less than 0, and the key node is significant, the key node is a high-low clustering pattern node. When the expression value of a key node is lower than the mean, the local Moran's I statistic is less than 0, and the key node is significant, the key node is a low-high clustering pattern node.
[0019] To achieve the above objectives, a third aspect of the present invention provides a pathway enrichment analysis device based on the mask Moran index, comprising: Data acquisition module: used to acquire omics data from biological samples and determine each molecular identifier and its corresponding expression change value; Identifier mapping module: used to map the molecular identifier to nodes of the pathway network in the preset pathway database; Weight matrix construction module: used to construct a spatial weight matrix based on the topology of the pathway network; Masking module: used to mask nodes in the path network that are not mapped to the expressed data; The statistics calculation module is used to calculate the global spatial autocorrelation of the pathway using the masked Moran's I statistic. Significance testing module: used to assess the significance level of the global spatial autocorrelation using Monte Carlo permutation tests; Local Analysis Module: Used to calculate local Moran's I statistics and identify key nodes in enriched pathways; Results output module: used to determine the list of enriched pathways based on the significance level.
[0020] The methods, systems, and storage media of this invention take into account pathway topology, integrate multi-omics data, employ innovative missing data processing, identify key regulatory nodes, and are easy to use and promote. By calculating local Moran's I statistics, this invention can identify key nodes (hotspots) in pathway networks. These nodes and their neighboring nodes exhibit significant co-expression change patterns, providing more accurate and targeted candidate targets for subsequent experimental verification.
[0021] Additional aspects and advantages of the invention will be set forth in part in the description which follows, and in part will be obvious from the description, or may be learned by practice of the invention. Attached Figure Description
[0022] The above and / or additional aspects and advantages of the present invention will become apparent and readily understood from the following description of the embodiments taken in conjunction with the accompanying drawings, wherein: Figure 1 A schematic flowchart of the pathway enrichment analysis method based on Moran's I provided in an embodiment of the present invention; Figure 2 A schematic diagram of the method for analyzing key regulatory nodes of enrichment pathways based on local Moran's I provided in an embodiment of the present invention; Figure 3 A schematic diagram illustrating the calculation of Moran's I statistic for a mask provided in an embodiment of the present invention; Figure 4 A comparative schematic diagram of the global and local Moran's I analysis flow provided in the embodiments of the present invention; Figure 5 This is a schematic diagram of the pathway enrichment analysis device based on Moran's I provided in an embodiment of the present invention. Detailed Implementation
[0023] It should be noted that, unless otherwise specified, the embodiments and features described in the present invention can be combined with each other. The present invention will now be described in detail with reference to the accompanying drawings and embodiments.
[0024] To enable those skilled in the art to better understand the present invention, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort should fall within the scope of protection of the present invention.
[0025] The following description, with reference to the accompanying drawings, describes a method and apparatus for path enrichment analysis based on the mask Moran index, according to an embodiment of the present invention.
[0026] Example 1 Figure 1 This is a flowchart of a pathway enrichment analysis method based on Moran's I according to an embodiment of the present invention.
[0027] like Figure 1 As shown, the pathway enrichment analysis method based on Moran's I includes the following steps: S1, acquire proteomics and / or metabolomics data of biological samples, and extract molecular identifiers and their corresponding expression change values.
[0028] Specifically, in some implementations, acquiring proteomics and / or metabolomics data from biological samples and extracting molecular identifiers and their corresponding expression change values is the initial key step of the method of this invention. Its technical implementation is based on a standardized processing flow for high-throughput omics data. The core of this step lies in extracting biologically significant molecular identifiers (such as UNIPROT, ENTREZID, ENSEMBL, or SYMBOL for proteins, and HMDB, KEGG, ChEBI, or CID for metabolites) and their expression differential values (logFC or normalized expression levels) from the raw experimental data, providing basic input for subsequent pathway mapping and spatial autocorrelation analysis.
[0029] As one implementation method, for proteomics data, protein identifiers (such as UNIPROT ID, ENTREZID, ENSEMBL ID, or gene SYMBOL) and their corresponding expression differential values are extracted. Expression differential values are typically log2 fold change (logFC), representing the expression change in the experimental group relative to the control group.
[0030] For example, proteomics data can be represented as: Protein_IDlogFC P046371.25 P00533-0.85 Q9Y6K90.67 ... For metabolomics data, extract metabolite identifiers (such as HMDB ID, KEGG ID, ChEBI ID, or CID) and their corresponding differential expression values.
[0031] For example, metabolomics data can be represented as: Metabolite_IDlogFC HMDB00001220.52 HMDB0000256-0.38 HMDB00000941.15 ... In this embodiment, logFC is preferably used to express the difference value because it has symmetry (the numerical range of upregulation and downregulation is consistent) and comparability (logFC or standardized logFC values under different experimental conditions can be directly compared).
[0032] At the technical implementation level, this step is typically accomplished through a data parsing module, supporting various input file formats such as CSV, Excel, or TSV. The module first identifies the identifier and expression value columns in the data, using regular expressions to match standard identifier formats, ensuring data accuracy and consistency. The expression difference value is typically log2 Fold Change (logFC), calculated based on the ratio of expression levels between the experimental and control groups. It possesses symmetry and comparability, making it suitable for spatial statistical analysis.
[0033] In practical applications, this step is widely applicable to multi-omics integrated analysis in biomedical research, such as tumor marker screening, drug target discovery, and the elucidation of metabolic disorder mechanisms. By extracting molecular identifiers and expression change values, researchers can precisely match experimental data with pathway networks, laying the foundation for subsequent Moran's I statistic calculations.
[0034] The technical advantage of this step lies in providing structured and standardized input data for pathway enrichment analysis, ensuring the accuracy of subsequent spatial weight matrix construction and statistical testing. Through reasonable screening and mapping, this invention effectively improves the biological interpretability and computational stability of pathway analysis.
[0035] S2, the molecular identifier is mapped to the node of the pathway network in the preset pathway database, and a spatial weight matrix is constructed based on the topology of the pathway network. The spatial weight matrix quantifies the proximity relationship by calculating the shortest path distance between nodes and using inverse distance weight or Gaussian kernel weight.
[0036] Further, step S2 includes: S21, Obtain pathway network data from a preset pathway database and convert molecular identifiers of proteins and metabolites; When the molecular identifier is a protein identifier, use a gene annotation database to convert the identifier type and convert the input protein identifier into a standard identifier used by the pathway database; When the molecular identifier is a metabolite identifier, use a metabolite mapping table to convert the identifier type and the metabolite mapping table contains the correspondence between metabolite identifiers in different databases. S22, Match the converted identifier with the identifier of the path network node to establish the correspondence between the expressed data and the network node; S23, determine whether the transformed identifier matches the node successfully. For nodes that match successfully, assign the corresponding expression difference value to the node; for nodes that do not match successfully, mark them as missing data nodes. S24, obtain the adjacency matrix of the path network, the adjacency matrix indicates whether there is a direct connection between each node, and calculate the spatial weight between nodes based on the network topology distance.
[0037] In actual omics experiments, due to limitations in detection technology and differences in molecular abundance in samples, only a small number of nodes in a pathway can often be successfully identified and quantified. For example, in a metabolic pathway containing 50 nodes, only 5-10 metabolites may be identified; in a signaling pathway containing 100 nodes, only 20-30 proteins may be detected.
[0038] If the adjacency matrix (i.e., only directly adjacent nodes have a weight of 1, and others have a weight of 0) is used directly as the spatial weight matrix, the following problems will occur: 1. Network connectivity disruption: When the observed nodes are sparsely distributed in the path, after extracting the subgraph of the observed nodes, the originally connected path may become multiple isolated nodes or small connected components, resulting in the total weight S0 of the spatial weight matrix being extremely small or even 0, making it impossible to calculate Moran's I statistic.
[0039] 2. Loss of long-range regulatory information: Biological pathways contain multi-step regulatory relationships. For example, upstream signaling molecules influence downstream effector molecules through multiple intermediate nodes. If only directly adjacent nodes are considered, these important long-range regulatory patterns will be overlooked.
[0040] Therefore, this invention employs a distance-based weight calculation method, which ensures that all node pairs have non-zero weights, and the weight values decay as the distance increases, thereby maintaining the global connectivity of the network.
[0041] As one implementation method, the weight calculation method used in this invention is as follows: First, graph theory algorithms are used to compute the shortest path distance matrix D between any two nodes in the path network, where D... ij This represents the length of the shortest path from node i to node j (the number of edges traversed). The shortest path distance can be calculated using Dijkstra's algorithm, Floyd-Warshall's algorithm, or breadth-first search (BFS) algorithm.
[0042] Method 1: Inverse Distance Weighting (Power-Law Decay) Inverse distance weights are calculated using a power-law function to determine the weights between nodes.
[0043] Among them, D ij Let be the shortest path distance between node i and node j; p is the distance decay parameter, controlling the rate at which the weight decays with distance. The larger p is, the faster the weight decays, and the smaller the influence of distant nodes. When p=1, it is a simple inverse distance weight, with nodes 2 steps apart having a weight of 1 / 2 and nodes 3 steps apart having a weight of 1 / 3. When p=2, it is a squared inverse distance weight, with even faster weight decay.
[0044] Furthermore, to accommodate most biological pathway analysis scenarios, it is recommended to use p=1 as the default value.
[0045] Method 2: Gaussian kernel weights (exponential decay) Gaussian kernel weights use an exponential function to calculate the weights between nodes:
[0046] Among them, D ij Let be the shortest path distance between node i and node j; σ is the bandwidth parameter, controlling the rate of weight decay. The larger σ is, the slower the weight decays, and the greater the influence of distant nodes. Gaussian kernel weights have the characteristic of smooth decay, with the weight approaching 1 when the distance is small, and decaying exponentially with increasing distance. When σ equals the average path length of the network, a good balance can be achieved between maintaining local proximity and global connectivity.
[0047] Furthermore, it is recommended that σ be set to 1 / 2 to 1 times the average path length of the pathway network.
[0048] To further improve the processing performance, post-processing of the weight matrix is performed: 1. Self-join processing: Set the diagonal elements of the weight matrix to 0, i.e., W[i,i] = 0, because nodes do not have spatial dependencies on themselves. When calculating the shortest path distance, first set the diagonal of the distance matrix to infinity (Inf), and then replace the infinity weights with 0 after calculating the weights.
[0049] 2. Symmetry treatment: To satisfy the symmetry assumptions of certain spatial statistical tests, the weight matrix can be symmetricized:
[0050] Among them, W T This represents the transpose of the weight matrix. After symmetry, the weight of node i with respect to node j is equal to the weight of node j with respect to node i, i.e., W. ij = W ji In embodiments of the invention, symmetry is recommended because a symmetric weight matrix ensures a more stable permutation distribution of Moran's I statistic and simplifies certain matrix operations. 3. Standardization: Perform row standardization on the weight matrix so that the sum of the elements in each row equals 1.
[0051] After standardization, the weight matrix represents the relative influence weights of neighboring nodes on the central node.
[0052] Specifically, in some implementations, mapping the molecular identifiers to nodes in a pathway network within a predefined pathway database and constructing a spatial weight matrix based on the topology of the pathway network is a key preliminary step for spatial autocorrelation analysis in the method of this invention. This step provides structured input for subsequent Moran's I statistic calculation by precisely matching molecular identifiers in the omics data with nodes in the pathway network.
[0053] At the technical implementation level, this step first extracts pathway network data from pre-defined pathway databases (such as KEGG, Reactome, and WikiPathways), where each pathway consists of nodes (representing biomolecules) and edges (representing regulatory or interaction relationships between molecules). For protein identifiers, if their type is inconsistent with the standard identifiers used in the pathway database (such as ENTREZID), identifier conversion is performed using a gene annotation database (such as org.Hs.eg.db); for metabolite identifiers, if their type is HMDB, ChEBI, etc., cross-database identifier conversion is performed using a metabolite mapping table. After mapping, successfully matched nodes are assigned their corresponding expression differential values (such as logFC), while unmatched nodes are marked as missing data nodes. This process ensures the structural consistency between omics data and pathway networks, laying the foundation for the subsequent construction of spatial weight matrices.
[0054] At the application level, this step is widely used in proteomics, metabolomics, and multi-omics integrated analysis. For example, in the KEGG metabolic pathway, if only some metabolites are detected, constructing a distance-based weight matrix can preserve the global structural information of the pathway and avoid network connectivity disruption caused by sparse nodes. In signaling pathways, if only some proteins are detected, this method can still quantify their potential regulatory relationships through topological distance, thereby improving the biological interpretability of enrichment analysis.
[0055] In terms of technical effectiveness, by mapping molecular identifiers to pathway network nodes and constructing a spatial weight matrix, this invention effectively integrates omics data with pathway topology, overcoming the shortcomings of traditional methods that ignore intermolecular regulatory relationships. This step ensures the computability and stability of Moran's I statistic under sparse data conditions, providing a structured basis for subsequent identification of key regulatory nodes in pathways, and significantly improving the accuracy and biological relevance of pathway enrichment analysis.
[0056] S3 performs masking on nodes that are not mapped to the expressed data, using mask vectors to mark the observed nodes and missing data nodes.
[0057] Specifically, in the Moran's I-based pathway enrichment analysis method, masking nodes not mapped to expression data is a crucial step in ensuring the accuracy of spatial autocorrelation calculations and the validity of permutation tests. This step utilizes logical mask vectors... In the labeled pathway network, nodes with valid expressed data (observation nodes) and nodes without detected expressed information are marked as missing data nodes. Masking not only affects the calculation of global spatial autocorrelation but also determines the random allocation strategy of samples in permutation tests.
[0058] At the technical implementation level, masking primarily relies on the results of identifier mapping. For each node in the pathway network, if its corresponding molecular identifier has a matching differential expression value (e.g., logFC) in the input proteomics or metabolomics data, then the node is marked as an observation node, and its expression value is assigned to the network node; if no match is found, the node is marked as a missing data node. Mask vector It is a length of The logical vector, where This represents the total number of nodes in the pathway network. Represents a node To observe the nodes, Represents a node This is a missing node.
[0059] In global spatial autocorrelation calculations, the expression values of missing nodes are assigned to zero to avoid interference with the statistics. However, in the Monte Carlo permutation test, these nodes are completely excluded from the permutation pool, ensuring that only the expression values of observed nodes are redistributed during random rearrangement. Specifically, in each permutation, a random sample is taken from the expression value pool. Values ( (To observe the total number of nodes), and randomly assign them to the pathway network. Each node One position is filled with zeros, and the remaining positions are filled with zeros to form a permutation vector. Subsequently, based on Construct a new mask at non-zero positions Used to extract permutation subgraphs and calculate permutation statistics. .
[0060] The parameter settings for mask processing include the method for constructing the mask vector, the sampling strategy for permutation samples, and the logic for handling missing nodes. Mask vector The construction of the permutation matrix must ensure a one-to-one correspondence with the network nodes. Sampling of the permutation samples must employ sampling without replacement to maintain the independence and representativeness of the samples. During the permutation process, the exclusion mechanism of mask nodes ensures the spatial weight matrix... The construction is based solely on observation nodes, thus avoiding weight matrix degradation or computational instability caused by missing data.
[0061] This step has broad applicability in practical applications, especially for the integrated analysis of proteomics and metabolomics data. In multi-omics data, due to limitations in detection technology, only some nodes in a pathway often have expression information. Masking can effectively address the problem of data sparsity and ensure the biological significance of the analysis results. Through the masking mechanism, this invention maintains the integrity of the pathway topology while avoiding interference from invalid nodes on the statistics, improving the computational accuracy and reliability of significance assessment of global and local Moran's I statistics, and providing a solid data foundation for identifying key regulatory nodes and enriched pathways.
[0062] S4 extracts subgraphs based on mask vectors and calculates the global spatial autocorrelation of pathways using the global Moran's I statistic.
[0063] Further, step S4 includes: S41, using the observation nodes in the mask vector, extract the expression values of the observation nodes from the complete expression value vector to form expression value sub-vectors, extract the weight relationships between the observation nodes from the complete spatial weight matrix to form observation node weight sub-matrix, and use the observation nodes, expression value sub-vectors and observation node weight sub-matrix to form a subgraph; S42, based on the extracted subgraphs, calculate the sum of subgraph weights and the sum of squares of subgraph representation values, and determine the subgraph degradation based on the sum of subgraph weights and the sum of squares of subgraph representation values. When the sum of subgraph weights is greater than 0, calculate Moran's I statistic; otherwise, output an invalid value. S43. Based on the results of Moran's I statistic, the global spatial autocorrelation of the pathway is determined. A Moran's I value close to 1 indicates strong positive spatial autocorrelation, where adjacent nodes in the pathway tend to be upregulated or downregulated simultaneously. A Moran's I value close to -1 indicates strong negative spatial autocorrelation, where adjacent nodes in the pathway tend to have opposite expression changes. A Moran's I value close to 0 indicates no spatial autocorrelation, where expression changes are randomly distributed in the pathway.
[0064] Specifically, in some implementations, calculating the global spatial autocorrelation of pathways using the masked Moran's I statistic and evaluating its significance through a Monte Carlo permutation test is one of the core steps of the method in this invention, aiming to extract biologically meaningful enrichment patterns from the topology of the pathway network. This step is based on the principles of spatial statistics, treating changes in molecular expression in the pathway as "spatial attributes," and constructing a spatial weight matrix. This is used to quantify the proximity relationships between nodes, thereby assessing the spatial clustering of expression values in the pathway network.
[0065] At the technical implementation level, the observation node subgraph is first extracted from the pathway network. This subgraph consists of nodes with effective expression difference values (such as logFC), and its expression value vector is denoted as... The weight submatrix is denoted as Then, the total weight of the subgraph is calculated. and the sum of squares of the expression values .like or If the expression value is not variable, it indicates that the pathway subgraph is not connected or there is no variation in expression value, so Moran's I statistic cannot be calculated and an invalid value NA is returned.
[0066] Furthermore, the formula for calculating the global Moran's I statistic is as follows:
[0067] Where, n eff To observe the number of nodes, This represents the quadratic form of the expression value under spatial weights, used to measure the weighted sum of the covariances of the expression values of a node and its neighboring nodes. The value range of this statistic is [-1, 1], with positive values indicating positive spatial autocorrelation, negative values indicating negative spatial autocorrelation, and 0 indicating no spatial autocorrelation.
[0068] Figure 3This diagram illustrates subgraph extraction and Moran's I statistic calculation. In the pathway network, green nodes represent observation nodes where expressed data was detected, and gray nodes represent nodes with missing data. Only the subgraph consisting of green nodes is used in the calculation.
[0069] S5. The significance was assessed by the Monte Carlo permutation test, and the p-value corresponding to the global Moran's I statistic was calculated.
[0070] Further, step S5 includes: S51, construct a permutation vector of length n, and randomly sample n from all expression values. eff The extracted expression values are then used to fill the first n values of the permutation vector. eff There are n positions, and the remaining positions are filled with 0, where n eff To observe the number of rows or columns in the weight submatrix of the nodes; S52, randomly sort the elements in the permutation vector to obtain the reordered permutation vector, and construct the permutation mask based on the non-zero positions of the reordered permutation vector. S53, calculate the permutation mask statistic. Based on the positive and negative values of the permutation mask statistic and the relationship between the permutation mask statistic and the global Moran's I statistic, perform extreme value statistics. When the permutation mask statistic is ≥0 and the global Moran's I statistic is ≥ the permutation mask statistic, increment the first counter by one. When the permutation mask statistic is <0 and the global Moran's I statistic is ≤ the permutation mask statistic, increment the second counter by one. S54, calculate the one-sided p-value based on the sign of the permutation mask statistic. When the permutation mask statistic is ≥ 0, the formula for calculating the p-value is: p = (ge + 1) / (N + 1); Where ge is the count of the first counter, and N is the set total number of replacements; When the permutation mask statistic is less than 0, the formula for calculating the p-value is: p = (le + 1) / (N + 1); Where le is the number of times the second counter counts; When the one-sided p-value is less than the preset significance threshold, the spatial autocorrelation of the pathway is considered significant.
[0071] Subsequently, steps S2-S5 were repeated to obtain the global Moran's I statistic and its corresponding p-value for all pathways. The pathways were then sorted according to their p-values, and pathways with p-values less than the significance threshold were selected as enriched pathways.
[0072] This permutation strategy uses a dynamic masking mechanism to simulate the random distribution of expression values under missing data conditions, thereby effectively evaluating the significance of pathway enrichment.
[0073] In application scenarios, this step is suitable for proteomics, metabolomics, and multi-omics integrated analysis, especially when there is a large amount of missing data in the pathway, it can maintain network connectivity and accurately identify enriched pathways. By introducing a spatial weight matrix and a masking mechanism, this invention overcomes the shortcomings of traditional enrichment analysis methods that ignore topological structure and cannot handle sparse data, providing a more accurate statistical basis for identifying pathways with co-expression patterns.
[0074] The pathway enrichment analysis method based on topological proximity in this invention can comprehensively consider the spatial correlation between pathway topology and molecular expression changes, thereby improving the biological interpretability and regulatory relationship identification ability of pathway enrichment analysis.
[0075] Example 2 Figure 2 This is a flowchart of an embodiment of the present invention for analyzing key regulatory nodes of enriched pathways based on local Moran's I.
[0076] Global Moran's I assesses the overall spatial autocorrelation of the entire pathway, while Local Moran's I assesses the local spatial autocorrelation of each node and its neighboring nodes. Local Moran's I can identify "hot spots" and "cold spots" in the network, i.e., regions where expression changes are significantly clustered.
[0077] like Figure 2 As shown, the method for analyzing key regulatory nodes in enrichment pathways based on local Moran's I includes the following steps: A1. Obtain proteomics and / or metabolomics data from biological samples, and extract molecular identifiers and their corresponding expression values.
[0078] A2, the molecular identifier is mapped to the node of the pathway network in the preset pathway database, and a spatial weight matrix is constructed based on the topology of the pathway network. The node in the pathway network represents a biomolecule, the edge between each node represents the regulation or interaction between biomolecules, and the spatial weight matrix represents the proximity relationship between nodes in the pathway network.
[0079] A3 performs masking on nodes that are not mapped to the expressed data, and uses the mask vector to mark the observation nodes and missing data nodes, and uses the observation nodes as candidate key nodes.
[0080] A4. Calculate the local Moran's I statistic for each candidate key node, and evaluate the significance of the local Moran's I for each candidate key node through permutation test. Specifically, the formula for calculating the local Moran's I statistic for each observation node is as follows: ; Where, n eff S0 is the number of observation nodes; S0 is the spatial weight matrix W. sub The sum of weights; z i Let z be the centered expression variance value of node i, which is the expression variance value of that node minus the mean of the expression values of the observed nodes; m T This represents the transpose of a vector.
[0081] In the local Moran's I statistic, positive I i The value indicates that the expression changes of node i and its neighboring nodes are in the same direction; a negative value indicates a negative value. I i The value indicates that the expression change of node i is opposite to that of its neighboring nodes.
[0082] A5. Based on the expression value, local Moran's I statistic, and significance of candidate key nodes, determine whether they are key nodes and classify them. The key node classification includes high-high clustering pattern, low-low clustering pattern, high-low clustering pattern, and low-high clustering pattern.
[0083] Specifically, the method for classifying key nodes based on the expression value of node i, local Moran's I statistic, and significance is as follows: When the expression value of a key node is above the mean, the local Moran's I statistic is >0, and the key node is significant, it is considered a high-high aggregation pattern node. High-high aggregation pattern nodes and their neighboring nodes all exhibit high expression, representing core regions of co-upregulation in the pathway. Biologically, high-high aggregation pattern nodes may correspond to key proteins or metabolites that play a positive regulatory role.
[0084] When the expression level of a key node is below the mean, the local Moran's I statistic is >0, and the key node is significant, it is considered a low-low aggregation pattern node. Low-low aggregation pattern nodes and their neighboring nodes all exhibit low expression, representing core regions of co-downregulation within the pathway. Biologically, low-low aggregation pattern nodes may represent repressed pathway branches.
[0085] When the expression value of a key node is higher than the mean, the local Moran's I statistic is <0, and the key node is significant, the key node is a high-low aggregation pattern node. A high-low aggregation pattern node is highly expressed, but its neighboring nodes are lowly expressed. Biologically, this node may be a key point for regulatory mode switching or a core node of negative feedback regulation.
[0086] When the expression value of a key node is below the mean, the local Moran's I statistic is <0, and the key node is significant, the key node is a low-high aggregation pattern node. Low-high aggregation pattern nodes are low in expression but their neighboring nodes are high in expression. Biologically, this node may be a key target of inhibition or a pathway branch point.
[0087] By identifying these key nodes, researchers can gain a deeper understanding of the biological mechanisms of pathway enrichment, provide priority candidate targets for subsequent experimental validation, and discover potential drug targets or biomarkers.
[0088] Figure 4 This demonstrates a comparison between global and local Moran's I analysis workflows. Global Moran's I assesses the enrichment of the entire pathway, while local Moran's I further locates key nodes within the pathway.
[0089] Example 3 Figure 5 This is a structural diagram of a pathway enrichment analysis device based on Moran's I according to an embodiment of the present invention.
[0090] like Figure 5 As shown, the pathway enrichment analysis device based on Moran's I includes: Data acquisition module 301: This module is designed to acquire omics data from biological samples and determine molecular identifiers and their corresponding expression changes. It supports multiple data input formats, including CSV and Excel, and can automatically identify data columns and extract molecular identifiers and expression values.
[0091] Identifier mapping module 302: This module is configured to map the molecular identifiers to nodes in the pathway network of a preset pathway database. It has multiple built-in annotation database interfaces (AnnotationDbi, AnnotationHub), supports various identifier type conversions for proteins and metabolites, and can automatically handle one-to-many or many-to-one mapping relationships.
[0092] Weight matrix construction module 303: This module is configured to construct a spatial weight matrix based on the topology of the path network. It implements various weight calculation methods, including adjacency matrix, distance decay, and random walk, which users can select as needed.
[0093] Masking module 304: Configured to mask nodes in the pathway network that are not mapped to expressed data. This module records which nodes are missing data nodes, providing information for subsequent Moran's I mask calculation and permutation test.
[0094] Statistical calculation module 305: set to calculate the global spatial autocorrelation of the pathway using mask Moran's I statistic.
[0095] Significance testing module 306: Configured to evaluate the significance level of the global spatial autocorrelation using Monte Carlo permutation tests. This module implements parallel permutation computation, which can utilize multi-core CPUs to accelerate the permutation testing process.
[0096] Local Analysis Module 307: Set to calculate the local Moran's I statistic to identify key nodes in enriched pathways.
[0097] Results output module 308: Configured to determine a list of enriched pathways based on the significance level. This module supports multiple output formats, including data tables, visualizations, and HTML reports, facilitating user viewing and interpretation of the results.
[0098] Optionally, the device further includes: Visualization module 309: Set to generate pathway network diagrams, use different colors to mark the expression changes of nodes and local Moran's I values, and intuitively display the enrichment patterns of pathways.
[0099] Parameter configuration module 310: Set to allow users to customize analysis parameters, such as the number of permutations, significance threshold, weight matrix type, etc.
[0100] In the description of this specification, the references to terms such as "one embodiment," "some embodiments," "example," "specific example," or "some examples," etc., indicate that a specific feature, structure, material, or characteristic described in connection with that embodiment or example is included in at least one embodiment or example of the present invention. In this specification, the illustrative expressions of the above terms do not necessarily refer to the same embodiment or example. Furthermore, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in one or more embodiments or examples. Moreover, without contradiction, those skilled in the art can combine and integrate the different embodiments or examples described in this specification, as well as the features of different embodiments or examples.
[0101] Furthermore, the terms "first" and "second" are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the number of technical features indicated. Thus, a feature defined as "first" or "second" may explicitly or implicitly include at least one of that feature. In the description of this invention, "a plurality of" means at least two, such as two, three, etc., unless otherwise explicitly specified.
Claims
1. A pathway enrichment analysis method based on Moran's I, characterized in that, include: S1, acquire proteomics and / or metabolomics data of biological samples, and extract molecular identifiers and their corresponding expression values; S2, map the molecular identifier to the nodes of the pathway network in the preset pathway database, and construct a spatial weight matrix based on the topology of the pathway network. In the pathway network, the nodes represent biomolecules, the edges between the nodes represent the regulatory or interaction relationships between biomolecules, and the spatial weight matrix represents the proximity relationship between the nodes in the pathway network. S3, masking is performed on nodes that are not mapped to the expressed data, and the observed nodes and missing data nodes are marked by the mask vector; S4, extract sub-graphs based on mask vectors, and calculate the global spatial autocorrelation of pathways using the global Moran's I statistic; S5. The significance was assessed by the Monte Carlo permutation test, and the p-value corresponding to the global Moran's I statistic was calculated. S6. Repeat steps S2-S5 to obtain the global Moran's I statistic and its corresponding p-value for all pathways. Sort the pathways according to the p-value and select pathways with p-values less than the significance threshold as enriched pathways.
2. The method as described in claim 1, characterized in that, S2 further includes: S21, Obtain pathway network data from a preset pathway database and convert molecular identifiers of proteins and metabolites; When the molecular identifier is a protein identifier, use a gene annotation database to convert the identifier type and convert the input protein identifier into a standard identifier used by the pathway database; When the molecular identifier is a metabolite identifier, use a metabolite mapping table to convert the identifier type and the metabolite mapping table contains the correspondence between metabolite identifiers in different databases. S22, Match the converted identifier with the identifier of the path network node to establish the correspondence between the expressed data and the network node; S23, determine whether the transformed identifier matches the node successfully. For nodes that match successfully, assign the corresponding expression difference value to the node; for nodes that do not match successfully, mark them as missing data nodes. S24, obtain the adjacency matrix of the path network, the adjacency matrix indicates whether there is a direct connection relationship between each node, and calculate the spatial weight between nodes based on the network topology distance; The method for calculating the spatial weights between nodes based on network topology distance is as follows: Calculate the spatial weights between nodes using the shortest path distance weighting method, the inverse distance weighting method, or the Gaussian kernel weighting method; Before calculating the spatial weights, the spatial weight matrix is self-connected.
3. The method as described in claim 1, characterized in that, S4 further includes: S41, using the observation nodes in the mask vector, extract the expression values of the observation nodes from the complete expression value vector to form expression value sub-vectors, extract the weight relationships between the observation nodes from the complete spatial weight matrix to form observation node weight sub-matrix, and use the observation nodes, expression value sub-vectors and observation node weight sub-matrix to form a subgraph; S42, based on the extracted subgraphs, calculate the sum of subgraph weights and the sum of squares of subgraph representation values, and determine the subgraph degradation based on the sum of subgraph weights and the sum of squares of subgraph representation values. When the sum of subgraph weights is greater than 0, calculate Moran's I statistic; otherwise, output an invalid value. S43. Based on the results of Moran's I statistic, the global spatial autocorrelation of the pathway is determined. A Moran's I value close to 1 indicates strong positive spatial autocorrelation, where adjacent nodes in the pathway tend to be upregulated or downregulated simultaneously. A Moran's I value close to -1 indicates strong negative spatial autocorrelation, where adjacent nodes in the pathway tend to have opposite expression changes. A Moran's I value close to 0 indicates no spatial autocorrelation, where expression changes are randomly distributed in the pathway.
4. The method as described in claim 3, characterized in that, The formula for calculating Moran's I statistic is as follows: ; Where, n eff To observe the number of rows or columns in the node weight submatrix, This represents the sum of the subgraph weights. The subgraph represents the sum of squared values. To express the value sub-vector, This is the weight submatrix for the observation nodes.
5. The method as described in claim 3, characterized in that, The S5 also includes: S51, construct a permutation vector of length n, and randomly sample n values from all expression values. eff The extracted expression values are then used to fill the first n values of the permutation vector. eff There are n positions, and the remaining positions are filled with 0, where n eff To observe the number of rows or columns in the weight submatrix of the nodes; S52, randomly sort the elements in the permutation vector to obtain the reordered permutation vector, and construct the permutation mask based on the non-zero positions of the reordered permutation vector. S53, calculate the permutation mask statistic. Based on the positive and negative values of the permutation mask statistic and the relationship between the permutation mask statistic and the global Moran's I statistic, perform extreme value statistics. When the permutation mask statistic is ≥0 and the global Moran's I statistic is ≥ the permutation mask statistic, increment the first counter by one. When the permutation mask statistic is <0 and the global Moran's I statistic is ≤ the permutation mask statistic, increment the second counter by one. S54, calculate the one-sided p-value based on the sign of the permutation mask statistic. When the permutation mask statistic is ≥ 0, the formula for calculating the p-value is: p = (ge + 1) / (N + 1); Where ge is the count of the first counter, and N is the set total number of replacements; When the permutation mask statistic is less than 0, the formula for calculating the p-value is: p = (le + 1) / (N + 1); Where le is the number of times the second counter counts; When the one-sided p-value is less than the preset significance threshold, the spatial autocorrelation of the pathway is considered significant.
6. A method for analyzing key regulatory nodes in enrichment pathways based on local Moran's I, characterized in that, include: A1. Obtain proteomics and / or metabolomics data from biological samples, and extract molecular identifiers and their corresponding expression values. A2, map the molecular identifier to the nodes of the pathway network in the preset pathway database, and construct a spatial weight matrix based on the topology of the pathway network. In the pathway network, the nodes represent biomolecules, the edges between the nodes represent the regulatory or interaction relationships between biomolecules, and the spatial weight matrix represents the proximity relationship between the nodes in the pathway network. A3, performs masking on nodes that are not mapped to the expressed data, marks the observation nodes and missing data nodes with mask vectors, and uses the observation nodes as candidate key nodes; A4. Calculate the local Moran's I statistic for each candidate key node, and evaluate the significance of the local Moran's I for each candidate key node through permutation test. A5. Based on the expression value, local Moran's I statistic, and significance of candidate key nodes, determine whether they are key nodes and classify them. The key node classification includes high-high clustering pattern, low-low clustering pattern, high-low clustering pattern, and low-high clustering pattern.
7. The method as described in claim 6, characterized in that, The formula for calculating the local Moran's I statistic for each observation node is as follows: ; Where, n eff S0 is the number of observation nodes; S0 is the spatial weight matrix W. sub The sum of weights; z i Let z be the centered expression variance value of node i, which is the expression variance value of that node minus the mean of the expression values of the observed nodes; m T This represents the transpose of a vector.
8. The method as described in claim 6, characterized in that, The A4 also includes: A41, construct a permutation vector of length n, and randomly sample n values from all representation values. eff The extracted expression values are then used to fill the first n values of the permutation vector. eff There are n positions, and the remaining positions are filled with 0, where n eff To observe the number of rows or columns in the weight submatrix of the nodes; A42, randomly sort the elements in the permutation vector to obtain the reordered permutation vector, and construct the permutation mask based on the non-zero positions of the reordered permutation vector. A43. Calculate the local Moran's I statistic and local permutation mask statistic for each node. When the local Moran's I statistic of a node is a finite value, increment the local first counter by one. Perform extreme value statistics based on the positive and negative values of the local permutation mask statistic and the relationship between the local permutation mask statistic and the local Moran's I statistic. When the local permutation mask statistic is ≥0 and the local Moran's I statistic is ≥ the local permutation mask statistic, increment the local second counter by one. When the local permutation mask statistic is <0 and the local Moran's I statistic is ≤ the local permutation mask statistic, increment the local third counter by one. A44, calculate one-sided p based on the sign of the local permutation mask statistic. local When the local permutation mask statistic is ≥0, p local The formula for calculating the value is: p local [i] = (cnt ge [i] + 1) / (cnt included [i] + 1); Among them, cnt included The local first counter counts the number of times, cnt ge This is the number of times the local second counter counts. When the local permutation mask statistic is < 0, p local The formula for calculating the value is: p local [i] = (cnt le [i] + 1) / (cnt included [i] + 1); Among them, cnt le This is the number of times the local third counter counts. When p local When the value is less than the preset significance threshold, the node is considered significant.
9. The method as described in claim 6, characterized in that, The method for determining the classification of key nodes in step A5 is as follows: When the expression value of the key node is higher than the mean, the local Moran's I statistic is greater than 0, and the key node is significant, the key node is a high-high clustering pattern node. When the expression value of the key node is lower than the mean, the local Moran's I statistic is greater than 0, and the key node is significant, the key node is a low-low clustering pattern node. When the expression value of the key node is higher than the mean, the local Moran's I statistic is less than 0, and the key node is significant, the key node is a high-low clustering pattern node. When the expression value of a key node is lower than the mean, the local Moran's I statistic is less than 0, and the key node is significant, the key node is a low-high clustering pattern node.
10. A pathway enrichment analysis device based on Moran's I, characterized in that, include: Data acquisition module: used to acquire omics data from biological samples and determine each molecular identifier and its corresponding expression change value; Identifier mapping module: used to map the molecular identifier to nodes of the pathway network in a preset pathway database; Weight matrix construction module: used to construct a spatial weight matrix based on the topology of the pathway network; Masking module: used to mask nodes in the path network that are not mapped to the expressed data; The statistics calculation module is used to calculate the global spatial autocorrelation of the pathway using the masked Moran's I statistic. Significance testing module: used to assess the significance level of the global spatial autocorrelation using Monte Carlo permutation tests; Local Analysis Module: Used to calculate local Moran's I statistics and identify key nodes in enriched pathways; Results output module: used to determine the list of enriched pathways based on the significance level.