A cell-cell communication prediction method based on heterogeneous graph integration and dynamic fusion

Abstract

The application belongs to the field of bioinformatics and relates to a cell communication prediction method based on heterogeneous graph integration and dynamic fusion. First, ligand receptor features and labels are obtained and a heterogeneous graph is constructed, outer layer stratified cross-validation and training fold stratified validation are adopted, and only training feature standardization is based; second, a deep neural network branch outputting a prediction logarithm value is trained, and a first prediction probability is obtained through Sigmoid; third, a multi-random seed heterogeneous graph auto-encoding sub-model is trained, an AUC based on a current outer layer reserved evaluation set is monitored, and the optimal parameters are rolled back, the sub-model edge prediction probability is averaged to obtain a second prediction probability; then, the two branch results are ranked and normalized, and dynamic gating or weight search fusion is performed based on the AUC of the current evaluation set; finally, the final interaction prediction score or ranking result is output. The application can reduce the risk of information leakage, alleviate the output scale difference and weak branch interference, and significantly improve the model stability and generalization ability.

Description

Technical Field

[0001] This invention belongs to the field of bioinformatics, specifically relating to a method for predicting intercellular communication based on heterogeneous graph integration and dynamic fusion, and more particularly to a method that combines deep neural network branches and heterogeneous graph autoencoder branches, and achieves ligand-receptor interaction prediction through rank normalization and dynamic gating fusion based on the current evaluation set AUC. Background Technology

[0002] Intercellular communication is a crucial molecular mechanism for maintaining growth, differentiation, immune responses, and disease progression in multicellular organisms, with ligand-receptor interactions typically serving as its core carrier. With the development of single-cell sequencing technology, accurately identifying potential ligand-receptor interactions from high-dimensional expression data has become a key challenge in cell communication analysis, tumor microenvironment research, and precision medicine.

[0003] Existing methods for predicting cell communication mainly include those based on expression co-occurrence statistics, those based on sequence feature learning, and those based on graph neural network topology modeling. Among these, single sequence feature methods can model the local properties of ligands and receptors, but they are difficult to fully utilize the global network topology. Single graph structure methods, while able to capture the association patterns between ligands and receptors, still have shortcomings in node attribute discrimination.

[0004] To address the aforementioned issues, existing ensemble prediction frameworks attempt to jointly model deep neural network branches with heterogeneous graph autoencoder branches, thereby simultaneously utilizing node attribute features and network topology features. However, existing ensemble methods still have the following shortcomings: First, the data partitioning method is not sufficiently adaptable to class imbalance, leading to significant fluctuations in evaluation results between different folds; second, the training objectives and output formats of deep neural network branches are not well-suited for binary classification AUC optimization, easily causing unstable probability characterization; third, the lack of dynamic learning rate adjustment and optimal parameter rollback mechanisms based on validation set performance limits the model's generalization ability; fourth, heterogeneous graph branches typically use single random initialization for training, resulting in sensitivity to initialization and significant result fluctuations; fifth, when directly fusing the prediction results of two branches with fixed weights, the fusion effect is easily limited due to inconsistencies in the output scales of the two branches and differences in the performance of individual branches.

[0005] Therefore, there is an urgent need in this field for a technical solution that can further improve the prediction accuracy and generalization stability of cell communication ligand receptor interactions by maintaining the advantages of the dual-branch integrated architecture, through hierarchical data partitioning, stable deep neural network training, heterogeneous graph integration, and AUC-based dynamic fusion mechanism. Summary of the Invention

[0006] The purpose of this invention is to provide a method for predicting intercellular communication based on heterogeneous graph integration and dynamic fusion, which solves the problems of unstable data partitioning, insufficient generalization of deep neural network training, sensitivity of heterogeneous graph models to initialization, and limited effect of dual-branch fixed weight fusion in the prior art.

[0007] To achieve the above objectives, the present invention adopts the following technical solution: a method for predicting intercellular communication based on heterogeneous graph integration and dynamic fusion, comprising the following steps: Step 1: Data Acquisition, Hierarchical Partitioning, and Standardization. Obtain the raw expression feature data and known interaction labels of the ligand-receptor pairs. Construct a heterogeneous graph based on the ligand-receptor pairs, where the heterogeneous graph data includes at least ligand nodes, receptor nodes, and their corresponding interaction edges. Using an outer-layer hierarchical cross-validation method, while maintaining a consistent positive and negative sample distribution ratio, divide the ligand-receptor pairs into training folds and test folds. Further partition the training folds into validation sets, and perform standardization processing on the training features, validation features, and test features input to the deep neural network branches based on the feature statistics corresponding to the training folds to obtain standardized feature data.

[0008] Step 2: Deep Neural Network Branch Training and First Prediction Result Generation. The standardized feature data obtained in Step 1 is input into the deep neural network branch, and the branch is trained based on the corresponding interaction labels and validation set. The output layer of the deep neural network branch does not have a Sigmoid activation function; it directly outputs the predicted logarithmic value of the ligand-receptor pair. Based on the training features and validation set, a binary cross-entropy loss function including Sigmoid transformation is used, combined with weight decay, gradient pruning, and dynamic learning rate scheduling based on the validation set AUC, to train the deep neural network branch. During training, the model parameters corresponding to the optimal validation set AUC are recorded. When the validation set AUC does not improve within a preset number of consecutive training rounds, training stops, and the deep neural network branch is rolled back to the specified model parameters. Subsequently, a Sigmoid transformation is performed on the predicted logarithmic value to obtain the first prediction probability result.

[0009] Step 3: Training the Heterogeneous Graph Autoencoder Ensemble Branch and Generating the Second Prediction Result. The heterogeneous graph data obtained in Step 1 is input into the heterogeneous graph autoencoder branch, and training is monitored based on the current outer cross-validation retained evaluation set. Multiple heterogeneous graph autoencoder sub-models are initialized using multiple different and reproducible random seeds. For each heterogeneous graph autoencoder sub-model, training is monitored based on the retained evaluation set AUC in the current outer cross-validation. When the retained evaluation set AUC does not improve within a preset number of consecutive training rounds, the training of the corresponding sub-model is stopped, and the parameters of the corresponding sub-model are rolled back to the training round with the optimal retained evaluation set AUC. The edge prediction probabilities output by multiple heterogeneous graph autoencoder sub-models are averaged according to the sample dimension to obtain the second prediction probability result.

[0010] Step 4: Rank normalization and dynamic fusion based on the retained evaluation set AUC in the current outer layer cross-validation. Rank normalization is performed on the first and second prediction probability results respectively, and the AUC of the deep neural network branch and the heterogeneous graph autoencoder branch on the current evaluation set is calculated respectively. When the current evaluation set AUC corresponding to the first prediction probability result is lower than the current evaluation set AUC corresponding to the second prediction probability result, and the difference between the two exceeds a preset threshold, the second prediction probability result is used as the final interactive prediction result. Otherwise, with the goal of maximizing the current evaluation set AUC, a dynamic fusion coefficient is searched within a preset fusion weight range where the weight of the second prediction probability result is greater than the weight of the first prediction probability result. The normalized first and second prediction probability results are then linearly weighted and fused according to the dynamic fusion coefficient to obtain the final interactive prediction score.

[0011] Step 5: Output the interaction prediction results. Based on the final interaction prediction results, output the interaction prediction score, interaction prediction probability, or ranking result of the candidate ligand-receptor pairs.

[0012] The outer-layer stratified cross-validation method in Step 1 is StratifiedKFold, a stratified K-fold cross-validation algorithm used to maintain a consistent ratio of positive to negative samples when dividing training and testing folds. Within each training fold, a stratified partitioning method based on label distribution consistency is further used to construct the internal training and validation sets of the deep neural network branch. The number of folds in the outer-layer stratified cross-validation is 5, and the internal validation set partition ratio is 15%. The standardization process in Step 1 calculates the mean and standard deviation of each feature dimension based only on the features of the internal training set, and uses the mean and standard deviation to perform zero-mean, unit-variance transformation on the features of the internal training set, validation set, and outer-layer test fold to avoid leakage of validation set or test fold information into the training process.

[0013] In step 2, the deep neural network branch employs a multi-layer fully connected neural network to learn the representation of the standardized ligand-receptor pair features. The output layer of this deep neural network branch does not use a sigmoid activation function; instead, it directly outputs the predicted logarithmic value. During training, a binary cross-entropy loss function (BCEWithLogitsLoss) incorporating a sigmoid transform and log-probability is used, along with explicit weight decay, gradient pruning, dynamic scheduling of the learning rate based on the validation set AUC using an evaluation metric-based adaptive learning rate, and early stopping and optimal parameter rollback mechanisms based on the validation set AUC. After training, the deep neural network branch is restored to the model parameters corresponding to the highest validation set AUC, and a sigmoid transform is performed on the predicted logarithmic value to obtain the first predicted probability result.

[0014] The heterogeneous graph autoencoder branch in step 3 consists of a graph encoder and a bilinear decoder, used to learn the heterogeneous graph topological representation formed by ligand nodes, receptor nodes, and their interaction edges. Multiple heterogeneous graph autoencoder sub-models are initialized and trained independently using different and reproducible random seeds. Each sub-model has a patience value mechanism that retains the evaluation set AUC based on the current outer layer cross-validation. When the AUC does not improve within a preset number of consecutive rounds, the training of the corresponding sub-model is terminated, and the parameters of that sub-model are rolled back to the parameter state corresponding to its highest AUC round. After multiple sub-models have completed training and rollback, the arithmetic mean of their output edge prediction probabilities is performed to obtain the second prediction probability result.

[0015] The ranking rank normalization process in step 4 involves generating rank values ​​based on the ranking positions of the first and second prediction probability results in their respective result sets, and mapping these rank values ​​to rank ratio values ​​within a preset interval to reduce the impact of the output scale difference between the two branches on the fusion result. During dynamic fusion, the AUC of the deep neural network branch and the heterogeneous graph autoencoder branch on the current evaluation set is compared first. When the AUC of the deep neural network branch is significantly lower than that of the heterogeneous graph autoencoder branch, dynamic gating is triggered, and the output of the heterogeneous graph autoencoder branch is directly used as the final prediction score. Otherwise, candidate fusion coefficients are searched within a preset weight interval where the heterogeneous graph autoencoder branch dominates, and the optimal fusion coefficient is determined with the maximum AUC of the current evaluation set as the optimization objective. Then, the rank normalization results of the two branches are linearly weighted and fused according to the optimal fusion coefficient to obtain the final interactive prediction score.

[0016] Compared with the prior art, the present invention has at least the following beneficial effects: First, by using hierarchical cross-validation and hierarchical validation set partitioning, the distribution of positive and negative samples during the training and evaluation processes is made more stable, reducing performance fluctuations between different experimental rounds from the source.

[0017] Secondly, by adjusting the output form of the deep neural network branches to predict log values, and combining BCEWithLogitsLoss, dynamic learning rate scheduling, gradient pruning, and optimal parameter rollback mechanisms, the training stability and generalization ability of the deep neural network branches in binary classification tasks are improved.

[0018] Third, by adopting a multi-random seed heterogeneous graph autoencoder model integration and an optimal loop rollover mechanism, the sensitivity of heterogeneous graph branches to single random initialization is reduced, and the stability of graph topology learning results is improved.

[0019] Fourth, by performing rank normalization on the dual-branch output and introducing a dynamic gating and dynamic weight fusion strategy based on the current evaluation set AUC, the performance loss caused by fixed weight fusion and inconsistent output scale is effectively avoided, thereby further improving the final prediction performance.

[0020] Fifth, in a specific example, the average AUC of the present invention on the standard dataset can reach 0.881, which is a further improvement compared to the AUC of about 0.861 of the original ensemble method, indicating that the present invention has better predictive performance and better experimental stability in the cell communication ligand receptor interaction prediction task. Attached Figure Description

[0021] Figure 1 This is a flowchart illustrating the overall process of an intercellular communication prediction method based on heterogeneous graph integration and dynamic fusion according to the present invention.

[0022] Figure 2 This is a flowchart of the data partitioning and preprocessing process.

[0023] Figure 3 This is a schematic diagram of the training process for a branch of a deep neural network.

[0024] Figure 4 This is a schematic diagram of the training and rollback process for the heterogeneous graph autoencoder ensemble branch.

[0025] Figure 5 This is a schematic diagram of the fusion process of rank normalization and dynamic gating based on the current evaluation set AUC.

[0026] Figure 6 This is a schematic diagram of a heterogeneous graph data structure. Detailed Implementation

[0027] The present invention will be described in detail below with reference to the accompanying drawings and specific embodiments, but the scope of protection of the present invention is not limited thereto.

[0028] This invention provides a method for predicting intercellular communication based on heterogeneous graph ensemble and dynamic fusion. This method uses ligand-receptor pairs as the basic prediction objects, and combines node attribute features and network topology features for bi-branch modeling. Finally, it obtains the predicted score or ranking result of candidate ligand-receptor interactions through rank normalization and dynamic gating fusion based on the current evaluation set AUC. The overall process can be found in [link to documentation]. Figure 1 For the data partitioning and preprocessing process, please refer to [link / reference]. Figure 2 For the deep neural network branch training process, please refer to [link / reference]. Figure 3 For the heterogeneous graph autoencoder ensemble branch training and rollback process, please refer to [link to relevant documentation]. Figure 4 For integrated decision-making processes, please refer to Figure 5 The heterogeneous graph data structure can be found in [reference]. Figure 6 .

[0029] This invention proposes a method for predicting intercellular communication based on heterogeneous graph integration and dynamic fusion. The overall system consists of a cascaded input layer, a deep neural network sequence learning layer, multiple sub-heterogeneous graph topology learning layers, a dynamic fusion and evaluation layer, and a prediction output layer. The specific steps are as follows: Step 1: Data Acquisition, Hierarchical Partitioning, and Standardization. This step aims to provide a stable, information-leak-free training and validation data stream for the two-branch model, eliminating distribution oscillations during model evaluation by strictly controlling the ratio of positive to negative samples.

[0030] The data input stage of this embodiment includes two complementary data sources: the first is the node feature matrix, where the ligand feature matrix is ​​denoted as... The receptor feature matrix is ​​denoted as Original feature dimension The first category consists of 1790 dimensions, which are mapped to fixed low-dimensional continuous vectors after nonlinear dimensionality reduction; the second category consists of ligand-receptor interaction-associated tag data, represented as a one-dimensional array of binary tags. .

[0031] To prevent model evaluation oscillations during cross-validation training of the multimodal joint model due to positive and negative sample distribution shifts, this method designs a double-layered nested sampling topology at the starting point of the data flow: Outer Layer Stratified Five-Fold Cross-Validation Loop: Instantiate the outer layer scheduling operator StratifiedKFold(n_splits=5, shuffle=True, random_state=random_seed). This operator, when splitting the full association matrix samples, enforces a constraint on the ratio of positive to negative samples in each outer training fold and outer test fold, making it consistent with the ratio of positive to negative samples in the global dataset. The distribution is strictly consistent.

[0032] Inner-layer hierarchical validation set loop interception: Within each outer training layer, the features and labels are split a second time using the inner-layer hierarchical scheduling function: train_test_split(feature_train, target_train, test_size=0.15, stratify=target_train, random_state=seed). By configuring class distribution consistency constraints, it is ensured that the split independent validation set and the new training set have the same distribution pattern in the class imbalance scenario. This serves as the internal validation benchmark for deep neural network branch learning rate scheduling, early stopping, and parameter rollback. The heterogeneous graph autoencoder branch and fusion module monitor AUC based on the evaluation set retained by the current outer-layer cross-validation.

[0033] At the input of a multilayer fully connected network, to eliminate dimensional differences between different biometric features and prevent future information leakage, standardized operations are subject to strict sequential constraints on the data flow: This method calculates the mean values ​​of each dimension of the features by performing state fitting and transformation, relying solely on the inner layer's new training set data stream. and standard deviation Subsequently, the statistical metric is frozen in memory, and stateless zero-mean and unit-variance scaling is performed on the validation set data stream and the external test set data stream using this frozen metric. This code-level temporal isolation completely prevents the reverse infiltration of test-phase distribution information into the training phase.

[0034] Step 2: Deep Neural Network Branch Training and First Prediction Result Generation. This step aims to enhance the generalization ability of the neural network in extracting local sequence features in binary classification tasks by reconstructing the deep network output mechanism, introducing gradient protection, and high-level learning rate scheduling.

[0035] In this embodiment, the deep neural network branch can use an existing multi-layer fully connected network structure as the basic model. Unlike simply adjusting the network layer structure, this invention further introduces stabilizing training mechanisms for this branch, including hierarchical validation, input feature standardization, predicted logarithmic output, BCEWithLogitsLoss loss function, learning rate scheduling based on validation set AUC, early stopping and rollback, gradient pruning, and explicit weight decay. These mechanisms aim to improve the training stability and generalization ability of the deep neural network branch in the ligand-receptor interaction binary classification task.

[0036] The sequence feature extraction branch of this invention relies on a deep neural network model. In the underlying network structure declaration, this network eliminates the nonlinear probability mapping activation layer at the end. The forward function of the network directly outputs the original predicted logits after dimensionality compression of the multi-layer fully connected network block.

[0037] For this unconstrained Logits output stream, the target loss function nn.BCEWithLogitsLoss() with linked modification is configured in the training script. At the underlying level, this loss function performs a unified mathematical combination of the exponential power transformation of Sigmoid and the cross-entropy formula, eliminating the numerical collapse risk of logarithmic calculation overflow caused by extremely small probability values in deep backpropagation.

[0038] To effectively suppress the overfitting risk of the deep neural network branch during the training life cycle, the following control flow is embedded in the iterative optimization loop of the model: Explicit regularization and gradient clipping: The deep neural network branch uses the Adam optimizer for parameter updates with a learning rate of 0.001, and the explicit weight decay parameter weight_decay = . After the loss function performs backpropagation loss.backward() and before the optimizer takes a step, by calling nn.utils.clip_grad_norm_(model0.parameters(), max_norm = 5.0), global norm constraint clipping is performed on the parameter gradients of each layer of the network, forcing the sudden abnormal gradient norm to be truncated within the preset upper threshold (such as 5.0).

[0039] Learning rate plateau scheduling based on validation metrics: Deploy the ReduceLROnPlateau learning rate adaptive scheduler, set mode='max', factor = 0.5, patience = 15, and use the validation set AUC as the scheduling metric. If the validation set evaluation metric does not show a positive improvement within the preset consecutive rounds, the scheduler automatically multiplies the current learning rate by the preset decay factor in memory for stepwise compression.

[0040] Adaptive parameter early stopping and peak weight rollback: Set the hard state counter wait and the early stopping threshold patience = 40. During the model operation, if the AUC metric of the current round of the validation set is better than the historical peak, that is, it satisfies the condition val_auc best_val_auc , then reset the counter wait = 0, and temporarily store the current optimal network state best_state through the underlying copy instruction copy.deepcopy(model0.state_dict()); if there is no performance improvement after reaching the early stopping threshold for consecutive rounds, the training loop is interrupted in advance. At the code boundary of exiting the loop, forcefully call the parameters in the temporary storage area to perform an accurate rollback on the neural network, and restore its parameters to the state of the round with the best generalization performance of the validation set.

[0041] After the model state rollback is complete, the external test fold feature vector is input into the network branch. After performing an explicit Sigmoid() nonlinear transformation on the original logarithmic values ​​Logits, the first set of predicted probability results is continuously distributed in the interval (0,1).

[0042] Step 3: Heterogeneous graph autoencoder ensemble branch training and second prediction result generation. This step aims to capture the topological relationships between ligands and receptors in macroscopic biological networks, and eliminate the sensitivity of graph aggregation operations to network structure initialization and random walk paths through an ensemble learning mechanism. To address the pain point of heterogeneous graph neural networks containing graph encoders (GraphEncoder) and bilinear decoders (BilinearDecoder) being sensitive to random evolution paths and message aggregation network topology initialization, this method constructs a parallel ensemble architecture that includes multiple models, with the preferred number of ensembles being hgae_ensemble=3.

[0043] To ensure that each subnetwork has a completely different search path and to make the experiment reproducible, each independent subnetwork is assigned a distinguishable and controllable independent offset random seed at the code level: ens_seed = random_seed + i 100+ ens_idx, where random_seed is the base seed. i `ens_idx` is the fold index for cross-validation, and `ens_idx` is the ensemble index for the current sub-network. At the beginning of each sub-network training loop, strict differentiation and reproducibility constraints on random walks and neighborhood feature aggregation paths are achieved by resetting the underlying random states: `mx.random.seed(ens_seed)`, `np.random.seed(ens_seed)`, and `random.seed(ens_seed)`.

[0044] After instantiating a network instance containing GraphSAGE aggregation operators, explicitly declare the regularization control term {'learning_rate'=0.001, 'wd'= in the optimizer configuration dictionary. }, using parameter value The high-intensity weight decay mechanism suppresses the overfitting tendency in graph topology learning. The feature flow of each subnetwork when updating the node state is as follows: the node's own representation is concatenated with the longitude-normalized neighbor aggregation representation, and then the updated node embedding is obtained through linear mapping, LeakyReLU and Dropout.

[0045] Within the set long training cycle, each sub-model independently initiates stop monitoring based on the AUC of the current outer layer cross-validation edge: It allocates an independent state dictionary structure `best_hgae_params` in memory; when the AUC of the retained evaluation set in the current round reaches a new high, it temporarily stores a copy of the weights of each network parameter item in real time using dictionary derivation, and sets an early stopping patience value `hgae_patience=30` specifically for heterogeneous graph autoencoders. If the AUC of the retained evaluation set does not show a significant improvement for 30 consecutive rounds, a break is triggered to interrupt the sub-model optimization. After exiting the current sub-network training stream, it iterates through the parameter items using pointers and calls `param.set_data(best_hgae_params[name])` to perform optimal loop rollback, accurately restoring the state of each sub-network to its own historical best performance node.

[0046] After all subnetworks have rolled back to their respective extreme generalization states, test edge feature matrix inference is performed in parallel, generating multiple independent one-dimensional arrays of graph probabilities. Finally, the mean ensemble operator np.mean(np.vstack(hgae_scores),axis=0) is used on the backbone data stream to perform element-wise summation and mean averaging of the prediction output matrices of the multiple subnetworks along the test pair dimension, outputting a second set of prediction probability results after eliminating topological perturbations.

[0047] Step 4: Ranking rank normalization and dynamic fusion based on the retained evaluation set AUC in the current outer cross-validation. This step aims to address the technical bottleneck of severe modal heterogeneity in the probability sets output by the sequence learning branch and the graph topology branch, and how direct weighting can easily lead to large numerical branches suppressing small numerical branches.

[0048] Because the model architecture, loss objective, and physical logic of the two-branch model are fundamentally different, the absolute numerical distributions of the first and second prediction probability result sets typically exhibit significant modal heterogeneity. To eliminate the strong suppression of the smaller numerical branch by the larger numerical branch during weighting, a normalized rank transformation algorithm is deployed at the entry point of the fusion module. Its execution logic is as follows: Stable sorting indexing: Receives the original real number array, calls the underlying stable merge sort instruction np.argsort(scores,kind='mergesort'), and outputs a sorted index array of sample points arranged in ascending order while ensuring that the original relative positions of elements with equal values ​​remain unchanged.

[0049] Spatial rank proportionalization mapping: Allocate a new floating-point memory container ranks, and let ranks[order]=np.arange(scores.shape[0]) to record the absolute relative ranking of each sample in its set; finally, execute the normalization mapping equation: , the discrete predicted absolute scores are losslessly transformed into relative rank proportion features within the range of 0.0 to 1.0. Thus, score0 is transformed into score0_norm, and score1 is transformed into score1_norm, completing the same-scale transformation of heterogeneous features based on dimensionless and distribution alignment.

[0050] Subsequently, the evaluation metrics AUC_DNN and AUC_HGAE of the two branches on the current evaluation set are calculated, and the fusion decision-maker triggers a performance-driven dominant gating control flow based on the relative difference in their effectiveness: Trigger the strong blocking gate: Execute the conditional branch check if auc0 + 0.01 < auc1. If the condition holds, the system determines that the local sequence features extracted by the deep fully connected layer introduce serious overfitting noise in the data distribution of this fold, that is, the performance lags far behind the topological graph features. At this time, the system directly triggers the gating blocking mechanism, cuts off and discards the data stream output of the deep neural network branch; forcibly locks the system fusion weight of the heterogeneous graph autoencoder branch to 1.0, and the final cell communication prediction probability directly adopts the second prediction probability result set.

[0051] Adaptive weight grid search and self-optimizing fusion: If the above blocking conditions are not met, the system determines that both bimodal features have collaborative decision-making value under the current distribution. The system starts a one-dimensional grid search loop with a fixed step size (such as 0.01) within the preset dominant weight window : Simulate the execution of weighted feature synthesis in memory and call the ROC evaluation module to dynamically calculate the composite evaluation metric under the current evaluation set; the system automatically captures and locks the optimal fusion weight coefficient when the composite evaluation metric reaches the global maximum. Finally, use this optimal fusion coefficient to perform linear weighted composition on the rank proportion features of the two-branch data on the main data stream and output the final interaction prediction score. This control strategy ensures the absolute control of topological information in the fusion system on the one hand, and eliminates the suboptimal fusion performance loss caused by artificially fixed parameters on the other hand.

[0052] In the final full-scale prediction stage, after this method repeats five-fold cross-validation in the preset custom rounds (such as 20 rounds), the system screens the running round with the optimal average performance according to the average AUC of each five-fold, and saves the parameters of the deep neural network branch model and the heterogeneous graph autoencoder branch model with the highest single-fold AUC in this round. After the end of the large loop, the system solidifies and exports the parameters of the specific round with the most robust generalization performance, completing the storage and deployment of the highly reliable cell communication prediction model.

[0053] Step 5: Output of Interactive Prediction Results. Based on the final interactive prediction results output in Step 4, the system performs a scan based on the optimal dynamic evaluation index (or directly outputs a list of prediction scores), and outputs the interactive prediction scores or ranking results of candidate ligand-receptor pairs according to the final judgment label, for use in downstream intercellular communication pathway analysis or targeted drug screening.

[0054] The above description, in conjunction with specific preferred embodiments, provides a further detailed explanation of the present invention. It should not be construed that the specific implementation of the present invention is limited to these descriptions. For those skilled in the art, various simple deductions or substitutions can be made without departing from the concept of the present invention, and all such modifications and substitutions should be considered within the scope of protection of the present invention.

Claims

1. A method for predicting intercellular communication based on heterogeneous graph ensemble and dynamic fusion, characterized in that, Includes the following steps: Step 1: Data Acquisition, Stratification, and Standardization: We obtain the expression features and interaction labels of ligand-receptor pairs, construct a heterogeneous graph containing ligand nodes, receptor nodes and interaction edges, use outer layer hierarchical cross-validation to divide the training fold and test fold, and divide the validation set within the training fold. We also standardize the input features of the deep neural network branches based on the statistics of the training features. Step 2: Training deep neural network branches and generating the first prediction result: The deep neural network branch that trains the output predicts the logarithmic value is trained using a binary cross-entropy loss with sigmoid transformation, combined with weight decay, gradient pruning, learning rate scheduling based on validation set AUC, early stopping, and parameter rollback to obtain the first prediction probability result. Step 3: Training of the heterogeneous graph autoencoder ensemble branch and generation of the second prediction result: The heterogeneous graph autoencoder ensemble branch is trained, multiple sub-models are initialized with multiple random seeds, the evaluation set AUC is retained based on the current outer layer cross-validation for training monitoring and optimal round parameter rollback, and the edge prediction probabilities output by multiple heterogeneous graph autoencoder sub-models are averaged according to the sample dimension to obtain the second prediction probability result. Step 4: Dynamic fusion of ordinal rank normalization and the retained evaluation set AUC based on the current outer cross-validation: The first and second prediction probability results are sorted and rank-normalized. Dynamic gating is performed based on the current evaluation set AUC, or the fusion coefficient is searched within the weight interval where the heterogeneous graph autoencoder branch is dominant, to obtain the final interactive prediction score. Step 5: Output the interactive prediction results: Output the interaction prediction score or ranking result of candidate ligand-receptor pairs.

2. The method according to claim 1, characterized in that, The outer layer stratified cross-validation described in step 1 is implemented using a stratified K-fold cross-validation algorithm. The algorithm divides the data into layers according to the ratio of ligand receptor to interactive label categories, ensuring that the ratio of positive and negative samples in each training fold and test fold is consistent with the overall sample distribution. Within each training fold, the deep neural network branch is further divided into internal training and validation sets according to the principle of label distribution consistency. The outer layer stratified cross-validation has 5 folds, and the internal validation set is divided at a ratio of 15%.

3. The method according to claim 1, characterized in that, The standardization process in step 1 calculates the mean and standard deviation of each feature based only on the features of the internal training set, and uses the mean and standard deviation to perform zero-mean and unit variance transformations on the features of the internal training set, validation set, and outer test set, respectively. The validation set and test set are not included in the calculation of standardized statistics to avoid information leakage from the evaluation data into the training process.

4. The method according to claim 1, characterized in that, The deep neural network branch in step 2 is a multi-layer fully connected neural network. Its output layer does not have a Sigmoid activation function, but directly outputs the predicted log value. During training, a binary cross-entropy loss function that combines log odds is used, so that the Sigmoid transformation and the binary cross-entropy calculation are completed together inside the loss function. In the prediction stage, the Sigmoid transformation is performed on the predicted log value to obtain the first prediction probability result.

5. The method according to claim 4, characterized in that, The deep neural network branch in step 2 uses the Adam optimizer for parameter updates with a learning rate of 0.001 and an explicit weight decay parameter of [value missing]. After backpropagation of the loss function and before the optimizer updates the parameters, global norm pruning is performed on the gradients of the network parameters with a pruning threshold of 5.0; and dynamic scheduling of the learning rate is performed based on the validation set AUC, where the scheduling mode is to maximize AUC, the learning rate decay factor is 0.5, and the scheduling patience value is 15.

6. The method according to claim 4, characterized in that, In step 2, when the validation set AUC exceeds the historical best value, the model parameters of the current deep neural network branch are temporarily stored; when the validation set AUC fails to improve continuously and reaches a preset early stopping threshold, training is stopped, and the deep neural network branch is restored to the model parameters corresponding to the highest validation set AUC. The early stopping threshold is 40 rounds, and the AUC improvement judgment threshold is... .

7. The method according to claim 1, characterized in that, The heterogeneous graph autoencoder branch in step 3 includes a graph encoder and a bilinear decoder. The graph encoder uses GraphSAGE aggregation to learn the topological representation of the heterogeneous graph composed of ligand nodes, receptor nodes and their interaction edges. The bilinear decoder calculates the edge prediction probability based on the ligand node embedding and the receptor node embedding.

8. The method according to claim 7, characterized in that, In step 3, multiple heterogeneous graph autoencoder sub-models are initialized and trained independently using different and reproducible random seeds. The random seeds are determined by the base random seed, the current outer layer cross-validation fold index, and the current sub-model index. Each heterogeneous graph autoencoder sub-model sets a patience value mechanism based on the current outer layer cross-validation retention evaluation set AUC. When the retention evaluation set AUC does not improve for a preset number of consecutive rounds, the training is terminated, and the sub-model parameters are rolled back to the parameter state corresponding to the highest retention evaluation set AUC of that round. In addition, the number of heterogeneous graph autoencoder sub-models in step 3 is 3; after each sub-model completes training and rollback of the optimal round parameters, it outputs the corresponding edge prediction probability, and performs an arithmetic mean on the edge prediction probabilities output by multiple sub-models according to the sample dimension to obtain the second prediction probability result.

9. The method according to claim 1, characterized in that, The ranking rank normalization process in step 4 includes: performing stable ranking on the first prediction probability result and the second prediction probability result respectively, generating a rank value according to the ranking position of each sample in the corresponding result set, and mapping the rank value to the rank ratio value within a preset interval to obtain the first rank normalization result and the second rank normalization result, thereby reducing the impact of the difference in the output scale of the two branches on the fusion result.

10. The method according to claim 9, characterized in that, The dynamic fusion in step 4 includes dynamic gating and dynamic weight search: when the AUC of the deep neural network branch on the current evaluation set is lower than the AUC of the heterogeneous graph autoencoder branch on the current evaluation set and the difference exceeds a preset threshold, the output of the heterogeneous graph autoencoder branch is directly used as the final interactive prediction score; when dynamic gating is not triggered, the candidate fusion coefficients are traversed within a preset fusion weight range where the weight of the heterogeneous graph autoencoder branch is greater than the weight of the deep neural network branch, and the optimal fusion coefficient is determined with the maximum AUC of the current evaluation set as the optimization objective. The first-rank normalization result and the second-rank normalization result are linearly weighted and fused according to the optimal fusion coefficient to obtain the final interactive prediction score; the weight of the heterogeneous graph autoencoder branch is searched between 0.80 and 1.00 with a search step size of 0.01.

Images