A multi-phenotype-based genotype and phenotype association analysis method

Abstract

The application provides a kind of genotypes and phenotypic association analysis method based on multiple phenotypes, belongs to the technical field of bioinformatics, and establishes SNP-gene-phenotype three-layer network to explore the relationship between genotypes and phenotypes using statistical data, and the internal correlation between genomics and genomics, such as genes and genes, and phenotypes and phenotypes, is considered in the three-layer network construction process, which is more consistent with biological reality;Solve the problem that clinical data is difficult to obtain and SNP and phenotype correlation cannot be predicted;The biological pathway correlation between SNP-gene-phenotype different omics layers is analyzed.Through the analysis of various database data, it is found that the internal correlation of each database omics has quantitative value, while the correlation between omics generally has qualitative value, i.e. the correlation is 1 and the non-correlation is 0.Through the establishment of three-layer model, the key correlation genes can be predicted using the quantitative value of the internal omics and the qualitative value of the inter-omeric relationship, and the pathway relationship score between each layer is calculated to analyze the pathway relationship between SNP-gene-phenotype.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of bioinformatics, and particularly relates to a genotype and phenotype correlation analysis method based on multiple phenotypes. BACKGROUND

[0002] In the genotype and phenotype correlation analysis based on multi-omics data, multi-omics analysis is mainly performed on each omics data in clinical data. In clinical data, the horizontal and vertical coordinates of the matrix form data represent samples and omics data characteristics, such as SNP data, gene expression data, methylation data and the like in each clinical data. The data source is mainly obtained through clinical collection. Due to the protection of patient privacy and the management requirements of each institution on data, and the relatively large characteristics of each omics data, in the process of analyzing clinical data, the number of omics data types is often insufficient, and the sample size is relatively small. In view of the problem of small sample size, a clustering grouping method is proposed to reduce the characteristics and achieve the purpose of adapting to small sample application[1]. In view of the problem of insufficient number of omics, a gene hidden layer is introduced into the deep neural network, and the existing prior statistical information is used to assist the multi-omics fusion algorithm for analysis, and the effect of multi-omics data fusion is achieved with the aid of statistical data. Here, the statistical data refers to the association data between each omics obtained from the existing database. The data in each database is generally collected through biopsy experiments, algorithm verification, statistical induction and the like. This provides us with a new research idea. Since the statistical data is easy to obtain and has a fast update speed, it is considered in the present application whether the statistical data can be used to analyze the genotype and phenotype correlation. In the past multi-omics method research, the clinical data is used as the basis, and the statistical data is only used as verification data or only plays an auxiliary role. When the clinical data is analyzed, the phenotype data only contains the specific disease or trait to be studied, and the correlation between phenotypes is almost not used for analysis. It is thus conceived whether the statistical data can be used to establish a genotype and phenotype correlation network for analysis to achieve an effect similar to that of clinical multi-omics analysis.

[0003] The statistical data is used to analyze the genotype and phenotype correlation, and two problems are mainly considered: 1) whether the correlation between phenotypes can be used to analyze the genotype and phenotype relationship; and 2) how to use the statistical data to construct a correlation network to analyze the genotype and phenotype relationship pathway. SUMMARY

[0004] The application provides a genotype and phenotype correlation analysis method based on multiple phenotypes. The application uses statistical data to establish a SNP-gene-phenotype three-layer network to explore the genotype and phenotype pathway relationship, and the internal correlation between genes and genes, phenotypes and phenotypes and other omics is considered in the three-layer network construction process, which is more in line with the biological reality; the problem that clinical data is difficult to obtain and the SNP and phenotype correlation relationship cannot be predicted is solved; and the biological pathway correlation relationship between SNP-gene-phenotype different omics layers is analyzed. Through analysis of various database data, it is found that the internal correlation relationship of each database omics has a quantitative value, and the correlation relationship between omics generally has only a qualitative value, that is, the correlation is 1 and the non-correlation is 0, and through the establishment of the three-layer model, the key correlation genes can be predicted by using the internal quantitative value of the omics and the qualitative value of the correlation between the omics, and the pathway relationship score between each layer is calculated, so as to analyze the pathway relationship between SNP-gene-phenotype.

[0005] The technical scheme adopted by the application is:

[0006] A genotype and phenotype correlation analysis method based on multiple phenotypes comprises the following steps:

[0007] Step 1: establishing a SNP-gene-phenotype three-layer network model;

[0008] Step 2: dividing the SNP-gene-phenotype three-layer network model into two double-layer networks, respectively calculating the interlayer weight of the double-layer networks, combining the weights of the two double-layer networks through a hyperparameter β to obtain a final weight score, and taking the weight score as a pathway score;

[0009] Step 3: based on the calculation result of step 2, establishing a horizontal coordinate as a SNP and phenotype correlation relationship and a vertical coordinate as a gene pathway relationship matrix;

[0010] Step 4: combining the correlation relationship verified in each public database to verify the reliability of the pathway relationship.

[0011] In one specific embodiment of the application, the SNP-gene-phenotype three-layer network model in step 1 comprises:

[0012] The first layer is an SNP layer,

[0013] The second layer is a gene expression layer, the intralayer relationship of the gene layer is mapped by using a PPI network, and the interlayer correlation relationship between the SNP and the gene is embodied by eQTL data, and the upper two double-layer networks are constructed by using the data;

[0014] The third layer is a phenotype layer, and the phenotype layer includes intralayer relationship and gene correlation data related to the phenotype, and then the lower two double-layer networks are constructed by using the data.

[0015] In a specific embodiment of the present application, preferably, the upper two-layer double-layer network in step 2 only contains the intra-layer relationship of the gene layer and the inter-layer relationship between the SNP layer and the gene layer, and when solving the inter-layer weight of the upper two-layer double-layer network, the PBMDA algorithm is used to solve the inter-layer relationship value first, and the inter-layer weight of the upper two-layer double-layer network is solved by taking the value as the initial edge weight of the RWR algorithm.

[0016] In a more preferred specific embodiment of the present application, the specific calculation process of the inter-layer weight of the upper two-layer double-layer network includes:

[0017] (1) The PBMDA algorithm is used to solve the inter-layer relationship value between the SNP layer and the gene layer;

[0018] There are multiple path relationships between any SNP and any gene, assuming that P is a set of paths connecting m and d, m represents an SNP, and d represents a gene, then P = {p1, p2, …, p n}, wherein p i represents the weight of each edge of the path (0 < p i ≤ 1), P is converted into a relationship value P s , which can be represented as follows:

[0019]

[0020] The index F decay (P) is a decay factor, which is represented as follows:

[0021] F decay (P) = α × len(P) (2)

[0022] Wherein, the parameter α is a decay coefficient, which is set to 2.26, and len(P) represents the length of the path P, and in the algorithm implementation process, len(P) is set to be less than or equal to 3; P s represents the influence of a path on the SNP-gene relationship, and for all paths, the SNP-gene relationship score(m, d) can be represented as:

[0023]

[0024] Wherein, P si represents the P s value of the i-th path in all paths, and the greater the value of score(m, d) is, the greater the correlation between m and d is;

[0025] (2) The value of score(m, d) is taken as the initial edge weight of the RWR algorithm to solve the inter-layer weight of the upper two-layer double-layer network;

[0026] The RWR algorithm can be represented as:

[0027] P t+1 = (1 - γ)W1P t + γP 0 (4)

[0028] where γ ∈ (0, 1) is the restart probability, representing the probability of returning to the seed node at each step of the random walk, P 0 is the initial probability vector, i.e., score(m, d); P t represents the probability vector of the random walk at time t; then the transition matrix of the two-layer network constructed by the SNP layer and the gene layer is where W S and W G represent the intra-layer transition matrices of the SNP layer and the gene layer, respectively, W GS and W SG are the inter-layer transition matrices of the SNP layer and the gene layer, respectively, and W GS is the transpose matrix of W SG . In the process of solving W GS , the inter-layer network association relationship is divided into two cases: associated and not associated, which belongs to a bipartite network; for example, in the W GS matrix, if the SNP and the gene have an association relationship, the value is 1, otherwise it is 0. The W GS matrix is normalized, and the formula is as follows:

[0029]

[0030] Similarly,

[0031]

[0032] where w′ ij represents the element in the i-th row and the j-th column of the adjacency matrix of the SNP layer and the gene layer.

[0033] In one specific embodiment of the present application, the lower two-layer network in step 2 is calculated by the CIPHER algorithm + RWR algorithm to obtain the inter-layer weight of the two-layer network, and the specific calculation process includes:

[0034] (1) The CIPHER algorithm is based on the assumption that genes with similar functions lead to similar phenotypes, and a regression model is used in the algorithm to incorporate this assumption. The similarity vector between phenotypes can be defined as:

[0035]

[0036] where S pp′ is the similarity value between phenotypes p and p'; L gg′is the topological distance between genes g and g' in the PPI network; G(p) is the set of all genes associated with phenotype p; G(p') is the set of genes associated with phenotype p'; C p is a constant, and is understood that if a gene in G(p') has no interaction with a gene in G(p) in the gene layer, then S pp′ is C p ; β pg is the coefficient of the entire regression model, and β pg represents the contribution of gene g in the similarity calculation of phenotype p and p'; in order to quantify the association between phenotype and gene, the similarity of gene g and phenotype p' can be defined as the sum of the correlation of all genes associated with phenotype p' and gene g, i.e.,

[0037]

[0038] where Φ gp′ represents the relationship value between gene g and phenotype p', and substituting the formula into the formula, it can be expressed as:

[0039] S pp′ = C p +∑ g∈G(p) β pg Φ gp′ (9)

[0040] If there are n phenotypes in the phenotype layer that are associated with gene g and phenotype p, the relationship vector between phenotype p and all other phenotypes in the phenotype layer can be obtained from the phenotype layer Similarly, the relationship vector between gene g and all phenotypes can be calculated using formula (8) Thus, formula (9) is extended to the following form:

[0041] S p = C p +∑ g∈G(p) β pg Φ g (10)

[0042] In the linear regression model of the above formula, the relationship value between phenotype p and gene g can be calculated using the Pearson linear correlation coefficient, which is defined as follows:

[0043]

[0044] where cov(S p ,Φ g ) represents the covariance between S p and Φ g , and σ(S p ) represents the standard deviation of vector S p , and σ(Φ g) represents the standard deviation of the vector Φ g The correlation coefficient CS pg reflects the correlation between the gene g in the gene layer and the phenotype p in the phenotype layer.

[0045] (2) The relationship value calculated above is used as the initial value of the edge weight of the RWR algorithm, and the RWR algorithm can be expressed as:

[0046] P t+1 = (1-γ) W2P t + γP 0 (12)

[0047] wherein γ ∈ (0, 1) is a restart probability, representing the probability of returning to the seed node at each step of random walk, P 0 is an initial probability vector, i.e., the correlation coefficient CS pg ; P t represents the probability vector of random walk at time t; and the transition matrix of the double-layer network constructed by the gene layer and the phenotype layer is wherein W G and W P represent the intra-layer transition matrices of the gene layer and the phenotype layer, W PG and W GP are the inter-layer transition matrices of the gene layer and the phenotype layer, W PG is the transpose matrix of W GP , and in the process of solving W PG , the inter-layer network correlation is divided into two cases: associated and not associated, which belongs to a binary network; for example, in the W PG matrix, if the gene and the phenotype have an association relationship, the value is 1, otherwise, the value is 0, and the W PG matrix is normalized, and the formula is as follows:

[0048]

[0049] Similarly,

[0050]

[0051] wherein w i ' j ' represents the i-th row and j-th column element in the adjacency matrix of the gene layer and the phenotype layer.

[0052] In a specific embodiment of the present application, the restart probability γ is preferably 0.7, and the hyperparameter β is 0.7.

[0053] Advantages of the present application:

[0054] 1. The present application establishes a SNP-gene-phenotype three-layer network model by using statistical data, considers the internal correlation of omics, and considers the correlation between omics, so as to more accurately reflect the biological mechanism.

[0055] 2. The present application divides the SNP-gene-phenotype three-layer network model into two double-layer networks, then calculates the interlayer weight of each double-layer network, combines the weights of the two double-layer networks by using a hyperparameter β, obtains the final weight score, judges the tightness of the pathway relationship according to the score, and obtains the pathway relationship between SNP-gene-phenotype, which is more accurate.

[0056] 3. The present application adopts the PCRWR algorithm, and after increasing the initial weight, the AUC value of the random walk algorithm is increased by about 7 percentage points, and the PCRWR algorithm performs best in terms of algorithm performance and stability. BRIEF DESCRIPTION OF DRAWINGS

[0057] In order to more clearly illustrate the technical solutions in the embodiments of the present application or the prior art, the drawings needed to be used in the embodiments or prior art description will be briefly introduced below. Obviously, the drawings in the following description are only some embodiments of the present application, and other drawings can be obtained by those skilled in the art without creative labor.

[0058] Figure 1 The flowchart of the present application;

[0059] Figure 2 The phenotype and miRNA pathway relationship diagram;

[0060] Figure 3 The gene-phenotype double-layer network regression flowchart;

[0061] Figure 4 The schematic diagram for comparing the performance of each algorithm by using the ROC curve;

[0062] Figure 5 The statistical diagram of the influence of the number of diseases on the performance of the algorithm. DETAILED DESCRIPTION

[0063] In order to make the purpose, technical scheme and advantages of the embodiments of the present application more clear, the technical scheme in the embodiments of the present application will be described clearly and completely below in combination with the drawings in the embodiments of the present application. Obviously, the described embodiments are part of the embodiments of the present application, not all the embodiments. The components of the embodiments of the present application described and shown in the drawings here can be arranged and designed in various different configurations.

[0064] The following detailed description of embodiments of the application provided in the accompanying drawings is not intended to limit the scope of the application as claimed, but merely represents selected embodiments of the application. Based upon the embodiments in the present application, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of the present application.

[0065] Referring to Figure 1 The present application provides a multi-phenotype-based genotype and phenotype correlation analysis method, comprising the following steps:

[0066] Step 1: Establish a SNP-gene-phenotype three-layer network model;

[0067] Step 2: Divide the SNP-gene-phenotype three-layer network model into two double-layer networks, calculate the inter-layer weights of the double-layer networks respectively, combine the weights of the two double-layer networks through a hyperparameter β to obtain a final weight score, and use the weight score as a pathway score;

[0068] Step 3: Based on the calculation results of Step 2, establish a SNP-phenotype correlation relationship as the horizontal coordinate and a gene pathway relationship matrix as the vertical coordinate;

[0069] Step 4: Verify the reliability of the above pathway relationship in combination with the correlation relationships confirmed in various public databases.

[0070] The SNP-gene-phenotype three-layer network model is as follows: Figure 1As shown in the middle part of a, the first layer is SNP layer. The association between SNPs is generally referred to as "epistatic association", and epistasis is an important genetic component of phenotypic variation and a key mechanism to explain the missing heritability. Identifying epistatic interactions in genetic association studies can help better understand the genetic structure of complex traits and diseases. Many statistical methods have been used to model and identify epistatic interactions between genetic variations. There are methods that use SNP position information to integrate summary data from GWAS and eQTL studies to predict complex trait gene targets. And use the relationship between SNPs and gene positions to model the causal relationship chain of genetic variations affecting chromatin and thus affecting gene expression, to better identify causal variations and chromatin markers that affect gene expression. However, this model does not consider epistatic association, because epistatic association is mainly established through SNP position relationship or linkage disequilibrium analysis, for example, multiple SNP sites are in the same gene segment, then it is considered that there is an association between these SNPs, and there are many epistatic relationships that are reflected through the association between genes. Such relationships have been reflected in eQTL data and gene association relationships (such as PPI networks). In addition, the SNP data that people are interested in is in the tens of thousands or even millions, if the SNP association relationship is added in this network, the SNP relationship matrix formed is not only too large, but also too sparse relative to the characteristic quantity of SNP data. For the above reasons, the SNP layer association is not added in the model building process.

[0071] The second layer is the gene expression layer, referred to as the gene layer. The intra-layer relationship of the gene layer is mapped using the PPI network, and the inter-layer association between SNPs and genes is reflected through eQTL data. Through such data, a two-layer network (SNP layer and gene layer) can be constructed, as shown in the middle part of b) of the upper half. Figure 1 Such a network can be analyzed using a two-layer network algorithm, thereby generating a relationship matrix between SNPs and genes, as shown in the middle part of c) of the upper half. Figure 1

[0072] The third layer is the phenotype layer. The Human Phenotype Ontology (HPO) provides relationships between phenotypes, and also contains gene association data related to the phenotype, which can also be downloaded from other public databases. With the intra-layer and inter-layer association relationship data, a two-layer network similar to the construction process of the upper two layers can also be used to solve the association matrix between genes and phenotypes using a two-layer network algorithm.

[0073] ​For the phenotype network, we can focus on a specific phenotype in the network, establish a phenotype-related network for the concerned phenotype to assist in analysis, and obtain the relationship path of the phenotype. We can also cluster the entire phenotype and analyze the phenotype network formed by a cluster. In the implementation process of the algorithm in this paper, the entire phenotype network is first clustered. Phenotypes in the same cluster are considered more likely to be related, and then the phenotypes in this cluster are analyzed in detail. The specific method is to use the K-means++ clustering algorithm (see K-means algorithm for details) to find closely related phenotypes.

[0074] K-means algorithm

[0075] Due to the large number of disease types, the number of associated relationships generated in the three-layer network increases exponentially, and processing such a large amount of data will affect the efficiency and accuracy of the algorithm. Therefore, the algorithm in this paper simplifies the network by clustering and analyzes the disease traits with stronger correlation together. Specifically, in the implementation process, the data is first filtered. The entire phenotype network is clustered using a clustering algorithm to generate different disease sets or trait clusters, and the concerned disease set or trait cluster is selected. Combined with the three-layer network, the relevant genes and SNPs are filtered layer by layer through the concerned disease set or trait cluster for pathway research. This paper uses the K-means clustering algorithm to cluster the nodes of the disease association network.

[0076] The K-means clustering algorithm can be described as follows: Let X = {x1,..., xn} be a sample set containing n samples, which is divided into k categories, denoted as C = {c1,..., ck}. The cluster centers of these k categories are set as A = {a1,..., ak}. n}。 k k

[0077] First, randomly select k from the n samples as the initial cluster center points, and calculate the distance d of each point in the data to the k cluster center points, denoted as:

[0078]

[0079] Second, classify each sample point. The classification principle is to attribute the sample point to the class in which the nearest cluster center point is located.

[0080] Third, reselect the cluster center points according to the k classes divided. The new cluster center point selection formula is:

[0081]

[0082] where N j represents the number of samples in class c j . ​​

[0083] Fourthly, check if the termination condition is reached, if yes, stop the operation, if not, repeat the second and third steps until the termination condition is reached.

[0084] The termination condition can be any of the following:

[0085] (1) the objects in the cluster remain unchanged or within the allowed range;

[0086] (2) the number of iterations reaches the preset maximum number;

[0087] (3) the error sum of squares is locally minimized. The error sum of squares formula is:

[0088]

[0089] The data processed by the K-means clustering algorithm is divided into several categories, and only the data of a certain category needs to be selected for analysis, so as to achieve the purpose of data simplification and stronger network correlation. Since the disease data obtained in this paper is mainly the correlation data between diseases, when using the clustering algorithm, this paper adopts the enhanced K-means clustering algorithm-K-means++ clustering algorithm. Compared with the K-means clustering algorithm, the only difference between the K-means++ clustering algorithm is that the distance relationship between sample points is mainly relied on when selecting the clustering center, and its principle is the same as that of the K-means clustering algorithm.

[0090] Double-layer network algorithm

[0091] The PBMDA algorithm, CIPHER algorithm and RWR algorithm are used to analyze the double-layer network.

[0092] 1. PBMDA algorithm

[0093] The path-based double-layer network algorithm (PBMDA algorithm) is an algorithm for predicting the relationship between miRNA and disease proposed by You Z H et al

[17] .

[0094] There are multiple path relationships between any miRNA m and any disease d in the miRNA-disease double-layer network, and these paths are composed of intra-layer networks such as disease association network and miRNA association network and inter-layer network miRNA-disease. The intra-layer network has both association relationship and association weight, while the inter-layer network only has association relationship (association is 1 and non-association is 0). The model schematic diagram is shown in Figure 2 .

[0095] Since there can be multiple paths for any m and d, we assume P is a set of paths connecting m and d, P = {p1, p2,..., p n}. Where p i represents the weight of each edge that makes up the path (0 < p i ≤ 1), then P is converted into a relationship value P s which can be expressed as follows:

[0096]

[0097] The index F decay (P) is a decay factor, which is expressed as follows:

[0098] F decay (P) = a x len(P) (5)

[0099] Where the parameter a is the decay coefficient, set to 2.26. len(P) represents the length of the path P. As can be seen, the longer the path, the greater the decay factor, and the more distant the relationship between m and d. When the path length is greater than 3, since the value obtained is too small, it is not considered to be of much significance, so in the subsequent algorithm implementation process, len(P) is set to be less than or equal to 3. P s represents the influence of a path on the miRNA-disease relationship. For all paths, the miRNA-disease relationship score(m, d) can be expressed as:

[0100]

[0101] Where P si i represents the P s value of the i-th path in all paths. The greater the value of score(m, d), the greater the correlation between m and d. Take Figure 2 for example, there are multiple paths between miRNA m1 and disease d2. For example: m1→d2, m1→d1→d2, m1→m3→d2, m1→m4→m3→d2, and the remaining paths are greater than 3 and are not considered, such as: m1→m4→d4→m3→d2. Therefore, according to formula (6), the score value is calculated as follows:

[0102] score(m1, d2) = 1.0 2.26*1 + (1.0 * 0.7) 2.26*2 + (0.6 * 1.0) 2.26*2 + (0.6 * 0.4 * 1) 2.26*3

[0103] ≈1+0.1995+0.0994+0.00006

[0104] = 1.29896 (7)

[0105] 2. CIPHER algorithm

[0106] CIPHER algorithm is also based on the hypothesis that genes with similar functions lead to similar phenotypes, and incorporates this hypothesis in the algorithm using a regression model. Finally, the score from the Pearson linear correlation coefficient is used to assess the likelihood of a gene being associated with a specific phenotype.

[0107] To build this regression model, the quantified similarity between phenotypes and the quantified interaction between genes and proteins are needed in the layer, as well as the complete list of known disease gene-phenotype associations between layers. The detailed data sources are shown in the data sources and preprocessing section.

[0108] The idea of the algorithm is: first, combine the human phenotype network, protein network and gene-phenotype association into a complete network, and then solve for a specific phenotype p and gene g. The specific process is divided into three steps: 1) find the relationship vector of phenotype p with other phenotypes in the phenotype layer; 2) calculate the relationship vector of gene g with other phenotypes through the double-layer network; 3) calculate the linear correlation of the two vectors as the association score between phenotype p and gene g, as shown in the following equation: Figure 3

[0109] The similarity vector between phenotypes can be defined as:

[0110]

[0111] where S pp′ is the similarity value between phenotype p and p'; L gg′ is the topological distance between genes g and g' in the PPI network, and the shortest path distance is used in the experiment; G(p) is the set of all genes related to phenotype p; G(p') is the set of genes related to phenotype p'; C p is a constant, which can be understood as if the genes in G(p') and the genes in G(p) have no interaction in the gene layer, then the S pp ' relationship value is C p ; β pg is the coefficient of the entire regression model, and β pg represents the contribution of gene g in the similarity calculation of phenotype p and p'. As can be seen from the formula, in this regression model, the similarity between phenotypes can be reflected by the close relationship between the genes related to the phenotypes in the PPI network. The similarity between phenotypes is converted into the interaction between the genes related to the phenotypes in the PPI network. Similarly, to quantify the association between phenotypes and genes, the similarity of gene g and phenotype p' can be defined as the sum of the relevance of all genes related to phenotype p' and gene g, i.e.:​

[0112]

[0113] where Φ gp′ represents the relationship value between the gene g and the phenotype p', and formula (9) is substituted into formula (8), which can be expressed as:

[0114] S pp′ = C p +∑ g∈G(p) β pg Φ gp′ (10)

[0115] If there are n phenotypes related to the gene g and the phenotype p in the phenotype layer, the relationship vector between the phenotype and all other phenotypes can be obtained from the phenotype layer Similarly, the relationship vector between the gene g and all phenotypes is calculated by formula (9) Thus, formula (10) is extended to the following form:

[0116]

[0117] In the linear regression model of the above formula, the relationship value between the phenotype p and the gene g can be calculated by using the Pearson linear correlation coefficient, which is defined as follows:

[0118]

[0119] where cov(S p ,Φ g ) represents the covariance between S p and Φ g , σ(S p ) represents the standard deviation of the vector S p , and σ(Φ g ) represents the standard deviation of the vector Φ g . The correlation coefficient CS pg reflects the degree of association between the gene g in the gene layer and the phenotype p in the phenotype layer. The greater the value of CS pg , the closer the relationship between the gene g and the phenotype p.

[0120] 3. RWR algorithm

[0121] The random walk algorithm with restart (RWR algorithm) is improved on the basis of the random walk algorithm. The random walk algorithm with restart can be expressed by the formula:

[0122] P t+1 = (1-γ) WP t + γP 0 (13)

[0123] where γ ∈ (0, 1) is the restart probability, representing the probability of returning to the seed node at each step of the random walk, P0is the initial probability vector, i.e., the correlation coefficient CS pg ; P t represents the probability vector of the random walk at time t; the transition matrix of the two-layer network is where W U and W V represent the intra-layer transition matrices of the first and second layers, respectively, W UV and W VU are the inter-layer transition matrices, W UV is the transpose matrix of W VU . In the process of solving W VU , the inter-layer network association relationship is divided into two cases: associated and not associated. If there is an association relationship, the value is 1, otherwise it is 0. The W VU matrix is normalized, and the formula is as follows:

[0124]

[0125] Similarly,

[0126]

[0127] where w ij represents the element in the i-th row and j-th column of the adjacency matrix.

[0128] As can be seen from the formula, in the use of RWR algorithm, the current node in the inter-layer network is assigned to all adjacent nodes with the same weight, i.e., the transition probability of each adjacent node is the same, which does not conform to the actual situation. In actual situations, genes in the gene layer have different degrees of association with different related SNPs, and similarly, phenotypes in the phenotype layer should also have different degrees of association with different related genes. Therefore, the transition probability of different nodes in the random walk should be considered differently.

[0129] Algorithm implementation

[0130] In the process of constructing a three-layer network, the upper two layers (SNP layer and gene expression layer) only contain the intra-layer relationship of the gene layer and the inter-layer relationship between the SNP layer and the gene layer, while the CIPHER algorithm requires the SNP intra-layer relationship to complete the inter-layer weight solving of the two-layer network. Therefore, the upper two layers only select the PBMDA algorithm to solve their inter-layer relationship values, and use them as the initial edge weight values of the random walk for association analysis. The lower two layers can have multiple choices, which can be PBMDA algorithm + RWR algorithm or CIPHER algorithm + RWR algorithm.

[0131] Considering that the influence of genes on SNPs and phenotypes may be inconsistent, a hyperparameter β should be added when establishing associations. For example, after solving the interlayer weights of two two-layer networks, if the weight of the relationship between SNP1 and gene1 is m, and the weight of the relationship between gene1 and phenotype1 is n, then the weight of the relationship between SNP1 and phenotype1 is considered to be m + βn, and this is used as the pathway score (PS). The hyperparameter β can be obtained through enumeration testing. This establishes a pathway relationship matrix with the SNP-phenotype association on the x-axis and the gene-phenotype association on the y-axis. The SNP-phenotype association can be verified by searching public databases such as PhenoScanner and DisGeNET (an association in the database is defined as 1, and no association is defined as 0). Figure 1 As shown in section d). Subsequent calculations have transformed the biological mechanism research into a mathematical problem through modeling, requiring only the following... Figure 1 The matrix in part d) can be classified binary using logistic regression. To ensure a balance between positive and negative samples, we use the same number of positive and negative samples, and employ K-fold cross-validation to avoid overfitting. During the solution process, key genes can be identified using PS values, and pathway relationships between SNPs, genes, and phenotypes can be established, facilitating a deeper understanding of disease mechanisms.

[0132] Data Sources and Preprocessing

[0133] As described in the algorithm introduction above, to verify the effectiveness of the algorithm, the following types of relational data are required: SNP-gene association data, gene-gene association data, gene-phenotype association data, phenotype-phenotype association data, and SNP-gene-phenotype association data. Here, phenotype is primarily reflected through disease phenotype data.

[0134] SNP-gene association data comes from the BioMart website, official website www.biomart.org. This portal provides access to numerous database projects supported by 30 scientific organizations, including over 800 different biological datasets covering genomics, proteomics, model biology, cancer data, ontology information, etc. The website also provides better support and scalability for data analysis and visualization tools. Using the biomaRt package in R to connect the hsapiens_snp and hsapiens_gene_ensembl sub-databases for data matching and naming transformation, only SNP data with MAF greater than 0.05 were used. The final processed SNP-gene association data consisted of 383,668 records.

[0135] The gene-gene association data is mapped from PPI network data. The PPI network data is from HPRD database, which contains not only the association information between proteins, but also the specific weight of the association. There are 402815 records in total.

[0136] The gene-phenotype association data is from DisGeNET database. DisGeNET is one of the largest databases containing human diseases and related genes. The data version used in this paper is DisGeNET(v7.0). It contains 21,671 genes and 30,170 diseases or abnormal phenotypes, and 1,134,942 associations between them. The obtained data contains multiple columns, such as geneSymbol, diseaseName, etc.

[0137] The phenotype-phenotype association data is from the literature

[27] , which contains 383 diseases.

[0138] The SNP-gene-phenotype data is from PheGenI database

[28] , which is used for the final data verification. The pathway relationships obtained by this algorithm are compared with the pathway relationships in the database to illustrate the accuracy of the algorithm. In the implementation process of the algorithm, the validation set can be used to screen the features of the gene layer network and the phenotype layer network, reducing the number of nodes in each layer network.

[0139] Analysis of experimental results

[0140] Analysis of prediction results

[0141] To analyze the genotype-phenotype pathway relationship using statistical data, if only the double-layer network constructed by genotype and phenotype is analyzed, the pathway information cannot be obtained, so in the comparative experiment of the results, all models are compared based on the SNP-gene-phenotype three-layer network architecture.

[0142] In the comparative experiment, first, the PBMDA algorithm+RWR algorithm is used in the upper two layers, and the PBMDA algorithm+RWR algorithm is used in the lower two layers (this method is marked as PPRWR), and the PBMDA algorithm+RWR algorithm is used in the upper two layers, and the CIPHER algorithm+RWR algorithm is used in the lower two layers (this method is marked as PCRWR). In order to verify the effectiveness of the method, the random walk algorithm with restart (this method is marked as RWR) is also added, the PBMDA algorithm is used in the upper and lower two layers, and the PBMDA+PBMDA method is used without using the random walk algorithm, and the PBMDA algorithm is used in the upper two layers, and the CIPHER algorithm is used in the lower two layers, and the PBMDA+CIPHER method is used without using the random walk algorithm.

[0143] Such binary classification problem mostly uses cross-validation method to evaluate. In this paper, K-fold cross-validation method is mainly used, through the solution of true positive rate (TP) and false positive rate (FP) indicators, the ROC curve is drawn and the AUC value under the curve is calculated to evaluate the algorithm. There are multiple parameters in the algorithm that need to be determined, such as the number of disease categories in the disease network, the weight distribution β between the two double-layer networks and the restart probability γ. For specific parameter value analysis, see the parameter analysis section, and in Figure 4 In the comparison of the performance of each algorithm using the ROC curve, we use the optimal value of each parameter for analysis, in which the number of disease categories is 112, the weight distribution β is 0.7, and the restart probability γ is 0.7.

[0144] From the analysis of the ROC curve results in Figure 0-4, it can be seen that the AUC values of PCRWR, PPRWR and RWR algorithms are higher than those of PBMDA+PBMDA algorithm and PBMDA+CIPHER algorithm, that is, the performance of algorithms containing restart random walk is generally higher than that of algorithms not containing the algorithm, which shows that the restart random walk algorithm has better effect in the analysis of the three-layer network. In addition, whether the initial weight is calculated by PBMDA algorithm or CIPHER algorithm, the results of PCRWR and PPRWR two weight random walk algorithms are better than that of RWR algorithm which does not contain weight, but overall, the weight calculated by CIPHER algorithm makes the effect of the whole algorithm better. In order to prove the superiority and robustness of PCRWR algorithm, we solve the AUC value of ROC curve for each algorithm 10 times, and the results are shown in Table 1. The AUC value of PBMDA+PBMDA algorithm is 0.675±0.041, the AUC value of PBMDA+CIPHER algorithm is 0.731±0.036, the AUC value of RWR algorithm is 0.802±0.049, the AUC value of PPRWR algorithm is 0.848±0.034, and the AUC value of PCRWR algorithm is 0.874±0.018. It can be found that PCRWR algorithm performs best in terms of algorithm performance and stability.

[0145] Table 1 AUC value of ROC curve of each algorithm 10 times statistical results

[0146]

[0147] Parameter analysis

[0148] There are three hyperparameters in the whole algorithm that need to be determined: the number of diseases, the weight distribution β between the two double-layer networks and the restart probability γ.

[0149] Firstly, the number of disease categories directly affects the complexity of the disease-related network, and then affects the entire three-layer network structure. Therefore, we control the number of diseases by K-means++ clustering to verify the influence of the number of diseases on the algorithm. In the original data of the disease, there are 383 diseases, and after preprocessing the data, the number of diseases is finally 201. In K-means++ clustering, we set the number of clusters to 2 and 3, respectively, and generate multiple disease sets with disease numbers of 38, 69, 112, 168, etc. The performance of the algorithm is tested with these disease numbers, and the results are shown in Figure 5 As the number of diseases changes, the performance of the algorithm also changes. When the number of diseases gradually increases, the performance of the algorithm improves, but when the number of diseases is greater than a certain number, its performance tends to decline, indicating that the number of diseases has a certain influence on the algorithm, and the size of the influence is related to the specific number of diseases. When the number of diseases is small, it will cause poor performance due to its inability to fully reflect the pathway relationship, but when the number of diseases reaches a certain level, even too much, it will produce some diseases with not enough close association, thereby affecting the performance of the algorithm.

[0150] The weight distribution β between the two double-layer networks is to control the proportion of the influence of the two double-layer networks on the results. In theory, its value range should be between (0, 1), but in the experiment, in order to facilitate calculation, we lock its value in [0.1, 0.9], and try to enumerate with an increment of 0.1. Finally, it is found that when β = 0.7, the result is the best, that is, the influence of the gene and disease double-layer network on the result is greater.

[0151] Like the weight distribution parameter β, the value range of the restart probability γ is also set to [0.1, 0.9]. Through experiments, it is found that when the incremental interval is set to 0.2, it can be better distinguished and compared, as shown in Table 2, as γ increases, the performance of the algorithm also improves, but when γ is greater than 0.7, the performance tends to decline, so in this experiment, the value of γ is set to 0.7.

[0152] Table 2 Influence of restart probability on PCRWR algorithm

[0153]

[0154] The present application mainly solves the problem of analyzing multi-omics pathway relationship using statistical data when clinical data is difficult to obtain. The main innovation lies in using statistical data in public databases that are easy to obtain to achieve the effect of analyzing clinical multi-omics data. At the same time, the phenotype association relationship that is difficult to integrate with clinical data is integrated in the analysis process.

[0155] In order to make the path relationship among SNP-gene-phenotype clearer, the intra-layer correlation and its quantitative relationship weight and the inter-layer correlation and its qualitative relationship weight are required in terms of data requirements. With the data as support, a three-layer heterogeneous network can be established. In order to facilitate analysis, the three-layer network is split into two double-layer networks in the algorithm, the inter-layer weight is predicted respectively, and the weights of the two double-layer networks are combined by increasing the super parameter β to obtain the final weight score. Then, the accuracy of the method is illustrated by using the K-fold cross-validation method in combination with the correlation relationship verified in each public database.

[0156] In the solution of the inter-layer correlation weight in the double-layer network, the main idea of the PBMDA algorithm, the CIPHER algorithm and the RWR algorithm is to map the intra-layer correlation relationship to the inter-layer relationship. For the methods with the same solution idea, the methods are combined for use, that is, in the original equal probability random walk algorithm, the PBMDA algorithm and the CIPHER algorithm are used to generate the biased weight of the random walk, and the PCRWR method with higher accuracy is obtained. The results show that after increasing the initial weight, the AUC value of the random walk algorithm is increased by about 7 percentage points.

[0157] The above description is only used to illustrate the technical solutions of the present application but not limit the present application. Other modifications or equivalent replacements to the technical solutions of the present application made by those skilled in the art should be covered in the scope of the claims of the present application as long as they do not depart from the spirit and scope of the present application.

Claims

1. A method for genotype-phenotype association analysis based on multiple phenotypes, characterized in that, The method comprises the following steps: Step 1: establishing a SNP-gene-phenotype three-layer network model; Step 2: The SNP-gene-phenotype three-layer network model is divided into two double-layer networks, and the inter-layer weights of the double-layer networks are calculated respectively The weights of the two double-layer networks are combined to obtain the final weight score, which is used as the pathway score; Step 3: based on the calculation results of step 2, a coordinate for SNP and phenotype association relationship and a vertical coordinate for gene pathway relationship matrix are established; Step 4: the reliability of the above pathway relationship is verified by combining the verified association relationship in each public database; The SNP-gene-phenotype three-layer network model in step 1 comprises: The first layer is the SNP layer, The second layer is the gene expression layer, the intra-layer relationship of the gene layer is mapped using the PPI network, and the inter-layer association relationship between the SNP and the gene is embodied through the eQTL data, and the above two layers of double-layer network are constructed through the data; The third layer is the phenotype layer, which contains the intra-layer relationship of the phenotype layer and the gene association data related to the phenotype, and then the lower two layers of double-layer network are constructed through the data; In step 2, the upper two layers of double-layer network only contain the intra-layer relationship of the gene layer and the inter-layer relationship between the SNP layer and the gene layer, when solving the inter-layer weight of the upper two layers of double-layer network, PBMDA algorithm is used to solve the inter-layer relationship value, and the initial value of the edge weight of RWR algorithm is used to complete the inter-layer weight solving of the upper two layers of double-layer network; The lower two layers of double-layer network are calculated by CIPHER algorithm+RWR algorithm to obtain the inter-layer weight of the double-layer network, and the inter-layer relationship value is solved by CIPHER algorithm, and the initial value of the edge weight of RWR algorithm is used to complete the inter-layer weight solving of the lower two layers of double-layer network.

2. The method of claim 1, wherein the method is a multi-phenotype based genotypic and phenotypic association analysis method. The specific calculation process comprises: (1) the inter-layer relationship value of the SNP layer and the gene layer is solved by PBMDA algorithm; There are multiple pathway relationships between any SNP and any gene, and the assumption Let Pm,d be the set of paths connecting m and d, where m represents a SNP and d represents a gene, then where represents the weight of each edge that makes up the path , then is converted to a relationship value can be represented as follows: （1） Exponent is a decay factor, which is expressed as follows: （2） Wherein, the parameter a is the attenuation coefficient, set to 2.26, Indicates The length of the path, in the process of algorithm implementation, Set less than or equal to 3; Indicates the influence of a path on SNP-gene relationship, for all paths, SNP-gene relationship Can be expressed as: （3） wherein, represents the value of the mth path of the dth path, The greater the value of m, the greater the correlation between m and d. (2) take the value of as the initial value of the edge weight of the RWR algorithm to complete the weight solving of the upper two layers of the double-layer network. RWR algorithm can be expressed as: （4） where, is the restart probability, which represents the probability of returning to the seed node at each step of the random walk, is the initial probability vector, i.e. ; represents the probability vector of the random walk at time t; the transition matrix of the two-layer network constructed by the SNP layer and the gene layer is , where and are the intra-layer transition matrices of the SNP layer and the gene layer, respectively, and are the inter-layer transition matrices of the SNP layer and the gene layer, is the transpose matrix of , and in the process of solving , the inter-layer network correlation relationship is divided into two cases: associated and not associated, which belongs to a binary network; for example, in the matrix, if the SNP and the gene have an association relationship, the value is 1, otherwise it is 0. The matrix is normalized, and the formula is as follows: （5） Similarly, （6） wherein, represents the element in the i-th row and j-th column of the SNP layer and gene layer adjacency matrix.

3. The method of claim 1, wherein the method is a multi-phenotype based genotypic and phenotypic association analysis method. In step 2, the lower two layers of double-layer network are calculated by CIPHER algorithm+RWR algorithm to obtain the inter-layer weight of the double-layer network, and the specific calculation process comprises: (1) CIPHER algorithm is based on the assumption that genes with similar functions lead to similar phenotypes, and the assumption is integrated into the algorithm by using a regression model, and the similarity vector between phenotypes can be defined as: （7） wherein, is a phenotype and a similarity value between them; is the topological distance between genes and in the PPI network; set is the set of all genes associated with phenotype ; set is the set of genes associated with phenotype ; is a constant, understood to be if the genes in have no interaction in the gene layer, then the relationship value is ; is the coefficient of the entire regression model, representing the contribution of gene to phenotype and in the similarity calculation; in order to quantify the association between the phenotype and the gene, the similarity of the gene to the phenotype can be defined as the sum of the relevance of all genes associated with phenotype to gene , that is: （8） wherein representative genes and phenotypes between the relationship values, formula (8) is substituted into formula (7), which can be expressed as: （9） If the phenotype layer has n phenotypes, the relationship vector between the phenotype and all other phenotypes can be obtained from the phenotype layer . The phenotypes associated with the gene are obtained from the phenotype layer . The relationship vector between the gene and all other phenotypes is calculated using formula (8) . Thus, formula (9) is extended as follows: （10） In the linear regression model of the above formula, the relationship value between the phenotype and the gene can be calculated by using the Pearson linear correlation coefficient, and the formula is defined as follows: （11） in express and Covariance between Representing vectors standard deviation Representing vectors Standard deviation, correlation coefficient This reflects the genes in the gene layer. Phenotype in phenotypic layer The degree of correlation between them; (2) the relationship value calculated above is used as the initial value of the edge weight of RWR algorithm, and RWR algorithm can be expressed as: （12） in, Let be the restart probability, representing the probability that a random walk can return to the seed node at each step. The initial probability vector, i.e., the correlation coefficient. ; Let represent the probability vector of the random walk at time t; then the transition matrix of the two-layer network constructed from the gene layer and the phenotype layer is: ,in and These represent the intra-layer transfer matrices of the gene layer and the phenotype layer, respectively. and This is the transfer matrix between the gene layer and the phenotype layer. yes The transpose matrix, in solving During the process, the inter-layer network relationships are divided into two cases: related and unrelated, belonging to a binary network; for example, in In the matrix, a value of 1 indicates a correlation between a gene and its phenotype; otherwise, a value of 0 indicates a correlation. The matrix is ​​standardized using the following formula: （13） Similarly, （14） wherein, represents the i-th row, j-th column element in the adjacency matrix of the genotype layer and the phenotype layer.

4. The method according to claim 2 or 3, wherein, restart probability is 0.

7.

5. The method of claim 1, wherein the method is a multi-phenotype based genotypic and phenotypic association analysis method, characterized in that, Hyperparameters in step 2 is 0.7.

Images