Peripheral gene regulation feature extraction method based on graph representation learning and attention

Abstract

The present invention discloses a method for extracting peripheral gene regulation features based on graph representation learning and attention mechanism, which solves the problem that the existing method cannot extract the diffusion pattern of peripheral genes to core gene modules in the network, and cannot reflect the biological regulation pattern of peripheral genes to any core gene. Its technical scheme is: determine the core gene set and the peripheral gene set in the construction of the interaction network, sample and splice the energy of the peripheral genes diffused to the core genes in the interaction network obtained by the graph representation learning algorithm to form topological regulation features, construct a graph attention neural network, and extract the biological regulation features of peripheral genes to each core gene. The present invention extracts the topological regulation features of peripheral genes through the graph representation learning method, builds and trains an attention neural network, so that the present invention can have a higher accuracy in the sub-classification task of samples.

Description

Technical Field

[0001] The present invention belongs to the field of physical technology, and further relates to a peripheral gene regulatory feature extraction method based on graph representation learning and attention mechanism in the field of information technology for genetic or protein-related data processing in computational molecular biology. The present invention can extract the topological regulatory features of peripheral genes in the interaction network through a graph representation learning algorithm, and then construct a graph attention neural network to extract the biological regulatory features of peripheral genes according to the biological functions involved in the peripheral genes. Background Art

[0002] At present, with the development of technology and the improvement of the understanding of traits, the whole genome graph model is a new exploration of the genetic structure of complex phenotypes. According to the relationship between the variant site and the phenotype, the genes are divided into core genes and peripheral genes. The core genes directly regulate the phenotype, and the peripheral genes are indirectly regulated by superimposing on the core genes through different trans-regulations of the phenotype. However, the contribution of core genes to the total heritability of the disease is generally small, while most peripheral genes expressed in related cell types contribute a large part to the heritability. The superposition of the weak effects of many peripheral genes plays a more critical regulatory role in the phenotype. Nowadays, phenotypes are analyzed based on the whole genome graph model, but the focus is on analyzing the genetic structure, energy control, etc. of the phenotype, and there is no mining of graph patterns to explain its peripheral effects. From the perspective of peripheral effects, new insights can be provided for understanding phenotypes, and quantifying how peripheral genes regulate core genes is a key step in understanding peripheral effects. With the rapid development of graph representation learning methods and machine learning, their applications in biology have also become mature. Graph representation learning methods can obtain feature representations of any node in the network based on biological molecular networks, while machine learning methods often require a large amount of sample training when performing prediction tasks, but are not suitable for few-sample prediction tasks and have poor interpretability.

[0003] Claire Donnat et al. published a method for extracting unsupervised peripheral gene regulation features based on structural similarity in their paper "Learning Structural Node Embeddings via Diffusion Wavelets, (in Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining. 2018, Association for Computing Machinery: London, United Kingdom. p. 1320–1329.)". The implementation steps of this method are: Step 1, calculate the Laplace matrix of the input network; Step 2, obtain the diffusion pattern of any node for the entire network node according to the spectral wavelet; Step 3, embed the spectral wavelet coefficient distribution into two-dimensional space by calculating the characteristic function of each node; Step 4, unsupervised uniform sampling of nodes in the network and splicing the values ​​to obtain the structural characteristics of the nodes. The shortcoming of this method is that although it uses a large number of mathematical formulas to illustrate that nodes with similar roles in any network can be extracted as similar features, which can be well applied to classification tasks, this feature is oriented to the entire network and cannot obtain the regulatory characteristics of peripheral genes on a core gene module.

[0004] Tianle Ma et al. published a method for extracting peripheral gene regulatory features based on deep learning in their paper "Affinitynet: semi-supervised few-shot learning for disease type prediction. (In Proceedings of the AAAI conference on artificial intelligence, vol. 33, no. 01, pp. 1069-1076. 2019)". The implementation steps of this method are: step 1, generating a sample training set; step 2, constructing a feature selection network to screen key features; step 3, constructing a K-nearest neighbor attention pooling layer to aggregate similar features; step 4, using aggregated features to establish a sample similarity graph; step 5, using a small number of labeled samples to train the model and identify samples. The disadvantage of this method is that although it can extract important features with a small number of training samples, it can be used in biomedicine to distinguish disease samples well, but this feature can only extract key features for distinguishing peripheral genes, and cannot reflect the biologically significant regulation of peripheral genes on any core gene. Summary of the invention

[0005] The purpose of the present invention is to address the deficiencies of the above-mentioned prior art and propose a method for extracting peripheral gene regulatory features based on graph representation learning and attention mechanism, so as to solve the problem that the existing method cannot extract the regulatory features of peripheral genes on a core gene module in the network, and the extracted features cannot reflect the regulation of peripheral genes on any core gene in a biological sense.

[0006] The idea of ​​achieving the purpose of the present invention is that the present invention extracts the regulatory features of peripheral genes from two aspects: topological structure and biological function. From the aspect of topological structure, the topological regulatory features of peripheral genes in the gene interaction network are extracted based on the graph representation learning method. The topological regulatory features reflect the diffusion pattern of peripheral genes in the network towards the core gene module, thereby solving the problem that the existing methods cannot extract the regulatory features of peripheral genes on a core module in the network. From the aspect of biological function, the present invention constructs a graph attention neural network to extract the biological regulatory features of peripheral genes on each core gene. The biological regulatory features include the biological regulatory pattern of peripheral genes on each core gene, thereby solving the problem that the features extracted by the existing methods cannot reflect the regulation of peripheral genes on each core gene in a biological sense.

[0007] The implementation steps of the present invention are as follows:

[0008] Step 1: Determine the core genes and peripheral genes in the constructed network:

[0009] Step 1.1, construct a human gene interaction network consisting of N nodes, each network node in the human gene interaction network corresponds to a gene, establish an edge between every two genes that interact, and model the network with an undirected graph G = (V, E), where N represents the total number of gene nodes in the human gene interaction network, V represents the total number of gene nodes in the undirected graph G, the value of V is equal to N, and E represents the total number of edges in the undirected graph G;

[0010] Step 1.2, add up all the regulatory values ​​of each gene in the human gene interaction network to obtain the regulatory factor of the gene on the phenotype, select genes with regulatory factors greater than or equal to 4 on the phenotype in the human gene interaction network to form the core gene set, and select all gene nodes with regulatory coefficients of core genes on the phenotype less than or equal to 0.05 in the human gene interaction network to form the peripheral gene set;

[0011] Step 2: Extract the topological regulatory features of each gene in the peripheral gene set on the core gene set:

[0012] Step 2.1, calculate the energy of each peripheral gene node propagating to the core gene set according to the following formula:

[0013]

[0014] Among them, δ j represents the energy propagated from the jth gene node in the peripheral gene set to the core gene set, C represents the size of the core gene set, ∑ represents the summation operation, r represents the sequence number of the core gene in the core gene set, exp(·) represents the exponential operation with the natural constant e as the base, i represents the symbol of the imaginary part of the complex number, and λ o represents the oth eigenvalue in the eigenvalue diagonal matrix decomposed by the Laplace matrix, o = 1, 2, ..., O, O represents the number of eigenvalues ​​in the eigenvalue diagonal matrix decomposed by the Laplace matrix, s represents the diffusion factor, the value of s is equal to the average value of all eigenvalues ​​in the eigenvalue diagonal matrix decomposed by the Laplace matrix, U rl represents the energy value propagated from the lth gene node in the eigenvector matrix of Laplace matrix decomposition to the rth gene node in the core gene set, U lj represents the energy value propagated from the jth gene node in the eigenvector matrix decomposed by the Laplace matrix to the lth gene node in the peripheral gene set, l = 1, 2, ..., L, L represents the size of the peripheral gene set, and the value of l is equal to that of o;

[0015] Step 2.2, sample the propagation of each gene in the peripheral gene set to each core gene in the core gene set to obtain the topological regulation characteristics of the gene node:

[0016] γ j =concat(Re(δ j ),Im(δ j ))

[0017] Among them, γ j represents the topological regulatory features of the j-th gene node in the peripheral gene set, Re(·) represents the real part operation, Im(·) represents the imaginary part operation, and concat(·) represents the concatenation operation;

[0018] Step 3, generate the initial selection vector for each gene node;

[0019] Step 3.1, select q biological function entries, each of which contains the coding gene nodes involved in the biological function; the total number of coding genes contained in the kth biological function entry is n kThe k-th position of each coding gene is set to 1, and the k-th position of the coding gene not included in the k-th biological function entry is set to 0, thereby representing the situation that each coding gene in the node set V participates in q biological function entries as a p-dimensional initial biological feature vector, and the initial selection vector of each coding gene node is obtained by mean normalizing the initial biological feature of each coding gene node, where p = q, q represents the number of all biological function entries;

[0020] Step 3.2, select the set of coding gene nodes X from the set of first-order neighbor gene nodes of each non-coding gene node in the undirected graph G; k The intersection of the coding gene set and the coding gene node set X is taken to obtain x genes; the number of randomly selected genes in the set of gene nodes is n k The intersection of the gene set and the coding gene node set X is x. r The overlapping node sets are randomly selected 1000 times to obtain 1000 overlapping node sets. The mean μ and variance σ of the selected 1000 overlapping node sets are calculated. Calculate the significance score z_score. If z_score is greater than 1.645, set the kth position of the non-coding gene node to 1. If z_score is less than 1.645, set the kth position of the non-coding gene node to 0. Thus, each non-coding gene is represented as a p-dimensional initial biological feature vector. The initial biological feature of each non-coding gene node is mean-normalized to obtain the initial selection vector of the non-coding gene node.

[0021] Step 4: Generate training set:

[0022] Step 4.1, set all peripheral genes in the peripheral gene set as positive samples, and set the label of the positive sample to 0; set the remaining genes in the gene node set excluding the peripheral gene set and the core gene set as negative samples, and set the label of the negative sample to 1;

[0023] Step 4.2, randomly select 20% of positive samples and 20% of negative samples and their corresponding labels from the set of gene nodes after setting sample labels to form a training set, and remove samples other than the training set samples and their corresponding labels from the set of gene nodes after setting sample labels to form a test set;

[0024] Step 4.3, calculate the feature average weight vector of the initial selection vector of all gene nodes in the training set according to the following formula:

[0025]

[0026] in, represents a feature average weight vector composed of all positive samples in the training set, F train represents the total number of samples in the training set, h j Represents the initial selection vector of the j-th gene node in the training set;

[0027] Step 5: Build a graph attention neural network to extract the biological regulation features of peripheral genes:

[0028] Step 5.1: Build a 5-layer graph attention neural network to extract the biological regulatory features of peripheral genes. Its structure is as follows: feature selection module, regulatory transformation group, splicing layer, KNN-graph attention module, output layer; set the number of neurons in the feature selection module to p s , p s The value of is equal to the dimension of the initial selection vector of the sample, the output layer is implemented by a K-nearest neighbor classifier, and the parameters of the K-nearest neighbor classifier are set to 64*2;

[0029] Step 5.2: The control conversion layer is composed of 20 control conversion layers with the same structure connected in parallel. The structure of each control conversion layer is: fully connected layer, first activation layer, conversion layer, second activation layer. The parameters of the fully connected layer are set to The first activation layer is implemented using the softmax function, and the number of neurons in the conversion layer is set to 32. The second activation layer is implemented using the sigmoid function, where: The value of is equal to the dimension of the initial selection vector of the sample;

[0030] Step 6: Train the graph attention neural network:

[0031] Step 6.1, input each sample in the training set into the graph attention neural network in turn, and output 20 regulatory features of the sample through 20 parallel regulatory conversion layers in the splicing layer of the graph attention neural network. After splicing the 20 regulatory features, the biological regulatory features of the sample are obtained, and the Euclidean distance between the biological regulatory features of every two samples in all samples is calculated. The dis samples with the smallest Euclidean distance to each sample are used as the neighbor set of the sample, thereby constructing an undirected graph G'=(V',E'); the undirected graph G' and the biological regulatory features of each sample are simultaneously input into the KNN attention pooling layer; wherein V' represents all sample nodes, and E' represents the set of edges established between each sample and each sample in its neighbor set;

[0032] Step 6.2, using the stochastic gradient descent method, iteratively update the parameters until the loss function converges to obtain the trained graph attention neural network;

[0033] Step 7: Extract the regulatory features of each peripheral gene:

[0034] Step 7.1, using the same method as step 2, extracting the topological regulatory features of each gene in the peripheral gene set whose topological regulatory features are to be extracted;

[0035] In step 7.2, the peripheral genes whose biological regulatory features are to be extracted are processed using the same method as in step 4, and the processed peripheral genes are input into the trained graph attention neural network to output the biological regulatory features of the peripheral genes; the biological regulatory features of the peripheral genes are concatenated with the topological regulatory features of the corresponding peripheral genes to obtain the regulatory features of the peripheral genes.

[0036] Compared with the prior art, the present invention has the following advantages:

[0037] First, since the present invention proposes a graph representation learning method for extracting the regulatory features of peripheral genes, it overcomes the defect that the existing methods cannot extract the regulatory features of peripheral genes in the network for a core module. The present invention can better reflect the regulatory form of peripheral genes on core genes in the topological structure for the core genes determined in any constructed network, thereby obtaining more accurate classification results.

[0038] Second, because the present invention utilizes the construction and training of a graph attention neural network to extract the biological regulatory features of peripheral genes, it overcomes the defect that the features extracted by the existing methods cannot reflect the regulation of peripheral genes on any core gene in a biological sense. After the biological regulatory features of peripheral genes on any core gene are extracted, the present invention can better reflect the regulatory mechanism of peripheral genes on core genes in the whole gene graph model, and can enhance the interpretability of the peripheral gene regulation process.

[0039] Third, since the present invention extracts the regulatory characteristics of peripheral genes on core genes in topological structure and biological function based on the whole genome model, it overcomes the defect that the application of the whole genome graph model in the existing methods is still in the conceptualization stage. The present invention can extract the regulatory characteristics of peripheral genes for any phenotype based on the whole genome graph model, thereby better supporting the phenotypic sub-classification problem. BRIEF DESCRIPTION OF THE DRAWINGS

[0040] Figure 1 is a flow chart of the present invention;

[0041] Figure 2 It is a result diagram of the simulation experiment of the present invention. DETAILED DESCRIPTION

[0042] The present invention is further described in detail below in conjunction with the accompanying drawings and embodiments.

[0043] Reference Figure 1, the implementation steps of the present invention are further described.

[0044] Step 1: Identify peripheral genes and core genes in the constructed network.

[0045] Step 1.1, construct human interaction network.

[0046] The network used in the embodiment of the present invention is a human gene interaction network composed of 24215 genes. Each network node in the human gene interaction network corresponds to a gene, and an edge is established between every two genes that interact with each other.

[0047] All genes in the interaction network are collected from eight databases, including RNAInter, miRecords, LncRNADisease, miRTarBase, BIOGRID, OncoBase, LncACTdb, and PPI, as nodes in the network. An edge is established between every two genes that interact with each other. Finally, a genetic interaction network (GenicInteractome Network, GIN) with 24,215 nodes, 314,748 edges, and no isolated nodes is obtained, which is composed of multiple node types.

[0048] The data sources in the interaction network refer to:

[0049] In their paper "RNAInter in 2020: RNA interactome repository with increased coverage and annotation. (Nucleic Acids Research, Oxford University Press, 2020, 48(D1): D189–D197.)", Yunqing Lind et al. disclosed 26,568 interaction edges supported by 67 strong experimental evidences, involving a total of 10,530 genes.

[0050] Feifei Xiao et al. published a paper titled “miRecords: An integrated resource for microRNA-target interactions. (Nucleic Acids Research, 2009, 37 (SUPPL. 1): 105–110.)” that disclosed 1,992 miRNA-target interaction links, covering a total of 1,356 miRNA genes.

[0051] Geng Chen et al. published a paper titled “LncRNADisease: A database for long-non-coding RNA-associated diseases. (Nucleic Acids Research, 2013, 41(D1): 983–986.)” which revealed 148 RNA interaction links, including 177 lncRNA genes in total.

[0052] Hsi-Yuan Huang et al. published a paper titled "MiRTarBase 2020: Updates to the experimentally validated microRNA-target interaction database[J].(Nucleic Acids Research,Oxford University Press,2020,48(D1):D148–D154.)" which disclosed 8489 miRNA-target interaction links supported by strong experimental evidence, covering a total of 3589 genes.

[0053] Andrew Chatr-aryamontri et al. published 8,335 human-related genetic interaction edges in their paper “The BioGRID interaction database: 2015 update [J]. (Nucleic Acids Research, 2015, 43 (D1): D470–D478.)”, which contain a total of 3,303 coding gene nodes.

[0054] Yun Zhu et al. published a paper titled “Constructing 3D interaction maps from 1Depigenomes[J]. (Nature Communications, Nature Publishing Group, 2016, 7: 1–11.)” that disclosed 65,578 promoter-enhancer edges, which contained a total of 22,820 gene nodes.

[0055] Peng Wang WANG et al. published a paper titled “LncACTdb 2.0: an updated database of experimentally supported ceRNA interactions curated from low-and high-throughput experiments[J]. (Nucleic Acids Research, Oxford University Press, 2019, 47(D1): D121–D127.)” in which they disclosed 2,681 human gene interaction edges, involving a total of 1,668 coding genes.

[0056] Feixiong Cheng et al. published a paper titled “Network-based prediction of drug combinations[J].(Nature Communications,2019,10(1).)” which disclosed 234,714 human protein interaction edges, involving a total of 16,348 coding genes.

[0057] Step 1.2, determine the core genes and peripheral genes in the constructed interaction network.

[0058] Step 1.2.1, each network node in the human gene interaction network corresponds to a gene, and an edge is established between every two genes that interact with each other, and the model is modeled using an undirected graph G = (V, E), where V = {v 1 ,v 2 ,…,v N} represents the set of all gene nodes in the gene interaction network, E represents the set of all edges between two gene nodes in the gene interaction network, N represents the total number of network nodes in the human gene interaction network, and in the embodiment of the present invention, N=24215, E=314748.

[0059] Step 1.2.2, the whole gene graph model proposes core genes and peripheral genes, and different phenotypes have different core genes and peripheral genes. The peripheral gene regulation feature extraction method described in the present invention is to extract the regulatory features of the peripheral genes determined for each phenotype. The present invention takes the schizophrenia phenotype as an application example of the peripheral gene regulation feature extraction method, and adds all the regulatory values ​​of each gene in the human gene interaction network in 8 databases to obtain the regulatory factor of the gene on the phenotype, and selects genes with a regulatory factor greater than or equal to 4 for the phenotype to form a core gene set Core, and selects all gene nodes with a regulatory coefficient of the core gene of the phenotype less than or equal to 0.05 in the human gene interaction network to form a peripheral gene set Per, wherein the size of the core gene set in the embodiment of the present invention is 20, and the size of the peripheral gene set is Peripheral 7477.

[0060] Step 2: Extract the topological regulatory features of each peripheral gene in the peripheral gene set on the core gene set.

[0061] Step 2.1, calculate the energy of each peripheral gene node propagating to the core gene set according to the following formula:

[0062]

[0063] Among them, δ j represents the energy propagated from the jth gene node in the peripheral gene set to the core gene set, C represents the size of the core gene set, ∑ represents the summation operation, r represents the sequence number of the core gene in the core gene set, exp(·) represents the exponential operation with the natural constant e as the base, i represents the symbol of the imaginary part of the complex number, and λ o represents the oth eigenvalue in the eigenvalue diagonal matrix decomposed by the Laplace matrix, o = 1, 2, ..., O, O represents the number of eigenvalues ​​in the eigenvalue diagonal matrix decomposed by the Laplace matrix, s represents the diffusion factor, and setting s too small or too large will make the diffusion energy unclear. The value of s is equal to the average value of all eigenvalues ​​in the eigenvalue diagonal matrix decomposed by the Laplace matrix, U rl represents the energy value propagated from the lth gene node in the eigenvector matrix of Laplace matrix decomposition to the rth gene node in the core gene set, U lj represents the energy value propagated from the j-th gene node in the eigenvector matrix decomposed by the Laplace matrix to the l-th gene node in the peripheral gene set, l=1, 2, ..., L, L represents the size of the peripheral gene set, the value of l is equal to that of o, and L=7477 in the embodiment of the present invention.

[0064] Step 2.2, sample the propagation of each gene in the peripheral gene set to each core gene in the core gene set to obtain the topological regulation characteristics of the gene node:

[0065] γ j =concat(Re(δ j ),Im(δ j ))

[0066] Among them, γ j represents the topological regulatory features of the j-th gene node in the peripheral gene set, Re(·) represents the real part operation, Im(·) represents the imaginary part operation, and concat(·) represents the concatenation operation.

[0067] Step 3: construct the initial selection vector for each gene node.

[0068] Step 3.1, Omer Weissbrod et al. published a paper titled “Functionally informed fine-mapping and polygenic localization of complex trait heritability[J].(Nature Genetics,2020:1-9.)” in which q biological function entries are disclosed. Each biological function entry contains coding gene nodes involved in the biological function. In the embodiment of the present invention, the total number of coding genes contained in the kth biological function entry is n. k The k-th position of each coding gene is set to 1, and the k-th position of the coding gene not included in the k-th biological function entry is set to 0. Thus, the situation that each coding gene in the node set V participates in q biological function entries is represented as a p-dimensional initial biological feature vector. The initial selection vector f of the coding gene node is obtained by mean normalizing the initial biological feature of each coding gene node, where represents a vector space, p=q, q represents the number of all biological function entries, and in the embodiment of the present invention, q=1555, p=1555.

[0069] Step 3.2, in the present example, in the undirected graph G, a coding gene node set X is selected from the first-order neighbor gene node set of each non-coding gene node; the number of nodes n contained in the kth biological function entry is k The intersection of the coding gene set and the coding gene node set X is taken to obtain x genes. Then the number of randomly selected genes in the set of gene nodes is n k The intersection of the gene set and the coding gene node set X is x. r The overlapping node sets are randomly selected 1000 times to obtain 1000 overlapping node sets. The mean μ and variance σ of the selected 1000 overlapping node sets are calculated. The significance score z_score is calculated. If z_score is greater than 1.645, the kth position of the non-coding gene node is set to 1. If z_score is less than 1.645, the kth position of the non-coding gene node is set to 0. Thus, each non-coding gene is represented as an initial biological feature vector of p' dimension. The initial biological feature of each non-coding gene node is mean-normalized to obtain the initial selection vector of the non-coding gene node. In the embodiment of the present invention, p'=1555.

[0070] Step 4: Generate training set and test set.

[0071] Step 4.1, set all peripheral genes in the peripheral gene set as positive samples, and set the label of the positive sample to 0; set the remaining genes in the set of gene nodes excluding the peripheral gene set and the core gene set as negative samples, and set the label of the negative sample to 1.

[0072] Step 4.2, randomly select 20% of positive samples and 20% of negative samples and their corresponding labels from the set of gene nodes after setting sample labels to form a training set, and remove samples other than the training set samples and their corresponding labels from the set of gene nodes after setting sample labels to form a test set; Step 4.3, calculate the feature average weight vector of the initial selection vector of all gene nodes in the training set according to the following formula:

[0073]

[0074] in, Represents a set consisting of all positive samples in the training set The feature average weight vector, F train represents the total number of samples in the training set, h j represents the initial selection vector of the j-th gene node in the training set, and the initial selection vector is equal to p. In the embodiment of the present invention

[0075] Step 5: Construct a graph attention neural network to extract the biological regulatory features of peripheral genes.

[0076] Build a 5-layer graph attention neural network to extract the biological regulatory features of peripheral genes. Its structure is as follows: feature selection module, regulatory transformation group, splicing layer, KNN-graph attention module, output layer. Set the number of neurons in the feature selection module to p s , where p s The value of is equal to the dimension of the initial selection vector of the sample. In the embodiment of the present invention, p s = 1555. The output layer is implemented by a K-nearest neighbor classifier, and the parameters of the K-nearest neighbor classifier are set to 64*2.

[0077] The control conversion group is composed of 20 control conversion layers with the same structure connected in parallel.

[0078] The structure of each control conversion layer is: fully connected layer, first activation layer, conversion layer, second activation layer. The parameters of the fully connected layer are set to 1555*32, the first activation layer is implemented using the softmax function, the number of neurons in the conversion layer is set to 32, and the second activation layer is implemented using the sigmoid function.

[0079] The KNN-graph attention module is composed of a KNN attention pooling layer and a feed-forward layer connected in series. The feed-forward layer is implemented using the perceptron affine kernel method.

[0080] The structure of the KNN attention pooling layer is: convolution layer, attention layer, pooling layer. The convolution kernel size of the convolution layer is set to 1×1, the attention layer is implemented using the activation function sigmoid, and the pooling layer uses average pooling.

[0081] An embodiment of the present invention adopts the structure and parameter setting of the KNN attention pooling layer disclosed in the paper “Affinitynet: semi-supervised few-shot learning for disease type prediction. (In Proceedings of the AAAI conference on artificial intelligence, vol. 33, no. 01, pp. 1069-1076. 2019)” published by Tianle Ma et al.

[0082] Step 6: Train the graph attention neural network.

[0083] Step 6.1, input each sample in the training set into the graph attention neural network in turn, and output 20 regulatory features of the sample through 20 parallel regulatory conversion layers in the splicing layer of the graph attention neural network. After splicing the 20 regulatory features, the biological regulatory features of the sample are obtained, and the Euclidean distance between the biological regulatory features of every two samples in all samples is calculated. The dis samples with the smallest Euclidean distance to each sample are used as the neighbor set of the sample, thereby constructing an undirected graph G'=(V',E'); the undirected graph G' and the biological regulatory features of each sample are simultaneously input into the KNN attention pooling layer, where V' represents all sample nodes, and E' represents the set of edges established between each sample and each sample in its neighbor set. In the embodiment of the present invention, dis=5.

[0084] Step 6.2, use the stochastic gradient descent method to iteratively update the parameters until the loss function converges to obtain the trained graph attention neural network.

[0085] The loss function is:

[0086]

[0087] Among them, L represents the loss function, y j represents the true label of the j-th gene node in the training set, log represents the logarithm operation with base 10, and p j represents the probability that the j-th gene node in the training set is predicted to be a peripheral gene, ||·|| represents the modulo operation, Represents the sum of the initial selection vectors of all positive samples in the training sample.

[0088] Step 7: Extract the regulatory features of each peripheral gene.

[0089] Step 7.1, using the same method as step 2, extract the topological regulatory features of each gene in the peripheral gene set whose topological regulatory features are to be extracted.

[0090] In step 7.2, the peripheral genes whose biological regulatory features are to be extracted are processed using the same method as in step 4, and the processed peripheral genes are input into the trained graph attention neural network to output the biological regulatory features of the peripheral genes; the biological regulatory features of the peripheral genes are concatenated with the topological regulatory features of the corresponding peripheral genes to obtain the regulatory features of the peripheral genes.

[0091] The effects of the present invention are further described below in conjunction with simulation experiments.

[0092] 1. Simulation conditions:

[0093] The hardware platform of the simulation experiment of the present invention has a processor of Intel(R) Core(TM) i7-8550U CPU@1.80GHz and a memory of 128GB.

[0094] The software platform of the simulation experiment of the present invention is: Windows 10 operating system, Pytorch 1.5.1, python 3.6.

[0095] The data set used in the simulation experiment of the present invention is: peripheral gene expression data of schizophrenia. The gene expression data is an expression matrix, which represents the expression values ​​of 15,906 genes on 59 samples. The 59 samples consist of 30 disease samples and 29 normal samples. It is the data disclosed by Huan Luo et al. in their published paper "Prognostic Implications of Novel Ten-Gene Signature in Uveal Melanoma[J].(Frontiers in Oncology,2020,10.)".

[0096] 2. Simulation experiment content and results analysis:

[0097] The simulation experiment of the present invention adopts seven methods, namely, the present invention and the existing feature extraction methods (ALL PEEP, Per PEEP, CorePEEP, Per Exp, Core Exp, DEG Exp), to extract the control features of the samples in the input data set respectively, classify the samples according to the control features, and plot the results of the seven methods into Figure 2 .

[0098] The prior art method for extracting PEEP features refers to: Menche et al. proposed a classification method for obtaining sample characteristics by establishing personalized perturbation spectra for samples in their paper "Integrating personalized gene expression profiles into predictive disease-associated gene pools[J].(Npj Systems Biology&Applications,2017,3(1):10.)", which is referred to as PEEP. The method for obtaining sample characteristics based on the personalized perturbation spectra of all genes in the interaction network is referred to as ALL PEEP, the method for obtaining sample characteristics based on the personalized perturbation spectra of peripheral genes is referred to as Per PEEP, and the method for obtaining sample characteristics based on the personalized perturbation spectra of core genes is referred to as Core PEEP.

[0099] The prior art Exp feature extraction method refers to the classification method of obtaining sample features by establishing expression profiles for samples proposed by Bin Tang et al. in their paper “Differential age-and disease-related effects on the expression of genesrelated to the arachidonic acid signaling pathway in schizophrenia. (Psychiatry Res. 2012 Apr 30; 196 (2-3): 201-6.)”, referred to as Exp. The method of obtaining sample features based on peripheral gene expression profiles is referred to as Per Exp, the method of obtaining sample features based on core gene expression profiles is referred to as Core Exp, and the method of obtaining sample features based on differentially expressed gene expression profiles is referred to as DEG Exp.

[0100] In order to verify the effect of the simulation experiment of the present invention, the classification results of each sample in the input data set are obtained according to the seven methods, and the true positive rate (TPR) and false positive rate (FPR) of each sample are obtained. The receiver operating characteristic curve is plotted with TPR as the ordinate and FPR as the abscissa to evaluate the classification results of the seven methods.

[0101] Combine the following Figure 2 The simulation diagram of the present invention is further described in the simulation experiment.

[0102] Figure 2 The horizontal axis represents the false positive rate, the vertical axis represents the true positive rate, the seven technical methods correspond to seven different types of curves, and each point on each curve represents the false positive rate and true positive rate obtained by each threshold of the method, among which PToC-GW represents the simulation experiment results of the present invention.

[0103] Combined with Figure 2 It can be seen that the peripheral gene features extracted by the present invention surpass the existing technology in sample classification, which proves that the present invention can learn the regulatory information of peripheral genes on core genes, and then apply it to disease subclassification to obtain higher accuracy.

[0104] The above description is only a specific example of the present invention and does not constitute any limitation to the present invention. Obviously, for professionals in this field, after understanding the content and principle of the present invention, it is possible to make various modifications and changes in form and details without departing from the principle and structure of the present invention, but these modifications and changes based on the idea of ​​the present invention are still within the scope of protection of the claims of the present invention.

Claims

1. A method for extracting peripheral gene regulation features based on graph representation learning and attention, characterized in that: The energy sampling and splicing of the diffusion of peripheral genes to core genes in the interaction network obtained by the graph representation learning algorithm are combined into topological regulatory features to construct a graph attention neural network. The steps of the feature extraction method include the following: Step 1: Determine the core genes and peripheral genes in the constructed network: Step 1.1, construct a human gene interaction network consisting of N nodes, each network node in the human gene interaction network corresponds to a gene, establish an edge between every two genes that interact, and model the network with an undirected graph G = (V, E), where N represents the total number of gene nodes in the human gene interaction network, V represents the total number of gene nodes in the undirected graph G, the value of V is equal to N, and E represents the total number of edges in the undirected graph G; Step 1.2, add up all the regulatory values ​​of each gene in the human gene interaction network to obtain the regulatory factor of the gene on the phenotype, select genes with regulatory factors greater than or equal to 4 on the phenotype in the human gene interaction network to form the core gene set, and select all gene nodes with regulatory coefficients of core genes on the phenotype less than or equal to 0.05 in the human gene interaction network to form the peripheral gene set; Step 2: Extract the topological regulatory features of each gene in the peripheral gene set on the core gene set: Step 2.1, calculate the energy of each peripheral gene node propagating to the core gene set according to the following formula: Among them, δ j represents the energy propagated from the jth gene node in the peripheral gene set to the core gene set, C represents the size of the core gene set, ∑ represents the summation operation, r represents the sequence number of the core gene in the core gene set, exp(·) represents the exponential operation with the natural constant e as the base, i represents the symbol of the imaginary part of the complex number, and λ o represents the oth eigenvalue in the eigenvalue diagonal matrix decomposed by the Laplace matrix, o = 1, 2, ..., O, O represents the number of eigenvalues ​​in the eigenvalue diagonal matrix decomposed by the Laplace matrix, s represents the diffusion factor, the value of s is equal to the average value of all eigenvalues ​​in the eigenvalue diagonal matrix decomposed by the Laplace matrix, U rl represents the energy value propagated from the lth gene node in the eigenvector matrix of Laplace matrix decomposition to the rth gene node in the core gene set, U lj represents the energy value propagated from the jth gene node in the eigenvector matrix decomposed by the Laplace matrix to the lth gene node in the peripheral gene set, l = 1, 2, ..., L, L represents the size of the peripheral gene set, and the value of l is equal to that of o; Step 2.2, sample the propagation of each gene in the peripheral gene set to each core gene in the core gene set to obtain the topological regulation characteristics of the gene node: c j =concat(Re(δ j ),Im(δ j )) Among them, γ j represents the topological regulatory features of the j-th gene node in the peripheral gene set, Re(·) represents the real part operation, Im(·) represents the imaginary part operation, and concat(·) represents the concatenation operation; Step 3, generate the initial selection vector for each gene node; Step 3.1, select q biological function entries, each of which contains coding gene nodes involved in the biological function; the total number of coding genes contained in the kth biological function entry is n k The k-th position of each coding gene is set to 1, and the k-th position of the coding gene not included in the k-th biological function entry is set to 0, thereby representing the situation that each coding gene in the node set V participates in q biological function entries as a p-dimensional initial biological feature vector, and the initial selection vector of each coding gene node is obtained by mean normalizing the initial biological feature of each coding gene node, where p = q, q represents the number of all biological function entries; Step 3.2, select the set of coding gene nodes X from the set of first-order neighbor gene nodes of each non-coding gene node in the undirected graph G; k The intersection of the coding gene set and the coding gene node set X is taken to obtain x genes; the number of randomly selected genes in the set of gene nodes is n k The intersection of the gene set and the coding gene node set X is x. r The overlapping node sets are randomly selected 1000 times to obtain 1000 overlapping node sets. The mean μ and variance σ of the selected 1000 overlapping node sets are calculated. Calculate the significance score z_score. If z_score is greater than 1.645, set the kth position of the non-coding gene node to 1. If z_score is less than 1.645, set the kth position of the non-coding gene node to 0. Thus, each non-coding gene is represented as a p-dimensional initial biological feature vector. The initial biological feature of each non-coding gene node is mean-normalized to obtain the initial selection vector of the non-coding gene node. Step 4: Generate training set: Step 4.1, set all peripheral genes in the peripheral gene set as positive samples, and set the label of the positive sample to 0; set the remaining genes in the gene node set excluding the peripheral gene set and the core gene set as negative samples, and set the label of the negative sample to 1; Step 4.2, randomly select 20% of positive samples and 20% of negative samples and their corresponding labels from the set of gene nodes after setting sample labels to form a training set, and remove samples other than the training set samples and their corresponding labels from the set of gene nodes after setting sample labels to form a test set; Step 4.3, calculate the feature average weight vector of the initial selection vector of all gene nodes in the training set according to the following formula: in, represents a feature average weight vector composed of all positive samples in the training set, F train represents the total number of samples in the training set, h j Represents the initial selection vector of the j-th gene node in the training set; Step 5: Build a graph attention neural network to extract the biological regulation features of peripheral genes: Step 5.1: Build a 5-layer graph attention neural network to extract the biological regulatory features of peripheral genes. Its structure is as follows: feature selection module, regulatory transformation group, splicing layer, KNN-graph attention module, output layer; set the number of neurons in the feature selection module to p s , p s The value of is equal to the dimension of the initial selection vector of the sample, the output layer is implemented by a K-nearest neighbor classifier, and the parameters of the K-nearest neighbor classifier are set to 64*2; Step 5.2: The control conversion group is composed of 20 control conversion layers with the same structure connected in parallel. The structure of each control conversion layer is: fully connected layer, first activation layer, conversion layer, second activation layer. The parameters of the fully connected layer are set to The first activation layer is implemented using the softmax function, and the number of neurons in the conversion layer is set to 32. The second activation layer is implemented using the sigmoid function, where the value of p is equal to the dimension of the initial selection vector of the sample; Step 6: Train the graph attention neural network: Step 6.1, input each sample in the training set into the graph attention neural network in turn, and output 20 regulatory features of the sample through 20 parallel regulatory conversion layers in the splicing layer of the graph attention neural network. After splicing the 20 regulatory features, the biological regulatory features of the sample are obtained, and the Euclidean distance between the biological regulatory features of every two samples in all samples is calculated. The dis samples with the smallest Euclidean distance to each sample are used as the neighbor set of the sample, thereby constructing an undirected graph G'=(V',E'); the undirected graph G' and the biological regulatory features of each sample are simultaneously input into the KNN attention pooling layer; wherein V' represents all sample nodes, and E' represents the set of edges established between each sample and each sample in its neighbor set; Step 6.2, using the stochastic gradient descent method, iteratively update the parameters until the loss function converges to obtain the trained graph attention neural network; Step 7: Extract the regulatory features of each peripheral gene: Step 7.1, using the same method as step 2, extracting the topological regulatory features of each gene in the peripheral gene set whose topological regulatory features are to be extracted; In step 7.2, the peripheral genes whose biological regulatory features are to be extracted are processed using the same method as in step 4, and the processed peripheral genes are input into the trained graph attention neural network to output the biological regulatory features of the peripheral genes; the biological regulatory features of the peripheral genes are concatenated with the topological regulatory features of the corresponding peripheral genes to obtain the regulatory features of the peripheral genes.

2. The peripheral gene regulation feature extraction method based on graph representation learning and attention according to claim 1 is characterized in that: The KNN-graph attention module described in step 5.1 is composed of a KNN attention pooling layer and a feedforward layer in series; the feedforward layer is implemented using the perceptron affine kernel method; the structure of the KNN attention pooling layer is: convolution layer, attention layer, pooling layer; the convolution kernel size of the convolution layer is set to 1×1, the attention layer is implemented using the activation function sigmoid, and the pooling layer uses average pooling.

3. The peripheral gene regulation feature extraction method based on graph representation learning and attention according to claim 1 is characterized in that: The loss function described in step 6.2 is as follows: Among them, L represents the loss function, y j represents the true label of the j-th gene node in the training set, log represents the logarithm operation with base 10, and p j represents the probability that the j-th gene node in the training set is predicted to be a peripheral gene, ||·|| represents the modulo operation, Represents the sum of the initial selection vectors of all positive samples in the training sample.

Images