Cell communication prediction method based on boosting, deep forest and single cell sequencing data and application

Abstract

The application discloses a cell communication prediction method and application based on Boosting, deep forest and single-cell sequencing data, and is characterized in that: on the basis of extracting ligand and receptor biological characteristics, a limit gradient boosting algorithm is designed to select the characteristics of ligand-receptor pairs; then, based on a category characteristic gradient boosting algorithm, a natural gradient boosting algorithm and a deep forest model, an integrated framework is designed to predict ligand-receptor interaction; and combined with single-cell sequencing data of tumor tissues, known and predicted ligand-receptor interaction data are filtered; and then, according to the filtered ligand-receptor interaction and single-cell sequencing data, combined with an expression product method and an expression threshold method, cell communication in a tumor microenvironment is predicted. The method can improve the prediction effect of cell communication, can be applied to cell communication prediction in human tumor tissues, and solves the problem that the accuracy of predicting cell communication intensity based on ligand-receptor interaction in the prior art is not high.

Description

Technical Field

[0001] This invention relates to the technical field of bioinformatics, and more specifically, to cell communication prediction methods and applications based on Boosting, deep forests, and single-cell sequencing data. Background Technology

[0002] In multicellular organisms, cell communication coordinates the activities of various cell types, leading to the formation of tissues, organs, and systems, and ultimately, the completion of various biological functions. Cell communication is also essential for complex bodily processes, such as immune responses, growth, and homeostasis under healthy or diseased conditions. To understand the biological functions of each cell type in its tissues, we need to understand the protein information transmitted by each type of cell.

[0003] Single-cell sequencing technology can accurately quantify the number of gene copies in a single cell nucleus. Since cancer cells often have parts of their genome deleted or amplified, leading to the loss or overexpression of key genes and interfering with normal cell growth, this method can be used to analyze gene copy numbers, thus having wide applications in cancer diagnosis. Single-cell sequencing often provides a large amount of gene data; identifying key intercellular relationships helps reveal the regulatory mechanisms between communicating cells, improving researchers' accuracy in predicting tissue function at homeostasis and disease changes. CN202011620086.X discloses a cell communication analysis method and system that uses cell communication prediction and ligand-target gene regulation prediction. Cell communication prediction includes ligand-receptor pair expression abundance analysis, ligand-receptor pair number analysis, significant enrichment of ligand-receptor pair number analysis, and construction of cell interaction network diagrams. Ligand-target gene regulation prediction includes ligand activity analysis and ligand-target gene regulatory potential analysis to describe intercellular relationships. Although this patent's cell communication analysis process is relatively efficient and comprehensive... However, this method has low performance, fails to visualize the prediction results, and lacks analysis of the tumor microenvironment. For the regulation of the interaction between secretory ligands and plasma membrane receptors in intracellular communication, i.e., ligand-receptor interaction, the accuracy of ligand-receptor interaction prediction is still limited. Summary of the Invention

[0004] The technical problem to be solved by this invention is to address the insufficient accuracy of existing cell communication prediction methods mediated by ligand-receptor interactions, and to provide a cell communication prediction method based on Boosting, deep forest, and single-cell sequencing data.

[0005] Another technical problem of the present invention is to provide an application of a cell communication prediction method based on Boosting, deep forest and single-cell sequencing data.

[0006] The objective of this invention is achieved through the following technical solution:

[0007] A cell communication prediction method based on Boosting, deep forest, and single-cell sequencing data, comprising the following steps:

[0008] S1. Extract biomarkers from the sequences of ligands and receptors, and use the limiting gradient algorithm to select the biomarkers for each ligand-receptor pair;

[0009] S2. The gradient boosting algorithm LRI-CatBoost is used to classify ligand-receptor pairs based on their biological characteristics.

[0010] S3. The LRI-NGBoost model based on natural gradient boosting is used to classify ligand-receptor pairs according to their biological characteristics;

[0011] S4. Using the deep forest algorithm, the biological characteristics of ligand-receptor are divided into positive and negative classes. The classes with higher probabilities are calculated and selected as the final classes.

[0012] S5. Filter known and predicted ligand-receptor interaction data;

[0013] S6. The final cell communication strength is calculated based on the filtered ligand-receptor interactions, single-cell sequencing data, and scoring methods.

[0014] Furthermore, the biometrics include 400-dimensional monoMono, 8000-dimensional monoDi, 8000-dimensional diMono, 147-dimensional CTD, and 80-dimensional PseudoAAC.

[0015] Furthermore, the limiting gradient algorithm is as follows:

[0016]

[0017] Where i is the i-th sample, I L g represents the number of samples in the left-hand node space. i h is the first-order partial derivative. i λ is the second-order partial derivative, and γ represents the regularization parameters.

[0018] Furthermore, the LRI-CatBoost gradient boosting algorithm classification steps include:

[0019] S21. Implement a symmetric decision tree using a top-down greedy algorithm. Each decision rule r consists of a feature i∈{1,..,l} and a threshold v∈R. At each level of the tree, decision rule r partitions k disjoint sets into 2k disjoint subsets. For a complete binary tree with level k', k=2k...k′ A set of feature vectors X∈R is divided into two completely independent subsets (Xi, Xj, X ... L and X R For each x∈X, LRI-CatBoost determines its category based on these two subsets:

[0020]

[0021] S22. Given a set and an objective function t:R l →R, the segmentation rule is defined as:

[0022]

[0023] Where M is used to evaluate X1,...,X k The optimality of the partitioning rule r on the surface;

[0024] S23. Obtain the prediction model M i,j M i,j (i) represents the permutation σ r The result of the i-th sample among the first j samples, in each iteration t, is from {σ1,...,σ...} S Construct a tree T t And calculate its gradient:

[0025]

[0026] S24. Calculate the gradient grad for each sample i. r,σ(i)-1 (i) After all possible action pairs have been predicted, the leaf value of sample i is calculated by taking the gradient grad of samples that previously belonged to the same leaf as sample i. r,σ(i)-1 The average value of (i) is used to construct the tree structure T. t Then, the unknown ligand-receptor pairs were classified.

[0027] Furthermore, M can be defined as:

[0028]

[0029] in Indicates about X i The target score set for the middle sample.

[0030] Furthermore, the LRI-NGBoost model consists of three parts: a base learner, a parameter probability distribution, and a prediction rule. For a sample x, LRI-NGBoost uses the conditional distribution P... θ Predict its label y, where the parameter θ is determined by the initial θ (0)This is achieved by combining the outputs of M base classifiers. For a normal distribution with parameters μ and logσ, each stage has two base classifiers. and therefore

[0031] Furthermore, for a sample x, LRI-NGBoost uses the conditional distribution P θ Predict its label y, where the parameter θ is determined by the initial θ (0) This is achieved by combining the outputs of M base classifiers. For a normal distribution with parameters μ and logσ, there are two base classifiers at each stage. and The predicted output is determined by a phased scaling factor p. (m) The learning rate η is used for evaluation, where the scaling factor ρ is the scaling factor. (m) It is a single scalar:

[0032]

[0033] Furthermore, random forest and additional trees are selected as base classifiers. For a ligand-receptor interaction feature, each predictor calculates the ratio of feature samples corresponding to the positive and negative classes in each layer. A class vector is generated from the class probabilities obtained from all predictors. This vector is concatenated with the original ligand-receptor interaction feature vector and used as the input to the next layer of deep forest.

[0034] When the prediction performance is better than all the previous layers, a new layer is added to the model; when the performance of the last two layers does not improve, training terminates. Finally, the average interaction probability is calculated for each ligand-receptor pair belonging to the positive and negative classes, and the class with the larger average interaction probability is taken as the final class.

[0035] Furthermore, the scoring method is a combination of the expression product method and the expression threshold method, and its cell communication score is calculated as follows:

[0036]

[0037] Where f1(k1,k2) is the cell communication score calculated based on the expression product method, and g1(k1,k2) is the cell communication score calculated based on the expression threshold method.

[0038] Furthermore, the cell communication score calculated based on the expression product method is:

[0039]

[0040] The cell communication score calculated based on the expression threshold method is:

[0041]

[0042] in, Cell types mediated by ligand i-receptor j interaction, calculated based on the expression product method. and The communication strength score, Cell types mediated by ligand i-receptor j interaction, calculated based on expression threshold method. and The communication strength score.

[0043] Furthermore, the present invention can also visualize the results of cell communication prediction.

[0044] The cell communication prediction method based on Boosting, deep forest, and single-cell sequencing data described above was applied to predict cell communication in human tumor tissues.

[0045] Compared with existing technologies, the beneficial effects are:

[0046] This invention employs a limiting gradient boosting algorithm to select features for ligand-receptor pairs based on extracted ligand and receptor biomarkers. Then, an ensemble framework is designed to predict ligand-receptor interactions based on categorical feature gradient boosting, natural gradient boosting, and a deep forest model. Furthermore, known and predicted ligand-receptor interactions are filtered using single-cell sequencing data, and expression product and expression threshold methods are combined to predict cell communication in the tumor microenvironment. The method described in this invention improves the predictive performance of cell communication. Attached Figure Description

[0047] Figure 1 This is a flowchart of cell communication prediction;

[0048] Figure 2 This is a flowchart of the ligand-receptor interaction prediction framework;

[0049] Figure 3 This is the AUC curve of the method described in this invention on datasets 1-4;

[0050] Figure 4 This is the AUPR curve of the method described in this invention in datasets 1-4;

[0051] Where a is dataset 1, b is dataset 2, c is dataset 3, and d is dataset 4;

[0052] Figure 5 This is a heatmap of cell communication ligand-receptor interactions in human head and neck squamous cell carcinoma tissue;

[0053] Figure 6This is a heatmap of cell communication intensity in human head and neck squamous cell carcinoma tissue;

[0054] Figure 7 This is a network diagram of cell communication intensity in human head and neck squamous cell carcinoma tissue;

[0055] Figure 8 This is a heatmap of cell communication ligand-receptor interactions in human breast cancer tissue.

[0056] Figure 9 This is a heatmap of cell communication intensity in human breast cancer tissue;

[0057] Figure 10 This is a network diagram of cell communication intensity in human breast cancer tissue. Detailed Implementation

[0058] The following examples further explain and clarify the invention, but the specific examples do not limit the invention in any way. Unless otherwise specified, the methods and equipment used in the examples are conventional methods and equipment, and the raw materials used are all commercially available conventional raw materials.

[0059] Example 1

[0060] like Figure 1-2 This embodiment provides a cell communication prediction method based on Boosting, deep forest, and single-cell sequencing data. The specific steps include:

[0061] S1. Biometrics are extracted from the ligand and receptor sequences, resulting in 400-dimensional monoMono, 8000-dimensional monoDi, 8000-dimensional diMono, 147-dimensional CTD, and 80-dimensional PseudoAAC. Each ligand or receptor can be described as a 16,627-dimensional vector, and a ligand-receptor pair can be represented as a 33,254-dimensional vector. The limiting gradient algorithm is used to select the biometrics for each ligand-receptor pair. The limiting gradient algorithm is as follows:

[0062]

[0063] Among them, I L This represents the number of samples in the left-hand node space. λ and γ represent the regularization parameters.

[0064] Higher feature gain implies more effective and more important features. After feature selection, each ligand-receptor pair is described as a d-dimensional vector.

[0065] S2. The gradient boosting algorithm LRI-CatBoost is used for classification based on the biological characteristics of ligand-receptor interactions;

[0066] Let D = (X, Y) denote a dataset with n ligand-receptor pairs, where x represents a training sample with a d-dimensional feature vector, and y ∈ Y represents its label. For the i-th ligand-receptor pair x... i If they interact, y i =1, otherwise y i =0.

[0067] A symmetric decision tree is implemented using a top-down greedy algorithm. Each decision rule r consists of a feature i∈{1,..,l} and a threshold v∈R. At each level of the tree, decision rule r partitions k disjoint sets into 2k disjoint subsets. Specifically, for a complete binary tree with level k', k=2... k′ A set of feature vectors X∈R is divided into two completely independent subsets (Xi, Xj, X ... L and X R For each x∈X, LRI-CatBoost can determine its category based on these two subsets:

[0068]

[0069] Therefore, based on the partitioning rule, any k mutually independent sets Both can be used to implement 2k independent sets.

[0070] When given a set and an objective function t:R l →R, the segmentation rule is defined as:

[0071]

[0072] Where M is used to evaluate X1,...,X k The optimality of the partitioning rule r on M. M can be defined as:

[0073]

[0074] in Indicates about X i The target score set for the middle sample.

[0075] Obtain the prediction model M i,j M i,j (i) represents the permutation σ r The result of the i-th sample from the first j samples. In each iteration t, from {σ1,...,σ...} S Construct a tree T t And calculate its gradient:

[0076]

[0077] For each sample i, its gradient grad r,σ(i)-1 (i) can be computed. Once all possible action pairs have been predicted, the leaf value of sample i can be calculated by computeding the gradient grad of samples that previously belonged to the same leaf as sample i. r,σ(i)-1 The average value of (i) is obtained. When the tree structure T t Once established, it allows for the classification of unknown ligand-receptor interaction data.

[0078] S3. The LRI-NGBoost natural gradient boosting model is used to predict the interaction probability of each ligand-receptor pair;

[0079] The LRI-NGBoost model consists of three parts: a base classifier (f), a parameter probability distribution (P), and a base classifier (f). θ The LRI-NGBoost algorithm uses the conditional distribution P and the prediction rule (S) to predict a given sample x. θ Predict its label y, where the parameter θ is determined by the initial θ (0) This is achieved by combining the outputs of M base classifiers. For a normal distribution with parameters μ and logσ, each stage has two base classifiers. and therefore The predicted output is determined by a phased scaling factor p. (m) The learning rate η is used for evaluation, where the scaling factor ρ is the scaling factor. (m) It is a single scalar:

[0080]

[0081] S4. Using the deep forest algorithm, the biological characteristics of ligand-receptor are divided into positive and negative classes. The class with the larger average interaction probability is calculated and selected as the final class.

[0082] Random forests and additional trees are chosen as base classifiers, with each cascaded layer consisting of two random forests and two additional trees. Each predictor consists of 100 decision trees. For a ligand-receptor interaction feature, each predictor calculates the ratio of feature samples corresponding to the positive and negative classes in each layer. The class probabilities obtained from all predictors generate a class vector. This vector is concatenated with the original ligand-receptor interaction feature vector and used as input to the next deep forest layer.

[0083] When the prediction performance surpasses that of all previous layers, a new layer is added to the model. Training terminates when the performance of the last two layers does not improve. Finally, the average interaction probability of each ligand-receptor pair belonging to the positive and negative classes is calculated separately. The class with the larger average interaction probability is selected as the final class.

[0084] Finally, we obtained the final classification of each ligand-receptor pair by integrating the results of LRI-CatBoost, LRI-NGBoost, and LRI-DF.

[0085] S5. Filter known and identified ligand-receptor interactions. If the ligand or receptor in a ligand-receptor interaction is not expressed in the cells of the single-cell sequencing data, the interaction is excluded from the corresponding cell communication.

[0086] S6. The final communication score is calculated based on the filtered ligand-receptor interactions, single-cell sequencing data, and scoring methods.

[0087] The scoring method employs a combination of the expression product method and the expression threshold method.

[0088] (1) Expression product method: predicting ligand i and receptor j with two cell types and The score of the interaction, where This indicates that ligand i and receptor j are in different cell types. The arithmetic mean in the middle:

[0089]

[0090] and The cell communication score f1(k1,k2) between cells can be calculated:

[0091]

[0092] (2) Expression threshold method: Predicting ligand i and receptor j with two cell types and The interaction score, where σ i and σ j Indicates standard deviation:

[0093]

[0094] and The cell communication score g1(k1,k2) between cells can be calculated:

[0095]

[0096] For calculations based on expression product method and expression threshold method and The cell communication scores f1(k1,k2) and g1(k1,k2) are combined to obtain the final cell communication score. In other words, and The cell communication score between cells can be calculated using the following formula:

[0097]

[0098] Example 2

[0099] This embodiment provides the cell communication prediction algorithm described in this invention and four representative protein interaction prediction methods, namely the extreme gradient boosting algorithm, support vector machine, distributed gradient boosting framework based on decision tree algorithm, and recurrent convolutional neural network algorithm based on ordinal regression. The performance was evaluated by 20 five-fold cross-validations, with AUC and AUPR as evaluation indicators. The higher the AUC and AUPR values, the better the algorithm performance.

[0100] The parameters of the extreme gradient boosting algorithm, support vector machine, and decision tree-based distributed gradient boosting framework are set to their default values. For the ordinal regression-based recurrent convolutional neural network algorithm, the parameters are set as follows: learning_rate = 0.01, n_estimators = 20, max_depth = 3, criterion = friedman_mse, loss = deviance, min_samples_split = 2. For the cell communication prediction algorithm proposed in this invention, the boosting type, max_depth, and n_estimators in LRI-CatBoost are set to Ordered, 10, and 2000, respectively; the learning rate, natural gradient, frac, and eval in LRI-NGBoost are set to 0.01, True, 1.0, and 100, respectively; n_trees in LRI-DF is set to 100, and the predictor is set to forest. The dimension of the ligand-receptor interaction feature vector after dimensionality reduction is set to 300.

[0101] In this experiment, we collected four different ligand-receptor interaction datasets. Datasets 1 and 2 were both from the CellTalk database. Dataset 3 was constructed by Skelly et al., and Dataset 4 was constructed by Ximerakis et al. The specific dataset details are shown in Table 1 below:

[0102] Table 1

[0103] Dataset 1 812 780 3390 Dataset 2 650 588 2031 Dataset 3 574 559 2006 Dataset 4 1129 1335 6585

[0104] The performance obtained using the different methods described above is shown in Table 2 below:

[0105] Table 2

[0106]

[0107] From Table 2 above and Figure 3-4 As can be seen, the ligand-receptor interaction prediction algorithm described in this invention achieved the best AUC and AUPR on the four datasets. The best AUCs on the four datasets were 0.8533, 0.8316, 0.8150, and 0.8434, respectively, which are 1.39%, 3.29%, 3.59%, and 1.89% higher than the second-ranked distributed gradient boosting framework based on the decision tree algorithm, respectively. Simultaneously, this invention also achieved the best AUPR on the four datasets, which were 0.8681, 0.8442, 0.8259, and 0.8632, respectively, which are 1.11%, 2.19%, 2.11%, and 1.54% higher than the second-ranked distributed gradient boosting framework based on the decision tree algorithm, respectively. The ligand-receptor interaction prediction algorithm LRI-CNbDP described in this invention achieved the best AUC and AUPR on the four datasets used in this experiment, demonstrating its powerful ligand-receptor interaction prediction performance.

[0108] Example 3

[0109] This embodiment provides a practical predictive application of the scheme described in this invention. Sequencing data of cells from human head and neck squamous cell carcinoma tissue, including cell types such as head and neck squamous cell carcinoma cells, fibroblasts, B cells, myocytes, macrophages, endothelial cells, T cells, dendritic cells, and mast cells, are downloaded from the GEO database. Combined with the filtered ligand-receptor interaction and single-cell sequencing data from this invention, a cell communication network related to breast cancer is established to predict cell communication in human tissues. For example... Figure 5-7 As shown, the method of the present invention has found that the communication strength between fibroblasts and human head and neck squamous cell carcinoma cells is higher in human head and neck squamous cells.

[0110] Example 4

[0111] This embodiment provides a practical predictive application of the scheme described in this invention. Relevant sequencing data from human breast cancer tissue are downloaded from the GEO database and combined with the filtered ligand-receptor interaction and single-cell sequencing data from this invention to establish a cell communication network related to breast cancer and predict cell communication in human tissues. For example... Figure 8-10As shown, the probability of communication between immune cells and breast cancer cells is higher in human breast cancer tissue.

[0112] Obviously, the above embodiments of the present invention are merely examples for clearly illustrating the present invention, and are not intended to limit the implementation of the present invention. Those skilled in the art can make other variations or modifications based on the above description. It is neither necessary nor possible to exhaustively describe all embodiments here. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included within the scope of protection of the claims of the present invention.

Claims

1. A cell communication prediction method based on Boosting, deep forest, and single-cell sequencing data, characterized in that, The method is applied to predict cell communication in the human tumor microenvironment, and the steps include: S1. Biometrics are extracted from the sequences of the ligand and receptor. The limiting gradient algorithm is used to select the biometrics for each ligand-receptor pair. The limiting gradient algorithm is as follows: in, For the first One sample, This represents the number of samples in the left-hand node space. It is the first-order partial derivative. It is the second-order partial derivative. and Represents the regularization parameter; S2. The gradient boosting algorithm LRI-CatBoost is used to classify ligand-receptor pairs based on their biological characteristics. S21. Use a top-down greedy algorithm to perform decision tree induction, for each decision rule. By a feature and a threshold Composition, at each level of the tree, decision rules Will A disjoint set is divided into A set of disjoint subsets, for a given set of subsets. A complete binary tree of level ' A set of feature vectors Divided into two completely independent subsets and Subset, for each LRI-CatBoost determines its category based on these two subsets: ； S22. When given a set and an objective function : The segmentation rule is defined as follows: in Used for evaluation Segmentation rules Optimality; S23. Obtain the prediction model ,in Indicates based on permutation Center front The first sample The results for each sample, in each iteration From Build a tree And calculate its gradient: ； S24. Calculate each sample gradient After all possible action pairs have been predicted, the results are calculated by comparing the previous results with the sample. Gradient of samples belonging to the same leaf The average value of the sample is obtained Use the values ​​of the leaf nodes to build a tree structure. Then, unknown ligand-receptor interactions are classified; S3. The LRI-NGBoost model based on natural gradient boosting is used to classify ligand-receptor pairs according to their biological characteristics; S4. Using the deep forest algorithm, the biological characteristics of ligand-receptor are divided into positive and negative classes. The classes with higher probabilities are calculated and selected as the final classes. S5. Filter known and predicted ligand-receptor interaction data; S6. The final cell communication strength is calculated based on the filtered ligand-receptor interactions, single-cell sequencing data, and scoring methods. The scoring method is a combination of the expression product method and the expression threshold method, and its cell communication score is calculated as follows: in, To express the cell communication fraction of the product method, Cell communication score based on expression threshold method; The cell communication score calculated based on the expression product method is: The cell communication score calculated based on the expression threshold method is: in, For calculation based on the expression product method, derived from ligands -receptor Interaction-mediated cell types and The communication strength score, For expression thresholding, calculated by ligands -receptor Interaction-mediated cell types and The communication strength score.

2. The cell communication prediction method based on Boosting, deep forest, and single-cell sequencing data according to claim 1, characterized in that, The biometrics include 400-dimensional monoMono, 8000-dimensional monoDi, 8000-dimensional diMono, 147-dimensional CTD, and 80-dimensional PseudoAAC.

3. The cell communication prediction method based on Boosting, deep forest, and single-cell sequencing data according to claim 1, characterized in that, It can be defined as: in express The target score set for the middle sample.

4. The cell communication prediction method based on Boosting, deep forest, and single-cell sequencing data according to claim 1, characterized in that, For a sample LRI-NGBoost uses conditional distribution Predict its label , where parameters From the initial and The combined implementation of the outputs of the base classifiers, for parameters of... and The normal distribution is followed, and there are two base classifiers at each stage. and , The predicted output is determined by a phased scaling factor. and learning rate To evaluate, where the scaling factor It is a single scalar: 。 5. The cell communication prediction method based on Boosting, deep forest, and single-cell sequencing data according to claim 1, characterized in that, Random forest and additional trees are chosen as the basic classifiers. For a ligand-receptor interaction feature, each predictor calculates the ratio of feature samples corresponding to the positive and negative classes in each layer. The class probabilities obtained from all predictors generate a class vector, which is concatenated with the original ligand-receptor interaction feature vector and used as the input to the next layer of deep forest. When the prediction performance is better than all the previous layers, a new layer is added to the model; when the performance of the last two layers does not improve, training terminates. Finally, the average interaction probability value is calculated for each ligand-receptor pair belonging to the positive and negative classes, and the class with the larger average interaction probability is taken as the final class.

Images