A protein signal network learning method based on parallel discrete artificial bee colony

Abstract

The application discloses a protein signal network learning method based on parallel discrete artificial bee colony, which combines an artificial bee colony algorithm and a K2 scoring method. The K2 scoring method is used to firstly perform discrete processing on multiple multi-parameter single-cell data sets input into the algorithm. Then, multiple bee colonies are initialized to perform parallel search on the multiple data sets after the discrete processing, and the causal action relationship between different types of phosphorylated protein biomolecules is learned. After each iteration of the parallel bee colonies is completed, the pheromone of the artificial bee colony is fused and updated to avoid falling into a local optimum. Finally, when a preset number of iterations is reached, according to a designed information fusion and extraction rule, multiple protein signal networks learned by all the bee colonies from the multiple data sets are fused together, and after screening, an optimal protein signal network is extracted. The application can effectively deduce a protein signal network, and provide a reference for medical researchers to analyze a signal communication mechanism in cells.

Description

TECHNICAL FIELD

[0001] The present application relates to the causal discovery and information fusion method of human immune system multi-parameter single-cell data, aiming to discover unknown effects in molecular biology, deepen the understanding of the interaction between different types of phosphorylated protein biomolecules, and assist in the diagnosis of cell diseases, and designs a protein signaling network learning method based on parallel discrete artificial bee colony. BACKGROUND

[0002] Protein signaling network is a network structure formed by the complex interaction between proteins involved in intracellular signal transduction. Protein signaling network plays a crucial role in cells, and understanding protein signaling network can help better understand the signal mechanism in cells. However, due to the complexity of the signal pathway, it is difficult to achieve this goal only through biological experiments. Fortunately, protein signaling network can be obtained in another way by analyzing variables such as protein or phospholipid expression levels in multi-parameter single-cell data. Learning protein signaling network from multi-parameter single-cell data set can be represented as a problem of detecting the interaction between different types of phosphorylated protein biomolecules. Specifically, protein signaling network is a directed acyclic graph composed of a set of nodes and directed edges, where each node represents a phosphorylated protein biomolecule, and each directed edge represents the causal relationship between two biomolecules. Therefore, constructing protein signaling network using biological data sets has become a frontier hotspot of current bioinformatics research.

[0003] In recent years, a large number of methods for learning protein signaling network from multi-parameter single-cell data set have emerged, including Bayesian network (BN) based method, linear non-Gaussian acyclic model (LiNGAM) based method, reinforcement learning (RL) based method and NoTears based on structure learning optimization method, etc.

[0004] Among the many methods, the Bayesian network-based method has been considered as one of the most promising frameworks for modeling complex systems (such as cell signaling cascades) and is a classic and effective method for inferring protein signaling networks. Bayesian networks are a form of graphical models and are one of the most effective theoretical models in the field of uncertain knowledge representation and reasoning at present, because they can represent the probabilistic dependency relationships between multiple interacting components. Since Bayesian networks are effective tools for building such multivariate relationships, the construction of protein signaling networks can be seen as a structure learning problem of Bayesian networks. Compared with other methods, the existing Bayesian network-based methods have the advantages of relatively simple model, not relying on a large amount of sample data, relatively low time cost, and strong interpretability. However, there are also many shortcomings. The existing Bayesian network-based protein signaling network learning methods are prone to fall into local optimum, and it is difficult to learn an accurate protein signaling network. Moreover, the time performance and learning accuracy of learning protein signaling networks from multiple data sets at the same time still need to be improved.

[0005] Therefore, how to accurately and quickly identify protein signaling networks from multiple multi-parameter single-cell data sets at the same time is still a challenging topic in this field. SUMMARY

[0006] In view of the problems of high time cost and low identification accuracy of the existing protein signaling network learning methods for learning protein signaling networks from multiple multi-parameter single-cell data sets at the same time, the present application proposes a protein signaling network learning method based on parallel discrete artificial bee colony (Parallel Discrete Artificial Bee Colony, PDABC). The method combines the artificial bee colony algorithm and the K2 scoring method. Since the K2 scoring method is used, like other Bayesian structure learning methods, the multiple multi-parameter single-cell data sets input into the algorithm need to be discretized first. Then multiple bee colonies are initialized to search the discretized multiple data sets in parallel, and the causal relationship between different types of phosphorylated protein biomolecules is learned. When each iteration of the parallel bee colony is completed, the pheromone of the artificial bee colony is fused and updated to avoid falling into local optimum. Finally, when the preset number of iterations is reached, according to the designed information fusion and extraction rules, the multiple protein signaling networks learned from multiple data sets by all bee colonies are fused together, and the final optimal protein signaling network is extracted after screening.

[0007] The main idea of realizing the present application is that the Bayesian network has been considered as the most promising method in inferring the protein signaling network, but the existing Bayesian network learning method has the disadvantages of low identification accuracy and high time cost when processing multiple large-scale data sets at the same time. Through the designed parallel artificial bee colony learning algorithm, firstly, multiple threads can be started, and multiple multi-parameter single-cell data sets can be searched in parallel, which can greatly improve the time performance of the algorithm. Secondly, through the designed information fusion strategy, the protein signaling network information learned by different bee colonies from different data sets can be exchanged and shared, which describes the complex relationship between different phosphorylated protein biomolecules. Through information sharing, the probability of the algorithm falling into a local optimal solution can be greatly reduced. Thirdly, after the algorithm iteration is completed, the optimal protein signaling network searched by each bee colony is fused together, and the final optimal protein signaling network is extracted according to the designed screening strategy, which can effectively solve the problem that small sample data cannot provide enough information, thereby further improving the accuracy of protein signaling network learning.

[0008] A protein signaling network learning method based on parallel discrete artificial bee colony includes the following steps:

[0009] Step one: data acquisition; in order to verify the effectiveness of the method proposed in the present application, experiments are carried out on the protein signaling network real data set provided by Sachs et al. The performance of the method for learning the protein signaling network is evaluated.

[0010] Step two: data discretization processing and initialization; K2 score is used to first discretize the multi-parameter single-cell data set, and then initialize multiple bee colonies and related parameters, each bee colony corresponding to a discretized data set.

[0011] Step three: parallel artificial bee colony searches for protein signaling network; K2 score is used to guide multiple artificial bee colonies to simultaneously search for and solve the protein signaling network of different data sets.

[0012] Step four: pheromone fusion and update; after each iteration is completed, the updated pheromone matrix of all bee colonies is fused and averaged, and the averaged pheromone matrix is updated to each population again, which is used to guide the search of the bee colony in the next iteration.

[0013] Step five: fusion and extraction of protein signaling network; steps four and five are repeated until the set number of iterations is completed, then the protein signaling network searched by each bee colony is fused into a total connection matrix. Finally, the final optimal protein signaling network is extracted according to the designed extraction strategy.

[0014] Compared with the prior art, the present application has the following obvious advantages and beneficial effects:

[0015] (1) combined with K2 score, a new protein signal network learning method based on parallel discrete artificial bee colony is proposed, and the designed parallel strategy can effectively improve the time performance of the algorithm.

[0016] (2) The pheromone fusion strategy designed by the application can share and exchange all the information learned by the bee colony from the whole data set, reduces the probability of falling into a local optimal solution, and effectively alleviates the problem that small sample data is difficult to provide enough data.

[0017] (3) The protein signal network fusion and extraction strategy designed by the application can optimize the obtained protein signal network, and improve the accuracy and performance of the algorithm learning.

[0018] (4) The experimental results on the real data set of protein signal network show that the application can effectively infer the protein signal network, and can provide reference for medical researchers to analyze the signal exchange mechanism in cells. BRIEF DESCRIPTION OF DRAWINGS

[0019] Figure 1 : The flowchart of the algorithm involved in the method.

[0020] Figure 2 : The protein signal network learned by PDABC on the multi-parameter single cell real data set. DETAILED DESCRIPTION

[0021] The specific implementation and detailed steps of the application will be described below, and the flow of the specific implementation of the application is shown in Figure 1 , which specifically includes:

[0022] (Step 1) Data acquisition.

[0023] In order to verify the effectiveness of the model proposed in the present application, in order to verify the effectiveness of the method proposed in the present application, experiments are carried out on the protein signal network real data set provided by Sachs et al. The performance of the method for learning the protein signal network is evaluated to verify the effectiveness of the method. The protein signal network real data set provided by Sachs et al. is widely recognized in the field of biological research, and the data set continuously measures the expression levels of a variety of phosphorylated proteins and phospholipids in human immune cells. The data set has 14 sub-data sets, involving 14 different biochemical experiments. The number of data points in each sub-data set is 723-917. Each sub-data set has 11 signal nodes, and each node represents a phosphorylated protein biomolecule. In the past two decades, a protein signal network containing 17 connection edges has been constructed for this data set through a large number of previous studies and classical biochemical and genetic analysis, which can be used as a standard protein signal network on this data set. The data set can be downloaded for free at https: / / www.science.org / doi / 10.1126 / science.1105809#supplementary-materials.

[0024] (Step two) Data discretization processing and algorithm initialization.

[0025] First, the data is discretized: data points that differ from the average value by more than three standard deviations are removed. Then, the data is discretized into three levels (low, medium or high level of phosphorylated protein) using clustering methods to minimize the loss of pairwise mutual information under chemical intervention conditions. After discretization, each part contains the same number of phosphorylated protein biomolecules, and each part contains the same number of data sample points.

[0026] Then, the algorithm is initialized. First, initialize N bee swarms containing K bees for parallel search of N multi-parameter single-cell sub-data sets, and then initialize an empty directed acyclic graph G(0) with the same number of nodes and a pheromone matrix M for each bee swarm, which are used to record the structure of the protein signal network and the change of pheromone in the search process, respectively. A directed acyclic graph can be represented by a joint probability distribution of a set of variables:

[0027]

[0028] In addition, the K2 score is used to evaluate the pros and cons of the learned protein signal network, and the initial expression of the K2 score is:

[0029]

[0030] where U is a given data set, G is a directed acyclic graph structure, c iIt is variable X i The number of possible values ​​that q can take. i It is Π(X) i The number of variables that can be configured in T) ijk It is U in X i The kth value, Π(X) i ) is the number of cases that are instantiated as the j-th value.

[0031] Finally, to prevent bees from getting stuck in an infinite loop while searching for a particular food source, a maximum digging limit and an algorithm iteration limit (Iter) are set for each population. max .

[0032] Step 3: Parallel artificial bee colony search for protein signaling networks.

[0033] After initialization, the parallel bee colony search for protein signaling networks begins. Each bee starts with an empty graph G(0) and gradually builds a connected solution by adding one arc at a time until the solution achieves the highest K2 score. The probability transition rule for bee k to select a directed edge from the current candidate edges at runtime t is defined as follows:

[0034]

[0035] Where, ε ij (t) and Θ ij (t) represents the pheromone concentration and the directional arc b, respectively. i，j The heuristic information, ξ, is the control Θ ij (t) Select the weighting coefficients for the arc; DA k (t) is the set of all candidate arcs with heuristic information greater than zero; q0 is the initial parameter that determines the relative importance of mining; I and J are a pair of nodes randomly selected with the following probabilities:

[0036]

[0037] Where parameter φ is the concentration ε of pheromones left by bees. rl To determine the relative importance of (t) and obtain the protein signaling network structure with the highest K2 score, the heuristic function of the directed arc is defined as follows:

[0038] Θ ij (t)=f(X i , Π(X i )∪X j )-f(X i , Π(X i (4)

[0039] Bees are divided into three types: mercenary bees, observer bees, and scout bees. The specific search process is as follows:

[0040] First, each employed bee is associated with an empty graph G(0) and then the employed bee searches for a new solution with a higher K2 score within the neighborhood of the current solution by testing and executing four basic evolutionary operators, which can be described as follows:

[0041] • Insert: The Insert operator inserts an arc between two randomly selected nodes in G(0) while ensuring that G(0) is a directed acyclic graph.

[0042] • Delete: The Delete operator deletes an arc between two randomly selected nodes.

[0043] • Reverse: The Reverse operator reverses the direction of an arc between two randomly selected nodes while ensuring that G(0) is a directed acyclic graph.

[0044] • Move: The Move operator deletes an arc between two randomly selected nodes, randomly selects a node from each of the parent node sets, and inserts a new arc between the selected nodes while ensuring that G(0) is a directed acyclic graph.

[0045] After all of the employed bees have completed the above process, they begin to share information about the solutions with the onlooker bees. The onlooker bees select their favorite solution according to the following probability:

[0046]

[0047] where f(G k ) is the K2 metric of the solution and K is the total number of solutions. According to this selection method, the higher the K2 value of a solution, the greater the probability of selection. Therefore, solutions with high K2 values can attract more onlooker bees than solutions with low K2 values.

[0048] After the onlooker bees have selected the solutions, the onlooker bees also enter the neighborhood search phase. In addition to the four basic evolutionary operators, each onlooker bee can also execute two knowledge-driven operators to search for a best solution with a higher K2 score, which will be the new position of the employed bee. The two basic operators can be described as follows:

[0049] • KD_Insert: This operator is an extended version of the Insert basic evolutionary operator. This operator inserts a new arc while ensuring that G(0) is a directed acyclic graph and that the K2 score will be higher after the graph structure is changed.

[0050] • KD_Delete: This operator is an extended version of the Delete basic evolutionary operator. This operator deletes a new arc and needs to guarantee that G(0) is a directed acyclic graph and the K2 score will be higher after the graph structure is changed.

[0051] If the current solution does not improve after limit iterations, the solution can be trapped in a local optimum. The employed bees will abandon the solution and turn into scout bees and then reconstruct new solutions according to the pheromone and heuristic information left on the arcs and become employed bees again. When the new states of all employed bees' solutions are determined, the best solution is the optimal protein signaling network for the current dataset searched by the swarm.

[0052] Let G + denote the current best solution found, f(G + : U) denote the K2 score of G + . At this time, Θ ij (t) is updated to preserve the relevant information of the best solution found in this iteration.

[0053] Θ ij (t) = (1 - p) Θ ij (t - 1) + p AΘ ij (t) (6)

[0054]

[0055] When all N parallel swarms complete the above process, the search phase of this iteration is completed and the algorithm enters the pheromone fusion and update phase.

[0056] Step Four: Pheromone Fusion and Update.

[0057] After Step Three is completed, the pheromone matrix of all swarms is updated, at this time, the pheromone matrix of all swarms is updated according to the following formula.

[0058]

[0059] The algorithm accumulates and sums all the pheromone matrices of the swarms and then takes the average matrix as the new pheromone matrix of each swarm to continue to guide the search of the swarm in the next iteration. At this time, the iteration process is completed. Repeat Step Three until Iter max iterations are completed, when Iter max iterations are completed, the algorithm performs Step Five.

[0060] Step Five: Fusion and Extraction of Protein Signaling Network.

[0061] After the iteration is completed, all the bee colonies obtain the optimal protein signal network corresponding to the searched data set. Each protein signal network is represented by an adjacency matrix containing 0 and 1. N adjacency matrices are accumulated to obtain an N*N matrix MG, and then MG is processed by designing protein signal network extraction rules. The specific extraction rules are as follows:

[0062] If the numerical element MG ij ≥ 40%*N, set to 1, otherwise, set to 0. Finally, the optimal protein signal network is obtained, and the algorithm ends.

[0063] In order to fully verify the superiority of the method, the PDABC of the present application is compared with six existing causal learning methods, namely DAG-GNN, ECDRL, PC, GAE and DABC. Six evaluation indexes, namely Precision, Recall, F1, Accuracy, SHD and Time, are used to evaluate the results. The six evaluation indexes are widely used in the evaluation of protein signal network learning. Among them, the values of Precision, Recall, F1 and Accuracy are between 0 and 1, and the larger the value, the higher the accuracy of the algorithm in learning the protein signal network. SHD is the structure hamming code distance, which represents the sum of the added edges, the missing edges and the reversed edges of the learned protein signal network compared with the standard protein signal network. The lower the value, the higher the accuracy of the model. Time represents the time required for the algorithm to run, and the unit is second. The less the time, the higher the performance of the algorithm.

[0064] Sachs et al. selected five sub-data sets from the data set as the input of the above six algorithms. The algorithm was repeatedly run for 10 times, and the average value of all index results was taken, as shown in Table 1.

[0065] Table 1: Comparison of each method on Sachs' experimental data set

[0066]

[0067] From Table 1, it can be found that the performance of the proposed algorithm is significantly higher than that of the other six comparison methods in the six evaluation indexes. Specifically, among the above algorithms, the deep learning method DAG-GNN and ECDRL learn the accuracy of the protein signal network significantly higher than PC and other traditional machine learning algorithms, but due to the use of too complex model structure, the time consumed is significantly more than other algorithms. The running time of GAE which also uses a deep learning model is 54.22 seconds, although the Precision reaches 0.5, but the Recall and F1 are 0.18 and 0.26 respectively, which shows that the algorithm identifies a large number of redundant edges, resulting in low overall accuracy. The performance of the PC algorithm in each evaluation index is the worst, and the SHD is as high as 28, which may be due to the setting of the acyclic constraint condition in the algorithm and the existence of the ring in the standard protein signal network, thereby leading to poor algorithm performance. Compared with the DABC algorithm which does not use parallel strategy and information fusion mechanism and all other algorithms, the running time of PDABC is the shortest, only 0.43 seconds, and the performance in each evaluation index is the best. The protein signal network learned by PDABC is shown in Figure 2 As shown in the figure, PDABC learns 9 correct signal connection edges, which is higher than all other algorithms, fully proving that the parallel strategy and information interaction mechanism of PDABC can effectively improve the accuracy and time performance of the algorithm.

[0068] The above experiments show that the method PDABC proposed by the present application has better performance and higher accuracy in learning protein signal networks on multiple sample multi-parameter single cell data sets compared with other methods, and therefore has great application prospect in computer-aided diagnosis of cell diseases.

Claims

1. A protein signal network learning method based on parallel discrete artificial bee colony, characterized by, The method comprises the following steps: Step one: data acquisition; experiments are carried out on the real data set of protein signaling network to evaluate the performance of learning protein signaling network; Step two: data discretization processing and initialization; K2 score is used to firstly discretize the multi-parameter single-cell data set, and then initialize multiple bee swarms and related parameters, each bee swarm corresponds to a discretized data set; Step three: parallel artificial bee colony search for protein signaling network; K2 score is used to guide multiple artificial bee colonies to simultaneously optimize and solve protein signaling network for different data sets; Step four: pheromone fusion and update; After each iteration, the updated pheromone matrix of all bee swarms is fused and averaged, and the averaged pheromone matrix is updated to each swarm again for guiding the search of the bee swarm in the next iteration; Step five: fusion and extraction of protein signaling network; steps four and five are repeated until the set number of iterations is completed, then the protein signaling networks searched by each bee swarm are fused into a total connection matrix; finally, the final optimal protein signaling network is extracted according to the designed extraction strategy; In step one, the real data set of protein signaling network has 14 sub-data sets, involving 14 different biochemical experiments; the number of data points in each sub-data set is 723-917; each sub-data set has 11 signal nodes, each node representing a phosphorylated protein biomolecule; In step two, firstly, the data is discretized: data points with a difference of more than three standard deviations from the average value are removed; then the data is discretized into three levels using clustering method; after discretization, each part contains the same number of phosphorylated protein biomolecules, and each part contains the same number of data sample points; Firstly, N bee swarms each containing K bees are initialized for parallel search of N multi-parameter single-cell sub-data sets, then an empty directed acyclic graph G(0) with the same number of nodes and a pheromone matrix M are initialized for each bee swarm to record the structure of the protein signaling network and the change of the pheromone in the search process respectively; the directed acyclic graph is represented by the joint probability distribution of a set of variables: In addition, K2 score is used to evaluate the pros and cons of the learned protein signaling network, and the initial expression of K2 score is: where U is a given dataset, G is a directed acyclic graph structure, c i is the number of possible values of the variable X i , q i is the number of possible configurations of variables in Π(X i ), T ijk is the kth value of X i in U, and ∏(X i ) is the number of cases that are instantiated to the jth value. Set a maximum number of digging times limit and the number of iterations of the algorithm Iter for each population max ; In step three, each bee starts from an empty graph G(0), gradually constructs a connected solution by adding an arc each time, until the solution obtains the highest K2 score; the probability transition rule for bee k to select a directed edge from the current candidate edge at running time t is defined as: where ε ij (t) and Θ ij (t) are the pheromone concentration and the heuristic information of the arc b i,j , respectively, ξ is the weighting factor that controls the selection of Θ ij (t) ; DA k (t) is the set of candidate arcs for which the heuristic information is greater than zero; q0 is an initial parameter that determines the relative importance of exploitation versus exploration, and I and J are a pair of nodes chosen randomly with probabilities given by where the parameter φ is the pheromone concentration ε left by the bees rl The relative importance of (t) for obtaining the highest K2 score for the protein signal network structure is defined by the heuristic function of the directed arcs as follows: Θ ij (t) = f(X i ,∏(X i )∪X j )-f(X i ,∏(X i )). (4) The bees are divided into three types: employed bees, onlookers and scouts, and the specific search process is as follows: Firstly, each employed bee is associated with an empty graph G(0), and then the employed bee finds a new solution with a higher K2 score in the neighborhood of the current solution by testing and executing four basic evolution operators, which can be described as: Insert: under the premise that G(0) is a directed acyclic graph, the Insert operator randomly selects two different nodes in G(0), and then inserts an arc between the two nodes; • Delete: Delete operator selects a random arc between two random nodes and deletes it; • Reverse: Reverse operator selects a random arc between two random nodes and reverses the direction of the arc, provided that G(0) is a DAG; • Move: Move operator selects two random nodes, deletes the arc between them, and inserts a new arc between two randomly selected nodes from their parent sets, provided that G(0) is a DAG; After all the employed bees complete the above process, they start to share the information of the solutions with the onlooker bees; the onlooker bees choose their favorite solution according to the following probability where f(G k ) is the K2 metric of the solution, K is the total number of solutions; according to this selection method, the higher the K2 value of a solution, the greater the probability of being selected; solutions with high K2 values attract more bystander bees than solutions with low K2 values. 2.The parallel discrete artificial bee colony based protein signaling network learning method according to claim 1, wherein, After the onlooker bees choose the solutions, the onlooker bees enter the neighborhood search phase; each onlooker bee can perform two knowledge-driven operators in addition to the four basic evolutionary operators to search for a better solution with a higher K2 score; the better solution will be the new position of the employed bee; the two basic operators are described as follows: • KD_Insert: is an extended version of the Insert basic evolutionary operator; this operator inserts a new arc, which not only needs to ensure that G(0) is a DAG but also needs to ensure that the K2 score will be higher after the graph structure is changed; • KD_Delete: is an extended version of the Delete basic evolutionary operator; this operator deletes a new arc, which not only needs to ensure that G(0) is a DAG but also needs to ensure that the K2 score will be higher after the graph structure is changed; If the current solution does not improve after limit iterations, the solution may be trapped in a local optimum; the employed bee will abandon the solution and become a scout bee, then reconstruct a new solution according to the pheromone and heuristic information left on the arcs, and then become an employed bee again; when the new states of all the solutions of the employed bees are determined, the best solution is the optimal protein signaling network corresponding to the current search data set of the bee colony; Let G + denote the set of all possible solutions; f(G + ) denotes the best solution found so far; and f(G + :U) denotes the K2 score of G . At this time, update Θ ij (t) with the relevant information that preserves the optimal solution found in this iteration; Θ ij (t) = (1 - p) Θ ij (t - 1) + p ΔΘ ij (t) (6) When the N parallel search bee colonies complete the above process, the search phase of this iteration ends, and the algorithm enters the pheromone fusion and update phase. 3.The parallel discrete artificial bee colony based method for learning protein signaling networks according to claim 1, wherein, In step four, the pheromone matrix of all bee colonies is updated; the pheromone matrix of all bee colonies is updated according to the following formula; The algorithm adds up all the pheromone matrices of the population, then takes the average matrix and uses it as the new pheromone matrix for each bee colony to continue guiding the search in the next iteration of the bee colony; at this point, the current iteration process is complete; repeat step three until Iter max max When the Iter max When the Iter max When the Iter 4.The parallel discrete artificial bee colony based protein signaling network learning method according to claim 2, wherein, In step five, after the iteration is completed, all bee colonies obtain the optimal protein signaling network corresponding to the search data set; each protein signaling network is represented by an adjacency matrix containing 0 and 1; N adjacency matrices are accumulated to obtain an N x N matrix MG, and then the protein signaling network extraction rule is designed to process MG; the extraction rule is as follows: If the numerical element in MG is MG ij ≥40%×N, set to 1; otherwise, set to 0; finally, the optimal protein signaling network was obtained, and the process ended.

Images