Method for identifying cancer driver modules based on graph embedding and hierarchical clustering

Abstract

This invention discloses a method for identifying cancer-driving modules based on graph embedding and hierarchical clustering, comprising the following steps: 1) setting a model; 2) constructing a weighted network; 3) extracting features of the weighted protein-protein interaction network (PPI) using graph embedding technology; 4) clustering the feature vectors; and 5) constructing a set of driving modules. This method identifies driving modules more effectively and accurately, and identifies driving modules with higher coverage and mutual exclusivity, while also detecting missed cancer-related genes.

Description

Technical Field

[0001] This invention relates to the fields of machine learning and bioinformatics, specifically a method for identifying cancer-driven modules based on graph embedding and hierarchical clustering. Background Technology

[0002] Cancer involves multiple omics, and current research indicates that cancer is caused by clonal proliferation resulting from mutations that give cells a selective growth advantage. These mutations, which promote cancer development, are called driver mutations. Driver mutations are the cause of tumorigenesis, capable of positive selection, driving the transformation of normal cells into proliferating cancer cells. The genes with these mutated mutations are called driver genes. Among a large number of somatic mutations, most occur during DNA replication and cell proliferation in both normal and cancer cells. Such mutations are called passenger mutations, and the genes with these mutations are called passenger genes. To mine somatic mutation data, several large-scale cancer genome sequencing projects have been launched, such as The Cancer Genome Atlas (TCGA) and the International Consortium on Cancer Genomes, providing a wide variety of high-throughput data for various cancer types. A significant challenge in cancer genome research is identifying driver mutations and driver genes. Functionally related driver genomes are also called driver modules or pathways. Recent research has shown that identifying and analyzing a single driver gene alone cannot effectively explain the complex genetic mechanisms of cancer. Therefore, driver modules composed of genes with driver mutations play a more important role in cancer development than individual driver genes. To address the shortcomings and low recognition rates of previous technologies, machine learning methods are being used to improve efficiency and accuracy.

[0003] Vandin et al. proposed the Dendrix method to identify driving pathways from somatic mutation data. This method introduces a weight W and maximizes a function W to achieve the objective, known as the maximum weighted submatrix problem. The Dendrix method can obtain gene sets with a large weight function W, but because it employs a stochastic search process, it may only find local optima.

[0004] The iMCMC method integrates somatic mutation, copy number variation, and gene expression data. However, iMCMC only calculates the mutual exclusivity of gene mutations within a module, and does not calculate the mutual exclusivity of genes to mutations.

[0005] The ComMDP and SpeMDP methods (Zhang et al., 2017) utilize integer linear programming to compute the maximum exclusion submatrix problem, focusing on identifying mutually exclusive gene sets. Due to the large number of combinations of mutated genes, de novo identification methods often reduce the inherent computational complexity by using pre-filtering based on mutation frequency, and may overlook some oncogenes or pathways containing rare mutations.

[0006] The HotNet2 method, proposed by Leiserson et al., is based on a network model and identifies driver modules and protein complexes through network analysis. HotNet2 uses an insulating thermal diffusion process that considers the mutation frequency of individual genes and the known topology of interactions between corresponding proteins to obtain a diffusion matrix from which cancer driver modules are identified. While many identified driver modules show significant co-occurrence in pan-cancer or single-cancer studies, mutual exclusion only occurs within modules and not between them.

[0007] The MEMCover method, as proposed by Kim et al., categorizes mutual exclusion into three types: intra-tissue mutual exclusion, cross-tissue type mutual exclusion, and inter-tissue type mutual exclusion. The subnetworks identified by this method not only include known pan-cancer dysregulation modules but also novel modules with unexplored cross-cancer effects. However, the selection of subnetworks in the MEMCover method is influenced by three factors: mutual exclusion, the connectivity of interacting networks, and sample coverage. Therefore, the subnetworks identified by this method must contain mutually exclusive gene sets.

[0008] The MEXCOwalk method, proposed by Ahmed et al., is an edge-weighted random walk approach. This method integrates mutation data with a protein-protein interaction network, combining gene mutual exclusion and coverage with protein interactions to form a vertex- and edge-weighted network. Vertex weights correspond to gene coverage, while edge weights correspond to gene mutual exclusion and the connectivity information of the protein-protein interaction network. MEXCOwalk can identify not only known significant oncogenes but also rare mutated oncogenes that are seldom mutated in pan-cancer data.

[0009] Wu et al. proposed ECSWalk, which introduces a vertex structure similarity based on JS divergence through network topology, improving the splitting and expansion methods in MEXCOwalk and enabling it to identify candidate gene sets with higher biological relevance and statistical significance. Summary of the Invention

[0010] The purpose of this invention is to address the shortcomings and low recognition rates of existing technologies by proposing a method for identifying cancer-driving modules based on graph embedding and hierarchical clustering. This method can identify driving modules more effectively and accurately, and can identify driving modules with higher coverage and mutual exclusivity, while also detecting missed cancer-related genes.

[0011] The technical solution to achieve the objective of this invention is:

[0012] The cancer-driven module identification method based on graph embedding and hierarchical clustering includes the following steps:

[0013] 1) Setting up the model:

[0014] Suppose there is a type L Mutation matrix of somatic cell data in cancer The row represents a group of cancer samples. , column represents a group of genes Each item in the matrix A value of 1 or 0 indicates a gene. In the sample Whether a value mutates is indicated by 1 for mutation and 0 for non-mutation. ,make Represents a sample set, in which genes These are all mutated genes, making The module set is as follows:

[0015] (1)

[0016] For any pair of genes ,in , It represents a selected set of genes, or a module, a subset of genes. Mutual exclusivity is usually expressed as:

[0017] (2)

[0018] like Then subset The genes in the sample are mutually exclusive, meaning that a subset of each sample... There is at least one gene mutation, and the gene subset Coverage is usually expressed as:

[0019] (3)

[0020] like Then the subset of genes Complete coverage of all patients, meaning that each sample has at least one subset. At least one gene within the gene has mutated;

[0021] Assumption It is a subset of genes. Although many previous studies have provided definitions of mutual exclusion and coverage for this, none of them have considered the impact of module size. Representation module Relative size, i.e. , The mutual exclusion score is defined as:

[0022] (4)

[0023] The coverage score is defined as:

[0024] (5);

[0025] For a graph and a group Genes that allow Indicates by The vertex corresponding to the gene in the middle The subgraph of cancer-driven module identification problem: given as input network Each gene of Total number of genes and minimum module size Find a set of disjoint modules. The module set score is usually defined as:

[0026] (6)

[0027] And meet the following conditions:

[0028] 1. They are connected.

[0029]

[0030] ;

[0031] 2) Construct a weighted network:

[0032] Given a network Representing the interconnected protein-protein interaction network, construct an undirected weighted graph. , where each vertex Indicates from genes The protein is expressed, and each undirected edge... Indicates gene and genes The corresponding interactions between proteins, in order to clarify the biological origin of the protein, which is derived from the gene. In network identification and citation, the name of the coding gene is used to express the encoding. To index and refer to the protein vertex, therefore, This represents both genes and corresponding vertices in a protein-protein interaction network; considering the close relationship between genes and their direct neighbor genes, the vertices are... The set of vertices of a vertex and its direct neighbors is defined as follows: As shown below:

[0033] (7)

[0034] To balance the mutual exclusion between genes and the chance of a gene coexisting with its surrounding genes, utilize and The average value is used as the mutual exclusivity of gene pairs in the network. As shown below:

[0035] (8)

[0036] To reduce the over- or under-coverage of a single gene, the coverage between gene pairs is represented by the product of the coverage of two genes. As shown below:

[0037] (9)

[0038] By introducing the confidence scores of genes and miRNAs from the mirDIP database, using Let be an expression, where 'e' represents an edge, and calculate the gene. and confidence levels between :

[0039] (10)

[0040] express and Interacting networks express The set of vertices, Contains only vertex, express The set of miRNA edges, , Contains only vertex, express and The set of edges between vertex sets. express and The set of edges between vertex sets. express and The number of edges between vertex sets. express and The number of edges between vertex sets, when and When there is no common miRNA vertex between them, that is ,but Set to 0.2;

[0041] A weighted undirected graph is calculated by combining three features: mutual exclusion, coverage, and confidence between gene pairs. The edge weights are as follows:

[0042] (11)

[0043] To further enhance the impact of mutual exclusion on weights, a threshold is introduced. ,for Rating less than The edge weight is assigned a value of 0, and finally for the edge weight in the edge weight is assigned a value of 0. Chinese gene pair The actual edge weights are as follows:

[0044] (12)

[0045] in Represents vertices The sum of the edge weights between it and its directly adjacent vertices;

[0046] 3) Extracting features of the weighted protein-protein interaction network (PPI) using graph embedding techniques:

[0047] Step 2) Constructed an undirected weighted graph ,exist The above applies vertex embedding, also known as graph embedding. Graph embedding is a broad technique that maps the entire graph to a low-dimensional space, converting graph data into feature vectors that can be understood by machine learning models. Vertex embedding is a more specific technique that maps vertex features to feature vectors that can be understood by machine learning models, thus better representing the similarity between vertices. Using the Node2vec framework proposed by Aditya Grover and Jure Leskovec for learning continuous feature representations of vertices in a network, the above applies... Feature learning in this context can be formulated as a maximum likelihood optimization problem, optimizing the following objective function:

[0048] (13)

[0049] in ,function Representing vertices respectively Mapping function to feature representation, Represents vertices generated through a neighborhood sampling strategy. The sampling strategy for the network neighborhood is a second-order biased random walk, meaning that reaching the next vertex depends not only on the current vertex but also on the previous vertex.

[0050] (14)

[0051] in, This represents the starting vertex of a second-order biased random walk. This indicates a second-order biased random walk that moves to the previous vertex of the current vertex, i.e. The previous vertex, Let Z represent the current vertex of a second-order biased random walk, and Z be the normalization constant. express and The non-normalized transition probabilities between them, i.e.:

[0052] (15)

[0053] As vertices With vertex The right to the border, This is represented as a search bias, i.e.:

[0054] (16)

[0055] Represents vertices To the top The shortest path distance, p and q, are the key parameters for deciding whether to use DFS or BFS. p controls the probability of immediately revisiting a vertex during the traversal. The random walk then tends to go closer to the vertex. The vertices of a sample are composed of vertices within a small local area. The walk tends to visit vertices that are far from the vertex. For distant vertices, use the Node2vec graph embedding method in After performing a second-order random walk using a neighborhood sampling strategy, a feature vector that can be understood by a machine learning model is obtained.

[0056] 4) Clustering of feature vectors:

[0057] In step 3), a feature vector was obtained, which can now be processed using machine learning methods. Hierarchical clustering is a cluster analysis method used to discover latent structures in a dataset and to extract useful information. In cluster analysis, the k-value (the number of groups in a cluster) refers to the number of clusters into which the dataset is divided. It affects the results and quality of clustering. If the k-value is too small, the intra-cluster variance will be too large, losing the meaning of clustering; if the k-value is too large, the intra-cluster variance will be too small, increasing computational complexity. Therefore, choosing a reasonable k-value is a key issue in cluster analysis. To solve this problem, the k-value is determined based on the number of sample gene sets. The minimum module size is used to determine the k-value for clustering. In studies of multiple cancer-driven modules, the minimum module size is mostly set to 3, such as in MEXCOwalk and ECSWalk. Hotnet2, however, sets the minimum module size to 2. However, for the size of the largest module However, this is rarely mentioned. Therefore, we propose an assumption: under reasonable and computationally convenient conditions, we set a maximum module size value. Of course, in reality, the maximum module size is obtained after clustering and processing, and the maximum module size is uncertain, so it will not directly affect the maximum module size in the actual result. A sequence of k values ​​can be obtained according to the following formula:

[0058]

[0059] According to settings Assuming the maximum module size is equal to the size of a single module, and ensuring that no two numbers are multiples of each other, we select the nearest prime number upwards to obtain the sequence. Finally, a suitable value of k is selected from the experiments;

[0060] In cluster analysis, the choice of distance to measure the similarity or dissimilarity between samples affects the clustering results and effectiveness. The impact of distance on clustering is a broad topic, and different clustering algorithms and distance metrics may have different effects. Generally speaking, the smaller the distance, the higher the similarity between data points, and the easier it is for them to be assigned to the same cluster; conversely, the larger the distance, the greater the dissimilarity between data points, and the easier it is for them to be assigned to different clusters. Since step 3) yields a multidimensional feature vector, Mahalanobis distance is chosen. It corrects the Euclidean distance. To address the issue of inconsistent and correlated scales across different dimensions, and to better represent the similarity between data, we assume... , The eigenvectors in the feature matrix are as follows: ,but and The Mahalanobis distance between them is:

[0061]

[0062] in It is the covariance matrix of multidimensional random variables. Mahalanobis distance eliminates the interference of correlation between variables, making the distance reasonable. Using the hierarchical clustering analysis above, a cluster is obtained.

[0063] 5) Construct the driver module set:

[0064] For the clusters obtained in step 4), the following cluster processing algorithm is proposed, which requires three main steps: namely, sorting and scoring the gene set, deleting genes, and extracting sub-connected components from each cluster to obtain the final cancer-driving module; firstly, a set of candidate modules is created, based on... The network topology selects vertices from each cluster with a number greater than or equal to 1. The sub-connected components are added as candidate modules to the initial module set P, and then based on... The vertex strength and sample coverage are used to score and rank each gene in the initial module set P, resulting in a gene sorting sequence. In the third step, the last few genes in the gene sorting sequence are removed from P, and the process returns to the previous two steps to update the module set P and the gene sorting sequence, until the total number of genes in the module set P is reduced to a certain value. The purpose of this step is to dynamically delete genes with low scores and then, based on... The connection structure in the image is re-acquired to obtain sub-connected components. After the above process, the final cancer driver module gene set is obtained. The cluster processing algorithm for constructing the driver module set is as follows:

[0065] Algorithm: Cluster Processing

[0066] Input: cluster Weighted undirected graph ;parameter Gene size parameters ,

[0067] Step length Minimum module size ;

[0068] Output: A set of driver modules

[0069] Foreach do

[0070] The connected subgraphs in are stored middle;

[0071] Foreach do

[0072] like but

[0073] Will Store to module set middle;

[0074] Finish

[0075] Finish

[0076] Finish

[0077] use function pairs Gene scoring yields a gene sequencing sequence. ;

[0078] when do

[0079] From the set Delete The last one One vertex;

[0080] ;

[0081] Set to empty;

[0082] foreach do

[0083] The connected subgraphs in the data are stored in middle;

[0084] foreach do

[0085] like but

[0086] Will Store to module set middle;

[0087] Finish

[0088] Finish

[0089] Finish

[0090] like and but

[0091] Output

[0092]

[0093] Finish

[0094] Finish

[0095] function:

[0096] Input: Weighted undirected graph ;vertex

[0097] Return: Vertex Score

[0098] foreach do

[0099]

[0100] end

[0101] ,in for coverage

[0102] return .

[0103] This technical solution introduces a miRNA-gene interaction network to reduce the negative impact of noise. miRNAs are important regulators of gene expression, which are regulated by binding to the genes to be regulated. This technical solution considers factors such as high mutual exclusion and high coverage between genes, as well as topological structure, and uses graph embedding technology to extract more hidden information between vertices. This not only enables more effective and accurate identification of driving modules, but also allows for the identification of driving modules with higher coverage and mutual exclusion. Attached Figure Description

[0104] Figure 1 Flowchart for an embodiment;

[0105] Figure 2 This is a schematic diagram of parameter settings in the embodiment;

[0106] Figure 3 The results are shown in the example.

[0107] Figure 4 This is a diagram illustrating the effect when the gene size is 100, as shown in the example. Detailed Implementation

[0108] The present invention will now be described in further detail with reference to the accompanying drawings and specific embodiments, but this is not intended to limit the scope of the invention.

[0109] Example:

[0110] 1. Data preparation for this example:

[0111] First, somatic mutation data, miRNA-gene interaction network data, and protein interaction network data of the samples were downloaded from the TCGA and mirDIP databases, respectively.

[0112] 2. Parameter settings in this example: Figure 2 As shown:

[0113] Set the mutual exclusion threshold to 0.7; the confidence threshold to 0.2;

[0114] Dimension 32; p-value 4; q-value 1; number of groups k is 127;

[0115] Each time, the number of vertices deleted is del_num(dn) = 1; the minimum module size is 3; the gene size is 2500-100, and the step size is 100.

[0116] Reference Figure 1 The cancer-driven module identification method based on graph embedding and hierarchical clustering includes the following steps:

[0117] 1) Setting up the model:

[0118] Suppose there is a type L Mutation matrix of somatic cell data in cancer The row represents a group of cancer samples. , column represents a group of genes Each item in the matrix A value of 1 or 0 indicates a gene. In the sample Whether a value mutates is indicated by 1 for mutation and 0 for non-mutation. ,make Represents a sample set, in which genes These are all mutated genes, making The module set is as follows:

[0119] (1)

[0120] For any pair of genes ,in , Represents a selected set of genes, signifying a module, or a subset of genes. Mutual exclusivity is usually expressed as:

[0121] (2)

[0122] like Then subset The genes in the sample are mutually exclusive, meaning that a subset of each sample... There is at least one gene mutation, and the gene subset Coverage is usually expressed as:

[0123] (3)

[0124] like Then the subset of genes Complete coverage of all patients, meaning that each sample has at least one subset. At least one gene within the gene has mutated;

[0125] Assumption It is a subset of genes. Although many previous studies have provided definitions of mutual exclusion and coverage for this, none of them have considered the impact of module size. Representation module Relative size, i.e. , The mutual exclusion score is defined as:

[0126] (4)

[0127] The coverage score is defined as:

[0128] (5);

[0129] For a graph and a group Genes that allow Indicates by The vertex corresponding to the gene in the middle The subgraph of cancer-driven module identification problem: given as input network Each gene of Total number of genes and minimum module size Find a set of disjoint modules. The module set score is usually defined as:

[0130] (6)

[0131] And meet the following conditions:

[0132] 1. They are connected.

[0133]

[0134] ;

[0135] 2) Construct a weighted network;

[0136] Given a network Representing the interconnected protein-protein interaction network, construct an undirected weighted graph. , where each vertex Indicates from genes The protein is expressed, and each undirected edge... Indicates gene and genes The corresponding interactions between proteins, in order to clarify the biological origin of the protein, which is derived from the gene. In network identification and citation, the name of the coding gene is used to express the encoding. To index and refer to the protein vertex, therefore, This represents both genes and corresponding vertices in a protein-protein interaction network; considering the close relationship between genes and their direct neighbor genes, the vertices are... The set of vertices of a vertex and its direct neighbors is defined as follows: As shown below:

[0137] (7)

[0138] To balance the mutual exclusion between genes and the chance of a gene coexisting with its surrounding genes, utilize and The average value is used as the mutual exclusivity of gene pairs in the network. As shown below:

[0139] (8)

[0140] To reduce the over- or under-coverage of a single gene, the coverage between gene pairs is represented by the product of the coverage of two genes. As shown below:

[0141] (9)

[0142] This example introduces the confidence scores of genes and miRNAs from the mirDIP database, using... Let be an expression, where 'e' represents an edge, and calculate the gene. and confidence levels between :

[0143] (10)

[0144] express and Interacting networks express The set of vertices, Contains only vertex, express The set of miRNA edges, , Contains only vertex, express and The set of edges between vertex sets. express and The set of edges between vertex sets. express and The number of edges between vertex sets. express and The number of edges between vertex sets, when and When there is no common miRNA vertex between them, that is ,but Set to 0.2;

[0145] This example combines three features—mutual exclusion, coverage, and confidence—between gene pairs to calculate a weighted undirected graph. The edge weights are as follows:

[0146] (11)

[0147] To further enhance the impact of mutual exclusion on weights, a threshold is introduced. ,for Rating less than The edge weight is assigned a value of 0, and finally for the edge weight in the edge weight is assigned a value of 0. Chinese gene pair The actual edge weights are as follows:

[0148] (12)

[0149] in Represents vertices The sum of the edge weights between it and its directly adjacent vertices;

[0150] 3) Extracting features of the weighted protein-protein interaction network (PPI) using graph embedding techniques:

[0151] Step 2 constructs an undirected weighted graph. ,exist The above applies vertex embedding, also known as graph embedding. Graph embedding is a broad technique that maps the entire graph to a low-dimensional space, converting graph data into feature vectors that can be understood by machine learning models. Vertex embedding is a more specific technique that maps vertex features to feature vectors that can be understood by machine learning models, thus better representing the similarity between vertices. Using the Node2vec framework proposed by Aditya Grover and Jure Leskovec for learning continuous feature representations of vertices in a network, the above applies... Feature learning in this context can be formulated as a maximum likelihood optimization problem, optimizing the following objective function:

[0152] (13)

[0153] in ,function Representing vertices respectively Mapping function to feature representation, Represents vertices generated through a neighborhood sampling strategy. The sampling strategy for the network neighborhood is a second-order biased random walk, meaning that reaching the next vertex depends not only on the current vertex but also on the previous vertex.

[0154] (14)

[0155] in, This represents the starting vertex of a second-order biased random walk. This indicates a second-order biased random walk that moves to the previous vertex of the current vertex, i.e. The previous vertex, Let Z represent the current vertex of a second-order biased random walk, and Z be the normalization constant. express and The non-normalized transition probabilities between them, i.e.:

[0156] (15)

[0157] As vertices With vertex The right to the border, This is represented as a search bias, i.e.:

[0158] (16)

[0159] Represents vertices To the top The shortest path distance, p and q, are the key parameters for deciding whether to use DFS or BFS. p controls the probability of immediately revisiting a vertex during the traversal. The random walk then tends to go closer to the vertex. The vertices of a sample are composed of vertices within a small local area. The walk tends to visit vertices that are far from the vertex. For distant vertices, use the Node2vec graph embedding method in After performing a second-order random walk using a neighborhood sampling strategy, a feature vector that can be understood by a machine learning model is obtained.

[0160] 4) Clustering of feature vectors:

[0161] In step 3), a feature vector was obtained, which can now be processed using machine learning methods. Hierarchical clustering is a cluster analysis method used to discover latent structures in a dataset and to extract useful information. In cluster analysis, the k-value (the number of groups in a cluster) refers to the number of clusters into which the dataset is divided. It affects the results and quality of clustering. If the k-value is too small, the intra-cluster variance will be too large, losing the meaning of clustering; if the k-value is too large, the intra-cluster variance will be too small, increasing computational complexity. Therefore, choosing a reasonable k-value is a key issue in cluster analysis. To solve this problem, this example uses the number of sample gene sets... The minimum module size is used to determine the k-value for clustering. In studies of multiple cancer-driven modules, the minimum module size is mostly set to 3, such as in MEXCOwalk and ECSWalk. Hotnet2, however, sets the minimum module size to 2. However, for the size of the largest module However, this is rarely mentioned. Therefore, this example proposes an assumption: under reasonable and computationally convenient conditions, a maximum module size is set. Of course, in reality, the maximum module size is obtained after clustering and processing, and the maximum module size is uncertain, so it will not directly affect the maximum module size in the actual result. A sequence of k values ​​can be obtained according to the following formula:

[0162]

[0163] According to settings Assuming the maximum module size is equal to the size of a single module, and ensuring that no two numbers are multiples of each other, we select the nearest prime number upwards to obtain the sequence. Finally, a suitable value of k is selected from the experimental results;

[0164] In cluster analysis, the choice of distance to measure the similarity or dissimilarity between samples affects the clustering results and effectiveness. The impact of distance on clustering is a broad topic, and different clustering algorithms and distance metrics may have different effects. Generally speaking, the smaller the distance, the higher the similarity between data points, and the easier it is for them to be assigned to the same cluster; conversely, the larger the distance, the greater the dissimilarity between data points, and the easier it is for them to be assigned to different clusters. Since step 3) yields a multidimensional feature vector, Mahalanobis distance is chosen. It corrects the Euclidean distance. To address the issue of inconsistent and correlated scales across different dimensions, and to better represent the similarity between data, we assume... , The eigenvectors in the feature matrix are as follows: ,but and The Mahalanobis distance between them is:

[0165]

[0166] in It is the covariance matrix of multidimensional random variables. Mahalanobis distance eliminates the interference of correlation between variables, making the distance reasonable.

[0167] After using the hierarchical clustering analysis described above, a cluster set is obtained;

[0168] 5) Construct the driver module set:

[0169] For the cluster set obtained in step 4), this example proposes the following cluster processing algorithm, which requires three main steps: namely, sorting and scoring the gene set, deleting genes, and extracting sub-connected components from each cluster to obtain the final cancer-driving module; firstly, a set of candidate modules is created, based on... The network topology selects vertices from each cluster with a number greater than or equal to 1. The sub-connected components are added as candidate modules to the initial module set P, and then based on... The vertex strength and sample coverage are used to score and rank each gene in the initial module set P, resulting in a gene sorting sequence. In the third step, the last few genes in the gene sorting sequence are removed from P, and the process returns to the previous two steps to update the module set P and the gene sorting sequence, until the total number of genes in the module set P is reduced to a certain value. The purpose of this step is to dynamically delete genes with low scores and then, based on... The connection structure in the image is re-acquired to obtain sub-connected components. After the above process, the final cancer driver module gene set is obtained. The cluster processing algorithm for constructing the driver module set is as follows:

[0170] Algorithm: Cluster Processing

[0171] Input: cluster Weighted undirected graph ;parameter Gene size parameters ,

[0172] Step length Minimum module size ;

[0173] Output: A set of driver modules

[0174] Foreach do

[0175] The connected subgraphs in are stored middle;

[0176] Foreach do

[0177] like but

[0178] Will Store to module set middle;

[0179] Finish

[0180] Finish

[0181] Finish

[0182] use function pairs Gene scoring yields a gene sequencing sequence. ;

[0183] when do

[0184] From the set Delete The last one One vertex;

[0185] ;

[0186] Set to empty;

[0187] foreach do

[0188] The connected subgraphs in the data are stored in middle;

[0189] foreach do

[0190] like but

[0191] Will Store to module set middle;

[0192] Finish

[0193] Finish

[0194] Finish

[0195] like and but

[0196] Output

[0197]

[0198] Finish

[0199] Finish

[0200] function:

[0201] Input: Weighted undirected graph ;vertex

[0202] Return: Vertex Score

[0203] foreach do

[0204]

[0205] end

[0206] ,in for coverage

[0207] return .

[0208] like Figure 3 The results are shown below.

[0209] When input gene size ( The range is 2500 to 100. When the step size is 100, an "out" folder will be output. Inside the "out" folder, a "txt" folder contains 25 "txt" files, which are the final result files. Each "txt" file records the 2500, 2400, 2300, ..., 100 genes respectively. The other files in the "out" folder are process files, storing information such as gene coverage, mutual exclusion, confidence, edge weight information, feature matrices, and clustering results generated during the process. During the run, the number of modules and gene size information will be continuously printed. The terminal displays a running time of 7068 seconds.

[0210] like Figure 4 The following is a demonstration of the effect when the gene size is 100:

[0211] When the total number of genes In this example, the method outputs a module. The vertices of the square correspond to genes. The size of a vertex is proportional to the mutation frequency of the corresponding gene. The edge weights reflect the thickness of the line segment. A module corresponds to the cancer type with the highest gene mutation enrichment in that module, and each module is named after the gene with the highest enrichment in that module.

Claims

1. A method for identifying cancer driver modules based on graph embedding and hierarchical clustering, characterized in that, Comprising the following steps: 1) Set the model: Suppose there is a mutation matrix from somatic data of L kinds of cancers where rows represent a set of cancer samples , columns represent a set of genes , and each item in the matrix has a value of 1 or 0, indicating whether the gene is mutated in the sample , 1 indicating mutation and 0 indicating no mutation, for each , let represent a sample set in which the genes are all mutated genes, and let be a module set, as follows:​ (1) For any pair of genes where , represents a selected set of genes, also a module, a subset of genes The mutual exclusivity of a set of genes is represented as: (2) If , then the genes in the subset are mutually exclusive, i.e., at least one gene in the subset is mutated in each sample, the coverage of the gene subset is represented as: （3） If then the gene subset completely covers all patients, i.e. each sample has at least one subset in which at least one gene is mutated; Assume is a subset of genes, with denotes a module relative size, i.e. , The mutual exclusivity score of two genes is defined as: （4） The coverage score of a coverage is defined as: （5）； For a graph and a set of genes let denote the subgraph of vertices corresponding to the genes The cancer driver module identification problem: Given as input a network , the number of genes , the number of genes , the total number of genes and the minimum module size , find a set of disjoint modules The module set score is typically defined as: （6） And meet the following conditions: (1) are connected; （2） ； （3） ； 2) Build a weighted network: Given a network representing a protein interaction network of connected proteins, construct an undirected weighted graph where each vertex represents a protein produced from expression of a gene each undirected edge represents an interaction between the proteins corresponding to genes and to identify the biological source of a protein encoded for expression by a gene in the network, use the name of the encoding gene to index and refer to the protein vertex in the network's identification and references, thus, both representing the gene and the protein with the corresponding vertex in the protein interaction network; considering that a gene has a close connection with its immediate neighbors, define the set of vertices and their immediate neighbor vertices as as follows: （7） To balance the exclusivity between genes and the chance of coexistence of a gene with its surrounding genes, the average of and is used as the exclusivity of a gene pair in the network as follows: （8） To reduce the coverage of individual genes being either too large or too small, the coverage product of two genes is used to represent the coverage between a pair of genes As follows: （9） By introducing the confidence of genes and miRNAs on the mirDIP database, the confidence between genes and miRNAs is calculated as follows: wherein e represents an edge, and the confidence between genes and is calculated as follows: : (10) denotes interaction network, denotes a set of vertices, contains only vertices, denotes a set of miRNA edges, , contains only vertices, denotes a set of edges between vertices, denotes a set of edges between vertices, denotes the number of edges between vertices, denotes the number of edges between vertices, when there is no common miRNA vertex between and , is set to 0.2; The three features of exclusivity, coverage and confidence between pairs of genes are integrated to calculate the edge weights of the weighted undirected graph as follows: (11) To further strengthen the influence of mutual exclusivity on weights, a threshold is introduced For edges with score less than 0 is assigned to the edge weight, and finally for the actual edge weight between genes in pairs is as follows: (12) wherein denotes a vertex the sum of the edge weights between its directly adjacent vertices; 3) Extract the features of the weighted protein-protein interaction network PPI by graph embedding technology: Step 2 constructs an undirected weighted graph. ,exist The above applies vertex embedding, also known as graph embedding. Graph embedding is a broad technique that maps the entire graph to a low-dimensional space, converting graph data into feature vectors that can be understood by machine learning models. Vertex embedding is a more specific technique that maps vertex features to feature vectors that can be understood by machine learning models, thus better representing the similarity between vertices. Using the Node2vec framework proposed by Aditya Grover and Jure Leskovec for learning continuous feature representations of vertices in a network, the above applies... Feature learning in this context can be formulated as a maximum likelihood optimization problem, optimizing the following objective function: (13) wherein , the function denotes a mapping function from the vertices to the feature representation, denotes the network neighborhood of a vertex generated by a neighborhood sampling strategy, which is a second order biased random walk, i.e. the next vertex reached depends not only on the current vertex but also on the previous vertex, i.e.: (14) wherein, denotes the starting vertex of the second order biased random walk, denotes the previous vertex of the current vertex of the second order biased random walk, i.e. the previous vertex of Z, denotes the current vertex of the second order biased random walk, Z is a normalization constant, denotes the non-normalized transition probability between and (15) for vertex edge weight, is represented as a search bias, i.e.: (16) representing a vertex to vertex , p and q are the key parameters that determine whether to use DFS or BFS, p controls the possibility of immediately re-visiting a vertex during the walk, , then the random walk is biased towards vertices that are close to vertex , i.e. the sample consists of vertices within a small local neighborhood; , then the walk is biased towards visiting vertices that are far from vertex , Node2vec graph embedding method performs a second-order random walk on by a neighborhood sampling strategy and obtains a feature vector that can be understood by machine learning models; 4) Cluster the feature vectors: Hierarchical clustering analysis was used to perform cluster analysis on the feature vectors obtained in step 3). In cluster analysis, the k-value refers to the number of clusters into which the dataset is divided. It affects the results and quality of clustering. If the k-value is too small, the intra-cluster variance will be too large, losing the meaning of clustering; if the k-value is too large, the intra-cluster variance will be too small, increasing computational complexity. Therefore, choosing a reasonable k-value is a key issue in cluster analysis. The k-value is determined based on the number of gene sets in the sample. The value of k is determined by the module size. In studies involving multiple cancer-driving modules, the k-value is set... Under reasonable and computationally convenient conditions, a maximum module size sequence is set; a sequence of k values ​​can be obtained according to the following formula: According to the setting , assuming the maximum module size is the size of a single module, and also making the two numbers not have a multiple relationship, the nearest prime number is selected upwards to obtain the sequence , and finally a suitable k value is selected according to the result; In cluster analysis, the choice of distance to measure the similarity or difference between samples affects the clustering results and effectiveness. Since step 3) yields a multidimensional feature vector, Mahalanobis distance is chosen, assuming... , The eigenvectors in the feature matrix are as follows: ,but and The Mahalanobis distance between them is: wherein is the covariance matrix of the multidimensional random variables, the Mahalanobis distance excludes the interference of correlation between variables, so that the distance is reasonable; After using the above hierarchical clustering analysis, a cluster set is obtained; 5) Build a driving module set: For the cluster set obtained in step 4), the following cluster processing algorithm is proposed, which needs to adopt three main steps: that is, scoring and sorting the gene set, deleting genes, and taking sub-connected components of each cluster to obtain the final cancer driving module; first, a set of candidate modules is created, and according to the network topology structure of , the sub-connected components with the number of vertices greater than or equal to are selected from each cluster, and the sub-connected components are added to the initial module set P as candidate modules, then each gene in the initial module set P is scored and sorted according to the vertex intensity and sample coverage in , and a gene sorting sequence is obtained, and in the third step, the last several genes in the gene sorting sequence are removed from P and the module set P and the gene sorting sequence are updated back to the above two steps until the total number of genes in the module set P is reduced to , the purpose of this step is to dynamically delete genes with low scores and then re-obtain sub-connected components according to the connection structure in , and after the above process, the final cancer driving module gene set is obtained, and the cluster processing algorithm process for constructing the driving module set is as follows: Algorithm: cluster set processing Input: cluster set ; weighted undirected graph ; parameters ; gene size parameters , , step size ; minimum module size ; Output: a set of driver modules Foreach do the connected sub-graphs in the are stored in the connected sub-graphs in the Foreach do If then To store into a module set ; End End End With the function on the gene scores in a gene ranking sequence When do from the set delete last vertex; ； Set to null; foreach do the connected subgraph in the graph G is stored in the graph G; foreach do If then To store into a module set ; End End End If and then Output End End Function: Input: Weighted undirected graph ; vertices Return: vertex score foreach do End wherein is coverage Return .

Images