Method and system for predicting protein drug binding sites based on multi-modal dynamic graph

Abstract

The application discloses a method and system for predicting protein drug binding sites based on a multi-modal dynamic graph, comprising: obtaining amino acid sequence data and three-dimensional structure data of a protein, and generating evolutionary conservation features, structure graph topology features and sequence features respectively; inputting the features into a multi-modal fusion encoder to generate first fusion features; inputting the first fusion features into a prediction decoder to generate initial prediction probabilities; iteratively updating the structure graph topology features according to a preset three-graph update rule according to the initial prediction probabilities and residue dynamic communication scores; re-inputting the updated structure graph topology features into the multi-modal fusion encoder and the prediction decoder, and repeatedly executing until a preset iteration number is reached, to generate final prediction probabilities; and generating a residue importance heat map and outputting a prediction report based on the final prediction probabilities through a gradient weighted class activation mapping algorithm. The application realizes the cooperative optimization of prediction and graph structure, and significantly improves the accuracy of binding site prediction.

Description

Technical Field

[0001] This invention relates to the field of bioinformatics, specifically to a method and system for predicting protein drug binding sites based on multimodal dynamic graphs. Background Technology

[0002] Accurate prediction of protein drug binding sites is a core component of computer-aided drug design. Currently, deep learning methods based on graph neural networks have become mainstream, such as constructing protein residue contact maps to extract structural topological features and combining them with sequence information for prediction. However, the protein graph structures constructed by existing methods are usually one-time static graphs based on fixed spatial distance thresholds, which remain unchanged throughout the prediction process. This static graph structure cannot be adaptively adjusted based on preliminary prediction results and is difficult to effectively utilize the inherent dynamic characteristics of proteins, such as functional communication and cooperative movement between residues, resulting in insufficient identification of key binding sites, especially allosteric sites. Furthermore, existing methods often employ fixed weight strategies when fusing multimodal information, failing to adaptively adjust the contribution of each modality feature according to different protein folding types, further limiting the generalization ability and accuracy of the prediction model. Summary of the Invention

[0003] In view of the above problems, the present invention provides a method and system for predicting protein drug binding sites based on multimodal dynamic graphs. By enabling the graph structure to iteratively evolve according to the prediction results and the dynamic communication information of the protein, the prediction and the graph structure are synergistically optimized.

[0004] To achieve the above objectives, in a first aspect, this application provides a method for predicting protein drug binding sites based on multimodal dynamic graphs, comprising:

[0005] Obtain the amino acid sequence data and three-dimensional structure data of the target protein;

[0006] The amino acid sequence data is used to construct a multiple sequence alignment file using a multiple sequence alignment tool, and the site-specific scoring matrix and residue conservation entropy value are calculated to generate evolutionary conservation features.

[0007] The three-dimensional structural data are used to calculate the dynamic communication score of residues through an elastic network model, extract the spatial contact relationship of residues, and generate the topological features of the structural diagram.

[0008] Amino acid sequence data is input into a protein language model for encoding to generate sequence features;

[0009] Evolutionary conservation features, structural graph topological features, and sequence features are input into a pre-trained multimodal fusion encoder to generate the first fusion feature;

[0010] The first fused feature is input into a pre-trained prediction decoder to generate the initial binding site prediction probability;

[0011] The initial binding site prediction probability and residue dynamic communication score are used to iteratively update the topological features of the structural graph according to the three preset graph update rules to generate the updated topological features of the structural graph. The three graph update rules include the high confidence residue edge enhancement rule, the low confidence far residue edge decay rule, and the high dynamic communication score residue edge addition rule.

[0012] The updated topological features of the structure graph are re-input into the multimodal fusion encoder and predictive decoder, and the feature fusion and prediction steps are repeated until the preset number of iterations is reached to generate the final binding site prediction probability.

[0013] Based on the predicted probability of the final binding site, the contribution weight of each residue to the prediction result is calculated by the gradient weighted class activation mapping algorithm, a residue importance heatmap is generated, and a prediction report containing the predicted probability of the final binding site and the residue importance heatmap is output.

[0014] Furthermore, the amino acid sequence data is used to construct a multiple sequence alignment file using a multiple sequence alignment tool, and the site-specific scoring matrix and residue conservation entropy value are calculated to generate evolutionary conservation features, including:

[0015] Obtain the amino acid sequence data of the target protein, input the amino acid sequence data into the HHblits tool, use the UniRef30 database as the search target, perform a multiple sequence alignment search, and generate a multiple sequence alignment file;

[0016] Read the multiple sequence alignment file, count the frequency of amino acid types appearing at each residue position during the evolution process, use a pseudo-counting smoothing strategy to smooth the frequency distribution, and generate a site-specific scoring matrix;

[0017] Based on the site-specific scoring matrix, the information entropy of each residue position is calculated using the Shannon entropy formula to generate the residue conservation entropy value;

[0018] The site-specific scoring matrix is ​​concatenated with the residue conservation entropy value to form the initial evolutionary feature matrix;

[0019] The initial evolutionary feature matrix is ​​input into a pre-trained Transformer encoder, which contains a multi-layer multi-head self-attention mechanism and a feedforward neural network. The multi-head self-attention mechanism captures the long-range evolutionary dependencies between residues in the initial evolutionary feature matrix and generates evolutionarily conserved features.

[0020] Furthermore, the three-dimensional structural data is used to calculate the dynamic communication score of residues through an elastic network model, including:

[0021] Obtain the three-dimensional structural data of the target protein, and extract the C-value of each residue from the three-dimensional structural data. β Atomic coordinates, for glycine residues, C is extracted. α Atomic coordinates are used as a substitute;

[0022] C based on all residues α Or C β Atomic coordinates are used to calculate the spatial Euclidean distance between residue pairs. A distance truncation threshold is set. When the spatial Euclidean distance between residue pairs is less than the distance truncation threshold, a spring connection is constructed between the residue pairs to generate the Kirchhoff matrix of the elastic network model.

[0023] Perform a pseudo-inverse operation on the Kirchhoff matrix to generate a pseudo-inverse matrix. The diagonal elements of the pseudo-inverse matrix represent the theoretical mean square fluctuations of each residue under the elastic network model, and the off-diagonal elements of the pseudo-inverse matrix represent the fluctuation covariance between different residue pairs.

[0024] Based on the off-diagonal and diagonal elements of the pseudo-inverse matrix, the Pearson correlation coefficient between each pair of residues is calculated using the normalized covariance formula, generating a dynamic correlation matrix of residues.

[0025] The residue dynamic correlation matrix is ​​row normalized, and the average value of the summation of the dynamic correlation coefficients of each residue with all other residues is taken to generate the residue dynamic communication score of each residue. The residue dynamic communication score is used to quantify the degree of communication hub of each residue in the global dynamic network of the protein.

[0026] Furthermore, the spatial contact relationships of residues are extracted to generate structural graph topological features, including:

[0027] Obtain the three-dimensional structural data of the target protein, and extract the C-value of each residue from the three-dimensional structural data. β Atomic coordinates, for glycine residues, C is extracted. α Atomic coordinates are used as a substitute;

[0028] C based on all residues α Atomic coordinates or C β Atomic coordinates are used to calculate the spatial Euclidean distance between residue pairs. A spatial contact distance threshold is set. When the spatial Euclidean distance between residue pairs is less than the spatial contact distance threshold, it is determined that there is a spatial contact relationship between the residue pairs, and a set of spatial contact edges is generated.

[0029] Obtain the sequence position index of residues in the amino acid sequence data. When the absolute value of the difference between the sequence position indices of two residues is a preset value, it is determined that there is a sequence adjacency relationship between the residue pairs, and a set of sequence adjacency edges is generated.

[0030] The set of spatial contact edges and the set of sequential adjacent edges are merged to generate a heterogeneous edge set, which contains two edge types: spatial contact edges and sequential adjacent edges.

[0031] Construct an edge feature vector for each edge in the heterogeneous edge set. The edge feature vector shall contain at least the normalized distance, the direction cosine component, and the edge type one-hot encoding.

[0032] An initial heterogeneous graph structure is constructed using residues as nodes, heterogeneous edge sets as edges, and edge feature vectors as edge attributes.

[0033] Obtain the Laplacian matrix of the initial heterogeneous graph structure, perform eigenvalue decomposition on the Laplacian matrix, extract the eigenvectors corresponding to the minimum non-zero eigenvalues ​​of the first preset number, generate the Laplacian position code, and attach the Laplacian position code as the node position feature to each node of the initial heterogeneous graph structure to generate the topological features of the structure graph.

[0034] Furthermore, the amino acid sequence data is input into a protein language model for encoding to generate sequence features, including:

[0035] The amino acid sequence data of the target protein is obtained, and the amino acid sequence data is standardized by replacing non-standard amino acids with corresponding standard amino acids to generate a standardized amino acid sequence.

[0036] The standardized amino acid sequence is input into the ESM-2 protein language model, which contains a multi-layer Transformer encoder. The network parameters of the first preset number of Transformer encoders of the ESM-2 protein language model are frozen, and the network parameters of the last preset number of Transformer encoders are fine-tuned to generate the ESM-2 hidden layer representation.

[0037] The ESM-2 hidden layer representation is input into a multi-scale one-dimensional convolution module, which contains multiple parallel one-dimensional convolution kernels, each with a different kernel size. Multi-scale local pattern extraction is performed on the ESM-2 hidden layer representation through multiple parallel one-dimensional convolution kernels to generate multi-scale convolution features.

[0038] The convolutional features at each scale in the multi-scale convolutional features are concatenated to generate concatenated convolutional features.

[0039] The concatenated convolutional features are input into a layer normalization layer and a fully connected layer for dimensionality transformation to generate sequence features.

[0040] Furthermore, the evolutionary conservation features, structural graph topological features, and sequence features are input into a pre-trained multimodal fusion encoder to generate the first fusion feature, including:

[0041] Evolutionary conservation features, structural graph topology features, and sequence features are obtained. These features are then input into their respective self-attention modules. Each self-attention module contains a multi-head self-attention mechanism. The multi-head self-attention mechanism is used to strengthen the long-range dependencies between residues within each modality feature, thereby generating self-attention refined evolutionary conservation features, self-attention refined structural graph topology features, and self-attention refined sequence features.

[0042] Obtain the SCOP fold type classification label of the target protein. If the SCOP fold type classification label cannot be obtained, use the preset general prior gating weight vector or predict its most likely fold type through an auxiliary classifier to generate the corresponding prior gating weight vector.

[0043] Based on the SCOP fold type classification label, the corresponding prior gating weight vector is found from the preset SCOP prior gating matrix. The rows of the SCOP prior gating matrix correspond to different SCOP fold types, and the columns of the SCOP prior gating matrix correspond to three modes: evolutionary conservation features, structural graph topology features, and sequence features.

[0044] The evolutionary conservatism features refined by self-attention, the topological features of the structure graph refined by self-attention, and the sequence features refined by self-attention are used as key vectors, and the globally learnable query vector is used as the query vector. The dot product attention score of the key vector and the query vector is calculated to generate the initial dynamic weights.

[0045] The initial dynamic weights are multiplied element-wise with the prior gated weight vector to generate the modulated dynamic weights. The modulated dynamic weights are then normalized using softmax to generate normalized dynamic weights.

[0046] Based on normalized dynamic weights, the evolutionary conservation features refined by self-attention, the topological features of the structure graph refined by self-attention, and the sequence features refined by self-attention are weighted and summed to generate residue-level fusion features.

[0047] Using residue-level fusion features as query vectors, the evolutionary conservation features refined by self-attention, the topological features of the structure graph refined by self-attention, and the sequence features refined by self-attention are concatenated as key-value pairs, and multi-head cross-attention operation is performed to generate cross-attention refined features.

[0048] The cross-attention refined features are input into the feedforward neural network for nonlinear transformation to generate the first fused features.

[0049] Furthermore, the first fused feature is input into a pre-trained prediction decoder to generate initial binding site prediction probabilities, including:

[0050] The first fusion feature is input into a bidirectional gated recurrent unit, which includes a forward-gated recurrent unit and a backward-gated recurrent unit. The forward-gated recurrent unit captures the contextual dependencies of residues along the amino acid sequence from the N-terminus to the C-terminus, and the backward-gated recurrent unit captures the contextual dependencies of residues along the amino acid sequence from the C-terminus to the N-terminus. The hidden outputs of the forward-gated recurrent unit and the backward-gated recurrent unit are spliced ​​together to generate sequence context features.

[0051] The sequence context features are input into a multilayer perceptron, which contains multiple sequentially connected fully connected layers and activation function layers. The sequence context features are then subjected to a layer-by-layer nonlinear transformation through the multiple sequentially connected fully connected layers and activation function layers to generate log odds values.

[0052] The log-odds value is input into the sigmoid activation function, which maps the log-odds value to a probability range between 0 and 1, generating the initial predicted probability of each residue belonging to the ortho-binding site.

[0053] Furthermore, the initial predicted binding site probabilities and residue dynamic communication scores are iteratively updated according to three preset graph update rules to generate updated structural graph topological features, including:

[0054] The predicted probability of the initial binding site, the dynamic communication score of the residue, and the topological features of the structural graph are obtained. The topological features of the structural graph include the set of nodes, the set of edges, and the edge weight matrix.

[0055] Set a first confidence threshold and a second confidence threshold. If the first confidence threshold is higher than the second confidence threshold, set a far-end distance threshold.

[0056] Traverse each edge in the edge set and obtain the initial binding site prediction probability of each of the two residue nodes connected by each edge and the spatial Euclidean distance between the two residue nodes;

[0057] When the initial binding site prediction probabilities of both residue nodes are higher than the first confidence threshold, the edge weights corresponding to the edges are multiplied by the enhancement coefficient to generate updated edge weights. The enhancement coefficient is positively correlated with the product of the initial binding site prediction probabilities of the two residue nodes.

[0058] When the initial binding site prediction probabilities of two residue nodes are both lower than the second confidence threshold and the spatial Euclidean distance between the two residue nodes is greater than the far-end distance threshold, the edge weight corresponding to the edge is multiplied by the decay coefficient to generate the updated edge weight. The decay coefficient is a positive number less than 1.

[0059] Traverse the dynamic communication scores of residues, obtain the dynamic correlation coefficient between each pair of residue nodes, set a communication score threshold, and when the dynamic correlation coefficient between a pair of residue nodes is higher than the communication score threshold and there is no edge connecting a pair of residue nodes in the edge set, add an edge to the edge set and assign an initial edge weight to the added edge.

[0060] The edge set and edge weight matrix after edge weight update and edge addition operations are used to replace the original edge set and original edge weight matrix in the topological features of the structure graph, thereby generating the updated topological features of the structure graph.

[0061] Furthermore, based on the predicted probability of the final binding site, the contribution weight of each residue to the prediction result is calculated using a gradient-weighted class activation mapping algorithm, generating a residue importance heatmap, including:

[0062] Obtain the final binding site prediction probability and the first fusion feature, which contains the fusion feature vector corresponding to each residue;

[0063] The predicted score corresponding to the target category in the final binding site prediction probability is used as the target score. The target category is the ortho-binding site category or the allosteric binding site category. Backpropagation is performed on the target score to calculate the gradient value of the target score relative to the fusion feature vector corresponding to each residue in the first fusion feature.

[0064] Global average pooling is performed on the gradient values ​​along the feature channel dimension to generate the neuron importance weights corresponding to each residue. The neuron importance weights represent the degree of contribution of each feature channel to the target score.

[0065] The fusion feature vector corresponding to each residue is weighted and summed with the neuron importance weight to generate the initial importance score for each residue;

[0066] The initial importance score is processed using the ReLU activation function, retaining the positive value part and suppressing the negative value part to generate the importance score after activation;

[0067] The importance score after activation is subjected to min-max normalization, which maps the importance score after activation to the interval between 0 and 1, generating a normalized residue importance score for each residue.

[0068] The normalized residue importance score is written as a temperature factor into the protein's three-dimensional structure file to generate a residue importance heatmap.

[0069] In a second aspect, the present invention also provides a system for predicting protein drug binding sites based on a multimodal dynamic graph, applicable to the method described in the first aspect. The system includes a data acquisition module, an evolutionary feature extraction module, a structural feature extraction module, a sequence feature extraction module, a multimodal fusion encoder, a prediction decoder, a graph topology iterative update module, an iteration control module, and an interpretable output module. The data acquisition module acquires the amino acid sequence data and three-dimensional structural data of the target protein. The evolutionary feature extraction module constructs a multi-sequence alignment file from the amino acid sequence data using a multi-sequence alignment tool, calculates the site-specific scoring matrix and residue conservation entropy, and generates evolutionarily conserved features. The structural feature extraction module calculates the residue dynamic communication score from the three-dimensional structural data using an elastic network model, extracts the spatial contact relationships of residues, and generates structural graph topology features. The sequence feature extraction module inputs the amino acid sequence data into a protein language model for encoding, generating sequence features. The multimodal fusion encoder integrates the evolutionarily conserved features, structural graph topology features, and sequence features. Adaptive weighted fusion generates the first fused feature; the predictive decoder performs sequence context modeling and nonlinear transformation on the first fused feature to generate the initial binding site prediction probability; the graph topology iterative update module iteratively updates the structural graph topology feature according to three preset graph update rules, including a high-confidence residue edge enhancement rule, a low-confidence far-end residue edge decay rule, and a high-dynamic-communication-score residue edge addition rule; the iterative control module re-inputs the updated structural graph topology feature into the multimodal fusion encoder and predictive decoder, repeating the feature fusion and prediction steps until a preset number of iterations is reached to generate the final binding site prediction probability; the interpretable output module calculates the contribution weight of each residue to the prediction result based on the final binding site prediction probability using a gradient-weighted class activation mapping algorithm, generates a residue importance heatmap, and outputs a prediction report containing the final binding site prediction probability and the residue importance heatmap.

[0070] Unlike existing technologies, the above technical solution provides a method and system for predicting protein drug binding sites based on multimodal dynamic graphs. The method includes: acquiring amino acid sequence data and three-dimensional structural data of the protein, generating evolutionarily conserved features, structural graph topological features, and sequence features, respectively; inputting these features into a multimodal fusion encoder to generate a first fusion feature; inputting the first fusion feature into a prediction decoder to generate an initial prediction probability; iteratively updating the structural graph topological features according to three preset graph update rules using the initial prediction probability and residue dynamic communication scores; re-inputting the updated structural graph topological features into the multimodal fusion encoder and prediction decoder, repeating this process until a preset number of iterations is reached to generate a final prediction probability; and generating a residue importance heatmap based on the final prediction probability using a gradient-weighted class activation mapping algorithm and outputting a prediction report. This invention achieves synergistic optimization of prediction and graph structure, significantly improving the accuracy of binding site prediction.

[0071] The above description of the invention is merely an overview of the technical solution of this application. In order to enable those skilled in the art to better understand the technical solution of this application and to implement it based on the description and drawings, and to make the above-mentioned objectives and other objectives, features and advantages of this application easier to understand, the following description is provided in conjunction with the specific embodiments and drawings of this application. Attached Figure Description

[0072] The accompanying drawings are only used to illustrate the principles, implementation methods, applications, features, and effects of specific embodiments of the present invention and other related contents, and should not be considered as limitations on this application.

[0073] In the accompanying drawings of the instruction manual:

[0074] Figure 1 This is a schematic diagram of the overall process of the method described in the specific implementation;

[0075] Figure 2 This is a schematic diagram illustrating steps S101 to S105 of the method described in the specific implementation embodiment;

[0076] Figure 3 This is a schematic diagram illustrating steps S201 to S207 of the method described in a specific implementation.

[0077] Figure 4 This is a schematic diagram illustrating steps S301 to S306 of the method described in a specific implementation.

[0078] Figure 5 This is a schematic diagram of the system described in a specific implementation. Detailed Implementation

[0079] To illustrate the possible application scenarios, technical principles, implementable specific solutions, and achievable objectives and effects of this application in detail, the following description, in conjunction with the listed specific embodiments and accompanying drawings, provides a detailed explanation. The embodiments described herein are merely illustrative of the technical solutions of this application and are therefore intended to limit the scope of protection of this application.

[0080] In this document, the term "embodiment" means that a specific feature, structure, or characteristic described in connection with an embodiment may be included in at least one embodiment of this application. The term "embodiment" appearing in various places throughout the specification does not necessarily refer to the same embodiment, nor does it specifically limit its independence or connection with other embodiments. In principle, in this application, as long as there are no technical contradictions or conflicts, the technical features mentioned in each embodiment can be combined in any way to form corresponding implementable technical solutions.

[0081] Unless otherwise defined, the technical terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this application pertains; the use of related terms herein is merely for the purpose of describing particular embodiments and is not intended to limit this application.

[0082] In the description of this application, the term "and / or" is used to describe the logical relationship between objects, indicating that three relationships can exist. For example, A and / or B means: A exists, B exists, and A and B exist simultaneously. Additionally, the character " / " in this document generally indicates that the preceding and following objects have an "or" logical relationship.

[0083] In this application, terms such as “first” and “second” are used only to distinguish one entity or operation from another, and do not necessarily require or imply any actual quantity, hierarchy or order relationship between these entities or operations.

[0084] Without further limitations, the use of terms such as “comprising,” “including,” “having,” or other similar open-ended expressions in this application is intended to cover non-exclusive inclusion, which does not exclude the presence of additional elements in a process, method, or product that includes the stated elements, such that a process, method, or product that includes a list of elements may include not only those defined elements but also other elements not expressly listed, or elements inherent to such a process, method, or product.

[0085] As understood in the Examination Guidelines, in this application, expressions such as "greater than," "less than," and "exceeding" are understood to exclude the stated number; expressions such as "above," "below," and "within" are understood to include the stated number. Furthermore, in the description of the embodiments in this application, "multiple" means two or more (including two), and similar expressions related to "multiple" are also understood in this way, such as "multiple groups" and "multiple times," unless otherwise explicitly specified.

[0086] The processor described in the embodiments of this application can be implemented by hardware, firmware, software, or a combination thereof. It can be a circuit, one or more of an application-specific integrated circuit (ASIC), a digital signal processor (DSP), a digital signal processing device (DSPD), a programmable logic device (PLD), a field-programmable gate array (FPGA), a central processing unit (CPU), a controller, a microcontroller, or a microprocessor. It also includes other physical, biological, or chemical structures that can implement the same or equivalent functions as the processors listed above, such as biological neurons, quantum computing units, DNA computing units, etc., so that the processor can execute some or all of the steps in the computer program or method involved in the various embodiments of this application, or any combination of the steps mentioned therein.

[0087] The computer program involved in the embodiments can be stored in a computer device readable storage medium, which includes, but is not limited to, disks, magnetic tapes, magnetic cards, floppy disks, flash memory, optical disks, optical cards, read-only memory (ROM), random access memory (RAM), erasable programmable ROM (EPROM), and electrically erasable programmable ROM (EEPROM), etc., and also includes other biological, physical, or chemical structures that can achieve the same or equivalent functions as the storage media listed above, such as DNA, RNA, proteins, and other units with information storage capabilities. In specific embodiments, the storage medium involved can be one of the above-mentioned media types, or a combination of the above-mentioned media types. In different embodiments, the computer program involved in the embodiments can be centrally stored in a single medium, or distributed and stored in multiple media. The memory containing the computer device readable storage medium can be non-volatile memory or random access memory. These computer device readable storage media can be built into the device, or can be connected to the device involved in the embodiments as an external device or part of an external device. In some embodiments, the memory having a computer device readable storage medium is deployed locally; in other embodiments, the memory may be deployed remotely from the processor, for example, as a network-attached memory accessed via RF circuitry or an external port and a communication network, wherein the communication network may be the Internet, one or more intranets, a local area network (LAN), a wide area network (WLAN), a storage area network (SAN), or a suitable combination thereof, as long as computer device access to the memory is enabled. Furthermore, the computer program involved in the embodiments may be stored in plaintext / ciphertext form, or it may be designed as training data, integrated and recombined through model training and implicitly stored in the parameter states of a deep neural network or other machine learning model.

[0088] Please see Figure 1 In a first aspect, this embodiment provides a method for predicting protein drug binding sites based on multimodal dynamic graphs, comprising:

[0089] Obtain the amino acid sequence data and three-dimensional structure data of the target protein;

[0090] The amino acid sequence data is used to construct a multiple sequence alignment file using a multiple sequence alignment tool, and the site-specific scoring matrix and residue conservation entropy value are calculated to generate evolutionary conservation features.

[0091] The three-dimensional structural data are used to calculate the dynamic communication score of residues through an elastic network model, extract the spatial contact relationship of residues, and generate the topological features of the structural diagram.

[0092] Amino acid sequence data is input into a protein language model for encoding to generate sequence features;

[0093] Evolutionary conservatism features, structural graph topology features, and sequence features are input into a pre-trained multimodal fusion encoder to generate the first fusion feature. The multimodal fusion encoder includes a gated attention fusion module with SCOP fold type as a priori. The gated attention fusion module performs adaptive weighted fusion of evolutionary conservatism features, structural graph topology features, and sequence features.

[0094] The first fused feature is input into a pre-trained prediction decoder to generate the initial binding site prediction probability;

[0095] The initial binding site prediction probability and residue dynamic communication score are used to iteratively update the topological features of the structural graph according to the three preset graph update rules to generate the updated topological features of the structural graph. The three graph update rules include the high confidence residue edge enhancement rule, the low confidence far residue edge decay rule, and the high dynamic communication score residue edge addition rule.

[0096] The updated topological features of the structure graph are re-input into the multimodal fusion encoder and predictive decoder, and the feature fusion and prediction steps are repeated until the preset number of iterations is reached to generate the final binding site prediction probability.

[0097] Based on the predicted probability of the final binding site, the contribution weight of each residue to the prediction result is calculated by the gradient weighted class activation mapping algorithm, a residue importance heatmap is generated, and a prediction report containing the predicted probability of the final binding site and the residue importance heatmap is output.

[0098] In this embodiment, the amino acid sequence data is derived from the primary structure information of the target protein, typically recorded in FASTA format as a linear sequence of amino acids; the three-dimensional structure data describes the spatial conformation of the protein, storing the three-dimensional coordinates of each atom in PDB format. This data can be obtained through searches of public databases or generated experimentally or using structure prediction tools, serving as the initial input for subsequent multimodal feature extraction.

[0099] Multiple sequence alignment (MSA) files are constructed using a multiple sequence alignment tool. A homology sequence search algorithm is then used to retrieve a set of homologous sequences with evolutionary relationships to the target sequence from a protein sequence database. These sequences are then arranged in a multi-alignment manner according to their evolutionary relationships. A site-specific scoring matrix is ​​calculated from the frequency of amino acids at each residue position in the MSA file after pseudo-counting smoothing, characterizing the preference of each residue position for different amino acid types. The residue conservation entropy value is calculated based on the aforementioned frequency distribution using the information entropy formula; a lower entropy value indicates greater conservation of the site during evolution. The two matrices are then concatenated and further extracted by an encoder to form evolutionary conservation features reflecting the evolutionary conservation pattern of the protein.

[0100] The elastic network model abstracts the three-dimensional structure of proteins into C-type structures.α Or C β An elastic network with atoms as nodes and spring connections between spatially adjacent residues as edges is used to calculate the degree of mutual influence of each residue in the global dynamic motion of the protein by analyzing the vibration modes of this network. Preferably, the residue dynamic communication score is obtained by normalizing the covariance after pseudo-inverse operation of the Kirchhoff matrix of the elastic network model, and is used to quantify the degree of communication hub of each residue in the protein dynamic network; the spatial contact relationship of residues is extracted by setting a spatial distance threshold to determine whether residue pairs are adjacent in three-dimensional space. The topological features of the structural graph construct a heterogeneous graph structure with residues as nodes and spatial contact relationships and sequence adjacency relationships as edges, and add Laplacian position encoding as node position features to comprehensively characterize the spatial topology and dynamic coupling information of the protein.

[0101] Protein language models learn the semantic representation of amino acids and sequence context dependency patterns through self-supervised pre-training on large-scale protein sequence data. When amino acid sequence data is input into the model, each residue is mapped to a high-dimensional continuous vector representation. The set of these vectors is the sequence feature, which is used to capture long-range semantic dependencies between residues from the sequence dimension.

[0102] In the multimodal fusion encoder, SCOP fold types are classified according to the three-dimensional folding patterns of protein domains. Proteins with different fold types have fundamentally different dependence patterns on evolutionary, structural, and sequence information. The gated attention fusion module introduces this prior information to dynamically adjust the contribution weights of the three modal features in the fusion process, generating a first fusion feature that integrates complementary evolutionary, structural, and sequence information.

[0103] The predictive decoder is used to perform sequence context modeling and nonlinear transformation on the first fusion feature, and outputs a probability value between 0 and 1 for each residue. This probability value represents the confidence that the residue belongs to the ortho-binding site, which is the initial binding site prediction probability.

[0104] The three graph update rules correspond to different biological motivations: the high-confidence residue pair enhancement rule targets residue pairs with high confidence in the initial prediction, amplifying the synergistic effect of high-confidence regions by increasing their connection edge weights; the low-confidence distant residue pair attenuation rule targets residue pairs with low confidence and large spatial distances, suppressing noise interference by reducing their connection edge weights; and the high-dynamic-communication-score residue pair addition rule targets residue pairs with high communication strength in the elastic network model but lacking direct connections in the current graph structure, capturing non-spatial functional coupling relationships by adding edges. The updated graph topology features are the graph structure optimized by the above rules.

[0105] The updated topological features of the structure map are re-input into the multimodal fusion encoder and predictor decoder, and the feature fusion and prediction steps are repeated. The preset number of iterations determines the number of times the prediction-feedback loop is executed. Through multiple iterations, the prediction results gradually converge to a stable state. Finally, the predicted probability of the binding site is the binding probability of each residue output after the preset number of iterations.

[0106] The gradient-weighted activation mapping algorithm calculates the gradient of the final predicted probability relative to the feature vector corresponding to each residue in the fused features through backpropagation. The gradient value reflects the contribution of each residue to the prediction result. The residue importance heatmap is generated by normalizing the above contribution weights, which can intuitively show the importance distribution of each residue in the binding site prediction. The prediction report integrates the final binding site prediction probability and the residue importance heatmap in a structured format, providing drug development practice with prediction results that are both accurate and interpretable.

[0107] This embodiment achieves high-precision prediction of protein drug binding sites by constructing a complete process from multimodal data acquisition to interpretable predictive output. It employs multiple sequence alignment technology to extract evolutionarily conserved features, elastic network models to extract structural dynamic features, and protein language models to extract sequence semantic features, comprehensively characterizing the biological properties of proteins from three complementary dimensions. An adaptive fusion mechanism with SCOP fold type prior gating enables the integration of protein type-aware multimodal information. A prediction-feedback dynamic graph topology iterative update mechanism is introduced, allowing the graph structure to be adaptively refined based on prediction results, forming a closed-loop optimized prediction framework. Combined with a gradient-weighted class activation mapping algorithm, it provides residue-level interpretable output, offering computational support with both accuracy and reliability for drug development practice.

[0108] Combination Figure 1 and Figure 5 As shown, the raw protein data (FASTA sequence + PDB structure file) consists of the amino acid sequence and three-dimensional structure data of the target protein. Evolutionary conservation features, structural topology features, and sequence features are obtained through four-modal parallel feature extraction (M-FEM). Optimal extraction can also extract physicochemical MLP (additional physicochemical property features). ESM-2 corresponds to the protein language model, MT-GAT heterogeneity graph corresponds to the extraction of structural topology features, and evolutionary Transformer corresponds to the extraction of evolutionary conservation features. The output dimension of each branch is uniformly N×256 (N is the number of residues, 256 is the feature dimension). The first fused feature is obtained through H-CMAF:SCOP prior-gated hierarchical fusion. SCOP gating weights correspond to a gating attention fusion mechanism based on SCOP folding type as a priori. L1 / L2 / L3 three-level fusion corresponds to hierarchical fusion at three granularities: residue level, domain level, and global protein level. Cross-attention refinement is a cross-modal cross-attention operation. TSL: thermodynamic stability-aware feature modulation is a further optimized scheme, using the three-source residue flexibility index (RFI) to modulate the amplitude of the fused features residue by residue, and outputting the modulated features. ,in The modulation coefficient, The weight matrix is ​​a learnable matrix; IPGR: Prediction-Feedback Dynamic Iteration iteratively updates the topological features of the structure graph, and re-inputs the updated topological features into the multimodal fusion encoder and prediction decoder, repeating the feature fusion and prediction steps until the preset number of iterations (T=3) is reached, generating the final binding site prediction probability, where R1 corresponds to the high-confidence residue edge enhancement rule, R2 corresponds to the low-confidence far-end residue edge decay rule, and R3 corresponds to the high-dynamic communication score residue edge addition rule. For the first The median probability of the round prediction, The updated graph structure is shown below. Initial and final binding site prediction probabilities are generated using DH-PD (Hyper-Ortho- / Hyper-Ortho ... The predicted probability is for allosteric sites; the final IOM is: Grad-CAM multidimensional interpretable output prediction report, Grad-CAM heatmap corresponding to residue importance heatmap, B factor format is the temperature factor column of the protein three-dimensional structure file with normalized residue importance scores written, five-dimensional confidence CC is an optional multi-angle confidence assessment, and JSON structured report is the structured output format of the prediction results.

[0109] In some embodiments, amino acid sequence data are used to construct a multiple sequence alignment file using a multiple sequence alignment tool, and a site-specific scoring matrix and residue conservation entropy value are calculated to generate evolutionary conservation features, including:

[0110] Obtain the amino acid sequence data of the target protein, input the amino acid sequence data into the HHblits tool, use the UniRef30 database as the search target, perform a multiple sequence alignment search, and generate a multiple sequence alignment file;

[0111] Read the multiple sequence alignment file, count the frequency of amino acid types appearing at each residue position during the evolution process, use a pseudo-counting smoothing strategy to smooth the frequency distribution, and generate a site-specific scoring matrix;

[0112] Based on the site-specific scoring matrix, the information entropy of each residue position is calculated using the Shannon entropy formula to generate the residue conservation entropy value;

[0113] The site-specific scoring matrix is ​​concatenated with the residue conservation entropy value to form the initial evolutionary feature matrix;

[0114] The initial evolutionary feature matrix is ​​input into a pre-trained Transformer encoder, which contains a multi-layer multi-head self-attention mechanism and a feedforward neural network. The multi-head self-attention mechanism captures the long-range evolutionary dependencies between residues in the initial evolutionary feature matrix and generates evolutionarily conserved features.

[0115] In this embodiment, HHblits is a homology sequence retrieval tool based on iterative search using a Hidden Markov Model (HMM). It transforms the amino acid sequence data of the target protein into an HMM profile and performs multiple rounds of iterative searches in a protein sequence database. Each iteration updates the HMM profile using matching sequences obtained from the previous round, thereby progressively improving the sensitivity and depth of the search. The UniRef30 database is a non-redundant protein sequence database that aggregates protein sequences from sources such as UniProtKB and performs clustering and merging, significantly reducing data redundancy while ensuring sequence diversity and providing a broad set of reference sequences for homology searches. Multiple sequence alignment (MSE) searches use the HMM profile of the target protein as the query, retrieving sequences with evolutionary homology relationships from the UniRef30 database. The retrieved sequences are then arranged according to their evolutionary relationship with the target sequence, generating a MSE file recording the amino acid variation patterns at each residue position.

[0116] The proportion of each of the twenty standard amino acids appearing in each residue column of the multiple sequence alignment file is calculated relative to the total number of sequences. A pseudo-counting smoothing strategy is used to introduce a small prior count when calculating frequencies, avoiding the problem of some amino acid types having zero frequencies due to limited sample size, thus making the frequency distribution more robust and reliable. The site-specific scoring matrix is ​​the log-probability score matrix of the amino acid at each residue position after pseudo-counting smoothing, reflecting the preference for different amino acid types at each position and the strength of evolutionary constraints.

[0117] Based on the site-specific scoring matrix, the information entropy of each residue position is calculated using the Shannon entropy formula. Specifically, the logarithm of the probability of occurrence of each of the twenty amino acids at each residue position is multiplied by that probability, the sums are taken, and the negative number is taken. The resulting entropy value is used to quantify the degree of disorder in the amino acid distribution at that site. The residue conservation entropy value is the information entropy sequence calculated above. The lower the entropy value, the more concentrated and conserved the amino acid type distribution at that residue position is during evolution; the higher the entropy value, the more dispersed and variable the amino acid type distribution is at that position.

[0118] Each row of the site-specific scoring matrix (corresponding to a residue position) and the corresponding residue conservation entropy value are merged along the feature dimension to form an initial evolutionary feature matrix in which each row contains twenty amino acid preference scores and one conservation entropy value, which serves as the input representation for subsequent deep learning models.

[0119] The pre-trained Transformer encoder incorporates a multi-layered multi-head self-attention mechanism and a feedforward neural network. During pre-training, a large amount of protein multiple sequence alignment data is used as training samples. A masking modeling approach is employed as the pre-training task, randomly masking a portion of the initial evolutionary feature matrix. The encoder is required to predict the values ​​at the masked positions based on contextual information, learning statistical regularities and residue dependencies in evolutionary features by minimizing the difference between the predicted and true values. After pre-training, the multi-head self-attention mechanism maps the input features to multiple different attention subspaces, independently calculating attention weights between residue pairs in each subspace. This allows the encoder to simultaneously capture long-range evolutionary dependencies between different residue positions in the initial evolutionary feature matrix from multiple perspectives, such as co-evolutionary patterns and functional coupling. The feedforward neural network performs a position-by-position nonlinear transformation on the output of the self-attention mechanism, further enhancing the feature expressiveness. The evolutionary conservatism features are the set of feature vectors output after processing by the Transformer encoder, integrating the original evolutionary statistics and deep dependencies.

[0120] This embodiment constructs a complete evolutionary feature extraction workflow from multiple sequence alignment search to Transformer encoder refinement, achieving in-depth mining of protein evolutionary conservation information. The combination of the HHblits tool and the UniRef30 database ensures the depth and breadth of homology search, while a pseudo-counting smoothing strategy guarantees the robustness of frequency statistics. The Shannon entropy formula provides an intuitive quantitative indicator of conservation, and the Transformer encoder captures long-range evolutionary dependencies between residues through a self-attention mechanism, providing high-quality evolutionary conservation feature input for subsequent multimodal fusion.

[0121] Please see Figure 2In some embodiments, the three-dimensional structural data is used to calculate the residue dynamic communication score through an elastic network model, including:

[0122] S101. Obtain the three-dimensional structural data of the target protein, and extract the C-value of each residue from the three-dimensional structural data. β Atomic coordinates, for glycine residues, C is extracted. α Atomic coordinates are used as a substitute;

[0123] S102, C based on all residues α Or C β Atomic coordinates are used to calculate the spatial Euclidean distance between residue pairs. A distance truncation threshold is set. When the spatial Euclidean distance between residue pairs is less than the distance truncation threshold, a spring connection is constructed between the residue pairs to generate the Kirchhoff matrix of the elastic network model.

[0124] S103. Perform pseudo-inverse operation on the Kirchhoff matrix to generate a pseudo-inverse matrix. The diagonal elements of the pseudo-inverse matrix represent the theoretical mean square fluctuation of each residue under the elastic network model, and the off-diagonal elements of the pseudo-inverse matrix represent the fluctuation covariance between different residue pairs.

[0125] S104. Based on the off-diagonal and diagonal elements of the pseudo-inverse matrix, the Pearson correlation coefficient between each pair of residues is calculated using the normalized covariance formula to generate a residue dynamic correlation matrix.

[0126] S105. Perform row normalization on the residue dynamic correlation matrix, sum the dynamic correlation coefficients of each residue with all other residues and take the mean to generate the residue dynamic communication score of each residue. The residue dynamic communication score is used to quantify the degree of communication hub of each residue in the global dynamic network of the protein.

[0127] In step S101, the three-dimensional structure data of the target protein is stored in a PDB format file, in which the three-dimensional spatial coordinates of each atom are arranged line by line in a standardized recording format; the PDB file marked C is parsed. β The atomic entries are read and their three-dimensional spatial coordinates are obtained. For glycine residues, since their side chains consist of only one hydrogen atom and therefore do not contain C atoms... β Atoms, therefore extract C α Atomic coordinates are used as a substitute, C α The atom is the main chain carbon atom common to all amino acids, and its coordinates are always present in the PDB file. The obtained set of coordinates serves as the input for the node positions in the subsequent construction of the elastic network model.

[0128] In step S102, the three-dimensional Euclidean distance between any two residues is obtained by taking the square root of the sum of the squares of the differences in the atomic coordinates of the two residues. The distance truncation threshold can be set according to the modeling conventions of protein elastic network models to determine whether there is an effective elastic interaction between two residues. When the spatial Euclidean distance between two residues is less than the threshold, it is considered that there is a direct mechanical coupling relationship between them, and a virtual spring connection is constructed between them. The stiffness coefficient of the spring can be set according to the distance. The Kirchhoff matrix is ​​an N×N symmetric matrix (N is the number of residues). Its diagonal elements are the sum of the number of springs connected to each residue, and the off-diagonal elements are negative spring stiffness coefficients at residue pairs with spring connections and zero values ​​at unconnected residue pairs. This matrix completely describes the topology and connection strength distribution of the protein elastic network.

[0129] In step S103, the Kirchhoff matrix is ​​a singular matrix (its row sum is zero), and its generalized inverse matrix can be calculated using singular value decomposition or the Moore-Penrose pseudo-inverse algorithm. The diagonal elements of the pseudo-inverse matrix correspond to the theoretical mean square fluctuation of each residue under the elastic network model. This value reflects the vibration amplitude of the residue around its equilibrium position under thermodynamic equilibrium. The larger the value, the higher the flexibility of the residue. The off-diagonal elements of the pseudo-inverse matrix represent the fluctuation covariance between different residue pairs. Positive values ​​indicate that the two residues tend to move in the same direction, and negative values ​​indicate that they tend to move in opposite directions. The absolute value reflects the coupling strength.

[0130] In step S104, when calculating the Pearson correlation coefficient between each pair of residues using the normalized covariance formula, the Pearson correlation coefficient is obtained by dividing the fluctuation covariance between the residue pair by the product of the square roots of the theoretical mean square fluctuations of the two residues. The resulting value is between negative one and positive one. The residue dynamic correlation matrix is ​​an N×N symmetric matrix composed of the Pearson correlation coefficients between all residue pairs, where the diagonal elements are 1 and the off-diagonal elements reflect the consistency of the motion direction and the coupling strength of the corresponding residue pairs in the global dynamic motion of the protein.

[0131] In step S105, for each residue, the absolute value or original value of its dynamic correlation coefficient with all other residues is summed, and then divided by the total number of residues minus one to obtain the mean value, which is the residue dynamic communication score of that residue. The residue dynamic communication score is used to quantify the degree of communication hub of each residue in the global dynamic network of the protein. The higher the value, the stronger the overall dynamic coupling of the residue with other residues, and the more likely it is to play the role of a key communication node in the functional conformational changes of the protein.

[0132] This embodiment establishes an elastic network model by extracting the atomic coordinates of residues, obtains residue fluctuation information using Kirchhoff matrix pseudo-inverse operations, and further calculates the dynamic correlation between residues through normalized covariance, ultimately aggregating the communication hub index for each residue. This process transforms the three-dimensional structural information of the protein into quantitative features reflecting the dynamic coupling relationship between residues, providing a basis for communication strength based on the intrinsic dynamics of the protein for subsequent dynamic updates of the graph structure.

[0133] In some embodiments, the spatial contact relationships of residues are extracted to generate structural graph topological features, including:

[0134] Obtain the three-dimensional structural data of the target protein, and extract the C-value of each residue from the three-dimensional structural data. β Atomic coordinates, for glycine residues, C is extracted. α Atomic coordinates are used as a substitute;

[0135] C based on all residues α Atomic coordinates or C β Atomic coordinates are used to calculate the spatial Euclidean distance between residue pairs. A spatial contact distance threshold is set. When the spatial Euclidean distance between residue pairs is less than the spatial contact distance threshold, it is determined that there is a spatial contact relationship between the residue pairs, and a set of spatial contact edges is generated.

[0136] Obtain the sequence position index of residues in the amino acid sequence data. When the absolute value of the difference between the sequence position indices of two residues is a preset value, it is determined that there is a sequence adjacency relationship between the residue pairs, and a set of sequence adjacency edges is generated.

[0137] The set of spatial contact edges and the set of sequential adjacent edges are merged to generate a heterogeneous edge set, which contains two edge types: spatial contact edges and sequential adjacent edges.

[0138] Construct an edge feature vector for each edge in the heterogeneous edge set. The edge feature vector shall contain at least the normalized distance, the direction cosine component, and the edge type one-hot encoding.

[0139] An initial heterogeneous graph structure is constructed using residues as nodes, heterogeneous edge sets as edges, and edge feature vectors as edge attributes.

[0140] Obtain the Laplacian matrix of the initial heterogeneous graph structure, perform eigenvalue decomposition on the Laplacian matrix, extract the eigenvectors corresponding to the minimum non-zero eigenvalues ​​of the first preset number, generate the Laplacian position code, and attach the Laplacian position code as the node position feature to each node of the initial heterogeneous graph structure to generate the topological features of the structure graph.

[0141] In this embodiment, the coordinate set obtained from the three-dimensional structural data of the target protein is used to subsequently determine the proximity of residue pairs in three-dimensional space, providing a geometric data basis for constructing a graph structure that reflects the spatial folding pattern of the protein.

[0142] The spatial contact distance threshold is set based on the effective range of non-bonded interactions between residues in the protein structure. When the spatial Euclidean distance between two residues is less than the threshold, they are considered to have a direct physical proximity relationship in the three-dimensional folded space, and this determination result is recorded as a spatial contact edge. After performing the above determination on all residue pairs one by one, the edges that satisfy the spatial proximity condition constitute a set of spatial contact edges. This set encodes the local spatial adjacency pattern between residues in the three-dimensional structure of the protein.

[0143] The sequence number of each residue along the polypeptide chain from the amino terminus to the carboxyl terminus is obtained from the amino acid sequence data. This number is the sequence position index. The preset value is set according to the natural property of covalent linkage in the protein backbone. When the absolute value of the difference between the sequence position indices of two residues is equal to the preset value, they are considered to be directly connected by peptide bonds on the polypeptide chain, and this relationship is recorded as a sequence adjacent edge. After performing the above judgment on all residue pairs one by one, the edges that satisfy the sequence adjacent condition constitute the sequence adjacent edge set. This set preserves the covalent linkage topology between residues in the primary structure of the protein.

[0144] The set of spatial contact edges and the set of sequential adjacent edges are merged into a unified edge set, where each edge is accompanied by a label identifying its origin, which is used to distinguish whether the edge belongs to a spatial proximity relationship or a sequential proximity relationship. The merged set is the heterogeneous edge set, which provides a structural basis for the subsequent graph neural network to distinguish and process different types of residue connections.

[0145] For each edge in the heterogeneous edge set, an edge feature vector is constructed. The normalized distance is obtained by dividing the spatial Euclidean distance between the two residues connected by the edge by the maximum value among all residue pairs, so that the distance value is mapped to a uniform numerical range. The direction cosine component is composed of the cosine values ​​of the angle between the atomic coordinate difference vector of the two residues and the three axes in three-dimensional space, which is used to describe the relative orientation of the two residues in space. The edge type one-hot encoding is represented by a vector with a length equal to the number of edge types. It takes a value of 1 at the position of the corresponding edge type and a value of zero at other positions, which is used to distinguish between spatial contact edges and sequential adjacent edges.

[0146] Using each residue as a node in the graph structure, each edge in the heterogeneous edge set as the connection between nodes, and the edge feature vector corresponding to each edge as the attribute information of the edge, an initial heterogeneous graph structure is constructed. This structure simultaneously encodes the protein's sequence covalent connection topology and three-dimensional spatial proximity topology.

[0147] The Laplacian matrix is ​​calculated from the initial heterogeneous graph structure. This matrix, obtained by the difference between the degree matrix and the adjacency matrix, reflects the connectivity degree of each node in the graph structure and the connection relationships between nodes. Eigenvalue decomposition is performed on the Laplacian matrix to extract the eigenvectors corresponding to the smallest non-zero eigenvalues ​​of the first preset number. These eigenvectors constitute the Laplacian position code, which is used to assign a low-dimensional vector representation of each node reflecting its global position in the graph structure. The Laplacian position code is then appended as a node position feature to each node in the initial heterogeneous graph structure to generate a structural graph topological feature. This feature integrates the connection topological relationships between nodes, the multidimensional geometric properties of edges, and the global position encoding information of nodes.

[0148] This embodiment extracts various residue relationships from protein 3D structure data and amino acid sequence data, constructing a rich graph structure representation that includes heterogeneous edge types, multidimensional edge feature vectors, and Laplace position encoding. The spatial contact edge set captures local spatial adjacency patterns in protein 3D folding, the sequence adjacent edge set preserves the main chain covalent connection topology of the polypeptide chain, the heterogeneous edge structure provides a basis for distinguishing different types of residue relationships, the normalized distance and direction cosine components in the edge feature vectors enrich the geometric information representation of the edges, and the Laplace position encoding endows nodes with global position awareness through graph theory. This structural graph topological feature provides a graph structure input containing rich topological and geometric information for subsequent multimodal fusion and dynamic graph iterative updates.

[0149] In some embodiments, amino acid sequence data is input into a protein language model for encoding to generate sequence features, including:

[0150] The amino acid sequence data of the target protein is obtained, and the amino acid sequence data is standardized by replacing non-standard amino acids with corresponding standard amino acids to generate a standardized amino acid sequence.

[0151] The standardized amino acid sequence is input into the ESM-2 protein language model, which contains a multi-layer Transformer encoder. The network parameters of the first preset number of Transformer encoders of the ESM-2 protein language model are frozen, and the network parameters of the last preset number of Transformer encoders are fine-tuned to generate the ESM-2 hidden layer representation.

[0152] The ESM-2 hidden layer representation is input into a multi-scale one-dimensional convolution module, which contains multiple parallel one-dimensional convolution kernels, each with a different kernel size. Multi-scale local pattern extraction is performed on the ESM-2 hidden layer representation through multiple parallel one-dimensional convolution kernels to generate multi-scale convolution features.

[0153] The convolutional features at each scale in the multi-scale convolutional features are concatenated to generate concatenated convolutional features.

[0154] The concatenated convolutional features are input into a layer normalization layer and a fully connected layer for dimensionality transformation to generate sequence features.

[0155] In this embodiment, the standardization process is performed according to a preset mapping table between non-standard amino acids and standard amino acids. This mapping table is determined based on the similarity of the side chain chemical structure of amino acids and their correlation with biosynthetic pathways. For example, selenocysteine ​​is mapped to cysteine ​​and pyrrolidone is mapped to lysine. By scanning each residue identifier in the amino acid sequence data one by one, if a residue not found in the twenty standard amino acids is found, it is replaced with the corresponding standard amino acid to generate a standardized amino acid sequence.

[0156] The ESM-2 protein language model consists of multiple stacked Transformer encoder layers. During its pre-training phase, protein sequences from the UniRef50 database are used as training data, employing a masked language modeling approach to learn the contextual representation of each residue in the sequence. Before freezing, the network parameters of the pre-set number of Transformer encoder layers are obtained by setting the gradients of the weight matrices and bias terms of these layers to zero during backpropagation; these layers do not participate in parameter updates. When fine-tuning the network parameters of the subsequent pre-set number of Transformer encoder layers, labeled data from the binding site prediction task are used as training samples. The prediction error is calculated using the cross-entropy loss function, and the parameters of these layers are updated via backpropagation to generate the ESM-2 hidden layer representation.

[0157] In the multi-scale one-dimensional convolution module, multiple parallel one-dimensional convolution kernels independently perform sliding window convolution operations on the ESM-2 hidden layer representation along the sequence dimension. The size of different convolution kernels determines the number of consecutive residues covered by each convolution operation. Smaller convolution kernels focus on local interaction patterns between neighboring residues, while larger convolution kernels cover a wider sequence window to capture short-range dependencies between slightly more distant residues. The output of each convolution kernel maintains the same length as the input in the sequence length dimension through padding operations, generating multi-scale convolution features.

[0158] The convolutional features at each scale in the multi-scale convolutional features are merged along the feature channel dimension, so that the feature vector corresponding to each residue position contains multi-scale local pattern information extracted from different convolutional kernels, generating spliced ​​convolutional features. The spliced ​​convolutional features are then input to a layer normalization layer, which calculates the mean and variance along the feature channel dimension and standardizes the features. The features are then input to a fully connected layer and linearly transformed to map the feature dimension to a preset uniform dimension, generating sequence features.

[0159] This embodiment extracts global semantic representations using the ESM-2 protein language model, and captures local sequence patterns of different ranges by combining a multi-scale one-dimensional convolution module. The partial parameter freezing strategy reduces the risk of overfitting in downstream tasks while preserving pre-training knowledge. Layer normalization and fully connected layers uniformly map multi-source features to a standard dimensional space, providing consistent sequence feature inputs for subsequent multimodal fusion.

[0160] Please see Figure 3 In some embodiments, evolutionary conservation features, structural graph topological features, and sequence features are input into a pre-trained multimodal fusion encoder to generate a first fusion feature, including:

[0161] S201. Obtain evolutionary conserved features, structural graph topology features, and sequence features. Input the evolutionary conserved features, structural graph topology features, and sequence features into the corresponding self-attention modules. Each self-attention module contains a multi-head self-attention mechanism. Through the multi-head self-attention mechanism, the long-range dependencies between residues within each modality feature are strengthened, generating self-attention refined evolutionary conserved features, self-attention refined structural graph topology features, and self-attention refined sequence features.

[0162] S202. Obtain the SCOP fold type classification label of the target protein. If the SCOP fold type classification label cannot be obtained, use the preset general prior gating weight vector or predict its most likely fold type through an auxiliary classifier to generate the corresponding prior gating weight vector. Based on the SCOP fold type classification label, find the corresponding prior gating weight vector from the preset SCOP prior gating matrix. The rows of the SCOP prior gating matrix correspond to different SCOP fold types, and the columns of the SCOP prior gating matrix correspond to three modalities: evolutionary conservation features, structural graph topology features, and sequence features.

[0163] S203. The evolutionary conservatism features refined by self-attention, the topological features of the structure graph refined by self-attention, and the sequence features refined by self-attention are respectively used as key vectors. The globally learnable query vector is used as the query vector. The dot product attention score of the key vector and the query vector is calculated to generate the initial dynamic weights.

[0164] S204. Multiply the initial dynamic weights element-wise with the prior gated weight vector to generate the modulated dynamic weights. Perform softmax normalization on the modulated dynamic weights to generate normalized dynamic weights.

[0165] S205. Based on normalized dynamic weights, the evolutionary conservation features refined by self-attention, the topological features of the structure graph refined by self-attention, and the sequence features refined by self-attention are weighted and summed to generate residue-level fusion features.

[0166] S206. Using residue-level fusion features as query vectors, concatenate the evolutionary conservation features refined by self-attention, the topological features of the structure graph refined by self-attention, and the sequence features refined by self-attention as key-value pairs, and perform multi-head cross-attention operation to generate cross-attention refined features.

[0167] S207. Input the cross-attention refined features into the feedforward neural network for nonlinear transformation to generate the first fused features.

[0168] In step S201, the self-attention module maps the input features to a query matrix, a key matrix, and a value matrix, respectively. After calculating the dot product of the query matrix and the key matrix, it normalizes the result using softmax to obtain the attention weight matrix, which is then weighted and aggregated with the value matrix to refine the features. The multi-head self-attention mechanism executes the above process in parallel multiple times, projecting the input features to different subspaces each time using different linear mapping parameters. The outputs of each head are concatenated and linearly transformed to obtain the final output, enabling the model to simultaneously capture long-range dependencies between residues from multiple different representation subspaces. For example, the association patterns of co-evolutionary residue pairs in evolutionary conservation features, the dependency strength between spatially adjacent residues in structural graph topology features, and the semantic associations between distant residues in sequence features generate evolutionary conservation features refined by self-attention, structural graph topology features refined by self-attention, and sequence features refined by self-attention, respectively.

[0169] In step S202, the SCOP fold type classification label of the target protein is obtained by querying the SCOP database. This database classifies proteins according to the three-dimensional folding pattern of protein domains, and each protein domain is assigned a unique fold type identifier. Each row of the SCOP prior gating matrix corresponds to a SCOP fold type, and each column corresponds to three modalities: evolutionary conservation features, structural graph topology features, and sequence features. The matrix element values ​​are initialized based on the statistical mutual information of the contribution of each modality to the prediction of binding sites in historical data for each fold type of protein, and are updated along with the model as learnable parameters during end-to-end training. The row vector corresponding to the SCOP fold type classification label is found in this matrix, which is the prior gating weight vector.

[0170] In step S203, the modal features refined by self-attention are compressed into key vectors through linear mapping. The globally learnable query vector is a vector parameter that is randomly initialized and updated through backpropagation during training. The attention score corresponding to each modality is obtained by calculating the dot product of the query vector and the key vector of each modality. This score reflects the relative importance of each modality in the current prediction task, and the initial dynamic weights are generated.

[0171] In step S204, the initial dynamic weights and the prior gating weight vector are multiplied element-wise at corresponding element positions, so that the dynamic importance of the data-driven evaluation is modulated by the prior knowledge of the SCOP folding type. The multiplication result is normalized by the softmax function, and the weight values ​​of each modality are converted into a probability distribution with a sum of one, generating normalized dynamic weights.

[0172] In step S205, the weight value corresponding to each mode in the normalized dynamic weights is multiplied by the self-attention refined feature of that mode, and then the three weighted feature matrices are added element by element along the feature dimension, so that the fusion feature vector corresponding to each residue position integrates the information of the three modes and the contribution ratio of each mode is adaptively determined by the normalized dynamic weights, thereby generating residue-level fusion features.

[0173] In step S206, the multi-head cross-attention operation maps the query vector and key-value pairs to multiple attention subspaces respectively. In each subspace, the attention weights of the query and the key are calculated independently and the corresponding values ​​are aggregated. This enables each residue position in the residue-level fusion feature to be located from the original refined features of each modality and extract semantically related complementary information. Different attention heads focus on different aspects of cross-modal association patterns. The outputs of each head are concatenated and linearly transformed to obtain the cross-attention refined features.

[0174] In step S207, the feedforward neural network consists of two fully connected layers and a nonlinear activation function between them. The cross-attention refined features are input into the first fully connected layer for dimensional expansion. After the nonlinear activation function introduces nonlinear transformation capability, the feature dimension is restored to the preset unified dimension through the second fully connected layer to generate the first fused feature.

[0175] This embodiment achieves protein type-aware adaptive multimodal fusion by combining SCOP prior gating with dynamic weights. The self-attention module strengthens the residue dependencies within each modality, the SCOP prior gating matrix introduces prior knowledge of protein folding types, and the globally learnable query vector provides data-driven dynamic weight evaluation. The fusion of prior and dynamic approaches ensures that modal weights consider both biological priors and input data characteristics. Residue-level weighted summation achieves information integration at the residue level, cross-modal cross-attention further eliminates semantic differences between modalities, and the feedforward neural network enhances feature expressiveness. This fusion strategy overcomes the limitation of fixed-weight fusion methods in adapting to the differences in information dependencies among different protein types, achieving protein type-aware adaptive multimodal fusion and providing high-quality multimodal fusion features for subsequent prediction and decoding.

[0176] In some embodiments, the first fused feature is input into a pre-trained prediction decoder to generate an initial binding site prediction probability, including:

[0177] The first fusion feature is input into a bidirectional gated recurrent unit, which includes a forward-gated recurrent unit and a backward-gated recurrent unit. The forward-gated recurrent unit captures the contextual dependencies of residues along the amino acid sequence from the N-terminus to the C-terminus, and the backward-gated recurrent unit captures the contextual dependencies of residues along the amino acid sequence from the C-terminus to the N-terminus. The hidden outputs of the forward-gated recurrent unit and the backward-gated recurrent unit are spliced ​​together to generate sequence context features.

[0178] The sequence context features are input into a multilayer perceptron, which contains multiple sequentially connected fully connected layers and activation function layers. The sequence context features are then subjected to a layer-by-layer nonlinear transformation through the multiple sequentially connected fully connected layers and activation function layers to generate log odds values.

[0179] The log-odds value is input into the sigmoid activation function, which maps the log-odds value to a probability range between 0 and 1, generating the initial predicted probability of each residue belonging to the ortho-binding site.

[0180] In this embodiment, the forward-gated loop unit and the backward-gated loop unit in the bidirectional gated loop unit each maintain independent gating parameters and hidden layer state dimensions. The forward-gated loop unit reads the first fusion feature residue by residue along the amino acid sequence from the N-terminus to the C-terminus. At each time step, the current hidden layer output is calculated through update gates and reset gates based on the input feature of the current residue and the hidden layer state of the previous time step. The backward-gated loop unit reads the first fusion feature residue by residue along the amino acid sequence from the C-terminus to the N-terminus in the same manner. The hidden layer outputs in both directions are aligned along the sequence length dimension according to the residue position and then spliced ​​along the feature dimension, so that the sequence context feature dimension corresponding to each residue is the sum of the forward and backward hidden layer dimensions.

[0181] The specific configuration of multiple sequentially connected fully connected layers and activation function layers in the multilayer perceptron is determined based on the dimension of the first fused feature and the complexity of the prediction task. The first fully connected layer linearly transforms the dimension of the sequence context feature to a preset intermediate dimension. After the activation function layer introduces a nonlinear transformation, it is passed to the next fully connected layer. The last fully connected layer maps the feature dimension to a single numerical value, which is the log-odds value, representing the model's unnormalized confidence that the residue belongs to the ortho-binding site.

[0182] The sigmoid activation function takes the logarithmic probability value as input and compresses it to a probability range between 0 and 1 through the ratio operation of the exponential function. The closer the output value is to 1, the higher the confidence that the residue belongs to the ortho-binding site. The closer it is to 0, the lower the confidence. The calculation results of all residue positions constitute the initial predicted probability of the binding site.

[0183] In this embodiment, context information is extracted and concatenated from two directions of the sequence through a bidirectional gated recurrent unit, so that the representation of each residue incorporates the complete sequence context dependency. The multilayer perceptron converts the context features into classification decision scores through layer-by-layer linear transformation and nonlinear mapping. The sigmoid activation function maps continuous scores into predicted values ​​with probabilistic significance, providing an initial prediction basis for subsequent dynamic graph iteration updates.

[0184] Please see Figure 4 In some embodiments, the initial binding site prediction probability and residue dynamic communication score are iteratively updated according to a preset three-graph update rule to generate updated structural graph topological features, including:

[0185] S301. Obtain the initial binding site prediction probability, residue dynamic communication score and structural graph topology features. The structural graph topology features include node set, edge set and edge weight matrix.

[0186] S302. Set a first confidence threshold and a second confidence threshold. The first confidence threshold is higher than the second confidence threshold. Set a far-end distance threshold.

[0187] S303. Traverse each edge in the edge set and obtain the initial binding site prediction probability of each of the two residue nodes connected by each edge and the spatial Euclidean distance between the two residue nodes.

[0188] S304. When the initial binding site prediction probabilities of two residue nodes are both higher than the first confidence threshold, the edge weights corresponding to the edges are multiplied by the enhancement coefficient to generate updated edge weights. The enhancement coefficient is positively correlated with the product of the initial binding site prediction probabilities of the two residue nodes.

[0189] When the initial binding site prediction probabilities of two residue nodes are both lower than the second confidence threshold and the spatial Euclidean distance between the two residue nodes is greater than the far-end distance threshold, the edge weight corresponding to the edge is multiplied by the decay coefficient to generate the updated edge weight. The decay coefficient is a positive number less than 1.

[0190] S305. Traverse the dynamic communication scores of residues, obtain the dynamic correlation coefficient between each pair of residue nodes, set a communication score threshold, and when the dynamic correlation coefficient between a pair of residue nodes is higher than the communication score threshold and there is no edge connecting a pair of residue nodes in the edge set, add an edge to the edge set and assign an initial edge weight to the added edge.

[0191] S306. Replace the original edge set and original edge weight matrix in the topological features of the structure graph with the edge set and edge weight matrix after edge weight update and edge addition operations, and generate the updated topological features of the structure graph.

[0192] In step S301, the topological features of the structural graph are stored in a graph data structure. The node set records the index identifiers of all residues in the protein, the edge set records all residue pairs that are determined to have a connection relationship, and the edge weight matrix stores the weight value corresponding to each edge in the form of a two-dimensional array. The magnitude of the weight value reflects the strength of the mutual influence between the two residues connected by the edge.

[0193] In step S302, the first confidence threshold and the second confidence threshold are determined based on the distribution characteristics of the predicted probabilities on the training set. The first confidence threshold is set in the high quantile interval of the predicted probability distribution to screen high-confidence positive predicted residues, and the second confidence threshold is set in the low quantile interval of the predicted probability distribution to screen low-confidence negative predicted residues. The far-end distance threshold is set based on the upper limit of the typical range of effective spatial interactions between residues in the protein structure to distinguish between spatially adjacent and spatially distant residue pairs.

[0194] In step S303, when traversing each edge in the edge set, the corresponding probability value is read from the initial binding site prediction probability vector through the two residue node indices of the edge record, and the spatial Euclidean distance between the two residue nodes is read from the pre-calculated three-dimensional spatial distance matrix.

[0195] In step S304, when the initial binding site prediction probabilities of both residue nodes are higher than the first confidence threshold, the enhancement coefficient is calculated by multiplying the product of these two probability values ​​by a preset scaling factor. This scaling factor can be set according to the control requirements of the growth rate of edge weights during the iteration process, so that the enhanced edge weights can effectively amplify the signal propagation in the high confidence region without excessively dominating the graph structure during the iteration process. When the initial binding site prediction probabilities of both residue nodes are lower than the second confidence threshold and the spatial Euclidean distance between the two residue nodes is greater than the far-end distance threshold, the attenuation coefficient is a positive number less than 1. Its value can be determined according to the noise suppression rate expected in each iteration, so that the connection strength between low-confidence far-end residue pairs gradually attenuates during the iteration process, reducing the propagation of noise information in subsequent predictions.

[0196] In step S305, when traversing the dynamic communication scores of residues, the communication score threshold is set according to the distribution of the dynamic communication scores of residues across the entire protein, typically taking the mean or median of the communication scores of all residue pairs. When the dynamic correlation coefficient between a pair of residue nodes is higher than the communication score threshold, it is checked whether an edge connecting these two nodes already exists in the edge set. If not, a new edge is added to the edge set, and a small initial edge weight is assigned to the new edge, so that the newly introduced connection has a relatively mild impact on the prediction results in the initial stage, gradually taking effect as the number of iterations increases. The dynamic communication score of residues is used to quantify the degree of communication hub status of each residue in the global dynamic network of the protein. In subsequent graph update steps, the dynamic correlation coefficients between paired residues are directly derived from the off-diagonal elements of the dynamic correlation matrix of residues generated in this step. This mechanism provides a key graph structure foundation for allosteric site prediction, enabling residue pairs that are not spatially in direct contact but are functionally coupled through the internal dynamic network of the protein to establish information transmission channels.

[0197] In step S306, the edge set and edge weight matrix after edge weight update and edge addition operations are used to replace the original data structure in the topological features of the structure graph. This makes the updated topological features of the structure graph contain the adjusted edge connection relationships and edge weight distribution, providing the graph structure input optimized based on the current prediction results for the next round of feature fusion and prediction.

[0198] This embodiment achieves feedback optimization of the graph structure based on the prediction results through three graph update rules. The high-confidence residue-edge enhancement rule amplifies the propagation of the synergistic effect between residues at high-confidence binding sites, the low-confidence distal residue-edge decay rule suppresses noise interference in low-confidence regions, and the high-dynamic-communication-score residue-edge addition rule establishes new information transmission channels for functionally coupled but spatially non-contact residue pairs. The synergistic effect of these three rules gradually refines the graph structure during the iteration process, providing a topological foundation for continuous optimization in multiple rounds of prediction.

[0199] In some embodiments, based on the predicted probability of the final binding site, a gradient-weighted class activation mapping algorithm is used to calculate the contribution weight of each residue to the prediction result, generating a residue importance heatmap, including:

[0200] Obtain the final binding site prediction probability and the first fusion feature, which contains the fusion feature vector corresponding to each residue;

[0201] The predicted score corresponding to the target category in the final binding site prediction probability is used as the target score. The target category is the ortho-binding site category or the allosteric binding site category. Backpropagation is performed on the target score to calculate the gradient value of the target score relative to the fusion feature vector corresponding to each residue in the first fusion feature.

[0202] Global average pooling is performed on the gradient values ​​along the feature channel dimension to generate the neuron importance weights corresponding to each residue. The neuron importance weights represent the degree of contribution of each feature channel to the target score.

[0203] The fusion feature vector corresponding to each residue is weighted and summed with the neuron importance weight to generate the initial importance score for each residue;

[0204] The initial importance score is processed using the ReLU activation function, retaining the positive value part and suppressing the negative value part to generate the importance score after activation;

[0205] The importance score after activation is subjected to min-max normalization, which maps the importance score after activation to the interval between 0 and 1, generating a normalized residue importance score for each residue.

[0206] The normalized residue importance score is written as a temperature factor into the protein's three-dimensional structure file to generate a residue importance heatmap.

[0207] In this embodiment, the first fusion feature is stored in the form of a two-dimensional matrix, where the rows of the matrix correspond to the position of each residue in the protein sequence and the columns of the matrix correspond to the dimensions of the fusion feature vector; the final binding site prediction probability is stored in the form of a vector, where the length of the vector is equal to the number of residues and each element is the prediction probability value of the corresponding residue belonging to the orthogenetic binding site.

[0208] The predicted score corresponding to the orthogonal binding site category in the final binding site prediction probability is used as the target score. This target score is obtained by taking the maximum value of the prediction probability of all residue positions. When backpropagating the target score, the gradient is backpropagated layer by layer along the computation graph starting from the output layer of the prediction decoder. The gradient value of the target score relative to the fusion feature vector corresponding to each residue in the first fusion feature is calculated. The magnitude of the gradient value reflects the sensitivity of the fusion feature of the residue to the prediction result.

[0209] Global average pooling is performed on the gradient values ​​along the feature channel dimension. That is, the gradient values ​​of all residue positions in each feature channel are summed and then divided by the total number of residues to obtain the average gradient value corresponding to that channel, which is used as the neuron importance weight. The larger the weight value, the more significant the contribution of that feature channel to the prediction of the orthotype binding site category.

[0210] The fusion feature vector corresponding to each residue is weighted and summed with the neuron importance weight. That is, for each residue position, the value of each dimension in its fusion feature vector is multiplied by the neuron importance weight of the corresponding feature channel, and then the product of all dimensions is added together to obtain the initial importance score of the residue. The higher the score, the more important the residue plays in the prediction.

[0211] The initial importance scores are processed using the ReLU activation function, setting negative values ​​to zero and leaving positive values ​​unchanged to generate post-activation importance scores. This ensures that only residues that positively contribute to the prediction results are included in subsequent analyses. The post-activation importance scores are then subjected to min-max normalization, which involves finding the minimum and maximum scores of all residues, subtracting the minimum from the score of each residue, and dividing by the range to map all scores to the interval between 0 and 1, generating normalized residue importance scores.

[0212] The normalized residue importance score is written as a temperature factor column into the protein's 3D structure file. This involves parsing the original protein 3D structure file, locating the atomic record row corresponding to each residue, replacing the original value of the temperature factor field in that atomic record row with the normalized residue importance score, and saving the modified file to generate a residue importance heatmap. This file can be directly loaded into molecular visualization software, where the residue importance is visually mapped to the protein's 3D structure surface using a temperature factor coloring mode. This allows drug developers to directly observe which residue regions contribute most to binding site prediction, providing a visual basis for subsequent mutation experiment design or drug molecule docking site selection.

[0213] This embodiment applies the gradient-weighted activation mapping algorithm to protein binding site prediction. Through backpropagation, the prediction results of the deep neural network are traced back to the input residue level, generating a biologically interpretable residue importance heatmap. The gradient-weighted activation mapping algorithm, originally used for image classification interpretability, is adapted to the protein residue sequence prediction task. The first fusion feature is used as the target feature layer for gradient backpropagation, allowing the importance score of each residue to comprehensively reflect its contribution to the multimodal fusion feature space. Simultaneously, the normalized residue importance score is written into the protein's 3D structure file in temperature factor format, fully utilizing the native support of temperature factor coloring in molecular visualization software. This enables the prediction results to be directly presented in mainstream molecular visualization tools with zero additional development cost, providing residue-level decision-making basis for target screening and experimental verification in drug development practice.

[0214] Please see Figure 5In a second aspect, this embodiment also provides a system for predicting protein drug binding sites based on multimodal dynamic graphs, applicable to the method described in the first aspect. The system includes a data acquisition module, an evolutionary feature extraction module, a structural feature extraction module, a sequence feature extraction module, a multimodal fusion encoder, a prediction decoder, a graph topology iterative update module, an iteration control module, and an interpretable output module. The data acquisition module acquires the amino acid sequence data and three-dimensional structural data of the target protein. The evolutionary feature extraction module constructs a multi-sequence alignment file from the amino acid sequence data using a multi-sequence alignment tool, calculates the site-specific scoring matrix and residue conservation entropy, and generates evolutionarily conserved features. The structural feature extraction module calculates the residue dynamic communication score from the three-dimensional structural data using an elastic network model, extracts the spatial contact relationships of residues, and generates structural graph topology features. The sequence feature extraction module inputs the amino acid sequence data into a protein language model for encoding, generating sequence features. The multimodal fusion encoder integrates the evolutionarily conserved features, structural graph topology features, and sequence features... The algorithm employs an adaptive weighted fusion method to generate a first fused feature. A predictive decoder performs sequence context modeling and nonlinear transformation on this first fused feature to generate an initial binding site prediction probability. A graph topology iterative update module iteratively updates the structural graph topology feature with the initial binding site prediction probability and residue dynamic communication score according to three preset graph update rules, generating an updated structural graph topology feature. These three rules include a high-confidence residue edge enhancement rule, a low-confidence far-end residue edge attenuation rule, and a high-dynamic-communication-score residue edge addition rule. An iterative control module re-inputs the updated structural graph topology feature into the multimodal fusion encoder and predictive decoder, repeating the feature fusion and prediction steps until a preset number of iterations is reached, generating the final binding site prediction probability. An interpretable output module calculates the contribution weight of each residue to the prediction result based on the final binding site prediction probability using a gradient-weighted class activation mapping algorithm, generating a residue importance heatmap and outputting a prediction report containing the final binding site prediction probability and the residue importance heatmap.

[0215] Combination Figure 5 As shown, this system takes the amino acid sequence and three-dimensional structure file of the target protein (data acquisition module) as the whole input and enters the data preprocessing module (a sub-function of the data acquisition module) to perform preprocessing operations such as sequence standardization and structural data parsing. Preferably, the preprocessed data flows through a three-stage domain-specific self-supervised pre-training module (injecting biological prior knowledge of binding sites into the backbone network) and then enters the multimodal feature extraction module (M-FEM).

[0216] The Multimodal Feature Extraction Module (M-FEM) contains four parallel coding branches: the ESM-2 protein language model sequence coding branch (sequence feature extraction module), the MT-GAT heterogeneous graph structure coding branch (structural feature extraction module), the evolutionary Transformer conservation coding branch (evolutionary feature extraction module), and the physicochemical property MLP coding branch (preferred, which can provide additional physicochemical property features). After extracting features in parallel, the feature vectors output by each branch are passed to the next module.

[0217] The features are fed into the SCOP-gated hierarchical cross-modal attention fusion module (H-CMAF, Multimodal Fusion Encoder). This module sequentially includes an intramodal self-attention refinement unit, an SCOP prior gating matrix unit, a three-level hierarchical feature fusion unit, and a cross-modal cross-attention refinement unit, ultimately outputting the first fused feature. .

[0218] First fusion feature Thermodynamic stability sensing layering (TSL, which integrates protein thermodynamic dynamics information in a differentiable manner) then proceeds; TSL uses the three-source residue flexibility index (RFI) as input to modulate the amplitude of the fusion feature residue by residue, and outputs the modulated first fusion feature. .

[0219] First fusion feature after modulation The system then enters the Predict-Feedback Dynamic Graph Topology Iterative Update (IPGR) module. This module incorporates R1 high-confidence residue edge enhancement rules, R2 low-confidence far-end residue edge decay rules, and R3 high-dynamic communication score residue edge addition rules (graph topology iterative update module). Through an iterative control mechanism (iterative control module), the updated graph topology features are re-inputted to the multimodal fusion encoder (multimodal fusion encoder) and the predictor decoder (predictor decoder). The feature fusion and prediction steps are repeated until the preset number of iterations is reached, achieving closed-loop iterative update of the graph topology.

[0220] After the iterative update is completed, the modulated first fused feature generated in the last iteration enters the ENM subgraph morphological dual-head prediction decoding module (DH-PD, prediction decoder). This module contains a normal site prediction head and a morphological site prediction head, and uses a bidirectional GRU and multilayer perceptron structure to perform prediction decoding, generating the initial binding site prediction probability and the final binding site prediction probability.

[0221] Finally, the prediction results are fed into the Grad-CAM multidimensional interpretable output module (IOM). This module calculates the contribution weight of each residue to the prediction result based on the gradient-weighted class activation mapping algorithm (Grad-CAM), generates a residue importance heatmap, and outputs the final prediction results including ortho- / alteral site coordinates, five-dimensional confidence scores, and a structured JSON report.

[0222] Preferably, the Physically Constrained Virtual Pocket Augmentation (PC-VPA, enabled only during the training phase) performs sample balancing during model training to improve the imbalance between positive and negative samples. This module does not participate in the prediction phase. These modules are connected sequentially to form a complete end-to-end prediction system architecture.

[0223] By adopting the above technical solutions, this invention differs from existing technologies and possesses the following beneficial effects: It extracts protein features from three complementary dimensions—evolution, structure, and sequence—using multiple sequence alignment tools, elastic network models, and protein language models, providing a comprehensive and complete information foundation for binding site prediction; the gated attention fusion module in the multimodal fusion encoder uses SCOP folding type as a priori and adaptively adjusts the contribution weights of each modality for different protein types, achieving dynamic fusion with protein type awareness; the prediction-feedback dynamic graph topology iterative update mechanism refines the graph structure gradually during iteration through three graph update rules, forming a closed-loop optimization framework where prediction and graph structure mutually promote each other; the gradient-weighted class activation mapping algorithm backtracks the prediction results to the residue level and writes them into the protein's three-dimensional structure file in temperature factor format, providing prediction outputs that combine accuracy and interpretability. These technical solutions achieve high-precision prediction of protein drug binding sites, significantly improving the predictive model's adaptability to different protein types and the interpretability of the prediction results.

[0224] Finally, it should be noted that although the above embodiments have been described in the text and drawings of this application, this should not limit the scope of patent protection of this application. Any technical solutions that are based on the essential concept of this application and utilize the content described in the text and drawings of this application, resulting in equivalent structural or procedural substitutions or modifications, as well as the direct or indirect application of the technical solutions of the above embodiments to other related technical fields, are all included within the scope of patent protection of this application.

Claims

1. A method for predicting protein drug binding sites based on multimodal dynamic graphs, characterized in that, include: Obtain the amino acid sequence data and three-dimensional structure data of the target protein; The amino acid sequence data is used to construct a multiple sequence alignment file using a multiple sequence alignment tool, and the site-specific scoring matrix and residue conservation entropy value are calculated to generate evolutionary conservation features. The three-dimensional structural data are used to calculate the dynamic communication score of residues through an elastic network model, extract the spatial contact relationship of residues, and generate the topological features of the structural diagram. Amino acid sequence data is input into a protein language model for encoding to generate sequence features; Evolutionary conservation features, structural graph topological features, and sequence features are input into a pre-trained multimodal fusion encoder to generate the first fusion feature; The first fused feature is input into a pre-trained prediction decoder to generate the initial binding site prediction probability; The initial binding site prediction probability and residue dynamic communication score are used to iteratively update the topological features of the structural graph according to the three preset graph update rules to generate the updated topological features of the structural graph. The three graph update rules include the high confidence residue edge enhancement rule, the low confidence far residue edge decay rule, and the high dynamic communication score residue edge addition rule. The updated topological features of the structure graph are re-input into the multimodal fusion encoder and predictive decoder, and the feature fusion and prediction steps are repeated until the preset number of iterations is reached to generate the final binding site prediction probability. Based on the predicted probability of the final binding site, the contribution weight of each residue to the prediction result is calculated by the gradient weighted class activation mapping algorithm, a residue importance heatmap is generated, and a prediction report containing the predicted probability of the final binding site and the residue importance heatmap is output.

2. The method for predicting protein drug binding sites based on multimodal dynamic graphs according to claim 1, characterized in that, Amino acid sequence data are used to construct multiple sequence alignment files using multiple sequence alignment tools. Site-specific scoring matrices and residue conservation entropy values ​​are calculated to generate evolutionary conservation features, including: Obtain the amino acid sequence data of the target protein, input the amino acid sequence data into the HHblits tool, use the UniRef30 database as the search target, perform a multiple sequence alignment search, and generate a multiple sequence alignment file; Read the multiple sequence alignment file, count the frequency of amino acid types appearing at each residue position during the evolution process, use a pseudo-counting smoothing strategy to smooth the frequency distribution, and generate a site-specific scoring matrix; Based on the site-specific scoring matrix, the information entropy of each residue position is calculated using the Shannon entropy formula to generate the residue conservation entropy value; The site-specific scoring matrix is ​​concatenated with the residue conservation entropy value to form the initial evolutionary feature matrix; An initial evolutionary feature matrix is ​​input into a pre-trained Transformer encoder, which includes a multi-layer multi-head self-attention mechanism and a feedforward neural network. The multi-head self-attention mechanism captures the long-range evolutionary dependencies between residues in the initial evolutionary feature matrix to generate evolutionarily conserved features.

3. The method for predicting protein drug binding sites based on multimodal dynamic graphs according to claim 1, characterized in that, The dynamic communication score of residues is calculated using a resilient network model based on three-dimensional structural data, including: Obtain the three-dimensional structural data of the target protein, and extract the C-value of each residue from the three-dimensional structural data. β Atomic coordinates, for glycine residues, C is extracted. α Atomic coordinates are used as a substitute; C based on all residues α Or C β Atomic coordinates are used to calculate the spatial Euclidean distance between residue pairs. A distance truncation threshold is set. When the spatial Euclidean distance between residue pairs is less than the distance truncation threshold, a spring connection is constructed between the residue pairs to generate the Kirchhoff matrix of the elastic network model. The Kirchhoff matrix is ​​pseudo-inverseed to generate a pseudo-inverse matrix. The diagonal elements of the pseudo-inverse matrix represent the theoretical mean square fluctuations of each residue under the elastic network model, and the off-diagonal elements of the pseudo-inverse matrix represent the fluctuation covariance between different residue pairs. Based on the off-diagonal and diagonal elements of the pseudo-inverse matrix, the Pearson correlation coefficient between each pair of residues is calculated using the normalized covariance formula to generate a residue dynamic correlation matrix. The residue dynamic correlation matrix is ​​row normalized, and the average value of the summation of the dynamic correlation coefficients of each residue with all other residues is taken to generate the residue dynamic communication score of each residue. The residue dynamic communication score is used to quantify the degree of communication hub of each residue in the global dynamic network of the protein.

4. The method for predicting protein drug binding sites based on multimodal dynamic graphs according to claim 3, characterized in that, Extracting spatial contact relationships of residues to generate topological features of the structural diagram, including: Obtain the three-dimensional structural data of the target protein, and extract the C-value of each residue from the three-dimensional structural data. β Atomic coordinates, for glycine residues, C is extracted. α Atomic coordinates are used as a substitute; C based on all residues α Atomic coordinates or C β Atomic coordinates are used to calculate the spatial Euclidean distance between residue pairs. A spatial contact distance threshold is set. When the spatial Euclidean distance between residue pairs is less than the spatial contact distance threshold, it is determined that there is a spatial contact relationship between the residue pairs, and a set of spatial contact edges is generated. Obtain the sequence position index of the residues in the amino acid sequence data. When the absolute value of the difference between the sequence position indices of two residues is a preset value, it is determined that there is a sequence adjacency relationship between the residue pairs, and a set of sequence adjacency edges is generated. The set of spatial contact edges and the set of sequential adjacent edges are merged to generate a heterogeneous edge set, which contains two edge types: spatial contact edges and sequential adjacent edges. For each edge in the heterogeneous edge set, an edge feature vector is constructed, and the edge feature vector includes at least normalized distance, direction cosine component, and edge type one-hot encoding; Using the residues as nodes, the heterogeneous edge set as edges, and the edge feature vectors as edge attributes, an initial heterogeneous graph structure is constructed. Obtain the Laplacian matrix of the initial heterogeneous graph structure, perform eigenvalue decomposition on the Laplacian matrix, extract the feature vectors corresponding to the minimum non-zero eigenvalues ​​of the first preset number, generate Laplacian position codes, and attach the Laplacian position codes as node position features to each node of the initial heterogeneous graph structure to generate structural graph topological features.

5. The method for predicting protein drug binding sites based on multimodal dynamic graphs according to claim 1, characterized in that, Amino acid sequence data is input into a protein language model for encoding to generate sequence features, including: The amino acid sequence data of the target protein is obtained, and the amino acid sequence data is standardized by replacing non-standard amino acids with corresponding standard amino acids to generate a standardized amino acid sequence. The standardized amino acid sequence is input into the ESM-2 protein language model, which contains a multi-layer Transformer encoder. The network parameters of the first preset number of Transformer encoders of the ESM-2 protein language model are frozen, and the network parameters of the last preset number of Transformer encoders are fine-tuned and trained to generate the ESM-2 hidden layer representation. The ESM-2 hidden layer representation is input into a multi-scale one-dimensional convolution module, which contains multiple parallel one-dimensional convolution kernels, each with a different kernel size. Multi-scale local pattern extraction is performed on the ESM-2 hidden layer representation through multiple parallel one-dimensional convolution kernels to generate multi-scale convolution features. The convolutional features at each scale in the multi-scale convolutional features are concatenated to generate concatenated convolutional features. The concatenated convolutional features are input into a layer normalization layer and a fully connected layer for dimensionality transformation to generate sequence features.

6. The method for predicting protein drug binding sites based on multimodal dynamic graphs according to claim 1, characterized in that, Evolutionary conservation features, structural graph topological features, and sequence features are input into a pre-trained multimodal fusion encoder to generate the first fusion feature, including: The evolutionary conservation features, the structural graph topology features, and the sequence features are obtained. The evolutionary conservation features, the structural graph topology features, and the sequence features are respectively input into the corresponding self-attention modules. Each self-attention module contains a multi-head self-attention mechanism. The multi-head self-attention mechanism is used to strengthen the long-range dependencies between residues within each modality feature, thereby generating the self-attention refined evolutionary conservation features, the self-attention refined structural graph topology features, and the self-attention refined sequence features. Obtain the SCOP fold type classification label of the target protein. If the SCOP fold type classification label cannot be obtained, use a preset general prior gating weight vector or predict its most likely fold type through an auxiliary classifier to generate the corresponding prior gating weight vector. Based on the SCOP fold type classification label, the corresponding prior gating weight vector is found from the preset SCOP prior gating matrix. The rows of the SCOP prior gating matrix correspond to different SCOP fold types, and the columns of the SCOP prior gating matrix correspond to the three modalities of evolutionary conservation features, structural graph topology features, and sequence features. The evolutionary conservatism features refined by self-attention, the topological features of the structure graph refined by self-attention, and the sequence features refined by self-attention are used as key vectors, and the globally learnable query vector is used as the query vector. The dot product attention score of the key vector and the query vector is calculated to generate the initial dynamic weights. The initial dynamic weights are multiplied element-wise with the prior gated weight vector to generate modulated dynamic weights. The modulated dynamic weights are then normalized using softmax to generate normalized dynamic weights. Based on the normalized dynamic weights, the evolutionary conservation features refined by self-attention, the topological features of the structure graph refined by self-attention, and the sequence features refined by self-attention are weighted and summed to generate residue-level fusion features. Using the residue-level fusion feature as the query vector, the evolutionary conservation feature refined by self-attention, the topological feature of the structure graph refined by self-attention, and the sequence feature refined by self-attention are concatenated as key-value pairs, and multi-head cross-attention operation is performed to generate cross-attention refined features. The cross-attention refined features are input into a feedforward neural network for nonlinear transformation to generate the first fused features.

7. The method for predicting protein drug binding sites based on multimodal dynamic graphs according to claim 1, characterized in that, The first fused feature is input into a pre-trained prediction decoder to generate initial binding site prediction probabilities, including: The first fusion feature is input into a bidirectional gated loop unit, which includes a forward gated loop unit and a backward gated loop unit. The forward gated loop unit captures the contextual dependencies of residues along the amino acid sequence from the N-terminus to the C-terminus, and the backward gated loop unit captures the contextual dependencies of residues along the amino acid sequence from the C-terminus to the N-terminus. The hidden outputs of the forward gated loop unit and the backward gated loop unit are spliced ​​together to generate a sequence context feature. The sequence context features are input into a multilayer perceptron, which includes multiple sequentially connected fully connected layers and activation function layers. The sequence context features are then subjected to a layer-by-layer nonlinear transformation through the multiple sequentially connected fully connected layers and activation function layers to generate a log-odds value. The log-odds value is input into the sigmoid activation function, which maps the log-odds value to a probability range between 0 and 1, generating the initial predicted probability of each residue belonging to the ortho-binding site.

8. The method for predicting protein drug binding sites based on multimodal dynamic graphs according to claim 1, characterized in that, The initial binding site prediction probability and residue dynamic communication score are used to iteratively update the topological features of the structural graph according to a preset three-graph update rule, generating updated topological features of the structural graph, including: The initial binding site prediction probability, the residue dynamic communication score, and the structural graph topology features are obtained, wherein the structural graph topology features include a node set, an edge set, and an edge weight matrix. Set a first confidence threshold and a second confidence threshold, where the first confidence threshold is higher than the second confidence threshold, and set a far-end distance threshold; Traverse each edge in the edge set and obtain the initial binding site prediction probability of each of the two residue nodes connected by each edge and the spatial Euclidean distance between the two residue nodes; When the initial binding site prediction probabilities of both residue nodes are higher than the first confidence threshold, the edge weight corresponding to the edge is multiplied by an enhancement coefficient to generate an updated edge weight. The enhancement coefficient is positively correlated with the product of the initial binding site prediction probabilities of the two residue nodes. When the initial binding site prediction probabilities of two residue nodes are both lower than the second confidence threshold and the spatial Euclidean distance between the two residue nodes is greater than the far-end distance threshold, the edge weight corresponding to the edge is multiplied by the decay coefficient to generate the updated edge weight, wherein the decay coefficient is a positive number less than 1. Traverse the residue dynamic communication scores, obtain the residue dynamic correlation coefficient between each pair of residue nodes, set a communication score threshold, and when the residue dynamic correlation coefficient between a pair of residue nodes is higher than the communication score threshold and there is no edge connecting the pair of residue nodes in the edge set, add an edge to the edge set and assign an initial edge weight to the added edge. The edge set and edge weight matrix after edge weight update and edge addition operations are used to replace the original edge set and original edge weight matrix in the topological feature of the structure graph to generate the updated topological feature of the structure graph.

9. The method for predicting protein drug binding sites based on multimodal dynamic graphs according to claim 1, characterized in that, Based on the predicted probability of the final binding site, the contribution weight of each residue to the prediction result is calculated using a gradient-weighted class activation mapping algorithm, generating a residue importance heatmap, including: The predicted probability of the final binding site and the first fusion feature are obtained, wherein the first fusion feature contains a fusion feature vector corresponding to each residue; The predicted score corresponding to the target category in the final binding site prediction probability is used as the target score. The target category is either an orthogonal binding site category or an allosteric binding site category. The target score is backpropagated to calculate the gradient value of the target score relative to the fusion feature vector corresponding to each residue in the first fusion feature. The gradient value is subjected to global average pooling along the feature channel dimension to generate the neuron importance weight corresponding to each residue. The neuron importance weight represents the degree of contribution of each feature channel to the target score. The fusion feature vector corresponding to each residue is weighted and summed with the importance weight of the neuron to generate an initial importance score for each residue; The initial importance score is processed using the ReLU activation function, retaining the positive value part and suppressing the negative value part to generate the activated importance score; The activated importance score is subjected to min-max normalization to map the activated importance score to a range between 0 and 1, generating a normalized residue importance score for each residue. The normalized residue importance score is written as a temperature factor column into the protein's three-dimensional structure file to generate a residue importance heatmap.

10. A system for predicting protein drug binding sites based on multimodal dynamic graphs, characterized in that, The system applicable to the method of any one of claims 1 to 9 comprises: The data acquisition module is used to acquire the amino acid sequence data and three-dimensional structure data of the target protein; The evolutionary feature extraction module is used to construct a multiple sequence alignment file from amino acid sequence data using a multiple sequence alignment tool, calculate the site-specific scoring matrix and residue conservation entropy value, and generate evolutionary conservation features. The structural feature extraction module is used to calculate the dynamic communication score of residues through an elastic network model, extract the spatial contact relationship of residues, and generate the topological features of the structural diagram. The sequence feature extraction module is used to input amino acid sequence data into a protein language model for encoding and to generate sequence features; A multimodal fusion encoder is used to adaptively weight and fuse evolutionary conservatism features, structural graph topology features, and sequence features to generate the first fused feature; The predictive decoder is used to perform sequence context modeling and nonlinear transformation on the first fused features to generate the initial binding site prediction probability. The graph topology iterative update module is used to iteratively update the topological features of the structural graph according to the initial binding site prediction probability and the residue dynamic communication score according to three preset graph update rules, and generate updated structural graph topological features. The three graph update rules include the high confidence residue edge enhancement rule, the low confidence far residue edge decay rule, and the high dynamic communication score residue edge addition rule. The iterative control module is used to re-input the updated topological features of the structure graph into the multimodal fusion encoder and predictive decoder, and repeat the feature fusion and prediction steps until the preset number of iterations is reached to generate the final binding site prediction probability. An interpretable output module is used to calculate the contribution weight of each residue to the prediction result based on the final binding site prediction probability using a gradient-weighted class activation mapping algorithm, generate a residue importance heatmap, and output a prediction report containing the final binding site prediction probability and the residue importance heatmap.

Images