A pan-cancer common driver pathway identification method based on GAN sample balance

Abstract

The present invention discloses a pan-cancer common driver pathway identification method based on GAN sample balance, comprising the following steps: 1) generating somatic mutation data of the corresponding cancer that conforms to the real data distribution; 2) minimizing the model CDP-HA of the discreteness between the total weights of each cancer; 3) introducing a parthenogenetic algorithm to solve the model CDP-HA. This method becomes a useful tool for identifying pan-cancer common driver pathways, with strong scalability and practicality.

Description

Technical Field

[0001] The present invention relates to the field of bioinformatics and is used for identifying cancer driver pathways, in particular to a pan-cancer public driver pathway identification method based on GAN sample balance. Background Art

[0002] Cancer is a very complex disease that threatens human health, and its etiology includes multiple genetic and environmental factors. Understanding the carcinogenic mechanism at the molecular level is a huge challenge, which is beneficial to the diagnosis, treatment and drug design of cancer in medicine. With the rapid development of next-generation sequencing technology (NGS), researchers can better describe the molecular characteristics of cancer. At present, several large-scale cancer genome projects (The Cancer Genome Atlas (TCGA), the International Cancer Genome Consortium (ICGC), and the Cancer Cell Line Encyclopedia (CCLE)) have generated and analyzed massive amounts of data, providing unprecedented opportunities for further understanding of the molecular mechanisms and carcinogenic mechanisms of cancer. Existing studies have shown that only functional driver mutations promote the development of cancer, while passenger mutations have little effect on cancer development. Distinguishing between functional driver mutations and passenger mutations has become an important task in studying the pathogenesis of cancer.

[0003] Most early studies were devoted to designing methods that can effectively identify individual driver genes with significantly higher mutation rates. However, cancer has huge mutational heterogeneity, that is, different driver gene mutations in the same cancer can lead to cancer. Therefore, identifying a single driver gene cannot effectively understand the mechanism of cancer development. Further studies have shown that the occurrence of cancer is usually caused by the destruction of some pathways, and different combinations of driver mutations may disrupt these pathways (cell signaling or regulatory pathways). Therefore, identifying driver pathways has become a key way to understand the carcinogenic mechanism of cancer at the pathway level. At present, the problem of driver pathway identification can be divided into three directions: identifying a single driver pathway, identifying collaborative driver pathways, and identifying pan-cancer common driver pathways and specific driver pathways. This article mainly studies the problem of identifying pan-cancer common driver pathways.

[0004] The purpose of identifying common driver pathways on a pan-cancer scale is to study the commonalities that may exist between different cancer types, which is conducive to enhancing the understanding of cancer pathogenic mechanisms. The TCGA Pan-Cancer Project has collected multi-platform mutation data generated by thousands of cancer patients with 12 cancer types, providing an opportunity to further study such issues. Recently, a class of research methods based on prior knowledge has been proposed. These methods usually use gene-gene interaction (GGI) networks, protein-protein interaction (PPI) networks, and pathway-pathway interaction (PaPaI) networks. Although they can achieve good identification results, relying on prior knowledge will miss the discovery of better mutation gene combinations on the one hand, and limit the scope of finding pathways on the other hand, because the existing prior knowledge is not perfect and contains partial pathway information. For example, Leiserson et al. proposed the HotNet2 method based on the directed heat diffusion model, which attempts to obtain pathways and protein complexes by combining protein-protein interaction networks for analysis; Kim et al. studied different types of mutual exclusivity between multiple cancer types and proposed the MEMCover method for identifying subnetworks / pathways based on the HumanNet network. Hajkarim et al. proposed the DAMOKLE algorithm based on a large gene-gene interaction network, which attempts to identify subnetworks with significantly different mutation frequencies in samples of two cancers. The other type is the de novo identification method. Zhang et al. proposed the ComMDP method, which utilizes two characteristics of the driving pathway: high mutual exclusivity and high coverage, and then directly extends the maximum weight submatrix problem model applicable to a single cancer to multiple cancer types, attempting to identify common driving pathways by accumulating absolute weight values. Wu et al. introduced the CDP-V model, which uses relative proportions instead of absolute numbers and uses variance to minimize the discreteness of each proportion. At the same time, they proposed the CDP-H model that uses the harmonic mean to minimize the discreteness of each proportion, reducing the use of parameters. Try to identify common driving pathways through relative weight values.

[0005] In the above research methods, the ComMDP method uses absolute weight values. It accumulates the absolute weight values ​​of each cancer type to obtain a set of genes with the largest weight values. It does not take into account the problem of unbalanced sample sizes between cancer types. Therefore, when the sample size difference of each cancer type is large, the identification results will be more inclined to the cancer type with a larger sample size, thereby missing some driving pathways with higher commonality. In addition, the PGA-V method uses relative weight calculation to solve the problem of unbalanced sample size, but it introduces an artificial parameter. The determination of the parameter requires a lot of experiments, and the experimental process is cumbersome, which is not conducive to expansion in practical applications. However, the problem of unbalanced sample size has not been actually solved.

[0006] With the rapid development of deep learning, solving the problem of data imbalance has become more urgent. In 2012, the data augmentation strategy was proposed, which is to generate additional data items that conform to the real data distribution by converting existing real data. But there was no better way to generate data in the same period. It was not until 2014 that Goodfellow et al. proposed a powerful generative model based on game theory - Generative Adversarial Networks (GANs). Although a large number of generative models based on deep learning were born in the early days, GANs is one of the most successful generative models and has been successfully applied to enhance data in various fields. Summary of the invention

[0007] The purpose of the present invention is to address the deficiencies of the prior art and provide a pan-cancer common driver pathway identification method based on GAN sample balance. This method becomes a useful tool for identifying pan-cancer common driver pathways with strong scalability and practicality.

[0008] The technical solution for achieving the purpose of the present invention is:

[0009] A pan-cancer common driver pathway identification method based on GAN sample balance includes the following steps:

[0010] 1) Generate somatic mutation data corresponding to cancer that conforms to the real data distribution:

[0011] 1.1) Set up the adversarial generation network framework:

[0012] Assume a system with m r samples, n r The example training set of genes, the generator network of SNV-GANs is defined as G(z), the input of the generator is: z~norm(0,1), and the definition of the generator is as follows:

[0013] 1.1.1) The input layer uses GFC1 to map the noise vector z into a tensor zn of dimension (1,128);

[0014] 1.1.2) The hidden layer puts the tensor zn in step 1.1.1) into GFC2 for mapping, and the result is put into GFC3 for mapping. The same is true for GFC4, and finally it is mapped into a tensor zn′ of dimension (1,1024);

[0015] 1.1.3) The output layer maps the tensor zn′ to dimension (1,m) through GFC5 r *n r ) tensor gn, and then reset gn to a tensor of dimension (m r ,n r ) Tensor is the output of the generator;

[0016] The input layer and the hidden layer both use the dropout function to freeze some neurons, and use the activation function ReLU defined by formula (1). The output layer uses the activation function Sigmod defined by formula (2). The discriminator network is defined as D(x), and the input of the discriminator is the real data x~P real Or generate data x represents a set of somatic mutation data samples, and the discriminator is defined as follows:

[0017] 1.1.4) The input layer uses DFC1 to map x to dimension (m r ,256) tensor xn;

[0018] 1.1.5) The hidden layer puts the tensor xn from step 1.1.4) into DFC2 for mapping, and the result is put into DFC3 for mapping. The same is true for DFC4, and the final mapping is dimension (m r ,16) tensor xn′;

[0019] 1.1.6) The output layer maps the tensor xn′ to dimension (m) through DFC5 r ,1) tensor dn, which is the output of the discriminator; the input layer and the hidden layer both use the activation function ReLU defined by formula (1), and the output layer uses the activation function Sigmod defined by formula (2);

[0020] ReLU(x)=max(0,x) (1),

[0021]

[0022] 1.2) Training process of SNV-GANs:

[0023] 1.2.1) Given a somatic mutation matrix A r (m r ×n r ) and a random sample ratio parameter ∝ (∝ < 1), according to the ratio parameter ∝ from the matrix A r The number of randomly selected samples is m, m = m r *∝'s submatrix M r , size is m×n r , a total of 64 sub-matrices M are extracted r Construct a training set X and input it into the generative adversarial network for training;

[0024] 1.2.2) Initialize the parameters θ of the discriminator D(.) d , the parameters θ of the generator G(.) g ;

[0025] 1.2.3) Set the current round epoch = 1, and randomly generate a 1×100 Gaussian distributed noise vector z;

[0026] 1.2.4) Take z from step 1.2.3) as the input of the generator to get a size of m×n r Vector 1.2.5) Calculate the generator loss value according to formula (3), and then update the generator parameter θ g :

[0027]

[0028] Where G(z (i) ) indicates that the generator passes the noise vector z (i) The generated data, D(G(z (i) )) represents the probability that the discriminator generates data that is judged to be real data, loss G The smaller the better;

[0029] 1.2.6) Randomly select a sample group x from the training set X;

[0030] 1.2.7) Randomly generate a 1×100 Gaussian distributed noise vector z;

[0031] 1.2.8) Take z from step 1.2.7) as the input of the generator to get a size of m×n r Vector 1.2.9) Calculate the generator loss value according to formula (4), and then update the generator parameter θ g :

[0032]

[0033] Where D(x (i) ) represents the probability that the discriminator determines the generated data as real data, 1-D(G(z (i) )) represents the probability that D determines the generated data as generated data, loss D Bigger is better;

[0034] 1.2.10) Determine whether the current epoch has reached the set maximum epoch: if so, stop training; otherwise, return to step 1.2.3) to finally obtain the trained generator G(.);

[0035] 1.3) Data processing:

[0036] 1.3.1) Randomly generate a 1×100 Gaussian distributed noise vector z;

[0037] 1.3.2) Input the vector z in step 1.3.1) into the generator G(.) obtained through training to get the generated data G data = G(z);

[0038] 1.3.3) Set the values greater than or equal to 0.85 in G data to 1 and the values less than 0.85 to 0 to obtain a new binary matrix A fakedata ;

[0039] 1.3.4) Take the sample size of the cancer with the largest number of samples in the R cancer somatic mutation matrices as m max , 0 < max < R; then insert the matrix A r into the amplification matrix At this time, the sample size of the amplification matrix is

[0040] 1.3.5) If the number of samples to be amplified is greater than the sample size m of the matrix A r , execute step 1.3.6), if the number of samples to be amplified r is less than the sample size m of the matrix A , execute step 1.3.8); r of the matrix A r ;

[0041] Step 1.3.6) Randomly select m samples from A in step 1.3.3 fakedata respectively calculate the mutation rate of each gene in the matrix r and the matrix A to obtain two corresponding mutation probability sets V and Q; and the matrix A r ;

[0042] 1.3.7) Input the two sets V and Q obtained in step 1.3.6) into the JS divergence formula defined by equation (5) to obtain a distribution value. The smaller the distribution value, the more similar the mutation rates in the extraction matrix and the matrix A r . Therefore, when the distribution value is less than or equal to 0.09, insert the matrix into the amplification matrix and update otherwise, repeat step 1.3.6);

[0043]

[0044] 1.3.8) Randomly select a sample from A in step 1.3.3 fakedata and directly add it to the amplification matrix Update

[0045] 1.3.9) If the current amplification matrix Sample size Equal to the maximum sample size m max , sample replenishment is completed; otherwise, execute step 1.3.5) to finally obtain a sample size equal to the maximum sample size m max Equal amplification matrix Then make the matrix 2) Model CDP-HA that minimizes the dispersion between the total weights of each cancer:

[0046] Suppose R, R ≥ 2 cancer types, for each cancer type, a binary somatic mutation matrix is ​​represented as It is used to record whether the gene in the sample is mutated. It has m r Row and n r Columns, rows represent samples or patients, columns represent genes, r = 1, 2, 3, .., R., a i- Represents the matrix A r The i-th sample in a -j Represents the matrix A r When the jth gene in the i-th sample in the mutation matrix of the r-th cancer mutates, otherwise Given a gene set S of size k, Represents the corresponding matrix A r Medium size is m r ×k submatrices, Represents a submatrix Samples with mutations in gene a, express The total number of samples covered in is used to measure the coverage of gene set S. The sum of samples with overlapping coverage measures the mutual exclusivity of gene set S;

[0047] According to the definition of the previous paragraph, the nonlinear maximization weight function model CDP-HA is constructed: given R types of cancer m r Row, n r The binary somatic mutation matrix A r , a parameter K, let W C (S) is the maximum weight and function, which determines an m×K submatrix The specific formula (6) is as follows:

[0048]

[0049] in represents the absolute weight value of gene set S in the rth cancer type;

[0050] 3) Introduce the parthenogenetic algorithm to solve the CDP-HA model:

[0051] 3.1) Set the fitness function:

[0052] Assume that chromosome E is given, let M E Represents a submatrix corresponding to the chromosome, matrix M E The scale is m×K, and the definition of the individual fitness function Fitness(E) is as shown in the following formula (7). The larger the value of the individual fitness function, the better the feasible solution;

[0053] Fitness(E)=W C (M E ) (7);

[0054] 3.2) Set the selection operator:

[0055] Roulette wheel selection and elite strategy are used to generate a new generation of population. The individuals with the highest fitness are directly inherited from the parent generation to the offspring, and then the roulette wheel selection operator is used to generate the remaining N-1 individuals.

[0056] 3.3) Set the reorganization operator:

[0057] A recombination operator based on a greedy strategy is adopted. The steps are as follows: First, given a parent chromosome E = {e 1 ,e 2 ,...,e k}(e i =1,2,...,n), where e i Represents a gene number, so E is also called a gene set, which determines the candidate gene set Secondly, randomly delete a gene from gene set E to obtain gene set E ′ ; Finally, based on the greedy strategy, the best gene is selected from the candidate set Right now And eventually produce new offspring

[0058] 3.4) Set parameters:

[0059] Input R cancer enhanced somatic mutation matrix A r , gene number g number , parameter k, population size N, number of algorithm executions t, maximum evolution generations maxg;

[0060] 3.5) Construct the initial population:

[0061] The chromosomes are encoded in decimal code, and one chromosome represents one individual, which is used to represent the solution vector of the problem. In the parthenogenetic algorithm, a set of K genes is used as a solution to the problem, that is, E = {e 1 ,e 2 ,...,e k}(e i =1,2,...,n), the method of initializing individuals in the population is: randomly generate natural data sets from 1 to n, each number represents a gene in the mutation matrix, group the n genes in order, and obtain n / k gene sets S 1 ,S 2 ,...,S n / k .make Select gene set S max The genes of the initial chromosomes are used to generate the initial population by generating N initial chromosomes; the first K numbers are selected as an initial chromosome to generate an initial population pop 0 , the population size is N, calculate pop 0 The fitness value of the population chromosome, pop 0 Compared with the best chromosome in the set, the best individual is saved in the variable best, and the initial iteration number step = 0;

[0062] 3.6) Perform iterative operations:

[0063] 3.6.1) If step>maxg, go to step 3.6.5) to obtain a common driving path of size K, otherwise go to step 3.6.2);

[0064] 3.6.2) For the population pop step First, pop step The best chromosome with the highest fitness value is placed in pop step+1 , and then execute the roulette wheel selection operator to select the remaining N-1 chromosomes and put them into pop step+1 ;

[0065] 3.6.3) If step<0.7*maxg or Fitness(E′)>Fitness(E), update chromosome E=E′, otherwise do not update, keep E, step=step+1;

[0066] 3.6.4) Take pop step+1 The chromosome with the highest fitness value among all the chromosomes is updated if its fitness value is greater than that of the best chromosome, that is, best = pop step+1 The optimal chromosome of

[0067] 3.6.5) The best chromosome is converted into a gene set, thereby obtaining a submatrix M, and the submatrix M is output. The output M is the common driving pathway S of size K.

[0068] This technical solution proposes to use Generative Adversarial Networks (GAN) to generate samples of somatic mutation data of small sample cancers, so that the samples between multiple cancer types are balanced, and then use the mathematical model CDP-HA of the absolute weight value of each cancer type to accumulate the harmonic mean, and finally use the parthenogenetic algorithm to solve the model. For several types of cancer, the method proposed in this technical solution can effectively supplement the small sample number of cancers and solve the problem of sample quantity differences. At the same time, the gene set identified based on the proposed model not only mutates in most of the samples of these cancers, but also the proportion of mutant samples in individual cancers is very close. In addition, this method detects some biologically significant gene sets that are missing in other methods. Therefore, it is possible to narrow the difference in sample data volume between different types of cancers based on the adversarial generative network proposed in this technical solution, which provides a new idea for identifying common driving pathways for cancer.

[0069] Compared with the existing technology, this technical solution has the following advantages:

[0070] (1) A nonlinear weight maximization function that minimizes the discreteness is designed to measure the relative weights of multiple cancer types.

[0071] (2) A data augmentation method SNV-GANs suitable for cancer somatic mutation data is proposed, which has high practical value and is the first time that GANs have been applied to cancer somatic mutation data.

[0072] (3) The pan-cancer common driver pathways found by the entire technical solution include more genes enriched in the same important signaling pathway, and the identified genes are enriched in more important signaling pathways.

[0073] This method can not only be a useful tool for enhancing biological cancer somatic mutation data and identifying pan-cancer common driver pathways, but also provide more biological information, have strong scalability and practicality, and can identify more genes enriched in important signaling pathways, as well as identify genes enriched in more important signaling pathways. BRIEF DESCRIPTION OF THE DRAWINGS

[0074] Figure 1 This is an example diagram of a common driving path in the embodiment;

[0075] Figure 2 An example diagram of a partheno-genetic algorithm in an embodiment;

[0076] Figure 3 It is the network model of SNV-GANs in the example;

[0077] Figure 4 Schematic diagram of the training process of SNV-GANs in the example;

[0078] Figure 5 Pseudo code for the SNV-GANs training process in the example;

[0079] Figure 6 This is a sample diagram of the actual effect of generating data in the example. DETAILED DESCRIPTION

[0080] The present invention will be further described in detail below in conjunction with the accompanying drawings and specific embodiments, but the present invention is not limited thereto.

[0081] Example:

[0082] In this experiment, step 1) was performed on a Linux server (Intel(R) Xeon(R) Gold 6230 2.10GHz CPU, 256G memory, 32G video memory), with Python 3.7.9 as the compilation and execution environment. Steps 2) and 3) were performed on a computer (Intel(R) Core(TM) i5-6500 3.20GHz CPU, 32G memory), with Windows 10 as the operating system, Eclipse 4.23 as the compilation and execution tool, and Java 1.8.0 as the compilation environment.

[0083] This example illustrates the problem of identifying pan-cancer common driver pathways.

[0084] A pan-cancer common driver pathway identification method based on GAN sample balance includes the following steps:

[0085] Assume there are somatic mutation data A of three cancer types r : A 1 (COADCORE has 95 samples and 211 genes), A 2 (BLCA has 95 samples and 211 genes), A 3 (LUAD has 95 samples, 211 genes).

[0086] 1) Generate somatic mutation data corresponding to cancer that conforms to the real data distribution:

[0087] 1.1) Set up the adversarial generation network framework:

[0088] Assume that BLCA cancer somatic mutation data A is used 2As an example training set, the generator network of SNV-GANs is defined as G(z). The input of the generator is: z~norm(0,1). The definition of the generator is as follows: Figure 3 As shown:

[0089] 1.1.1) The input layer uses GFC1 to map the noise vector z into a tensor zn of dimension (1,128);

[0090] 1.1.2) The hidden layer puts the tensor zn in step 1.1.1) into GFC2 for mapping, and the result is put into GFC3 for mapping. The same is true for GFC4, and finally it is mapped into a tensor zn′ of dimension (1,1024);

[0091] 1.1.3) The output layer maps the tensor zn′ to a tensor gn of dimension (1,95*211) through GFC5, and then resets gn to a tensor of dimension (95,211) Tensor is the output of the generator;

[0092] The input layer and the hidden layer both use the dropout function to freeze some neurons, and use the activation function ReLU defined by formula (1). The output layer uses the activation function Sigmod defined by formula (2). The discriminator network is defined as D(x), and the input of the discriminator is the real data x~P real Or generate data x represents a set of somatic mutation data samples, and the discriminator is defined as follows: Figure 3 As shown:

[0093] 1.1.4) The input layer uses DFC1 to map x into a tensor xn of dimension (95,256);

[0094] 1.1.5) The hidden layer puts the tensor xn in step 1.1.4) into DFC2 for mapping, and the result is put into DFC3 for mapping. The same is true for DFC4, and finally it is mapped into a tensor xn′ of dimension (95,16);

[0095] 1.1.6) The output layer maps the tensor xn′ to a tensor dn of dimension (95,1) through DFC5, which is the output of the discriminator; the input layer and the hidden layer both use the activation function ReLU defined by formula (1), and the output layer uses the activation function Sigmod defined by formula (2);

[0096] ReLU(x)=max(0,x) (1),

[0097]

[0098] 1.2) Training process of SNV-GANs:

[0099] 1.2.1) Given a somatic mutation matrix A 2 (95×211) and a random sample ratio parameter ∝ = 0.7, according to the ratio parameter ∝ from the matrix A r Randomly select a submatrix M with a sample size of m = 95*0.7≈67 2 , the size is 67×211, and a total of 64 sub-matrices M are extracted 2 Construct a training set X and input it into the generative adversarial network for training, such as Figure 4 As shown;

[0100] 1.2.2) Initialize the parameters θ of the discriminator D(.) d , the parameters θ of the generator G(.) g ;

[0101] 1.2.3) Set the current round epoch = 1, and randomly generate a 1×100 Gaussian distributed noise vector z;

[0102] 1.2.4) Take z from step 1.2.3) as the input of the generator to get a vector of size 67×211

[0103] 1.2.5) Calculate the generator loss value according to formula (3), and then update the generator parameter θ g :

[0104]

[0105] Where G(z (i) ) indicates that the generator passes the noise vector z (i) The generated data, D(G(z (i) )) represents the probability that the generator determines the generated data as real data, loss G The smaller the better;

[0106] 1.2.6) Randomly select a sample group x from the training set X;

[0107] 1.2.7) Randomly generate a 1×100 Gaussian distributed noise vector z;

[0108] 1.2.8) Take z from step 1.2.7) as the input of the generator to get a vector of size 95×211

[0109] 1.2.9) Calculate the generator loss value according to formula (4), and then update the generator parameter θ g :

[0110]

[0111] Where D(x (i) ) represents the probability that the discriminator determines the generated data as real data, 1-D(G(z (i) )) represents the probability that D determines the generated data as generated data, loss D Bigger is better;

[0112] 1.2.10) Determine whether the current epoch has reached 10,000: If so, stop training; otherwise, return to step 1.2.3) to finally obtain the trained generator G(.);

[0113] 1.3) Data processing:

[0114] 1.3.1) Randomly generate a 1×100 Gaussian distributed noise vector z;

[0115] 1.3.2) Input the vector z in step 1.3.1) into the generator G(.) obtained through training to obtain the generated data G data =G(z);

[0116] 1.3.3) G data The values ​​greater than or equal to 0.85 are set to 1, and the values ​​less than 0.85 are set to 0, and a new binary matrix A is obtained. fakedata ;

[0117] 1.3.4) Take the sample size of the cancer with the largest number of samples in the somatic mutation matrix of the three cancers as m max =m 1 =489; then matrix A 2 Insertion Amplification Matrix At this point, the amplification matrix The sample size is

[0118] 1.3.5) If the number of samples to be amplified Greater than matrix A 2 The sample size m 2 =95, proceed to step 1.3.6), if the number of samples to be amplified Smaller than matrix A 2 The sample size m 2 =95, proceed to step 1.3.8);

[0119] 1.3.6) From step 1.3.3) fakedata The number of random draws is m 2 Sample Calculate the matrix separately and the matrix A2 The mutation rate of each gene in , and two corresponding mutation probability sets V and Q are obtained respectively; 1.3.7) The two sets V and Q obtained in step 1.3.6) are input into the JS divergence formula defined by formula (5) to obtain a distribution value. The smaller the distribution value, the smaller the extraction matrix With the matrix A r The more similar the mutation rates are, the distribution value is less than or equal to 0.09, that is, the matrix Insertion Amplification Matrix and update Otherwise, repeat step 1.3.6);

[0120]

[0121] 1.3.8) From step 1.3.3) fakedata Randomly select a sample and add it directly to matrix A r Update

[0122] 1.3.9) If the current amplification matrix Sample size Equal to the maximum sample size m max =489, the sample supplementation is completed; otherwise, execute step 1.3.5) to finally obtain a sample size equal to the maximum sample size m max =489 equal matrices Get a new amplification matrix The mutation rate is similar to the original matrix A 2 The mutation rate is basically the same, and then the matrix like Figure 6 As shown; the pseudo code of the training process is as follows Figure 5 ;

[0123] 2) Model CDP-HA that minimizes the dispersion between the total weights of each cancer:

[0124] Suppose R = 3 cancer types. For each cancer type, a binary somatic mutation matrix is ​​represented as It is used to record whether the gene in the sample is mutated. It has m r Row and n r Columns, rows represent samples or patients, columns represent genes, r = 1, 2, 3, .., R., a i- Represents the matrix A r The i-th sample in a -j Represents the matrix A r When the jth gene in the i-th sample in the mutation matrix of the r-th cancer mutates, otherwise Given a gene set S of size k, Represents the corresponding matrix A r with size m r ×k sub-matrix in Represents the sub-matrix Samples in which gene a has mutated Represents The total number of samples covered in, used to measure the coverage of gene set S The sum of overlapping covered samples, measuring the mutual exclusivity of gene set S

[0125] According to the above definitions, a non-linear maximum weight function model CDP-HA is constructed: Given an m r -row, n r -column binary somatic mutation matrix A r , a parameter K (0 < K ≤ 10), let W C (S) be the maximum weight sum function, and determine an m×K sub-matrix The specific formula (6) is as follows:

[0126]

[0127] Where Represents the absolute weight value of gene set S in the r-th cancer type; as Figure 1 Shown is the mutation matrix of three cancers, with two sub-matrices S 1 and S 2 , the scale size K is 3, according to the model proposed in this example, we can get It is found that S 2 has a higher weight, and S 2 is more in line with the commonality required by the common driver pathway;

[0128] 3) Introduce a single-parent genetic algorithm to solve the model CDP-HA:

[0129] 3.1) Set the fitness function:

[0130] Assume a given chromosome E, let M E represent a sub-matrix corresponding to the chromosome, and the scale size of matrix M E is m×K. The definition of the individual fitness function Fitness(E) is shown in the following formula (7). The larger the individual fitness function value, the better the feasible solution plan;

[0131] Fitness(E) = W C (M E ) (7);

[0132] 3.2) Set the selection operator:

[0133] Roulette wheel selection and elite strategy are used to generate a new generation of population. The individuals with the highest fitness are directly inherited from the parent generation to the offspring, and then the roulette wheel selection operator is used to generate the remaining N-1 individuals.

[0134] 3.3) Set the reorganization operator:

[0135] A recombination operator based on a greedy strategy is adopted. The steps are as follows: First, given a parent chromosome E = {e 1 ,e 2 ,...,e k}(e i =1,2,...,n), where represents a gene number, so E is also called a gene set, and the candidate gene set is determined from this Secondly, randomly delete a gene from gene set E to obtain gene set E ′ ; Finally, based on the greedy strategy, the best gene is selected from the candidate set Right now And eventually produce new offspring

[0136] 3.4) Set parameters:

[0137] Input the enhanced somatic mutation matrix A of 3 cancers r , gene number g number =211, parameter k=3, population size N=20, algorithm execution times step=10, maximum evolution generation maxg=1000;

[0138] 3.5) Construct the initial population:

[0139] The chromosomes are encoded in decimal, and one chromosome represents one individual, which is used to represent the solution vector of the problem. In the one-parent genetic algorithm, a set of K = 3 genes is used as a problem solution, that is, E = {e 1 ,e 2 ,...,e k}(e i =1,2,...,n), the method of initializing individuals in the population is: randomly generate 1 to 20 natural data sets, each number represents a gene in the mutation matrix, group the 20 genes in order, and obtain n / k=20 / 3≈7 gene sets S 1 ,S 2 ,...,S n / k .make Select gene set S max The genes of the initial chromosomes are used to generate the initial population by generating N initial chromosomes; the first K = 3 numbers are selected as an initial chromosome to generate an initial population pop0 , the population size is N, calculate pop 0 The fitness value of the population chromosome, pop 0 Compared with the best chromosome in , the best individual is saved in the variable best, and the initial iteration number step = 0, such as Figure 2 As shown;

[0140] 3.6) Perform iterative operations:

[0141] 3.6.1) If step>maxg, go to step 3.6.5) to obtain a common driving path of size K, otherwise go to step 3.6.2);

[0142] 3.6.2) For the population pop step First, pop step The best chromosome with the highest fitness value is placed in pop step+1 , and then execute the roulette wheel selection operator to select the remaining N-1=20-1=19 chromosomes and put them into pop step+1 ; 3.6.3) If step<700 or Fitness(E′)>Fitness(E), update chromosome E=E′, otherwise do not update, keep E, step=step+1;

[0143] 3.6.4) Take pop step+1 If the chromosome with the highest fitness value is greater than the fitness value of the best chromosome, the best chromosome is updated, that is, best = pop step+1 The optimal chromosome of

[0144] 3.6.5) The best chromosome is converted into a gene set, thereby obtaining a submatrix M, and the submatrix M is output. The output M is a common driving pathway S of size K=3.

Claims

1. A pan-cancer common driver pathway identification method based on GAN sample balance, characterized by comprising the following steps: 1) Generate somatic mutation data corresponding to cancer that conforms to the real data distribution: 1.1) Set up the adversarial generation network framework: Assume a system with m r samples, n r The example training set of genes, the generator network of SNV-GANs is defined as G(z), the input of the generator is: z~norm(0,1), and the definition of the generator is as follows: 1.1.1) The input layer uses GFC1 to map the noise vector z into a tensor zn of dimension (1,128); 1.1.2) The hidden layer puts the tensor zn in step 1.1.1) into GFC2 for mapping, and the result is put into GFC3 for mapping. The same is true for GFC4, and finally it is mapped into a tensor zn′ of dimension (1,1024); 1.1.3) The output layer maps the tensor zn′ to dimension (1,m) through GFC5 r *n r ) tensor gn, and then reset gn to a tensor of dimension (m r ,n r ) Tensor is the output of the generator; The input layer and the hidden layer both use the dropout function to freeze some neurons, and use the activation function ReLU defined by formula (1). The output layer uses the activation function Sigmod defined by formula (2). The discriminator network is defined as D(x), and the input of the discriminator is the real data x~P real Or generate data x represents a set of somatic mutation data samples, and the discriminator is defined as follows: 1.1.4) The input layer uses DFC1 to map x to dimension (m r ,256) tensor xn; 1.1.5) The hidden layer puts the tensor xn from step 1.1.4) into DFC2 for mapping, and the result is put into DFC3 for mapping. The same is true for DFC4, and the final mapping is dimension (m r ,16) tensor xn′; 1.1.6) The output layer maps the tensor xn′ to dimension (m) through DFC5 r ,1) tensor dn, which is the output of the discriminator; The input layer and hidden layer both use the activation function ReLU defined by formula (1), and the output layer uses the activation function Sigmod defined by formula (2); ReLU(x)=max(0,x) (1), 1.2) Training process of SNV-GANs: 1.2.1) Given a somatic mutation matrix A r (m r ×n r ) and a random sample ratio parameter ∝, ∝ < 1, according to the ratio parameter ∝ from the matrix A r The number of randomly selected samples is m, m = m r *∝'s submatrix M r , size is m×n r , a total of 64 sub-matrices M are extracted r Construct a training set X and input it into the generative adversarial network for training; 1.2.2) Initialize the parameters θ of the discriminator D(.) d , the parameters θ of the generator G(.) g ; 1.2.3) Set the current round epoch = 1, and randomly generate a 1×100 Gaussian distributed noise vector z; 1.2.4) Take z from step 1.2.3) as the input of the generator to get a size of m×n r Vector 1.2.5) Calculate the generator loss value according to formula (3), and then update the generator parameter θ g : Where G(z (i) ) indicates that the generator passes the noise vector z (i) The generated data, D(G(z (i) )) represents the probability that the discriminator determines the generated data as real data, loss G The smaller the better; 1.2.6) Randomly select a sample group x from the training set X; 1.2.7) Randomly generate a 1×100 Gaussian distributed noise vector z; 1.2.8) Take z from step 1.2.7) as the input of the generator to get a size of m×n r Vector 1.2.9) Calculate the generator loss value according to formula (4), and then update the generator parameter θ g : Where D(x (i) ) represents the probability that the discriminator determines the generated data as real data, 1-D(G(z (i) )) represents the probability that D determines the generated data as generated data, loss D Bigger is better; 1.2.10) Determine whether the current epoch has reached the set maximum epoch: if so, stop training; otherwise, return to step 1.2.3) to finally obtain the trained generator G(.); 1.3) Data processing: 1.3.1) Randomly generate a 1×100 Gaussian distributed noise vector z; 1.3.2) Input the vector z in step 1.3.1) into the generator G(.) obtained through training to obtain the generated data G data =G(z); 1.3.3) G data The values ​​greater than or equal to 0.85 are set to 1, and the values ​​less than 0.85 are set to 0, and a new binary matrix A is obtained. fakedata ; 1.3.4) Take the sample size of the cancer with the largest number of samples in the somatic mutation matrix of R cancers as m max , 0 < max < R; Then insert matrix A r into the amplification matrix At this time, the amplification matrix has a sample quantity of 1.3.5) If the number of samples to be amplified Greater than matrix A r The sample size m r , proceed to step 1.3.6), if the number of samples to be amplified Smaller than matrix A r The sample size m r , execute step 1.3.8); 1.3.6) From step 1.3.3) fakedata The number of random draws is m r Sample Calculate the matrix separately and the matrix A r The mutation rate of each gene in , and two corresponding mutation probability sets V and Q are obtained respectively; 1.3.7) Input the two sets V and Q obtained in step 1.3.6) into the JS divergence formula defined by formula (5) to obtain a distribution value. The smaller the distribution value, the smaller the extraction matrix With the matrix A r The more similar the mutation rates are, the distribution value is less than or equal to 0.09, that is, the matrix Insertion Amplification Matrix and update Otherwise, repeat step 1.3.6); 1.3.8) From step 1.3.3) fakedata Randomly select a sample and add it directly to the amplification matrix Update 1.3.9) If the current amplification matrix Sample size Equal to the maximum sample size m max , sample replenishment is completed; otherwise, execute step 1.3.5) to finally obtain a sample size equal to the maximum sample size m max Equal amplification matrix Then make the matrix 2) Model CDP-HA that minimizes the dispersion between the total weights of each cancer: Suppose R, R ≥ 2 cancer types, for each cancer type, a binary somatic mutation matrix is ​​represented as It is used to record whether the gene in the sample is mutated. It has m r Row and n r Columns, rows represent samples or patients, columns represent genes, r = 1, 2, 3, .., R., a i- Represents the matrix A r The i-th sample in a -j Represents the matrix A r When the jth gene in the i-th sample in the mutation matrix of the r-th cancer mutates, otherwise Given a gene set S of size k, Represents the corresponding matrix A r Medium size is m r ×k submatrices, Represents a submatrix Samples with mutations in gene a, express The total number of samples covered in is used to measure the coverage of gene set S. The sum of samples with overlapping coverage measures the mutual exclusivity of gene set S; According to the definition of the previous paragraph, the nonlinear maximization weight function model CDP-HA is constructed: given R types of cancer m r Row, n r The binary somatic mutation matrix A r , a parameter K, let W C (S) is the maximum weight and function, which determines an m×K submatrix The specific formula (6) is as follows: in represents the absolute weight value of gene set S in the rth cancer type; 3) introducing a parthenogenetic algorithm to solve the model CDP-HA: 3.1) Set the fitness function: Assume that chromosome E is given, let M E Represents a submatrix corresponding to the chromosome, matrix M E The scale is m×K, and the definition of the individual fitness function Fitness(E) is as shown in the following formula (7). The larger the value of the individual fitness function, the better the feasible solution; Fitness(E)=W C (M E ) (7); 3.2) Set the selection operator: Roulette wheel selection and elite strategy are used to generate a new generation of population. The individuals with the highest fitness are directly inherited from the parent generation to the offspring, and then the roulette wheel selection operator is used to generate the remaining N-1 individuals. 3.3) Set the reorganization operator: A recombination operator based on a greedy strategy is adopted. The steps are as follows: First, given a parent chromosome E = {e1, e2, ..., e k },e i =1,2,...,n, where e i Represents a gene number, so E is also called a gene set, which determines the candidate gene set Secondly, randomly delete a gene from gene set E to obtain gene set E ′ ; Finally, based on the greedy strategy, the best gene is selected from the candidate set Right now And eventually produce new offspring 3.4) Set parameters: Input R cancer enhanced somatic mutation matrix A r , gene number g number , parameter k, population size N, number of algorithm executions t, maximum evolution generations maxg; 3.5) Construct the initial population: The chromosome is encoded in decimal. One chromosome represents one individual and is used to represent the solution vector of the problem. In the parthenogenetic algorithm, a set of K genes is used as a solution to the problem, that is, E = {e1, e2, ..., e k },e i =1,2,...,n, the method of initializing individuals in the population is: randomly generate natural data sets from 1 to n, each number represents a gene in the mutation matrix, group the n genes in order, and obtain n / k gene sets S1,S2,...,S n / k ,make Select gene set S max The genes of constitute the initial chromosome, and the initial population is generated by generating N initial chromosomes; the first K numbers are selected as an initial chromosome to generate an initial population pop0 with a population size of N, the fitness value of the chromosome of the pop0 population is calculated, the best chromosome in pop0 is compared, and the best individual is saved in the variable best. The initial iteration number step = 0; 3.6) Perform iterative operations: 3.6.1) If step>maxg, go to step 3.6.5) to obtain a common driving path of size K, otherwise go to step 3.6.2); 3.6.2) For the population pop step First, pop step The best chromosome with the highest fitness value is placed in pop step+1 , and then execute the roulette wheel selection operator to select the remaining N-1 chromosomes and put them into pop step+1 ; 3.6.3) If step<0.7*maxg or Fitness(E′)>Fitness(E), update chromosome E=E′, otherwise do not update, keep X, step=step+1; 3.6.4) Take pop step+1 If the chromosome with the highest fitness value is greater than the fitness value of the best chromosome, the best chromosome is updated, that is, best = pop step+1 The optimal chromosome of 3.6.5) The best chromosome is converted into a gene set, thereby obtaining a submatrix M, and the submatrix M is output. The output M is the common driving pathway S of size K.

Images